internal_step.f90

Go to the documentation of this file.

!/===========================================================================/
! Copyright (c) 2007, The University of Massachusetts Dartmouth 
! Produced at the School of Marine Science & Technology 
! Marine Ecosystem Dynamics Modeling group
! All rights reserved.
!
! FVCOM has been developed by the joint UMASSD-WHOI research team. For 
! details of authorship and attribution of credit please see the FVCOM
! technical manual or contact the MEDM group.
!
! 
! This file is part of FVCOM. For details, see http://fvcom.smast.umassd.edu 
! The full copyright notice is contained in the file COPYRIGHT located in the 
! root directory of the FVCOM code. This original header must be maintained
! in all distributed versions.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
! AND ANY EXPRESS OR  IMPLIED WARRANTIES, INCLUDING,  BUT NOT  LIMITED TO,
! THE IMPLIED WARRANTIES OF MERCHANTABILITY AND  FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED.  
!
!/---------------------------------------------------------------------------/
! CVS VERSION INFORMATION
! $Id$
! $Name$
! $Revision$
!/===========================================================================/

SUBROUTINE internal_step
  USE mod_nesting
  USE all_vars
  USE mod_utils
  USE mod_obcs
  USE mod_time
  USE eqs_of_state
  USE mod_wd
  USE mod_assim
  USE mod_par
  USE mod_ice
  USE mod_ice2d
  USE mod_northpole



  IMPLICIT NONE
  INTEGER :: I,J,K
  REAL(SP) :: UTMP,VTMP
  integer :: i1,i2
!------------------------------------------------------------------------------|

  if(dbg_set(dbg_sbr)) write(ipt,*)&
       & "Start: internal_step"

!----SET RAMP FACTOR TO EASE SPINUP--------------------------------------------!
  ramp = 0.0_sp
  IF(iramp /= 0) THEN
     ramp=tanh(real(iint,sp)/real(iramp,sp))
  ELSE
     ramp = 1.0_sp
  END IF

!----OFFLINE SEDIMENT UPDATE
!----loop end for offline sediment



!----SET UP WATER QUALITY MODEL COEFFICIENTS-----------------------------------!


  

  !----ADJUST CONSISTENCY BETWEEN 3-D VELOCITY AND VERT-AVERAGED VELOCITIES------!
    CALL adjust2d3d(1)

  !----SPECIFY THE SOLID BOUNDARY CONDITIONS OF U&V INTERNAL MODES---------------!
    CALL bcond_gcn(5,0)

! end defined semi-implicit && nh

  !----SPECIFY THE SURFACE FORCING OF INTERNAL MODES-----------------------------!
  CALL bcond_gcn(8,0)

  ! New Open Boundary Condition ----3

  !end !defined (TWO_D_MODEL)
  
  !----SPECIFY THE BOTTOM ROUGHNESS AND CALCULATE THE BOTTOM STRESSES------------!
  CALL bottom_roughness


