Changeset 14665
- Timestamp:
- 2021-03-31T19:05:59+02:00 (3 years ago)
- Location:
- NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE
- Files:
-
- 14 edited
-
TRA/eosbn2.F90 (modified) (4 diffs)
-
ZDF/zdfddm.F90 (modified) (5 diffs)
-
ZDF/zdfdrg.F90 (modified) (4 diffs)
-
ZDF/zdfevd.F90 (modified) (6 diffs)
-
ZDF/zdfgls.F90 (modified) (20 diffs)
-
ZDF/zdfiwm.F90 (modified) (6 diffs)
-
ZDF/zdfmfc.F90 (modified) (7 diffs)
-
ZDF/zdfmxl.F90 (modified) (5 diffs)
-
ZDF/zdfosm.F90 (modified) (4 diffs)
-
ZDF/zdfphy.F90 (modified) (7 diffs)
-
ZDF/zdfric.F90 (modified) (3 diffs)
-
ZDF/zdfsh2.F90 (modified) (1 diff)
-
ZDF/zdfswm.F90 (modified) (1 diff)
-
ZDF/zdftke.F90 (modified) (27 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/eosbn2.F90
r14131 r14665 577 577 578 578 SUBROUTINE eos_insitu_pot_2d( pts, prhop ) 579 !! 580 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pts ! 1 : potential temperature [Celsius] 581 ! ! 2 : salinity [psu] 582 REAL(wp), DIMENSION(:,:) , INTENT( out) :: prhop ! potential density (surface referenced) 583 !! 584 CALL eos_insitu_pot_2d_t( pts, is_tile(pts), prhop, is_tile(prhop) ) 585 END SUBROUTINE eos_insitu_pot_2d 586 587 588 SUBROUTINE eos_insitu_pot_2d_t( pts, ktts, prhop, ktrhop ) 579 589 !!---------------------------------------------------------------------- 580 590 !! *** ROUTINE eos_insitu_pot *** … … 589 599 !! 590 600 !!---------------------------------------------------------------------- 591 REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in ) :: pts ! 1 : potential temperature [Celsius] 601 INTEGER , INTENT(in ) :: ktts, ktrhop 602 REAL(wp), DIMENSION(A2D_T(ktts),JPTS), INTENT(in ) :: pts ! 1 : potential temperature [Celsius] 592 603 ! ! 2 : salinity [psu] 593 REAL(wp), DIMENSION( jpi,jpj), INTENT( out) :: prhop ! potential density (surface referenced)604 REAL(wp), DIMENSION(A2D_T(ktrhop) ), INTENT( out) :: prhop ! potential density (surface referenced) 594 605 ! 595 606 INTEGER :: ji, jj, jk, jsmp ! dummy loop indices … … 606 617 CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==! 607 618 ! 608 DO_2D( 1, 1, 1, 1)609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 619 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 620 ! 621 zt = pts (ji,jj,jp_tem) * r1_T0 ! temperature 622 zs = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity 623 ztm = tmask(ji,jj,1) ! tmask 624 ! 625 zn0 = (((((EOS060*zt & 626 & + EOS150*zs+EOS050)*zt & 627 & + (EOS240*zs+EOS140)*zs+EOS040)*zt & 628 & + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt & 629 & + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt & 630 & + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt & 631 & + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 632 ! 633 ! 634 prhop(ji,jj) = zn0 * ztm ! potential density referenced at the surface 635 ! 636 END_2D 626 637 627 638 CASE( np_seos ) !== simplified EOS ==! 628 639 ! 629 DO_2D( 1, 1, 1, 1)640 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 630 641 zt = pts (ji,jj,jp_tem) - 10._wp 631 642 zs = pts (ji,jj,jp_sal) - 35._wp … … 646 657 IF( ln_timing ) CALL timing_stop('eos-pot') 647 658 ! 648 END SUBROUTINE eos_insitu_pot_2d 659 END SUBROUTINE eos_insitu_pot_2d_t 649 660 650 661 -
NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfddm.F90
r14636 r14665 83 83 REAL(dp) :: zavfs ! - - 84 84 REAL(wp) :: zavdt, zavds ! - - 85 REAL(wp), DIMENSION( jpi,jpj) :: zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd385 REAL(wp), DIMENSION(A2D(nn_hls)) :: zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3 86 86 !!---------------------------------------------------------------------- 87 87 ! … … 96 96 97 97 ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 ) !== R=zrau = (alpha / beta) (dk[t] / dk[s]) ==! 98 DO_2D( nn_hls , nn_hls, nn_hls, nn_hls) !== R=zrau = (alpha / beta) (dk[t] / dk[s]) ==!98 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) !== R=zrau = (alpha / beta) (dk[t] / dk[s]) ==! 99 99 zrw = ( gdepw(ji,jj,jk ,Kmm) - gdept(ji,jj,jk,Kmm) ) & 100 100 !!gm please, use e3w at Kmm below … … 113 113 114 114 ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 ) !== indicators ==! 115 DO_2D( nn_hls , nn_hls, nn_hls, nn_hls) !== indicators ==!115 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) !== indicators ==! 116 116 ! stability indicator: msks=1 if rn2>0; 0 elsewhere 117 117 IF( rn2(ji,jj,jk) + 1.e-12 <= 0. ) THEN ; zmsks(ji,jj) = 0._wp … … 137 137 END_2D 138 138 ! mask zmsk in order to have avt and avs masked 139 zmsks(:,:) = zmsks(:,:) * wmask( :,:,jk)139 zmsks(:,:) = zmsks(:,:) * wmask(A2D(nn_hls),jk) 140 140 141 141 … … 144 144 ! Constant eddy coefficient: reset to the background value 145 145 ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 ) 146 DO_2D ( nn_hls, nn_hls, nn_hls, nn_hls)146 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 147 147 zinr = 1._wp / zrau(ji,jj) 148 148 ! salt fingering -
NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfdrg.F90
r14636 r14665 118 118 IF( l_log_not_linssh ) THEN !== "log layer" ==! compute Cd and -Cd*|U| 119 119 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 120 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )120 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 121 121 imk = k_mk(ji,jj) ! ocean bottom level at t-points 122 122 zut = uu(ji,jj,imk,Kmm) + uu(ji-1,jj,imk,Kmm) ! 2 x velocity at t-point … … 131 131 ELSE !== standard Cd ==! 132 132 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 133 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )133 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 134 134 imk = k_mk(ji,jj) ! ocean bottom level at t-points 135 135 zut = uu(ji,jj,imk,Kmm) + uu(ji-1,jj,imk,Kmm) ! 2 x velocity at t-point … … 179 179 180 180 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 181 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )181 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 182 182 ikbu = mbku(ji,jj) ! deepest wet ocean u- & v-levels 183 183 ikbv = mbkv(ji,jj) … … 193 193 IF( ln_isfcav ) THEN ! ocean cavities 194 194 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 195 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )195 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 196 196 ikbu = miku(ji,jj) ! first wet ocean u- & v-levels 197 197 ikbv = mikv(ji,jj) -
NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfevd.F90
r14636 r14665 62 62 ! 63 63 INTEGER :: ji, jj, jk ! dummy loop indices 64 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zavt_evd, zavm_evd 64 ! NOTE: [tiling] use a SAVE array to store diagnostics, then send after all tiles are finished. This is necessary because p_avt/p_avm are modified on adjacent tiles when using nn_hls > 1. zavt_evd/zavm_evd are then zero on some points when subsequently calculated for these tiles. 65 REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: zavt_evd, zavm_evd 65 66 !!---------------------------------------------------------------------- 66 67 ! 67 IF( kt == nit000 ) THEN 68 IF(lwp) WRITE(numout,*) 69 IF(lwp) WRITE(numout,*) 'zdf_evd : Enhanced Vertical Diffusion (evd)' 70 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 71 IF(lwp) WRITE(numout,*) 68 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 69 IF( kt == nit000 ) THEN 70 IF(lwp) WRITE(numout,*) 71 IF(lwp) WRITE(numout,*) 'zdf_evd : Enhanced Vertical Diffusion (evd)' 72 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 73 IF(lwp) WRITE(numout,*) 74 ENDIF 75 76 ALLOCATE( zavt_evd(jpi,jpj,jpk) ) 77 IF( nn_evdm == 1 ) ALLOCATE( zavm_evd(jpi,jpj,jpk) ) 72 78 ENDIF 73 79 ! 74 80 ! 75 zavt_evd(:,:,:) = p_avt(:,:,:) ! set avt prior to evd application 81 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk ) 82 zavt_evd(ji,jj,jk) = p_avt(ji,jj,jk) ! set avt prior to evd application 83 END_3D 76 84 ! 77 85 SELECT CASE ( nn_evdm ) … … 79 87 CASE ( 1 ) !== enhance tracer & momentum Kz ==! (if rn2<-1.e-12) 80 88 ! 81 zavm_evd(:,:,:) = p_avm(:,:,:) ! set avm prior to evd application 89 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk ) 90 zavm_evd(ji,jj,jk) = p_avm(ji,jj,jk) ! set avm prior to evd application 91 END_3D 82 92 ! 83 93 !! change last digits results … … 88 98 ! 89 99 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 90 DO_3D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )100 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 91 101 IF( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) THEN 92 102 p_avt(ji,jj,jk) = rn_evd * wmask(ji,jj,jk) … … 95 105 END_3D 96 106 ! 97 zavm_evd(:,:,:) = p_avm(:,:,:) - zavm_evd(:,:,:) ! change in avm due to evd 98 CALL iom_put( "avm_evd", zavm_evd ) ! output this change 107 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk ) 108 zavm_evd(ji,jj,jk) = p_avm(ji,jj,jk) - zavm_evd(ji,jj,jk) ! change in avm due to evd 109 END_3D 110 IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN ! Do only on the last tile 111 CALL iom_put( "avm_evd", zavm_evd ) ! output this change 112 DEALLOCATE( zavm_evd ) 113 ENDIF 99 114 ! 100 115 CASE DEFAULT !== enhance tracer Kz ==! (if rn2<-1.e-12) … … 105 120 106 121 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 107 DO_3D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )122 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 108 123 IF( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) & 109 124 p_avt(ji,jj,jk) = rn_evd * wmask(ji,jj,jk) … … 112 127 END SELECT 113 128 ! 114 zavt_evd(:,:,:) = p_avt(:,:,:) - zavt_evd(:,:,:) ! change in avt due to evd 115 CALL iom_put( "avt_evd", zavt_evd ) ! output this change 129 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk ) 130 zavt_evd(ji,jj,jk) = p_avt(ji,jj,jk) - zavt_evd(ji,jj,jk) ! change in avt due to evd 131 END_3D 132 IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN ! Do only on the last tile 133 CALL iom_put( "avt_evd", zavt_evd ) ! output this change 134 DEALLOCATE( zavt_evd ) 135 ENDIF 116 136 IF( l_trdtra ) CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_evd, zavt_evd ) 117 137 ! -
NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfgls.F90
r14636 r14665 137 137 USE zdf_oce , ONLY : en, avtb, avmb ! ocean vertical physics 138 138 !! 139 INTEGER , INTENT(in ) :: kt ! ocean time step140 INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices141 REAL(wp), DIMENSION( :,:,:), INTENT(in ) :: p_sh2 ! shear production term142 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: p_avm, p_avt ! momentum and tracer Kz (w-points)139 INTEGER , INTENT(in ) :: kt ! ocean time step 140 INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices 141 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in ) :: p_sh2 ! shear production term 142 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: p_avm, p_avt ! momentum and tracer Kz (w-points) 143 143 ! 144 144 INTEGER :: ji, jj, jk ! dummy loop arguments … … 151 151 REAL(wp) :: gh, gm, shr, dif, zsqen, zavt, zavm ! - - 152 152 REAL(wp) :: zmsku, zmskv ! - - 153 REAL(wp), DIMENSION( jpi,jpj) :: zdep154 REAL(wp), DIMENSION( jpi,jpj) :: zkar155 REAL(wp), DIMENSION( jpi,jpj) :: zflxs! Turbulence fluxed induced by internal waves156 REAL(wp), DIMENSION( jpi,jpj) :: zhsro! Surface roughness (surface waves)157 REAL(wp), DIMENSION( jpi,jpj) :: zice_fra! Tapering of wave breaking under sea ice158 REAL(wp), DIMENSION( jpi,jpj,jpk) :: eb! tke at time before159 REAL(wp), DIMENSION( jpi,jpj,jpk) :: hmxl_b! mixing length at time before160 REAL(wp), DIMENSION( jpi,jpj,jpk) :: eps! dissipation rate161 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zwall_psi! Wall function use in the wb case (ln_sigpsi)162 REAL(wp), DIMENSION( jpi,jpj,jpk) :: psi! psi at time now163 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zd_lw, zd_up, zdiag ! lower, upper and diagonal of the matrix164 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zstt, zstm! stability function on tracer and momentum153 REAL(wp), DIMENSION(A2D(nn_hls)) :: zdep 154 REAL(wp), DIMENSION(A2D(nn_hls)) :: zkar 155 REAL(wp), DIMENSION(A2D(nn_hls)) :: zflxs ! Turbulence fluxed induced by internal waves 156 REAL(wp), DIMENSION(A2D(nn_hls)) :: zhsro ! Surface roughness (surface waves) 157 REAL(wp), DIMENSION(A2D(nn_hls)) :: zice_fra ! Tapering of wave breaking under sea ice 158 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: eb ! tke at time before 159 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: hmxl_b ! mixing length at time before 160 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: eps ! dissipation rate 161 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zwall_psi ! Wall function use in the wb case (ln_sigpsi) 162 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: psi ! psi at time now 163 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zd_lw, zd_up, zdiag ! lower, upper and diagonal of the matrix 164 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zstt, zstm ! stability function on tracer and momentum 165 165 !!-------------------------------------------------------------------- 166 166 ! 167 167 ! Preliminary computing 168 169 ustar2_surf(:,:) = 0._wp ; psi(:,:,:) = 0._wp 170 ustar2_top (:,:) = 0._wp ; zwall_psi(:,:,:) = 0._wp 171 ustar2_bot (:,:) = 0._wp 168 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 169 ustar2_surf(ji,jj) = 0._wp ; ustar2_top(ji,jj) = 0._wp ; ustar2_bot(ji,jj) = 0._wp 170 END_2D 171 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk ) 172 psi(ji,jj,jk) = 0._wp ; zwall_psi(ji,jj,jk) = 0._wp 173 END_3D 172 174 173 175 SELECT CASE ( nn_z0_ice ) 174 176 CASE( 0 ) ; zice_fra(:,:) = 0._wp 175 CASE( 1 ) ; zice_fra(:,:) = TANH( fr_i( :,:) * 10._wp )176 CASE( 2 ) ; zice_fra(:,:) = fr_i( :,:)177 CASE( 3 ) ; zice_fra(:,:) = MIN( 4._wp * fr_i( :,:) , 1._wp )177 CASE( 1 ) ; zice_fra(:,:) = TANH( fr_i(A2D(nn_hls)) * 10._wp ) 178 CASE( 2 ) ; zice_fra(:,:) = fr_i(A2D(nn_hls)) 179 CASE( 3 ) ; zice_fra(:,:) = MIN( 4._wp * fr_i(A2D(nn_hls)) , 1._wp ) 178 180 END SELECT 179 181 180 182 ! Compute surface, top and bottom friction at T-points 181 183 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) !== surface ocean friction 182 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) !== surface ocean friction184 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) !== surface ocean friction 183 185 ustar2_surf(ji,jj) = r1_rho0 * taum(ji,jj) * tmask(ji,jj,1) ! surface friction 184 186 END_2D … … 188 190 IF( .NOT.ln_drg_OFF ) THEN !== top/bottom friction (explicit before friction) 189 191 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) ! bottom friction (explicit before friction) 190 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! bottom friction (explicit before friction)192 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! bottom friction (explicit before friction) 191 193 zmsku = 0.5_wp * ( 2._wp - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 192 194 zmskv = 0.5_wp * ( 2._wp - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) ! (CAUTION: CdU<0) … … 196 198 IF( ln_isfcav ) THEN 197 199 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) ! top friction 198 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! top friction200 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! top friction 199 201 zmsku = 0.5_wp * ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) ) 200 202 zmskv = 0.5_wp * ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) ) ! (CAUTION: CdU<0) … … 209 211 zhsro(:,:) = rn_hsro 210 212 CASE ( 1 ) ! Standard Charnock formula 211 zhsro(:,:) = MAX( rsbc_zs1 * ustar2_surf( :,:) , rn_hsro )213 zhsro(:,:) = MAX( rsbc_zs1 * ustar2_surf(A2D(nn_hls)) , rn_hsro ) 212 214 CASE ( 2 ) ! Roughness formulae according to Rascle et al., Ocean Modelling (2008) 213 215 !!gm faster coding : the 2 comment lines should be used 214 216 !!gm zcof = 2._wp * 0.6_wp / 28._wp 215 217 !!gm zdep(:,:) = 30._wp * TANH( zcof/ SQRT( MAX(ustar2_surf(:,:),rsmall) ) ) ! Wave age (eq. 10) 216 zdep (:,:) = 30.*TANH( 2.*0.3/(28.*SQRT(MAX(ustar2_surf(:,:),rsmall))) ) ! Wave age (eq. 10) 217 zhsro(:,:) = MAX(rsbc_zs2 * ustar2_surf(:,:) * zdep(:,:)**1.5, rn_hsro) ! zhsro = rn_frac_hs * Hsw (eq. 11) 218 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 219 zcof = 30.*TANH( 2.*0.3/(28.*SQRT(MAX(ustar2_surf(ji,jj),rsmall))) ) ! Wave age (eq. 10) 220 zhsro(ji,jj) = MAX(rsbc_zs2 * ustar2_surf(ji,jj) * zcof**1.5, rn_hsro) ! zhsro = rn_frac_hs * Hsw (eq. 11) 221 END_2D 218 222 CASE ( 3 ) ! Roughness given by the wave model (coupled or read in file) 219 zhsro(:,:) = MAX(rn_frac_hs * hsw( :,:), rn_hsro) ! (rn_frac_hs=1.6 see Eq. (5) of Rascle et al. 2008 )223 zhsro(:,:) = MAX(rn_frac_hs * hsw(A2D(nn_hls)), rn_hsro) ! (rn_frac_hs=1.6 see Eq. (5) of Rascle et al. 2008 ) 220 224 END SELECT 221 225 ! 222 226 ! adapt roughness where there is sea ice 223 zhsro(:,:) = ( (1._wp-zice_fra(:,:)) * zhsro(:,:) + zice_fra(:,:) * rn_hsri )*tmask(:,:,1) + (1._wp - tmask(:,:,1))*rn_hsro 227 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 228 zhsro(ji,jj) = ( (1._wp-zice_fra(ji,jj)) * zhsro(ji,jj) + zice_fra(ji,jj) * rn_hsri )*tmask(ji,jj,1) + & 229 & (1._wp - tmask(ji,jj,1))*rn_hsro 230 END_2D 224 231 ! 225 232 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 ) !== Compute dissipation rate ==! … … 229 236 230 237 ! Save tke at before time step 231 eb (:,:,:) = en (:,:,:) 232 hmxl_b(:,:,:) = hmxl_n(:,:,:) 238 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk ) 239 eb (ji,jj,jk) = en (ji,jj,jk) 240 hmxl_b(ji,jj,jk) = hmxl_n(ji,jj,jk) 241 END_3D 233 242 234 243 IF( nn_clos == 0 ) THEN ! Mellor-Yamada 235 244 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 236 DO_3D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )245 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 237 246 zup = hmxl_n(ji,jj,jk) * gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) 238 247 zdown = vkarmn * gdepw(ji,jj,jk,Kmm) * ( -gdepw(ji,jj,jk,Kmm) + gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) ) … … 256 265 257 266 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 258 DO_3D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )267 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 259 268 ! 260 269 buoy = - p_avt(ji,jj,jk) * rn2(ji,jj,jk) ! stratif. destruction … … 309 318 ! 310 319 CASE ( 0 ) ! Dirichlet boundary condition (set e at k=1 & 2) 311 ! First level 312 en (:,:,1) = MAX( rn_emin , rc02r * ustar2_surf(:,:) * (1._wp + (1._wp-zice_fra(:,:))*rsbc_tke1)**r2_3 ) 313 zd_lw(:,:,1) = en(:,:,1) 314 zd_up(:,:,1) = 0._wp 315 zdiag(:,:,1) = 1._wp 316 ! 317 ! One level below 318 en (:,:,2) = MAX( rc02r * ustar2_surf(:,:) * ( 1._wp + (1._wp-zice_fra(:,:))*rsbc_tke1 * ((zhsro(:,:)+gdepw(:,:,2,Kmm)) & 319 & / zhsro(:,:) )**(1.5_wp*ra_sf) )**(2._wp/3._wp) , rn_emin ) 320 zd_lw(:,:,2) = 0._wp 321 zd_up(:,:,2) = 0._wp 322 zdiag(:,:,2) = 1._wp 323 ! 324 ! 320 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 321 ! First level 322 en (ji,jj,1) = MAX( rn_emin , rc02r * ustar2_surf(ji,jj) * (1._wp + (1._wp-zice_fra(ji,jj))*rsbc_tke1)**r2_3 ) 323 zd_lw(ji,jj,1) = en(ji,jj,1) 324 zd_up(ji,jj,1) = 0._wp 325 zdiag(ji,jj,1) = 1._wp 326 ! 327 ! One level below 328 en (ji,jj,2) = MAX( rn_emin , rc02r * ustar2_surf(ji,jj) * (1._wp + (1._wp-zice_fra(ji,jj))*rsbc_tke1 & 329 & * ((zhsro(ji,jj)+gdepw(ji,jj,2,Kmm)) / zhsro(ji,jj) )**(1.5_wp*ra_sf) )**r2_3 ) 330 zd_lw(ji,jj,2) = 0._wp 331 zd_up(ji,jj,2) = 0._wp 332 zdiag(ji,jj,2) = 1._wp 333 END_2D 334 ! 335 ! 325 336 CASE ( 1 ) ! Neumann boundary condition (set d(e)/dz) 326 !327 ! Dirichlet conditions at k=1328 en (:,:,1) = MAX( rc02r * ustar2_surf(:,:) * (1._wp + (1._wp-zice_fra(:,:))*rsbc_tke1)**r2_3 , rn_emin )329 zd_lw(:,:,1) = en(:,:,1)330 zd_up(:,:,1) = 0._wp331 zdiag(:,:,1) = 1._wp332 !333 ! at k=2, set de/dz=Fw334 !cbr335 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) ! zdiag zd_lw not defined/used on the halo336 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )! zdiag zd_lw not defined/used on the halo337 zdiag(ji,jj,2) = zdiag(ji,jj,2) + zd_lw(ji,jj,2) ! Remove zd_lw from zdiag338 zd_lw(ji,jj,2) = 0._wp339 END_2D340 zkar (:,:) = (rl_sf + (vkarmn-rl_sf)*(1.-EXP(-rtrans*gdept(:,:,1,Kmm)/zhsro(:,:)) ))341 zflxs(:,:) = rsbc_tke2 * (1._wp-zice_fra(:,:)) * ustar2_surf(:,:)**1.5_wp * zkar(:,:) &342 & * ( ( zhsro(:,:)+gdept(:,:,1,Kmm) ) / zhsro(:,:) )**(1.5_wp*ra_sf)337 ! 338 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 339 ! Dirichlet conditions at k=1 340 en (ji,jj,1) = MAX( rn_emin , rc02r * ustar2_surf(ji,jj) * (1._wp + (1._wp-zice_fra(ji,jj))*rsbc_tke1)**r2_3 ) 341 zd_lw(ji,jj,1) = en(ji,jj,1) 342 zd_up(ji,jj,1) = 0._wp 343 zdiag(ji,jj,1) = 1._wp 344 ! 345 ! at k=2, set de/dz=Fw 346 !cbr 347 ! zdiag zd_lw not defined/used on the halo 348 zdiag(ji,jj,2) = zdiag(ji,jj,2) + zd_lw(ji,jj,2) ! Remove zd_lw from zdiag 349 zd_lw(ji,jj,2) = 0._wp 350 ! 351 zkar (ji,jj) = (rl_sf + (vkarmn-rl_sf)*(1.-EXP(-rtrans*gdept(ji,jj,1,Kmm)/zhsro(ji,jj)) )) 352 zflxs(ji,jj) = rsbc_tke2 * (1._wp-zice_fra(ji,jj)) * ustar2_surf(ji,jj)**1.5_wp * zkar(ji,jj) & 353 & * ( ( zhsro(ji,jj)+gdept(ji,jj,1,Kmm) ) / zhsro(ji,jj) )**(1.5_wp*ra_sf) 343 354 !!gm why not : * ( 1._wp + gdept(:,:,1,Kmm) / zhsro(:,:) )**(1.5_wp*ra_sf) 344 en(:,:,2) = en(:,:,2) + zflxs(:,:) / e3w(:,:,2,Kmm) 345 ! 346 ! 355 en(ji,jj,2) = en(ji,jj,2) + zflxs(ji,jj) / e3w(ji,jj,2,Kmm) 356 END_2D 357 ! 358 ! 347 359 END SELECT 348 360 … … 356 368 ! ! Balance between the production and the dissipation terms 357 369 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 358 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )370 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 359 371 !!gm This means that bottom and ocean w-level above have a specified "en" value. Sure ???? 360 372 !! With thick deep ocean level thickness, this may be quite large, no ??? … … 373 385 END_2D 374 386 ! 387 ! NOTE: [tiling] not tested- ctl_stop with ln_isfcav! 375 388 IF( ln_isfcav) THEN ! top boundary (ocean cavity) 376 389 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 377 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )390 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 378 391 itop = mikt(ji,jj) ! k top w-point 379 392 itopp1 = mikt(ji,jj) + 1 ! k+1 1st w-point below the top one … … 394 407 ! 395 408 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 396 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )409 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 397 410 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 398 411 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 … … 408 421 en (ji,jj,ibot) = z_en 409 422 END_2D 423 ! NOTE: [tiling] not tested- ctl_stop with ln_isfcav! 410 424 IF( ln_isfcav) THEN ! top boundary (ocean cavity) 411 425 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 412 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )426 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 413 427 itop = mikt(ji,jj) ! k top w-point 414 428 itopp1 = mikt(ji,jj) + 1 ! k+1 1st w-point below the top one … … 431 445 ! ---------------------------------------------------------- 432 446 ! 433 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 447 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 434 448 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 435 449 zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) … … 440 454 END_3D 441 455 ! [comm_cleanup] ! DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 ) ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 442 DO_3DS ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpkm1, 2, -1 ) ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk456 DO_3DS_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpkm1, 2, -1 ) ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 443 457 en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 444 458 END_3D 445 459 ! ! set the minimum value of tke 446 en(:,:,:) = MAX( en(:,:,:), rn_emin ) 460 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk ) 461 en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) 462 END_3D 447 463 448 464 !!----------------------------------------!! … … 535 551 CASE ( 0 ) ! Dirichlet boundary conditions 536 552 ! 537 ! Surface value 538 zdep (:,:) = zhsro(:,:) * rl_sf ! Cosmetic 539 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 540 zd_lw(:,:,1) = psi(:,:,1) 541 zd_up(:,:,1) = 0._wp 542 zdiag(:,:,1) = 1._wp 543 ! 544 ! One level below 545 zkar (:,:) = (rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdepw(:,:,2,Kmm)/zhsro(:,:) ))) 546 zdep (:,:) = (zhsro(:,:) + gdepw(:,:,2,Kmm)) * zkar(:,:) 547 psi (:,:,2) = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 548 zd_lw(:,:,2) = 0._wp 549 zd_up(:,:,2) = 0._wp 550 zdiag(:,:,2) = 1._wp 553 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 554 ! Surface value 555 zdep (ji,jj) = zhsro(ji,jj) * rl_sf ! Cosmetic 556 psi (ji,jj,1) = rc0**rpp * en(ji,jj,1)**rmm * zdep(ji,jj)**rnn * tmask(ji,jj,1) 557 zd_lw(ji,jj,1) = psi(ji,jj,1) 558 zd_up(ji,jj,1) = 0._wp 559 zdiag(ji,jj,1) = 1._wp 560 ! 561 ! One level below 562 zkar (ji,jj) = (rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdepw(ji,jj,2,Kmm)/zhsro(ji,jj) ))) 563 zdep (ji,jj) = (zhsro(ji,jj) + gdepw(ji,jj,2,Kmm)) * zkar(ji,jj) 564 psi (ji,jj,2) = rc0**rpp * en(ji,jj,2)**rmm * zdep(ji,jj)**rnn * tmask(ji,jj,1) 565 zd_lw(ji,jj,2) = 0._wp 566 zd_up(ji,jj,2) = 0._wp 567 zdiag(ji,jj,2) = 1._wp 568 END_2D 551 569 ! 552 570 CASE ( 1 ) ! Neumann boundary condition on d(psi)/dz 553 571 ! 554 ! Surface value: Dirichlet 555 zdep (:,:) = zhsro(:,:) * rl_sf 556 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 557 zd_lw(:,:,1) = psi(:,:,1) 558 zd_up(:,:,1) = 0._wp 559 zdiag(:,:,1) = 1._wp 560 ! 561 ! Neumann condition at k=2 562 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) ! zdiag zd_lw not defined/used on the halo 563 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! zdiag zd_lw not defined/used on the halo 572 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 573 ! Surface value: Dirichlet 574 zdep (ji,jj) = zhsro(ji,jj) * rl_sf 575 psi (ji,jj,1) = rc0**rpp * en(ji,jj,1)**rmm * zdep(ji,jj)**rnn * tmask(ji,jj,1) 576 zd_lw(ji,jj,1) = psi(ji,jj,1) 577 zd_up(ji,jj,1) = 0._wp 578 zdiag(ji,jj,1) = 1._wp 579 ! 580 ! Neumann condition at k=2, zdiag zd_lw not defined/used on the halo 564 581 zdiag(ji,jj,2) = zdiag(ji,jj,2) + zd_lw(ji,jj,2) ! Remove zd_lw from zdiag 565 582 zd_lw(ji,jj,2) = 0._wp 583 ! 584 ! Set psi vertical flux at the surface: 585 zkar (ji,jj) = rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdept(ji,jj,1,Kmm)/zhsro(ji,jj) )) ! Lengh scale slope 586 zdep (ji,jj) = ((zhsro(ji,jj) + gdept(ji,jj,1,Kmm)) / zhsro(ji,jj))**(rmm*ra_sf) 587 zflxs(ji,jj) = (rnn + (1._wp-zice_fra(ji,jj))*rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(ji,jj)) & 588 & *(1._wp + (1._wp-zice_fra(ji,jj))*rsbc_tke1*zdep(ji,jj))**(2._wp*rmm/3._wp-1_wp) 589 zdep (ji,jj) = rsbc_psi1 * (zwall_psi(ji,jj,1)*p_avm(ji,jj,1)+zwall_psi(ji,jj,2)*p_avm(ji,jj,2)) * & 590 & ustar2_surf(ji,jj)**rmm * zkar(ji,jj)**rnn * (zhsro(ji,jj) + gdept(ji,jj,1,Kmm))**(rnn-1.) 591 zflxs(ji,jj) = zdep(ji,jj) * zflxs(ji,jj) 592 psi (ji,jj,2) = psi(ji,jj,2) + zflxs(ji,jj) / e3w(ji,jj,2,Kmm) 566 593 END_2D 567 !568 ! Set psi vertical flux at the surface:569 zkar (:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdept(:,:,1,Kmm)/zhsro(:,:) )) ! Lengh scale slope570 zdep (:,:) = ((zhsro(:,:) + gdept(:,:,1,Kmm)) / zhsro(:,:))**(rmm*ra_sf)571 zflxs(:,:) = (rnn + (1._wp-zice_fra(:,:))*rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:)) &572 & *(1._wp + (1._wp-zice_fra(:,:))*rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp)573 zdep (:,:) = rsbc_psi1 * (zwall_psi(:,:,1)*p_avm(:,:,1)+zwall_psi(:,:,2)*p_avm(:,:,2)) * &574 & ustar2_surf(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + gdept(:,:,1,Kmm))**(rnn-1.)575 zflxs(:,:) = zdep(:,:) * zflxs(:,:)576 psi (:,:,2) = psi(:,:,2) + zflxs(:,:) / e3w(:,:,2,Kmm)577 594 ! 578 595 END SELECT … … 638 655 ! ---------------- 639 656 ! 640 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 657 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 641 658 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 642 659 zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) … … 688 705 ! -------------------------------------------------- 689 706 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! Note that this set boundary conditions on hmxl_n at the same time 690 DO_3D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) ! Note that this set boundary conditions on hmxl_n at the same time707 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) ! Note that this set boundary conditions on hmxl_n at the same time 691 708 ! limitation 692 709 eps (ji,jj,jk) = MAX( eps(ji,jj,jk), rn_epsmin ) 693 710 hmxl_n(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk) 694 ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated) 695 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 696 IF( ln_length_lim ) hmxl_n(ji,jj,jk) = MIN( rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), hmxl_n(ji,jj,jk) ) 697 END_3D 711 END_3D 712 IF( ln_length_lim ) THEN ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated) 713 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 714 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 715 hmxl_n(ji,jj,jk) = MIN( rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), hmxl_n(ji,jj,jk) ) 716 END_3D 717 ENDIF 698 718 699 719 ! … … 760 780 END_2D 761 781 762 zstt(:,:, 1) = wmask( :,:, 1) ! default value not needed but avoid a bug when looking for undefined values (-fpe0)763 zstt(:,:,jpk) = wmask( :,:,jpk) ! default value not needed but avoid a bug when looking for undefined values (-fpe0)782 zstt(:,:, 1) = wmask(A2D(nn_hls), 1) ! default value not needed but avoid a bug when looking for undefined values (-fpe0) 783 zstt(:,:,jpk) = wmask(A2D(nn_hls),jpk) ! default value not needed but avoid a bug when looking for undefined values (-fpe0) 764 784 765 785 !!gm should be done for ISF (top boundary cond.) … … 772 792 ! for zd_lw and zd_up but they have to be defined to avoid a bug when looking for undefined values (-fpe0) 773 793 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpk ) 774 DO_3D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk )794 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk ) 775 795 zsqen = SQRT( 2._wp * en(ji,jj,jk) ) * hmxl_n(ji,jj,jk) 776 796 zavt = zsqen * zstt(ji,jj,jk) … … 779 799 p_avm(ji,jj,jk) = MAX( zavm, avmb(jk) ) ! Note that avm is not masked at the surface and the bottom 780 800 END_3D 781 p_avt( :,:,1) = 0._wp801 p_avt(A2D(nn_hls),1) = 0._wp 782 802 ! 783 803 IF(sn_cfctl%l_prtctl) THEN -
NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfiwm.F90
r14636 r14665 125 125 ! 126 126 INTEGER :: ji, jj, jk ! dummy loop indices 127 REAL(wp) :: zztmp, ztmp1, ztmp2 ! scalar workspace 128 REAL(wp), DIMENSION(jpi,jpj) :: zfact ! Used for vertical structure 129 REAL(wp), DIMENSION(jpi,jpj) :: zhdep ! Ocean depth 130 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwkb ! WKB-stretched height above bottom 131 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zweight ! Weight for high mode vertical distribution 132 REAL(wp), DIMENSION(jpi,jpj,jpk) :: znu_t ! Molecular kinematic viscosity (T grid) 133 REAL(wp), DIMENSION(jpi,jpj,jpk) :: znu_w ! Molecular kinematic viscosity (W grid) 134 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zReb ! Turbulence intensity parameter 135 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zemx_iwm ! local energy density available for mixing (W/kg) 136 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zav_ratio ! S/T diffusivity ratio (only for ln_tsdiff=T) 137 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zav_wave ! Internal wave-induced diffusivity 127 ! NOTE: [tiling] declare zztmp with SAVE for mpp_sum 128 REAL(wp), SAVE :: zztmp 129 REAL(wp) :: ztmp1, ztmp2 ! scalar workspace 130 REAL(wp), DIMENSION(A2D(nn_hls)) :: zfact ! Used for vertical structure 131 REAL(wp), DIMENSION(A2D(nn_hls)) :: zhdep ! Ocean depth 132 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zwkb ! WKB-stretched height above bottom 133 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zweight ! Weight for high mode vertical distribution 134 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: znu_t ! Molecular kinematic viscosity (T grid) 135 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: znu_w ! Molecular kinematic viscosity (W grid) 136 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zReb ! Turbulence intensity parameter 137 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zemx_iwm ! local energy density available for mixing (W/kg) 138 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zav_ratio ! S/T diffusivity ratio (only for ln_tsdiff=T) 139 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zav_wave ! Internal wave-induced diffusivity 138 140 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z3d ! 3D workspace used for iom_put 139 141 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z2d ! 2D - - - - … … 337 339 ! 338 340 IF( kt == nit000 ) THEN !* Control print at first time-step: diagnose the energy consumed by zav_wave 339 zztmp = 0._wp341 IF( .NOT. l_istiled .OR. ntile == 1 ) zztmp = 0._wp ! Do only on the first tile 340 342 !!gm used of glosum 3D.... 341 343 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 ) … … 344 346 & * MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj) 345 347 END_3D 346 CALL mpp_sum( 'zdfiwm', zztmp ) 347 zztmp = rho0 * zztmp ! Global integral of rauo * Kz * N^2 = power contributing to mixing 348 ! 349 IF(lwp) THEN 350 WRITE(numout,*) 351 WRITE(numout,*) 'zdf_iwm : Internal wave-driven mixing (iwm)' 352 WRITE(numout,*) '~~~~~~~ ' 353 WRITE(numout,*) 354 WRITE(numout,*) ' Total power consumption by av_wave = ', zztmp * 1.e-12_wp, 'TW' 348 349 IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN ! Do only on the last tile 350 CALL mpp_sum( 'zdfiwm', zztmp ) 351 zztmp = rho0 * zztmp ! Global integral of rauo * Kz * N^2 = power contributing to mixing 352 ! 353 IF(lwp) THEN 354 WRITE(numout,*) 355 WRITE(numout,*) 'zdf_iwm : Internal wave-driven mixing (iwm)' 356 WRITE(numout,*) '~~~~~~~ ' 357 WRITE(numout,*) 358 WRITE(numout,*) ' Total power consumption by av_wave = ', zztmp * 1.e-12_wp, 'TW' 359 ENDIF 355 360 ENDIF 356 361 ENDIF … … 373 378 CALL iom_put( "av_ratio", zav_ratio ) 374 379 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 ) !* update momentum & tracer diffusivity with wave-driven mixing 375 DO_3D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) !* update momentum & tracer diffusivity with wave-driven mixing380 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) !* update momentum & tracer diffusivity with wave-driven mixing 376 381 p_avs(ji,jj,jk) = p_avs(ji,jj,jk) + zav_wave(ji,jj,jk) * zav_ratio(ji,jj,jk) 377 382 p_avt(ji,jj,jk) = p_avt(ji,jj,jk) + zav_wave(ji,jj,jk) … … 381 386 ELSE !* update momentum & tracer diffusivity with wave-driven mixing 382 387 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 383 DO_3D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )388 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 384 389 p_avs(ji,jj,jk) = p_avs(ji,jj,jk) + zav_wave(ji,jj,jk) 385 390 p_avt(ji,jj,jk) = p_avt(ji,jj,jk) + zav_wave(ji,jj,jk) … … 394 399 ! vertical integral of rho0 * Kz * N^2 , energy density (zemx_iwm) 395 400 IF( iom_use("bflx_iwm") .OR. iom_use("pcmap_iwm") ) THEN 396 ALLOCATE( z2d( jpi,jpj) , z3d(jpi,jpj,jpk) )401 ALLOCATE( z2d(A2D(nn_hls)) , z3d(A2D(nn_hls),jpk) ) 397 402 ! Initialisation for iom_put 398 403 DO_2D( 0, 0, 0, 0 ) 399 404 z3d(ji,jj,1) = 0._wp ; z3d(ji,jj,jpk) = 0._wp 400 405 END_2D 401 z3d( 1:nn_hls,:,:) = 0._wp ; z3d(:, 1:nn_hls,:) = 0._wp402 z3d(jpi-nn_hls+1:jpi ,:,:) = 0._wp ; z3d(:,jpj-nn_hls+1: jpj,:) = 0._wp403 z2d( 1:nn_hls,: ) = 0._wp ; z2d(:, 1:nn_hls ) = 0._wp404 z2d(jpi-nn_hls+1:jpi ,: ) = 0._wp ; z2d(:,jpj-nn_hls+1: jpj ) = 0._wp405 406 406 407 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) -
NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfmfc.F90
r14636 r14665 96 96 INTEGER , INTENT(in) :: Kmm, Krhs ! time level indices 97 97 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers and RHS of tracer equation 98 REAL(wp), DIMENSION( jpi,jpj,jpk,2) :: ztsp ! T/S of the plume99 REAL(wp), DIMENSION( jpi,jpj,jpk,2) :: ztse ! T/S at W point100 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zrwp !101 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zrwp2 !102 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zapp !103 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zedmf !104 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zepsT, zepsW !105 ! 106 REAL(wp), DIMENSION( jpi,jpj) :: zustar, zustar2 !107 REAL(wp), DIMENSION( jpi,jpj) :: zuws, zvws, zsws, zfnet !108 REAL(wp), DIMENSION( jpi,jpj) :: zfbuo, zrautbm1, zrautb, zraupl109 REAL(wp), DIMENSION( jpi,jpj) :: zwpsurf !110 REAL(wp), DIMENSION( jpi,jpj) :: zop0 , zsp0 !111 REAL(wp), DIMENSION( jpi,jpj) :: zrwp_0, zrwp2_0 !112 REAL(wp), DIMENSION( jpi,jpj) :: zapp0 !113 REAL(wp), DIMENSION( jpi,jpj) :: zphp, zph, zphpm1, zphm1, zNHydro114 REAL(wp), DIMENSION( jpi,jpj) :: zhcmo !115 ! 116 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zn2 ! N^2117 REAL(wp), DIMENSION( jpi,jpj,2 ) :: zab, zabm1, zabp ! alpha and beta98 REAL(wp), DIMENSION(A2D(nn_hls),jpk,2) :: ztsp ! T/S of the plume 99 REAL(wp), DIMENSION(A2D(nn_hls),jpk,2) :: ztse ! T/S at W point 100 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zrwp ! 101 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zrwp2 ! 102 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zapp ! 103 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zedmf ! 104 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zepsT, zepsW ! 105 ! 106 REAL(wp), DIMENSION(A2D(nn_hls)) :: zustar, zustar2 ! 107 REAL(wp), DIMENSION(A2D(nn_hls)) :: zuws, zvws, zsws, zfnet ! 108 REAL(wp), DIMENSION(A2D(nn_hls)) :: zfbuo, zrautbm1, zrautb, zraupl 109 REAL(wp), DIMENSION(A2D(nn_hls)) :: zwpsurf ! 110 REAL(wp), DIMENSION(A2D(nn_hls)) :: zop0 , zsp0 ! 111 REAL(wp), DIMENSION(A2D(nn_hls)) :: zrwp_0, zrwp2_0 ! 112 REAL(wp), DIMENSION(A2D(nn_hls)) :: zapp0 ! 113 REAL(wp), DIMENSION(A2D(nn_hls)) :: zphp, zph, zphpm1, zphm1, zNHydro 114 REAL(wp), DIMENSION(A2D(nn_hls)) :: zhcmo ! 115 ! 116 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zn2 ! N^2 117 REAL(wp), DIMENSION(A2D(nn_hls),2 ) :: zab, zabm1, zabp ! alpha and beta 118 118 119 119 REAL(wp), PARAMETER :: zepsilon = 1.e-30 ! local small value … … 136 136 zcd = 1._wp 137 137 138 !------------------------------------------------------------------ 139 ! Surface boundary condition 140 !------------------------------------------------------------------ 141 ! surface Stress 142 !-------------------- 143 zuws(:,:) = utau(:,:) * r1_rho0 144 zvws(:,:) = vtau(:,:) * r1_rho0 145 zustar2(:,:) = SQRT(zuws(:,:)*zuws(:,:)+zvws(:,:)*zvws(:,:)) 146 zustar(:,:) = SQRT(zustar2(:,:)) 147 148 ! Heat Flux 149 !-------------------- 150 zfnet(:,:) = qns(:,:) + qsr(:,:) 151 zfnet(:,:) = zfnet(:,:) / (rho0 * rcp) 152 153 ! Water Flux 154 !--------------------- 155 zsws(:,:) = emp(:,:) 156 157 !------------------------------------------- 158 ! Initialisation of prognostic variables 159 !------------------------------------------- 160 zrwp (:,:,:) = 0._wp ; zrwp2(:,:,:) = 0._wp ; zedmf(:,:,:) = 0._wp 161 zph (:,:) = 0._wp ; zphm1(:,:) = 0._wp ; zphpm1(:,:) = 0._wp 162 ztsp(:,:,:,:)= 0._wp 163 164 ! Tracers inside plume (ztsp) and environment (ztse) 165 ztsp(:,:,1,jp_tem) = pts(:,:,1,jp_tem,Kmm) * tmask(:,:,1) 166 ztsp(:,:,1,jp_sal) = pts(:,:,1,jp_sal,Kmm) * tmask(:,:,1) 167 ztse(:,:,1,jp_tem) = pts(:,:,1,jp_tem,Kmm) * tmask(:,:,1) 168 ztse(:,:,1,jp_sal) = pts(:,:,1,jp_sal,Kmm) * tmask(:,:,1) 138 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 139 !------------------------------------------------------------------ 140 ! Surface boundary condition 141 !------------------------------------------------------------------ 142 ! surface Stress 143 !-------------------- 144 zuws(ji,jj) = utau(ji,jj) * r1_rho0 145 zvws(ji,jj) = vtau(ji,jj) * r1_rho0 146 zustar2(ji,jj) = SQRT(zuws(ji,jj)*zuws(ji,jj)+zvws(ji,jj)*zvws(ji,jj)) 147 zustar(ji,jj) = SQRT(zustar2(ji,jj)) 148 149 ! Heat Flux 150 !-------------------- 151 zfnet(ji,jj) = qns(ji,jj) + qsr(ji,jj) 152 zfnet(ji,jj) = zfnet(ji,jj) / (rho0 * rcp) 153 154 ! Water Flux 155 !--------------------- 156 zsws(ji,jj) = emp(ji,jj) 157 158 !------------------------------------------- 159 ! Initialisation of prognostic variables 160 !------------------------------------------- 161 zrwp (ji,jj,:) = 0._wp ; zrwp2(ji,jj,:) = 0._wp ; zedmf(ji,jj,:) = 0._wp 162 zph (ji,jj) = 0._wp ; zphm1(ji,jj) = 0._wp ; zphpm1(ji,jj) = 0._wp 163 ztsp(ji,jj,:,:)= 0._wp 164 165 ! Tracers inside plume (ztsp) and environment (ztse) 166 ztsp(ji,jj,1,jp_tem) = pts(ji,jj,1,jp_tem,Kmm) * tmask(ji,jj,1) 167 ztsp(ji,jj,1,jp_sal) = pts(ji,jj,1,jp_sal,Kmm) * tmask(ji,jj,1) 168 ztse(ji,jj,1,jp_tem) = pts(ji,jj,1,jp_tem,Kmm) * tmask(ji,jj,1) 169 ztse(ji,jj,1,jp_sal) = pts(ji,jj,1,jp_sal,Kmm) * tmask(ji,jj,1) 170 END_2D 169 171 170 172 CALL eos( ztse(:,:,1,:) , zrautb(:,:) ) … … 174 176 ! Boundary Condition of Mass Flux (plume velo.; convective area, entrain/detrain) 175 177 !------------------------------------------- 176 zhcmo(:,:) = e3t( :,:,1,Kmm)178 zhcmo(:,:) = e3t(A1Di(nn_hls),A1Dj(nn_hls),1,Kmm) 177 179 zfbuo(:,:) = 0._wp 178 180 WHERE ( ABS(zrautb(:,:)) > 1.e-20 ) zfbuo(:,:) = & 179 & grav * ( 2.e-4_wp *zfnet(:,:) - 7.6E-4_wp*pts(:,:,1,jp_sal,Kmm)*zsws(:,:)/zrautb(:,:)) * zhcmo(:,:) 181 & grav * ( 2.e-4_wp *zfnet(:,:) & 182 & - 7.6E-4_wp*pts(A2D(nn_hls),1,jp_sal,Kmm) & 183 & * zsws(:,:)/zrautb(:,:)) * zhcmo(:,:) 180 184 181 185 zedmf(:,:,1) = -0.065_wp*(ABS(zfbuo(:,:)))**(1._wp/3._wp)*SIGN(1.,zfbuo(:,:)) … … 211 215 CALL eos( ztsp(:,:,jk-1,: ) , zraupl(:,:) ) 212 216 213 zphm1(:,:) = zphm1(:,:) + grav * zrautbm1(:,:) * e3t(:,:,jk-1, Kmm) 214 zphpm1(:,:) = zphpm1(:,:) + grav * zraupl(:,:) * e3t(:,:,jk-1, Kmm) 215 zph(:,:) = zphm1(:,:) + grav * zrautb(:,:) * e3t(:,:,jk , Kmm) 216 zph(:,:) = MAX( zph(:,:), zepsilon) 217 DO_2D( 0, 0, 0, 0 ) 218 zphm1(ji,jj) = zphm1(ji,jj) + grav * zrautbm1(ji,jj) * e3t(ji,jj,jk-1, Kmm) 219 zphpm1(ji,jj) = zphpm1(ji,jj) + grav * zraupl(ji,jj) * e3t(ji,jj,jk-1, Kmm) 220 zph(ji,jj) = zphm1(ji,jj) + grav * zrautb(ji,jj) * e3t(ji,jj,jk , Kmm) 221 zph(ji,jj) = MAX( zph(ji,jj), zepsilon) 222 END_2D 217 223 218 224 WHERE(zrautbm1 .NE. 0.) zfbuo(:,:) = grav * (zraupl(:,:) - zrautbm1(:,:)) / zrautbm1(:,:) 219 225 220 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 221 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 226 DO_2D( 0, 0, 0, 0 ) 222 227 223 228 ! Compute Environment of Plume. Interpolation T/S (before time step) on W-points … … 323 328 324 329 ! Compute Mass Flux on T-point 325 DO jk=1,jpk-1 326 edmfm(:,:,jk) = (zedmf(:,:,jk+1) + zedmf(:,:,jk) )*0.5_wp 327 END DO 328 edmfm(:,:,jpk) = zedmf(:,:,jpk) 330 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 331 edmfm(ji,jj,jk) = (zedmf(ji,jj,jk+1) + zedmf(ji,jj,jk) )*0.5_wp 332 END_3D 333 DO_2D( 0, 0, 0, 0 ) 334 edmfm(ji,jj,jpk) = zedmf(ji,jj,jpk) 335 END_2D 329 336 330 337 ! Save variable (on T point) … … 339 346 ! Computation of a tridiagonal matrix and right hand side terms of the linear system 340 347 !================================================================================= 341 edmfa(:,:,:) = 0._wp 342 edmfb(:,:,:) = 0._wp 343 edmfc(:,:,:) = 0._wp 344 edmftra(:,:,:,:) = 0._wp 348 DO_3D( 0, 0, 0, 0, 1, jpk ) 349 edmfa(ji,jj,jk) = 0._wp 350 edmfb(ji,jj,jk) = 0._wp 351 edmfc(ji,jj,jk) = 0._wp 352 edmftra(ji,jj,jk,:) = 0._wp 353 END_3D 345 354 346 355 !--------------------------------------------------------------- 347 356 ! Diagonal terms 348 357 !--------------------------------------------------------------- 349 DO jk=1,jpk-1 350 edmfa(:,:,jk) = 0._wp 351 edmfb(:,:,jk) = -edmfm(:,:,jk ) / e3w(:,:,jk+1,Kmm) 352 edmfc(:,:,jk) = edmfm(:,:,jk+1) / e3w(:,:,jk+1,Kmm) 353 END DO 354 edmfa(:,:,jpk) = -edmfm(:,:,jpk-1) / e3w(:,:,jpk,Kmm) 355 edmfb(:,:,jpk) = edmfm(:,:,jpk ) / e3w(:,:,jpk,Kmm) 356 edmfc(:,:,jpk) = 0._wp 358 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 359 edmfa(ji,jj,jk) = 0._wp 360 edmfb(ji,jj,jk) = -edmfm(ji,jj,jk ) / e3w(ji,jj,jk+1,Kmm) 361 edmfc(ji,jj,jk) = edmfm(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm) 362 END_3D 363 DO_2D( 0, 0, 0, 0 ) 364 edmfa(ji,jj,jpk) = -edmfm(ji,jj,jpk-1) / e3w(ji,jj,jpk,Kmm) 365 edmfb(ji,jj,jpk) = edmfm(ji,jj,jpk ) / e3w(ji,jj,jpk,Kmm) 366 edmfc(ji,jj,jpk) = 0._wp 367 END_2D 357 368 358 369 !--------------------------------------------------------------- 359 370 ! right hand side term for Temperature 360 371 !--------------------------------------------------------------- 361 DO jk=1,jpk-1 362 edmftra(:,:,jk,1) = - edmfm(:,:,jk ) * ztsp(:,:,jk ,jp_tem) / e3w(:,:,jk+1,Kmm) & 363 & + edmfm(:,:,jk+1) * ztsp(:,:,jk+1,jp_tem) / e3w(:,:,jk+1,Kmm) 364 END DO 365 edmftra(:,:,jpk,1) = - edmfm(:,:,jpk-1) * ztsp(:,:,jpk-1,jp_tem) / e3w(:,:,jpk,Kmm) & 366 & + edmfm(:,:,jpk ) * ztsp(:,:,jpk ,jp_tem) / e3w(:,:,jpk,Kmm) 367 372 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 373 edmftra(ji,jj,jk,1) = - edmfm(ji,jj,jk ) * ztsp(ji,jj,jk ,jp_tem) / e3w(ji,jj,jk+1,Kmm) & 374 & + edmfm(ji,jj,jk+1) * ztsp(ji,jj,jk+1,jp_tem) / e3w(ji,jj,jk+1,Kmm) 375 END_3D 376 DO_2D( 0, 0, 0, 0 ) 377 edmftra(ji,jj,jpk,1) = - edmfm(ji,jj,jpk-1) * ztsp(ji,jj,jpk-1,jp_tem) / e3w(ji,jj,jpk,Kmm) & 378 & + edmfm(ji,jj,jpk ) * ztsp(ji,jj,jpk ,jp_tem) / e3w(ji,jj,jpk,Kmm) 379 END_2D 380 368 381 !--------------------------------------------------------------- 369 382 ! Right hand side term for Salinity 370 383 !--------------------------------------------------------------- 371 DO jk=1,jpk-1 372 edmftra(:,:,jk,2) = - edmfm(:,:,jk ) * ztsp(:,:,jk ,jp_sal) / e3w(:,:,jk+1,Kmm) & 373 & + edmfm(:,:,jk+1) * ztsp(:,:,jk+1,jp_sal) / e3w(:,:,jk+1,Kmm) 374 END DO 375 edmftra(:,:,jpk,2) = - edmfm(:,:,jpk-1) * ztsp(:,:,jpk-1,jp_sal) / e3w(:,:,jpk,Kmm) & 376 & + edmfm(:,:,jpk ) * ztsp(:,:,jpk ,jp_sal) / e3w(:,:,jpk,Kmm) 377 ! 378 ! 379 ! [comm_cleanup] 380 IF (nn_hls.eq.1) CALL lbc_lnk( 'zdfmfc', edmfm,'T',1., edmfa,'T',1., edmfb,'T',1., edmfc,'T',1., edmftra(:,:,:,jp_tem),'T',1., edmftra(:,:,:,jp_sal),'T',1.) 384 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 385 edmftra(ji,jj,jk,2) = - edmfm(ji,jj,jk ) * ztsp(ji,jj,jk ,jp_sal) / e3w(ji,jj,jk+1,Kmm) & 386 & + edmfm(ji,jj,jk+1) * ztsp(ji,jj,jk+1,jp_sal) / e3w(ji,jj,jk+1,Kmm) 387 END_3D 388 DO_2D( 0, 0, 0, 0 ) 389 edmftra(ji,jj,jpk,2) = - edmfm(ji,jj,jpk-1) * ztsp(ji,jj,jpk-1,jp_sal) / e3w(ji,jj,jpk,Kmm) & 390 & + edmfm(ji,jj,jpk ) * ztsp(ji,jj,jpk ,jp_sal) / e3w(ji,jj,jpk,Kmm) 391 END_2D 381 392 ! 382 393 END SUBROUTINE tra_mfc … … 385 396 SUBROUTINE diag_mfc( zdiagi, zdiagd, zdiags, p2dt, Kaa ) 386 397 387 REAL(wp), DIMENSION( jpi,jpj,jpk), INTENT(inout) :: zdiagi, zdiagd, zdiags ! inout: tridaig. terms388 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step389 INTEGER , INTENT(in ) :: Kaa ! ocean time level indices398 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(inout) :: zdiagi, zdiagd, zdiags ! inout: tridaig. terms 399 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 400 INTEGER , INTENT(in ) :: Kaa ! ocean time level indices 390 401 391 402 INTEGER :: ji, jj, jk ! dummy loop arguments -
NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfmxl.F90
r14636 r14665 26 26 PRIVATE 27 27 28 PUBLIC zdf_mxl ! called by zdfphy.F9028 PUBLIC zdf_mxl, zdf_mxl_turb ! called by zdfphy.F90 29 29 30 30 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: nmln !: number of level in the mixed layer (used by LDF, ZDF, TRD, TOP) … … 65 65 !! *** ROUTINE zdfmxl *** 66 66 !! 67 !! ** Purpose : Compute the turbocline depth and the mixed layer depth 68 !! with density criteria. 67 !! ** Purpose : Compute the mixed layer depth with density criteria. 69 68 !! 70 69 !! ** Method : The mixed layer depth is the shallowest W depth with 71 70 !! the density of the corresponding T point (just bellow) bellow a 72 71 !! given value defined locally as rho(10m) + rho_c 73 !! The turbocline depth is the depth at which the vertical74 !! eddy diffusivity coefficient (resulting from the vertical physics75 !! alone, not the isopycnal part, see trazdf.F) fall below a given76 !! value defined locally (avt_c here taken equal to 5 cm/s2 by default)77 72 !! 78 !! ** Action : nmln, hml d, hmlp, hmlpt73 !! ** Action : nmln, hmlp, hmlpt 79 74 !!---------------------------------------------------------------------- 80 75 INTEGER, INTENT(in) :: kt ! ocean time-step index … … 82 77 ! 83 78 INTEGER :: ji, jj, jk ! dummy loop indices 84 INTEGER :: iik n, iiki, ikt! local integer79 INTEGER :: iik, ikt ! local integer 85 80 REAL(wp) :: zN2_c ! local scalar 86 INTEGER, DIMENSION(jpi,jpj) :: imld ! 2D workspace87 81 !!---------------------------------------------------------------------- 88 82 ! 89 IF( kt == nit000 ) THEN 90 IF(lwp) WRITE(numout,*) 91 IF(lwp) WRITE(numout,*) 'zdf_mxl : mixed layer depth' 92 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 93 ! ! allocate zdfmxl arrays 94 IF( zdf_mxl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_mxl : unable to allocate arrays' ) 83 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 84 IF( kt == nit000 ) THEN 85 IF(lwp) WRITE(numout,*) 86 IF(lwp) WRITE(numout,*) 'zdf_mxl : mixed layer depth' 87 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 88 ! ! allocate zdfmxl arrays 89 IF( zdf_mxl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_mxl : unable to allocate arrays' ) 90 ENDIF 95 91 ENDIF 96 92 ! 97 93 ! w-level of the mixing and mixed layers 98 nmln(:,:) = nlb10 ! Initialization to the number of w ocean point 99 hmlp(:,:) = 0._wp ! here hmlp used as a dummy variable, integrating vertically N^2 94 DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls ) 95 nmln(ji,jj) = nlb10 ! Initialization to the number of w ocean point 96 hmlp(ji,jj) = 0._wp ! here hmlp used as a dummy variable, integrating vertically N^2 97 END_2D 100 98 zN2_c = grav * rho_c * r1_rho0 ! convert density criteria into N^2 criteria 101 99 ! [comm_cleanup] ! DO_3D( 1, 1, 1, 1, nlb10, jpkm1 ) ! Mixed layer level: w-level 102 DO_3D ( nn_hls, nn_hls, nn_hls, nn_hls, nlb10, jpkm1 ) ! Mixed layer level: w-level100 DO_3D_OVR( nn_hls, nn_hls, nn_hls, nn_hls, nlb10, jpkm1 ) ! Mixed layer level: w-level 103 101 ikt = mbkt(ji,jj) 104 102 hmlp(ji,jj) = & … … 106 104 IF( hmlp(ji,jj) < zN2_c ) nmln(ji,jj) = MIN( jk , ikt ) + 1 ! Mixed layer level 107 105 END_3D 108 ! 109 ! w-level of the turbocline and mixing layer (iom_use) 110 imld(:,:) = mbkt(:,:) + 1 ! Initialization to the number of w ocean point 111 ! [comm_cleanup] ! DO_3DS( 1, 1, 1, 1, jpkm1, nlb10, -1 ) ! from the bottom to nlb10 112 DO_3DS( nn_hls, nn_hls, nn_hls, nn_hls, jpkm1, nlb10, -1 ) ! from the bottom to nlb10 113 IF( avt (ji,jj,jk) < avt_c * wmask(ji,jj,jk) ) imld(ji,jj) = jk ! Turbocline 114 END_3D 115 ! depth of the mixing and mixed layers 106 ! depth of the mixed layer 116 107 ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 ) 117 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 118 iiki = imld(ji,jj) 119 iikn = nmln(ji,jj) 120 hmld (ji,jj) = gdepw(ji,jj,iiki ,Kmm) * ssmask(ji,jj) ! Turbocline depth 121 hmlp (ji,jj) = gdepw(ji,jj,iikn ,Kmm) * ssmask(ji,jj) ! Mixed layer depth 122 hmlpt(ji,jj) = gdept(ji,jj,iikn-1,Kmm) * ssmask(ji,jj) ! depth of the last T-point inside the mixed layer 108 DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls ) 109 iik = nmln(ji,jj) 110 hmlp (ji,jj) = gdepw(ji,jj,iik ,Kmm) * ssmask(ji,jj) ! Mixed layer depth 111 hmlpt(ji,jj) = gdept(ji,jj,iik-1,Kmm) * ssmask(ji,jj) ! depth of the last T-point inside the mixed layer 123 112 END_2D 124 113 ! 125 IF( .NOT.l_offline ) THEN 126 IF( iom_use("mldr10_1") ) THEN 127 IF( ln_isfcav ) THEN ; CALL iom_put( "mldr10_1", hmlp - risfdep) ! mixed layer thickness 128 ELSE ; CALL iom_put( "mldr10_1", hmlp ) ! mixed layer depth 129 END IF 114 IF( .NOT.l_offline .AND. iom_use("mldr10_1") ) THEN 115 IF( ln_isfcav ) THEN ; CALL iom_put( "mldr10_1", hmlp - risfdep) ! mixed layer thickness 116 ELSE ; CALL iom_put( "mldr10_1", hmlp ) ! mixed layer depth 130 117 END IF 131 IF( iom_use("mldkz5") ) THEN132 IF( ln_isfcav ) THEN ; CALL iom_put( "mldkz5" , hmld - risfdep ) ! turbocline thickness133 ELSE ; CALL iom_put( "mldkz5" , hmld ) ! turbocline depth134 END IF135 ENDIF136 118 ENDIF 137 119 ! … … 140 122 END SUBROUTINE zdf_mxl 141 123 124 125 SUBROUTINE zdf_mxl_turb( kt, Kmm ) 126 !!---------------------------------------------------------------------- 127 !! *** ROUTINE zdf_mxl_turb *** 128 !! 129 !! ** Purpose : Compute the turbocline depth. 130 !! 131 !! ** Method : The turbocline depth is the depth at which the vertical 132 !! eddy diffusivity coefficient (resulting from the vertical physics 133 !! alone, not the isopycnal part, see trazdf.F) fall below a given 134 !! value defined locally (avt_c here taken equal to 5 cm/s2 by default) 135 !! 136 !! ** Action : hmld 137 !!---------------------------------------------------------------------- 138 INTEGER, INTENT(in) :: kt ! ocean time-step index 139 INTEGER, INTENT(in) :: Kmm ! ocean time level index 140 ! 141 INTEGER :: ji, jj, jk ! dummy loop indices 142 INTEGER :: iik ! local integer 143 INTEGER, DIMENSION(A2D(nn_hls)) :: imld ! 2D workspace 144 !!---------------------------------------------------------------------- 145 ! 146 ! w-level of the turbocline and mixing layer (iom_use) 147 imld(:,:) = mbkt(A2D(nn_hls)) + 1 ! Initialization to the number of w ocean point 148 DO_3DS( 0, 0, 0, 0, jpkm1, nlb10, -1 ) ! from the bottom to nlb10 149 IF( avt (ji,jj,jk) < avt_c * wmask(ji,jj,jk) ) imld(ji,jj) = jk ! Turbocline 150 END_3D 151 ! depth of the mixing layer 152 DO_2D( 0, 0, 0, 0 ) 153 iik = imld(ji,jj) 154 hmld (ji,jj) = gdepw(ji,jj,iik ,Kmm) * ssmask(ji,jj) ! Turbocline depth 155 END_2D 156 ! 157 IF( .NOT.l_offline .AND. iom_use("mldkz5") ) THEN 158 IF( ln_isfcav ) THEN ; CALL iom_put( "mldkz5" , hmld - risfdep ) ! turbocline thickness 159 ELSE ; CALL iom_put( "mldkz5" , hmld ) ! turbocline depth 160 END IF 161 ENDIF 162 ! 163 END SUBROUTINE zdf_mxl_turb 142 164 !!====================================================================== 143 165 END MODULE zdfmxl -
NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfosm.F90
r14636 r14665 1214 1214 ! Lateral boundary conditions on final outputs for hbl, on T-grid (sign unchanged) 1215 1215 ! [comm_cleanup] ! no need lbc_lnk for output 1216 ! CALL lbc_lnk( 'zdfosm', hbl, 'T', 1., dh, 'T', 1., hmle, 'T', 1. ) 1216 ! NOTE: [tiling] still needed (at least for hmle in tra_mle_trp) 1217 CALL lbc_lnk( 'zdfosm', hbl, 'T', 1., dh, 'T', 1., hmle, 'T', 1. ) 1217 1218 ! Lateral boundary conditions on final outputs for gham[ts], on W-grid (sign unchanged) 1218 1219 ! Lateral boundary conditions on final outputs for gham[uv], on [UV]-grid (sign changed) … … 2874 2875 ! 2875 2876 IF( kt == nit000 ) THEN 2876 IF( ntile == 0.OR. ntile == 1 ) THEN ! Do only on the first tile2877 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 2877 2878 IF(lwp) WRITE(numout,*) 2878 2879 IF(lwp) WRITE(numout,*) 'tra_osm : OSM non-local tracer fluxes' … … 2886 2887 ENDIF 2887 2888 2888 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 2889 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 2889 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 2890 2890 pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) & 2891 2891 & - ( ghamt(ji,jj,jk ) & … … 2948 2948 ! 2949 2949 IF( kt == nit000 ) THEN 2950 IF(lwp) WRITE(numout,*) 2951 IF(lwp) WRITE(numout,*) 'dyn_osm : OSM non-local velocity' 2952 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 2950 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 2951 IF(lwp) WRITE(numout,*) 2952 IF(lwp) WRITE(numout,*) 'dyn_osm : OSM non-local velocity' 2953 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 2954 ENDIF 2953 2955 ENDIF 2954 2956 !code saving tracer trends removed, replace with trdmxl_oce 2955 2957 2956 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! add non-local u and v fluxes 2957 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) ! add non-local u and v fluxes 2958 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! add non-local u and v fluxes 2958 2959 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) & 2959 2960 & - ( ghamu(ji,jj,jk ) & -
NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfphy.F90
r14636 r14665 12 12 !!---------------------------------------------------------------------- 13 13 USE oce ! ocean dynamics and tracers variables 14 ! TEMP: [tiling] This change not necessary after finalisation of zdf_osm 15 USE domtile 14 16 USE zdf_oce ! vertical physics: shared variables 15 17 USE zdfdrg ! vertical physics: top/bottom drag coef. … … 56 58 LOGICAL, PUBLIC :: l_zdfsh2 ! shear production term flag (=F for CST, =T otherwise (i.e. TKE, GLS, RIC)) 57 59 60 REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: avm_k_n !: "Now" avm_k used for calculation of zsh2 with tiling 61 62 !! * Substitutions 63 # include "do_loop_substitute.h90" 58 64 !!---------------------------------------------------------------------- 59 65 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 210 216 ENDIF 211 217 ! ! shear production term flag 212 IF( ln_zdfcst ) THEN ; l_zdfsh2 = .FALSE. 213 ELSE ; l_zdfsh2 = .TRUE. 214 ENDIF 218 IF( ln_zdfcst .OR. ln_zdfosm ) THEN ; l_zdfsh2 = .FALSE. 219 ELSE ; l_zdfsh2 = .TRUE. 220 ENDIF 221 IF( ln_tile .AND. l_zdfsh2 ) ALLOCATE( avm_k_n(jpi,jpj,jpk) ) 215 222 ! !== Mass Flux Convectiive algorithm ==! 216 223 IF( ln_zdfmfc ) CALL zdf_mfc_init ! Convection computed with eddy diffusivity mass flux … … 246 253 ! 247 254 INTEGER :: ji, jj, jk ! dummy loop indice 248 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zsh2 ! shear production 255 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zsh2 ! shear production 256 ! TEMP: [tiling] This change not necessary after finalisation of zdf_osm 257 LOGICAL :: lskip 249 258 !! --------------------------------------------------------------------- 250 259 ! 251 260 IF( ln_timing ) CALL timing_start('zdf_phy') 261 262 ! TEMP: [tiling] These changes not necessary after finalisation of zdf_osm 263 lskip = .FALSE. 264 265 IF( ln_tile .AND. nzdf_phy == np_OSM ) THEN 266 IF( ntile == 1 ) THEN 267 CALL dom_tile_stop( ldhold=.TRUE., cstr='zdfphy' ) 268 ELSE 269 lskip = .TRUE. 270 ENDIF 271 ENDIF 252 272 ! 253 273 IF( l_zdfdrg ) THEN !== update top/bottom drag ==! (non-linear cases) … … 264 284 ENDIF 265 285 ! 286 ! TODO: [tiling] NOT TESTED- requires SI3 and isfcav 266 287 #if defined key_si3 267 288 IF ( ln_drgice_imp) THEN 268 289 IF ( ln_isfcav ) THEN 269 rCdU_top(:,:) = rCdU_top(:,:) + ssmask(:,:) * tmask(:,:,1) * rCdU_ice(:,:) 290 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 291 rCdU_top(ji,jj) = rCdU_top(ji,jj) + ssmask(ji,jj) * tmask(ji,jj,1) * rCdU_ice(ji,jj) 292 END_2D 270 293 ELSE 271 rCdU_top(:,:) = rCdU_ice(:,:) 294 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 295 rCdU_top(ji,jj) = rCdU_ice(ji,jj) 296 END_2D 272 297 ENDIF 273 298 ENDIF 274 299 #endif 275 300 ! 276 ! !== Kz from chosen turbulent closure ==! (avm_k, avt_k) 277 ! 278 IF( l_zdfsh2 ) & !* shear production at w-points (energy conserving form) 279 CALL zdf_sh2( Kbb, Kmm, avm_k, & ! <<== in 280 & zsh2 ) ! ==>> out : shear production 281 ! 282 SELECT CASE ( nzdf_phy ) !* Vertical eddy viscosity and diffusivity coefficients at w-points 283 CASE( np_RIC ) ; CALL zdf_ric( kt, Kmm, zsh2, avm_k, avt_k ) ! Richardson number dependent Kz 284 CASE( np_TKE ) ; CALL zdf_tke( kt, Kbb, Kmm, zsh2, avm_k, avt_k ) ! TKE closure scheme for Kz 285 CASE( np_GLS ) ; CALL zdf_gls( kt, Kbb, Kmm, zsh2, avm_k, avt_k ) ! GLS closure scheme for Kz 286 ! [comm_cleanup] ! modified but not tested - no ref config uses this scheme 287 CASE( np_OSM ) ; CALL zdf_osm( kt, Kbb, Kmm, Krhs, avm_k, avt_k ) ! OSMOSIS closure scheme for Kz 288 ! CASE( np_CST ) ! Constant Kz (reset avt, avm to the background value) 289 ! ! avt_k and avm_k set one for all at initialisation phase 301 ! NOTE: [tiling] zdf_mxl has been split into zdf_mxl (nmln, hmlp, hmlpt) and zdf_mxl_turb (hmld), which depend on rn2b and avt respectively. zdf_mxl has been moved here from below (which is now zdf_mxl_turb) because zdf_tke with nn_etau = 2 depends on nmln. The previous solution was to call zdf_mxl in zdf_tke_init, but rn2b is not yet calculated here (0 everywhere), so this specific option was not restartable and tiling changed the results. Furthermore, the MLD diagnostics sent from zdf_mxl were incorrect for the first timestep. This bug fix has therefore changed the results when using zdf_tke with nn_etau = 2. 302 CALL zdf_mxl( kt, Kmm ) !* mixed layer depth, and level 303 304 ! TEMP: [tiling] These changes not necessary after finalisation of zdf_osm 305 IF( .NOT. lskip ) THEN 306 ! !== Kz from chosen turbulent closure ==! (avm_k, avt_k) 307 ! 308 ! NOTE: [tiling] the closure schemes (zdf_tke etc) will update avm_k. With tiling, the calculation of zsh2 on adjacent tiles then uses both updated (next timestep) and non-updated (current timestep) values of avm_k. To preserve results, we save a read-only copy of the "now" avm_k to use in the calculation of zsh2. 309 IF( l_zdfsh2 ) THEN !* shear production at w-points (energy conserving form) 310 IF( ln_tile ) THEN 311 IF( ntile == 1 ) avm_k_n(:,:,:) = avm_k(:,:,:) ! Preserve "now" avm_k for calculation of zsh2 312 CALL zdf_sh2( Kbb, Kmm, avm_k_n, & ! <<== in 313 & zsh2 ) ! ==>> out : shear production 314 ELSE 315 CALL zdf_sh2( Kbb, Kmm, avm_k, & ! <<== in 316 & zsh2 ) ! ==>> out : shear production 317 ENDIF 318 ENDIF 319 ! 320 SELECT CASE ( nzdf_phy ) !* Vertical eddy viscosity and diffusivity coefficients at w-points 321 CASE( np_RIC ) ; CALL zdf_ric( kt, Kmm, zsh2, avm_k, avt_k ) ! Richardson number dependent Kz 322 CASE( np_TKE ) ; CALL zdf_tke( kt, Kbb, Kmm, zsh2, avm_k, avt_k ) ! TKE closure scheme for Kz 323 CASE( np_GLS ) ; CALL zdf_gls( kt, Kbb, Kmm, zsh2, avm_k, avt_k ) ! GLS closure scheme for Kz 324 CASE( np_OSM ) ; CALL zdf_osm( kt, Kbb, Kmm, Krhs, avm_k, avt_k ) ! OSMOSIS closure scheme for Kz 325 ! CASE( np_CST ) ! Constant Kz (reset avt, avm to the background value) 326 ! ! avt_k and avm_k set one for all at initialisation phase 290 327 !!gm avt(2:jpim1,2:jpjm1,1:jpkm1) = rn_avt0 * wmask(2:jpim1,2:jpjm1,1:jpkm1) 291 328 !!gm avm(2:jpim1,2:jpjm1,1:jpkm1) = rn_avm0 * wmask(2:jpim1,2:jpjm1,1:jpkm1) 292 END SELECT 329 END SELECT 330 331 IF( ln_tile .AND. .NOT. l_istiled ) CALL dom_tile_start( ldhold=.TRUE., cstr='zdfphy' ) 332 ENDIF 293 333 ! 294 334 ! !== ocean Kz ==! (avt, avs, avm) 295 335 ! 296 336 ! !* start from turbulent closure values 297 avt(:,:,2:jpkm1) = avt_k(:,:,2:jpkm1) 298 avm(:,:,2:jpkm1) = avm_k(:,:,2:jpkm1) 337 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 338 avt(ji,jj,jk) = avt_k(ji,jj,jk) 339 avm(ji,jj,jk) = avm_k(ji,jj,jk) 340 END_3D 299 341 ! 300 342 IF( ln_rnf_mouth ) THEN !* increase diffusivity at rivers mouths 301 DO jk = 2, nkrnf302 avt( :,:,jk) = avt(:,:,jk) + 2._wp * rn_avt_rnf * rnfmsk(:,:) * wmask(:,:,jk)303 END DO343 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, nkrnf ) 344 avt(ji,jj,jk) = avt(ji,jj,jk) + 2._wp * rn_avt_rnf * rnfmsk(ji,jj) * wmask(ji,jj,jk) 345 END_3D 304 346 ENDIF 305 347 ! … … 310 352 CALL zdf_ddm( kt, Kmm, avm, avt, avs ) 311 353 ELSE ! same mixing on all tracers 312 avs(2:jpim1,2:jpjm1,1:jpkm1) = avt(2:jpim1,2:jpjm1,1:jpkm1) 354 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 355 avs(ji,jj,jk) = avt(ji,jj,jk) 356 END_3D 313 357 ENDIF 314 358 ! 315 359 ! !* wave-induced mixing 360 ! TODO: [tiling] NOT TESTED- requires ln_wave and ln_sdw 316 361 IF( ln_zdfswm ) CALL zdf_swm( kt, Kmm, avm, avt, avs ) ! surface wave (Qiao et al. 2004) 317 362 IF( ln_zdfiwm ) CALL zdf_iwm( kt, Kmm, avm, avt, avs ) ! internal wave (de Lavergne et al 2017) … … 319 364 #if defined key_agrif 320 365 ! interpolation parent grid => child grid for avm_k ( ex : at west border: update column 1 and 2) 321 IF( l_zdfsh2 ) CALL Agrif_avm 366 ! TODO: [tiling] NOT TESTED- requires AGRIF 367 IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN ! Do only on the last tile 368 IF( l_zdfsh2 ) CALL Agrif_avm 369 ENDIF 322 370 #endif 323 371 324 372 ! !* Lateral boundary conditions (sign unchanged) 325 373 ! [comm_cleanup] ! lbc_lnk shifted in stp 326 IF(nn_hls.eq.1) THEN 327 IF( l_zdfsh2 ) THEN 328 CALL lbc_lnk( 'zdfphy', avm_k, 'W', 1.0_wp , avt_k, 'W', 1.0_wp, & 329 & avm , 'W', 1.0_wp , avt , 'W', 1.0_wp , avs , 'W', 1.0_wp ) 330 ELSE 331 CALL lbc_lnk( 'zdfphy', avm , 'W', 1.0_wp , avt , 'W', 1.0_wp , avs , 'W', 1.0_wp ) 332 ENDIF 333 ! 334 IF( l_zdfdrg ) THEN ! drag have been updated (non-linear cases) 335 IF( ln_isfcav ) THEN ; CALL lbc_lnk( 'zdfphy', rCdU_top, 'T', 1.0_wp , rCdU_bot, 'T', 1.0_wp ) ! top & bot drag 336 ELSE ; CALL lbc_lnk( 'zdfphy', rCdU_bot, 'T', 1.0_wp ) ! bottom drag only 374 IF(nn_hls.eq.1) THEN 375 ! NOTE: [tiling] this solution for lbc_lnk will not work if outer halo values are accessed before the end of the tiling loop 376 IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN ! Do only on the last tile 377 IF( l_zdfsh2 ) THEN 378 CALL lbc_lnk( 'zdfphy', avm_k, 'W', 1.0_wp , avt_k, 'W', 1.0_wp, & 379 & avm , 'W', 1.0_wp , avt , 'W', 1.0_wp , avs , 'W', 1.0_wp ) 380 ELSE 381 CALL lbc_lnk( 'zdfphy', avm , 'W', 1.0_wp , avt , 'W', 1.0_wp , avs , 'W', 1.0_wp ) 337 382 ENDIF 338 ENDIF 339 ENDIF 340 ! 341 CALL zdf_mxl( kt, Kmm ) !* mixed layer depth, and level 342 ! 343 IF( lrst_oce ) THEN !* write TKE, GLS or RIC fields in the restart file 344 IF( ln_zdftke ) CALL tke_rst( kt, 'WRITE' ) 345 IF( ln_zdfgls ) CALL gls_rst( kt, 'WRITE' ) 346 IF( ln_zdfric ) CALL ric_rst( kt, 'WRITE' ) 347 ! NB. OSMOSIS restart (osm_rst) will be called in step.F90 after ww has been updated 383 ! 384 IF( l_zdfdrg ) THEN ! drag have been updated (non-linear cases) 385 IF( ln_isfcav ) THEN ; CALL lbc_lnk( 'zdfphy', rCdU_top, 'T', 1.0_wp , rCdU_bot, 'T', 1.0_wp ) ! top & bot drag 386 ELSE ; CALL lbc_lnk( 'zdfphy', rCdU_bot, 'T', 1.0_wp ) ! bottom drag only 387 ENDIF 388 ENDIF 389 ENDIF 390 ENDIF 391 ! 392 CALL zdf_mxl_turb( kt, Kmm ) !* turbocline depth 393 ! 394 IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN ! Do only on the last tile 395 IF( lrst_oce ) THEN !* write TKE, GLS or RIC fields in the restart file 396 IF( ln_zdftke ) CALL tke_rst( kt, 'WRITE' ) 397 IF( ln_zdfgls ) CALL gls_rst( kt, 'WRITE' ) 398 IF( ln_zdfric ) CALL ric_rst( kt, 'WRITE' ) 399 ! NB. OSMOSIS restart (osm_rst) will be called in step.F90 after ww has been updated 400 ENDIF 348 401 ENDIF 349 402 ! -
NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfric.F90
r14636 r14665 51 51 !! * Substitutions 52 52 # include "do_loop_substitute.h90" 53 # include "domzgr_substitute.h90" 53 54 !!---------------------------------------------------------------------- 54 55 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 145 146 !! PFJ Lermusiaux 2001. 146 147 !!---------------------------------------------------------------------- 147 INTEGER , INTENT(in ) :: kt ! ocean time-step148 INTEGER , INTENT(in ) :: Kmm ! ocean time level index149 REAL(wp), DIMENSION( :,:,:), INTENT(in ) :: p_sh2 ! shear production term150 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: p_avm, p_avt ! momentum and tracer Kz (w-points)148 INTEGER , INTENT(in ) :: kt ! ocean time-step 149 INTEGER , INTENT(in ) :: Kmm ! ocean time level index 150 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in ) :: p_sh2 ! shear production term 151 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: p_avm, p_avt ! momentum and tracer Kz (w-points) 151 152 !! 152 153 INTEGER :: ji, jj, jk ! dummy loop indices 153 154 REAL(wp) :: zcfRi, zav, zustar, zhek ! local scalars 154 REAL(wp), DIMENSION( jpi,jpj) :: zh_ekm ! 2D workspace155 REAL(wp), DIMENSION(A2D(nn_hls)) :: zh_ekm ! 2D workspace 155 156 !!---------------------------------------------------------------------- 156 157 ! 157 158 ! !== avm and avt = F(Richardson number) ==! 158 159 ! [comm_cleanup] ! DO_3D( 1, 0, 1, 0, 2, jpkm1 ) ! coefficient = F(richardson number) (avm-weighted Ri) 159 DO_3D ( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 2, jpkm1 ) ! coefficient = F(richardson number) (avm-weighted Ri)160 DO_3D_OVR( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 2, jpkm1 ) ! coefficient = F(richardson number) (avm-weighted Ri) 160 161 zcfRi = 1._wp / ( 1._wp + rn_alp * MAX( 0._wp , avm(ji,jj,jk) * rn2(ji,jj,jk) / ( p_sh2(ji,jj,jk) + 1.e-20 ) ) ) 161 162 zav = rn_avmri * zcfRi**nn_ric … … 177 178 END_2D 178 179 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 ) !* minimum mixing coeff. within the Ekman layer 179 DO_3D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) !* minimum mixing coeff. within the Ekman layer180 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) !* minimum mixing coeff. within the Ekman layer 180 181 IF( gdept(ji,jj,jk,Kmm) < zh_ekm(ji,jj) ) THEN 181 182 p_avm(ji,jj,jk) = MAX( p_avm(ji,jj,jk), rn_wvmix ) * wmask(ji,jj,jk) -
NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfsh2.F90
r14636 r14665 55 55 !! References : Bruchard, OM 2002 56 56 !! --------------------------------------------------------------------- 57 INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices58 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: p_avm ! vertical eddy viscosity (w-points)59 REAL(wp), DIMENSION( :,:,:) , INTENT( out) :: p_sh2 ! shear production of TKE (w-points)57 INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices 58 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: p_avm ! vertical eddy viscosity (w-points) 59 REAL(wp), DIMENSION(A2D(nn_hls),jpk) , INTENT( out) :: p_sh2 ! shear production of TKE (w-points) 60 60 ! 61 61 INTEGER :: ji, jj, jk ! dummy loop arguments 62 REAL(wp), DIMENSION( jpi,jpj) :: zsh2u, zsh2v ! 2D workspace62 REAL(wp), DIMENSION(A2D(nn_hls)) :: zsh2u, zsh2v ! 2D workspace 63 63 !!-------------------------------------------------------------------- 64 64 ! 65 65 DO jk = 2, jpkm1 !* Shear production at uw- and vw-points (energy conserving form) 66 ! TODO: [tiling] NOT TESTED- requires key_oasis3 66 67 IF ( cpl_sdrftx .AND. ln_stshear ) THEN ! Surface Stokes Drift available ===>>> shear + stokes drift contibution 67 68 ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 ) -
NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfswm.F90
r14636 r14665 64 64 zcoef = 1._wp * 0.353553_wp 65 65 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 66 DO_3D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )66 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 67 67 zqb = zcoef * hsw(ji,jj) * tsd2d(ji,jj) * EXP( -3. * wnum(ji,jj) * gdepw(ji,jj,jk,Kmm) ) * wmask(ji,jj,jk) 68 68 ! -
NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdftke.F90
r14636 r14665 168 168 !! Bruchard OM 2002 169 169 !!---------------------------------------------------------------------- 170 INTEGER , INTENT(in ) :: kt ! ocean time step171 INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices172 REAL(wp), DIMENSION( :,:,:), INTENT(in ) :: p_sh2 ! shear production term173 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: p_avm, p_avt ! momentum and tracer Kz (w-points)170 INTEGER , INTENT(in ) :: kt ! ocean time step 171 INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices 172 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in ) :: p_sh2 ! shear production term 173 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: p_avm, p_avt ! momentum and tracer Kz (w-points) 174 174 !!---------------------------------------------------------------------- 175 175 ! … … 201 201 USE zdf_oce , ONLY : en ! ocean vertical physics 202 202 !! 203 INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices204 REAL(wp), DIMENSION( :,:,:) , INTENT(in ) :: p_sh2 ! shear production term205 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: p_avm, p_avt ! vertical eddy viscosity & diffusivity (w-points)203 INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices 204 REAL(wp), DIMENSION(A2D(nn_hls),jpk) , INTENT(in ) :: p_sh2 ! shear production term 205 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: p_avm, p_avt ! vertical eddy viscosity & diffusivity (w-points) 206 206 ! 207 207 INTEGER :: ji, jj, jk ! dummy loop arguments … … 216 216 REAL(wp) :: zzd_up, zzd_lw ! - - 217 217 REAL(wp) :: ztaui, ztauj, z1_norm 218 INTEGER , DIMENSION( jpi,jpj) :: imlc219 REAL(wp), DIMENSION( jpi,jpj) :: zice_fra, zhlc, zus3, zWlc2220 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zpelc, zdiag, zd_up, zd_lw218 INTEGER , DIMENSION(A2D(nn_hls)) :: imlc 219 REAL(wp), DIMENSION(A2D(nn_hls)) :: zice_fra, zhlc, zus3, zWlc2 220 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zpelc, zdiag, zd_up, zd_lw 221 221 !!-------------------------------------------------------------------- 222 222 ! … … 232 232 SELECT CASE ( nn_eice ) 233 233 CASE( 0 ) ; zice_fra(:,:) = 0._wp 234 CASE( 1 ) ; zice_fra(:,:) = TANH( fr_i( :,:) * 10._wp )235 CASE( 2 ) ; zice_fra(:,:) = fr_i( :,:)236 CASE( 3 ) ; zice_fra(:,:) = MIN( 4._wp * fr_i( :,:) , 1._wp )234 CASE( 1 ) ; zice_fra(:,:) = TANH( fr_i(A2D(nn_hls)) * 10._wp ) 235 CASE( 2 ) ; zice_fra(:,:) = fr_i(A2D(nn_hls)) 236 CASE( 3 ) ; zice_fra(:,:) = MIN( 4._wp * fr_i(A2D(nn_hls)) , 1._wp ) 237 237 END SELECT 238 238 ! … … 242 242 ! 243 243 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 244 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )244 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 245 245 en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) 246 246 zdiag(ji,jj,1) = 1._wp/en(ji,jj,1) … … 260 260 ! 261 261 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) ! bottom friction 262 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! bottom friction262 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! bottom friction 263 263 zmsku = ( 2. - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 264 264 zmskv = ( 2. - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) … … 270 270 IF( ln_isfcav ) THEN 271 271 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) ! top friction 272 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! top friction272 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! top friction 273 273 zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) ) 274 274 zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) ) … … 290 290 ! !* Langmuir velocity scale 291 291 ! 292 ! TODO: [tiling] NOT TESTED- requires key_oasis3 292 293 IF ( cpl_sdrftx ) THEN ! Surface Stokes Drift available 293 294 ! ! Craik-Leibovich velocity scale Wlc = ( u* u_s )^1/2 with u* = (taum/rho0)^1/2 … … 312 313 zWlc2(ji,jj) = 0.5_wp * z1_norm * ( MAX( ut0sd(ji,jj)*ztaui + vt0sd(ji,jj)*ztauj, 0._wp ) )**2 313 314 END_2D 314 ! [comm_cleanup]315 IF (nn_hls.eq.1) CALL lbc_lnk ( 'zdftke', zWlc2, 'T', 1. )316 !317 315 ELSE ! Surface Stokes drift deduced from surface stress 318 316 ! ! Wlc = u_s with u_s = 0.016*U_10m, the surface stokes drift (Axell 2002, Eq.44) … … 322 320 zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag ) ! to convert stress in 10m wind using a constant drag 323 321 ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 ) 324 DO_2D( nn_hls , nn_hls, nn_hls, nn_hls)322 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 325 323 zWlc2(ji,jj) = zcof * taum(ji,jj) 326 324 END_2D … … 330 328 ! !* Depth of the LC circulation (Axell 2002, Eq.47) 331 329 ! !- LHS of Eq.47 332 zpelc(:,:,1) = MAX( rn2b(:,:,1), 0._wp ) * gdepw(:,:,1,Kmm) * e3w(:,:,1,Kmm) 333 DO jk = 2, jpk 334 zpelc(:,:,jk) = zpelc(:,:,jk-1) + & 335 & MAX( rn2b(:,:,jk), 0._wp ) * gdepw(:,:,jk,Kmm) * e3w(:,:,jk,Kmm) 336 END DO 330 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 331 zpelc(ji,jj,1) = MAX( rn2b(ji,jj,1), 0._wp ) * gdepw(ji,jj,1,Kmm) * e3w(ji,jj,1,Kmm) 332 END_2D 333 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpk ) 334 zpelc(ji,jj,jk) = zpelc(ji,jj,jk-1) + & 335 & MAX( rn2b(ji,jj,jk), 0._wp ) * gdepw(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) 336 END_3D 337 337 ! 338 338 ! !- compare LHS to RHS of Eq.47 339 imlc(:,:) = mbkt( :,:) + 1 ! Initialization to the number of w ocean point (=2 over land)339 imlc(:,:) = mbkt(A2D(nn_hls)) + 1 ! Initialization to the number of w ocean point (=2 over land) 340 340 ! [comm_cleanup] ! DO_3DS( 1, 1, 1, 1, jpkm1, 2, -1 ) 341 DO_3DS( nn_hls , nn_hls, nn_hls, nn_hls, jpkm1, 2, -1 )341 DO_3DS( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpkm1, 2, -1 ) 342 342 IF( zpelc(ji,jj,jk) > zWlc2(ji,jj) ) imlc(ji,jj) = jk 343 343 END_3D 344 344 ! ! finite LC depth 345 345 ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 ) 346 DO_2D( nn_hls , nn_hls, nn_hls, nn_hls)346 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 347 347 zhlc(ji,jj) = gdepw(ji,jj,imlc(ji,jj),Kmm) 348 348 END_2D … … 355 355 END_2D 356 356 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 ) !* TKE Langmuir circulation source term added to en 357 DO_3D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) !* TKE Langmuir circulation source term added to en357 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) !* TKE Langmuir circulation source term added to en 358 358 IF ( zus3(ji,jj) /= 0._wp ) THEN 359 359 IF ( gdepw(ji,jj,jk,Kmm) - zhlc(ji,jj) < 0 .AND. wmask(ji,jj,jk) /= 0. ) THEN … … 377 377 IF( nn_pdl == 1 ) THEN !* Prandtl number = F( Ri ) 378 378 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 379 DO_3D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )379 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 380 380 ! ! local Richardson number 381 381 IF (rn2b(ji,jj,jk) <= 0.0_wp) then … … 390 390 ! 391 391 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 ) !* Matrix and right hand side in en 392 DO_3D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) !* Matrix and right hand side in en392 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) !* Matrix and right hand side in en 393 393 zcof = zfact1 * tmask(ji,jj,jk) 394 394 ! ! A minimum of 2.e-5 m2/s is imposed on TKE vertical … … 415 415 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 416 416 ! 417 ! TODO: [tiling] NOT TESTED- requires key_oasis3 and waves 417 418 IF ( cpl_phioc .and. ln_phioc ) THEN 418 419 SELECT CASE (nn_bc_surf) ! Boundary Condition using surface TKE flux from waves … … 420 421 CASE ( 0 ) ! Dirichlet BC 421 422 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) ! en(1) = rn_ebb taum / rho0 (min value rn_emin0) 422 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! en(1) = rn_ebb taum / rho0 (min value rn_emin0)423 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! en(1) = rn_ebb taum / rho0 (min value rn_emin0) 423 424 IF ( phioc(ji,jj) < 0 ) phioc(ji,jj) = 0._wp 424 425 en(ji,jj,1) = MAX( rn_emin0, .5 * ( 15.8 * phioc(ji,jj) / rho0 )**(2./3.) ) * tmask(ji,jj,1) … … 428 429 CASE ( 1 ) ! Neumann BC 429 430 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 430 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )431 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 431 432 IF ( phioc(ji,jj) < 0 ) phioc(ji,jj) = 0._wp 432 433 en(ji,jj,2) = en(ji,jj,2) + ( rn_Dt * phioc(ji,jj) / rho0 ) /e3w(ji,jj,2,Kmm) … … 455 456 END_3D 456 457 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 457 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk458 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 458 459 en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1) 459 460 END_2D 460 461 ! [comm_cleanup] ! DO_3DS( 0, 0, 0, 0, jpk-2, 2, -1 ) 461 DO_3DS ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpk-2, 2, -1 )462 DO_3DS_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpk-2, 2, -1 ) 462 463 en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 463 464 END_3D 464 465 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! set the minimum value of tke 465 DO_3D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! set the minimum value of tke466 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! set the minimum value of tke 466 467 en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk) 467 468 END_3D … … 477 478 IF( nn_etau == 1 ) THEN !* penetration below the mixed layer (rn_efr fraction) 478 479 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 479 DO_3D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )480 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 480 481 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) ) & 481 482 & * MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) … … 483 484 ELSEIF( nn_etau == 2 ) THEN !* act only at the base of the mixed layer (jk=nmln) (rn_efr fraction) 484 485 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 485 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )486 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 486 487 jk = nmln(ji,jj) 487 488 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) ) & 488 489 & * MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 489 490 END_2D 491 ! TODO: [tiling] NOT TESTED- requires ln_cpl 490 492 ELSEIF( nn_etau == 3 ) THEN !* penetration belox the mixed layer (HF variability) 491 493 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 492 DO_3D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )494 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 493 495 ztx2 = utau(ji-1,jj ) + utau(ji,jj) 494 496 zty2 = vtau(ji ,jj-1) + vtau(ji,jj) … … 547 549 REAL(wp) :: zdku, zdkv, zsqen ! - - 548 550 REAL(wp) :: zemxl, zemlm, zemlp, zmaxice ! - - 549 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zmxlm, zmxld ! 3D workspace551 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zmxlm, zmxld ! 3D workspace 550 552 !!-------------------------------------------------------------------- 551 553 ! … … 560 562 zmxld(:,:,:) = rmxl_min 561 563 ! 564 ! TODO: [tiling] NOT TESTED- requires waves 562 565 IF(ln_sdw .AND. ln_mxhsw) THEN 563 566 zmxlm(:,:,1)= vkarmn * MAX ( 1.6 * hsw(:,:) , 0.02 ) ! surface mixing length = F(wave height) … … 576 579 END_2D 577 580 #else 581 ! TODO: [tiling] NOT TESTED- requires sea ice 578 582 SELECT CASE( nn_mxlice ) ! Type of scaling under sea-ice 579 583 ! … … 698 702 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 699 703 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !* vertical eddy viscosity & diffivity at w-points 700 DO_3D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) !* vertical eddy viscosity & diffivity at w-points704 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) !* vertical eddy viscosity & diffivity at w-points 701 705 zsqen = SQRT( en(ji,jj,jk) ) 702 706 zav = rn_ediff * zmxlm(ji,jj,jk) * zsqen … … 709 713 IF( nn_pdl == 1 ) THEN !* Prandtl number case: update avt 710 714 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 711 DO_3D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )715 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 712 716 p_avt(ji,jj,jk) = MAX( apdlr(ji,jj,jk) * p_avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) 713 717 END_3D … … 825 829 ! 826 830 ! !* Check of some namelist values 827 IF( nn_mxl < 0 .OR. nn_mxl > 3 ) CALL ctl_stop( 'bad flag: nn_mxl is 0, 1 or 2' )828 IF( nn_pdl < 0 .OR. nn_pdl > 1 ) CALL ctl_stop( 'bad flag: nn_pdl is 0 or 1 829 IF( nn_htau < 0 .OR. nn_htau > 1 ) CALL ctl_stop( 'bad flag: nn_htau is 0 , 1 or 2' )831 IF( nn_mxl < 0 .OR. nn_mxl > 3 ) CALL ctl_stop( 'bad flag: nn_mxl is 0, 1, 2 or 3' ) 832 IF( nn_pdl < 0 .OR. nn_pdl > 1 ) CALL ctl_stop( 'bad flag: nn_pdl is 0 or 1' ) 833 IF( nn_htau < 0 .OR. nn_htau > 1 ) CALL ctl_stop( 'bad flag: nn_htau is 0 or 1' ) 830 834 IF( nn_etau == 3 .AND. .NOT. ln_cpl ) CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' ) 831 835 ! … … 835 839 rn_mxl0 = rmxl_min 836 840 ENDIF 837 838 IF( nn_etau == 2 ) CALL zdf_mxl( nit000, Kmm ) ! Initialization of nmln839 840 841 ! !* depth of penetration of surface tke 841 842 IF( nn_etau /= 0 ) THEN
Note: See TracChangeset
for help on using the changeset viewer.
