Changeset 12495 for NEMO/branches/2020/dev_r12472_ASINTER-05_Masson_CurrentFeedback/tests/BENCH/MY_SRC/zdfiwm.F90
- Timestamp:
- 2020-03-02T09:10:34+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r12472_ASINTER-05_Masson_CurrentFeedback/tests/BENCH/MY_SRC/zdfiwm.F90
r12377 r12495 48 48 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hcri_iwm ! decay scale for low-mode critical slope dissipation (m) 49 49 50 !! * Substitutions 51 # include "do_loop_substitute.h90" 50 52 !!---------------------------------------------------------------------- 51 53 !! NEMO/OCE 4.0 , NEMO Consortium (2018) 52 !! $Id: zdfiwm.F90 8093 2017-05-30 08:13:14Z gm$53 !! Software governed by the CeCILL licen ce (./LICENSE)54 !! $Id: zdfiwm.F90 12377 2020-02-12 14:39:06Z acc $ 55 !! Software governed by the CeCILL license (see ./LICENSE) 54 56 !!---------------------------------------------------------------------- 55 57 CONTAINS … … 85 87 !! This is divided into three components: 86 88 !! 1. Bottom-intensified low-mode dissipation at critical slopes 87 !! zemx_iwm(z) = ( ecri_iwm / r au0 ) * EXP( -(H-z)/hcri_iwm )89 !! zemx_iwm(z) = ( ecri_iwm / rho0 ) * EXP( -(H-z)/hcri_iwm ) 88 90 !! / ( 1. - EXP( - H/hcri_iwm ) ) * hcri_iwm 89 91 !! where hcri_iwm is the characteristic length scale of the bottom 90 92 !! intensification, ecri_iwm a map of available power, and H the ocean depth. 91 93 !! 2. Pycnocline-intensified low-mode dissipation 92 !! zemx_iwm(z) = ( epyc_iwm / r au0 ) * ( sqrt(rn2(z))^nn_zpyc )94 !! zemx_iwm(z) = ( epyc_iwm / rho0 ) * ( sqrt(rn2(z))^nn_zpyc ) 93 95 !! / SUM( sqrt(rn2(z))^nn_zpyc * e3w(z) ) 94 96 !! where epyc_iwm is a map of available power, and nn_zpyc … … 96 98 !! energy dissipation. 97 99 !! 3. WKB-height dependent high mode dissipation 98 !! zemx_iwm(z) = ( ebot_iwm / r au0 ) * rn2(z) * EXP(-z_wkb(z)/hbot_iwm)100 !! zemx_iwm(z) = ( ebot_iwm / rho0 ) * rn2(z) * EXP(-z_wkb(z)/hbot_iwm) 99 101 !! / SUM( rn2(z) * EXP(-z_wkb(z)/hbot_iwm) * e3w(z) ) 100 102 !! where hbot_iwm is the characteristic length scale of the WKB bottom … … 121 123 ! 122 124 INTEGER :: ji, jj, jk ! dummy loop indices 123 REAL(wp) :: zztmp ! scalar workspace125 REAL(wp) :: zztmp, ztmp1, ztmp2 ! scalar workspace 124 126 REAL(wp), DIMENSION(jpi,jpj) :: zfact ! Used for vertical structure 125 127 REAL(wp), DIMENSION(jpi,jpj) :: zhdep ! Ocean depth … … 147 149 ! !* Critical slope mixing: distribute energy over the time-varying ocean depth, 148 150 ! using an exponential decay from the seafloor. 149 DO jj = 1, jpj ! part independent of the level 150 DO ji = 1, jpi 151 zhdep(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1) ! depth of the ocean 152 zfact(ji,jj) = rau0 * ( 1._wp - EXP( -zhdep(ji,jj) / hcri_iwm(ji,jj) ) ) 153 IF( zfact(ji,jj) /= 0._wp ) zfact(ji,jj) = ecri_iwm(ji,jj) / zfact(ji,jj) 154 END DO 155 END DO 156 !!gm gde3w ==>>> check for ssh taken into account.... seem OK gde3w_n=gdept(Kmm) - ssh(Kmm) 157 DO jk = 2, jpkm1 ! complete with the level-dependent part 158 zemx_iwm(:,:,jk) = zfact(:,:) * ( EXP( ( gde3w(:,:,jk ) - zhdep(:,:) ) / hcri_iwm(:,:) ) & 159 & - EXP( ( gde3w(:,:,jk-1) - zhdep(:,:) ) / hcri_iwm(:,:) ) ) * wmask(:,:,jk) & 160 & / ( gde3w(:,:,jk) - gde3w(:,:,jk-1) ) 161 162 !!gm delta(gde3w) = e3t(Kmm) !! Please verify the grid-point position w versus t-point 151 DO_2D_11_11 152 zhdep(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1) ! depth of the ocean 153 zfact(ji,jj) = rho0 * ( 1._wp - EXP( -zhdep(ji,jj) / hcri_iwm(ji,jj) ) ) 154 IF( zfact(ji,jj) /= 0._wp ) zfact(ji,jj) = ecri_iwm(ji,jj) / zfact(ji,jj) 155 END_2D 156 !!gm gde3w ==>>> check for ssh taken into account.... seem OK gde3w_n=gdept(:,:,:,Kmm) - ssh(:,:,Kmm) 157 DO_3D_11_11( 2, jpkm1 ) 158 IF ( zfact(ji,jj) == 0._wp .OR. wmask(ji,jj,jk) == 0._wp ) THEN ! optimization 159 zemx_iwm(ji,jj,jk) = 0._wp 160 ELSE 161 zemx_iwm(ji,jj,jk) = zfact(ji,jj) * ( EXP( ( gde3w(ji,jj,jk ) - zhdep(ji,jj) ) / hcri_iwm(ji,jj) ) & 162 & - EXP( ( gde3w(ji,jj,jk-1) - zhdep(ji,jj) ) / hcri_iwm(ji,jj) ) ) & 163 & / ( gde3w(ji,jj,jk) - gde3w(ji,jj,jk-1) ) 164 ENDIF 165 END_3D 166 !!gm delta(gde3w) = e3t(:,:,:,Kmm) !! Please verify the grid-point position w versus t-point 163 167 !!gm it seems to me that only 1/hcri_iwm is used ==> compute it one for all 164 168 165 END DO166 169 167 170 ! !* Pycnocline-intensified mixing: distribute energy over the time-varying … … 177 180 END DO 178 181 ! 179 DO jj = 1, jpj 180 DO ji = 1, jpi 181 IF( zfact(ji,jj) /= 0 ) zfact(ji,jj) = epyc_iwm(ji,jj) / ( rau0 * zfact(ji,jj) ) 182 END DO 183 END DO 182 DO_2D_11_11 183 IF( zfact(ji,jj) /= 0 ) zfact(ji,jj) = epyc_iwm(ji,jj) / ( rho0 * zfact(ji,jj) ) 184 END_2D 184 185 ! 185 186 DO jk = 2, jpkm1 ! complete with the level-dependent part … … 194 195 END DO 195 196 ! 196 DO jj= 1, jpj 197 DO ji = 1, jpi 198 IF( zfact(ji,jj) /= 0 ) zfact(ji,jj) = epyc_iwm(ji,jj) / ( rau0 * zfact(ji,jj) ) 199 END DO 200 END DO 197 DO_2D_11_11 198 IF( zfact(ji,jj) /= 0 ) zfact(ji,jj) = epyc_iwm(ji,jj) / ( rho0 * zfact(ji,jj) ) 199 END_2D 201 200 ! 202 201 DO jk = 2, jpkm1 ! complete with the level-dependent part … … 223 222 !!gm 224 223 ! 225 DO jk = 2, jpkm1 226 DO jj = 1, jpj 227 DO ji = 1, jpi 228 IF( zfact(ji,jj) /= 0 ) zwkb(ji,jj,jk) = zhdep(ji,jj) * ( zfact(ji,jj) - zwkb(ji,jj,jk) ) & 229 & * wmask(ji,jj,jk) / zfact(ji,jj) 230 END DO 231 END DO 232 END DO 224 DO_3D_11_11( 2, jpkm1 ) 225 IF( zfact(ji,jj) /= 0 ) zwkb(ji,jj,jk) = zhdep(ji,jj) * ( zfact(ji,jj) - zwkb(ji,jj,jk) ) & 226 & * wmask(ji,jj,jk) / zfact(ji,jj) 227 END_3D 233 228 zwkb(:,:,1) = zhdep(:,:) * wmask(:,:,1) 234 229 ! 235 zweight(:,:,:) = 0._wp 236 DO jk = 2, jpkm1 237 zweight(:,:,jk) = MAX( 0._wp, rn2(:,:,jk) ) * hbot_iwm(:,:) * wmask(:,:,jk) & 238 & * ( EXP( -zwkb(:,:,jk) / hbot_iwm(:,:) ) - EXP( -zwkb(:,:,jk-1) / hbot_iwm(:,:) ) ) 239 END DO 230 DO_3D_11_11( 2, jpkm1 ) 231 IF ( rn2(ji,jj,jk) <= 0._wp .OR. wmask(ji,jj,jk) == 0._wp ) THEN ! optimization 232 zweight(ji,jj,jk) = 0._wp 233 ELSE 234 zweight(ji,jj,jk) = rn2(ji,jj,jk) * hbot_iwm(ji,jj) & 235 & * ( EXP( -zwkb(ji,jj,jk) / hbot_iwm(ji,jj) ) - EXP( -zwkb(ji,jj,jk-1) / hbot_iwm(ji,jj) ) ) 236 ENDIF 237 END_3D 240 238 ! 241 239 zfact(:,:) = 0._wp … … 244 242 END DO 245 243 ! 246 DO jj = 1, jpj 247 DO ji = 1, jpi 248 IF( zfact(ji,jj) /= 0 ) zfact(ji,jj) = ebot_iwm(ji,jj) / ( rau0 * zfact(ji,jj) ) 249 END DO 250 END DO 244 DO_2D_11_11 245 IF( zfact(ji,jj) /= 0 ) zfact(ji,jj) = ebot_iwm(ji,jj) / ( rho0 * zfact(ji,jj) ) 246 END_2D 251 247 ! 252 248 DO jk = 2, jpkm1 ! complete with the level-dependent part … … 259 255 ! Calculate molecular kinematic viscosity 260 256 znu_t(:,:,:) = 1.e-4_wp * ( 17.91_wp - 0.53810_wp * ts(:,:,:,jp_tem,Kmm) + 0.00694_wp * ts(:,:,:,jp_tem,Kmm) * ts(:,:,:,jp_tem,Kmm) & 261 & + 0.02305_wp * ts(:,:,:,jp_sal,Kmm) ) * tmask(:,:,:) * r1_r au0257 & + 0.02305_wp * ts(:,:,:,jp_sal,Kmm) ) * tmask(:,:,:) * r1_rho0 262 258 DO jk = 2, jpkm1 263 259 znu_w(:,:,jk) = 0.5_wp * ( znu_t(:,:,jk-1) + znu_t(:,:,jk) ) * wmask(:,:,jk) … … 276 272 ! 277 273 IF( ln_mevar ) THEN ! Variable mixing efficiency case : modify zav_wave in the 278 DO jk = 2, jpkm1 ! energetic (Reb > 480) and buoyancy-controlled (Reb <10.224 ) regimes 279 DO jj = 1, jpj 280 DO ji = 1, jpi 281 IF( zReb(ji,jj,jk) > 480.00_wp ) THEN 282 zav_wave(ji,jj,jk) = 3.6515_wp * znu_w(ji,jj,jk) * SQRT( zReb(ji,jj,jk) ) 283 ELSEIF( zReb(ji,jj,jk) < 10.224_wp ) THEN 284 zav_wave(ji,jj,jk) = 0.052125_wp * znu_w(ji,jj,jk) * zReb(ji,jj,jk) * SQRT( zReb(ji,jj,jk) ) 285 ENDIF 286 END DO 287 END DO 288 END DO 274 DO_3D_11_11( 2, jpkm1 ) 275 IF( zReb(ji,jj,jk) > 480.00_wp ) THEN 276 zav_wave(ji,jj,jk) = 3.6515_wp * znu_w(ji,jj,jk) * SQRT( zReb(ji,jj,jk) ) 277 ELSEIF( zReb(ji,jj,jk) < 10.224_wp ) THEN 278 zav_wave(ji,jj,jk) = 0.052125_wp * znu_w(ji,jj,jk) * zReb(ji,jj,jk) * SQRT( zReb(ji,jj,jk) ) 279 ENDIF 280 END_3D 289 281 ENDIF 290 282 ! … … 296 288 zztmp = 0._wp 297 289 !!gm used of glosum 3D.... 298 DO jk = 2, jpkm1 299 DO jj = 1, jpj 300 DO ji = 1, jpi 301 zztmp = zztmp + e3w(ji,jj,jk,Kmm) * e1e2t(ji,jj) & 302 & * MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj) 303 END DO 304 END DO 305 END DO 290 DO_3D_11_11( 2, jpkm1 ) 291 zztmp = zztmp + e3w(ji,jj,jk,Kmm) * e1e2t(ji,jj) & 292 & * MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj) 293 END_3D 306 294 CALL mpp_sum( 'zdfiwm', zztmp ) 307 zztmp = r au0 * zztmp ! Global integral of rauo * Kz * N^2 = power contributing to mixing295 zztmp = rho0 * zztmp ! Global integral of rauo * Kz * N^2 = power contributing to mixing 308 296 ! 309 297 IF(lwp) THEN … … 321 309 ! 322 310 IF( ln_tsdiff ) THEN !* Option for differential mixing of salinity and temperature 323 DO jk = 2, jpkm1 ! Calculate S/T diffusivity ratio as a function of Reb324 DO jj = 1, jpj325 DO ji = 1, jpi326 zav_ratio(ji,jj,jk) = ( 0.505_wp + 0.495_wp * &327 & TANH( 0.92_wp * ( LOG10( MAX( 1.e-20_wp, zReb(ji,jj,jk) * 5._wp * r1_6 ) ) - 0.60_wp ) ) &328 & ) * wmask(ji,jj,jk)329 END DO330 END DO331 END DO311 ztmp1 = 0.505_wp + 0.495_wp * TANH( 0.92_wp * ( LOG10( 1.e-20_wp ) - 0.60_wp ) ) 312 DO_3D_11_11( 2, jpkm1 ) 313 ztmp2 = zReb(ji,jj,jk) * 5._wp * r1_6 314 IF ( ztmp2 > 1.e-20_wp .AND. wmask(ji,jj,jk) == 1._wp ) THEN 315 zav_ratio(ji,jj,jk) = 0.505_wp + 0.495_wp * TANH( 0.92_wp * ( LOG10(ztmp2) - 0.60_wp ) ) 316 ELSE 317 zav_ratio(ji,jj,jk) = ztmp1 * wmask(ji,jj,jk) 318 ENDIF 319 END_3D 332 320 CALL iom_put( "av_ratio", zav_ratio ) 333 321 DO jk = 2, jpkm1 !* update momentum & tracer diffusivity with wave-driven mixing … … 349 337 !* output useful diagnostics: Kz*N^2 , 350 338 !!gm Kz*N2 should take into account the ratio avs/avt if it is used.... (see diaar5) 351 ! vertical integral of r au0 * Kz * N^2 , energy density (zemx_iwm)339 ! vertical integral of rho0 * Kz * N^2 , energy density (zemx_iwm) 352 340 IF( iom_use("bflx_iwm") .OR. iom_use("pcmap_iwm") ) THEN 353 341 ALLOCATE( z2d(jpi,jpj) , z3d(jpi,jpj,jpk) ) … … 357 345 z2d(:,:) = z2d(:,:) + e3w(:,:,jk,Kmm) * z3d(:,:,jk) * wmask(:,:,jk) 358 346 END DO 359 z2d(:,:) = r au0 * z2d(:,:)347 z2d(:,:) = rho0 * z2d(:,:) 360 348 CALL iom_put( "bflx_iwm", z3d ) 361 349 CALL iom_put( "pcmap_iwm", z2d ) … … 364 352 CALL iom_put( "emix_iwm", zemx_iwm ) 365 353 366 IF( ln_ctl) CALL prt_ctl(tab3d_1=zav_wave , clinfo1=' iwm - av_wave: ', tab3d_2=avt, clinfo2=' avt: ', kdim=jpk)354 IF(sn_cfctl%l_prtctl) CALL prt_ctl(tab3d_1=zav_wave , clinfo1=' iwm - av_wave: ', tab3d_2=avt, clinfo2=' avt: ', kdim=jpk) 367 355 ! 368 356 END SUBROUTINE zdf_iwm … … 395 383 !! de Lavergne et al. in prep., 2017 396 384 !!---------------------------------------------------------------------- 397 INTEGER :: ji, jj, jk ! dummy loop indices398 385 INTEGER :: inum ! local integer 399 386 INTEGER :: ios
Note: See TracChangeset
for help on using the changeset viewer.