!==============================================================================!
!  CALCULATE DISPERSION (GX/GY) AND BAROCLINIC PRESSURE GRADIENT TERMS         !
!==============================================================================!

    CALL advection_edge_gcn(advx,advy)          !Calculate 3-D Adv/Diff       !

  IF(recalculate_rho_mean) THEN
     IF(recalc_rho_mean .LT. inttime)THEN
        recalc_rho_mean = recalc_rho_mean + delt_rho_mean
        CALL rho_pmean    
     END IF
  END IF


  IF(.NOT. barotropic)THEN                    !Barotropic Flow ?            !
     SELECT CASE(baroclinic_pressure_gradient)
     CASE ("sigma levels")
        CALL baropg      !Sigma Level Pressure Gradient!
     CASE('z coordinates')
        CALL phy_baropg  !Z Level Pressure Gradient    !
     CASE DEFAULT
        CALL fatal_error("UNKNOW BAROCLINIC PRESURE GRADIENT TYPE",&
             & trim(baroclinic_pressure_gradient))
     END SELECT

  END IF                                      !                             !


  adx2d = 0.0_sp ; ady2d = 0.0_sp             !Initialize GX/GY Terms       !
  drx2d = 0.0_sp ; dry2d = 0.0_sp             !Initialize BCPG for Ext Mode !

  DO k=1,kbm1
     DO i=1, n
        adx2d(i)=adx2d(i)+advx(i,k)   !*DZ1(I,K)
        ady2d(i)=ady2d(i)+advy(i,k)   !*DZ1(I,K)
        drx2d(i)=drx2d(i)+drhox(i,k)  !*DZ1(I,K)
        dry2d(i)=dry2d(i)+drhoy(i,k)  !*DZ1(I,K)
     END DO
  END DO

  CALL advave_edge_gcn(advua,advva)           !Compute Ext Mode Adv/Diff
  adx2d = adx2d - advua                       !Subtract to Form GX
  ady2d = ady2d - advva                       !Subtract to Form GY

  !----INITIALIZE ARRAYS USED TO CALCULATE AVERAGE UA/E  OVER EXTERNAL STEPS-----!
  uard = 0.0_sp
  vard = 0.0_sp
  egf  = 0.0_sp


!!# if defined (1)
!!  UARDS = 0.0_SP
!!  VARDS = 0.0_SP
!!# endif

  IF(iobcn > 0) THEN
     uard_obcn(1:iobcn)=0.0_sp
  END IF
! end defined semi-implicit
! end defined (TWO_D_MODEL)


  ! New Open Boundary Condition ----4



  !==============================================================================!
  !  LOOP OVER EXTERNAL TIME STEPS                                               !
  !==============================================================================!
  DO iext=1,isplit

     IF (dbg_set(dbg_sbrio)) WRITE(ipt,*) "/// EXT SETP: ",iext

     exttime = exttime + imdte

     CALL external_step
     
     
  END DO


!  new update  !  jqi ggao 0730/2007  for E-P

!  new update  !  jqi ggao 0730/2007

! end defined semi-implicit

!==============================================================================!
!==============================================================================!
!                     BEGIN THREE D ADJUSTMENTS
!==============================================================================!
!==============================================================================!

  !==============================================================================!
  !    ADJUST INTERNAL VELOCITY FIELD TO CORRESPOND TO EXTERNAL                  !
  !==============================================================================!

    CALL adjust2d3d(2)

  ! New Open Boundary Condition ----9
  
  !==============================================================================!
  !     CALCULATE INTERNAL VELOCITY FLUXES                                       |
  !==============================================================================!
  !                                                    !
!!

!end !defined (TWO_D_MODEL)

!==============================================================================!
!     CALCULATE INTERNAL VELOCITY FLUXES                                       |
!==============================================================================!


    CALL vertvl_edge     ! Calculate/Update Sigma Vertical Velocity (Omega)   !

    IF(wetting_drying_on) CALL wd_update(2)

    CALL viscof_h        ! Calculate horizontal diffusion coefficient scalars !
  
    CALL adv_uv_edge_gcn ! Horizontal Advect/Diff + Vertical Advection        !

     CALL vdif_uv      ! Implicit Integration of Vertical Diffusion of U/V  !

    IF(adcor_on) THEN
      CALL adcor
     CALL vdif_uv      ! Implicit Integration of Vertical Diffusion of U/V  !
!      CALL VDIF_UV   ! Implicit Integration of Vertical Diffusion of U/V  !

    ENDIF

    DO i=1,n
      IF(h1(i) <= static_ssh_adj ) THEN
        DO k=1,kbm1
           uf(i,k)=ua(i)
           vf(i,k)=va(i)
        END DO
      END IF
    END DO

    CALL bcond_gcn(3,0)    ! Boundary Condition on U/V At River Input           !




    IF(nesting_on )THEN
      CALL set_var(inttime,u=uf)
      CALL set_var(inttime,v=vf)
    END IF  


