- Timestamp:
- 2021-05-06T19:58:34+02:00 (3 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/OCE/ZDF/zdfosm.F90
r14798 r14802 158 158 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: shol ! Stability parameter for boundary layer 159 159 ! 160 ! Layer averages: BL 161 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: av_t_bl ! Temperature average 162 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: av_s_bl ! Salinity average 163 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: av_u_bl ! Velocity average (u) 164 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: av_v_bl ! Velocity average (v) 165 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: av_b_bl ! Buoyancy 166 ! 167 ! Layer averages: ML 168 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: av_t_ml ! Temperature average 169 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: av_s_ml ! Salinity average 170 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: av_u_ml ! Velocity average (u) 171 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: av_v_ml ! Velocity average (v) 172 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: av_b_ml ! Buoyancy 173 ! 174 ! Layer averages: MLE 175 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: av_t_mle ! Temperature average 176 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: av_s_mle ! Salinity average 177 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: av_u_mle ! Velocity average (u) 178 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: av_v_mle ! Velocity average (v) 179 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: av_b_mle ! Buoyancy 180 ! 160 181 ! ** Namelist namzdf_osm ** 161 182 LOGICAL :: ln_use_osm_la ! Use namelist rn_osm_la … … 216 237 !!---------------------------------------------------------------------- 217 238 CONTAINS 218 ! 239 219 240 INTEGER FUNCTION zdf_osm_alloc() 220 241 !!---------------------------------------------------------------------- … … 246 267 zdf_osm_alloc = zdf_osm_alloc + ierr 247 268 ! 269 ALLOCATE( av_t_bl(jpi,jpj), av_s_bl(jpi,jpj), av_u_bl(jpi,jpj), av_v_bl(jpi,jpj), av_b_bl(jpi,jpj), STAT=ierr) 270 zdf_osm_alloc = zdf_osm_alloc + ierr 271 ! 272 ALLOCATE( av_t_ml(jpi,jpj), av_s_ml(jpi,jpj), av_u_ml(jpi,jpj), av_v_ml(jpi,jpj), av_b_ml(jpi,jpj), STAT=ierr) 273 zdf_osm_alloc = zdf_osm_alloc + ierr 274 ! 275 ALLOCATE( av_t_mle(jpi,jpj), av_s_mle(jpi,jpj), av_u_mle(jpi,jpj), av_v_mle(jpi,jpj), av_b_mle(jpi,jpj), STAT=ierr) 276 zdf_osm_alloc = zdf_osm_alloc + ierr 277 ! 248 278 CALL mpp_sum ( 'zdfosm', zdf_osm_alloc ) 249 279 IF( zdf_osm_alloc /= 0 ) CALL ctl_warn( 'zdf_osm_alloc: failed to allocate zdf_osm arrays' ) 250 280 ! 251 281 END FUNCTION zdf_osm_alloc 252 ! 282 253 283 SUBROUTINE zdf_osm( kt, Kbb, Kmm, Krhs, p_avm, & 254 284 & p_avt ) … … 287 317 !! the equation number. (LMD94, here after) 288 318 !!---------------------------------------------------------------------- 289 INTEGER , INTENT(in ) :: kt ! ocean time step 290 INTEGER , INTENT(in ) :: Kbb, Kmm, Krhs ! ocean time level indices 291 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: p_avm, p_avt ! momentum and tracer Kz (w-points) 292 !! 293 INTEGER :: ji, jj, jk, jkflt ! dummy loop indices 294 295 INTEGER :: jl ! dummy loop indices 296 297 INTEGER :: ikbot, jkm1, jkp2 ! 298 299 REAL(wp) :: ztx, zty, zflageos, zstabl, zbuofdep,zucube ! 300 REAL(wp) :: zbeta, zthermal ! 301 REAL(wp) :: zehat, zeta, zhrib, zsig, zscale, zwst, zws, zwm ! Velocity scales 302 REAL(wp) :: zwsun, zwmun, zcons, zconm, zwcons, zwconm ! 303 REAL(wp) :: zsr, zbw, ze, zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zcomp , zrhd,zrhdr,zbvzed ! In situ density 304 INTEGER :: jm ! dummy loop indices 305 REAL(wp) :: zr1, zr2, zr3, zr4, zrhop ! Compression terms 306 REAL(wp) :: zflag, zrn2, zdep21, zdep32, zdep43 307 REAL(wp) :: zesh2, zri, zfri ! Interior richardson mixing 308 REAL(wp) :: zdelta, zdelta2, zdzup, zdzdn, zdzh, zvath, zgat1, zdat1, zkm1m, zkm1t 309 REAL(wp) :: zt,zs,zu,zv,zrh ! variables used in constructing averages 319 INTEGER , INTENT(in ) :: kt ! Ocean time step 320 INTEGER , INTENT(in ) :: Kbb, Kmm, Krhs ! Ocean time level indices 321 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: p_avm, p_avt ! Momentum and tracer Kz (w-points) 322 ! 323 ! Local variables 324 INTEGER :: ji, jj, jk, jl, jm, jkflt ! Dummy loop indices 325 ! 326 REAL(wp) :: zthermal, zbeta 327 REAL(wp) :: zesh2, zri, zfri ! Interior Richardson mixing 328 ! 310 329 ! Scales 311 REAL(wp), DIMENSION(A2D(0)) :: zrad0 ! Surface solar temperature flux (deg m/s) 312 REAL(wp), DIMENSION(A2D(0)) :: zradh ! Radiative flux at bl base (Buoyancy units) 313 REAL(wp) :: zradav ! Radiative flux, bl average (Buoyancy Units) 314 REAL(wp) :: zvw0 ! Surface v-momentum flux 315 REAL(wp), DIMENSION(A2D(0)) :: zwb0tot ! Total surface buoyancy flux including insolation 316 REAL(wp), DIMENSION(jpi,jpj) :: zwb_ent ! Buoyancy entrainment flux 317 REAL(wp), DIMENSION(jpi,jpj) :: zwb_min 318 319 REAL(wp), DIMENSION(jpi,jpj) :: zwb_fk_b ! MLE buoyancy flux averaged over OSBL 320 REAL(wp), DIMENSION(jpi,jpj) :: zwb_fk ! max MLE buoyancy flux 321 REAL(wp), DIMENSION(jpi,jpj) :: zdiff_mle ! extra MLE vertical diff 322 REAL(wp), DIMENSION(jpi,jpj) :: zvel_mle ! velocity scale for dhdt with stable ML and FK 323 330 REAL(wp), DIMENSION(A2D(0)) :: zrad0 ! Surface solar temperature flux (deg m/s) 331 REAL(wp), DIMENSION(A2D(0)) :: zradh ! Radiative flux at bl base (Buoyancy units) 332 REAL(wp) :: zradav ! Radiative flux, bl average (Buoyancy Units) 333 REAL(wp) :: zvw0 ! Surface v-momentum flux 334 REAL(wp), DIMENSION(A2D(0)) :: zwb0tot ! Total surface buoyancy flux including insolation 335 REAL(wp), DIMENSION(jpi,jpj) :: zwb_ent ! Buoyancy entrainment flux 336 REAL(wp), DIMENSION(jpi,jpj) :: zwb_min 337 REAL(wp), DIMENSION(jpi,jpj) :: zwb_fk_b ! MLE buoyancy flux averaged over OSBL 338 REAL(wp), DIMENSION(jpi,jpj) :: zwb_fk ! Max MLE buoyancy flux 339 REAL(wp), DIMENSION(jpi,jpj) :: zdiff_mle ! Extra MLE vertical diff 340 REAL(wp), DIMENSION(jpi,jpj) :: zvel_mle ! Velocity scale for dhdt with stable ML and FK 341 ! 324 342 ! mixed-layer variables 325 326 INTEGER, DIMENSION(jpi,jpj) :: jp_ext ! offset for external level 327 328 REAL(wp), DIMENSION(jpi,jpj) :: zhbl ! bl depth - grid 329 REAL(wp), DIMENSION(jpi,jpj) :: zhml ! ml depth - grid 330 331 REAL(wp), DIMENSION(jpi,jpj) :: zhmle ! MLE depth - grid 332 REAL(wp), DIMENSION(jpi,jpj) :: zmld ! ML depth on grid 333 334 REAL(wp), DIMENSION(jpi,jpj) :: zdh ! pycnocline depth - grid 335 REAL(wp), DIMENSION(A2D(0)) :: zdhdt ! BL depth tendency 336 REAL(wp), DIMENSION(A2D(0)) :: zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext ! external temperature/salinity and buoyancy gradients 337 REAL(wp), DIMENSION(jpi,jpj) :: zdtdx, zdtdy, zdsdx, zdsdy ! horizontal gradients for Fox-Kemper parametrization. 338 339 REAL(wp), DIMENSION(jpi,jpj) :: zt_bl,zs_bl,zu_bl,zv_bl,zb_bl ! averages over the depth of the blayer 340 REAL(wp), DIMENSION(jpi,jpj) :: zt_ml,zs_ml,zu_ml,zv_ml,zb_ml ! averages over the depth of the mixed layer 341 REAL(wp), DIMENSION(jpi,jpj) :: zt_mle,zs_mle,zu_mle,zv_mle,zb_mle ! averages over the depth of the MLE layer 342 REAL(wp), DIMENSION(jpi,jpj) :: zdt_bl,zds_bl,zdu_bl,zdv_bl,zdb_bl ! difference between blayer average and parameter at base of blayer 343 REAL(wp), DIMENSION(jpi,jpj) :: zdt_ml,zds_ml,zdu_ml,zdv_ml,zdb_ml ! difference between mixed layer average and parameter at base of blayer 344 ! REAL(wp), DIMENSION(jpi,jpj) :: zwth_ent,zws_ent ! heat and salinity fluxes at the top of the pycnocline 345 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdbdz_pyc ! parametrised gradient of buoyancy in the pycnocline 346 REAL(wp), DIMENSION(jpi,jpj) :: zdbds_mle ! Magnitude of horizontal buoyancy gradient. 343 INTEGER, DIMENSION(jpi,jpj) :: jp_ext ! Offset for external level 344 ! 345 REAL(wp), DIMENSION(jpi,jpj) :: zhbl ! BL depth - grid 346 REAL(wp), DIMENSION(jpi,jpj) :: zhml ! ML depth - grid 347 ! 348 REAL(wp), DIMENSION(jpi,jpj) :: zhmle ! MLE depth - grid 349 REAL(wp), DIMENSION(jpi,jpj) :: zmld ! ML depth on grid 350 ! 351 REAL(wp), DIMENSION(jpi,jpj) :: zdh ! Pycnocline depth - grid 352 REAL(wp), DIMENSION(A2D(0)) :: zdhdt ! BL depth tendency 353 REAL(wp), DIMENSION(A2D(0)) :: zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext ! External temperature/salinity and buoyancy gradients 354 REAL(wp), DIMENSION(jpi,jpj) :: zdtdx, zdtdy, zdsdx, zdsdy ! Horizontal gradients for Fox-Kemper parametrization. 355 ! 356 REAL(wp), DIMENSION(jpi,jpj) :: zdt_bl,zds_bl,zdu_bl,zdv_bl,zdb_bl ! Difference between blayer average and parameter at base of blayer 357 REAL(wp), DIMENSION(jpi,jpj) :: zdt_ml,zds_ml,zdu_ml,zdv_ml,zdb_ml ! Difference between mixed layer average and parameter at base of blayer 358 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdbdz_pyc ! parametrised gradient of buoyancy in the pycnocline 359 REAL(wp), DIMENSION(jpi,jpj) :: zdbds_mle ! Magnitude of horizontal buoyancy gradient. 347 360 ! Flux-gradient relationship variables 348 REAL(wp), DIMENSION(jpi, jpj) ::zshear ! Shear production349 350 REAL(wp), DIMENSION(jpi,jpj) :: zhbl_t ! holds boundary layer depth updated by full timestep351 361 REAL(wp), DIMENSION(jpi,jpj) :: zshear ! Shear production 362 ! 363 REAL(wp), DIMENSION(jpi,jpj) :: zhbl_t ! Holds boundary layer depth updated by full timestep 364 ! 352 365 ! For calculating Ri#-dependent mixing 353 366 REAL(wp), DIMENSION(jpi,jpj) :: z2du ! u-shear^2 354 367 REAL(wp), DIMENSION(jpi,jpj) :: z2dv ! v-shear^2 355 368 REAL(wp) :: zrimix ! Spatial form of ri#-induced diffusion 356 369 ! 357 370 ! Temporary variables 358 INTEGER :: inhml 359 REAL(wp) :: znd,znd_d,zznd_ml,zznd_pyc,zznd_d ! temporary non-dimensional depths used in various routines 360 REAL(wp) :: ztemp, zari, zpert, zzdhdt, zdb ! temporary variables 361 REAL(wp) :: zthick, zz0, zz1 ! temporary variables 362 REAL(wp) :: zvel_max, zhbl_s ! temporary variables 363 REAL(wp) :: zfac, ztmp ! temporary variable 364 REAL(wp) :: zus_x, zus_y ! temporary Stokes drift 365 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zviscos ! viscosity 366 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdiffut ! t-diffusivity 367 REAL(wp), DIMENSION(jpi,jpj) :: zalpha_pyc 368 INTEGER :: ibld_ext=0 ! does not have to be zero for modified scheme 369 REAL(wp) :: zgamma_b_nd, zgamma_b, zdhoh, ztau 370 REAL(wp) :: zzeta_s = 0._wp 371 REAL(wp) :: zzeta_v = 0.46 372 REAL(wp) :: zabsstke 373 REAL(wp) :: zsqrtpi, z_two_thirds, zproportion, ztransp, zthickness 374 REAL(wp) :: z2k_times_thickness, zsqrt_depth, zexp_depth, zdstokes0, zf, zexperfc 375 371 REAL(wp) :: znd ! Temporary non-dimensional depth 372 REAL(wp) :: zz0, zz1, zfac 373 REAL(wp) :: zus_x, zus_y ! Temporary Stokes drift 374 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zviscos ! Viscosity 375 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdiffut ! t-diffusivity 376 REAL(wp), DIMENSION(jpi,jpj) :: zalpha_pyc 377 REAL(wp) :: zabsstke 378 REAL(wp) :: zsqrtpi, z_two_thirds, zthickness 379 REAL(wp) :: z2k_times_thickness, zsqrt_depth, zexp_depth, zf, zexperfc 380 ! 376 381 ! For debugging 377 INTEGER :: ikt 378 REAL(wp) :: zlarge = -1.0e10_wp, zero = 0.0_wp 379 !!-------------------------------------------------------------------- 382 REAL(wp), PARAMETER :: pp_large = -1e10_wp 380 383 ! 381 384 IF( ln_timing ) CALL timing_start('zdf_osm') 382 nbld(:,:) = 0 ; nmld(:,:) = 0 383 sustar(:,:) = zlarge 384 swstrl(:,:) = zlarge ; svstr(:,:) = zlarge ; swstrc(:,:) = zlarge ; suw0(:,:) = zlarge 385 swth0(:,:) = zlarge ; sws0(:,:) = zlarge ; swb0(:,:) = zlarge 386 swthav(:,:) = zlarge ; swsav(:,:) = zlarge ; swbav(:,:) = zlarge 387 sustke(:,:) = zlarge ; sla(:,:) = zlarge ; scos_wind(:,:) = zlarge ; ssin_wind(:,:) = zlarge 388 shol(:,:) = zlarge ; zalpha_pyc(:,:) = zlarge 389 l_pyc(:,:) = .FALSE. ; l_flux(:,:) = .FALSE. ; l_mle(:,:) = .FALSE. 390 ! mixed layer 391 ! no initialization of zhbl or zhml (or zdh?) 392 zhbl(:,:) = zlarge ; zhml(:,:) = zlarge ; zdh(:,:) = zlarge 393 zt_bl(:,:) = zlarge ; zs_bl(:,:) = zlarge ; zu_bl(:,:) = zlarge 394 zv_bl(:,:) = zlarge ; zb_bl(:,:) = zlarge 395 zt_ml(:,:) = zlarge ; zs_ml(:,:) = zlarge ; zu_ml(:,:) = zlarge 396 zt_mle(:,:) = zlarge ; zs_mle(:,:) = zlarge ; zu_mle(:,:) = zlarge 397 zb_mle(:,:) = zlarge 398 zv_ml(:,:) = zlarge ; zdt_bl(:,:) = zlarge ; zds_bl(:,:) = zlarge 399 zdu_bl(:,:) = zlarge ; zdv_bl(:,:) = zlarge ; zdb_bl(:,:) = zlarge 400 zdt_ml(:,:) = zlarge ; zds_ml(:,:) = zlarge ; zdu_ml(:,:) = zlarge ; zdv_ml(:,:) = zlarge 401 zdb_ml(:,:) = zlarge 402 ! 403 zdbdz_pyc(:,:,:) = zlarge 385 ! 386 nbld(:,:) = 0 ; nmld(:,:) = 0 387 sustar(:,:) = pp_large 388 swstrl(:,:) = pp_large ; svstr(:,:) = pp_large ; swstrc(:,:) = pp_large ; suw0(:,:) = pp_large 389 swth0(:,:) = pp_large ; sws0(:,:) = pp_large ; swb0(:,:) = pp_large 390 swthav(:,:) = pp_large ; swsav(:,:) = pp_large ; swbav(:,:) = pp_large 391 sustke(:,:) = pp_large ; sla(:,:) = pp_large ; scos_wind(:,:) = pp_large ; ssin_wind(:,:) = pp_large 392 shol(:,:) = pp_large ; zalpha_pyc(:,:) = pp_large 393 l_pyc(:,:) = .FALSE. ; l_flux(:,:) = .FALSE. ; l_mle(:,:) = .FALSE. 394 ! Mixed layer 395 ! No initialization of zhbl or zhml (or zdh?) 396 zhbl(:,:) = pp_large ; zhml(:,:) = pp_large ; zdh(:,:) = pp_large 397 av_t_bl(:,:) = pp_large ; av_s_bl(:,:) = pp_large ; av_u_bl(:,:) = pp_large 398 av_v_bl(:,:) = pp_large ; av_b_bl(:,:) = pp_large ; av_t_ml(:,:) = pp_large 399 av_s_ml(:,:) = pp_large ; av_u_ml(:,:) = pp_large ; av_v_ml(:,:) = pp_large 400 av_b_ml(:,:) = pp_large ; av_t_mle(:,:) = pp_large ; av_s_mle(:,:) = pp_large 401 av_u_mle(:,:) = pp_large ; av_v_mle(:,:) = pp_large ; av_b_mle(:,:) = pp_large 402 zdt_bl(:,:) = pp_large ; zds_bl(:,:) = pp_large ; zdu_bl(:,:) = pp_large 403 zdv_bl(:,:) = pp_large ; zdb_bl(:,:) = pp_large ; zdt_ml(:,:) = pp_large 404 zds_ml(:,:) = pp_large ; zdu_ml(:,:) = pp_large ; zdv_ml(:,:) = pp_large 405 zdb_ml(:,:) = pp_large 406 ! 407 zdbdz_pyc(:,:,:) = pp_large 404 408 zdbdz_pyc(A2D(0),:) = 0.0_wp 405 409 ! 406 407 IF ( ln_osm_mle ) THEN ! only initialise arrays if needed 408 zdtdx(:,:) = zlarge ; zdtdy(:,:) = zlarge ; zdsdx(:,:) = zlarge 409 zdsdy(:,:) = zlarge ; dbdx_mle(:,:) = zlarge ; dbdy_mle(:,:) = zlarge 410 zwb_fk(:,:) = zlarge ; zvel_mle(:,:) = zlarge ; zdiff_mle(:,:) = zlarge 411 zhmle(:,:) = zlarge ; zmld(:,:) = zlarge 410 IF ( ln_osm_mle ) THEN ! Only initialise arrays if needed 411 zdtdx(:,:) = pp_large ; zdtdy(:,:) = pp_large ; zdsdx(:,:) = pp_large 412 zdsdy(:,:) = pp_large ; dbdx_mle(:,:) = pp_large ; dbdy_mle(:,:) = pp_large 413 zwb_fk(:,:) = pp_large ; zvel_mle(:,:) = pp_large ; zdiff_mle(:,:) = pp_large 414 zhmle(:,:) = pp_large ; zmld(:,:) = pp_large 412 415 ENDIF 413 zhbl_t(:,:) = zlarge414 415 zdiffut(:,:,:) = zlarge ; zviscos(:,:,:) = zlarge416 zdiffut(A2D(0),:) = 0.0_wp ; zviscos(A2D(0),:) = 0.0_wp417 ghamt(:,:,:) = zlarge ; ghams(:,:,:) = zlarge418 ghamt(A2D(0),:) = 0.0_wp ; ghams(A2D(0),:) = 0.0_wp419 ghamu(:,:,:) = zlarge ; ghamv(:,:,:) = zlarge420 ghamu(A2D(0),:) = 0.0_wp ; ghamv(A2D(0),:) = 0.0_wp416 zhbl_t(:,:) = pp_large 417 ! 418 zdiffut(:,:,:) = pp_large ; zviscos(:,:,:) = pp_large 419 zdiffut(A2D(0),:) = 0.0_wp ; zviscos(A2D(0),:) = 0.0_wp 420 ghamt(:,:,:) = pp_large ; ghams(:,:,:) = pp_large 421 ghamt(A2D(0),:) = 0.0_wp ; ghams(A2D(0),:) = 0.0_wp 422 ghamu(:,:,:) = pp_large ; ghamv(:,:,:) = pp_large 423 ghamu(A2D(0),:) = 0.0_wp ; ghamv(A2D(0),:) = 0.0_wp 421 424 zdiff_mle(A2D(0)) = 0.0_wp 422 423 425 ! 424 426 #ifdef key_osm_debug 425 427 IF(mi0(nn_idb)==mi1(nn_idb) .AND. mj0(nn_jdb)==mj1(nn_jdb) .AND. & … … 427 429 nn_narea_db = narea 428 430 iloc_db=mi0(nn_idb); jloc_db=mj0(nn_jdb) 429 430 431 WRITE(narea+100,*) 431 432 WRITE(narea+100,'(a,i7)')'timestep=',kt … … 441 442 END IF 442 443 #endif 443 444 ! 444 445 ! hbl = MAX(hbl,epsln) 445 446 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 451 452 zz1 = 1.0_wp - rn_abs 452 453 DO_2D( 0, 0, 0, 0 ) 453 zrad0(ji,jj) = qsr(ji,jj) * r1_rho0_rcp 454 zrad0(ji,jj) = qsr(ji,jj) * r1_rho0_rcp ! Surface downward irradiance (so always +ve) 454 455 zradh(ji,jj) = zrad0(ji,jj) * & ! Downwards irradiance at base of boundary layer 455 456 & ( zz0 * EXP( -1.0_wp * hbl(ji,jj) / rn_si0 ) + zz1 * EXP( -1.0_wp * hbl(ji,jj) / rn_si1 ) ) 456 zradav = zrad0(ji,jj) * & ! Downwards irradiance averaged over depth of the OSBL457 & ( zz0 * ( 1.0_wp - EXP( -hbl(ji,jj)/rn_si0 ) ) * rn_si0 + &458 & zz1 * ( 1.0_wp - EXP( -hbl(ji,jj)/rn_si1 ) ) * rn_si1 ) / hbl(ji,jj)459 swth0(ji,jj) = - qns(ji,jj) * r1_rho0_rcp * tmask(ji,jj,1) 460 swthav(ji,jj) = 0.5_wp * swth0(ji,jj) - & ! Turbulent heat flux averaged over depth of OSBL461 & ( 0.5_wp * ( zrad0(ji,jj) + zradh(ji,jj) ) - zradav )457 zradav = zrad0(ji,jj) * & ! Downwards irradiance averaged 458 & ( zz0 * ( 1.0_wp - EXP( -hbl(ji,jj)/rn_si0 ) ) * rn_si0 + & ! over depth of the OSBL 459 & zz1 * ( 1.0_wp - EXP( -hbl(ji,jj)/rn_si1 ) ) * rn_si1 ) / hbl(ji,jj) 460 swth0(ji,jj) = - qns(ji,jj) * r1_rho0_rcp * tmask(ji,jj,1) ! Upwards surface Temperature flux for non-local term 461 swthav(ji,jj) = 0.5_wp * swth0(ji,jj) - ( 0.5_wp * ( zrad0(ji,jj) + zradh(ji,jj) ) - & ! Turbulent heat flux averaged 462 & zradav ) ! over depth of OSBL 462 463 END_2D 463 464 DO_2D( 0, 0, 0, 0 ) 464 sws0(ji,jj) = -1.0_wp * & ! Upwards surface salinity flux for non-local term465 & ( ( emp(ji,jj) - rnf(ji,jj) ) * ts(ji,jj,1,jp_sal,Kmm) + sfx(ji,jj) ) * r1_rho0 * tmask(ji,jj,1)465 sws0(ji,jj) = -1.0_wp * ( ( emp(ji,jj) - rnf(ji,jj) ) * ts(ji,jj,1,jp_sal,Kmm) + & ! Upwards surface salinity flux 466 & sfx(ji,jj) ) * r1_rho0 * tmask(ji,jj,1) ! for non-local term 466 467 zthermal = rab_n(ji,jj,1,jp_tem) 467 468 zbeta = rab_n(ji,jj,1,jp_sal) 468 swb0(ji,jj) = grav * zthermal * swth0(ji,jj) - & ! Non radiative upwards surface buoyancy flux 469 & grav * zbeta * sws0(ji,jj) 470 zwb0tot(ji,jj) = swb0(ji,jj) - grav * zthermal * & ! Total upwards surface buoyancy flux 471 & ( zrad0(ji,jj) - zradh(ji,jj) ) 472 swsav(ji,jj) = 0.5 * sws0(ji,jj) ! Turbulent salinity flux averaged over depth of the OBSL 469 swb0(ji,jj) = grav * zthermal * swth0(ji,jj) - grav * zbeta * sws0(ji,jj) ! Non radiative upwards surface buoyancy flux 470 zwb0tot(ji,jj) = swb0(ji,jj) - grav * zthermal * ( zrad0(ji,jj) - zradh(ji,jj) ) ! Total upwards surface buoyancy flux 471 swsav(ji,jj) = 0.5_wp * sws0(ji,jj) ! Turbulent salinity flux averaged over depth of the OBSL 473 472 swbav(ji,jj) = grav * zthermal * swthav(ji,jj) - & ! Turbulent buoyancy flux averaged over the depth of the 474 473 & grav * zbeta * swsav(ji,jj) ! OBSBL 475 474 END_2D 476 475 DO_2D( 0, 0, 0, 0 ) 477 suw0(ji,jj) = - 0.5 * (utau(ji-1,jj) + utau(ji,jj)) * & ! Surface upward velocity fluxes 478 & r1_rho0 * tmask(ji,jj,1) 479 zvw0 = - 0.5 * (vtau(ji,jj-1) + vtau(ji,jj)) * r1_rho0 * tmask(ji,jj,1) 480 sustar(ji,jj) = MAX( SQRT( SQRT( suw0(ji,jj) * & ! Friction velocity (sustar), at T-point : LMD94 eq. 2 481 & suw0(ji,jj) + zvw0 * zvw0 ) ), 1.0e-8_wp ) 482 scos_wind(ji,jj) = -suw0(ji,jj) / ( sustar(ji,jj) * sustar(ji,jj) ) 483 ssin_wind(ji,jj) = -zvw0 / ( sustar(ji,jj) * sustar(ji,jj) ) 476 suw0(ji,jj) = -0.5_wp * (utau(ji-1,jj) + utau(ji,jj)) * r1_rho0 * tmask(ji,jj,1) ! Surface upward velocity fluxes 477 zvw0 = -0.5_wp * (vtau(ji,jj-1) + vtau(ji,jj)) * r1_rho0 * tmask(ji,jj,1) 478 sustar(ji,jj) = MAX( SQRT( SQRT( suw0(ji,jj) * suw0(ji,jj) + zvw0 * zvw0 ) ), & ! Friction velocity (sustar), at 479 & 1e-8_wp ) ! T-point : LMD94 eq. 2 480 scos_wind(ji,jj) = -1.0_wp * suw0(ji,jj) / ( sustar(ji,jj) * sustar(ji,jj) ) 481 ssin_wind(ji,jj) = -1.0_wp * zvw0 / ( sustar(ji,jj) * sustar(ji,jj) ) 484 482 #ifdef key_osm_debug 485 483 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN … … 503 501 CASE(0) 504 502 DO_2D( 0, 0, 0, 0 ) 505 zus_x = scos_wind(ji,jj) * sustar(ji,jj) / 0.3 **2506 zus_y = ssin_wind(ji,jj) * sustar(ji,jj) / 0.3 **2503 zus_x = scos_wind(ji,jj) * sustar(ji,jj) / 0.3_wp**2 504 zus_y = ssin_wind(ji,jj) * sustar(ji,jj) / 0.3_wp**2 507 505 ! Linearly 508 sustke(ji,jj) = MAX ( SQRT( zus_x*zus_x + zus_y*zus_y), 1.0e-8)506 sustke(ji,jj) = MAX( SQRT( zus_x * zus_x + zus_y * zus_y ), 1e-8_wp ) 509 507 dstokes(ji,jj) = rn_osm_dstokes 510 508 END_2D … … 513 511 DO_2D( 0, 0, 0, 0 ) 514 512 ! Use wind speed wndm included in sbc_oce module 515 sustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8)516 dstokes(ji,jj) = MAX ( 0.12 * wndm(ji,jj)**2 / grav, 5.e-1)513 sustke(ji,jj) = MAX ( 0.016_wp * wndm(ji,jj), 1e-8_wp ) 514 dstokes(ji,jj) = MAX ( 0.12_wp * wndm(ji,jj)**2 / grav, 5e-1_wp ) 517 515 END_2D 518 516 ! Use ECMWF wave fields as output from SBCWAVE 519 517 CASE(2) 520 518 zfac = 2.0_wp * rpi / 16.0_wp 521 519 ! 522 520 DO_2D( 0, 0, 0, 0 ) 523 IF ( hsw(ji,jj) > 1.e-4) THEN521 IF ( hsw(ji,jj) > 1e-4_wp ) THEN 524 522 ! Use wave fields 525 zabsstke = SQRT( ut0sd(ji,jj)**2 + vt0sd(ji,jj)**2)526 sustke(ji,jj) = MAX ( ( scos_wind(ji,jj) * ut0sd(ji,jj) + ssin_wind(ji,jj) * vt0sd(ji,jj) ), 1.0e-8)527 dstokes(ji,jj) = MAX (zfac * hsw(ji,jj)*hsw(ji,jj) / ( MAX(zabsstke * wmp(ji,jj), 1.0e-7 ) ), 5.0e-1)523 zabsstke = SQRT( ut0sd(ji,jj)**2 + vt0sd(ji,jj)**2 ) 524 sustke(ji,jj) = MAX( ( scos_wind(ji,jj) * ut0sd(ji,jj) + ssin_wind(ji,jj) * vt0sd(ji,jj) ), 1e-8_wp ) 525 dstokes(ji,jj) = MAX( zfac * hsw(ji,jj) * hsw(ji,jj) / ( MAX( zabsstke * wmp(ji,jj), 1e-7 ) ), 5e-1_wp ) 528 526 ELSE 529 527 ! Assume masking issue (e.g. ice in ECMWF reanalysis but not in model run) 530 528 ! .. so default to Pierson-Moskowitz 531 sustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8)532 dstokes(ji,jj) = MAX ( 0.12 * wndm(ji,jj)**2 / grav, 5.e-1)529 sustke(ji,jj) = MAX( 0.016_wp * wndm(ji,jj), 1e-8_wp ) 530 dstokes(ji,jj) = MAX( 0.12_wp * wndm(ji,jj)**2 / grav, 5e-1_wp ) 533 531 END IF 534 532 END_2D … … 541 539 END IF 542 540 #endif 543 541 ! 544 542 IF (ln_zdfosm_ice_shelter) THEN 545 543 ! Reduce both Stokes drift and its depth scale by ocean fraction to represent sheltering by ice 546 544 DO_2D( 0, 0, 0, 0 ) 547 sustke(ji,jj) = sustke(ji,jj) * (1.0_wp - fr_i(ji,jj))548 dstokes(ji,jj) = dstokes(ji,jj) * ( 1.0_wp - fr_i(ji,jj))545 sustke(ji,jj) = sustke(ji,jj) * ( 1.0_wp - fr_i(ji,jj) ) 546 dstokes(ji,jj) = dstokes(ji,jj) * ( 1.0_wp - fr_i(ji,jj) ) 549 547 END_2D 550 548 END IF 551 549 ! 552 550 SELECT CASE (nn_osm_SD_reduce) 553 ! Reduce surface Stokes drift by a constant factor or following Breivik (2016) + van 551 ! Reduce surface Stokes drift by a constant factor or following Breivik (2016) + van Roekel (2012) or Grant (2020). 554 552 CASE(0) 555 ! The Langmur number from the ECMWF model (or from PM) appears to give La<0.3 for wind-driven seas. 556 ! The coefficient rn_zdfosm_adjust_sd = 0.8 gives La=0.3 in this situation. 557 ! It could represent the effects of the spread of wave directions 558 ! around the mean wind. The effect of this adjustment needs to be tested. 553 ! The Langmur number from the ECMWF model (or from PM) appears to give La<0.3 for wind-driven seas. 554 ! The coefficient rn_zdfosm_adjust_sd = 0.8 gives La=0.3 in this situation. 555 ! It could represent the effects of the spread of wave directions around the mean wind. The effect of this adjustment needs to be tested. 559 556 IF(nn_osm_wave > 0) THEN 560 557 sustke(A2D(0)) = rn_zdfosm_adjust_sd * sustke(A2D(0)) 561 558 END IF 562 559 CASE(1) 563 ! vanRoekel (2012): consider average SD over top 10% of boundary layer564 ! assumes approximate depth profile of SD from Breivik (2016)560 ! Van Roekel (2012): consider average SD over top 10% of boundary layer 561 ! Assumes approximate depth profile of SD from Breivik (2016) 565 562 zsqrtpi = SQRT(rpi) 566 563 z_two_thirds = 2.0_wp / 3.0_wp 567 568 564 DO_2D( 0, 0, 0, 0 ) 569 565 zthickness = rn_osm_hblfrac*hbl(ji,jj) 570 z2k_times_thickness = zthickness * 2.0_wp / MAX( ABS( 5.97_wp * dstokes(ji,jj) ), 0.0000001_wp )571 zsqrt_depth = SQRT( z2k_times_thickness)572 zexp_depth = EXP( -z2k_times_thickness)573 sustke(ji,jj) = sustke(ji,jj) * ( 1.0_wp - zexp_depth&574 & - z_two_thirds * ( zsqrtpi*zsqrt_depth*z2k_times_thickness * ERFC(zsqrt_depth)&575 & + 1.0_wp - (1.0_wp + z2k_times_thickness)*zexp_depth ) ) / z2k_times_thickness576 566 z2k_times_thickness = zthickness * 2.0_wp / MAX( ABS( 5.97_wp * dstokes(ji,jj) ), 1e-7_wp ) 567 zsqrt_depth = SQRT( z2k_times_thickness ) 568 zexp_depth = EXP( -1.0_wp * z2k_times_thickness ) 569 sustke(ji,jj) = sustke(ji,jj) * ( 1.0_wp - zexp_depth - & 570 & z_two_thirds * ( zsqrtpi * zsqrt_depth * z2k_times_thickness * ERFC(zsqrt_depth) + & 571 & 1.0_wp - ( 1.0_wp + z2k_times_thickness ) * zexp_depth ) ) / & 572 & z2k_times_thickness 577 573 END_2D 578 574 CASE(2) 579 575 ! Grant (2020): Match to exponential with same SD and d/dz(Sd) at depth 10% of boundary layer 580 ! assumes approximate depth profile of SD from Breivik (2016)576 ! Assumes approximate depth profile of SD from Breivik (2016) 581 577 zsqrtpi = SQRT(rpi) 582 583 578 DO_2D( 0, 0, 0, 0 ) 584 579 zthickness = rn_osm_hblfrac*hbl(ji,jj) 585 z2k_times_thickness = zthickness * 2.0_wp / MAX( ABS( 5.97_wp * dstokes(ji,jj) ), 0.0000001_wp ) 586 587 IF(z2k_times_thickness < 50._wp) THEN 588 zsqrt_depth = SQRT(z2k_times_thickness) 589 zexperfc = zsqrtpi * zsqrt_depth * ERFC(zsqrt_depth) * EXP(z2k_times_thickness) 580 z2k_times_thickness = zthickness * 2.0_wp / MAX( ABS( 5.97_wp * dstokes(ji,jj) ), 1e-7_wp ) 581 IF( z2k_times_thickness < 50.0_wp ) THEN 582 zsqrt_depth = SQRT( z2k_times_thickness ) 583 zexperfc = zsqrtpi * zsqrt_depth * ERFC(zsqrt_depth) * EXP( z2k_times_thickness ) 590 584 ELSE 591 ! asymptotic expansion of sqrt(pi)*zsqrt_depth*EXP(z2k_times_thickness)*ERFC(zsqrt_depth) for large z2k_times_thickness585 ! Asymptotic expansion of sqrt(pi)*zsqrt_depth*EXP(z2k_times_thickness)*ERFC(zsqrt_depth) for large z2k_times_thickness 592 586 ! See Abramowitz and Stegun, Eq. 7.1.23 593 ! zexperfc = 1._wp - (1/2)/(z2k_times_thickness) + (3/4)/(z2k_times_thickness**2) - (15/8)/(z2k_times_thickness**3) 594 zexperfc = ((- 1.875_wp/z2k_times_thickness + 0.75_wp)/z2k_times_thickness - 0.5_wp)/z2k_times_thickness + 1.0_wp 587 ! zexperfc = 1._wp - (1/2)/(z2k_times_thickness) + (3/4)/(z2k_times_thickness**2) - (15/8)/(z2k_times_thickness**3) 588 zexperfc = ( ( -1.875_wp / z2k_times_thickness + 0.75_wp ) / z2k_times_thickness - 0.5_wp ) / & 589 & z2k_times_thickness + 1.0_wp 595 590 END IF 596 zf = z2k_times_thickness*(1.0_wp/zexperfc - 1.0_wp) 597 dstokes(ji,jj) = 5.97 * zf * dstokes(ji,jj) 598 sustke(ji,jj) = sustke(ji,jj) * EXP(z2k_times_thickness * ( 1.0_wp / (2. * zf) - 1.0_wp )) * ( 1.0_wp - zexperfc) 591 zf = z2k_times_thickness * ( 1.0_wp / zexperfc - 1.0_wp ) 592 dstokes(ji,jj) = 5.97_wp * zf * dstokes(ji,jj) 593 sustke(ji,jj) = sustke(ji,jj) * EXP( z2k_times_thickness * ( 1.0_wp / ( 2.0_wp * zf ) - 1.0_wp ) ) * & 594 & ( 1.0_wp - zexperfc ) 599 595 END_2D 600 596 END SELECT 601 597 ! 602 598 ! Langmuir velocity scale (swstrl), La # (sla) 603 ! mixed scale (svstr), convective velocity scale (swstrc)599 ! Mixed scale (svstr), convective velocity scale (swstrc) 604 600 DO_2D( 0, 0, 0, 0 ) 605 601 ! Langmuir velocity scale (swstrl), at T-point 606 602 swstrl(ji,jj) = ( sustar(ji,jj) * sustar(ji,jj) * sustke(ji,jj) )**pthird 607 sla(ji,jj) = MAX(MIN(SQRT ( sustar(ji,jj) / ( swstrl(ji,jj) + epsln ) )**3, 4.0), 0.2)608 IF (sla(ji,jj) > 0.45) dstokes(ji,jj) = MIN(dstokes(ji,jj), 0.5_wp*hbl(ji,jj))603 sla(ji,jj) = MAX( MIN( SQRT( sustar(ji,jj) / ( swstrl(ji,jj) + epsln ) )**3, 4.0_wp ), 0.2_wp ) 604 IF ( sla(ji,jj) > 0.45_wp ) dstokes(ji,jj) = MIN( dstokes(ji,jj), 0.5_wp * hbl(ji,jj) ) 609 605 ! Velocity scale that tends to sustar for large Langmuir numbers 610 svstr(ji,jj) = ( swstrl(ji,jj)**3 + & 611 & ( 1.0 - EXP( -0.5 * sla(ji,jj)**2 ) ) * sustar(ji,jj) * sustar(ji,jj) * sustar(ji,jj) )**pthird 612 613 ! limit maximum value of Langmuir number as approximate treatment for shear turbulence. 614 ! Note sustke and swstrl are not amended. 606 svstr(ji,jj) = ( swstrl(ji,jj)**3 + ( 1.0_wp - EXP( -0.5_wp * sla(ji,jj)**2 ) ) * sustar(ji,jj) * sustar(ji,jj) * & 607 & sustar(ji,jj) )**pthird 615 608 ! 616 ! get convective velocity (swstrc), stabilty scale (shol) and logical conection flag l_conv 617 IF ( swbav(ji,jj) > 0.0) THEN 618 swstrc(ji,jj) = ( 2.0 * swbav(ji,jj) * 0.9 * hbl(ji,jj) )**pthird 619 shol(ji,jj) = -0.9 * hbl(ji,jj) * 2.0 * swbav(ji,jj) / (svstr(ji,jj)**3 + epsln ) 609 ! Limit maximum value of Langmuir number as approximate treatment for shear turbulence 610 ! Note sustke and swstrl are not amended 611 ! 612 ! Get convective velocity (swstrc), stabilty scale (shol) and logical conection flag l_conv 613 IF ( swbav(ji,jj) > 0.0_wp ) THEN 614 swstrc(ji,jj) = ( 2.0_wp * swbav(ji,jj) * 0.9_wp * hbl(ji,jj) )**pthird 615 shol(ji,jj) = -0.9_wp * hbl(ji,jj) * 2.0_wp * swbav(ji,jj) / ( svstr(ji,jj)**3 + epsln ) 620 616 ELSE 621 617 swstrc(ji,jj) = 0.0_wp 622 shol(ji,jj) = -hbl(ji,jj) * 2.0 * swbav(ji,jj)/ (svstr(ji,jj)**3 + epsln )618 shol(ji,jj) = -1.0_wp * hbl(ji,jj) * 2.0_wp * swbav(ji,jj) / ( svstr(ji,jj)**3 + epsln ) 623 619 ENDIF 624 620 #ifdef key_osm_debug … … 632 628 #endif 633 629 END_2D 634 630 ! 635 631 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 636 632 ! Mixed-layer model - calculate averages over the boundary layer, and the change in the boundary layer depth … … 640 636 ! zdf_osm_zmld_horizontal_gradients need halo values for nbld, so must 641 637 ! previously exist for hbl also. 642 638 ! 643 639 ! agn 23/6/20: not clear all this is needed, as hbl checked after it is re-calculated anyway 644 640 ! ########################################################################## … … 651 647 END_3D 652 648 ! ########################################################################## 653 649 ! 654 650 DO_2D( 0, 0, 0, 0 ) 655 651 zhbl(ji,jj) = gdepw(ji,jj,nbld(ji,jj),Kmm) … … 665 661 &' zhbl =',zhbl(ji,jj) , ' zhml=', zhml(ji,jj), ' zdh=', zdh(ji,jj),& 666 662 &' imld=', nmld(ji,jj), ' ibld=', nbld(ji,jj) 667 668 663 WRITE(narea+100,'(a,g11.3,a,2g11.3)') 'Physics: ssh ',ssh(ji,jj,Kmm),' T S surface=',ts(ji,jj,1,jp_tem,Kmm),ts(ji,jj,1,jp_sal,Kmm) 669 664 jl = nmld(ji,jj) - 1; jm = MIN( nbld(ji,jj) + 2, mbkt(ji,jj) ) … … 679 674 END IF 680 675 #endif 681 676 ! 682 677 ! Averages over well-mixed and boundary layer, note BL averages use jp_ext=2 everywhere 683 678 jp_ext(:,:) = 1 ! ag 19/03 684 CALL zdf_osm_vertical_average( Kbb, Kmm, nbld, zt_bl, zs_bl, &685 & zb_bl, zu_bl, zv_bl, jp_ext, zdt_bl, &679 CALL zdf_osm_vertical_average( Kbb, Kmm, nbld, av_t_bl, av_s_bl, & 680 & av_b_bl, av_u_bl, av_v_bl, jp_ext, zdt_bl, & 686 681 & zds_bl, zdb_bl, zdu_bl, zdv_bl ) 687 682 jp_ext(:,:) = nbld(:,:) - nmld(:,:) + jp_ext(:,:) + 1 ! ag 19/03 688 CALL zdf_osm_vertical_average( Kbb, Kmm, nmld - 1, zt_ml, zs_ml, &689 & zb_ml, zu_ml, zv_ml,jp_ext, zdt_ml, &690 & zds_ml, zdb_ml, zdu_ml, 683 CALL zdf_osm_vertical_average( Kbb, Kmm, nmld - 1, av_t_ml, av_s_ml, & 684 & av_b_ml, av_u_ml, av_v_ml, jp_ext, zdt_ml, & 685 & zds_ml, zdb_ml, zdu_ml, zdv_ml ) 691 686 #ifdef key_osm_debug 692 687 IF(narea==nn_narea_db) THEN 693 688 ji=iloc_db; jj=jloc_db 694 689 WRITE(narea+100,'(4(3(a,g11.3),/), 2(4(a,g11.3),/))') & 695 & 'After averaging, with old hbl (& jp_ext==2), hml: zt_bl=', zt_bl(ji,jj),&696 & ' zs_bl=', zs_bl(ji,jj), ' zb_bl=', zb_bl(ji,jj),&690 & 'After averaging, with old hbl (& jp_ext==2), hml: zt_bl=', av_t_bl(ji,jj),& 691 & ' zs_bl=', av_s_bl(ji,jj), ' zb_bl=', av_b_bl(ji,jj),& 697 692 & 'zdt_bl=', zdt_bl(ji,jj), ' zds_bl=', zds_bl(ji,jj), ' zdb_bl=', zdb_bl(ji,jj),& 698 & 'zt_ml=', zt_ml(ji,jj), ' zs_ml=', zs_ml(ji,jj), ' zb_ml=', zb_ml(ji,jj),&693 & 'zt_ml=', av_t_ml(ji,jj), ' zs_ml=', av_s_ml(ji,jj), ' zb_ml=', av_b_ml(ji,jj),& 699 694 & 'zdt_ml=', zdt_ml(ji,jj), ' zds_ml=', zds_ml(ji,jj), ' zdb_ml=', zdb_ml(ji,jj),& 700 & 'zu_bl =', zu_bl(ji,jj) , ' zv_bl=', zv_bl(ji,jj), ' zdu_bl=', zdu_bl(ji,jj), ' zdv_bl=', zdv_bl(ji,jj),&701 & 'zu_ml =', zu_ml(ji,jj) , ' zv_ml=', zv_ml(ji,jj), ' zdu_ml=', zdu_ml(ji,jj), ' zdv_ml=', zdv_ml(ji,jj)695 & 'zu_bl =', av_u_bl(ji,jj) , ' zv_bl=', av_v_bl(ji,jj), ' zdu_bl=', zdu_bl(ji,jj), ' zdv_bl=', zdv_bl(ji,jj),& 696 & 'zu_ml =', av_u_ml(ji,jj) , ' zv_ml=', av_v_ml(ji,jj), ' zdu_ml=', zdu_ml(ji,jj), ' zdv_ml=', zdv_ml(ji,jj) 702 697 FLUSH(narea+100) 703 698 END IF 704 699 #endif 705 700 ! Velocity components in frame aligned with surface stress 706 CALL zdf_osm_velocity_rotation( zu_ml, zv_ml )701 CALL zdf_osm_velocity_rotation( av_u_ml, av_v_ml ) 707 702 CALL zdf_osm_velocity_rotation( zdu_ml, zdv_ml ) 708 CALL zdf_osm_velocity_rotation( zu_bl, zv_bl )703 CALL zdf_osm_velocity_rotation( av_u_bl, av_v_bl ) 709 704 CALL zdf_osm_velocity_rotation( zdu_bl, zdv_bl ) 710 705 #ifdef key_osm_debug … … 713 708 WRITE(narea+100,'(a,/, 2(4(a,g11.3),/))') & 714 709 & 'After rotation, with old hbl (& jp_ext==2), hml:', & 715 & 'zu_bl =', zu_bl(ji,jj) , ' zv_bl=', zv_bl(ji,jj), ' zdu_bl=', zdu_bl(ji,jj), ' zdv_bl=', zdv_bl(ji,jj),&716 & 'zu_ml =', zu_ml(ji,jj) , ' zv_ml=', zv_ml(ji,jj), ' zdu_ml=', zdu_ml(ji,jj), ' zdv_ml=', zdv_ml(ji,jj)710 & 'zu_bl =', av_u_bl(ji,jj) , ' zv_bl=', av_v_bl(ji,jj), ' zdu_bl=', zdu_bl(ji,jj), ' zdv_bl=', zdv_bl(ji,jj),& 711 & 'zu_ml =', av_u_ml(ji,jj) , ' zv_ml=', av_v_ml(ji,jj), ' zdu_ml=', zdu_ml(ji,jj), ' zdv_ml=', zdv_ml(ji,jj) 717 712 FLUSH(narea+100) 718 713 END IF 719 714 #endif 720 715 ! 721 716 ! Determine the state of the OSBL, stable/unstable, shear/no shear 722 717 CALL zdf_osm_osbl_state( Kmm, zwb_ent, zwb_min, zshear, zhbl, & 723 718 & zhml, zdh, zdb_bl, zdu_bl, zdv_bl, & 724 719 & zdb_ml, zdu_ml, zdv_ml ) 725 720 ! 726 721 #ifdef key_osm_debug 727 722 IF(narea==nn_narea_db) THEN … … 739 734 IF ( hmle(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk) 740 735 END_3D 741 CALL zdf_osm_vertical_average( Kbb, Kmm, mld_prof, zt_mle, zs_mle, &742 & zb_mle, zu_mle, zv_mle )743 736 CALL zdf_osm_vertical_average( Kbb, Kmm, mld_prof, av_t_mle, av_s_mle, & 737 & av_b_mle, av_u_mle, av_v_mle ) 738 ! 744 739 DO_2D( 0, 0, 0, 0 ) 745 740 zhmle(ji,jj) = gdepw(ji,jj,mld_prof(ji,jj),Kmm) … … 750 745 WRITE(narea+100,'(2(a,g11.3), a, i7,/,(3(a,g11.3),/),2(a,g11.3),/)') & 751 746 & 'Before updating hmle: hmle =',hmle(ji,jj) , ' zhmle=', zhmle(ji,jj), ' mld_prof=', mld_prof(ji,jj), & 752 & 'averaging over hmle: zt_mle=', zt_mle(ji,jj), ' zs_mle=', zs_mle(ji,jj), ' zb_mle=', zb_mle(ji,jj),&753 & 'zu_mle =', zu_mle(ji,jj), ' zv_mle=', zv_mle(ji,jj)747 & 'averaging over hmle: zt_mle=', av_t_mle(ji,jj), ' zs_mle=', av_s_mle(ji,jj), ' zb_mle=', av_b_mle(ji,jj),& 748 & 'zu_mle =', av_u_mle(ji,jj), ' zv_mle=', av_v_mle(ji,jj) 754 749 FLUSH(narea+100) 755 750 END IF 756 751 #endif 757 758 ! !Calculate fairly-well-mixed depth zmld & its index mld_prof + lateral zmld-averaged gradients759 CALL zdf_osm_zmld_horizontal_gradients( Kmm, zmld, zdtdx,zdtdy, zdsdx, &752 ! 753 ! Calculate fairly-well-mixed depth zmld & its index mld_prof + lateral zmld-averaged gradients 754 CALL zdf_osm_zmld_horizontal_gradients( Kmm, zmld, zdtdx, zdtdy, zdsdx, & 760 755 & zdsdy, dbdx_mle, dbdy_mle, zdbds_mle ) 761 !! Calculate max vertical FK flux zwb_fk & set logical descriptors 762 CALL zdf_osm_osbl_state_fk( Kmm, zwb_fk, zt_mle, zs_mle, zb_mle, & 763 & zhbl, zhmle, zwb_ent, zt_bl, zs_bl, & 764 & zb_bl, zdb_bl, zdbds_mle ) 765 !! recalculate hmle, zmle, zvel_mle, zdiff_mle & redefine mld_proc to be index for new hmle 766 CALL zdf_osm_mle_parameters( Kmm, mld_prof, zmld, zhmle, zvel_mle, & 767 & zdiff_mle, zdbds_mle, zhbl, zb_bl, zwb0tot ) 756 ! Calculate max vertical FK flux zwb_fk & set logical descriptors 757 CALL zdf_osm_osbl_state_fk( Kmm, zwb_fk, zhbl, zhmle, zwb_ent, & 758 & zdb_bl, zdbds_mle ) 759 ! Recalculate hmle, zmle, zvel_mle, zdiff_mle & redefine mld_proc to be index for new hmle 760 CALL zdf_osm_mle_parameters( Kmm, mld_prof, zmld, zhmle, zvel_mle, & 761 & zdiff_mle, zdbds_mle, zhbl, zwb0tot ) 768 762 #ifdef key_osm_debug 769 763 IF(narea==nn_narea_db) THEN … … 782 776 ELSE ! ln_osm_mle 783 777 ! FK not selected, Boundary Layer only. 784 l_pyc(:,:) = .TRUE.778 l_pyc(:,:) = .TRUE. 785 779 l_flux(:,:) = .FALSE. 786 l_mle(:,:) = .FALSE.780 l_mle(:,:) = .FALSE. 787 781 DO_2D( 0, 0, 0, 0 ) 788 782 IF ( l_conv(ji,jj) .AND. zdb_bl(ji,jj) < rn_osm_bl_thresh ) l_pyc(ji,jj) = .FALSE. 789 783 END_2D 790 784 ENDIF ! ln_osm_mle 791 785 ! 792 786 !! External gradient below BL needed both with and w/o FK 793 787 CALL zdf_osm_external_gradients( Kmm, nbld + 1, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext ) ! ag 19/03 794 788 ! 795 789 ! Test if pycnocline well resolved 796 ! DO_2D( 0, 0, 0, 0 ) Removed with ag 19/03 changes. A change in eddy diffusivity/viscosity797 ! IF (l_conv(ji,jj) ) THEN should account for this.798 ! ztmp = 0.2 * zhbl(ji,jj) / e3w(ji,jj,nbld(ji,jj),Kmm)799 ! IF ( ztmp > 6 ) THEN800 ! ! pycnocline well resolved801 ! jp_ext(ji,jj) = 1802 ! ELSE803 ! ! pycnocline poorly resolved804 ! jp_ext(ji,jj) = 0805 ! ENDIF806 ! ELSE807 ! ! Stable conditions808 ! jp_ext(ji,jj) = 0809 ! ENDIF810 ! END_2D790 ! DO_2D( 0, 0, 0, 0 ) Removed with ag 19/03 changes. A change in eddy diffusivity/viscosity 791 ! IF (l_conv(ji,jj) ) THEN should account for this. 792 ! ztmp = 0.2 * zhbl(ji,jj) / e3w(ji,jj,nbld(ji,jj),Kmm) 793 ! IF ( ztmp > 6 ) THEN 794 ! ! pycnocline well resolved 795 ! jp_ext(ji,jj) = 1 796 ! ELSE 797 ! ! pycnocline poorly resolved 798 ! jp_ext(ji,jj) = 0 799 ! ENDIF 800 ! ELSE 801 ! ! Stable conditions 802 ! jp_ext(ji,jj) = 0 803 ! ENDIF 804 ! END_2D 811 805 #ifdef key_osm_debug 812 806 IF(narea==nn_narea_db) THEN … … 819 813 END IF 820 814 #endif 821 815 ! 822 816 ! Recalculate bl averages using jp_ext & ml averages .... note no rotation of u & v here.. 823 817 jp_ext(:,:) = 1 ! ag 19/03 824 CALL zdf_osm_vertical_average( Kbb, Kmm, nbld, zt_bl, zs_bl, &825 & zb_bl, zu_bl, zv_bl,jp_ext, zdt_bl, &818 CALL zdf_osm_vertical_average( Kbb, Kmm, nbld, av_t_bl, av_s_bl, & 819 & av_b_bl, av_u_bl, av_v_bl, jp_ext, zdt_bl, & 826 820 & zds_bl, zdb_bl, zdu_bl, zdv_bl ) 827 821 jp_ext(:,:) = nbld(:,:) - nmld(:,:) + jp_ext(:,:) + 1 ! ag 19/03 828 CALL zdf_osm_vertical_average( Kbb, Kmm, nmld - 1, zt_ml, zs_ml, &829 & zb_ml, zu_ml, zv_ml,jp_ext, zdt_ml, &830 & zds_ml, zdb_ml, zdu_ml, 822 CALL zdf_osm_vertical_average( Kbb, Kmm, nmld - 1, av_t_ml, av_s_ml, & 823 & av_b_ml, av_u_ml, av_v_ml, jp_ext, zdt_ml, & 824 & zds_ml, zdb_ml, zdu_ml, zdv_ml ) ! ag 19/03 831 825 #ifdef key_osm_debug 832 826 IF(narea==nn_narea_db) THEN 833 827 ji=iloc_db; jj=jloc_db 834 828 WRITE(narea+100,'(4(3(a,g11.3),/), 2(4(a,g11.3),/))') & 835 & 'After averaging, with old hbl (&correct jp_ext), hml: zt_bl=', zt_bl(ji,jj),&836 & ' zs_bl=', zs_bl(ji,jj), ' zb_bl=', zb_bl(ji,jj),&829 & 'After averaging, with old hbl (&correct jp_ext), hml: zt_bl=', av_t_bl(ji,jj),& 830 & ' zs_bl=', av_s_bl(ji,jj), ' zb_bl=', av_b_bl(ji,jj),& 837 831 & 'zdt_bl=', zdt_bl(ji,jj), ' zds_bl=', zds_bl(ji,jj), ' zdb_bl=', zdb_bl(ji,jj),& 838 & 'zt_ml=', zt_ml(ji,jj), ' zs_ml=', zs_ml(ji,jj), ' zb_ml=', zb_ml(ji,jj),&832 & 'zt_ml=', av_t_ml(ji,jj), ' zs_ml=', av_s_ml(ji,jj), ' zb_ml=', av_b_ml(ji,jj),& 839 833 & 'zdt_ml=', zdt_ml(ji,jj), ' zds_ml=', zds_ml(ji,jj), ' zdb_ml=', zdb_ml(ji,jj),& 840 & 'zu_bl =', zu_bl(ji,jj) , ' zv_bl=', zv_bl(ji,jj), ' zdu_bl=', zdu_bl(ji,jj), ' zdv_bl=', zdv_bl(ji,jj),&841 & 'zu_ml =', zu_ml(ji,jj) , ' zv_ml=', zv_ml(ji,jj), ' zdu_ml=', zdu_ml(ji,jj), ' zdv_ml=', zdv_ml(ji,jj)834 & 'zu_bl =', av_u_bl(ji,jj) , ' zv_bl=', av_v_bl(ji,jj), ' zdu_bl=', zdu_bl(ji,jj), ' zdv_bl=', zdv_bl(ji,jj),& 835 & 'zu_ml =', av_u_ml(ji,jj) , ' zv_ml=', av_v_ml(ji,jj), ' zdu_ml=', zdu_ml(ji,jj), ' zdv_ml=', zdv_ml(ji,jj) 842 836 FLUSH(narea+100) 843 837 END IF 844 838 #endif 845 846 847 ! Rate of change of hbl 848 CALL zdf_osm_calculate_dhdt( zdhdt, zhbl, zdh, zwb_ent, zwb_min, & 849 & zdb_bl, zdbdz_bl_ext, zdb_ml, zwb_fk_b, zwb_fk, & 839 ! 840 ! Rate of change of hbl 841 CALL zdf_osm_calculate_dhdt( zdhdt, zhbl, zdh, zwb_ent, zwb_min, & 842 & zdb_bl, zdbdz_bl_ext, zdb_ml, zwb_fk_b, zwb_fk, & 850 843 & zvel_mle ) 851 844 ! Test if surface boundary layer coupled to bottom … … 853 846 DO_2D( 0, 0, 0, 0 ) 854 847 zhbl_t(ji,jj) = hbl(ji,jj) + ( zdhdt(ji,jj) - ww(ji,jj,nbld(ji,jj)) ) * rn_Dt ! Certainly need ww here, so subtract it 855 ! adjustment to represent limiting by ocean bottom856 IF ( mbkt(ji,jj) > 2 ) THEN ! to ensure mbkt(ji,jj) - 2 > 0 so no incorrect array access848 ! Adjustment to represent limiting by ocean bottom 849 IF ( mbkt(ji,jj) > 2 ) THEN ! To ensure mbkt(ji,jj) - 2 > 0 so no incorrect array access 857 850 IF ( zhbl_t(ji,jj) > gdepw(ji, jj,mbkt(ji,jj)-2,Kmm) ) THEN 858 851 zhbl_t(ji,jj) = MIN( zhbl_t(ji,jj), gdepw(ji,jj,mbkt(ji,jj)-2,Kmm) ) ! ht(:,:)) 859 l_pyc(ji,jj) = .FALSE.852 l_pyc(ji,jj) = .FALSE. 860 853 l_coup(ji,jj) = .TRUE. ! ag 19/03 861 854 END IF … … 870 863 #endif 871 864 END_2D 872 865 ! 873 866 nmld(:,:) = nbld(:,:) ! use nmld to hold previous blayer index 874 867 nbld(:,:) = 4 875 868 ! 876 869 DO_3D( 0, 0, 0, 0, 4, jpkm1 ) 877 870 IF ( zhbl_t(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN 878 871 nbld(ji,jj) = jk 879 END IF872 END IF 880 873 END_3D 881 874 ! 882 875 ! 883 876 ! Step through model levels taking account of buoyancy change to determine the effect on dhdt 884 877 ! 885 CALL zdf_osm_timestep_hbl( Kmm, zdhdt, zhbl, zhbl_t, zt_bl, &886 & z s_bl, zwb_ent, zwb_fk_b )887 ! is external level in bounds?888 889 ! 878 CALL zdf_osm_timestep_hbl( Kmm, zdhdt, zhbl, zhbl_t, zwb_ent, & 879 & zwb_fk_b ) 880 ! Is external level in bounds? 881 ! 882 ! Recalculate BL averages and differences using new BL depth 890 883 jp_ext(:,:) = 1 ! ag 19/03 891 CALL zdf_osm_vertical_average( Kbb, Kmm, nbld, zt_bl, zs_bl, &892 & zb_bl, zu_bl, zv_bl,jp_ext, zdt_bl, &884 CALL zdf_osm_vertical_average( Kbb, Kmm, nbld, av_t_bl, av_s_bl, & 885 & av_b_bl, av_u_bl, av_v_bl, jp_ext, zdt_bl, & 893 886 & zds_bl, zdb_bl, zdu_bl, zdv_bl ) 894 895 CALL zdf_osm_pycnocline_thickness( Kmm, zdh, zhml,zdhdt, zdb_bl, &887 ! 888 CALL zdf_osm_pycnocline_thickness( Kmm, zdh, zhml, zdhdt, zdb_bl, & 896 889 & zhbl, zwb_ent, zdbdz_bl_ext, zwb_fk_b ) 897 890 ! 898 891 ! Reset l_pyc before calculating terms in the flux-gradient relationship 899 900 892 DO_2D( 0, 0, 0, 0 ) 901 IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh . or. nbld(ji,jj) >= mbkt(ji,jj) - 2 .or. &902 & nbld(ji,jj) - nmld(ji,jj) == 1 .or. zdhdt(ji,jj) < 0.0_wp ) THEN! ag 19/03903 l_pyc(ji,jj) = .FALSE. 893 IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh .OR. nbld(ji,jj) >= mbkt(ji,jj) - 2 .OR. & 894 & nbld(ji,jj) - nmld(ji,jj) == 1 .OR. zdhdt(ji,jj) < 0.0_wp ) THEN ! ag 19/03 895 l_pyc(ji,jj) = .FALSE. ! ag 19/03 904 896 IF ( nbld(ji,jj) >= mbkt(ji,jj) -2 ) THEN 905 nmld(ji,jj) = nbld(ji,jj) - 1 ! ag 19/03906 zdh(ji,jj) = gdepw(ji,jj,nbld(ji,jj),Kmm) - gdepw(ji,jj,nmld(ji,jj),Kmm) ! ag 19/03907 zhml(ji,jj) = gdepw(ji,jj,nmld(ji,jj),Kmm) ! ag 19/03908 dh(ji,jj) = zdh(ji,jj)! ag 19/03909 hml(ji,jj) = hbl(ji,jj) - dh(ji,jj)! ag 19/03897 nmld(ji,jj) = nbld(ji,jj) - 1 ! ag 19/03 898 zdh(ji,jj) = gdepw(ji,jj,nbld(ji,jj),Kmm) - gdepw(ji,jj,nmld(ji,jj),Kmm) ! ag 19/03 899 zhml(ji,jj) = gdepw(ji,jj,nmld(ji,jj),Kmm) ! ag 19/03 900 dh(ji,jj) = zdh(ji,jj) ! ag 19/03 901 hml(ji,jj) = hbl(ji,jj) - dh(ji,jj) ! ag 19/03 910 902 #ifdef key_osm_debug 911 903 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN … … 919 911 ENDIF ! ag 19/03 920 912 END_2D 921 922 dstokes(:,:) = MIN ( dstokes(:,:), hbl(:,:)/3. ) ! Limit delta for shallow boundary layers for calculating flux-gradient terms. 913 ! 914 dstokes(:,:) = MIN ( dstokes(:,:), hbl(:,:)/ 3.0_wp ) ! Limit delta for shallow boundary layers for calculating 915 ! ! flux-gradient terms 923 916 ! 924 917 ! Average over the depth of the mixed layer in the convective boundary layer 925 ! jp_ext = nbld - nmld + 1926 ! 918 ! jp_ext = nbld - nmld + 1 919 ! Recalculate ML averages and differences using new ML depth 927 920 jp_ext(:,:) = nbld(:,:) - nmld(:,:) + jp_ext(:,:) + 1 ! ag 19/03 928 CALL zdf_osm_vertical_average( Kbb, Kmm, nmld - 1, zt_ml, zs_ml, &929 & zb_ml, zu_ml, zv_ml,jp_ext, zdt_ml, &930 & zds_ml, zdb_ml, zdu_ml, 931 921 CALL zdf_osm_vertical_average( Kbb, Kmm, nmld - 1, av_t_ml, av_s_ml, & 922 & av_b_ml, av_u_ml, av_v_ml, jp_ext, zdt_ml, & 923 & zds_ml, zdb_ml, zdu_ml, zdv_ml ) 924 ! 932 925 CALL zdf_osm_external_gradients( Kmm, nbld + 1, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext ) 933 926 #ifdef key_osm_debug … … 935 928 ji=iloc_db; jj=jloc_db 936 929 WRITE(narea+100,'(4(3(a,g11.3),/), 2(4(a,g11.3),/))') & 937 & 'After averaging, with new hbl (&correct jp_ext), hml: zt_bl=', zt_bl(ji,jj),&938 & ' zs_bl=', zs_bl(ji,jj), ' zb_bl=', zb_bl(ji,jj),&930 & 'After averaging, with new hbl (&correct jp_ext), hml: zt_bl=', av_t_bl(ji,jj),& 931 & ' zs_bl=', av_s_bl(ji,jj), ' zb_bl=', av_b_bl(ji,jj),& 939 932 & 'zdt_bl=', zdt_bl(ji,jj), ' zds_bl=', zds_bl(ji,jj), ' zdb_bl=', zdb_bl(ji,jj),& 940 & 'zt_ml=', zt_ml(ji,jj), ' zs_ml=', zs_ml(ji,jj), ' zb_ml=', zb_ml(ji,jj),&933 & 'zt_ml=', av_t_ml(ji,jj), ' zs_ml=', av_s_ml(ji,jj), ' zb_ml=', av_b_ml(ji,jj),& 941 934 & 'zdt_ml=', zdt_ml(ji,jj), ' zds_ml=', zds_ml(ji,jj), ' zdb_ml=', zdb_ml(ji,jj),& 942 & 'zu_bl =', zu_bl(ji,jj) , ' zv_bl=', zv_bl(ji,jj), ' zdu_bl=', zdu_bl(ji,jj), ' zdv_bl=', zdv_bl(ji,jj),&943 & 'zu_ml =', zu_ml(ji,jj) , ' zv_ml=', zv_ml(ji,jj), ' zdu_ml=', zdu_ml(ji,jj), ' zdv_ml=', zdv_ml(ji,jj)935 & 'zu_bl =', av_u_bl(ji,jj) , ' zv_bl=', av_v_bl(ji,jj), ' zdu_bl=', zdu_bl(ji,jj), ' zdv_bl=', zdv_bl(ji,jj),& 936 & 'zu_ml =', av_u_ml(ji,jj) , ' zv_ml=', av_v_ml(ji,jj), ' zdu_ml=', zdu_ml(ji,jj), ' zdv_ml=', zdv_ml(ji,jj) 944 937 FLUSH(narea+100) 945 938 END IF 946 939 #endif 947 940 ! Rotate mean currents and changes onto wind aligned co-ordinates 948 CALL zdf_osm_velocity_rotation( zu_ml, zv_ml)941 CALL zdf_osm_velocity_rotation( av_u_ml, av_v_ml ) 949 942 CALL zdf_osm_velocity_rotation( zdu_ml, zdv_ml ) 950 CALL zdf_osm_velocity_rotation( zu_bl, zv_bl)943 CALL zdf_osm_velocity_rotation( av_u_bl, av_v_bl ) 951 944 CALL zdf_osm_velocity_rotation( zdu_bl, zdv_bl ) 952 945 #ifdef key_osm_debug … … 955 948 WRITE(narea+100,'(a,/, 2(4(a,g11.3),/))') & 956 949 & 'After rotation, with new hbl (& correct jp_ext), hml:', & 957 & 'zu_bl =', zu_bl(ji,jj) , ' zv_bl=', zv_bl(ji,jj), ' zdu_bl=', zdu_bl(ji,jj), ' zdv_bl=', zdv_bl(ji,jj),&958 & 'zu_ml =', zu_ml(ji,jj) , ' zv_ml=', zv_ml(ji,jj), ' zdu_ml=', zdu_ml(ji,jj), ' zdv_ml=', zdv_ml(ji,jj)950 & 'zu_bl =', av_u_bl(ji,jj) , ' zv_bl=', av_v_bl(ji,jj), ' zdu_bl=', zdu_bl(ji,jj), ' zdv_bl=', zdv_bl(ji,jj),& 951 & 'zu_ml =', av_u_ml(ji,jj) , ' zv_ml=', av_v_ml(ji,jj), ' zdu_ml=', zdu_ml(ji,jj), ' zdv_ml=', zdv_ml(ji,jj) 959 952 FLUSH(narea+100) 960 953 END IF … … 963 956 ! Pycnocline gradients for scalars and velocity 964 957 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 965 966 958 jp_ext(:,:) = 1 ! ag 19/03 967 CALL zdf_osm_pycnocline_buoyancy_profiles( Kmm, jp_ext, zdbdz_pyc, zalpha_pyc, zdh,&968 & zdb_bl, zhbl, zdbdz_bl_ext, zhml,zdb_ml, &959 CALL zdf_osm_pycnocline_buoyancy_profiles( Kmm, jp_ext, zdbdz_pyc, zalpha_pyc, zdh, & 960 & zdb_bl, zhbl, zdbdz_bl_ext, zhml, zdb_ml, & 969 961 & zdhdt ) 970 962 #ifdef key_osm_debug … … 988 980 ! Eddy viscosity/diffusivity and non-gradient terms in the flux-gradient relationship 989 981 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 990 CALL zdf_osm_diffusivity_viscosity( Kbb, Kmm, zdiffut, zviscos, zhbl, & 991 & zhml, zdh, zdhdt, zshear, zb_bl, & 992 & zu_bl, zv_bl, zb_ml, zu_ml, zv_ml, & 993 & zwb_ent, zwb_min ) 982 CALL zdf_osm_diffusivity_viscosity( Kbb, Kmm, zdiffut, zviscos, zhbl, & 983 & zhml, zdh, zdhdt, zshear, zwb_ent, & 984 & zwb_min ) 994 985 #ifdef key_osm_debug 995 986 IF(narea==nn_narea_db) THEN … … 1002 993 END IF 1003 994 #endif 1004 1005 995 ! 1006 996 ! Calculate non-gradient components of the flux-gradient relationships 1007 997 ! -------------------------------------------------------------------- 1008 CALL zdf_osm_fgr_terms( Kmm, jp_ext, zhbl, zhml, zdh,&1009 & zdhdt, zshear, zdt_bl, zds_bl, zdb_bl,&1010 & zdu_bl, zdv_bl, zdt_ml, zds_ml, zdb_ml,&1011 & zdu_ml, zdv_ml,zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_pyc, &998 CALL zdf_osm_fgr_terms( Kmm, jp_ext, zhbl, zhml, zdh, & 999 & zdhdt, zshear, zdt_bl, zds_bl, zdb_bl, & 1000 & zdu_bl, zdv_bl, zdt_ml, zds_ml, zdb_ml, & 1001 & zdu_ml, zdv_ml, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_pyc, & 1012 1002 & zalpha_pyc, zdiffut, zviscos ) 1013 1003 ! 1014 1004 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 1015 1005 ! Need to put in code for contributions that are applied explicitly to 1016 1006 ! the prognostic variables 1017 1007 ! 1. Entrainment flux 1018 !1019 1008 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 1020 1009 ! 1021 1010 ! Rotate non-gradient velocity terms back to model reference frame 1022 1011 CALL zdf_osm_velocity_rotation( ghamu, ghamv, .FALSE., 2, nbld ) 1023 1012 ! 1024 1013 ! KPP-style Ri# mixing 1025 1014 IF ( ln_kpprimix ) THEN … … 1031 1020 ! Shear production at uw- and vw-points (energy conserving form) 1032 1021 DO_2D( 1, 0, 1, 0 ) 1033 z2du(ji,jj) = 0.5_wp * ( uu(ji,jj,jk-1,Kmm) - uu(ji,jj,jk,Kmm) ) * & 1034 & ( uu(ji,jj,jk-1,Kbb) - uu(ji,jj,jk,Kbb) ) * wumask(ji,jj,jk) / & 1035 & ( e3uw(ji,jj,jk,Kmm) * e3uw(ji,jj,jk,Kbb) ) 1036 z2dv(ji,jj) = 0.5_wp * ( vv(ji,jj,jk-1,Kmm) - vv(ji,jj,jk,Kmm) ) * & 1037 & ( vv(ji,jj,jk-1,Kbb) - vv(ji,jj,jk,Kbb) ) * wvmask(ji,jj,jk) / & 1038 & ( e3vw(ji,jj,jk,Kmm) * e3vw(ji,jj,jk,Kbb) ) 1022 z2du(ji,jj) = 0.5_wp * ( uu(ji,jj,jk-1,Kmm) - uu(ji,jj,jk,Kmm) ) * ( uu(ji,jj,jk-1,Kbb) - uu(ji,jj,jk,Kbb) ) * & 1023 & wumask(ji,jj,jk) / ( e3uw(ji,jj,jk,Kmm) * e3uw(ji,jj,jk,Kbb) ) 1024 z2dv(ji,jj) = 0.5_wp * ( vv(ji,jj,jk-1,Kmm) - vv(ji,jj,jk,Kmm) ) * ( vv(ji,jj,jk-1,Kbb) - vv(ji,jj,jk,Kbb) ) * & 1025 & wvmask(ji,jj,jk) / ( e3vw(ji,jj,jk,Kmm) * e3vw(ji,jj,jk,Kbb) ) 1039 1026 END_2D 1040 1027 DO_2D( 0, 0, 0, 0 ) 1041 1028 IF ( jk > nbld(ji,jj) ) THEN 1042 1029 ! Shear prod. at w-point weightened by mask 1043 zesh2 = ( z2du(ji-1,jj) + z2du(ji,jj) ) / MAX( 1.0_wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )&1044 & +( z2dv(ji,jj-1) + z2dv(ji,jj) ) / MAX( 1.0_wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )1030 zesh2 = ( z2du(ji-1,jj) + z2du(ji,jj) ) / MAX( 1.0_wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) + & 1031 & ( z2dv(ji,jj-1) + z2dv(ji,jj) ) / MAX( 1.0_wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) ) 1045 1032 ! Local Richardson number 1046 zri = MAX( rn2b(ji,jj,jk), 0.0_wp ) / MAX(zesh2, epsln)1047 zfri = MIN( zri / rn_riinfty, 1.0_wp )1048 zfri = ( 1.0_wp - zfri * zfri )1033 zri = MAX( rn2b(ji,jj,jk), 0.0_wp ) / MAX( zesh2, epsln ) 1034 zfri = MIN( zri / rn_riinfty, 1.0_wp ) 1035 zfri = ( 1.0_wp - zfri * zfri ) 1049 1036 zrimix = zfri * zfri * zfri * wmask(ji, jj, jk) 1050 1037 zdiffut(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), zrimix*rn_difri ) … … 1053 1040 END_2D 1054 1041 END DO 1055 END IF ! ln_kpprimix = .true.1056 1042 END IF ! ln_kpprimix = .true. 1043 ! 1057 1044 ! KPP-style set diffusivity large if unstable below BL 1058 IF ( ln_convmix) THEN1045 IF ( ln_convmix) THEN 1059 1046 DO_2D( 0, 0, 0, 0 ) 1060 1047 DO jk = nbld(ji,jj) + 1, jpkm1 1061 IF ( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12) zdiffut(ji,jj,jk) = MAX( rn_difconv, zdiffut(ji,jj,jk) )1048 IF ( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1e-12_wp ) zdiffut(ji,jj,jk) = MAX( rn_difconv, zdiffut(ji,jj,jk) ) 1062 1049 END DO 1063 1050 END_2D 1064 END IF ! ln_convmix = .true.1051 END IF ! ln_convmix = .true. 1065 1052 #ifdef key_osm_debug 1066 1053 IF(narea==nn_narea_db) THEN … … 1074 1061 END IF 1075 1062 #endif 1076 1077 1078 1079 IF ( ln_osm_mle ) THEN ! set up diffusivity and non-gradient mixing 1063 ! 1064 IF ( ln_osm_mle ) THEN ! Set up diffusivity and non-gradient mixing 1080 1065 DO_2D( 0, 0, 0, 0 ) 1081 IF ( l_flux(ji,jj) ) THEN ! MLE mixing extends below boundary layer1066 IF ( l_flux(ji,jj) ) THEN ! MLE mixing extends below boundary layer 1082 1067 ! Calculate MLE flux contribution from surface fluxes 1083 1068 DO jk = 1, nbld(ji,jj) 1084 znd = gdepw(ji,jj,jk,Kmm) / MAX( zhbl(ji,jj),epsln)1085 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) - ( swth0(ji,jj) - zrad0(ji,jj) + zradh(ji,jj) ) * ( 1.0 - znd )1086 ghams(ji,jj,jk) = ghams(ji,jj,jk) - sws0(ji,jj) * ( 1.0 - znd )1069 znd = gdepw(ji,jj,jk,Kmm) / MAX( zhbl(ji,jj), epsln ) 1070 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) - ( swth0(ji,jj) - zrad0(ji,jj) + zradh(ji,jj) ) * ( 1.0_wp - znd ) 1071 ghams(ji,jj,jk) = ghams(ji,jj,jk) - sws0(ji,jj) * ( 1.0_wp - znd ) 1087 1072 END DO 1088 1073 DO jk = 1, mld_prof(ji,jj) 1089 znd = gdepw(ji,jj,jk,Kmm) / MAX( zhmle(ji,jj),epsln)1090 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + ( swth0(ji,jj) - zrad0(ji,jj) + zradh(ji,jj) ) * ( 1.0 - znd )1091 ghams(ji,jj,jk) = ghams(ji,jj,jk) + sws0(ji,jj) * ( 1.0 -znd )1074 znd = gdepw(ji,jj,jk,Kmm) / MAX( zhmle(ji,jj), epsln ) 1075 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + ( swth0(ji,jj) - zrad0(ji,jj) + zradh(ji,jj) ) * ( 1.0_wp - znd ) 1076 ghams(ji,jj,jk) = ghams(ji,jj,jk) + sws0(ji,jj) * ( 1.0_wp -znd ) 1092 1077 END DO 1093 1078 ! Viscosity for MLEs 1094 1079 DO jk = 1, mld_prof(ji,jj) 1095 znd = -gdepw(ji,jj,jk,Kmm) / MAX(zhmle(ji,jj),epsln) 1096 zdiffut(ji,jj,jk) = zdiffut(ji,jj,jk) + zdiff_mle(ji,jj) * ( 1.0 - ( 2.0 * znd + 1.0 )**2 ) * ( 1.0 + 5.0 / 21.0 * ( 2.0 * znd + 1.0 )** 2 ) 1080 znd = -1.0_wp * gdepw(ji,jj,jk,Kmm) / MAX( zhmle(ji,jj), epsln ) 1081 zdiffut(ji,jj,jk) = zdiffut(ji,jj,jk) + zdiff_mle(ji,jj) * ( 1.0_wp - ( 2.0_wp * znd + 1.0_wp )**2 ) * & 1082 & ( 1.0_wp + 5.0_wp / 21.0_wp * ( 2.0_wp * znd + 1.0_wp )**2 ) 1097 1083 END DO 1098 ELSE 1099 ! Surface transports limited to OSBL. 1084 ELSE ! Surface transports limited to OSBL 1100 1085 ! Viscosity for MLEs 1101 1086 DO jk = 1, mld_prof(ji,jj) 1102 znd = -gdepw(ji,jj,jk,Kmm) / MAX(zhmle(ji,jj),epsln) 1103 zdiffut(ji,jj,jk) = zdiffut(ji,jj,jk) + zdiff_mle(ji,jj) * ( 1.0 - ( 2.0 * znd + 1.0 )**2 ) * ( 1.0 + 5.0 / 21.0 * ( 2.0 * znd + 1.0 )** 2 ) 1087 znd = -1.0_wp * gdepw(ji,jj,jk,Kmm) / MAX( zhmle(ji,jj), epsln ) 1088 zdiffut(ji,jj,jk) = zdiffut(ji,jj,jk) + zdiff_mle(ji,jj) * ( 1.0_wp - ( 2.0_wp * znd + 1.0_wp )**2 ) * & 1089 & ( 1.0_wp + 5.0_wp / 21.0_wp * ( 2.0_wp * znd + 1.0_wp )**2 ) 1104 1090 END DO 1105 END IF1091 END IF 1106 1092 END_2D 1107 1093 #ifdef key_osm_debug … … 1118 1104 #endif 1119 1105 ENDIF 1120 1106 ! 1121 1107 ! Lateral boundary conditions on zvicos (sign unchanged), needed to caclulate viscosities on u and v grids 1122 !CALL lbc_lnk( 'zdfosm', zviscos(:,:,:), 'W', 1.0_wp ) 1123 1108 ! CALL lbc_lnk( 'zdfosm', zviscos(:,:,:), 'W', 1.0_wp ) 1124 1109 ! GN 25/8: need to change tmask --> wmask 1125 1126 1110 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 1127 1111 p_avt(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), avtb(jk) ) * tmask(ji,jj,jk) 1128 1112 p_avm(ji,jj,jk) = MAX( zviscos(ji,jj,jk), avmb(jk) ) * tmask(ji,jj,jk) 1129 1113 END_3D 1130 ! Lateral boundary conditions on ghamu and ghamv, currently on W-grid (sign unchanged), needed to caclulate gham[uv] on u and v grids 1131 CALL lbc_lnk_multi( 'zdfosm', p_avt, 'W', 1.0_wp , p_avm, 'W', 1.0_wp, & 1132 & ghamu, 'W', 1.0_wp , ghamv, 'W', 1.0_wp ) 1114 ! Lateral boundary conditions on ghamu and ghamv, currently on W-grid (sign unchanged), needed to caclulate gham[uv] on u and 1115 ! v grids 1116 CALL lbc_lnk_multi( 'zdfosm', p_avt, 'W', 1.0_wp, & 1117 & p_avm, 'W', 1.0_wp, & 1118 & ghamu, 'W', 1.0_wp, & 1119 & ghamv, 'W', 1.0_wp ) 1133 1120 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 1134 ghamu(ji,jj,jk) = ( ghamu(ji,jj,jk) + ghamu(ji+1,jj,jk) ) & 1135 & / MAX( 1., tmask(ji,jj,jk) + tmask (ji + 1,jj,jk) ) * umask(ji,jj,jk) 1136 1137 ghamv(ji,jj,jk) = ( ghamv(ji,jj,jk) + ghamv(ji,jj+1,jk) ) & 1138 & / MAX( 1., tmask(ji,jj,jk) + tmask (ji,jj+1,jk) ) * vmask(ji,jj,jk) 1139 1121 ghamu(ji,jj,jk) = ( ghamu(ji,jj,jk) + ghamu(ji+1,jj,jk) ) / MAX( 1.0_wp, tmask(ji,jj,jk) + tmask (ji+1,jj,jk) ) * & 1122 & umask(ji,jj,jk) 1123 ghamv(ji,jj,jk) = ( ghamv(ji,jj,jk) + ghamv(ji,jj+1,jk) ) / MAX( 1.0_wp, tmask(ji,jj,jk) + tmask (ji,jj+1,jk) ) * & 1124 & vmask(ji,jj,jk) 1140 1125 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) * tmask(ji,jj,jk) 1141 1126 ghams(ji,jj,jk) = ghams(ji,jj,jk) * tmask(ji,jj,jk) 1142 1127 END_3D 1143 ! Lateral boundary conditions on final outputs for hbl, on T-grid (sign unchanged) 1144 CALL lbc_lnk_multi( 'zdfosm', hbl, 'T', 1., dh, 'T', 1., hmle, 'T', 1. ) 1145 ! Lateral boundary conditions on final outputs for gham[ts], on W-grid (sign unchanged) 1146 ! Lateral boundary conditions on final outputs for gham[uv], on [UV]-grid (sign changed) 1147 CALL lbc_lnk_multi( 'zdfosm', ghamt, 'W', 1.0_wp , ghams, 'W', 1.0_wp, & 1148 & ghamu, 'U', -1.0_wp , ghamv, 'V', -1.0_wp ) 1128 ! Lateral boundary conditions on final outputs for hbl, on T-grid (sign unchanged) 1129 CALL lbc_lnk_multi( 'zdfosm', hbl, 'T', 1.0_wp, & 1130 & dh, 'T', 1.0_wp, & 1131 & hmle, 'T', 1.0_wp ) 1132 ! Lateral boundary conditions on final outputs for gham[ts], on W-grid (sign unchanged) 1133 ! Lateral boundary conditions on final outputs for gham[uv], on [UV]-grid (sign changed) 1134 CALL lbc_lnk_multi( 'zdfosm', ghamt, 'W', 1.0_wp, & 1135 & ghams, 'W', 1.0_wp, & 1136 & ghamu, 'U', -1.0_wp, & 1137 & ghamv, 'V', -1.0_wp ) 1149 1138 #ifdef key_osm_debug 1150 1139 IF(narea==nn_narea_db) THEN … … 1162 1151 END IF 1163 1152 #endif 1164 1153 ! 1165 1154 IF(ln_dia_osm) THEN 1166 1155 SELECT CASE (nn_osm_wave) 1167 ! Stokes drift set by assumimg onstant La#=0.3 (=0) or Pierson-Moskovitz spectrum (=1).1156 ! Stokes drift set by assumimg onstant La#=0.3 (=0) or Pierson-Moskovitz spectrum (=1) 1168 1157 CASE(0:1) 1169 IF ( iom_use("us_x") ) CALL iom_put( "us_x", tmask(:,:,1)*sustke*scos_wind ) ! x surface Stokes drift 1170 IF ( iom_use("us_y") ) CALL iom_put( "us_y", tmask(:,:,1)*sustke*ssin_wind ) ! y surface Stokes drift 1171 IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rho0*tmask(:,:,1)*sustar**2*sustke ) 1158 IF ( iom_use("us_x") ) CALL iom_put( "us_x", tmask(:,:,1)*sustke*scos_wind ) ! x surface Stokes drift 1159 IF ( iom_use("us_y") ) CALL iom_put( "us_y", tmask(:,:,1)*sustke*ssin_wind ) ! y surface Stokes drift 1160 IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.0_wp * rho0 * tmask(:,:,1) * & 1161 & sustar**2 * sustke ) 1172 1162 ! Stokes drift read in from sbcwave (=2). 1173 1163 CASE(2:3) 1174 IF ( iom_use("us_x") ) CALL iom_put( "us_x", ut0sd*umask(:,:,1) ) ! x surface Stokes drift 1175 IF ( iom_use("us_y") ) CALL iom_put( "us_y", vt0sd*vmask(:,:,1) ) ! y surface Stokes drift 1176 IF ( iom_use("wmp") ) CALL iom_put( "wmp", wmp*tmask(:,:,1) ) ! wave mean period 1177 IF ( iom_use("hsw") ) CALL iom_put( "hsw", hsw*tmask(:,:,1) ) ! significant wave height 1178 IF ( iom_use("wmp_NP") ) CALL iom_put( "wmp_NP", (2.*rpi*1.026/(0.877*grav) )*wndm*tmask(:,:,1) ) ! wave mean period from NP spectrum 1179 IF ( iom_use("hsw_NP") ) CALL iom_put( "hsw_NP", (0.22/grav)*wndm**2*tmask(:,:,1) ) ! significant wave height from NP spectrum 1180 IF ( iom_use("wndm") ) CALL iom_put( "wndm", wndm*tmask(:,:,1) ) ! U_10 1181 IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rho0*tmask(:,:,1)*sustar**2* & 1182 & SQRT(ut0sd**2 + vt0sd**2 ) ) 1164 IF ( iom_use("us_x") ) CALL iom_put( "us_x", ut0sd*umask(:,:,1) ) ! x surface Stokes drift 1165 IF ( iom_use("us_y") ) CALL iom_put( "us_y", vt0sd*vmask(:,:,1) ) ! y surface Stokes drift 1166 IF ( iom_use("wmp") ) CALL iom_put( "wmp", wmp*tmask(:,:,1) ) ! Wave mean period 1167 IF ( iom_use("hsw") ) CALL iom_put( "hsw", hsw*tmask(:,:,1) ) ! Significant wave height 1168 IF ( iom_use("wmp_NP") ) CALL iom_put( "wmp_NP", ( 2.0_wp * rpi * 1.026_wp / & ! Wave mean period from NP 1169 & ( 0.877_wp * grav ) ) * & ! spectrum 1170 & wndm * tmask(:,:,1) ) 1171 IF ( iom_use("hsw_NP") ) CALL iom_put( "hsw_NP", ( 0.22_wp / grav ) * wndm**2 * & ! Significant wave 1172 & tmask(:,:,1) ) ! height from NP spectrum 1173 IF ( iom_use("wndm") ) CALL iom_put( "wndm", wndm*tmask(:,:,1) ) ! U_10 1174 IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", & 1175 & 1000.0_wp * rho0 * tmask(:,:,1) * sustar**2 * & 1176 & SQRT( ut0sd**2 + vt0sd**2 ) ) 1183 1177 END SELECT 1184 IF ( iom_use("ghamt") ) CALL iom_put( "ghamt", tmask*ghamt ) ! <Tw_NL>1185 IF ( iom_use("ghams") ) CALL iom_put( "ghams", tmask*ghams ) ! <Sw_NL>1186 IF ( iom_use("ghamu") ) CALL iom_put( "ghamu", umask*ghamu ) ! <uw_NL>1187 IF ( iom_use("ghamv") ) CALL iom_put( "ghamv", vmask*ghamv ) ! <vw_NL>1188 IF ( iom_use("zwth0") ) CALL iom_put( "zwth0", tmask(:,:,1)*swth0 ) ! <Tw_0>1189 IF ( iom_use("zws0") ) CALL iom_put( "zws0", tmask(:,:,1)*sws0 ) ! <Sw_0>1190 IF ( iom_use("zwb0") ) CALL iom_put( "zwb0", tmask(:,:,1)*swb0 ) ! <Sw_0>1191 IF ( iom_use("zwbav") ) CALL iom_put( "zwbav", tmask(:,:,1)*swth0 ) ! upward BL-avged turb buoyancy flux1192 IF ( iom_use("hbl") ) CALL iom_put( "hbl", tmask(:,:,1)*hbl ) ! boundary-layer depth1193 IF ( iom_use("ibld") ) CALL iom_put( "ibld", tmask(:,:,1)*nbld ) ! boundary-layer max k1194 IF ( iom_use("zdt_bl") ) CALL iom_put( "zdt_bl", tmask(:,:,1)*zdt_bl ) ! dt at ml base1195 IF ( iom_use("zds_bl") ) CALL iom_put( "zds_bl", tmask(:,:,1)*zds_bl ) ! ds at ml base1196 IF ( iom_use("zdb_bl") ) CALL iom_put( "zdb_bl", tmask(:,:,1)*zdb_bl ) ! db at ml base1197 IF ( iom_use("zdu_bl") ) CALL iom_put( "zdu_bl", tmask(:,:,1)*zdu_bl ) ! du at ml base1198 IF ( iom_use("zdv_bl") ) CALL iom_put( "zdv_bl", tmask(:,:,1)*zdv_bl ) ! dv at ml base1199 IF ( iom_use("dh") ) CALL iom_put( "dh", tmask(:,:,1)*dh ) ! Initial boundary-layer depth1200 IF ( iom_use("hml") ) CALL iom_put( "hml", tmask(:,:,1)*hml ) ! Initial boundary-layer depth1201 IF ( iom_use("zdt_ml") ) CALL iom_put( "zdt_ml", tmask(:,:,1)*zdt_ml ) ! dt at ml base1202 IF ( iom_use("zds_ml") ) CALL iom_put( "zds_ml", tmask(:,:,1)*zds_ml ) ! ds at ml base1203 IF ( iom_use("zdb_ml") ) CALL iom_put( "zdb_ml", tmask(:,:,1)*zdb_ml ) ! db at ml base1204 IF ( iom_use("dstokes") ) CALL iom_put( "dstokes", tmask(:,:,1)*dstokes ) ! Stokes drift penetration depth1178 IF ( iom_use("ghamt") ) CALL iom_put( "ghamt", tmask*ghamt ) ! <Tw_NL> 1179 IF ( iom_use("ghams") ) CALL iom_put( "ghams", tmask*ghams ) ! <Sw_NL> 1180 IF ( iom_use("ghamu") ) CALL iom_put( "ghamu", umask*ghamu ) ! <uw_NL> 1181 IF ( iom_use("ghamv") ) CALL iom_put( "ghamv", vmask*ghamv ) ! <vw_NL> 1182 IF ( iom_use("zwth0") ) CALL iom_put( "zwth0", tmask(:,:,1)*swth0 ) ! <Tw_0> 1183 IF ( iom_use("zws0") ) CALL iom_put( "zws0", tmask(:,:,1)*sws0 ) ! <Sw_0> 1184 IF ( iom_use("zwb0") ) CALL iom_put( "zwb0", tmask(:,:,1)*swb0 ) ! <Sw_0> 1185 IF ( iom_use("zwbav") ) CALL iom_put( "zwbav", tmask(:,:,1)*swth0 ) ! Upward BL-avged turb buoyancy flux 1186 IF ( iom_use("hbl") ) CALL iom_put( "hbl", tmask(:,:,1)*hbl ) ! Boundary-layer depth 1187 IF ( iom_use("ibld") ) CALL iom_put( "ibld", tmask(:,:,1)*nbld ) ! Boundary-layer max k 1188 IF ( iom_use("zdt_bl") ) CALL iom_put( "zdt_bl", tmask(:,:,1)*zdt_bl ) ! dt at ml base 1189 IF ( iom_use("zds_bl") ) CALL iom_put( "zds_bl", tmask(:,:,1)*zds_bl ) ! ds at ml base 1190 IF ( iom_use("zdb_bl") ) CALL iom_put( "zdb_bl", tmask(:,:,1)*zdb_bl ) ! db at ml base 1191 IF ( iom_use("zdu_bl") ) CALL iom_put( "zdu_bl", tmask(:,:,1)*zdu_bl ) ! du at ml base 1192 IF ( iom_use("zdv_bl") ) CALL iom_put( "zdv_bl", tmask(:,:,1)*zdv_bl ) ! dv at ml base 1193 IF ( iom_use("dh") ) CALL iom_put( "dh", tmask(:,:,1)*dh ) ! Initial boundary-layer depth 1194 IF ( iom_use("hml") ) CALL iom_put( "hml", tmask(:,:,1)*hml ) ! Initial boundary-layer depth 1195 IF ( iom_use("zdt_ml") ) CALL iom_put( "zdt_ml", tmask(:,:,1)*zdt_ml ) ! dt at ml base 1196 IF ( iom_use("zds_ml") ) CALL iom_put( "zds_ml", tmask(:,:,1)*zds_ml ) ! ds at ml base 1197 IF ( iom_use("zdb_ml") ) CALL iom_put( "zdb_ml", tmask(:,:,1)*zdb_ml ) ! db at ml base 1198 IF ( iom_use("dstokes") ) CALL iom_put( "dstokes", tmask(:,:,1)*dstokes ) ! Stokes drift penetration depth 1205 1199 IF ( iom_use("zustke") ) CALL iom_put( "zustke", tmask(:,:,1)*sustke ) ! Stokes drift magnitude at T-points 1206 IF ( iom_use("zwstrc") ) CALL iom_put( "zwstrc", tmask(:,:,1)*swstrc ) ! convective velocity scale1207 IF ( iom_use("zwstrl") ) CALL iom_put( "zwstrl", tmask(:,:,1)*swstrl ) ! Langmuir velocity scale1208 IF ( iom_use("zustar") ) CALL iom_put( "zustar", tmask(:,:,1)*sustar ) ! friction velocity scale1209 IF ( iom_use("zvstr") ) CALL iom_put( "zvstr", tmask(:,:,1)*svstr ) ! mixed velocity scale1210 IF ( iom_use("zla") ) CALL iom_put( "zla", tmask(:,:,1)*sla ) ! langmuir #1211 IF ( iom_use("wind_power") ) CALL iom_put( "wind_power", 1000. *rho0*tmask(:,:,1)*sustar**3 )! BL depth internal to zdf_osm routine1212 IF ( iom_use("wind_wave_power") ) CALL iom_put( "wind_wave_power", 1000.*rho0*tmask(:,:,1)*sustar**2*sustke )1213 IF ( iom_use(" zhbl") ) CALL iom_put( "zhbl", tmask(:,:,1)*zhbl ) ! BL depth internal to zdf_osm routine1214 IF ( iom_use("zh ml") ) CALL iom_put( "zhml", tmask(:,:,1)*zhml ) ! ML depth internal to zdf_osm routine1215 IF ( iom_use(" imld") ) CALL iom_put( "imld", tmask(:,:,1)*nmld ) ! index forML depth internal to zdf_osm routine1216 IF ( iom_use(" jp_ext") ) CALL iom_put( "jp_ext", tmask(:,:,1)*jp_ext ) ! =1 if pycnocline resolvedinternal to zdf_osm routine1217 IF ( iom_use("j _ddh") ) CALL iom_put( "j_ddh", tmask(:,:,1)*n_ddh ) ! index forpyc thicknesshinternal to zdf_osm routine1218 IF ( iom_use(" zshear") ) CALL iom_put( "zshear", tmask(:,:,1)*zshear ) ! shear production of TKEinternal to zdf_osm routine1219 IF ( iom_use("z dh") ) CALL iom_put( "zdh", tmask(:,:,1)*zdh ) ! pyc thicknesshinternal to zdf_osm routine1220 IF ( iom_use("z hol") ) CALL iom_put( "zhol", tmask(:,:,1)*shol ) ! ML depth internal to zdf_osm routine1221 IF ( iom_use("z wb_ent") ) CALL iom_put( "zwb_ent", tmask(:,:,1)*zwb_ent ) ! upward turb buoyancy entrainment flux1222 IF ( iom_use("z t_ml") ) CALL iom_put( "zt_ml", tmask(:,:,1)*zt_ml ) ! average T in ML1223 1224 IF ( iom_use("hmle") ) CALL iom_put( "hmle", tmask(:,:,1)*hmle ) ! FK layer depth1225 IF ( iom_use("zmld") ) CALL iom_put( "zmld", tmask(:,:,1)*zmld ) ! FK target layer depth1226 IF ( iom_use("zwb_fk") ) CALL iom_put( "zwb_fk", tmask(:,:,1)*zwb_fk ) ! FK b flux1227 IF ( iom_use("zwb_fk_b") ) CALL iom_put( "zwb_fk_b", tmask(:,:,1)*zwb_fk_b ) ! FK b flux averaged over ML1228 IF ( iom_use("mld_prof") ) CALL iom_put( "mld_prof", tmask(:,:,1)*mld_prof ) ! FK layer max k1229 IF ( iom_use("zdtdx") ) CALL iom_put( "zdtdx", umask(:,:,1)*zdtdx ) ! FK dtdx at u-pt1230 IF ( iom_use("zdtdy") ) CALL iom_put( "zdtdy", vmask(:,:,1)*zdtdy ) ! FK dtdy at v-pt1231 IF ( iom_use("zdsdx") ) CALL iom_put( "zdsdx", umask(:,:,1)*zdsdx ) ! FK dtdx at u-pt1232 IF ( iom_use("zdsdy") ) CALL iom_put( "zdsdy", vmask(:,:,1)*zdsdy ) ! FK dsdy at v-pt1233 IF ( iom_use("dbdx_mle") ) CALL iom_put( "dbdx_mle", umask(:,:,1)*dbdx_mle ) 1234 IF ( iom_use("dbdy_mle") ) CALL iom_put( "dbdy_mle", vmask(:,:,1)*dbdy_mle ) 1235 IF ( iom_use("zdiff_mle") ) CALL iom_put( "zdiff_mle", tmask(:,:,1)*zdiff_mle ) ! FK diff in MLE at t-pt1236 IF ( iom_use("zvel_mle") ) CALL iom_put( "zvel_mle", tmask(:,:,1)*zdiff_mle ) ! FK diff in MLE at t-pt1200 IF ( iom_use("zwstrc") ) CALL iom_put( "zwstrc", tmask(:,:,1)*swstrc ) ! Convective velocity scale 1201 IF ( iom_use("zwstrl") ) CALL iom_put( "zwstrl", tmask(:,:,1)*swstrl ) ! Langmuir velocity scale 1202 IF ( iom_use("zustar") ) CALL iom_put( "zustar", tmask(:,:,1)*sustar ) ! Friction velocity scale 1203 IF ( iom_use("zvstr") ) CALL iom_put( "zvstr", tmask(:,:,1)*svstr ) ! Mixed velocity scale 1204 IF ( iom_use("zla") ) CALL iom_put( "zla", tmask(:,:,1)*sla ) ! Langmuir # 1205 IF ( iom_use("wind_power") ) CALL iom_put( "wind_power", 1000.0_wp * rho0 * & ! BL depth internal to zdf_osm routine 1206 & tmask(:,:,1) * sustar**3 ) 1207 IF ( iom_use("wind_wave_power") ) CALL iom_put( "wind_wave_power", 1000.0_wp * rho0 * tmask(:,:,1) * sustar**2 * sustke ) 1208 IF ( iom_use("zhbl") ) CALL iom_put( "zhbl", tmask(:,:,1)*zhbl ) ! BL depth internal to zdf_osm routine 1209 IF ( iom_use("zhml") ) CALL iom_put( "zhml", tmask(:,:,1)*zhml ) ! ML depth internal to zdf_osm routine 1210 IF ( iom_use("imld") ) CALL iom_put( "imld", tmask(:,:,1)*nmld ) ! Index for ML depth internal to zdf_osm routine 1211 IF ( iom_use("jp_ext") ) CALL iom_put( "jp_ext", tmask(:,:,1)*jp_ext ) ! =1 if pycnocline resolved internal to zdf_osm routine 1212 IF ( iom_use("j_ddh") ) CALL iom_put( "j_ddh", tmask(:,:,1)*n_ddh ) ! Index forpyc thicknessh internal to zdf_osm routine 1213 IF ( iom_use("zshear") ) CALL iom_put( "zshear", tmask(:,:,1)*zshear ) ! Shear production of TKE internal to zdf_osm routine 1214 IF ( iom_use("zdh") ) CALL iom_put( "zdh", tmask(:,:,1)*zdh ) ! Pyc thicknessh internal to zdf_osm routine 1215 IF ( iom_use("zhol") ) CALL iom_put( "zhol", tmask(:,:,1)*shol ) ! ML depth internal to zdf_osm routine 1216 IF ( iom_use("zwb_ent") ) CALL iom_put( "zwb_ent", tmask(:,:,1)*zwb_ent ) ! Upward turb buoyancy entrainment flux 1217 IF ( iom_use("zt_ml") ) CALL iom_put( "zt_ml", tmask(:,:,1)*av_t_ml ) ! Average T in ML 1218 IF ( iom_use("hmle") ) CALL iom_put( "hmle", tmask(:,:,1)*hmle ) ! FK layer depth 1219 IF ( iom_use("zmld") ) CALL iom_put( "zmld", tmask(:,:,1)*zmld ) ! FK target layer depth 1220 IF ( iom_use("zwb_fk") ) CALL iom_put( "zwb_fk", tmask(:,:,1)*zwb_fk ) ! FK b flux 1221 IF ( iom_use("zwb_fk_b") ) CALL iom_put( "zwb_fk_b", tmask(:,:,1)*zwb_fk_b ) ! FK b flux averaged over ML 1222 IF ( iom_use("mld_prof") ) CALL iom_put( "mld_prof", tmask(:,:,1)*mld_prof ) ! FK layer max k 1223 IF ( iom_use("zdtdx") ) CALL iom_put( "zdtdx", umask(:,:,1)*zdtdx ) ! FK dtdx at u-pt 1224 IF ( iom_use("zdtdy") ) CALL iom_put( "zdtdy", vmask(:,:,1)*zdtdy ) ! FK dtdy at v-pt 1225 IF ( iom_use("zdsdx") ) CALL iom_put( "zdsdx", umask(:,:,1)*zdsdx ) ! FK dtdx at u-pt 1226 IF ( iom_use("zdsdy") ) CALL iom_put( "zdsdy", vmask(:,:,1)*zdsdy ) ! FK dsdy at v-pt 1227 IF ( iom_use("dbdx_mle") ) CALL iom_put( "dbdx_mle", umask(:,:,1)*dbdx_mle ) ! FK dbdx at u-pt 1228 IF ( iom_use("dbdy_mle") ) CALL iom_put( "dbdy_mle", vmask(:,:,1)*dbdy_mle ) ! FK dbdy at v-pt 1229 IF ( iom_use("zdiff_mle") ) CALL iom_put( "zdiff_mle", tmask(:,:,1)*zdiff_mle ) ! FK diff in MLE at t-pt 1230 IF ( iom_use("zvel_mle") ) CALL iom_put( "zvel_mle", tmask(:,:,1)*zdiff_mle ) ! FK diff in MLE at t-pt 1237 1231 1238 1232 END IF … … 1241 1235 END SUBROUTINE zdf_osm 1242 1236 1243 SUBROUTINE zdf_osm_vertical_average( Kbb, Kmm, knlev, pt, 1244 & pb, pu, pv,kp_ext, pdt, &1245 & pds, pdb, pdu, 1237 SUBROUTINE zdf_osm_vertical_average( Kbb, Kmm, knlev, pt, ps, & 1238 & pb, pu, pv, kp_ext, pdt, & 1239 & pds, pdb, pdu, pdv ) 1246 1240 !!--------------------------------------------------------------------- 1247 1241 !! *** ROUTINE zdf_vertical_average *** … … 1431 1425 END SUBROUTINE zdf_osm_velocity_rotation_3d 1432 1426 1433 SUBROUTINE zdf_osm_osbl_state( Kmm, pwb_ent, pwb_min, pshear, phbl,&1434 & phml, pdh, pdb_bl, pdu_bl, pdv_bl,&1435 & pdb_ml, pdu_ml, 1427 SUBROUTINE zdf_osm_osbl_state( Kmm, pwb_ent, pwb_min, pshear, phbl, & 1428 & phml, pdh, pdb_bl, pdu_bl, pdv_bl, & 1429 & pdb_ml, pdu_ml, pdv_ml ) 1436 1430 !!--------------------------------------------------------------------- 1437 1431 !! *** ROUTINE zdf_osm_osbl_state *** … … 1467 1461 REAL(wp) :: zwcor, zrf_conv, zrf_shear, zrf_langmuir, zr_stokes 1468 1462 ! 1469 REAL(wp), PARAMETER :: za_shr = 0.4_wp, zb_shr = 6.5_wp, za_wb_s = 0.8_wp 1470 REAL(wp), PARAMETER :: zalpha_c = 0.2_wp, zalpha_lc = 0.03_wp 1471 REAL(wp), PARAMETER :: zalpha_ls = 0.06_wp, zalpha_s = 0.15_wp 1472 REAL(wp), PARAMETER :: rn_ri_p_thresh = 27.0_wp 1473 REAL(wp), PARAMETER :: zri_c = 0.25_wp 1474 REAL(wp), PARAMETER :: zek = 4.0_wp 1475 REAL(wp), PARAMETER :: zrot = 0.0_wp ! Dummy rotation rate of surface stress 1476 REAL(wp), PARAMETER :: zlarge = -1e10_wp 1463 REAL(wp), PARAMETER :: pp_a_shr = 0.4_wp, pp_b_shr = 6.5_wp, pp_a_wb_s = 0.8_wp 1464 REAL(wp), PARAMETER :: pp_alpha_c = 0.2_wp, pp_alpha_lc = 0.03_wp 1465 REAL(wp), PARAMETER :: pp_alpha_ls = 0.06_wp, pp_alpha_s = 0.15_wp 1466 REAL(wp), PARAMETER :: pp_ri_p_thresh = 27.0_wp 1467 REAL(wp), PARAMETER :: pp_ri_c = 0.25_wp 1468 REAL(wp), PARAMETER :: pp_ek = 4.0_wp 1469 REAL(wp), PARAMETER :: pp_large = -1e10_wp 1477 1470 ! 1478 1471 IF( ln_timing ) CALL timing_start('zdf_osm_os') … … 1482 1475 l_shear(:,:) = .FALSE. 1483 1476 n_ddh(:,:) = 1 1484 pwb_ent(:,:) = zlarge1485 pwb_min(:,:) = zlarge1486 pshear(:,:) = zlarge1477 pwb_ent(:,:) = pp_large 1478 pwb_min(:,:) = pp_large 1479 pshear(:,:) = pp_large 1487 1480 ! 1488 1481 ! Determins stability and set flag l_conv … … 1496 1489 ! 1497 1490 pshear(A2D(0)) = 0._wp 1498 zekman(:,:) = EXP( -1.0_wp * zek * ABS( ff_t(A2D(0)) ) * phbl(A2D(0)) / MAX(sustar(A2D(0)), 1.e-8 ) )1491 zekman(:,:) = EXP( -1.0_wp * pp_ek * ABS( ff_t(A2D(0)) ) * phbl(A2D(0)) / MAX(sustar(A2D(0)), 1.e-8 ) ) 1499 1492 ! 1500 1493 #ifdef key_osm_debug … … 1520 1513 & MAX( pdv_ml(ji,jj), 1e-5_wp)**2 ) 1521 1514 END IF 1522 pshear(ji,jj) = za_shr * zekman(ji,jj) * &1515 pshear(ji,jj) = pp_a_shr * zekman(ji,jj) * & 1523 1516 & ( MAX( sustar(ji,jj)**2 * pdu_ml(ji,jj) / phbl(ji,jj), 0.0_wp ) + & 1524 & zb_shr * MAX( -1.0_wp * ff_t(ji,jj) * sustke(ji,jj) * dstokes(ji,jj) * &1517 & pp_b_shr * MAX( -1.0_wp * ff_t(ji,jj) * sustke(ji,jj) * dstokes(ji,jj) * & 1525 1518 & pdv_ml(ji,jj) / phbl(ji,jj), 0.0_wp ) ) 1526 1519 #ifdef key_osm_debug … … 1532 1525 #endif 1533 1526 ! Stability dependence 1534 pshear(ji,jj) = pshear(ji,jj) * EXP( -0.75_wp * MAX( 0.0_wp, ( zri_b(ji,jj) - zri_c ) / zri_c ) )1527 pshear(ji,jj) = pshear(ji,jj) * EXP( -0.75_wp * MAX( 0.0_wp, ( zri_b(ji,jj) - pp_ri_c ) / pp_ri_c ) ) 1535 1528 #ifdef key_osm_debug 1536 1529 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN … … 1544 1537 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1545 1538 IF ( pshear(ji,jj) > 1e-10 ) THEN 1546 IF ( zri_p(ji,jj) < rn_ri_p_thresh .AND. MIN( hu(ji,jj,Kmm), hu(ji-1,jj,Kmm), hv(ji,jj,Kmm), hv(ji,jj-1,Kmm) ) > 100.0_wp ) THEN1539 IF ( zri_p(ji,jj) < pp_ri_p_thresh .AND. MIN( hu(ji,jj,Kmm), hu(ji-1,jj,Kmm), hv(ji,jj,Kmm), hv(ji,jj-1,Kmm) ) > 100.0_wp ) THEN 1547 1540 ! Growing shear layer 1548 1541 n_ddh(ji,jj) = 0 … … 1585 1578 zr_stokes = 1.0 - EXP( -25.0_wp * dstokes(ji,jj) / hbl(ji,jj) * ( 1.0_wp + 4.0_wp * dstokes(ji,jj) / hbl(ji,jj) ) ) 1586 1579 END IF 1587 pwb_ent(ji,jj) = -2.0_wp * zalpha_c * zrf_conv * swbav(ji,jj) - &1588 & zalpha_s * zrf_shear * sustar(ji,jj)**3 / phml(ji,jj) + &1589 & zr_stokes * ( zalpha_s * EXP( -1.5_wp * sla(ji,jj) ) * zrf_shear * sustar(ji,jj)**3 - &1590 & zrf_langmuir * zalpha_lc * swstrl(ji,jj)**3 ) / phml(ji,jj)1580 pwb_ent(ji,jj) = -2.0_wp * pp_alpha_c * zrf_conv * swbav(ji,jj) - & 1581 & pp_alpha_s * zrf_shear * sustar(ji,jj)**3 / phml(ji,jj) + & 1582 & zr_stokes * ( pp_alpha_s * EXP( -1.5_wp * sla(ji,jj) ) * zrf_shear * sustar(ji,jj)**3 - & 1583 & zrf_langmuir * pp_alpha_lc * swstrl(ji,jj)**3 ) / phml(ji,jj) 1591 1584 #ifdef key_osm_debug 1592 1585 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN … … 1602 1595 IF ( l_conv(ji,jj) ) THEN 1603 1596 ! Unstable OSBL 1604 zwb_shr = -1.0_wp * za_wb_s * zri_b(ji,jj) * pshear(ji,jj)1597 zwb_shr = -1.0_wp * pp_a_wb_s * zri_b(ji,jj) * pshear(ji,jj) 1605 1598 #ifdef key_osm_debug 1606 1599 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN … … 1613 1606 1614 1607 ! pshear_u = MAX( zustar(ji,jj)**2 * MAX( pdu_ml(ji,jj), 0._wp ) / phbl(ji,jj), 0._wp ) 1615 ! pshear(ji,jj) = pshear(ji,jj) + pshear_u * ( 1.0 - MIN( zri_p(ji,jj) / rn_ri_p_thresh, 1.d0 )**2 )1608 ! pshear(ji,jj) = pshear(ji,jj) + pshear_u * ( 1.0 - MIN( zri_p(ji,jj) / pp_ri_p_thresh, 1.d0 )**2 ) 1616 1609 ! pshear(ji,jj) = MIN( pshear(ji,jj), pshear_u ) 1617 1610 1618 ! zwb_shr = zwb_shr - 0.25 * MAX ( pshear_u, 0._wp) * ( 1.0 - MIN( zri_p(ji,jj) / rn_ri_p_thresh, 1._wp )**2 )1611 ! zwb_shr = zwb_shr - 0.25 * MAX ( pshear_u, 0._wp) * ( 1.0 - MIN( zri_p(ji,jj) / pp_ri_p_thresh, 1._wp )**2 ) 1619 1612 ! zwb_shr = MAX( zwb_shr, -0.25 * pshear_u ) 1620 1613 #ifdef key_osm_debug … … 1666 1659 REAL(wp) :: zthermal, zbeta 1667 1660 ! 1668 REAL(wp), PARAMETER :: zlarge = -1e10_wp1661 REAL(wp), PARAMETER :: pp_large = -1e10_wp 1669 1662 ! 1670 1663 IF( ln_timing ) CALL timing_start('zdf_osm_eg') 1671 1664 ! 1672 pdtdz(:,:) = zlarge1673 pdsdz(:,:) = zlarge1674 pdbdz(:,:) = zlarge1665 pdtdz(:,:) = pp_large 1666 pdsdz(:,:) = pp_large 1667 pdbdz(:,:) = pp_large 1675 1668 ! 1676 1669 DO_2D( 0, 0, 0, 0 ) … … 1694 1687 END SUBROUTINE zdf_osm_external_gradients 1695 1688 1696 SUBROUTINE zdf_osm_calculate_dhdt( pdhdt, phbl, pdh, pwb_ent, pwb_min,&1697 & pdb_bl, pdbdz_bl_ext, pdb_ml, pwb_fk_b, pwb_fk, 1689 SUBROUTINE zdf_osm_calculate_dhdt( pdhdt, phbl, pdh, pwb_ent, pwb_min, & 1690 & pdb_bl, pdbdz_bl_ext, pdb_ml, pwb_fk_b, pwb_fk, & 1698 1691 & pvel_mle ) 1699 1692 !!--------------------------------------------------------------------- … … 1722 1715 REAL(wp) :: zvel_max, zddhdt 1723 1716 ! 1724 REAL(wp), PARAMETER :: zzeta_m = 0.3_wp 1725 REAL(wp), PARAMETER :: zgamma_c = 2.0_wp 1726 REAL(wp), PARAMETER :: zdhoh = 0.1_wp 1727 REAL(wp), PARAMETER :: zalpha_b = 0.3_wp 1728 REAL(wp), PARAMETER :: a_ddh = 2.5_wp, a_ddh_2 = 3.5_wp ! Also in pycnocline_depth 1729 REAL(wp), PARAMETER :: zlarge = -1e10_wp 1717 REAL(wp), PARAMETER :: pp_alpha_b = 0.3_wp 1718 REAL(wp), PARAMETER :: pp_ddh = 2.5_wp, pp_ddh_2 = 3.5_wp ! Also in pycnocline_depth 1719 REAL(wp), PARAMETER :: pp_large = -1e10_wp 1730 1720 ! 1731 1721 IF( ln_timing ) CALL timing_start('zdf_osm_cd') 1732 1722 ! 1733 pdhdt(:,:) = zlarge1734 pwb_fk_b(:,:) = zlarge1723 pdhdt(:,:) = pp_large 1724 pwb_fk_b(:,:) = pp_large 1735 1725 ! 1736 1726 DO_2D( 0, 0, 0, 0 ) … … 1756 1746 zpsi = zpsi + 1.75_wp * ( 1.0_wp - 0.5_wp * pdh(ji,jj) / phbl(ji,jj) ) * & 1757 1747 & ( pdh(ji,jj) / phbl(ji,jj) + zgamma_b_nd ) * MIN( ( pwb_min(ji,jj) + 2.0_wp * pwb_fk_b(ji,jj) ), 0.0_wp ) 1758 zpsi = zalpha_b * MAX( zpsi, 0.0_wp )1748 zpsi = pp_alpha_b * MAX( zpsi, 0.0_wp ) 1759 1749 pdhdt(ji,jj) = -1.0_wp * ( pwb_ent(ji,jj) + 2.0_wp * pwb_fk_b(ji,jj) ) / & 1760 1750 & ( zvel_max + MAX( pdb_bl(ji,jj), 1e-15_wp ) ) + & … … 1780 1770 ENDIF 1781 1771 ! Relaxation to dh_ref = zari * hbl 1782 zddhdt = -1.0_wp * a_ddh_2 * ( 1.0_wp - pdh(ji,jj) / ( zari * phbl(ji,jj) ) ) * pwb_ent(ji,jj) / &1772 zddhdt = -1.0_wp * pp_ddh_2 * ( 1.0_wp - pdh(ji,jj) / ( zari * phbl(ji,jj) ) ) * pwb_ent(ji,jj) / & 1783 1773 & ( zvel_max + MAX( pdb_bl(ji,jj), 1e-15_wp ) ) 1784 1774 #ifdef key_osm_debug … … 1789 1779 #endif 1790 1780 ELSE IF ( n_ddh(ji,jj) == 0 ) THEN ! Growing shear layer 1791 zddhdt = -1.0_wp * a_ddh * ( 1.0_wp - 1.6_wp * pdh(ji,jj) / phbl(ji,jj) ) * pwb_ent(ji,jj) / &1781 zddhdt = -1.0_wp * pp_ddh * ( 1.0_wp - 1.6_wp * pdh(ji,jj) / phbl(ji,jj) ) * pwb_ent(ji,jj) / & 1792 1782 & ( zvel_max + MAX( pdb_bl(ji,jj), 1e-15_wp ) ) 1793 1783 zddhdt = EXP( -4.0_wp * ABS( ff_t(ji,jj) ) * phbl(ji,jj) / MAX(sustar(ji,jj), 1e-8_wp ) ) * zddhdt … … 1795 1785 zddhdt = 0.0_wp 1796 1786 ENDIF ! n_ddh 1797 pdhdt(ji,jj) = pdhdt(ji,jj) + zalpha_b * ( 1.0_wp - 0.5_wp * pdh(ji,jj) / phbl(ji,jj) ) * &1787 pdhdt(ji,jj) = pdhdt(ji,jj) + pp_alpha_b * ( 1.0_wp - 0.5_wp * pdh(ji,jj) / phbl(ji,jj) ) * & 1798 1788 & pdb_ml(ji,jj) * MAX( zddhdt, 0.0_wp ) / & 1799 1789 & ( zvel_max + MAX( pdb_bl(ji,jj), 1e-15_wp ) ) … … 1889 1879 END SUBROUTINE zdf_osm_calculate_dhdt 1890 1880 1891 SUBROUTINE zdf_osm_timestep_hbl( Kmm, pdhdt, phbl, phbl_t, pt_bl, &1892 & p s_bl, pwb_ent, pwb_fk_b )1881 SUBROUTINE zdf_osm_timestep_hbl( Kmm, pdhdt, phbl, phbl_t, pwb_ent, & 1882 & pwb_fk_b ) 1893 1883 !!--------------------------------------------------------------------- 1894 1884 !! *** ROUTINE zdf_osm_timestep_hbl *** … … 1906 1896 REAL(wp), DIMENSION(:,:), INTENT(inout) :: phbl ! BL depth 1907 1897 REAL(wp), DIMENSION(:,:), INTENT(in ) :: phbl_t ! BL depth 1908 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pt_bl ! Temperature average over the depth of the blayer1909 REAL(wp), DIMENSION(:,:), INTENT(in ) :: ps_bl ! Salinity average over the depth of the blayer1910 1898 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pwb_ent ! Buoyancy entrainment flux 1911 1899 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pwb_fk_b ! MLE buoyancy flux averaged over OSBL … … 1950 1938 #endif 1951 1939 DO jk = nmld(ji,jj), nbld(ji,jj) 1952 zdb = MAX( grav * ( zthermal * ( pt_bl(ji,jj) - ts(ji,jj,jm,jp_tem,Kmm) ) - &1953 & zbeta * ( ps_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ), 0.0_wp ) + zvel_max1940 zdb = MAX( grav * ( zthermal * ( av_t_bl(ji,jj) - ts(ji,jj,jm,jp_tem,Kmm) ) - & 1941 & zbeta * ( av_s_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ), 0.0_wp ) + zvel_max 1954 1942 ! 1955 1943 IF ( ln_osm_mle ) THEN … … 1984 1972 #endif 1985 1973 DO jk = nmld(ji,jj), nbld(ji,jj) 1986 zdb = MAX( grav * ( zthermal * ( pt_bl(ji,jj) - ts(ji,jj,jm,jp_tem,Kmm) ) - &1987 & zbeta * ( ps_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ), 0.0_wp ) + &1974 zdb = MAX( grav * ( zthermal * ( av_t_bl(ji,jj) - ts(ji,jj,jm,jp_tem,Kmm) ) - & 1975 & zbeta * ( av_s_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ), 0.0_wp ) + & 1988 1976 & 2.0 * svstr(ji,jj)**2 / zhbl_s 1989 1977 ! … … 2032 2020 END SUBROUTINE zdf_osm_timestep_hbl 2033 2021 2034 SUBROUTINE zdf_osm_pycnocline_thickness( Kmm, pdh, phml,pdhdt, pdb_bl, &2022 SUBROUTINE zdf_osm_pycnocline_thickness( Kmm, pdh, phml, pdhdt, pdb_bl, & 2035 2023 & phbl, pwb_ent, pdbdz_bl_ext, pwb_fk_b ) 2036 2024 !!--------------------------------------------------------------------- … … 2063 2051 REAL(wp) :: ztmp ! Auxiliary variable 2064 2052 ! 2065 REAL, PARAMETER :: a_ddh = 2.5_wp, a_ddh_2 = 3.5_wp ! Also in pycnocline_depth2053 REAL, PARAMETER :: pp_ddh = 2.5_wp, pp_ddh_2 = 3.5_wp ! Also in pycnocline_depth 2066 2054 ! 2067 2055 IF( ln_timing ) CALL timing_start('zdf_osm_pt') … … 2077 2065 zvel_max = ( svstr(ji,jj)**3 + 0.5_wp * swstrc(ji,jj)**3 )**p2third / hbl(ji,jj) 2078 2066 ! ddhdt for pycnocline determined in osm_calculate_dhdt 2079 zddhdt = -1.0_wp * a_ddh * ( 1.0_wp - 1.6_wp * pdh(ji,jj) / phbl(ji,jj) ) * pwb_ent(ji,jj) / &2067 zddhdt = -1.0_wp * pp_ddh * ( 1.0_wp - 1.6_wp * pdh(ji,jj) / phbl(ji,jj) ) * pwb_ent(ji,jj) / & 2080 2068 & ( zvel_max + MAX( pdb_bl(ji,jj), 1e-15 ) ) 2081 2069 zddhdt = EXP( -4.0_wp * ABS( ff_t(ji,jj) ) * phbl(ji,jj) / MAX( sustar(ji,jj), 1e-8 ) ) * zddhdt … … 2091 2079 & pdb_bl(ji,jj)**2 / MAX( 4.5_wp * ztmp**2, & 2092 2080 & 1e-12_wp ) ) ), 0.2_wp ) 2093 ztau = MAX( pdb_bl(ji,jj) * ( zari * hbl(ji,jj) ) / ( a_ddh_2 * MAX( -1.0_wp * pwb_ent(ji,jj), 1e-12_wp ) ), &2081 ztau = MAX( pdb_bl(ji,jj) * ( zari * hbl(ji,jj) ) / ( pp_ddh_2 * MAX( -1.0_wp * pwb_ent(ji,jj), 1e-12_wp ) ), & 2094 2082 & 2.0_wp * rn_Dt ) 2095 2083 dh(ji,jj) = dh(ji,jj) * EXP( -1.0_wp * rn_Dt / ztau ) + & … … 2210 2198 END SUBROUTINE zdf_osm_pycnocline_thickness 2211 2199 2212 SUBROUTINE zdf_osm_pycnocline_buoyancy_profiles( Kmm, kp_ext, pdbdz, palpha, pdh,&2213 & pdb_bl, phbl, pdbdz_bl_ext, phml,pdb_ml, &2200 SUBROUTINE zdf_osm_pycnocline_buoyancy_profiles( Kmm, kp_ext, pdbdz, palpha, pdh, & 2201 & pdb_bl, phbl, pdbdz_bl_ext, phml, pdb_ml, & 2214 2202 & pdhdt ) 2215 2203 !!--------------------------------------------------------------------- … … 2240 2228 REAL(wp) :: ztmp ! Auxiliary variable 2241 2229 ! 2242 REAL(wp), PARAMETER :: pp gamma_b = 2.25_wp2230 REAL(wp), PARAMETER :: pp_gamma_b = 2.25_wp 2243 2231 ! 2244 2232 IF( ln_timing ) CALL timing_start('zdf_osm_pscp') … … 2253 2241 ! 2254 2242 zzeta_m = 0.1_wp + 0.3_wp / ( 1.0_wp + EXP( -3.5_wp * LOG10( -1.0_wp * shol(ji,jj) ) ) ) 2255 palpha(ji,jj) = 2.0_wp * ( 1.0_wp - ( 0.80_wp * zzeta_m + 0.5_wp * SQRT( 3.14159_wp / pp gamma_b ) ) * &2243 palpha(ji,jj) = 2.0_wp * ( 1.0_wp - ( 0.80_wp * zzeta_m + 0.5_wp * SQRT( 3.14159_wp / pp_gamma_b ) ) * & 2256 2244 & pdbdz_bl_ext(ji,jj) * pdh(ji,jj) / pdb_ml(ji,jj) ) / & 2257 & ( 0.723_wp + SQRT( 3.14159_wp / pp gamma_b ) )2245 & ( 0.723_wp + SQRT( 3.14159_wp / pp_gamma_b ) ) 2258 2246 palpha(ji,jj) = MAX( palpha(ji,jj), 0.0_wp ) 2259 2247 ztmp = 1.0_wp / MAX( pdh(ji,jj), epsln ) … … 2276 2264 & EXP( -6.0_wp * ( znd -zzeta_m )**2 ) 2277 2265 ELSE 2278 ! zdtdz(ji,jj,jk) = ztgrad * EXP( - zgamma_b * ( znd - zzeta_m )**2 )2279 ! zdsdz(ji,jj,jk) = zsgrad * EXP( - zgamma_b * ( znd - zzeta_m )**2 )2280 pdbdz(ji,jj,jk) = zbgrad * EXP( -1.0_wp * pp gamma_b * ( znd - zzeta_m )**2 )2266 ! zdtdz(ji,jj,jk) = ztgrad * EXP( -pp_gamma_b * ( znd - zzeta_m )**2 ) 2267 ! zdsdz(ji,jj,jk) = zsgrad * EXP( -pp_gamma_b * ( znd - zzeta_m )**2 ) 2268 pdbdz(ji,jj,jk) = zbgrad * EXP( -1.0_wp * pp_gamma_b * ( znd - zzeta_m )**2 ) 2281 2269 END IF 2282 2270 END DO … … 2335 2323 END SUBROUTINE zdf_osm_pycnocline_buoyancy_profiles 2336 2324 2337 SUBROUTINE zdf_osm_diffusivity_viscosity( Kbb, Kmm, pdiffut, pviscos, phbl, & 2338 & phml, pdh, pdhdt, pshear, pb_bl, & 2339 & pu_bl, pv_bl, pb_ml, pu_ml, pv_ml, & 2325 SUBROUTINE zdf_osm_diffusivity_viscosity( Kbb, Kmm, pdiffut, pviscos, phbl, & 2326 & phml, pdh, pdhdt, pshear, & 2340 2327 & pwb_ent, pwb_min ) 2341 2328 !!--------------------------------------------------------------------- … … 2356 2343 REAL(wp), DIMENSION(A2D(0)), INTENT(in ) :: pdhdt ! BL depth tendency 2357 2344 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pshear ! Shear production 2358 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pb_bl ! Buoyancy average over the depth of the boundary layer2359 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pu_bl ! Velocity average over the depth of the boundary layer2360 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pv_bl ! Velocity average over the depth of the boundary layer2361 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pb_ml ! Buoyancy average over the depth of the mixed layer2362 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pu_ml ! Velocity average over the depth of the mixed layer2363 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pv_ml ! Velocity average over the depth of the mixed layer2364 2345 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pwb_ent ! Buoyancy entrainment flux 2365 2346 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pwb_min … … 2381 2362 REAL(wp) :: zmsku, zmskv 2382 2363 ! 2383 REAL(wp), PARAMETER :: rn_dif_ml = 0.8_wp, rn_vis_ml = 0.375_wp2384 REAL(wp), PARAMETER :: rn_dif_pyc = 0.15_wp, rn_vis_pyc = 0.142_wp2385 REAL(wp), PARAMETER :: rn_vispyc_shr = 0.15_wp2364 REAL(wp), PARAMETER :: pp_dif_ml = 0.8_wp, pp_vis_ml = 0.375_wp 2365 REAL(wp), PARAMETER :: pp_dif_pyc = 0.15_wp, pp_vis_pyc = 0.142_wp 2366 REAL(wp), PARAMETER :: pp_vispyc_shr = 0.15_wp 2386 2367 ! 2387 2368 IF( ln_timing ) CALL timing_start('zdf_osm_dv') … … 2397 2378 & ( 1.4_wp - 0.4_wp / ( 1.0_wp + EXP(-3.5_wp * LOG10(-shol(ji,jj) ) ) )**1.25_wp ) )**2 2398 2379 ! 2399 zdifml_sc(ji,jj) = rn_dif_ml * phml(ji,jj) * zvel_sc_ml2400 zvisml_sc(ji,jj) = rn_vis_ml * zdifml_sc(ji,jj)2380 zdifml_sc(ji,jj) = pp_dif_ml * phml(ji,jj) * zvel_sc_ml 2381 zvisml_sc(ji,jj) = pp_vis_ml * zdifml_sc(ji,jj) 2401 2382 #ifdef key_osm_debug 2402 2383 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN … … 2408 2389 ! 2409 2390 IF ( l_pyc(ji,jj) ) THEN 2410 zdifpyc_n_sc(ji,jj) = rn_dif_pyc * zvel_sc_ml * pdh(ji,jj) 2411 zvispyc_n_sc(ji,jj) = 0.09_wp * zvel_sc_pyc * ( 1.0_wp - phbl(ji,jj) / pdh(ji,jj) )**2 * & 2412 & ( 0.005 * ( pu_ml(ji,jj)-pu_bl(ji,jj) )**2 + 0.0075 * ( pv_ml(ji,jj)-pv_bl(ji,jj) )**2 ) / & 2391 zdifpyc_n_sc(ji,jj) = pp_dif_pyc * zvel_sc_ml * pdh(ji,jj) 2392 zvispyc_n_sc(ji,jj) = 0.09_wp * zvel_sc_pyc * ( 1.0_wp - phbl(ji,jj) / pdh(ji,jj) )**2 * & 2393 & ( 0.005_wp * ( av_u_ml(ji,jj) - av_u_bl(ji,jj) )**2 + & 2394 & 0.0075_wp * ( av_v_ml(ji,jj) - av_v_bl(ji,jj) )**2 ) / & 2413 2395 & pdh(ji,jj) 2414 zvispyc_n_sc(ji,jj) = rn_vis_pyc * zvel_sc_ml * pdh(ji,jj) + zvispyc_n_sc(ji,jj) * zstab_fac2396 zvispyc_n_sc(ji,jj) = pp_vis_pyc * zvel_sc_ml * pdh(ji,jj) + zvispyc_n_sc(ji,jj) * zstab_fac 2415 2397 #ifdef key_osm_debug 2416 2398 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN … … 2421 2403 ! 2422 2404 IF ( l_shear(ji,jj) .AND. n_ddh(ji,jj) /= 2 ) THEN 2423 ztmp = rn_vispyc_shr * ( pshear(ji,jj) * phbl(ji,jj) )**pthird * phbl(ji,jj)2405 ztmp = pp_vispyc_shr * ( pshear(ji,jj) * phbl(ji,jj) )**pthird * phbl(ji,jj) 2424 2406 zdifpyc_n_sc(ji,jj) = zdifpyc_n_sc(ji,jj) + ztmp 2425 2407 zvispyc_n_sc(ji,jj) = zvispyc_n_sc(ji,jj) + ztmp … … 2433 2415 ! 2434 2416 zdifpyc_s_sc(ji,jj) = pwb_ent(ji,jj) + 0.0025_wp * zvel_sc_pyc * ( phbl(ji,jj) / pdh(ji,jj) - 1.0_wp ) * & 2435 & ( pb_ml(ji,jj) - pb_bl(ji,jj) )2417 & ( av_b_ml(ji,jj) - av_b_bl(ji,jj) ) 2436 2418 zvispyc_s_sc(ji,jj) = 0.09_wp * ( pwb_min(ji,jj) + 0.0025_wp * zvel_sc_pyc * & 2437 2419 & ( phbl(ji,jj) / pdh(ji,jj) - 1.0_wp ) * & 2438 & ( pb_ml(ji,jj) - pb_bl(ji,jj) ) )2420 & ( av_b_ml(ji,jj) - av_b_bl(ji,jj) ) ) 2439 2421 #ifdef key_osm_debug 2440 2422 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN … … 2459 2441 zbeta_v_sc(ji,jj) = 1.0_wp - 2.0_wp * ( zvispyc_n_sc(ji,jj) + zvispyc_s_sc(ji,jj) ) / ( zvisml_sc(ji,jj) + epsln ) 2460 2442 ELSE 2461 zdifpyc_n_sc(ji,jj) = rn_dif_pyc * zvel_sc_ml * pdh(ji,jj) ! ag 19/032443 zdifpyc_n_sc(ji,jj) = pp_dif_pyc * zvel_sc_ml * pdh(ji,jj) ! ag 19/03 2462 2444 zdifpyc_s_sc(ji,jj) = 0.0_wp ! ag 19/03 2463 zvispyc_n_sc(ji,jj) = rn_vis_pyc * zvel_sc_ml * pdh(ji,jj) ! ag 19/032445 zvispyc_n_sc(ji,jj) = pp_vis_pyc * zvel_sc_ml * pdh(ji,jj) ! ag 19/03 2464 2446 zvispyc_s_sc(ji,jj) = 0.0_wp ! ag 19/03 2465 2447 IF(l_coup(ji,jj) ) THEN ! ag 19/03 … … 2600 2582 END SUBROUTINE zdf_osm_diffusivity_viscosity 2601 2583 2602 SUBROUTINE zdf_osm_fgr_terms( Kmm, kp_ext, phbl, phml, pdh,&2603 & pdhdt, pshear, pdt_bl, pds_bl, pdb_bl,&2604 & pdu_bl, pdv_bl, pdt_ml, pds_ml, pdb_ml,&2605 & pdu_ml, pdv_ml,pdtdz_bl_ext, pdsdz_bl_ext, pdbdz_pyc, &2584 SUBROUTINE zdf_osm_fgr_terms( Kmm, kp_ext, phbl, phml, pdh, & 2585 & pdhdt, pshear, pdt_bl, pds_bl, pdb_bl, & 2586 & pdu_bl, pdv_bl, pdt_ml, pds_ml, pdb_ml, & 2587 & pdu_ml, pdv_ml, pdtdz_bl_ext, pdsdz_bl_ext, pdbdz_pyc, & 2606 2588 & palpha_pyc, pdiffut, pviscos ) 2607 2589 !!--------------------------------------------------------------------- … … 3225 3207 END SUBROUTINE zdf_osm_fgr_terms 3226 3208 3227 SUBROUTINE zdf_osm_zmld_horizontal_gradients( Kmm, pmld, pdtdx,pdtdy, pdsdx, &3209 SUBROUTINE zdf_osm_zmld_horizontal_gradients( Kmm, pmld, pdtdx, pdtdy, pdsdx, & 3228 3210 & pdsdy, pdbdx_mle, pdbdy_mle, pdbds_mle ) 3229 3211 !!---------------------------------------------------------------------- … … 3326 3308 END SUBROUTINE zdf_osm_zmld_horizontal_gradients 3327 3309 3328 SUBROUTINE zdf_osm_osbl_state_fk( Kmm, pwb_fk, pt_mle, ps_mle, pb_mle, & 3329 & phbl, phmle, pwb_ent, pt_bl, ps_bl, & 3330 & pb_bl, pdb_bl, pdbds_mle ) 3310 SUBROUTINE zdf_osm_osbl_state_fk( Kmm, pwb_fk, phbl, phmle, pwb_ent, & 3311 & pdb_bl, pdbds_mle ) 3331 3312 !!--------------------------------------------------------------------- 3332 3313 !! *** ROUTINE zdf_osm_osbl_state_fk *** … … 3349 3330 INTEGER, INTENT(in ) :: Kmm ! Time-level index 3350 3331 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pwb_fk 3351 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pt_mle ! Temperature average over the depth of the MLE layer3352 REAL(wp), DIMENSION(:,:), INTENT(inout) :: ps_mle ! Salinity average over the depth of the MLE layer3353 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pb_mle ! Buoyancy average over the depth of the MLE layer3354 3332 REAL(wp), DIMENSION(:,:), INTENT(in ) :: phbl ! BL depth 3355 3333 REAL(wp), DIMENSION(:,:), INTENT(in ) :: phmle ! MLE depth 3356 3334 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pwb_ent ! Buoyancy entrainment flux 3357 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pt_bl ! Temperature average over the depth of the blayer3358 REAL(wp), DIMENSION(:,:), INTENT(in ) :: ps_bl ! Salinity average over the depth of the blayer3359 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pb_bl ! Buoyancy average over the depth of the blayer3360 3335 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdb_bl ! Difference between blayer average and parameter at base of blayer 3361 3336 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdbds_mle ! Magnitude of horizontal buoyancy gradient … … 3384 3359 IF ( l_conv(ji,jj) ) THEN 3385 3360 IF ( phmle(ji,jj) > 1.2_wp * phbl(ji,jj) ) THEN 3386 pt_mle(ji,jj) = ( pt_mle(ji,jj) * phmle(ji,jj) - pt_bl(ji,jj) * phbl(ji,jj) ) / ( phmle(ji,jj) - phbl(ji,jj) )3387 ps_mle(ji,jj) = ( ps_mle(ji,jj) * phmle(ji,jj) - ps_bl(ji,jj) * phbl(ji,jj) ) / ( phmle(ji,jj) - phbl(ji,jj) )3388 pb_mle(ji,jj) = ( pb_mle(ji,jj) * phmle(ji,jj) - pb_bl(ji,jj) * phbl(ji,jj) ) / ( phmle(ji,jj) - phbl(ji,jj) )3389 zdbdz_mle_int = ( pb_bl(ji,jj) - ( 2.0_wp * pb_mle(ji,jj) - pb_bl(ji,jj) ) ) / ( phmle(ji,jj) - phbl(ji,jj) )3361 av_t_mle(ji,jj) = ( av_t_mle(ji,jj) * phmle(ji,jj) - av_t_bl(ji,jj) * phbl(ji,jj) ) / ( phmle(ji,jj) - phbl(ji,jj) ) 3362 av_s_mle(ji,jj) = ( av_s_mle(ji,jj) * phmle(ji,jj) - av_s_bl(ji,jj) * phbl(ji,jj) ) / ( phmle(ji,jj) - phbl(ji,jj) ) 3363 av_b_mle(ji,jj) = ( av_b_mle(ji,jj) * phmle(ji,jj) - av_b_bl(ji,jj) * phbl(ji,jj) ) / ( phmle(ji,jj) - phbl(ji,jj) ) 3364 zdbdz_mle_int = ( av_b_bl(ji,jj) - ( 2.0_wp * av_b_mle(ji,jj) - av_b_bl(ji,jj) ) ) / ( phmle(ji,jj) - phbl(ji,jj) ) 3390 3365 ! Calculate potential energies of actual profile and reference profile 3391 3366 zpe_mle_layer = 0.0_wp … … 3396 3371 zbuoy = grav * ( zthermal * ts(ji,jj,jk,jp_tem,Kmm) - zbeta * ts(ji,jj,jk,jp_sal,Kmm) ) 3397 3372 zpe_mle_layer = zpe_mle_layer + zbuoy * gdepw(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) 3398 zpe_mle_ref = zpe_mle_ref + ( pb_bl(ji,jj) - zdbdz_mle_int * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) ) * &3373 zpe_mle_ref = zpe_mle_ref + ( av_b_bl(ji,jj) - zdbdz_mle_int * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) ) * & 3399 3374 & gdepw(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) 3400 3375 END DO … … 3470 3445 END SUBROUTINE zdf_osm_osbl_state_fk 3471 3446 3472 SUBROUTINE zdf_osm_mle_parameters( Kmm, 3473 & pdiff_mle, pdbds_mle, phbl, p b_bl, pwb0tot )3447 SUBROUTINE zdf_osm_mle_parameters( Kmm, kmld_prof, pmld, phmle, pvel_mle, & 3448 & pdiff_mle, pdbds_mle, phbl, pwb0tot ) 3474 3449 !!---------------------------------------------------------------------- 3475 3450 !! *** ROUTINE zdf_osm_mle_parameters *** … … 3493 3468 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdbds_mle ! Magnitude of horizontal buoyancy gradient 3494 3469 REAL(wp), DIMENSION(:,:), INTENT(in ) :: phbl ! BL depth 3495 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pb_bl ! Buoyancy average over the depth of the blayer3496 3470 REAL(wp), DIMENSION(A2D(0)), INTENT(in ) :: pwb0tot ! Total surface buoyancy flux including insolation 3497 3471 ! … … 3526 3500 zbuoy = grav * ( zthermal * ts(ji,jj,kmld_prof(ji,jj)+2,jp_tem,Kmm) - & 3527 3501 & zbeta * ts(ji,jj,kmld_prof(ji,jj)+2,jp_sal,Kmm) ) 3528 zdb_mle = pb_bl(ji,jj) - zbuoy3502 zdb_mle = av_b_bl(ji,jj) - zbuoy 3529 3503 ! Timestep hmle 3530 3504 hmle(ji,jj) = hmle(ji,jj) + pwb0tot(ji,jj) * rn_Dt / zdb_mle … … 3563 3537 !! called at the first timestep (nit000) 3564 3538 !! 3565 !! ** input : Namlist namosm 3539 !! ** input : Namlists namzdf_osm and namosm_mle 3540 !! 3566 3541 !!---------------------------------------------------------------------- 3567 INTEGER, INTENT(in) :: Kmm ! time level 3568 INTEGER :: ios ! local integer 3569 INTEGER :: ji, jj, jk ! dummy loop indices 3570 REAL z1_t2 3571 REAL(wp) :: zlarge = -1.0e10_wp, zero = 0.0_wp 3572 !! 3573 NAMELIST/namzdf_osm/ ln_use_osm_la, rn_osm_la, rn_osm_dstokes, nn_ave & 3574 & ,nn_osm_wave, ln_dia_osm, rn_osm_hbl0, rn_zdfosm_adjust_sd & 3575 & ,ln_kpprimix, rn_riinfty, rn_difri, ln_convmix, rn_difconv, nn_osm_wave & 3576 #ifdef key_osm_debug 3577 & ,nn_osm_SD_reduce, ln_osm_mle, rn_osm_hblfrac, rn_osm_bl_thresh, ln_zdfosm_ice_shelter & 3578 & ,nn_idb, nn_jdb, nn_kdb, nn_narea_db 3579 #else 3580 & ,nn_osm_SD_reduce, ln_osm_mle, rn_osm_hblfrac, rn_osm_bl_thresh, ln_zdfosm_ice_shelter 3581 #endif 3582 ! Namelist for Fox-Kemper parametrization. 3583 NAMELIST/namosm_mle/ nn_osm_mle, rn_osm_mle_ce, rn_osm_mle_lf, rn_osm_mle_time, rn_osm_mle_lat,& 3584 & rn_osm_mle_rho_c, rn_osm_mle_thresh, rn_osm_mle_tau, ln_osm_hmle_limit, rn_osm_hmle_limit 3585 3586 !!---------------------------------------------------------------------- 3542 INTEGER, INTENT(in ) :: Kmm ! Time level 3543 ! 3544 ! Local variables 3545 INTEGER :: ios ! Local integer 3546 INTEGER :: ji, jj, jk ! Dummy loop indices 3547 REAL(wp) :: z1_t2 3548 ! 3549 REAL(wp), PARAMETER :: pp_large = -1e10_wp 3550 ! 3551 NAMELIST/namzdf_osm/ ln_use_osm_la, rn_osm_la, rn_osm_dstokes, nn_ave, nn_osm_wave, & 3552 & ln_dia_osm, rn_osm_hbl0, rn_zdfosm_adjust_sd, ln_kpprimix, rn_riinfty, & 3553 & rn_difri, ln_convmix, rn_difconv, nn_osm_wave, nn_osm_SD_reduce, & 3554 & ln_osm_mle, rn_osm_hblfrac, rn_osm_bl_thresh, ln_zdfosm_ice_shelter 3555 #ifdef key_osm_debug 3556 NAMELIST/namzdf_osm/ nn_idb, nn_jdb, nn_kdb, nn_narea_db 3557 #endif 3558 ! Namelist for Fox-Kemper parametrization 3559 NAMELIST/namosm_mle/ nn_osm_mle, rn_osm_mle_ce, rn_osm_mle_lf, rn_osm_mle_time, rn_osm_mle_lat, & 3560 & rn_osm_mle_rho_c, rn_osm_mle_thresh, rn_osm_mle_tau, ln_osm_hmle_limit, rn_osm_hmle_limit 3587 3561 ! 3588 3562 IF( ln_timing ) CALL timing_start('zdf_osm_init') … … 3599 3573 WRITE(numout,*) '~~~~~~~~~~~~' 3600 3574 WRITE(numout,*) ' Namelist namzdf_osm : set osm mixing parameters' 3601 WRITE(numout,*) ' Use rn_osm_la ln_use_osm_la= ', ln_use_osm_la3602 WRITE(numout,*) ' Use MLE in OBL, i.e. Fox-Kemper param ln_osm_mle= ', ln_osm_mle3603 WRITE(numout,*) ' Turbulent Langmuir number rn_osm_la= ', rn_osm_la3604 WRITE(numout,*) ' Stokes drift reduction factorrn_zdfosm_adjust_sd = ', rn_zdfosm_adjust_sd3605 WRITE(numout,*) ' Initial hbl for 1D runs rn_osm_hbl0= ', rn_osm_hbl03606 WRITE(numout,*) ' Depth scale of Stokes drift rn_osm_dstokes= ', rn_osm_dstokes3607 WRITE(numout,*) ' horizontal average flag nn_ave= ', nn_ave3608 WRITE(numout,*) ' Stokes drift nn_osm_wave= ', nn_osm_wave3575 WRITE(numout,*) ' Use rn_osm_la ln_use_osm_la = ', ln_use_osm_la 3576 WRITE(numout,*) ' Use MLE in OBL, i.e. Fox-Kemper param ln_osm_mle = ', ln_osm_mle 3577 WRITE(numout,*) ' Turbulent Langmuir number rn_osm_la = ', rn_osm_la 3578 WRITE(numout,*) ' Stokes drift reduction factor rn_zdfosm_adjust_sd = ', rn_zdfosm_adjust_sd 3579 WRITE(numout,*) ' Initial hbl for 1D runs rn_osm_hbl0 = ', rn_osm_hbl0 3580 WRITE(numout,*) ' Depth scale of Stokes drift rn_osm_dstokes = ', rn_osm_dstokes 3581 WRITE(numout,*) ' Horizontal average flag nn_ave = ', nn_ave 3582 WRITE(numout,*) ' Stokes drift nn_osm_wave = ', nn_osm_wave 3609 3583 SELECT CASE (nn_osm_wave) 3610 3584 CASE(0) 3611 WRITE(numout,*) ' calculated assuming constant La#=0.3'3585 WRITE(numout,*) ' Calculated assuming constant La#=0.3' 3612 3586 CASE(1) 3613 WRITE(numout,*) ' calculated from Pierson Moskowitz wind-waves'3587 WRITE(numout,*) ' Calculated from Pierson Moskowitz wind-waves' 3614 3588 CASE(2) 3615 WRITE(numout,*) ' calculated from ECMWF wave fields'3589 WRITE(numout,*) ' Calculated from ECMWF wave fields' 3616 3590 END SELECT 3617 WRITE(numout,*) ' Stokes drift reduction nn_osm_SD_reduce', nn_osm_SD_reduce3618 WRITE(numout,*) ' fraction of hbl to average SD over/fit'3619 WRITE(numout,*) ' exponential with nn_osm_SD_reduce = 1 or 2 rn_osm_hblfrac =', rn_osm_hblfrac3591 WRITE(numout,*) ' Stokes drift reduction nn_osm_SD_reduce = ', nn_osm_SD_reduce 3592 WRITE(numout,*) ' Fraction of hbl to average SD over/fit' 3593 WRITE(numout,*) ' Exponential with nn_osm_SD_reduce = 1 or 2 rn_osm_hblfrac = ', rn_osm_hblfrac 3620 3594 SELECT CASE (nn_osm_SD_reduce) 3621 3595 CASE(0) … … 3626 3600 WRITE(numout,*) ' Fit exponential to slope rn_osm_hblfrac of BL' 3627 3601 END SELECT 3628 WRITE(numout,*) ' reduce surface SD and depth scale under ice ln_zdfosm_ice_shelter=', ln_zdfosm_ice_shelter 3629 WRITE(numout,*) ' Output osm diagnostics ln_dia_osm = ', ln_dia_osm 3630 WRITE(numout,*) ' Threshold used to define BL rn_osm_bl_thresh = ', rn_osm_bl_thresh, 'm^2/s' 3631 WRITE(numout,*) ' Use KPP-style shear instability mixing ln_kpprimix = ', ln_kpprimix 3632 WRITE(numout,*) ' local Richardson Number limit for shear instability rn_riinfty = ', rn_riinfty 3633 WRITE(numout,*) ' maximum shear diffusivity at Rig = 0 (m2/s) rn_difri = ', rn_difri 3634 WRITE(numout,*) ' Use large mixing below BL when unstable ln_convmix = ', ln_convmix 3635 WRITE(numout,*) ' diffusivity when unstable below BL (m2/s) rn_difconv = ', rn_difconv 3602 WRITE(numout,*) ' Reduce surface SD and depth scale under ice ln_zdfosm_ice_shelter = ', ln_zdfosm_ice_shelter 3603 WRITE(numout,*) ' Output osm diagnostics ln_dia_osm = ', ln_dia_osm 3604 WRITE(numout,*) ' Threshold used to define BL rn_osm_bl_thresh = ', rn_osm_bl_thresh, & 3605 & 'm^2/s' 3606 WRITE(numout,*) ' Use KPP-style shear instability mixing ln_kpprimix = ', ln_kpprimix 3607 WRITE(numout,*) ' Local Richardson Number limit for shear instability rn_riinfty = ', rn_riinfty 3608 WRITE(numout,*) ' Maximum shear diffusivity at Rig = 0 (m2/s) rn_difri = ', rn_difri 3609 WRITE(numout,*) ' Use large mixing below BL when unstable ln_convmix = ', ln_convmix 3610 WRITE(numout,*) ' Diffusivity when unstable below BL (m2/s) rn_difconv = ', rn_difconv 3636 3611 #ifdef key_osm_debug 3637 3612 WRITE(numout,*) 'nn_idb', nn_idb, 'nn_jdb', nn_jdb, 'nn_kdb', nn_kdb, 'nn_narea_db', nn_narea_db 3638 3639 3613 iloc_db = mi0(nn_idb) 3640 3614 jloc_db = mj0(nn_jdb) … … 3642 3616 #endif 3643 3617 ENDIF 3644 3645 3618 ! 3646 3619 ! ! Check wave coupling settings ! 3647 3620 ! ! Further work needed - see ticket #2447 ! 3648 IF ( nn_osm_wave == 2 ) THEN3621 IF ( nn_osm_wave == 2 ) THEN 3649 3622 IF (.NOT. ( ln_wave .AND. ln_sdw )) & 3650 3623 & CALL ctl_stop( 'zdf_osm_init : ln_zdfosm and nn_osm_wave=2, ln_wave and ln_sdw must be true' ) 3651 3624 END IF 3652 3625 ! 3653 3626 ! Flags associated with diagnostic output 3654 3627 IF ( ln_dia_osm .AND. ( iom_use("zdudz_pyc") .OR. iom_use("zdvdz_pyc") ) ) ln_dia_pyc_shr = .TRUE. 3655 3628 IF ( ln_dia_osm .AND. ( iom_use("zdtdz_pyc") .OR. iom_use("zdsdz_pyc") .OR. iom_use("zdbdz_pyc" ) ) ) ln_dia_pyc_scl = .TRUE. 3656 3657 ! ! allocate zdfosm arrays 3658 IF( zdf_osm_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_osm_init : unable to allocate arrays' ) 3659 3660 3661 IF( ln_osm_mle ) THEN 3662 ! Initialise Fox-Kemper parametrization 3629 ! 3630 ! Allocate zdfosm arrays 3631 IF( zdf_osm_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_osm_init : unable to allocate arrays' ) 3632 ! 3633 IF( ln_osm_mle ) THEN ! Initialise Fox-Kemper parametrization 3663 3634 READ ( numnam_ref, namosm_mle, IOSTAT = ios, ERR = 903) 3664 903 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namosm_mle in reference namelist') 3665 3635 903 IF( ios /= 0 ) CALL ctl_nam( ios, 'namosm_mle in reference namelist' ) 3666 3636 READ ( numnam_cfg, namosm_mle, IOSTAT = ios, ERR = 904 ) 3667 904 IF( ios > 0 ) CALL ctl_nam ( ios , 'namosm_mle in configuration namelist')3637 904 IF( ios > 0 ) CALL ctl_nam( ios, 'namosm_mle in configuration namelist' ) 3668 3638 IF(lwm) WRITE ( numond, namosm_mle ) 3669 3670 IF(lwp) THEN 3639 ! 3640 IF(lwp) THEN ! Namelist print 3671 3641 WRITE(numout,*) 3672 3642 WRITE(numout,*) 'zdf_osm_init : initialise mixed layer eddy (MLE)' 3673 3643 WRITE(numout,*) '~~~~~~~~~~~~~' 3674 3644 WRITE(numout,*) ' Namelist namosm_mle : ' 3675 WRITE(numout,*) ' MLE type: =0 standard Fox-Kemper ; =1 new formulation nn_osm_mle = ', nn_osm_mle 3676 WRITE(numout,*) ' magnitude of the MLE (typical value: 0.06 to 0.08) rn_osm_mle_ce = ', rn_osm_mle_ce 3677 WRITE(numout,*) ' scale of ML front (ML radius of deformation) (nn_osm_mle=0) rn_osm_mle_lf = ', rn_osm_mle_lf, 'm' 3678 WRITE(numout,*) ' maximum time scale of MLE (nn_osm_mle=0) rn_osm_mle_time = ', rn_osm_mle_time, 's' 3679 WRITE(numout,*) ' reference latitude (degrees) of MLE coef. (nn_osm_mle=1) rn_osm_mle_lat = ', rn_osm_mle_lat, 'deg' 3680 WRITE(numout,*) ' Density difference used to define ML for FK rn_osm_mle_rho_c = ', rn_osm_mle_rho_c 3681 WRITE(numout,*) ' Threshold used to define MLE for FK rn_osm_mle_thresh = ', rn_osm_mle_thresh, 'm^2/s' 3682 WRITE(numout,*) ' Timescale for OSM-FK rn_osm_mle_tau = ', rn_osm_mle_tau, 's' 3683 WRITE(numout,*) ' switch to limit hmle ln_osm_hmle_limit = ', ln_osm_hmle_limit 3684 WRITE(numout,*) ' fraction of zmld to limit hmle to if ln_osm_hmle_limit =.T. rn_osm_hmle_limit = ', rn_osm_hmle_limit 3685 ENDIF ! 3686 ENDIF 3645 WRITE(numout,*) ' MLE type: =0 standard Fox-Kemper ; =1 new formulation nn_osm_mle = ', nn_osm_mle 3646 WRITE(numout,*) ' Magnitude of the MLE (typical value: 0.06 to 0.08) rn_osm_mle_ce = ', rn_osm_mle_ce 3647 WRITE(numout,*) ' Scale of ML front (ML radius of deform.) (nn_osm_mle=0) rn_osm_mle_lf = ', rn_osm_mle_lf, & 3648 & 'm' 3649 WRITE(numout,*) ' Maximum time scale of MLE (nn_osm_mle=0) rn_osm_mle_time = ', & 3650 & rn_osm_mle_time, 's' 3651 WRITE(numout,*) ' Reference latitude (deg) of MLE coef. (nn_osm_mle=1) rn_osm_mle_lat = ', rn_osm_mle_lat, & 3652 & 'deg' 3653 WRITE(numout,*) ' Density difference used to define ML for FK rn_osm_mle_rho_c = ', rn_osm_mle_rho_c 3654 WRITE(numout,*) ' Threshold used to define MLE for FK rn_osm_mle_thresh = ', & 3655 & rn_osm_mle_thresh, 'm^2/s' 3656 WRITE(numout,*) ' Timescale for OSM-FK rn_osm_mle_tau = ', rn_osm_mle_tau, 's' 3657 WRITE(numout,*) ' Switch to limit hmle ln_osm_hmle_limit = ', ln_osm_hmle_limit 3658 WRITE(numout,*) ' hmle limit (fraction of zmld) (ln_osm_hmle_limit = .T.) rn_osm_hmle_limit = ', rn_osm_hmle_limit 3659 END IF 3660 END IF 3687 3661 ! 3688 3662 IF(lwp) THEN 3689 3663 WRITE(numout,*) 3690 IF ( ln_osm_mle ) THEN3664 IF ( ln_osm_mle ) THEN 3691 3665 WRITE(numout,*) ' ==>>> Mixed Layer Eddy induced transport added to OSMOSIS BL calculation' 3692 3666 IF( nn_osm_mle == 0 ) WRITE(numout,*) ' Fox-Kemper et al 2010 formulation' … … 3694 3668 ELSE 3695 3669 WRITE(numout,*) ' ==>>> Mixed Layer induced transport NOT added to OSMOSIS BL calculation' 3696 END IF3697 END IF3698 ! 3699 IF( ln_osm_mle ) THEN 3670 END IF 3671 END IF 3672 ! 3673 IF( ln_osm_mle ) THEN ! MLE initialisation 3700 3674 ! 3701 rb_c = grav * rn_osm_mle_rho_c / rho0! Mixed Layer buoyancy criteria3675 rb_c = grav * rn_osm_mle_rho_c / rho0 ! Mixed Layer buoyancy criteria 3702 3676 IF(lwp) WRITE(numout,*) 3703 3677 IF(lwp) WRITE(numout,*) ' ML buoyancy criteria = ', rb_c, ' m/s2 ' 3704 3678 IF(lwp) WRITE(numout,*) ' associated ML density criteria defined in zdfmxl = ', rn_osm_mle_rho_c, 'kg/m3' 3705 3679 ! 3706 IF( nn_osm_mle == 0 ) THEN ! MLE array allocation & initialisation ! 3707 ! 3708 ELSEIF( nn_osm_mle == 1 ) THEN ! MLE array allocation & initialisation 3709 rc_f = rn_osm_mle_ce/ ( 5.e3_wp * 2._wp * omega * SIN( rad * rn_osm_mle_lat ) ) 3710 ! 3711 ENDIF 3712 ! ! 1/(f^2+tau^2)^1/2 at t-point (needed in both nn_osm_mle case) 3713 z1_t2 = 2.e-5 3680 IF( nn_osm_mle == 1 ) THEN 3681 rc_f = rn_osm_mle_ce / ( 5e3_wp * 2.0_wp * omega * SIN( rad * rn_osm_mle_lat ) ) 3682 END IF 3683 ! 1/(f^2+tau^2)^1/2 at t-point (needed in both nn_osm_mle case) 3684 z1_t2 = 2e-5_wp 3714 3685 DO_2D( 1, 1, 1, 1 ) 3715 r1_ft(ji,jj) = MIN( 1./( ABS(ff_t(ji,jj)) + epsln ), ABS(ff_t(ji,jj))/z1_t2**2)3686 r1_ft(ji,jj) = MIN( 1.0_wp / ( ABS( ff_t(ji,jj)) + epsln ), ABS( ff_t(ji,jj) ) / z1_t2**2 ) 3716 3687 END_2D 3717 3688 ! z1_t2 = 1._wp / ( rn_osm_mle_time * rn_osm_mle_timeji,jj ) 3718 3689 ! r1_ft(:,:) = 1._wp / SQRT( ff_t(:,:) * ff_t(:,:) + z1_t2 ) 3719 3690 ! 3720 ENDIF 3721 3722 call osm_rst( nit000, Kmm, 'READ' ) !* read or initialize hbl, dh, hmle 3723 3724 3725 IF( ln_zdfddm) THEN 3691 END IF 3692 ! 3693 CALL osm_rst( nit000, Kmm, 'READ' ) ! Read or initialize hbl, dh, hmle 3694 ! 3695 IF ( ln_zdfddm ) THEN 3726 3696 IF(lwp) THEN 3727 3697 WRITE(numout,*) 3728 3698 WRITE(numout,*) ' Double diffusion mixing on temperature and salinity ' 3729 3699 WRITE(numout,*) ' CAUTION : done in routine zdfosm, not in routine zdfddm ' 3730 ENDIF 3731 ENDIF 3732 3733 3734 !set constants not in namelist 3735 !----------------------------- 3736 3700 END IF 3701 END IF 3702 ! 3703 ! Set constants not in namelist 3704 ! ----------------------------- 3737 3705 IF(lwp) THEN 3738 3706 WRITE(numout,*) 3739 END IF3740 3741 dstokes(:,:) = zlarge3707 END IF 3708 ! 3709 dstokes(:,:) = pp_large 3742 3710 IF (nn_osm_wave == 0) THEN 3743 3711 dstokes(:,:) = rn_osm_dstokes 3744 3712 END IF 3745 3713 ! 3746 3714 ! Horizontal average : initialization of weighting arrays 3747 3715 ! ------------------- 3748 3749 3716 SELECT CASE ( nn_ave ) 3750 3751 3717 CASE ( 0 ) ! no horizontal average 3752 3718 IF(lwp) WRITE(numout,*) ' no horizontal average on avt' 3753 3719 IF(lwp) WRITE(numout,*) ' only in very high horizontal resolution !' 3754 ! weighting mean arrays etmean3720 ! Weighting mean arrays etmean 3755 3721 ! ( 1 1 ) 3756 3722 ! avt = 1/4 ( 1 1 ) 3757 3723 ! 3758 etmean(:,:,:) = 0. e03759 3724 etmean(:,:,:) = 0.0_wp 3725 ! 3760 3726 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 3761 etmean(ji,jj,jk) = tmask(ji,jj,jk) & 3762 & / MAX( 1., umask(ji-1,jj ,jk) + umask(ji,jj,jk) & 3763 & + vmask(ji ,jj-1,jk) + vmask(ji,jj,jk) ) 3727 etmean(ji,jj,jk) = tmask(ji,jj,jk) / MAX( 1.0_wp, umask(ji-1,jj, jk) + umask(ji,jj,jk) + & 3728 & vmask(ji, jj-1,jk) + vmask(ji,jj,jk) ) 3764 3729 END_3D 3765 3766 3730 CASE ( 1 ) ! horizontal average 3767 3731 IF(lwp) WRITE(numout,*) ' horizontal average on avt' 3768 ! weighting mean arrays etmean3732 ! Weighting mean arrays etmean 3769 3733 ! ( 1/2 1 1/2 ) 3770 3734 ! avt = 1/8 ( 1 2 1 ) 3771 3735 ! ( 1/2 1 1/2 ) 3772 etmean(:,:,:) = 0. e03773 3736 etmean(:,:,:) = 0.0_wp 3737 ! 3774 3738 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 3775 etmean(ji,jj,jk) = tmask(ji, jj,jk) & 3776 & / MAX( 1., 2.* tmask(ji,jj,jk) & 3777 & +.5 * ( tmask(ji-1,jj+1,jk) + tmask(ji-1,jj-1,jk) & 3778 & +tmask(ji+1,jj+1,jk) + tmask(ji+1,jj-1,jk) ) & 3779 & +1. * ( tmask(ji-1,jj ,jk) + tmask(ji ,jj+1,jk) & 3780 & +tmask(ji ,jj-1,jk) + tmask(ji+1,jj ,jk) ) ) 3739 etmean(ji,jj,jk) = tmask(ji, jj,jk) / MAX( 1.0_wp, 2.0_wp * tmask(ji,jj,jk) + & 3740 & 0.5_wp * ( tmask(ji-1,jj+1,jk) + tmask(ji-1,jj-1,jk) + & 3741 & tmask(ji+1,jj+1,jk) + tmask(ji+1,jj-1,jk) ) + & 3742 & 1.0_wp * ( tmask(ji-1,jj, jk) + tmask(ji, jj+1,jk) + & 3743 & tmask(ji, jj-1,jk) + tmask(ji+1,jj, jk) ) ) 3781 3744 END_3D 3782 3783 3745 CASE DEFAULT 3784 3746 WRITE(ctmp1,*) ' bad flag value for nn_ave = ', nn_ave 3785 3747 CALL ctl_stop( ctmp1 ) 3786 3787 3748 END SELECT 3788 3749 ! 3789 3750 ! Initialization of vertical eddy coef. to the background value 3790 3751 ! ------------------------------------------------------------- 3791 3752 DO jk = 1, jpk 3792 avt 3753 avt(:,:,jk) = avtb(jk) * tmask(:,:,jk) 3793 3754 END DO 3794 3795 ! zero the surface flux for non local term and osm mixed layer depth3755 ! 3756 ! Zero the surface flux for non local term and osm mixed layer depth 3796 3757 ! ------------------------------------------------------------------ 3797 ghamt(:,:,:) = 0. 3798 ghams(:,:,:) = 0. 3799 ghamu(:,:,:) = 0. 3800 ghamv(:,:,:) = 0. 3758 ghamt(:,:,:) = 0.0_wp 3759 ghams(:,:,:) = 0.0_wp 3760 ghamu(:,:,:) = 0.0_wp 3761 ghamv(:,:,:) = 0.0_wp 3801 3762 ! 3802 3763 IF( ln_timing ) CALL timing_stop('zdf_osm_init') 3764 ! 3803 3765 END SUBROUTINE zdf_osm_init 3804 3766 … … 3811 3773 !! ** Method : use of IOM library. If the restart does not contain 3812 3774 !! required fields, they are recomputed from stratification 3775 !! 3813 3776 !!---------------------------------------------------------------------- 3814 3815 INTEGER , INTENT(in) :: kt ! ocean time step index 3816 INTEGER , INTENT(in) :: Kmm ! ocean time level index (middle) 3817 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 3818 3819 INTEGER :: id1, id2, id3 ! iom enquiry index 3820 INTEGER :: ji, jj, jk ! dummy loop indices 3821 INTEGER :: iiki, ikt ! local integer 3822 REAL(wp) :: zhbf ! tempory scalars 3823 REAL(wp) :: zN2_c ! local scalar 3824 REAL(wp) :: rho_c = 0.01_wp !: density criterion for mixed layer depth 3825 INTEGER, DIMENSION(jpi,jpj) :: imld_rst ! level of mixed-layer depth (pycnocline top) 3826 !!---------------------------------------------------------------------- 3777 INTEGER , INTENT(in ) :: kt ! Ocean time step index 3778 INTEGER , INTENT(in ) :: Kmm ! Ocean time level index (middle) 3779 CHARACTER(len=*), INTENT(in ) :: cdrw ! "READ"/"WRITE" flag 3780 ! 3781 ! Local variables 3782 INTEGER :: id1, id2, id3 ! iom enquiry index 3783 INTEGER :: ji, jj, jk ! Dummy loop indices 3784 INTEGER :: iiki, ikt ! Local integer 3785 REAL(wp) :: zhbf ! Tempory scalars 3786 REAL(wp) :: zN2_c ! Local scalar 3787 REAL(wp) :: rho_c = 0.01_wp ! Density criterion for mixed layer depth 3788 INTEGER, DIMENSION(jpi,jpj) :: imld_rst ! Level of mixed-layer depth (pycnocline top) 3827 3789 ! 3828 3790 IF( ln_timing ) CALL timing_start('osm_rst') 3791 ! 3829 3792 !!----------------------------------------------------------------------------- 3830 3793 ! If READ/WRITE Flag is 'READ', try to get hbl from restart file. If successful then return 3831 3794 !!----------------------------------------------------------------------------- 3832 IF( TRIM(cdrw) == 'READ' .AND. ln_rstart) THEN3833 id1 = iom_varid( numror, 'wn' 3834 IF( id1 > 0 ) THEN 3795 IF( TRIM(cdrw) == 'READ' .AND. ln_rstart) THEN 3796 id1 = iom_varid( numror, 'wn', ldstop = .FALSE. ) 3797 IF( id1 > 0 ) THEN ! 'wn' exists; read 3835 3798 CALL iom_get( numror, jpdom_auto, 'wn', ww ) 3836 3799 WRITE(numout,*) ' ===>>>> : wn read from restart file' 3837 3800 ELSE 3838 ww(:,:,:) = 0. _wp3801 ww(:,:,:) = 0.0_wp 3839 3802 WRITE(numout,*) ' ===>>>> : wn not in restart file, set to zero initially' 3840 3803 END IF 3841 3842 id1 = iom_varid( numror, 'hbl' 3843 id2 = iom_varid( numror, 'dh' ,ldstop = .FALSE. )3844 IF( id1 > 0 .AND. id2 > 0 ) THEN! 'hbl' exists; read and return3845 CALL iom_get( numror, jpdom_auto, 'hbl' 3846 CALL iom_get( numror, jpdom_auto, 'dh', dh)3804 ! 3805 id1 = iom_varid( numror, 'hbl', ldstop = .FALSE. ) 3806 id2 = iom_varid( numror, 'dh', ldstop = .FALSE. ) 3807 IF( id1 > 0 .AND. id2 > 0 ) THEN ! 'hbl' exists; read and return 3808 CALL iom_get( numror, jpdom_auto, 'hbl', hbl ) 3809 CALL iom_get( numror, jpdom_auto, 'dh', dh ) 3847 3810 hml(:,:) = hbl(:,:) - dh(:,:) ! Initialise ML depth 3848 3811 WRITE(numout,*) ' ===>>>> : hbl & dh read from restart file' 3849 3812 IF( ln_osm_mle ) THEN 3850 id3 = iom_varid( numror, 'hmle' 3851 IF( id3 > 0 ) THEN3852 CALL iom_get( numror, jpdom_auto, 'hmle' 3813 id3 = iom_varid( numror, 'hmle', ldstop = .FALSE. ) 3814 IF( id3 > 0 ) THEN 3815 CALL iom_get( numror, jpdom_auto, 'hmle', hmle ) 3853 3816 WRITE(numout,*) ' ===>>>> : hmle read from restart file' 3854 3817 ELSE 3855 3818 WRITE(numout,*) ' ===>>>> : hmle not found, set to hbl' 3856 hmle(:,:) = hbl(:,:) ! Initialise MLE depth.3819 hmle(:,:) = hbl(:,:) ! Initialise MLE depth 3857 3820 END IF 3858 3821 END IF 3859 3822 RETURN 3860 ELSE 3823 ELSE ! 'hbl' & 'dh' not in restart file, recalculate 3861 3824 WRITE(numout,*) ' ===>>>> : previous run without osmosis scheme, hbl computed from stratification' 3862 3825 END IF 3863 3826 END IF 3864 3827 ! 3865 3828 !!----------------------------------------------------------------------------- 3866 3829 ! If READ/WRITE Flag is 'WRITE', write hbl into the restart file, then return 3867 3830 !!----------------------------------------------------------------------------- 3868 IF ( TRIM(cdrw) == 'WRITE') THEN !* Write hbl into the restart file, then return3831 IF ( TRIM(cdrw) == 'WRITE' ) THEN 3869 3832 IF(lwp) WRITE(numout,*) '---- osm-rst ----' 3870 CALL iom_rstput( kt, nitrst, numrow, 'wn' ,ww )3871 CALL iom_rstput( kt, nitrst, numrow, 'hbl' 3872 CALL iom_rstput( kt, nitrst, numrow, 'dh' ,dh )3873 IF ( ln_osm_mle ) THEN3833 CALL iom_rstput( kt, nitrst, numrow, 'wn', ww ) 3834 CALL iom_rstput( kt, nitrst, numrow, 'hbl', hbl ) 3835 CALL iom_rstput( kt, nitrst, numrow, 'dh', dh ) 3836 IF ( ln_osm_mle ) THEN 3874 3837 CALL iom_rstput( kt, nitrst, numrow, 'hmle', hmle ) 3875 3838 END IF 3876 3839 RETURN 3877 3840 END IF 3878 3841 ! 3879 3842 !!----------------------------------------------------------------------------- 3880 3843 ! Getting hbl, no restart file with hbl, so calculate from surface stratification … … 3883 3846 ! w-level of the mixing and mixed layers 3884 3847 CALL eos_rab( ts(:,:,:,:,Kmm), rab_n, Kmm ) 3885 CALL bn2( ts(:,:,:,:,Kmm), rab_n, rn2, Kmm)3886 imld_rst(:,:) = nlb10! Initialization to the number of w ocean point3887 hbl(:,:) = 0._wp ! here hbl used as a dummy variable, integrating vertically N^23888 zN2_c = grav * rho_c * r1_rho0 ! convert density criteria into N^2 criteria3889 ! 3890 hbl(:,:) = 0. _wp ! here hbl used as a dummy variable, integrating vertically N^23848 CALL bn2( ts(:,:,:,:,Kmm), rab_n, rn2, Kmm ) 3849 imld_rst(:,:) = nlb10 ! Initialization to the number of w ocean point 3850 hbl(:,:) = 0.0_wp ! Here hbl used as a dummy variable, integrating vertically N^2 3851 zN2_c = grav * rho_c * r1_rho0 ! Convert density criteria into N^2 criteria 3852 ! 3853 hbl(:,:) = 0.0_wp ! Here hbl used as a dummy variable, integrating vertically N^2 3891 3854 DO_3D( 1, 1, 1, 1, 1, jpkm1 ) 3892 3855 ikt = mbkt(ji,jj) 3893 hbl(ji,jj) = hbl(ji,jj) + MAX( rn2(ji,jj,jk) , 0. _wp ) * e3w(ji,jj,jk,Kmm)3894 IF ( hbl(ji,jj) < zN2_c )imld_rst(ji,jj) = MIN( jk , ikt ) + 1 ! Mixed layer level3856 hbl(ji,jj) = hbl(ji,jj) + MAX( rn2(ji,jj,jk) , 0.0_wp ) * e3w(ji,jj,jk,Kmm) 3857 IF ( hbl(ji,jj) < zN2_c ) imld_rst(ji,jj) = MIN( jk , ikt ) + 1 ! Mixed layer level 3895 3858 END_3D 3896 3859 ! 3897 3860 DO_2D( 1, 1, 1, 1 ) 3898 iiki = MAX( 4,imld_rst(ji,jj))3899 hbl(ji,jj) = gdepw(ji,jj,iiki,Kmm ) 3900 dh (ji,jj) = e3t(ji,jj,iiki-1,Kmm )! Turbocline depth3861 iiki = MAX( 4, imld_rst(ji,jj) ) 3862 hbl(ji,jj) = gdepw(ji,jj,iiki,Kmm ) ! Turbocline depth 3863 dh(ji,jj) = e3t(ji,jj,iiki-1,Kmm ) ! Turbocline depth 3901 3864 hml(ji,jj) = hbl(ji,jj) - dh(ji,jj) 3902 3865 END_2D 3903 3866 ! 3904 3867 WRITE(numout,*) ' ===>>>> : hbl computed from stratification' 3905 3868 ! 3906 3869 IF( ln_osm_mle ) THEN 3907 3870 hmle(:,:) = hbl(:,:) ! Initialise MLE depth. 3908 3871 WRITE(numout,*) ' ===>>>> : hmle set = to hbl' 3909 3872 END IF 3910 3873 ! 3911 3874 ww(:,:,:) = 0._wp 3912 3875 WRITE(numout,*) ' ===>>>> : wn not in restart file, set to zero initially' 3876 ! 3913 3877 IF( ln_timing ) CALL timing_stop('osm_rst') 3878 ! 3914 3879 END SUBROUTINE osm_rst 3915 3916 3880 3917 3881 SUBROUTINE tra_osm( kt, Kmm, pts, Krhs ) … … 3922 3886 !! 3923 3887 !! ** Method : ??? 3888 !! 3924 3889 !!---------------------------------------------------------------------- 3890 INTEGER , INTENT(in ) :: kt ! Time step index 3891 INTEGER , INTENT(in ) :: Kmm, Krhs ! Time level indices 3892 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! Active tracers and RHS of tracer equation 3893 ! 3894 ! Local variables 3895 INTEGER :: ji, jj, jk 3925 3896 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdt, ztrds ! 3D workspace 3926 !!----------------------------------------------------------------------3927 INTEGER , INTENT(in) :: kt ! time step index3928 INTEGER , INTENT(in) :: Kmm, Krhs ! time level indices3929 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers and RHS of tracer equation3930 !3931 INTEGER :: ji, jj, jk3932 3897 ! 3933 3898 IF( ln_timing ) CALL timing_start('tra_osm') 3934 IF( kt == nit000 ) THEN 3935 IF( ntile == 0 .OR. ntile == 1 ) THEN ! Do only on the first tile 3899 ! 3900 IF ( kt == nit000 ) THEN 3901 IF ( ntile == 0 .OR. ntile == 1 ) THEN ! Do only on the first tile 3936 3902 IF(lwp) WRITE(numout,*) 3937 3903 IF(lwp) WRITE(numout,*) 'tra_osm : OSM non-local tracer fluxes' 3938 3904 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 3939 ENDIF 3940 ENDIF 3941 3942 IF( l_trdtra ) THEN !* Save ta and sa trends 3943 ALLOCATE( ztrdt(jpi,jpj,jpk) ) ; ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) 3944 ALLOCATE( ztrds(jpi,jpj,jpk) ) ; ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) 3945 ENDIF 3946 3905 END IF 3906 END IF 3907 ! 3908 IF ( l_trdtra ) THEN ! Save ta and sa trends 3909 ALLOCATE( ztrdt(jpi,jpj,jpk), ztrds(jpi,jpj,jpk) ) 3910 ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) 3911 ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) 3912 END IF 3913 ! 3947 3914 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 3948 3915 pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) & … … 3953 3920 & - ghams(ji,jj,jk+1) ) / e3t(ji,jj,jk,Kmm) 3954 3921 END_3D 3955 3956 ! save the non-local tracer flux trends for diagnostics 3957 IF( l_trdtra ) THEN 3922 ! 3923 IF ( l_trdtra ) THEN ! Save the non-local tracer flux trends for diagnostics 3958 3924 ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:) 3959 3925 ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) - ztrds(:,:,:) 3960 3961 3926 CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_osm, ztrdt ) 3962 3927 CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_osm, ztrds ) 3963 DEALLOCATE( ztrdt ) ; DEALLOCATE(ztrds )3964 END IF3965 3966 IF (sn_cfctl%l_prtctl) THEN3928 DEALLOCATE( ztrdt, ztrds ) 3929 END IF 3930 ! 3931 IF ( sn_cfctl%l_prtctl ) THEN 3967 3932 CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' osm - Ta: ', mask1=tmask, & 3968 & 3969 END IF3933 & tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) 3934 END IF 3970 3935 ! 3971 3936 IF( ln_timing ) CALL timing_stop('tra_osm') 3937 ! 3972 3938 END SUBROUTINE tra_osm 3973 3939 3974 3975 SUBROUTINE trc_osm( kt ) ! Dummy routine 3940 SUBROUTINE trc_osm( kt ) ! Dummy routine 3976 3941 !!---------------------------------------------------------------------- 3977 3942 !! *** ROUTINE trc_osm *** … … 3982 3947 !! 3983 3948 !! ** Method : ??? 3984 !!---------------------------------------------------------------------- 3985 ! 3949 !! 3986 3950 !!---------------------------------------------------------------------- 3987 3951 INTEGER, INTENT(in) :: kt 3952 ! 3988 3953 IF( ln_timing ) CALL timing_start('trc_osm') 3954 ! 3989 3955 WRITE(*,*) 'trc_osm: Not written yet', kt 3956 ! 3990 3957 IF( ln_timing ) CALL timing_stop('trc_osm') 3958 ! 3991 3959 END SUBROUTINE trc_osm 3992 3993 3960 3994 3961 SUBROUTINE dyn_osm( kt, Kmm, puu, pvv, Krhs ) … … 4000 3967 !! 4001 3968 !! ** Method : ??? 3969 !! 4002 3970 !!---------------------------------------------------------------------- 4003 INTEGER , INTENT( in ) :: kt ! ocean time step index4004 INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices4005 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation3971 INTEGER , INTENT(in ) :: kt ! Ocean time step index 3972 INTEGER , INTENT(in ) :: Kmm, Krhs ! Ocean time level indices 3973 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! Ocean velocities and RHS of momentum equation 4006 3974 ! 4007 3975 INTEGER :: ji, jj, jk ! dummy loop indices … … 4009 3977 ! 4010 3978 IF( ln_timing ) CALL timing_start('dyn_osm') 4011 IF( kt == nit000 ) THEN 3979 ! 3980 IF ( kt == nit000 ) THEN 4012 3981 IF(lwp) WRITE(numout,*) 4013 3982 IF(lwp) WRITE(numout,*) 'dyn_osm : OSM non-local velocity' 4014 3983 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 4015 ENDIF 4016 !code saving tracer trends removed, replace with trdmxl_oce 4017 4018 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! add non-local u and v fluxes 4019 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) & 4020 & - ( ghamu(ji,jj,jk ) & 4021 & - ghamu(ji,jj,jk+1) ) / e3u(ji,jj,jk,Kmm) 4022 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) & 4023 & - ( ghamv(ji,jj,jk ) & 4024 & - ghamv(ji,jj,jk+1) ) / e3v(ji,jj,jk,Kmm) 3984 END IF 3985 ! 3986 ! Code saving tracer trends removed, replace with trdmxl_oce 3987 ! 3988 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! Add non-local u and v fluxes 3989 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( ghamu(ji,jj,jk ) - & 3990 & ghamu(ji,jj,jk+1) ) / e3u(ji,jj,jk,Kmm) 3991 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( ghamv(ji,jj,jk ) - & 3992 & ghamv(ji,jj,jk+1) ) / e3v(ji,jj,jk,Kmm) 4025 3993 END_3D 4026 3994 ! 4027 ! code for saving tracer trends removed3995 ! Code for saving tracer trends removed 4028 3996 ! 4029 3997 IF( ln_timing ) CALL timing_stop('dyn_osm') 3998 ! 4030 3999 END SUBROUTINE dyn_osm 4031 4000
Note: See TracChangeset
for help on using the changeset viewer.
