New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 8093 for branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfphy.F90 – NEMO

Ignore:
Timestamp:
2017-05-30T10:13:14+02:00 (7 years ago)
Author:
gm
Message:

#1880 (HPC-09) - step-6: prepare some forthcoming evolutions (ZDF modules mainly)

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfphy.F90

    r7990 r8093  
    1111   !!   zdf_phy       : upadate at each time-step the vertical mixing coeff.  
    1212   !!---------------------------------------------------------------------- 
    13    USE par_oce        ! ocean parameters 
     13   USE oce            ! ocean dynamics and tracers variables 
    1414   USE zdf_oce        ! vertical physics: shared variables          
     15   USE zdfdrg         ! vertical physics: top/bottom drag coef. 
     16!!gm old 
    1517   USE zdfbfr         ! vertical physics: bottom friction 
    16    USE zdfric         ! vertical physics: Richardson vertical mixing    
     18!!gm 
     19   USE zdfsh2         ! vertical physics: shear production term of TKE 
     20   USE zdfric         ! vertical physics: RIChardson dependent vertical mixing    
    1721   USE zdftke         ! vertical physics: TKE vertical mixing 
    1822   USE zdfgls         ! vertical physics: GLS vertical mixing 
     
    4448   INTEGER, PARAMETER ::   np_TKE = 3   ! Turbulente Kinetic Eenergy closure scheme for Kz 
    4549   INTEGER, PARAMETER ::   np_GLS = 4   ! Generic Length Scale closure scheme for Kz 
    46        
     50 
     51   LOGICAL ::   l_zdfsh2   ! shear production term flag (=F for CST, =T otherwise (i.e. TKE, GLS, RIC)) 
     52 
    4753   !!---------------------------------------------------------------------- 
    4854   !! NEMO/OPA 4.0 , NEMO Consortium (2017) 
     
    6167      !!                set horizontal shape and vertical profile of background mixing coef. 
    6268      !!---------------------------------------------------------------------- 
     69      INTEGER ::   jk            ! dummy loop indices 
    6370      INTEGER ::   ioptio, ios   ! local integers 
    6471      !! 
     
    9097         WRITE(numout,*) '      vertical closure scheme' 
    9198         WRITE(numout,*) '         constant vertical mixing coefficient    ln_zdfcst = ', ln_zdfcst 
    92          WRITE(numout,*) '         constant vertical mixing coefficient    ln_zdfric = ', ln_zdfric 
    93          WRITE(numout,*) '         constant vertical mixing coefficient    ln_zdftke = ', ln_zdftke 
    94          WRITE(numout,*) '         constant vertical mixing coefficient    ln_zdfgls = ', ln_zdfgls 
     99         WRITE(numout,*) '         Richardson number dependent closure     ln_zdfric = ', ln_zdfric 
     100         WRITE(numout,*) '         Turbulent Kinetic Energy closure (TKE)  ln_zdftke = ', ln_zdftke 
     101         WRITE(numout,*) '         Generic Length Scale closure (GLS)      ln_zdfgls = ', ln_zdfgls 
    95102         WRITE(numout,*) '      convection: ' 
    96103         WRITE(numout,*) '         enhanced vertical diffusion             ln_zdfevd = ', ln_zdfevd 
     
    137144      ENDIF 
    138145      ! 
    139       !                                   ! set 1st/last level Av to zero one for all 
    140       avt(:,:,1) = 0._wp   ;   avs(:,:,1) = 0._wp   ;   avm(:,:,1) = 0._wp 
     146      DO jk = 1, jpk                      ! set turbulent closure Kz to the background value (avt_k, avm_k) 
     147         avt_k(:,:,jk) = avtb_2d(:,:) * avtb(jk) * wmask (:,:,jk) 
     148         avm_k(:,:,jk) =                avmb(jk) * wmask (:,:,jk) 
     149      END DO 
     150!!gm  to be tested only the 1st & last levels 
     151!      avt  (:,:, 1 ) = 0._wp   ;   avs(:,:, 1 ) = 0._wp   ;   avm  (:,:, 1 ) = 0._wp 
     152!      avt  (:,:,jpk) = 0._wp   ;   avs(:,:,jpk) = 0._wp   ;   avm  (:,:,jpk) = 0._wp 
     153!!gm 
     154      avt  (:,:,:) = 0._wp   ;   avs(:,:,:) = 0._wp   ;   avm  (:,:,:) = 0._wp 
    141155 
    142156      !                          !==  Convection  ==! 
     
    170184      IF( ln_zdfric .OR. ln_zdfgls )    CALL ctl_stop( 'zdf_phy_init: zdfric and zdfgls never tested with ice shelves cavities ' ) 
    171185      ENDIF 
     186      !                                ! shear production term flag 
     187      IF( ln_zdfcst ) THEN   ;   l_zdfsh2 = .FALSE. 
     188      ELSE                   ;   l_zdfsh2 = .TRUE. 
     189      ENDIF 
    172190 
    173191      !                          !== gravity wave-driven mixing  ==! 
     
    175193      IF( ln_zdfswm )   CALL zdf_swm_init       ! surface  wave-driven mixing 
    176194 
    177       !                          !== bottom friction  ==! 
     195      !                          !== top/bottom friction  ==! 
     196      CALL zdf_drg_init 
     197!!gm old 
    178198      CALL zdf_bfr_init 
     199!!gm 
    179200      ! 
    180201      !                          !== time-stepping  ==! 
     
    200221      ! 
    201222      INTEGER ::   ji, jj, jk   ! dummy loop indice 
     223      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zsh2   ! shear production 
    202224      !! --------------------------------------------------------------------- 
    203225      ! 
     226!!gm old 
    204227      CALL zdf_bfr( kt )                        !* bottom friction (if quadratic) 
    205       !                        
     228!!gm 
     229      ! 
     230!      IF( l_zdfdrg ) THEN     !==  update top/bottom drag  ==!   (non-linear cases) 
     231!         ! 
     232!         !                       !* bottom drag 
     233!         CALL zdf_drg( kt, mbkt    , r_Cdmin_bot, r_Cdmax_bot,   &   ! <<== in  
     234!            &              r_z0_bot,   r_ke0_bot,    rCd0_bot,   & 
     235!            &                                        rCdU_bot  )     ! ==>> out : bottom drag [m/s] 
     236!         IF( ln_isfcav ) THEN    !* top drag   (ocean cavities) 
     237!            CALL zdf_drg( kt, mikt    , r_Cdmin_top, r_Cdmax_top,   &   ! <<== in  
     238!               &              r_z0_top,   r_ke0_top,    rCd0_top,   & 
     239!               &                                        rCdU_top  )     ! ==>> out : bottom drag [m/s] 
     240!         ENDIF 
     241!      ENDIF 
     242      ! 
     243      !                       !==  Kz from chosen turbulent closure  ==!   (avm_k, avt_k) 
     244      ! 
     245      IF( l_zdfsh2 )   &         !* shear production at w-points (energy conserving form) 
     246         CALL zdf_sh2( ub, vb, un, vn, avm_k,   &     ! <<== in 
     247            &                           zsh2    )     ! ==>> out : shear production 
     248      ! 
    206249      SELECT CASE ( nzdf_phy )                  !* Vertical eddy viscosity and diffusivity coefficients at w-points 
    207       CASE( np_RIC )   ;   CALL zdf_ric( kt )         ! Richardson number dependent Kz 
    208       CASE( np_TKE )   ;   CALL zdf_tke( kt )         ! TKE closure scheme for Kz 
    209       CASE( np_GLS )   ;   CALL zdf_gls( kt )         ! GLS closure scheme for Kz 
    210       CASE( np_CST )                                  ! Constant Kz (reset avt, avm to the background value) 
    211          avt(2:jpim1,2:jpjm1,1:jpkm1) = rn_avt0 * wmask(2:jpim1,2:jpjm1,1:jpkm1) 
    212          avm(2:jpim1,2:jpjm1,1:jpkm1) = rn_avm0 * wmask(2:jpim1,2:jpjm1,1:jpkm1) 
     250      CASE( np_RIC )   ;   CALL zdf_ric( kt, gdept_n, zsh2, avm_k, avt_k )    ! Richardson number dependent Kz 
     251      CASE( np_TKE )   ;   CALL zdf_tke( kt         , zsh2, avm_k, avt_k )    ! TKE closure scheme for Kz 
     252      CASE( np_GLS )   ;   CALL zdf_gls( kt         , zsh2, avm_k, avt_k )    ! GLS closure scheme for Kz 
     253!     CASE( np_CST )                                  ! Constant Kz (reset avt, avm to the background value) 
     254!         ! avt_k and avm_k set one for all at initialisation phase 
     255!!gm         avt(2:jpim1,2:jpjm1,1:jpkm1) = rn_avt0 * wmask(2:jpim1,2:jpjm1,1:jpkm1) 
     256!!gm         avm(2:jpim1,2:jpjm1,1:jpkm1) = rn_avm0 * wmask(2:jpim1,2:jpjm1,1:jpkm1) 
    213257      END SELECT 
     258      !   
     259      !                          !==  ocean Kz  ==!   (avt, avs, avm) 
     260      ! 
     261      !                                         !* start from turbulent closure values 
     262      avt(:,:,2:jpkm1) = avt_k(:,:,2:jpkm1) 
     263      avm(:,:,2:jpkm1) = avm_k(:,:,2:jpkm1) 
    214264      ! 
    215265      IF( ln_rnf_mouth ) THEN                   !* increase diffusivity at rivers mouths 
     
    219269      ENDIF 
    220270      ! 
    221       IF( ln_zdfevd )   CALL zdf_evd( kt )      !* convection: enhanced vertical eddy diffusivity 
     271      IF( ln_zdfevd )   CALL zdf_evd( kt, avm, avt )  !* convection: enhanced vertical eddy diffusivity 
    222272      ! 
    223273      !                                         !* double diffusive mixing 
    224274      IF( ln_zdfddm ) THEN                            ! update avt and compute avs 
    225                         CALL zdf_ddm( kt ) 
     275                        CALL zdf_ddm( kt, avm, avt, avs ) 
    226276      ELSE                                            ! same mixing on all tracers 
    227277         avs(2:jpim1,2:jpjm1,1:jpkm1) = avt(2:jpim1,2:jpjm1,1:jpkm1) 
     
    229279      ! 
    230280      !                                         !* wave-induced mixing  
    231       IF( ln_zdfswm )   CALL zdf_swm( kt )            ! surface  wave (Qiao et al. 2004)  
    232       IF( ln_zdfiwm )   CALL zdf_iwm( kt )            ! internal wave (de Lavergne et al 2017) 
     281      IF( ln_zdfswm )   CALL zdf_swm( kt, avm, avt, avs )   ! surface  wave (Qiao et al. 2004)  
     282      IF( ln_zdfiwm )   CALL zdf_iwm( kt, avm, avt, avs )   ! internal wave (de Lavergne et al 2017) 
    233283 
    234284 
    235285      !                                         !* Lateral boundary conditions (sign unchanged) 
    236       CALL lbc_lnk( avm, 'W', 1. ) 
    237       CALL lbc_lnk( avt, 'W', 1. ) 
    238       CALL lbc_lnk( avs, 'W', 1. ) 
    239       ! 
    240 !!gm  ===>>>   TO BE REMOVED  
    241       DO jk = 1, jpkm1            !* vertical eddy viscosity at wu- and wv-points 
    242          DO jj = 2, jpjm1 
    243             DO ji = 2, jpim1 
    244                avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj  ,jk) ) * wumask(ji,jj,jk) 
    245                avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji  ,jj+1,jk) ) * wvmask(ji,jj,jk) 
    246             END DO 
    247          END DO 
    248       END DO 
    249       CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! Lateral boundary conditions 
    250 !!gm end 
     286      CALL lbc_lnk( avm_k, 'W', 1. )                  ! needed to compute the shear production term 
     287      CALL lbc_lnk( avt_k, 'W', 1. )                  !!gm a priori useless ==>> to be tested 
     288      CALL lbc_lnk( avm  , 'W', 1. )                  ! needed to compute avm at u- and v-points 
     289      CALL lbc_lnk( avt  , 'W', 1. )                  !!gm  a priori only avm_k and avm are required 
     290      CALL lbc_lnk( avs  , 'W', 1. )                  !!gm  To be tested 
     291      ! 
    251292 
    252293 
    253294      CALL zdf_mxl( kt )                        !* mixed layer depth, and level 
    254295 
    255 !!gm  !==>> to be moved in zdftke & zdfgls     
    256                                                       ! write TKE or GLS information in the restart file 
    257       IF( lrst_oce .AND. ln_zdftke )   CALL tke_rst( kt, 'WRITE' ) 
    258       IF( lrst_oce .AND. ln_zdfgls )   CALL gls_rst( kt, 'WRITE' ) 
    259       !       
     296      IF( lrst_oce ) THEN                       !* write TKE, GLS or RIC fields in the restart file 
     297         IF( ln_zdftke )   CALL tke_rst( kt, 'WRITE' ) 
     298         IF( ln_zdfgls )   CALL gls_rst( kt, 'WRITE' ) 
     299         IF( ln_zdfric )   CALL ric_rst( kt, 'WRITE' )  
     300      ENDIF 
     301      ! 
    260302   END SUBROUTINE zdf_phy 
    261303 
Note: See TracChangeset for help on using the changeset viewer.