!if !defined (SEMI_IMPLICIT)

           
  CALL wreal           ! Calculate True Vertical Velocity (W)               !

!  CALL REPORT("before VISCOF_H")

  !==============================================================================!
  !    TURBULENCE MODEL SECTION                                                  |
  !==============================================================================!

  
  SELECT CASE(vertical_mixing_type)
  CASE('closure')
     !=================General Ocean Turbulence Model==========================!
     !===================Original FVCOM MY-2.5/Galperin 1988 Model=============!

     CALL adv_q(q2,q2f)       !!Advection of Q2 

     CALL adv_q(q2l,q2lf) 



     IF(scalar_positivity_control) CALL fct_q2             !Conservation Correction   !
     IF(scalar_positivity_control) CALL fct_q2l            !Conservation Correction   !
!end !defined (ONE_D_MODEL)



     CALL vdif_q                  !! Solve Q2,Q2*L eqns for KH/KM/KQ 
      q2  = q2f
      q2l = q2lf


! end if defined(GOTM)
  CASE('constant')
     km = umol
     kh = umol*vprnu
  END SELECT

  
  CALL n2e3d(km,km1)




!==============================================================================!
!    SEDIMENT MODEL SECTION  
!      Advance Sed Model (Erode/Deposit/Advect/Diffuse)      
!      Change bathymetry (if MORPHO_MODEL=T)  
!        Note:  morph array returned from sed:  Deposition:  morph > 0
!        0 < morpho_factor < 1  
!      morpho_model and morpho_factor are stored and set in mod_sed.F              
!==============================================================================!
  

  !==============================================================================!
  !    UPDATE TEMPERATURE IN NON-BAROTROPIC CASE                                 !
  !==============================================================================!
  IF(temperature_active)THEN 
     
     CALL adv_t                                     !Advection                 !



     !#                                                   if !defined (DOUBLE_PRECISION)
     IF(scalar_positivity_control) CALL fct_t            !Conservation Correction   !
     !#                                                   endif
        !end !defined (ONE_D_MODEL)                                    !                          !

     IF(casename(1:3) == 'gom')THEN
        CALL vdif_ts_gom(1,tf1)
     ELSE  
        CALL vdif_ts(1,tf1)                            !Vertical Diffusion        !
     END IF

    
     CALL bcond_ts(1)                               !Boundary Conditions       !

     IF(nesting_on )THEN
       CALL set_var(inttime,t1=tf1)
     END IF

!!$!QXU{
!!$#  if defined (SST_GRID_ASSIM)
!!$      IF(ASSIM_FLAG==0 .AND. .NOT. SST_ASSIM_GRD)CALL TEMP_NUDGING
!!$      IF(ASSIM_FLAG==1 .AND. SST_ASSIM_GRD)CALL TEMP_NUDGING
!!$#  else
!!$      IF(ASSIM_FLAG==0 .AND. .NOT. SST_ASSIM)CALL TEMP_NUDGING
!!$      IF(ASSIM_FLAG==1 .AND. SST_ASSIM)CALL TEMP_NUDGING
!!$#  endif
!!$!QXU}
     




!     IF(NESTING_ON )THEN
!       CALL SET_VAR(intTime,T1=TF1)
!     END IF

!J. Ge for tracer advection
     IF(backward_advection==.true.)THEN
       IF(backward_step==1)THEN
        t0 = t1
       ELSEIF(backward_step==2)THEN
         t2 = t0
         t0 = t1
       ENDIF
     ENDIF
