Changeset 15517 for NEMO/branches/2021
- Timestamp:
- 2021-11-16T15:48:25+01:00 (3 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/dev_r15388_updated_zdfiwm/src/OCE/ZDF/zdfiwm.F90
r15466 r15517 42 42 LOGICAL :: ln_tsdiff ! account for differential T/S wave-driven mixing (=T) or not (=F) 43 43 44 REAL(wp):: r1_6 = 1._wp / 6._wp45 REAL(wp):: rnu = 1.4e-6_wp ! molecular kinematic viscosity44 REAL(wp):: z1_6 = 1._wp / 6._wp 45 REAL(wp):: znu = 1.4e-6_wp ! molecular kinematic viscosity 46 46 47 47 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ebot_iwm ! bottom-intensified dissipation above abyssal hills (W/m2) … … 135 135 INTEGER :: ji, jj, jk ! dummy loop indices 136 136 REAL(wp), SAVE :: zztmp 137 REAL(wp) :: ztmp1, ztmp2 ! scalar workspace137 ! 138 138 REAL(wp), DIMENSION(A2D(nn_hls)) :: zfact ! Used for vertical structure 139 139 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zReb ! Turbulence intensity parameter … … 145 145 !!---------------------------------------------------------------------- 146 146 ! 147 ! !* Set to zero the 1st and last vertical levels of appropriate variables 148 IF( iom_use("emix_iwm") ) THEN 149 zemx_iwm(:,:,:) = 0._wp 150 ENDIF 151 IF( iom_use("av_ratio") ) THEN 152 zav_ratio(:,:,:) = 0._wp 153 ENDIF 154 IF( iom_use("av_wave") .OR. sn_cfctl%l_prtctl ) THEN 155 zav_wave(:,:,:) = 0._wp 156 ENDIF 147 ! !* Initialize appropriately certain variables 148 zav_ratio(:,:,:) = 1._wp * wmask(:,:,:) ! important to set it to 1 here 149 IF( iom_use("emix_iwm") ) zemx_iwm (:,:,:) = 0._wp 150 IF( iom_use("av_wave") .OR. sn_cfctl%l_prtctl ) zav_wave (:,:,:) = 0._wp 157 151 ! 158 152 ! ! ----------------------------- ! … … 163 157 ! !* ocean depth using an exponential decay from the seafloor. 164 158 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! part independent of the level 165 zfact(ji,jj) = rho0 * ( 1._wp - EXP( -ht(ji,jj) * hcri_iwm(ji,jj) ) ) 166 IF( zfact(ji,jj) /= 0._wp ) zfact(ji,jj) = ecri_iwm(ji,jj) / zfact(ji,jj) 159 IF( ht(ji,jj) /= 0._wp ) THEN ; zfact(ji,jj) = ecri_iwm(ji,jj) * r1_rho0 / ( 1._wp - EXP( -ht(ji,jj) * hcri_iwm(ji,jj) ) ) 160 ELSE ; zfact(ji,jj) = 0._wp 161 ENDIF 167 162 END_2D 168 163 169 164 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! complete with the level-dependent part 170 IF ( zfact(ji,jj) == 0._wp .OR. wmask(ji,jj,jk) == 0._wp ) THEN ! optimization 171 zemx_iwm(ji,jj,jk) = 0._wp 172 ELSE 173 zemx_iwm(ji,jj,jk) = zfact(ji,jj) * ( EXP( ( gdept(ji,jj,jk ,Kmm) - ht(ji,jj) ) * hcri_iwm(ji,jj) ) & 174 & - EXP( ( gdept(ji,jj,jk-1,Kmm) - ht(ji,jj) ) * hcri_iwm(ji,jj) ) ) & 175 & / e3w(ji,jj,jk,Kmm) 176 ENDIF 165 zemx_iwm(ji,jj,jk) = zfact(ji,jj) * ( EXP( ( gdept(ji,jj,jk ,Kmm) - ht(ji,jj) ) * hcri_iwm(ji,jj) ) & 166 & - EXP( ( gdept(ji,jj,jk-1,Kmm) - ht(ji,jj) ) * hcri_iwm(ji,jj) ) & 167 & ) * wmask(ji,jj,jk) / e3w(ji,jj,jk,Kmm) 177 168 END_3D 178 169 179 170 !* 'bot' component: distribute energy over the time-varying 180 171 !* ocean depth using an algebraic decay above the seafloor. 181 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )182 zfact(ji,jj) = 0._wp183 END_2D184 172 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! part independent of the level 185 IF( ht(ji,jj) /= 0._wp ) zfact(ji,jj) = ebot_iwm(ji,jj) * ( 1._wp + hbot_iwm(ji,jj) / ht(ji,jj) ) / rho0 173 IF( ht(ji,jj) /= 0._wp ) THEN ; zfact(ji,jj) = ebot_iwm(ji,jj) * ( 1._wp + hbot_iwm(ji,jj) / ht(ji,jj) ) * r1_rho0 174 ELSE ; zfact(ji,jj) = 0._wp 175 ENDIF 186 176 END_2D 187 177 188 178 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! complete with the level-dependent part 189 zemx_iwm(ji,jj,jk) = zemx_iwm(ji,jj,jk) + 190 & zfact(ji,jj) * ( 1._wp / ( 1._wp + ( ht(ji,jj) - gdept(ji,jj,jk ,Kmm) ) / hbot_iwm(ji,jj) )&191 & - 1._wp / ( 1._wp + ( ht(ji,jj) - gdept(ji,jj,jk-1,Kmm) ) / hbot_iwm(ji,jj) ) ) * wmask(ji,jj,jk) &192 & / e3w(ji,jj,jk,Kmm)179 zemx_iwm(ji,jj,jk) = zemx_iwm(ji,jj,jk) + & 180 & zfact(ji,jj) * ( 1._wp / ( 1._wp + ( ht(ji,jj) - gdept(ji,jj,jk ,Kmm) ) / hbot_iwm(ji,jj) ) & 181 & - 1._wp / ( 1._wp + ( ht(ji,jj) - gdept(ji,jj,jk-1,Kmm) ) / hbot_iwm(ji,jj) ) & 182 & ) * wmask(ji,jj,jk) / e3w(ji,jj,jk,Kmm) 193 183 END_3D 194 184 … … 203 193 ! 204 194 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 205 IF( zfact(ji,jj) /= 0._wp ) zfact(ji,jj) = ensq_iwm(ji,jj) / ( rho0 * zfact(ji,jj))195 IF( zfact(ji,jj) /= 0._wp ) zfact(ji,jj) = ensq_iwm(ji,jj) * r1_rho0 / zfact(ji,jj) 206 196 END_2D 207 197 ! … … 216 206 END_2D 217 207 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! part independent of the level 218 zfact(ji,jj) = zfact(ji,jj) + e3w(ji,jj,jk,Kmm) * SQRT( MAX( 0._wp, rn2(ji,jj,jk) ))208 zfact(ji,jj) = zfact(ji,jj) + e3w(ji,jj,jk,Kmm) * SQRT( MAX( 0._wp, rn2(ji,jj,jk) ) ) 219 209 END_3D 220 210 ! 221 211 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 222 IF( zfact(ji,jj) /= 0._wp ) zfact(ji,jj) = esho_iwm(ji,jj) / ( rho0 * zfact(ji,jj))212 IF( zfact(ji,jj) /= 0._wp ) zfact(ji,jj) = esho_iwm(ji,jj) * r1_rho0 / zfact(ji,jj) 223 213 END_2D 224 214 ! 225 215 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! complete with the level-dependent part 226 zemx_iwm(ji,jj,jk) = zemx_iwm(ji,jj,jk) + zfact(ji,jj) * SQRT( MAX( 0._wp, rn2(ji,jj,jk) ))216 zemx_iwm(ji,jj,jk) = zemx_iwm(ji,jj,jk) + zfact(ji,jj) * SQRT( MAX( 0._wp, rn2(ji,jj,jk) ) ) 227 217 END_3D 228 218 229 219 ! Calculate turbulence intensity parameter Reb 230 220 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 231 zReb(ji,jj,jk) = zemx_iwm(ji,jj,jk) / MAX( 1.e-20_wp, rnu * rn2(ji,jj,jk) )221 zReb(ji,jj,jk) = zemx_iwm(ji,jj,jk) / MAX( 1.e-20_wp, znu * rn2(ji,jj,jk) ) 232 222 END_3D 233 223 ! 234 224 ! Define internal wave-induced diffusivity 235 225 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 236 zav_wave(ji,jj,jk) = zReb(ji,jj,jk) * r1_6 * rnu ! This corresponds to a constant mixing efficiency of 1/6226 zav_wave(ji,jj,jk) = zReb(ji,jj,jk) * z1_6 * znu ! This corresponds to a constant mixing efficiency of 1/6 237 227 END_3D 238 228 ! … … 240 230 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! energetic (Reb > 480) and buoyancy-controlled (Reb <10.224) regimes 241 231 IF( zReb(ji,jj,jk) > 480.00_wp ) THEN 242 zav_wave(ji,jj,jk) = 3.6515_wp * rnu * SQRT( zReb(ji,jj,jk) )232 zav_wave(ji,jj,jk) = 3.6515_wp * znu * SQRT( zReb(ji,jj,jk) ) 243 233 ELSEIF( zReb(ji,jj,jk) < 10.224_wp ) THEN 244 zav_wave(ji,jj,jk) = 0.052125_wp * rnu * zReb(ji,jj,jk) * SQRT( zReb(ji,jj,jk) )234 zav_wave(ji,jj,jk) = 0.052125_wp * znu * zReb(ji,jj,jk) * SQRT( zReb(ji,jj,jk) ) 245 235 ENDIF 246 236 END_3D … … 248 238 ! 249 239 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! Bound diffusivity by molecular value and 100 cm2/s 250 zav_wave(ji,jj,jk) = MIN( MAX( 1.4e-7_wp, zav_wave(ji,jj,jk) ), 1.e-2_wp ) * wmask(ji,jj,jk) 251 END_3D 240 zav_wave(ji,jj,jk) = MIN( MAX( 1.4e-7_wp, zav_wave(ji,jj,jk) ), 1.e-2_wp ) * wmask(ji,jj,jk) 241 END_3D 242 ! 243 ! ! ----------------------- ! 244 ! ! Update mixing coefs ! 245 ! ! ----------------------- ! 246 ! 247 IF( ln_tsdiff ) THEN !* Option for differential mixing of salinity and temperature 248 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! Calculate S/T diffusivity ratio as a function of Reb (else it is set to 1) 249 zav_ratio(ji,jj,jk) = ( 0.505_wp + & 250 & 0.495_wp * TANH( 0.92_wp * ( LOG10( MAX( 1.e-20, zReb(ji,jj,jk) * 5._wp * z1_6 ) ) - 0.60_wp ) ) & 251 & ) * wmask(ji,jj,jk) 252 END_3D 253 ENDIF 254 CALL iom_put( "av_ratio", zav_ratio ) 255 ! 256 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 257 p_avs(ji,jj,jk) = p_avs(ji,jj,jk) + zav_wave(ji,jj,jk) * zav_ratio(ji,jj,jk) 258 p_avt(ji,jj,jk) = p_avt(ji,jj,jk) + zav_wave(ji,jj,jk) 259 p_avm(ji,jj,jk) = p_avm(ji,jj,jk) + zav_wave(ji,jj,jk) 260 END_3D 261 ! !* output internal wave-driven mixing coefficient 262 CALL iom_put( "av_wave", zav_wave ) 263 !* output useful diagnostics: Kz*N^2 , 264 ! vertical integral of rho0 * Kz * N^2 , energy density (zemx_iwm) 265 IF( iom_use("bflx_iwm") .OR. iom_use("pcmap_iwm") ) THEN 266 ALLOCATE( z2d(A2D(nn_hls)) , z3d(A2D(nn_hls),jpk) ) 267 z2d(:,:) = 0._wp ; z3d(:,:,:) = 0._wp ! Initialisation for iom_put 268 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 269 z3d(ji,jj,jk) = MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk) 270 z2d(ji,jj) = z2d(ji,jj) + rho0 * e3w(ji,jj,jk,Kmm) * z3d(ji,jj,jk) * wmask(ji,jj,jk) 271 END_3D 272 CALL iom_put( "bflx_iwm", z3d ) 273 CALL iom_put( "pcmap_iwm", z2d ) 274 DEALLOCATE( z2d , z3d ) 275 ENDIF 276 CALL iom_put( "emix_iwm", zemx_iwm ) 277 252 278 ! 253 279 IF( kt == nit000 ) THEN !* Control print at first time-step: diagnose the energy consumed by zav_wave … … 260 286 IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN ! Do only on the last tile 261 287 CALL mpp_sum( 'zdfiwm', zztmp ) 262 zztmp = rho0 * zztmp ! Global integral of r auo* Kz * N^2 = power contributing to mixing288 zztmp = rho0 * zztmp ! Global integral of rho0 * Kz * N^2 = power contributing to mixing 263 289 ! 264 290 IF(lwp) THEN … … 272 298 ENDIF 273 299 274 ! ! ----------------------- !275 ! ! Update mixing coefs !276 ! ! ----------------------- !277 !278 IF( ln_tsdiff ) THEN !* Option for differential mixing of salinity and temperature279 ztmp1 = 0.505_wp + 0.495_wp * TANH( 0.92_wp * ( LOG10( 1.e-20_wp ) - 0.60_wp ) )280 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! Calculate S/T diffusivity ratio as a function of Reb281 ztmp2 = zReb(ji,jj,jk) * 5._wp * r1_6282 IF ( ztmp2 > 1.e-20_wp .AND. wmask(ji,jj,jk) == 1._wp ) THEN283 zav_ratio(ji,jj,jk) = 0.505_wp + 0.495_wp * TANH( 0.92_wp * ( LOG10(ztmp2) - 0.60_wp ) )284 ELSE285 zav_ratio(ji,jj,jk) = ztmp1 * wmask(ji,jj,jk)286 ENDIF287 END_3D288 CALL iom_put( "av_ratio", zav_ratio )289 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) !* update momentum & tracer diffusivity with wave-driven mixing290 p_avs(ji,jj,jk) = p_avs(ji,jj,jk) + zav_wave(ji,jj,jk) * zav_ratio(ji,jj,jk)291 p_avt(ji,jj,jk) = p_avt(ji,jj,jk) + zav_wave(ji,jj,jk)292 p_avm(ji,jj,jk) = p_avm(ji,jj,jk) + zav_wave(ji,jj,jk)293 END_3D294 !295 ELSE !* update momentum & tracer diffusivity with wave-driven mixing296 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )297 p_avs(ji,jj,jk) = p_avs(ji,jj,jk) + zav_wave(ji,jj,jk)298 p_avt(ji,jj,jk) = p_avt(ji,jj,jk) + zav_wave(ji,jj,jk)299 p_avm(ji,jj,jk) = p_avm(ji,jj,jk) + zav_wave(ji,jj,jk)300 END_3D301 ENDIF302 ! !* output internal wave-driven mixing coefficient303 CALL iom_put( "av_wave", zav_wave )304 !* output useful diagnostics: Kz*N^2 ,305 ! vertical integral of rho0 * Kz * N^2 , energy density (zemx_iwm)306 IF( iom_use("bflx_iwm") .OR. iom_use("pcmap_iwm") ) THEN307 ALLOCATE( z2d(A2D(nn_hls)) , z3d(A2D(nn_hls),jpk) )308 ! Initialisation for iom_put309 z2d(:,:) = 0._wp ; z3d(:,:,:) = 0._wp310 311 DO_3D( 0, 0, 0, 0, 2, jpkm1 )312 z3d(ji,jj,jk) = MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk)313 z2d(ji,jj) = z2d(ji,jj) + e3w(ji,jj,jk,Kmm) * z3d(ji,jj,jk) * wmask(ji,jj,jk)314 END_3D315 DO_2D( 0, 0, 0, 0 )316 z2d(ji,jj) = rho0 * z2d(ji,jj)317 END_2D318 CALL iom_put( "bflx_iwm", z3d )319 CALL iom_put( "pcmap_iwm", z2d )320 DEALLOCATE( z2d , z3d )321 ENDIF322 CALL iom_put( "emix_iwm", zemx_iwm )323 324 300 IF(sn_cfctl%l_prtctl) CALL prt_ctl(tab3d_1=zav_wave , clinfo1=' iwm - av_wave: ', tab3d_2=avt, clinfo2=' avt: ', kdim=jpk) 325 301 ! … … 356 332 INTEGER :: inum ! local integer 357 333 INTEGER :: ios 358 REAL(wp) :: zbot, zcri, znsq, zsho ! local scalars359 334 ! 360 335 CHARACTER(len=256) :: cn_dir ! Root directory for location of ssr files … … 372 347 TYPE(FLD ), DIMENSION(jpiwm) :: sf_iwm ! structure of input fields (file informations, fields read) 373 348 ! 349 REAL(wp), DIMENSION(jpi,jpj,4) :: ztmp 350 REAL(wp), DIMENSION(4) :: zdia 351 ! 374 352 NAMELIST/namzdf_iwm/ ln_mevar, ln_tsdiff, & 375 353 & cn_dir, sn_mpb, sn_mpc, sn_mpn, sn_mps, sn_dsb, sn_dsc … … 395 373 ! ensures that avt remains larger than its molecular value (=1.4e-7). Therefore, avtb should 396 374 ! be set here to a very small value, and avmb to its (uniform) molecular value (=1.4e-6). 397 avmb(:) = rnu ! molecular value375 avmb(:) = znu ! molecular value 398 376 avtb(:) = 1.e-10_wp ! very small diffusive minimum (background avt is specified in zdf_iwm) 399 avtb_2d(:,:) = 1._wp ! uniform377 avtb_2d(:,:) = 1._wp ! uniform 400 378 IF(lwp) THEN ! Control print 401 379 WRITE(numout,*) … … 422 400 sf_iwm(jp_mpb)%fnow(:,:,1) = 1.e-10_wp 423 401 sf_iwm(jp_mpc)%fnow(:,:,1) = 1.e-10_wp 424 sf_iwm(jp_mpn)%fnow(:,:,1) = 1.e- 6_wp402 sf_iwm(jp_mpn)%fnow(:,:,1) = 1.e-5_wp 425 403 sf_iwm(jp_mps)%fnow(:,:,1) = 1.e-10_wp 426 404 sf_iwm(jp_dsb)%fnow(:,:,1) = 100._wp … … 439 417 hcri_iwm(:,:) = 1._wp / hcri_iwm(:,:) ! only the inverse height is used, hence calculated here once for all 440 418 441 zbot = glob_sum( 'zdfiwm', e1e2t(:,:) * ebot_iwm(:,:) ) 442 zcri = glob_sum( 'zdfiwm', e1e2t(:,:) * ecri_iwm(:,:) ) 443 znsq = glob_sum( 'zdfiwm', e1e2t(:,:) * ensq_iwm(:,:) ) 444 zsho = glob_sum( 'zdfiwm', e1e2t(:,:) * esho_iwm(:,:) ) 419 ! diags 420 ztmp(:,:,1) = e1e2t(:,:) * ebot_iwm(:,:) 421 ztmp(:,:,2) = e1e2t(:,:) * ecri_iwm(:,:) 422 ztmp(:,:,3) = e1e2t(:,:) * ensq_iwm(:,:) 423 ztmp(:,:,4) = e1e2t(:,:) * esho_iwm(:,:) 424 425 zdia(1:4) = glob_sum_vec( 'zdfiwm', ztmp(:,:,1:4) ) 445 426 446 427 IF(lwp) THEN 447 WRITE(numout,*) ' Dissipation above abyssal hills: ', z bot* 1.e-12_wp, 'TW'448 WRITE(numout,*) ' Dissipation along topographic slopes: ', z cri* 1.e-12_wp, 'TW'449 WRITE(numout,*) ' Dissipation scaling with N^2: ', z nsq* 1.e-12_wp, 'TW'450 WRITE(numout,*) ' Dissipation due to shoaling: ', z sho* 1.e-12_wp, 'TW'428 WRITE(numout,*) ' Dissipation above abyssal hills: ', zdia(1) * 1.e-12_wp, 'TW' 429 WRITE(numout,*) ' Dissipation along topographic slopes: ', zdia(2) * 1.e-12_wp, 'TW' 430 WRITE(numout,*) ' Dissipation scaling with N^2: ', zdia(3) * 1.e-12_wp, 'TW' 431 WRITE(numout,*) ' Dissipation due to shoaling: ', zdia(4) * 1.e-12_wp, 'TW' 451 432 ENDIF 452 433 !
Note: See TracChangeset
for help on using the changeset viewer.
