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 13899 for NEMO/branches/2020/tickets_icb_1900/src/OCE/ZDF/zdfiwm.F90 – NEMO

Ignore:
Timestamp:
2020-11-27T17:26:33+01:00 (4 years ago)
Author:
mathiot
Message:

ticket #1900: update branch to trunk and add ICB test case

Location:
NEMO/branches/2020/tickets_icb_1900
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/tickets_icb_1900

    • Property svn:externals
      •  

        old new  
        22^/utils/build/makenemo@HEAD   makenemo 
        33^/utils/build/mk@HEAD         mk 
        4 ^/utils/tools/@HEAD           tools 
         4^/utils/tools@HEAD            tools 
        55^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS      ext/AGRIF 
        66^/vendors/FCM@HEAD            ext/FCM 
         
        88 
        99# SETTE 
        10 ^/utils/CI/sette@12931        sette 
         10^/utils/CI/sette@13559        sette 
  • NEMO/branches/2020/tickets_icb_1900/src/OCE/ZDF/zdfiwm.F90

    r13237 r13899  
    140140      !!---------------------------------------------------------------------- 
    141141      ! 
    142       !                       !* Set to zero the 1st and last vertical levels of appropriate variables 
    143       zemx_iwm (:,:,1) = 0._wp   ;   zemx_iwm (:,:,jpk) = 0._wp 
    144       zav_ratio(:,:,1) = 0._wp   ;   zav_ratio(:,:,jpk) = 0._wp 
    145       zav_wave (:,:,1) = 0._wp   ;   zav_wave (:,:,jpk) = 0._wp 
     142      !                        
     143      ! Set to zero the 1st and last vertical levels of appropriate variables 
     144      IF( iom_use("emix_iwm") ) THEN 
     145         DO_2D( 0, 0, 0, 0 ) 
     146            zemx_iwm (ji,jj,1) = 0._wp   ;   zemx_iwm (ji,jj,jpk) = 0._wp 
     147         END_2D 
     148      ENDIF 
     149      IF( iom_use("av_ratio") ) THEN 
     150         DO_2D( 0, 0, 0, 0 ) 
     151            zav_ratio(ji,jj,1) = 0._wp   ;   zav_ratio(ji,jj,jpk) = 0._wp 
     152         END_2D 
     153      ENDIF 
     154      IF( iom_use("av_wave") .OR. sn_cfctl%l_prtctl ) THEN 
     155         DO_2D( 0, 0, 0, 0 ) 
     156            zav_wave (ji,jj,1) = 0._wp   ;   zav_wave (ji,jj,jpk) = 0._wp 
     157         END_2D 
     158      ENDIF 
    146159      ! 
    147160      !                       ! ----------------------------- ! 
     
    151164      !                       !* Critical slope mixing: distribute energy over the time-varying ocean depth, 
    152165      !                                                 using an exponential decay from the seafloor. 
    153       DO_2D_11_11 
     166      DO_2D( 0, 0, 0, 0 )             ! part independent of the level 
    154167         zhdep(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1)       ! depth of the ocean 
    155168         zfact(ji,jj) = rho0 * (  1._wp - EXP( -zhdep(ji,jj) / hcri_iwm(ji,jj) )  ) 
     
    157170      END_2D 
    158171!!gm gde3w ==>>>  check for ssh taken into account.... seem OK gde3w_n=gdept(:,:,:,Kmm) - ssh(:,:,Kmm) 
    159       DO_3D_11_11( 2, jpkm1 ) 
     172      DO_3D( 0, 0, 0, 0, 2, jpkm1 )   ! complete with the level-dependent part 
    160173         IF ( zfact(ji,jj) == 0._wp .OR. wmask(ji,jj,jk) == 0._wp ) THEN   ! optimization 
    161174            zemx_iwm(ji,jj,jk) = 0._wp 
     
    177190      CASE ( 1 )               ! Dissipation scales as N (recommended) 
    178191         ! 
    179          zfact(:,:) = 0._wp 
    180          DO jk = 2, jpkm1              ! part independent of the level 
    181             zfact(:,:) =   & 
    182                &  zfact(:,:) +   & 
    183                &  e3w(:,:,jk,Kmm) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk) 
    184          END DO 
    185          ! 
    186          DO_2D_11_11 
     192         DO_2D( 0, 0, 0, 0 ) 
     193            zfact(ji,jj) = 0._wp 
     194         END_2D 
     195         DO_3D( 0, 0, 0, 0, 2, jpkm1 )       ! part independent of the level 
     196            zfact(ji,jj) = zfact(ji,jj) + e3w(ji,jj,jk,Kmm) * SQRT(  MAX( 0._wp, rn2(ji,jj,jk) )  ) * wmask(ji,jj,jk) 
     197         END_3D 
     198         ! 
     199         DO_2D( 0, 0, 0, 0 ) 
    187200            IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rho0 * zfact(ji,jj) ) 
    188201         END_2D 
    189202         ! 
    190          DO jk = 2, jpkm1              ! complete with the level-dependent part 
    191             zemx_iwm(:,:,jk) = zemx_iwm(:,:,jk) + zfact(:,:) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk) 
    192          END DO 
     203         DO_3D( 0, 0, 0, 0, 2, jpkm1 )       ! complete with the level-dependent part 
     204            zemx_iwm(ji,jj,jk) = zemx_iwm(ji,jj,jk) + zfact(ji,jj) * SQRT(  MAX( 0._wp, rn2(ji,jj,jk) )  ) * wmask(ji,jj,jk) 
     205         END_3D 
    193206         ! 
    194207      CASE ( 2 )               ! Dissipation scales as N^2 
    195208         ! 
    196          zfact(:,:) = 0._wp 
    197          DO jk = 2, jpkm1              ! part independent of the level 
    198             zfact(:,:) = zfact(:,:) + e3w(:,:,jk,Kmm) * MAX( 0._wp, rn2(:,:,jk) ) * wmask(:,:,jk) 
    199          END DO 
    200          ! 
    201          DO_2D_11_11 
     209         DO_2D( 0, 0, 0, 0 ) 
     210            zfact(ji,jj) = 0._wp 
     211         END_2D 
     212         DO_3D( 0, 0, 0, 0, 2, jpkm1 )       ! part independent of the level 
     213            zfact(ji,jj) = zfact(ji,jj) + e3w(ji,jj,jk,Kmm) * MAX( 0._wp, rn2(ji,jj,jk) ) * wmask(ji,jj,jk) 
     214         END_3D 
     215         ! 
     216         DO_2D( 0, 0, 0, 0 ) 
    202217            IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rho0 * zfact(ji,jj) ) 
    203218         END_2D 
    204219         ! 
    205          DO jk = 2, jpkm1              ! complete with the level-dependent part 
    206             zemx_iwm(:,:,jk) = zemx_iwm(:,:,jk) + zfact(:,:) * MAX( 0._wp, rn2(:,:,jk) ) * wmask(:,:,jk) 
    207          END DO 
     220         DO_3D( 0, 0, 0, 0, 2, jpkm1 )       ! complete with the level-dependent part 
     221            zemx_iwm(ji,jj,jk) = zemx_iwm(ji,jj,jk) + zfact(ji,jj) * MAX( 0._wp, rn2(ji,jj,jk) ) * wmask(ji,jj,jk) 
     222         END_3D 
    208223         ! 
    209224      END SELECT 
     
    212227      !                        !* ocean depth as proportional to rn2 * exp(-z_wkb/rn_hbot) 
    213228      ! 
    214       zwkb (:,:,:) = 0._wp 
    215       zfact(:,:)   = 0._wp 
    216       DO jk = 2, jpkm1 
    217          zfact(:,:) = zfact(:,:) + e3w(:,:,jk,Kmm) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk) 
    218          zwkb(:,:,jk) = zfact(:,:) 
    219       END DO 
    220 !!gm even better: 
    221 !      DO jk = 2, jpkm1 
    222 !         zwkb(:,:) = zwkb(:,:) + e3w(:,:,jk,Kmm) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) 
    223 !      END DO 
    224 !      zfact(:,:) = zwkb(:,:,jpkm1) 
    225 !!gm or just use zwkb(k=jpk-1) instead of zfact... 
    226 !!gm 
    227       ! 
    228       DO_3D_11_11( 2, jpkm1 ) 
     229      DO_2D( 0, 0, 0, 0 ) 
     230         zwkb(ji,jj,1) = 0._wp 
     231      END_2D 
     232      DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     233         zwkb(ji,jj,jk) = zwkb(ji,jj,jk-1) + e3w(ji,jj,jk,Kmm) * SQRT(  MAX( 0._wp, rn2(ji,jj,jk) )  ) * wmask(ji,jj,jk) 
     234      END_3D 
     235      DO_2D( 0, 0, 0, 0 ) 
     236         zfact(ji,jj) = zwkb(ji,jj,jpkm1) 
     237      END_2D 
     238      ! 
     239      DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
    229240         IF( zfact(ji,jj) /= 0 )   zwkb(ji,jj,jk) = zhdep(ji,jj) * ( zfact(ji,jj) - zwkb(ji,jj,jk) )   & 
    230241            &                                     * wmask(ji,jj,jk) / zfact(ji,jj) 
    231242      END_3D 
    232       zwkb(:,:,1) = zhdep(:,:) * wmask(:,:,1) 
    233       ! 
    234       DO_3D_11_11( 2, jpkm1 ) 
    235          IF ( rn2(ji,jj,jk) <= 0._wp .OR. wmask(ji,jj,jk) == 0._wp ) THEN   ! optimization 
     243      DO_2D( 0, 0, 0, 0 ) 
     244         zwkb (ji,jj,1) = zhdep(ji,jj) * wmask(ji,jj,1) 
     245      END_2D 
     246      ! 
     247      DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     248         IF ( rn2(ji,jj,jk) <= 0._wp .OR. wmask(ji,jj,jk) == 0._wp ) THEN   ! optimization: EXP coast a lot 
    236249            zweight(ji,jj,jk) = 0._wp 
    237250         ELSE 
     
    241254      END_3D 
    242255      ! 
    243       zfact(:,:) = 0._wp 
    244       DO jk = 2, jpkm1              ! part independent of the level 
    245          zfact(:,:) = zfact(:,:) + zweight(:,:,jk) 
    246       END DO 
    247       ! 
    248       DO_2D_11_11 
     256      DO_2D( 0, 0, 0, 0 ) 
     257         zfact(ji,jj) = 0._wp 
     258      END_2D 
     259      DO_3D( 0, 0, 0, 0, 2, jpkm1 )       ! part independent of the level 
     260         zfact(ji,jj) = zfact(ji,jj) + zweight(ji,jj,jk) 
     261      END_3D 
     262      ! 
     263      DO_2D( 0, 0, 0, 0 ) 
    249264         IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = ebot_iwm(ji,jj) / ( rho0 * zfact(ji,jj) ) 
    250265      END_2D 
    251266      ! 
    252       DO jk = 2, jpkm1              ! complete with the level-dependent part 
    253          zemx_iwm(:,:,jk) = zemx_iwm(:,:,jk) + zweight(:,:,jk) * zfact(:,:) * wmask(:,:,jk)   & 
    254             &                                / ( gde3w(:,:,jk) - gde3w(:,:,jk-1) ) 
    255 !!gm  use of e3t(:,:,:,Kmm) just above? 
    256       END DO 
     267      DO_3D( 0, 0, 0, 0, 2, jpkm1 )       ! complete with the level-dependent part 
     268         zemx_iwm(ji,jj,jk) = zemx_iwm(ji,jj,jk) + zweight(ji,jj,jk) * zfact(ji,jj) * wmask(ji,jj,jk)   & 
     269            &                                                        / ( gde3w(ji,jj,jk) - gde3w(ji,jj,jk-1) ) 
     270!!gm  use of e3t(ji,jj,:,Kmm) just above? 
     271      END_3D 
    257272      ! 
    258273!!gm  this is to be replaced by just a constant value znu=1.e-6 m2/s 
    259274      ! Calculate molecular kinematic viscosity 
    260       znu_t(:,:,:) = 1.e-4_wp * (  17.91_wp - 0.53810_wp * ts(:,:,:,jp_tem,Kmm) + 0.00694_wp * ts(:,:,:,jp_tem,Kmm) * ts(:,:,:,jp_tem,Kmm)  & 
    261          &                                  + 0.02305_wp * ts(:,:,:,jp_sal,Kmm)  ) * tmask(:,:,:) * r1_rho0 
    262       DO jk = 2, jpkm1 
    263          znu_w(:,:,jk) = 0.5_wp * ( znu_t(:,:,jk-1) + znu_t(:,:,jk) ) * wmask(:,:,jk) 
    264       END DO 
     275      DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     276         znu_t(ji,jj,jk) = 1.e-4_wp * (  17.91_wp - 0.53810_wp * ts(ji,jj,jk,jp_tem,Kmm)   & 
     277            &                                     + 0.00694_wp * ts(ji,jj,jk,jp_tem,Kmm) * ts(ji,jj,jk,jp_tem,Kmm)  & 
     278            &                                     + 0.02305_wp * ts(ji,jj,jk,jp_sal,Kmm)  ) * tmask(ji,jj,jk) * r1_rho0 
     279      END_3D 
     280      DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     281         znu_w(ji,jj,jk) = 0.5_wp * ( znu_t(ji,jj,jk-1) + znu_t(ji,jj,jk) ) * wmask(ji,jj,jk) 
     282      END_3D 
    265283!!gm end 
    266284      ! 
    267285      ! Calculate turbulence intensity parameter Reb 
    268       DO jk = 2, jpkm1 
    269          zReb(:,:,jk) = zemx_iwm(:,:,jk) / MAX( 1.e-20_wp, znu_w(:,:,jk) * rn2(:,:,jk) ) 
    270       END DO 
     286      DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     287         zReb(ji,jj,jk) = zemx_iwm(ji,jj,jk) / MAX( 1.e-20_wp, znu_w(ji,jj,jk) * rn2(ji,jj,jk) ) 
     288      END_3D 
    271289      ! 
    272290      ! Define internal wave-induced diffusivity 
    273       DO jk = 2, jpkm1 
    274          zav_wave(:,:,jk) = znu_w(:,:,jk) * zReb(:,:,jk) * r1_6   ! This corresponds to a constant mixing efficiency of 1/6 
    275       END DO 
    276       ! 
    277       IF( ln_mevar ) THEN              ! Variable mixing efficiency case : modify zav_wave in the 
    278          DO_3D_11_11( 2, jpkm1 ) 
     291      DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     292         zav_wave(ji,jj,jk) = znu_w(ji,jj,jk) * zReb(ji,jj,jk) * r1_6   ! This corresponds to a constant mixing efficiency of 1/6 
     293      END_3D 
     294      ! 
     295      IF( ln_mevar ) THEN                ! Variable mixing efficiency case : modify zav_wave in the 
     296         DO_3D( 0, 0, 0, 0, 2, jpkm1 )   ! energetic (Reb > 480) and buoyancy-controlled (Reb <10.224 ) regimes 
    279297            IF( zReb(ji,jj,jk) > 480.00_wp ) THEN 
    280298               zav_wave(ji,jj,jk) = 3.6515_wp * znu_w(ji,jj,jk) * SQRT( zReb(ji,jj,jk) ) 
     
    285303      ENDIF 
    286304      ! 
    287       DO jk = 2, jpkm1                 ! Bound diffusivity by molecular value and 100 cm2/s 
    288          zav_wave(:,:,jk) = MIN(  MAX( 1.4e-7_wp, zav_wave(:,:,jk) ), 1.e-2_wp  ) * wmask(:,:,jk) 
    289       END DO 
     305      DO_3D( 0, 0, 0, 0, 2, jpkm1 )      ! Bound diffusivity by molecular value and 100 cm2/s 
     306         zav_wave(ji,jj,jk) = MIN(  MAX( 1.4e-7_wp, zav_wave(ji,jj,jk) ), 1.e-2_wp  ) * wmask(ji,jj,jk) 
     307      END_3D 
    290308      ! 
    291309      IF( kt == nit000 ) THEN        !* Control print at first time-step: diagnose the energy consumed by zav_wave 
    292310         zztmp = 0._wp 
    293311!!gm used of glosum 3D.... 
    294          DO_3D_11_11( 2, jpkm1 ) 
     312         DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
    295313            zztmp = zztmp + e3w(ji,jj,jk,Kmm) * e1e2t(ji,jj)   & 
    296314               &          * MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj) 
     
    312330      !                          ! ----------------------- ! 
    313331      !       
    314       IF( ln_tsdiff ) THEN          !* Option for differential mixing of salinity and temperature 
     332      IF( ln_tsdiff ) THEN                !* Option for differential mixing of salinity and temperature 
    315333         ztmp1 = 0.505_wp + 0.495_wp * TANH( 0.92_wp * ( LOG10( 1.e-20_wp ) - 0.60_wp ) ) 
    316          DO_3D_11_11( 2, jpkm1 ) 
     334         DO_3D( 0, 0, 0, 0, 2, jpkm1 )       ! Calculate S/T diffusivity ratio as a function of Reb 
    317335            ztmp2 = zReb(ji,jj,jk) * 5._wp * r1_6 
    318336            IF ( ztmp2 > 1.e-20_wp .AND. wmask(ji,jj,jk) == 1._wp ) THEN 
     
    323341         END_3D 
    324342         CALL iom_put( "av_ratio", zav_ratio ) 
    325          DO jk = 2, jpkm1           !* update momentum & tracer diffusivity with wave-driven mixing 
    326             p_avs(:,:,jk) = p_avs(:,:,jk) + zav_wave(:,:,jk) * zav_ratio(:,:,jk) 
    327             p_avt(:,:,jk) = p_avt(:,:,jk) + zav_wave(:,:,jk) 
    328             p_avm(:,:,jk) = p_avm(:,:,jk) + zav_wave(:,:,jk) 
    329          END DO 
    330          ! 
    331       ELSE                          !* update momentum & tracer diffusivity with wave-driven mixing 
    332          DO jk = 2, jpkm1 
    333             p_avs(:,:,jk) = p_avs(:,:,jk) + zav_wave(:,:,jk) 
    334             p_avt(:,:,jk) = p_avt(:,:,jk) + zav_wave(:,:,jk) 
    335             p_avm(:,:,jk) = p_avm(:,:,jk) + zav_wave(:,:,jk) 
    336          END DO 
    337       ENDIF 
    338  
    339       !                             !* output internal wave-driven mixing coefficient 
     343         DO_3D( 0, 0, 0, 0, 2, jpkm1 )    !* update momentum & tracer diffusivity with wave-driven mixing 
     344            p_avs(ji,jj,jk) = p_avs(ji,jj,jk) + zav_wave(ji,jj,jk) * zav_ratio(ji,jj,jk) 
     345            p_avt(ji,jj,jk) = p_avt(ji,jj,jk) + zav_wave(ji,jj,jk) 
     346            p_avm(ji,jj,jk) = p_avm(ji,jj,jk) + zav_wave(ji,jj,jk) 
     347         END_3D 
     348         ! 
     349      ELSE                                !* update momentum & tracer diffusivity with wave-driven mixing 
     350         DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     351            p_avs(ji,jj,jk) = p_avs(ji,jj,jk) + zav_wave(ji,jj,jk) 
     352            p_avt(ji,jj,jk) = p_avt(ji,jj,jk) + zav_wave(ji,jj,jk) 
     353            p_avm(ji,jj,jk) = p_avm(ji,jj,jk) + zav_wave(ji,jj,jk) 
     354         END_3D 
     355      ENDIF 
     356 
     357      !                                   !* output internal wave-driven mixing coefficient 
    340358      CALL iom_put( "av_wave", zav_wave ) 
    341                                     !* output useful diagnostics: Kz*N^2 ,  
     359                                          !* output useful diagnostics: Kz*N^2 ,  
    342360!!gm Kz*N2 should take into account the ratio avs/avt if it is used.... (see diaar5) 
    343                                     !  vertical integral of rho0 * Kz * N^2 , energy density (zemx_iwm) 
     361                                          !  vertical integral of rho0 * Kz * N^2 , energy density (zemx_iwm) 
    344362      IF( iom_use("bflx_iwm") .OR. iom_use("pcmap_iwm") ) THEN 
    345363         ALLOCATE( z2d(jpi,jpj) , z3d(jpi,jpj,jpk) ) 
    346          z3d(:,:,:) = MAX( 0._wp, rn2(:,:,:) ) * zav_wave(:,:,:) 
    347          z2d(:,:) = 0._wp 
    348          DO jk = 2, jpkm1 
    349             z2d(:,:) = z2d(:,:) + e3w(:,:,jk,Kmm) * z3d(:,:,jk) * wmask(:,:,jk) 
    350          END DO 
    351          z2d(:,:) = rho0 * z2d(:,:) 
    352          CALL iom_put( "bflx_iwm", z3d ) 
     364         ! Initialisation for iom_put 
     365         DO_2D( 0, 0, 0, 0 ) 
     366            z3d(ji,jj,1) = 0._wp   ;   z3d(ji,jj,jpk) = 0._wp 
     367         END_2D 
     368         z3d(           1:nn_hls,:,:) = 0._wp   ;   z3d(:,           1:nn_hls,:) = 0._wp 
     369         z3d(jpi-nn_hls+1:jpi   ,:,:) = 0._wp   ;   z3d(:,jpj-nn_hls+1:   jpj,:) = 0._wp 
     370         z2d(           1:nn_hls,:  ) = 0._wp   ;   z2d(:,           1:nn_hls  ) = 0._wp 
     371         z2d(jpi-nn_hls+1:jpi   ,:  ) = 0._wp   ;   z2d(:,jpj-nn_hls+1:   jpj  ) = 0._wp 
     372 
     373         DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     374            z3d(ji,jj,jk) = MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk) 
     375         END_3D 
     376         DO_2D( 0, 0, 0, 0 ) 
     377            z2d(ji,jj) = 0._wp 
     378         END_2D 
     379         DO_3D( 0, 0, 0, 0, 2, jpkm1 )  
     380            z2d(ji,jj) = z2d(ji,jj) + e3w(ji,jj,jk,Kmm) * z3d(ji,jj,jk) * wmask(ji,jj,jk) 
     381         END_3D 
     382         DO_2D( 0, 0, 0, 0 ) 
     383            z2d(ji,jj) = rho0 * z2d(ji,jj) 
     384         END_2D 
     385         CALL iom_put(  "bflx_iwm", z3d ) 
    353386         CALL iom_put( "pcmap_iwm", z2d ) 
    354387         DEALLOCATE( z2d , z3d ) 
Note: See TracChangeset for help on using the changeset viewer.