!J. Ge for tracer advection
     t1 = tf1                                       !Update to new time level  !

     CALL n2e3d(t1,t)                               !Shift to Elements         !


  END IF                                         !                          !

  !==============================================================================!
  !    UPDATE SALINITY IN NON-BAROTROPIC CASE                                    !
  !==============================================================================!

  IF(salinity_active)THEN                            !                          !   

     CALL adv_s                                     !Advection                 !

     

     !#                                                   if !defined (DOUBLE_PRECISION)
     IF(scalar_positivity_control) CALL fct_s       !Conservation Correction   !
     !#                                                   endif
     !end !defined (ONE_D_MODEL)                                    !                          !

     IF(casename(1:3) == 'gom')THEN
        CALL vdif_ts_gom(2,sf1)
     ELSE  
        CALL vdif_ts(2,sf1)                            !Vertical Diffusion        !
     END IF


     CALL bcond_ts(2)                               !Boundary Conditions       !

     IF(nesting_on )THEN
       CALL set_var(inttime,s1=sf1)
     END IF




!     IF(NESTING_ON )THEN
!       CALL SET_VAR(intTime,S1=SF1)
!     END IF

!J. Ge for tracer advection
     IF(backward_advection==.true.)THEN
       IF(backward_step==1)THEN
        s0 = s1
       ELSEIF(backward_step==2)THEN
         s2 = s0
         s0 = s1
       ENDIF
     ENDIF
!J. Ge for tracer advection
     s1 = sf1                                     !Update to new time level  !

     CALL n2e3d(s1,s)                               !Shift to Elements         !


  END IF                      

  !==============================================================================!
  !endif defined (DYE)

!==================================================================================!
!    ADJUST TEMPERATURE AND SALINITY AT RIVER MOUTHS
!==================================================================================!
  IF( river_ts_setting == 'calculated')THEN
     CALL adjust_ts
  END IF
  




  
  
  !==============================================================================!
  !     UPDATE THE DENSITY IN NON-BAROTROPIC CASE                                |
  !==============================================================================!
  IF(.NOT.barotropic)THEN
   SELECT CASE(sea_water_density_function)
   CASE(sw_dens1)
      CALL dens1
   CASE(sw_dens2)
      CALL dens2
   CASE(sw_dens3)
      CALL dens3
   CASE DEFAULT
      CALL fatal_error("INVALID DENSITY FUNCTION SELECTED:",&
           & "   "//trim(sea_water_density_function) )
   END SELECT
  END IF
  !==============================================================================!
  !     MIMIC CONVECTIVE OVERTURNING TO STABILIZE VERTICAL DENSITY PROFILE       |
  !==============================================================================!
  
  IF(convective_overturning)THEN
     CALL conv_over
     IF(.NOT.barotropic)THEN
        SELECT CASE(sea_water_density_function)
        CASE(sw_dens1)
           CALL dens1
        CASE(sw_dens2)
           CALL dens2
        CASE(sw_dens3)
           CALL dens3
        CASE DEFAULT
           CALL fatal_error("INVALID DENSITY FUNCTION SELECTED:",&
                & "   "//trim(sea_water_density_function) )
        END SELECT
        
     END IF
  END IF
!==============================================================================!
!==============================================================================!
  ! end if !defined (TWO_D_MODEL)
!==============================================================================!
!==============================================================================!
  
  
  
  
  !==============================================================================!
  !     UPDATE VELOCITY FIELD (NEEDED TO WAIT FOR SALINITY/TEMP/TURB/TRACER)     |
  !==============================================================================!
  u = uf
  v = vf
  !==============================================================================!
  !    PERFORM DATA EXCHANGE FOR ELEMENT BASED INFORMATION AT PROC BNDRIES       |
  !==============================================================================!
  
  !==============================================================================!
  !     PERFORM DATA EXCHANGE FOR WATER QUALITY VARIABLES                        |
  !==============================================================================!
  ! end if !defined (TWO_D_MODEL)
  


  !
  !----SHIFT SEA SURFACE ELEVATION AND DEPTH TO CURRENT TIME LEVEL---------------!
  !
  et  = el  
  dt  = d 
  et1 = el1
  dt1 = d1
  
  IF(wetting_drying_on) CALL wd_update(3)
  
  ! New Open Boundary Condition ----10


!!$  CALL DUMP_PROBE_DATA 


  if(dbg_set(dbg_sbr)) write(ipt,*)&
       & "End: internal_step"

END SUBROUTINE internal_step