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 12960 for NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/OCE/ZDF/zdfiwm.F90 – NEMO

Ignore:
Timestamp:
2020-05-22T09:05:34+02:00 (4 years ago)
Author:
smasson
Message:

Extra_Halo: additional bugfixes and developments, see #2366

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/OCE/ZDF/zdfiwm.F90

    r12510 r12960  
    139139      !!---------------------------------------------------------------------- 
    140140      ! 
    141       !                       !* Set to zero the 1st and last vertical levels of appropriate variables 
    142       zemx_iwm (:,:,1) = 0._wp   ;   zemx_iwm (:,:,jpk) = 0._wp 
    143       zav_ratio(:,:,1) = 0._wp   ;   zav_ratio(:,:,jpk) = 0._wp 
    144       zav_wave (:,:,1) = 0._wp   ;   zav_wave (:,:,jpk) = 0._wp 
     141      !                        
     142      ! Set to zero the 1st and last vertical levels of appropriate variables 
     143      IF( iom_use("emix_iwm") ) THEN 
     144         DO_2D_00_00 
     145            zemx_iwm (ji,jj,1) = 0._wp   ;   zemx_iwm (ji,jj,jpk) = 0._wp 
     146         END_2D 
     147         zemx_iwm (           1:nn_hls,:,:) = 0._wp   ;   zemx_iwm (:,           1:nn_hls,:) = 0._wp 
     148         zemx_iwm (jpi-nn_hls+1:jpi   ,:,:) = 0._wp   ;   zemx_iwm (:,jpj-nn_hls+1:   jpj,:) = 0._wp 
     149      ENDIF 
     150      IF( iom_use("av_ratio") ) THEN 
     151         DO_2D_00_00 
     152            zav_ratio(ji,jj,1) = 0._wp   ;   zav_ratio(ji,jj,jpk) = 0._wp 
     153         END_2D 
     154         zav_ratio(           1:nn_hls,:,:) = 0._wp   ;   zav_ratio(:,           1:nn_hls,:) = 0._wp 
     155         zav_ratio(jpi-nn_hls+1:jpi   ,:,:) = 0._wp   ;   zav_ratio(:,jpj-nn_hls+1:   jpj,:) = 0._wp 
     156      ENDIF 
     157      IF( iom_use("av_wave") ) THEN 
     158         DO_2D_00_00 
     159            zav_wave (ji,jj,1) = 0._wp   ;   zav_wave (ji,jj,jpk) = 0._wp 
     160         END_2D 
     161         zav_wave(           1:nn_hls,:,:) = 0._wp   ;   zav_wave(:,           1:nn_hls,:) = 0._wp 
     162         zav_wave(jpi-nn_hls+1:jpi   ,:,:) = 0._wp   ;   zav_wave(:,jpj-nn_hls+1:   jpj,:) = 0._wp 
     163      ENDIF 
    145164      ! 
    146165      !                       ! ----------------------------- ! 
     
    150169      !                       !* Critical slope mixing: distribute energy over the time-varying ocean depth, 
    151170      !                                                 using an exponential decay from the seafloor. 
    152       DO_2D_11_11 
     171      DO_2D_00_00 
    153172         zhdep(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1)       ! depth of the ocean 
    154173         zfact(ji,jj) = rho0 * (  1._wp - EXP( -zhdep(ji,jj) / hcri_iwm(ji,jj) )  ) 
     
    156175      END_2D 
    157176!!gm gde3w ==>>>  check for ssh taken into account.... seem OK gde3w_n=gdept(:,:,:,Kmm) - ssh(:,:,Kmm) 
    158       DO_3D_11_11( 2, jpkm1 ) 
     177      DO_3D_00_00( 2, jpkm1 ) 
    159178         IF ( zfact(ji,jj) == 0._wp .OR. wmask(ji,jj,jk) == 0._wp ) THEN   ! optimization 
    160179            zemx_iwm(ji,jj,jk) = 0._wp 
     
    176195      CASE ( 1 )               ! Dissipation scales as N (recommended) 
    177196         ! 
    178          zfact(:,:) = 0._wp 
    179          DO jk = 2, jpkm1              ! part independent of the level 
    180             zfact(:,:) = zfact(:,:) + e3w(:,:,jk,Kmm) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk) 
    181          END DO 
    182          ! 
    183          DO_2D_11_11 
     197         DO_2D_00_00 
     198            zfact(ji,jj) = 0._wp 
     199         END_2D 
     200         DO_3D_00_00( 2, jpkm1 )       ! part independent of the level 
     201            zfact(ji,jj) = zfact(ji,jj) + e3w(ji,jj,jk,Kmm) * SQRT(  MAX( 0._wp, rn2(ji,jj,jk) )  ) * wmask(ji,jj,jk) 
     202         END_3D 
     203         ! 
     204         DO_2D_00_00 
    184205            IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rho0 * zfact(ji,jj) ) 
    185206         END_2D 
    186207         ! 
    187          DO jk = 2, jpkm1              ! complete with the level-dependent part 
    188             zemx_iwm(:,:,jk) = zemx_iwm(:,:,jk) + zfact(:,:) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk) 
    189          END DO 
     208         DO_3D_00_00( 2, jpkm1 )       ! complete with the level-dependent part 
     209            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) 
     210         END_3D 
    190211         ! 
    191212      CASE ( 2 )               ! Dissipation scales as N^2 
    192213         ! 
    193          zfact(:,:) = 0._wp 
    194          DO jk = 2, jpkm1              ! part independent of the level 
    195             zfact(:,:) = zfact(:,:) + e3w(:,:,jk,Kmm) * MAX( 0._wp, rn2(:,:,jk) ) * wmask(:,:,jk) 
    196          END DO 
    197          ! 
    198          DO_2D_11_11 
     214         DO_2D_00_00 
     215            zfact(ji,jj) = 0._wp 
     216         END_2D 
     217         DO_3D_00_00( 2, jpkm1 )       ! part independent of the level 
     218            zfact(ji,jj) = zfact(ji,jj) + e3w(ji,jj,jk,Kmm) * MAX( 0._wp, rn2(ji,jj,jk) ) * wmask(ji,jj,jk) 
     219         END_3D 
     220         ! 
     221         DO_2D_00_00 
    199222            IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rho0 * zfact(ji,jj) ) 
    200223         END_2D 
    201224         ! 
    202          DO jk = 2, jpkm1              ! complete with the level-dependent part 
    203             zemx_iwm(:,:,jk) = zemx_iwm(:,:,jk) + zfact(:,:) * MAX( 0._wp, rn2(:,:,jk) ) * wmask(:,:,jk) 
    204          END DO 
     225         DO_3D_00_00( 2, jpkm1 )       ! complete with the level-dependent part 
     226            zemx_iwm(ji,jj,jk) = zemx_iwm(ji,jj,jk) + zfact(ji,jj) * MAX( 0._wp, rn2(ji,jj,jk) ) * wmask(ji,jj,jk) 
     227         END_3D 
    205228         ! 
    206229      END SELECT 
     
    209232      !                        !* ocean depth as proportional to rn2 * exp(-z_wkb/rn_hbot) 
    210233      ! 
    211       zwkb (:,:,:) = 0._wp 
    212       zfact(:,:)   = 0._wp 
    213       DO jk = 2, jpkm1 
    214          zfact(:,:) = zfact(:,:) + e3w(:,:,jk,Kmm) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk) 
    215          zwkb(:,:,jk) = zfact(:,:) 
    216       END DO 
    217 !!gm even better: 
    218 !      DO jk = 2, jpkm1 
    219 !         zwkb(:,:) = zwkb(:,:) + e3w(:,:,jk,Kmm) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) 
    220 !      END DO 
    221 !      zfact(:,:) = zwkb(:,:,jpkm1) 
    222 !!gm or just use zwkb(k=jpk-1) instead of zfact... 
    223 !!gm 
    224       ! 
    225       DO_3D_11_11( 2, jpkm1 ) 
     234      DO_2D_00_00 
     235         zwkb(ji,jj,1) = 0._wp 
     236      END_2D 
     237      DO_3D_00_00( 2, jpkm1 ) 
     238         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) 
     239      END_3D 
     240      DO_2D_00_00 
     241         zfact(ji,jj) = zwkb(ji,jj,jpkm1) 
     242      END_2D 
     243      ! 
     244      DO_3D_00_00( 2, jpkm1 ) 
    226245         IF( zfact(ji,jj) /= 0 )   zwkb(ji,jj,jk) = zhdep(ji,jj) * ( zfact(ji,jj) - zwkb(ji,jj,jk) )   & 
    227246            &                                     * wmask(ji,jj,jk) / zfact(ji,jj) 
    228247      END_3D 
    229       zwkb(:,:,1) = zhdep(:,:) * wmask(:,:,1) 
    230       ! 
    231       DO_3D_11_11( 2, jpkm1 ) 
    232          IF ( rn2(ji,jj,jk) <= 0._wp .OR. wmask(ji,jj,jk) == 0._wp ) THEN   ! optimization 
     248      DO_2D_00_00 
     249         zwkb (ji,jj,1) = zhdep(ji,jj) * wmask(ji,jj,1) 
     250      END_2D 
     251      ! 
     252      DO_3D_00_00( 2, jpkm1 ) 
     253         IF ( rn2(ji,jj,jk) <= 0._wp .OR. wmask(ji,jj,jk) == 0._wp ) THEN   ! optimization: EXP coast a lot 
    233254            zweight(ji,jj,jk) = 0._wp 
    234255         ELSE 
     
    238259      END_3D 
    239260      ! 
    240       zfact(:,:) = 0._wp 
    241       DO jk = 2, jpkm1              ! part independent of the level 
    242          zfact(:,:) = zfact(:,:) + zweight(:,:,jk) 
    243       END DO 
    244       ! 
    245       DO_2D_11_11 
     261      DO_2D_00_00 
     262         zfact(ji,jj) = 0._wp 
     263      END_2D 
     264      DO_3D_00_00( 2, jpkm1 )       ! part independent of the level 
     265         zfact(ji,jj) = zfact(ji,jj) + zweight(ji,jj,jk) 
     266      END_3D 
     267      ! 
     268      DO_2D_00_00 
    246269         IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = ebot_iwm(ji,jj) / ( rho0 * zfact(ji,jj) ) 
    247270      END_2D 
    248271      ! 
    249       DO jk = 2, jpkm1              ! complete with the level-dependent part 
    250          zemx_iwm(:,:,jk) = zemx_iwm(:,:,jk) + zweight(:,:,jk) * zfact(:,:) * wmask(:,:,jk)   & 
    251             &                                / ( gde3w(:,:,jk) - gde3w(:,:,jk-1) ) 
    252 !!gm  use of e3t(:,:,:,Kmm) just above? 
    253       END DO 
     272      DO_3D_00_00( 2, jpkm1 )       ! complete with the level-dependent part 
     273         zemx_iwm(ji,jj,jk) = zemx_iwm(ji,jj,jk) + zweight(ji,jj,jk) * zfact(ji,jj) * wmask(ji,jj,jk)   & 
     274            &                                                        / ( gde3w(ji,jj,jk) - gde3w(ji,jj,jk-1) ) 
     275!!gm  use of e3t(ji,jj,:,Kmm) just above? 
     276      END_3D 
    254277      ! 
    255278!!gm  this is to be replaced by just a constant value znu=1.e-6 m2/s 
    256279      ! Calculate molecular kinematic viscosity 
    257       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)  & 
    258          &                                  + 0.02305_wp * ts(:,:,:,jp_sal,Kmm)  ) * tmask(:,:,:) * r1_rho0 
    259       DO jk = 2, jpkm1 
    260          znu_w(:,:,jk) = 0.5_wp * ( znu_t(:,:,jk-1) + znu_t(:,:,jk) ) * wmask(:,:,jk) 
    261       END DO 
     280      DO_3D_00_00( 1, jpkm1 ) 
     281         znu_t(ji,jj,jk) = 1.e-4_wp * (  17.91_wp - 0.53810_wp * ts(ji,jj,jk,jp_tem,Kmm)   & 
     282            &                                     + 0.00694_wp * ts(ji,jj,jk,jp_tem,Kmm) * ts(ji,jj,jk,jp_tem,Kmm)  & 
     283            &                                     + 0.02305_wp * ts(ji,jj,jk,jp_sal,Kmm)  ) * tmask(ji,jj,jk) * r1_rho0 
     284      END_3D 
     285      DO_3D_00_00( 2, jpkm1 ) 
     286         znu_w(ji,jj,jk) = 0.5_wp * ( znu_t(ji,jj,jk-1) + znu_t(ji,jj,jk) ) * wmask(ji,jj,jk) 
     287      END_3D 
    262288!!gm end 
    263289      ! 
    264290      ! Calculate turbulence intensity parameter Reb 
    265       DO jk = 2, jpkm1 
    266          zReb(:,:,jk) = zemx_iwm(:,:,jk) / MAX( 1.e-20_wp, znu_w(:,:,jk) * rn2(:,:,jk) ) 
    267       END DO 
     291      DO_3D_00_00( 2, jpkm1 ) 
     292         zReb(ji,jj,jk) = zemx_iwm(ji,jj,jk) / MAX( 1.e-20_wp, znu_w(ji,jj,jk) * rn2(ji,jj,jk) ) 
     293      END_3D 
    268294      ! 
    269295      ! Define internal wave-induced diffusivity 
    270       DO jk = 2, jpkm1 
    271          zav_wave(:,:,jk) = znu_w(:,:,jk) * zReb(:,:,jk) * r1_6   ! This corresponds to a constant mixing efficiency of 1/6 
    272       END DO 
     296      DO_3D_00_00( 2, jpkm1 ) 
     297         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 
     298      END_3D 
    273299      ! 
    274300      IF( ln_mevar ) THEN              ! Variable mixing efficiency case : modify zav_wave in the 
    275          DO_3D_11_11( 2, jpkm1 ) 
     301         DO_3D_00_00( 2, jpkm1 ) 
    276302            IF( zReb(ji,jj,jk) > 480.00_wp ) THEN 
    277303               zav_wave(ji,jj,jk) = 3.6515_wp * znu_w(ji,jj,jk) * SQRT( zReb(ji,jj,jk) ) 
     
    282308      ENDIF 
    283309      ! 
    284       DO jk = 2, jpkm1                 ! Bound diffusivity by molecular value and 100 cm2/s 
    285          zav_wave(:,:,jk) = MIN(  MAX( 1.4e-7_wp, zav_wave(:,:,jk) ), 1.e-2_wp  ) * wmask(:,:,jk) 
    286       END DO 
     310      DO_3D_00_00( 2, jpkm1 )          ! Bound diffusivity by molecular value and 100 cm2/s 
     311         zav_wave(ji,jj,jk) = MIN(  MAX( 1.4e-7_wp, zav_wave(ji,jj,jk) ), 1.e-2_wp  ) * wmask(ji,jj,jk) 
     312      END_3D 
    287313      ! 
    288314      IF( kt == nit000 ) THEN        !* Control print at first time-step: diagnose the energy consumed by zav_wave 
    289315         zztmp = 0._wp 
    290316!!gm used of glosum 3D.... 
    291          DO_3D_11_11( 2, jpkm1 ) 
     317         DO_3D_00_00( 2, jpkm1 ) 
    292318            zztmp = zztmp + e3w(ji,jj,jk,Kmm) * e1e2t(ji,jj)   & 
    293319               &          * MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj) 
     
    311337      IF( ln_tsdiff ) THEN          !* Option for differential mixing of salinity and temperature 
    312338         ztmp1 = 0.505_wp + 0.495_wp * TANH( 0.92_wp * ( LOG10( 1.e-20_wp ) - 0.60_wp ) ) 
    313          DO_3D_11_11( 2, jpkm1 ) 
     339         DO_3D_00_00( 2, jpkm1 ) 
    314340            ztmp2 = zReb(ji,jj,jk) * 5._wp * r1_6 
    315341            IF ( ztmp2 > 1.e-20_wp .AND. wmask(ji,jj,jk) == 1._wp ) THEN 
     
    320346         END_3D 
    321347         CALL iom_put( "av_ratio", zav_ratio ) 
    322          DO jk = 2, jpkm1           !* update momentum & tracer diffusivity with wave-driven mixing 
    323             p_avs(:,:,jk) = p_avs(:,:,jk) + zav_wave(:,:,jk) * zav_ratio(:,:,jk) 
    324             p_avt(:,:,jk) = p_avt(:,:,jk) + zav_wave(:,:,jk) 
    325             p_avm(:,:,jk) = p_avm(:,:,jk) + zav_wave(:,:,jk) 
    326          END DO 
     348         DO_3D_00_00( 2, jpkm1 )    !* update momentum & tracer diffusivity with wave-driven mixing 
     349            p_avs(ji,jj,jk) = p_avs(ji,jj,jk) + zav_wave(ji,jj,jk) * zav_ratio(ji,jj,jk) 
     350            p_avt(ji,jj,jk) = p_avt(ji,jj,jk) + zav_wave(ji,jj,jk) 
     351            p_avm(ji,jj,jk) = p_avm(ji,jj,jk) + zav_wave(ji,jj,jk) 
     352         END_3D 
    327353         ! 
    328354      ELSE                          !* update momentum & tracer diffusivity with wave-driven mixing 
    329          DO jk = 2, jpkm1 
    330             p_avs(:,:,jk) = p_avs(:,:,jk) + zav_wave(:,:,jk) 
    331             p_avt(:,:,jk) = p_avt(:,:,jk) + zav_wave(:,:,jk) 
    332             p_avm(:,:,jk) = p_avm(:,:,jk) + zav_wave(:,:,jk) 
    333          END DO 
     355         DO_3D_00_00( 2, jpkm1 ) 
     356            p_avs(ji,jj,jk) = p_avs(ji,jj,jk) + zav_wave(ji,jj,jk) 
     357            p_avt(ji,jj,jk) = p_avt(ji,jj,jk) + zav_wave(ji,jj,jk) 
     358            p_avm(ji,jj,jk) = p_avm(ji,jj,jk) + zav_wave(ji,jj,jk) 
     359         END_3D 
    334360      ENDIF 
    335361 
     
    341367      IF( iom_use("bflx_iwm") .OR. iom_use("pcmap_iwm") ) THEN 
    342368         ALLOCATE( z2d(jpi,jpj) , z3d(jpi,jpj,jpk) ) 
    343          z3d(:,:,:) = MAX( 0._wp, rn2(:,:,:) ) * zav_wave(:,:,:) 
    344          z2d(:,:) = 0._wp 
    345          DO jk = 2, jpkm1 
    346             z2d(:,:) = z2d(:,:) + e3w(:,:,jk,Kmm) * z3d(:,:,jk) * wmask(:,:,jk) 
    347          END DO 
    348          z2d(:,:) = rho0 * z2d(:,:) 
    349          CALL iom_put( "bflx_iwm", z3d ) 
     369         ! Initialisation for iom_put 
     370         DO_2D_00_00 
     371            z3d(ji,jj,1) = 0._wp   ;   z3d(ji,jj,jpk) = 0._wp 
     372         END_2D 
     373         z3d(           1:nn_hls,:,:) = 0._wp   ;   z3d(:,           1:nn_hls,:) = 0._wp 
     374         z3d(jpi-nn_hls+1:jpi   ,:,:) = 0._wp   ;   z3d(:,jpj-nn_hls+1:   jpj,:) = 0._wp 
     375         z2d(           1:nn_hls,:  ) = 0._wp   ;   z2d(:,           1:nn_hls  ) = 0._wp 
     376         z2d(jpi-nn_hls+1:jpi   ,:  ) = 0._wp   ;   z2d(:,jpj-nn_hls+1:   jpj  ) = 0._wp 
     377 
     378         DO_3D_00_00( 2, jpkm1 ) 
     379            z3d(ji,jj,jk) = MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk) 
     380         END_3D 
     381         DO_2D_00_00 
     382            z2d(ji,jj) = 0._wp 
     383         END_2D 
     384         DO_3D_00_00( 2, jpkm1 )  
     385            z2d(ji,jj) = z2d(ji,jj) + e3w(ji,jj,jk,Kmm) * z3d(ji,jj,jk) * wmask(ji,jj,jk) 
     386         END_3D 
     387         DO_2D_00_00 
     388            z2d(ji,jj) = rho0 * z2d(ji,jj) 
     389         END_2D 
     390         CALL iom_put(  "bflx_iwm", z3d ) 
    350391         CALL iom_put( "pcmap_iwm", z2d ) 
    351392         DEALLOCATE( z2d , z3d ) 
Note: See TracChangeset for help on using the changeset viewer.