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 9939 for NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/ZDF – NEMO

Ignore:
Timestamp:
2018-07-13T09:28:50+02:00 (6 years ago)
Author:
gm
Message:

#1911 (ENHANCE-04): RK3 branche phased with MLF@9937 branche

Location:
NEMO/branches/2018/dev_r9838_ENHANCE04_RK3
Files:
8 edited
1 copied

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/ZDF/zdfddm.F90

    r9598 r9939  
    77   !!   NEMO     1.0  ! 2002-06  (G. Madec)  F90: Free form and module 
    88   !!            3.3  ! 2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase 
    9    !!            3.6  ! 2013-04  (G. Madec, F. Roquet) zrau compute locally using interpolation of alpha & beta 
     9   !!            3.6  ! 2013-04  (G. Madec, F. Roquet) zrho compute locally using interpolation of alpha & beta 
    1010   !!            4.0  !  2017-04  (G. Madec)  remove CPP ddm key & avm at t-point only  
    1111   !!---------------------------------------------------------------------- 
    … …  
    7979      REAL(wp) ::   zavft, zavfs    !   -      - 
    8080      REAL(wp) ::   zavdt, zavds    !   -      - 
    81       REAL(wp), DIMENSION(jpi,jpj) ::   zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3 
     81      REAL(wp), DIMENSION(jpi,jpj) ::   zrho, zmsks, zmskf, zmskd1, zmskd2, zmskd3 
    8282      !!---------------------------------------------------------------------- 
    8383      ! 
    … …  
    9191!!gm                            and many acces in memory 
    9292          
    93          DO jj = 1, jpj                !==  R=zrau = (alpha / beta) (dk[t] / dk[s])  ==! 
     93         DO jj = 1, jpj                !==  R=zrho = (alpha / beta) (dk[t] / dk[s])  ==! 
    9494            DO ji = 1, jpi 
    9595               zrw =   ( gdepw_n(ji,jj,jk  ) - gdept_n(ji,jj,jk) )   & 
    … …  
    105105               zds = zbw * ( tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal) )  
    106106               IF( ABS( zds) <= 1.e-20_wp )   zds = 1.e-20_wp 
    107                zrau(ji,jj) = MAX(  1.e-20, zdt / zds  )    ! only retains positive value of zrau 
     107               zrho(ji,jj) = MAX(  1.e-20, zdt / zds  )    ! only retains positive value of zrho 
    108108            END DO 
    109109         END DO 
    … …  
    116116               ENDIF 
    117117               ! salt fingering indicator: msksf=1 if R>1; 0 elsewhere             
    118                IF( zrau(ji,jj) <= 1.             ) THEN   ;   zmskf(ji,jj) = 0._wp 
     118               IF( zrho(ji,jj) <= 1.             ) THEN   ;   zmskf(ji,jj) = 0._wp 
    119119               ELSE                                       ;   zmskf(ji,jj) = 1._wp 
    120120               ENDIF 
    121121               ! diffusive layering indicators:  
    122122               !     ! mskdl1=1 if 0< R <1; 0 elsewhere 
    123                IF( zrau(ji,jj) >= 1.             ) THEN   ;   zmskd1(ji,jj) = 0._wp 
     123               IF( zrho(ji,jj) >= 1.             ) THEN   ;   zmskd1(ji,jj) = 0._wp 
    124124               ELSE                                       ;   zmskd1(ji,jj) = 1._wp 
    125125               ENDIF 
    126126               !     ! mskdl2=1 if 0< R <0.5; 0 elsewhere 
    127                IF( zrau(ji,jj) >= 0.5            ) THEN   ;   zmskd2(ji,jj) = 0._wp 
     127               IF( zrho(ji,jj) >= 0.5            ) THEN   ;   zmskd2(ji,jj) = 0._wp 
    128128               ELSE                                       ;   zmskd2(ji,jj) = 1._wp 
    129129               ENDIF 
    130130               !   mskdl3=1 if 0.5< R <1; 0 elsewhere 
    131                IF( zrau(ji,jj) <= 0.5 .OR. zrau(ji,jj) >= 1. ) THEN   ;   zmskd3(ji,jj) = 0._wp 
     131               IF( zrho(ji,jj) <= 0.5 .OR. zrho(ji,jj) >= 1. ) THEN   ;   zmskd3(ji,jj) = 0._wp 
    132132               ELSE                                                   ;   zmskd3(ji,jj) = 1._wp 
    133133               ENDIF 
    … …  
    143143         DO jj = 1, jpj 
    144144            DO ji = 1, jpi 
    145                zinr = 1._wp / zrau(ji,jj) 
     145               zinr = 1._wp / zrho(ji,jj) 
    146146               ! salt fingering 
    147                zrr = zrau(ji,jj) / rn_hsbfr 
     147               zrr = zrho(ji,jj) / rn_hsbfr 
    148148               zrr = zrr * zrr 
    149149               zavfs = rn_avts / ( 1 + zrr*zrr*zrr ) * zmsks(ji,jj) * zmskf(ji,jj) 
    … …  
    151151               ! diffusive layering 
    152152               zavdt = 1.3635e-6 * EXP(  4.6 * EXP( -0.54*(zinr-1.) )  ) * zmsks(ji,jj) * zmskd1(ji,jj) 
    153                zavds = zavdt * zmsks(ji,jj) * (  ( 1.85 * zrau(ji,jj) - 0.85 ) * zmskd3(ji,jj)   & 
    154                   &                             +  0.15 * zrau(ji,jj)          * zmskd2(ji,jj)  ) 
     153               zavds = zavdt * zmsks(ji,jj) * (  ( 1.85 * zrho(ji,jj) - 0.85 ) * zmskd3(ji,jj)   & 
     154                  &                             +  0.15 * zrho(ji,jj)          * zmskd2(ji,jj)  ) 
    155155               ! add to the eddy viscosity coef. previously computed 
    156156               p_avs(ji,jj,jk) = p_avt(ji,jj,jk) + zavfs + zavds 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/ZDF/zdfdrg.F90

    r9598 r9939  
    162162      INTEGER  ::   ji, jj       ! dummy loop indexes 
    163163      INTEGER  ::   ikbu, ikbv   ! local integers 
    164       REAL(wp) ::   zm1_2dt      ! local scalar 
    165164      REAL(wp) ::   zCdu, zCdv   !   -      - 
    166165      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdu, ztrdv 
    167166      !!--------------------------------------------------------------------- 
    168167      ! 
    169 !!gm bug : time step is only rdt (not 2 rdt if euler start !) 
    170       zm1_2dt = - 1._wp / ( 2._wp * rdt ) 
    171  
    172168      IF( l_trddyn ) THEN      ! trends: store the input trends 
    173169         ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) ) 
    … …  
    185181            zCdv = 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / e3v_n(ji,jj,ikbv) 
    186182            ! 
    187             pua(ji,jj,ikbu) = pua(ji,jj,ikbu) + MAX(  zCdu , zm1_2dt  ) * pub(ji,jj,ikbu) 
    188             pva(ji,jj,ikbv) = pva(ji,jj,ikbv) + MAX(  zCdv , zm1_2dt  ) * pvb(ji,jj,ikbv) 
     183            pua(ji,jj,ikbu) = pua(ji,jj,ikbu) + MAX(  zCdu , - r1_Dt  ) * pub(ji,jj,ikbu) 
     184            pva(ji,jj,ikbv) = pva(ji,jj,ikbv) + MAX(  zCdv , - r1_Dt  ) * pvb(ji,jj,ikbv) 
    189185         END DO 
    190186      END DO 
    … …  
    200196               zCdv = 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) / e3v_n(ji,jj,ikbv) 
    201197               ! 
    202                pua(ji,jj,ikbu) = pua(ji,jj,ikbu) + MAX(  zCdu , zm1_2dt  ) * pub(ji,jj,ikbu) 
    203                pva(ji,jj,ikbv) = pva(ji,jj,ikbv) + MAX(  zCdv , zm1_2dt  ) * pvb(ji,jj,ikbv) 
     198               pua(ji,jj,ikbu) = pua(ji,jj,ikbu) + MAX(  zCdu , - r1_Dt  ) * pub(ji,jj,ikbu) 
     199               pva(ji,jj,ikbv) = pva(ji,jj,ikbv) + MAX(  zCdv , - r1_Dt  ) * pvb(ji,jj,ikbv) 
    204200           END DO 
    205201         END DO 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/ZDF/zdfgls.F90

    r9598 r9939  
    170170            ! 
    171171            ! surface friction 
    172             ustar2_surf(ji,jj) = r1_rau0 * taum(ji,jj) * tmask(ji,jj,1) 
     172            ustar2_surf(ji,jj) = r1_rho0 * taum(ji,jj) * tmask(ji,jj,1) 
    173173            !    
    174174!!gm Rq we may add here r_ke0(_top/_bot) ?  ==>> think about that... 
    … …  
    280280               zd_up(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) ) / ( e3t_n(ji,jj,jk  ) * e3w_n(ji,jj,jk) ) 
    281281               !                                               ! diagonal 
    282                zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk)  + rdt * zdiss * wmask(ji,jj,jk)  
     282               zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk)  + rn_Dt * zdiss * wmask(ji,jj,jk)  
    283283               !                                               ! right hand side in en 
    284                en(ji,jj,jk) = en(ji,jj,jk) + rdt * zesh2 * wmask(ji,jj,jk) 
     284               en(ji,jj,jk) = en(ji,jj,jk) + rn_Dt * zesh2 * wmask(ji,jj,jk) 
    285285            END DO 
    286286         END DO 
    … …  
    530530               zd_up(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) ) / ( e3t_n(ji,jj,jk  ) * e3w_n(ji,jj,jk) ) 
    531531               !                                               ! diagonal 
    532                zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) + rdt * zdiss * wmask(ji,jj,jk) 
     532               zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) + rn_Dt * zdiss * wmask(ji,jj,jk) 
    533533               !                                               ! right hand side in psi 
    534                psi(ji,jj,jk) = psi(ji,jj,jk) + rdt * zesh2 * wmask(ji,jj,jk) 
     534               psi(ji,jj,jk) = psi(ji,jj,jk) + rn_Dt * zesh2 * wmask(ji,jj,jk) 
    535535            END DO 
    536536         END DO 
    … …  
    11051105      rc04  = rc03 * rc0 
    11061106      rsbc_tke1 = -3._wp/2._wp*rn_crban*ra_sf*rl_sf                      ! Dirichlet + Wave breaking 
    1107       rsbc_tke2 = rdt * rn_crban / rl_sf                                 ! Neumann + Wave breaking  
     1107      rsbc_tke2 = rn_Dt * rn_crban / rl_sf                               ! Neumann + Wave breaking  
    11081108      zcr = MAX(rsmall, rsbc_tke1**(1./(-ra_sf*3._wp/2._wp))-1._wp ) 
    11091109      rtrans = 0.2_wp / zcr                                              ! Ad. inverse transition length between log and wave layer  
    11101110      rsbc_zs1  = rn_charn/grav                                          ! Charnock formula for surface roughness 
    11111111      rsbc_zs2  = rn_frac_hs / 0.85_wp / grav * 665._wp                  ! Rascle formula for surface roughness  
    1112       rsbc_psi1 = -0.5_wp * rdt * rc0**(rpp-2._wp*rmm) / rsc_psi 
    1113       rsbc_psi2 = -0.5_wp * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking  
    1114       ! 
    1115       rfact_tke = -0.5_wp / rsc_tke * rdt                                ! Cst used for the Diffusion term of tke 
    1116       rfact_psi = -0.5_wp / rsc_psi * rdt                                ! Cst used for the Diffusion term of tke 
     1112      rsbc_psi1 = -0.5_wp * rn_Dt * rc0**(rpp-2._wp*rmm) / rsc_psi 
     1113      rsbc_psi2 = -0.5_wp * rn_Dt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking  
     1114      ! 
     1115      rfact_tke = -0.5_wp / rsc_tke * rn_Dt                              ! Cst used for the Diffusion term of tke 
     1116      rfact_psi = -0.5_wp / rsc_psi * rn_Dt                              ! Cst used for the Diffusion term of tke 
    11171117      ! 
    11181118      !                                !* Wall proximity function 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/ZDF/zdfiwm.F90

    r9598 r9939  
    8787      !!              This is divided into three components: 
    8888      !!                 1. Bottom-intensified low-mode dissipation at critical slopes 
    89       !!                     zemx_iwm(z) = ( ecri_iwm / rau0 ) * EXP( -(H-z)/hcri_iwm ) 
     89      !!                     zemx_iwm(z) = ( ecri_iwm / rho0 ) * EXP( -(H-z)/hcri_iwm ) 
    9090      !!                                   / ( 1. - EXP( - H/hcri_iwm ) ) * hcri_iwm 
    9191      !!              where hcri_iwm is the characteristic length scale of the bottom  
    9292      !!              intensification, ecri_iwm a map of available power, and H the ocean depth. 
    9393      !!                 2. Pycnocline-intensified low-mode dissipation 
    94       !!                     zemx_iwm(z) = ( epyc_iwm / rau0 ) * ( sqrt(rn2(z))^nn_zpyc ) 
     94      !!                     zemx_iwm(z) = ( epyc_iwm / rho0 ) * ( sqrt(rn2(z))^nn_zpyc ) 
    9595      !!                                   / SUM( sqrt(rn2(z))^nn_zpyc * e3w(z) ) 
    9696      !!              where epyc_iwm is a map of available power, and nn_zpyc 
    … …  
    9898      !!              energy dissipation. 
    9999      !!                 3. WKB-height dependent high mode dissipation 
    100       !!                     zemx_iwm(z) = ( ebot_iwm / rau0 ) * rn2(z) * EXP(-z_wkb(z)/hbot_iwm) 
     100      !!                     zemx_iwm(z) = ( ebot_iwm / rho0 ) * rn2(z) * EXP(-z_wkb(z)/hbot_iwm) 
    101101      !!                                   / SUM( rn2(z) * EXP(-z_wkb(z)/hbot_iwm) * e3w(z) ) 
    102102      !!              where hbot_iwm is the characteristic length scale of the WKB bottom  
    … …  
    151151         DO ji = 1, jpi 
    152152            zhdep(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1)       ! depth of the ocean 
    153             zfact(ji,jj) = rau0 * (  1._wp - EXP( -zhdep(ji,jj) / hcri_iwm(ji,jj) )  ) 
     153            zfact(ji,jj) = rho0 * (  1._wp - EXP( -zhdep(ji,jj) / hcri_iwm(ji,jj) )  ) 
    154154            IF( zfact(ji,jj) /= 0._wp )   zfact(ji,jj) = ecri_iwm(ji,jj) / zfact(ji,jj) 
    155155         END DO 
    … …  
    180180         DO jj = 1, jpj 
    181181            DO ji = 1, jpi 
    182                IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rau0 * zfact(ji,jj) ) 
     182               IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rho0 * zfact(ji,jj) ) 
    183183            END DO 
    184184         END DO 
    … …  
    197197         DO jj= 1, jpj 
    198198            DO ji = 1, jpi 
    199                IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rau0 * zfact(ji,jj) ) 
     199               IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rho0 * zfact(ji,jj) ) 
    200200            END DO 
    201201         END DO 
    … …  
    247247      DO jj = 1, jpj 
    248248         DO ji = 1, jpi 
    249             IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = ebot_iwm(ji,jj) / ( rau0 * zfact(ji,jj) ) 
     249            IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = ebot_iwm(ji,jj) / ( rho0 * zfact(ji,jj) ) 
    250250         END DO 
    251251      END DO 
    … …  
    260260      ! Calculate molecular kinematic viscosity 
    261261      znu_t(:,:,:) = 1.e-4_wp * (  17.91_wp - 0.53810_wp * tsn(:,:,:,jp_tem) + 0.00694_wp * tsn(:,:,:,jp_tem) * tsn(:,:,:,jp_tem)  & 
    262          &                                  + 0.02305_wp * tsn(:,:,:,jp_sal)  ) * tmask(:,:,:) * r1_rau0 
     262         &                                  + 0.02305_wp * tsn(:,:,:,jp_sal)  ) * tmask(:,:,:) * r1_rho0 
    263263      DO jk = 2, jpkm1 
    264264         znu_w(:,:,jk) = 0.5_wp * ( znu_t(:,:,jk-1) + znu_t(:,:,jk) ) * wmask(:,:,jk) 
    … …  
    306306         END DO 
    307307         IF( lk_mpp )   CALL mpp_sum( zztmp ) 
    308          zztmp = rau0 * zztmp ! Global integral of rauo * Kz * N^2 = power contributing to mixing  
     308         zztmp = rho0 * zztmp ! Global integral of rhoo * Kz * N^2 = power contributing to mixing  
    309309         ! 
    310310         IF(lwp) THEN 
    … …  
    350350                                    !* output useful diagnostics: Kz*N^2 ,  
    351351!!gm Kz*N2 should take into account the ratio avs/avt if it is used.... (see diaar5) 
    352                                     !  vertical integral of rau0 * Kz * N^2 , energy density (zemx_iwm) 
     352                                    !  vertical integral of rho0 * Kz * N^2 , energy density (zemx_iwm) 
    353353      IF( iom_use("bflx_iwm") .OR. iom_use("pcmap_iwm") ) THEN 
    354354         ALLOCATE( z2d(jpi,jpj) , z3d(jpi,jpj,jpk) ) 
    … …  
    358358            z2d(:,:) = z2d(:,:) + e3w_n(:,:,jk) * z3d(:,:,jk) * wmask(:,:,jk) 
    359359         END DO 
    360          z2d(:,:) = rau0 * z2d(:,:) 
     360         z2d(:,:) = rho0 * z2d(:,:) 
    361361         CALL iom_put( "bflx_iwm", z3d ) 
    362362         CALL iom_put( "pcmap_iwm", z2d ) 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/ZDF/zdfmxl.F90

    r9598 r9939  
    9393      nmln(:,:)  = nlb10               ! Initialization to the number of w ocean point 
    9494      hmlp(:,:)  = 0._wp               ! here hmlp used as a dummy variable, integrating vertically N^2 
    95       zN2_c = grav * rho_c * r1_rau0   ! convert density criteria into N^2 criteria 
     95      zN2_c = grav * rho_c * r1_rho0   ! convert density criteria into N^2 criteria 
    9696      DO jk = nlb10, jpkm1 
    9797         DO jj = 1, jpj                ! Mixed layer level: w-level  

  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/ZDF/zdfosm.F90

    r9598 r9939  
    298298        DO ji = 2, jpim1 
    299299           ! Surface downward irradiance (so always +ve) 
    300            zrad0(ji,jj) = qsr(ji,jj) * r1_rau0_rcp 
     300           zrad0(ji,jj) = qsr(ji,jj) * r1_rho0_rcp 
    301301           ! Downwards irradiance at base of boundary layer 
    302302           zradh(ji,jj) = zrad0(ji,jj) * ( zz0 * EXP( -hbl(ji,jj)/rn_si0 ) + zz1 * EXP( -hbl(ji,jj)/rn_si1) ) 
    … …  
    312312           zbeta    = rab_n(ji,jj,1,jp_sal) 
    313313           ! Upwards surface Temperature flux for non-local term 
    314            zwth0(ji,jj) = - qns(ji,jj) * r1_rau0_rcp * tmask(ji,jj,1) 
     314           zwth0(ji,jj) = - qns(ji,jj) * r1_rho0_rcp * tmask(ji,jj,1) 
    315315           ! Upwards surface salinity flux for non-local term 
    316            zws0(ji,jj) = - ( ( emp(ji,jj)-rnf(ji,jj) ) * tsn(ji,jj,1,jp_sal)  + sfx(ji,jj) ) * r1_rau0 * tmask(ji,jj,1) 
     316           zws0(ji,jj) = - ( ( emp(ji,jj)-rnf(ji,jj) ) * tsn(ji,jj,1,jp_sal)  + sfx(ji,jj) ) * r1_rho0 * tmask(ji,jj,1) 
    317317           ! Non radiative upwards surface buoyancy flux 
    318318           zwb0(ji,jj) = grav * zthermal * zwth0(ji,jj) -  grav * zbeta * zws0(ji,jj) 
    … …  
    324324           zwbav(ji,jj) = grav  * zthermal * zwthav(ji,jj) - grav  * zbeta * zwsav(ji,jj) 
    325325           ! Surface upward velocity fluxes 
    326            zuw0(ji,jj) = -utau(ji,jj) * r1_rau0 * tmask(ji,jj,1) 
    327            zvw0(ji,jj) = -vtau(ji,jj) * r1_rau0 * tmask(ji,jj,1) 
     326           zuw0(ji,jj) = -utau(ji,jj) * r1_rho0 * tmask(ji,jj,1) 
     327           zvw0(ji,jj) = -vtau(ji,jj) * r1_rho0 * tmask(ji,jj,1) 
    328328           ! Friction velocity (zustar), at T-point : LMD94 eq. 2 
    329329           zustar(ji,jj) = MAX( SQRT( SQRT( zuw0(ji,jj) * zuw0(ji,jj) + zvw0(ji,jj) * zvw0(ji,jj) ) ), 1.0e-8 ) 
    … …  
    455455                           &            + 0.135 * zla(ji,jj) * zwstrl(ji,jj)**3/hbl(ji,jj) ) 
    456456 
    457                       zvel_max =  - ( 1.0 + 1.0 * ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / hbl(ji,jj) ) & 
     457                      zvel_max =  - ( 1.0 + 1.0 * ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_Dt / hbl(ji,jj) ) & 
    458458                           &   * zwb_ent(ji,jj) / ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
    459459! Entrainment including component due to shear turbulence. Modified Langmuir component, but gives same result for La=0.3 For testing uncomment. 
    … …  
    461461!                           &            + ( 0.15 * ( 1.0 - EXP( -0.5 * zla(ji,jj) ) ) + 0.03 / zla(ji,jj)**2 ) * zustar(ji,jj)**3/hbl(ji,jj) ) 
    462462 
    463 !                      zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / zhbl(ji,jj) ) * zwb_ent(ji,jj) / & 
     463!                      zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_Dt / zhbl(ji,jj) ) * zwb_ent(ji,jj) / & 
    464464!                           &       ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
    465465                      zzdhdt = - zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj),0.0) ) 
    … …  
    472472                      IF ( zzdhdt < 0._wp ) THEN 
    473473                      ! For long timsteps factor in brackets slows the rapid collapse of the OSBL 
    474                          zpert   = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_rdt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj) 
     474                         zpert   = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_Dt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj) 
    475475                      ELSE 
    476                          zpert   = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_rdt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj) & 
     476                         zpert   = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_Dt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj) & 
    477477                              &  + MAX( zdb_bl(ji,jj), 0.0 ) 
    478478                      ENDIF 
    … …  
    487487      ibld(:,:) = 3 
    488488 
    489       zhbl_t(:,:) = hbl(:,:) + (zdhdt(:,:) - wn(ji,jj,ibld(ji,jj)))* rn_rdt ! certainly need wb here, so subtract it 
     489      zhbl_t(:,:) = hbl(:,:) + (zdhdt(:,:) - wn(ji,jj,ibld(ji,jj)))* rn_Dt ! certainly need wb here, so subtract it 
    490490      zhbl_t(:,:) = MIN(zhbl_t(:,:), ht_n(:,:)) 
    491       zdhdt(:,:) = MIN(zdhdt(:,:), (zhbl_t(:,:) - hbl(:,:))/rn_rdt + wn(ji,jj,ibld(ji,jj))) ! adjustment to represent limiting by ocean bottom 
     491      zdhdt(:,:) = MIN(zdhdt(:,:), (zhbl_t(:,:) - hbl(:,:))/rn_Dt + wn(ji,jj,ibld(ji,jj))) ! adjustment to represent limiting by ocean bottom 
    492492 
    493493      DO jk = 4, jpkm1 
    … …  
    516516               IF ( lconv(ji,jj) ) THEN 
    517517!unstable 
    518                   zvel_max =  - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / hbl(ji,jj) ) & 
     518                  zvel_max =  - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_Dt / hbl(ji,jj) ) & 
    519519                       &   * zwb_ent(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
    520520 
    … …  
    523523                          & - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ), 0.0 ) + zvel_max 
    524524 
    525                      zhbl_s = zhbl_s + MIN( - zwb_ent(ji,jj) / zdb * rn_rdt / FLOAT(ibld(ji,jj)-imld(ji,jj) ), e3w_n(ji,jj,jk) ) 
     525                     zhbl_s = zhbl_s + MIN( - zwb_ent(ji,jj) / zdb * rn_Dt / FLOAT(ibld(ji,jj)-imld(ji,jj) ), e3w_n(ji,jj,jk) ) 
    526526                     zhbl_s = MIN(zhbl_s, ht_n(ji,jj)) 
    527527 
    … …  
    13271327            IF ( iom_use("us_x") ) CALL iom_put( "us_x", tmask(:,:,1)*zustke*zcos_wind )   ! x surface Stokes drift 
    13281328            IF ( iom_use("us_y") ) CALL iom_put( "us_y", tmask(:,:,1)*zustke*zsin_wind )  ! y surface Stokes drift 
    1329             IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rau0*tmask(:,:,1)*zustar**2*zustke ) 
     1329            IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rho0*tmask(:,:,1)*zustar**2*zustke ) 
    13301330         ! Stokes drift read in from sbcwave  (=2). 
    13311331         CASE(2) 
    13321332            IF ( iom_use("us_x") ) CALL iom_put( "us_x", ut0sd )               ! x surface Stokes drift 
    13331333            IF ( iom_use("us_y") ) CALL iom_put( "us_y", vt0sd )               ! y surface Stokes drift 
    1334             IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rau0*tmask(:,:,1)*zustar**2* & 
     1334            IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rho0*tmask(:,:,1)*zustar**2* & 
    13351335                 & SQRT(ut0sd**2 + vt0sd**2 ) ) 
    13361336         END SELECT 
    … …  
    13481348         IF ( iom_use("zwstrl") ) CALL iom_put( "zwstrl", tmask(:,:,1)*zwstrl )         ! Langmuir velocity scale 
    13491349         IF ( iom_use("zustar") ) CALL iom_put( "zustar", tmask(:,:,1)*zustar )         ! friction velocity scale 
    1350          IF ( iom_use("wind_power") ) CALL iom_put( "wind_power", 1000.*rau0*tmask(:,:,1)*zustar**3 ) ! BL depth internal to zdf_osm routine 
    1351          IF ( iom_use("wind_wave_power") ) CALL iom_put( "wind_wave_power", 1000.*rau0*tmask(:,:,1)*zustar**2*zustke ) 
     1350         IF ( iom_use("wind_power") ) CALL iom_put( "wind_power", 1000.*rho0*tmask(:,:,1)*zustar**3 ) ! BL depth internal to zdf_osm routine 
     1351         IF ( iom_use("wind_wave_power") ) CALL iom_put( "wind_wave_power", 1000.*rho0*tmask(:,:,1)*zustar**2*zustke ) 
    13521352         IF ( iom_use("zhbl") ) CALL iom_put( "zhbl", tmask(:,:,1)*zhbl )               ! BL depth internal to zdf_osm routine 
    13531353         IF ( iom_use("zhml") ) CALL iom_put( "zhml", tmask(:,:,1)*zhml )               ! ML depth internal to zdf_osm routine 
    … …  
    15841584     imld_rst(:,:)  = nlb10         ! Initialization to the number of w ocean point 
    15851585     hbl(:,:)  = 0._wp              ! here hbl used as a dummy variable, integrating vertically N^2 
    1586      zN2_c = grav * rho_c * r1_rau0 ! convert density criteria into N^2 criteria 
     1586     zN2_c = grav * rho_c * r1_rho0 ! convert density criteria into N^2 criteria 
    15871587     ! 
    15881588     hbl(:,:)  = 0._wp              ! here hbl used as a dummy variable, integrating vertically N^2 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/ZDF/zdfric.F90

    r9598 r9939  
    181181         DO jj = 2, jpjm1        !* Ekman depth 
    182182            DO ji = 2, jpim1 
    183                zustar = SQRT( taum(ji,jj) * r1_rau0 ) 
     183               zustar = SQRT( taum(ji,jj) * r1_rho0 ) 
    184184               zhek   = rn_ekmfc * zustar / ( ABS( ff_t(ji,jj) ) + rsmall )   ! Ekman depth 
    185185               zh_ekm(ji,jj) = MAX(  rn_mldmin , MIN( zhek , rn_mldmax )  )   ! set allowed range 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/ZDF/zdftke.F90

    r9598 r9939  
    195195      REAL(wp) ::   zrhoa  = 1.22              ! Air density kg/m3 
    196196      REAL(wp) ::   zcdrag = 1.5e-3            ! drag coefficient 
    197       REAL(wp) ::   zbbrau, zri                ! local scalars 
     197      REAL(wp) ::   zbbrho, zri                ! local scalars 
    198198      REAL(wp) ::   zfact1, zfact2, zfact3     !   -         - 
    199199      REAL(wp) ::   ztx2  , zty2  , zcof       !   -         - 
    … …  
    206206      !!-------------------------------------------------------------------- 
    207207      ! 
    208       zbbrau = rn_ebb / rau0       ! Local constant initialisation 
    209       zfact1 = -.5_wp * rdt  
    210       zfact2 = 1.5_wp * rdt * rn_ediss 
    211       zfact3 = 0.5_wp       * rn_ediss 
     208      zbbrho = rn_ebb * r1_rho0       ! Local constant initialisation 
     209      zfact1 = -.5_wp * rn_Dt  
     210      zfact2 = 1.5_wp * rn_Dt * rn_ediss 
     211      zfact3 = 0.5_wp         * rn_ediss 
    212212      ! 
    213213      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    … …  
    215215      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    216216       
    217       DO jj = 2, jpjm1            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0) 
     217      DO jj = 2, jpjm1            ! en(1)   = rn_ebb taum / rho0  (min value rn_emin0) 
    218218         DO ji = fs_2, fs_jpim1   ! vector opt. 
    219             en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 
     219            en(ji,jj,1) = MAX( rn_emin0, zbbrho * taum(ji,jj) ) * tmask(ji,jj,1) 
    220220         END DO 
    221221      END DO 
    … …  
    232232      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    233233      ! 
    234       !   en(bot)   = (ebb0/rau0)*0.5*sqrt(u_botfr^2+v_botfr^2) (min value rn_emin) 
     234      !   en(bot)   = (ebb0/rho0)*0.5*sqrt(u_botfr^2+v_botfr^2) (min value rn_emin) 
    235235      ! where ebb0 does not includes surface wave enhancement (i.e. ebb0=3.75) 
    236236      ! Note that stress averaged is done using an wet-only calculation of u and v at t-point like in zdfsh2 
    … …  
    242242               zmsku = ( 2. - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 
    243243               zmskv = ( 2. - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) 
    244                !                       ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0) 
     244               !                       ! where 0.001875 = (rn_ebb0/rho0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0) 
    245245               zebot = - 0.001875_wp * rCdU_bot(ji,jj) * SQRT(  ( zmsku*( ub(ji,jj,mbkt(ji,jj))+ub(ji-1,jj,mbkt(ji,jj)) ) )**2  & 
    246246                  &                                           + ( zmskv*( vb(ji,jj,mbkt(ji,jj))+vb(ji,jj-1,mbkt(ji,jj)) ) )**2  ) 
    … …  
    253253                  zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) ) 
    254254                  zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) ) 
    255                   !                             ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000.  (CAUTION CdU<0) 
     255                  !                             ! where 0.001875 = (rn_ebb0/rho0) * 0.5 = 3.75*0.5/1000.  (CAUTION CdU<0) 
    256256                  zetop = - 0.001875_wp * rCdU_top(ji,jj) * SQRT(  ( zmsku*( ub(ji,jj,mikt(ji,jj))+ub(ji-1,jj,mikt(ji,jj)) ) )**2  & 
    257257                     &                                           + ( zmskv*( vb(ji,jj,mikt(ji,jj))+vb(ji,jj-1,mikt(ji,jj)) ) )**2  ) 
    … …  
    298298                  zwlc = zind * rn_lc * zus * SIN( rpi * pdepw(ji,jj,jk) / zhlc(ji,jj) ) 
    299299                  !                                           ! TKE Langmuir circulation source term 
    300                   en(ji,jj,jk) = en(ji,jj,jk) + rdt * MAX(0.,1._wp - 4.*fr_i(ji,jj) ) * ( zwlc * zwlc * zwlc )   & 
    301                      &                              / zhlc(ji,jj) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
     300                  en(ji,jj,jk) = en(ji,jj,jk) + rn_Dt * MAX(0.,1._wp - 4.*fr_i(ji,jj) ) * ( zwlc * zwlc * zwlc )   & 
     301                     &                                / zhlc(ji,jj) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    302302               END DO 
    303303            END DO 
    … …  
    342342               ! 
    343343               !                                   ! right hand side in en 
    344                en(ji,jj,jk) = en(ji,jj,jk) + rdt * (  p_sh2(ji,jj,jk)                          &   ! shear 
    345                   &                                 - p_avt(ji,jj,jk) * rn2(ji,jj,jk)          &   ! stratification 
    346                   &                                 + zfact3 * dissl(ji,jj,jk) * en(ji,jj,jk)  &   ! dissipation 
    347                   &                                ) * wmask(ji,jj,jk) 
     344               en(ji,jj,jk) = en(ji,jj,jk) + rn_Dt * (  p_sh2(ji,jj,jk)                          &   ! shear 
     345                  &                                   - p_avt(ji,jj,jk) * rn2(ji,jj,jk)          &   ! stratification 
     346                  &                                   + zfact3 * dissl(ji,jj,jk) * en(ji,jj,jk)  &   ! dissipation 
     347                  &                                  ) * wmask(ji,jj,jk) 
    348348            END DO 
    349349         END DO 
    … …  
    422422                  zdif = taum(ji,jj) - ztau                            ! mean of modulus - modulus of the mean  
    423423                  zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications... 
    424                   en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) )   & 
     424                  en(ji,jj,jk) = en(ji,jj,jk) + zbbrho * zdif * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) )   & 
    425425                     &                        * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    426426               END DO 
    … …  
    473473      ! 
    474474      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    475       REAL(wp) ::   zrn2, zraug, zcoef, zav   ! local scalars 
     475      REAL(wp) ::   zrn2, zrhog, zcoef, zav   ! local scalars 
    476476      REAL(wp) ::   zdku,   zdkv, zsqen       !   -      - 
    477477      REAL(wp) ::   zemxl, zemlm, zemlp       !   -      - 
    … …  
    489489      zmxld(:,:,:)  = rmxl_min 
    490490      ! 
    491       IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g) 
    492          zraug = vkarmn * 2.e5_wp / ( rau0 * grav ) 
     491      IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rho0*g) 
     492         zrhog = vkarmn * 2.e5_wp / ( rho0 * grav ) 
    493493         DO jj = 2, jpjm1 
    494494            DO ji = fs_2, fs_jpim1 
    495                zmxlm(ji,jj,1) = MAX( rn_mxl0, zraug * taum(ji,jj) * tmask(ji,jj,1) ) 
     495               zmxlm(ji,jj,1) = MAX( rn_mxl0, zrhog * taum(ji,jj) * tmask(ji,jj,1) ) 
    496496            END DO 
    497497         END DO 
Note: See TracChangeset for help on using the changeset viewer.