- Timestamp:
- 2019-12-11T12:02:38+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/ZDF/zdfosm.F90
r11250 r12178 25 25 !! (12) Replace zwstrl with zvstr in calculation of eddy viscosity. 26 26 !! 27/09/2017 (13) Calculate Stokes drift and Stokes penetration depth from wave information 27 !! (14) B uoyancy flux due to entrainment changed to include contribution from shear turbulence.27 !! (14) Bouyancy flux due to entrainment changed to include contribution from shear turbulence (for testing commented out). 28 28 !! 28/09/2017 (15) Calculation of Stokes drift moved into separate do-loops to allow for different options for the determining the Stokes drift to be added. 29 29 !! (16) Calculation of Stokes drift from windspeed for PM spectrum (for testing, commented out) 30 30 !! (17) Modification to Langmuir velocity scale to include effects due to the Stokes penetration depth (for testing, commented out) 31 !! ??/??/2018 (18) Revision to code structure, selected using key_osmldpth1. Inline code moved into subroutines. Changes to physics made,32 !! (a) Pycnocline temperature and salinity profies changed for unstable layers33 !! (b) The stable OSBL depth parametrization changed.34 !! 16/05/2019 (19) Fox-Kemper parametrization of restratification through mixed layer eddies added to revised code.35 !! 23/05/19 (20) Old code where key_osmldpth1` is *not* set removed, together with the key key_osmldpth136 31 !!---------------------------------------------------------------------- 37 32 … … 45 40 !! trc_osm : compute and add to the passive tracer trend the non-local flux (TBD) 46 41 !! dyn_osm : compute and add to u & v trensd the non-local flux 47 !!48 !! Subroutines in revised code.49 42 !!---------------------------------------------------------------------- 50 43 USE oce ! ocean dynamics and active tracers … … 58 51 USE traqsr ! details of solar radiation absorption 59 52 USE zdfddm ! double diffusion mixing (avs array) 60 USE zdfmxl ! mixed layer depth61 53 USE iom ! I/O library 62 54 USE lib_mpp ! MPP library … … 85 77 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: etmean !: averaging operator for avt 86 78 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hbl !: boundary layer depth 87 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: dh ! depth of pycnocline 88 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hml ! ML depth 79 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hbli !: intial boundary layer depth for stable blayer 89 80 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: dstokes !: penetration depth of the Stokes drift. 81 90 82 ! !!** Namelist namzdf_osm ** 91 83 LOGICAL :: ln_use_osm_la ! Use namelist rn_osm_la … … 105 97 106 98 ! !!! ** General constants ** 107 REAL(wp) :: epsln = 1.0e-20_wp ! a small positive number to ensure no div by zero 108 REAL(wp) :: depth_tol = 1.0e-6_wp ! a small-ish positive number to give a hbl slightly shallower than gdepw 99 REAL(wp) :: epsln = 1.0e-20_wp ! a small positive number 109 100 REAL(wp) :: pthird = 1._wp/3._wp ! 1/3 110 101 REAL(wp) :: p2third = 2._wp/3._wp ! 2/3 … … 123 114 !! *** FUNCTION zdf_osm_alloc *** 124 115 !!---------------------------------------------------------------------- 125 ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk),ghams(jpi,jpj,jpk), & 126 & hbl(jpi,jpj), dh(jpi,jpj), hml(jpi,jpj), dstokes(jpi, jpj), & 127 & etmean(jpi,jpj,jpk), STAT= zdf_osm_alloc ) 128 ! ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk),ghams(jpi,jpj,jpk), & ! would ths be better ? 129 ! & hbl(jpi,jpj), dh(jpi,jpj), hml(jpi,jpj), dstokes(jpi, jpj), & 130 ! & etmean(jpi,jpj,jpk), STAT= zdf_osm_alloc ) 131 ! IF( zdf_osm_alloc /= 0 ) CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm arrays') 132 ! IF ( ln_osm_mle ) THEN 133 ! Allocate( hmle(jpi,jpj), r1_ft(jpi,jpj), STAT= zdf_osm_alloc ) 134 ! IF( zdf_osm_alloc /= 0 ) CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm mle arrays') 135 ! ENDIF 136 116 ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk), ghams(jpi,jpj,jpk), & 117 & hbl(jpi,jpj) , hbli(jpi,jpj) , dstokes(jpi, jpj) , & 118 & etmean(jpi,jpj,jpk), STAT= zdf_osm_alloc ) 137 119 IF( zdf_osm_alloc /= 0 ) CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm arrays') 138 120 CALL mpp_sum ( 'zdfosm', zdf_osm_alloc ) … … 184 166 REAL(wp) :: zbeta, zthermal ! 185 167 REAL(wp) :: zehat, zeta, zhrib, zsig, zscale, zwst, zws, zwm ! Velocity scales 186 REAL(wp) :: zwsun, zwmun, zcons, zconm, zwcons, zwconm ! 187 168 REAL(wp) :: zwsun, zwmun, zcons, zconm, zwcons, zwconm ! 188 169 REAL(wp) :: zsr, zbw, ze, zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zcomp , zrhd,zrhdr,zbvzed ! In situ density 189 170 INTEGER :: jm ! dummy loop indices … … 210 191 REAL(wp), DIMENSION(jpi,jpj) :: zwbav ! Buoyancy flux - bl average 211 192 REAL(wp), DIMENSION(jpi,jpj) :: zwb_ent ! Buoyancy entrainment flux 212 213 193 REAL(wp), DIMENSION(jpi,jpj) :: zustke ! Surface Stokes drift 214 194 REAL(wp), DIMENSION(jpi,jpj) :: zla ! Trubulent Langmuir number … … 216 196 REAL(wp), DIMENSION(jpi,jpj) :: zsin_wind ! Sin angle of surface stress 217 197 REAL(wp), DIMENSION(jpi,jpj) :: zhol ! Stability parameter for boundary layer 218 LOGICAL, DIMENSION( jpi,jpj) :: lconv! unstable/stable bl198 LOGICAL, DIMENSION(:,:), ALLOCATABLE :: lconv ! unstable/stable bl 219 199 220 200 ! mixed-layer variables … … 222 202 INTEGER, DIMENSION(jpi,jpj) :: ibld ! level of boundary layer base 223 203 INTEGER, DIMENSION(jpi,jpj) :: imld ! level of mixed-layer depth (pycnocline top) 204 224 205 REAL(wp) :: ztgrad,zsgrad,zbgrad ! Temporary variables used to calculate pycnocline gradients 225 206 REAL(wp) :: zugrad,zvgrad ! temporary variables for calculating pycnocline shear … … 229 210 REAL(wp), DIMENSION(jpi,jpj) :: zdh ! pycnocline depth - grid 230 211 REAL(wp), DIMENSION(jpi,jpj) :: zdhdt ! BL depth tendency 231 REAL(wp), DIMENSION(jpi,jpj) :: zdhdt_2 ! correction to dhdt due to internal structure.232 REAL(wp), DIMENSION(jpi,jpj) :: zdtdz_ext,zdsdz_ext,zdbdz_ext ! external temperature/salinity and buoyancy gradients233 212 REAL(wp), DIMENSION(jpi,jpj) :: zt_bl,zs_bl,zu_bl,zv_bl,zrh_bl ! averages over the depth of the blayer 234 213 REAL(wp), DIMENSION(jpi,jpj) :: zt_ml,zs_ml,zu_ml,zv_ml,zrh_ml ! averages over the depth of the mixed layer … … 259 238 ! Temporary variables 260 239 INTEGER :: inhml 240 INTEGER :: i_lconv_alloc 261 241 REAL(wp) :: znd,znd_d,zznd_ml,zznd_pyc,zznd_d ! temporary non-dimensional depths used in various routines 262 242 REAL(wp) :: ztemp, zari, zpert, zzdhdt, zdb ! temporary variables … … 268 248 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdiffut ! t-diffusivity 269 249 270 INTEGER :: ibld_ext=0 ! does not have to be zero for modified scheme271 REAL(wp) :: zwb_min, zgamma_b_nd, zgamma_b, zdhoh, ztau, zddhdt272 REAL(wp) :: zzeta_s = 0._wp273 REAL(wp) :: zzeta_v = 0.46274 REAL(wp) :: zabsstke275 276 250 ! For debugging 277 251 INTEGER :: ikt 278 252 !!-------------------------------------------------------------------- 279 253 ! 254 ALLOCATE( lconv(jpi,jpj), STAT= i_lconv_alloc ) 255 IF( i_lconv_alloc /= 0 ) CALL ctl_warn('zdf_osm: failed to allocate lconv') 256 280 257 ibld(:,:) = 0 ; imld(:,:) = 0 281 258 zrad0(:,:) = 0._wp ; zradh(:,:) = 0._wp ; zradav(:,:) = 0._wp ; zustar(:,:) = 0._wp … … 291 268 zt_bl(:,:) = 0._wp ; zs_bl(:,:) = 0._wp ; zu_bl(:,:) = 0._wp ; zv_bl(:,:) = 0._wp 292 269 zrh_bl(:,:) = 0._wp ; zt_ml(:,:) = 0._wp ; zs_ml(:,:) = 0._wp ; zu_ml(:,:) = 0._wp 293 294 270 zv_ml(:,:) = 0._wp ; zrh_ml(:,:) = 0._wp ; zdt_bl(:,:) = 0._wp ; zds_bl(:,:) = 0._wp 295 271 zdu_bl(:,:) = 0._wp ; zdv_bl(:,:) = 0._wp ; zdrh_bl(:,:) = 0._wp ; zdb_bl(:,:) = 0._wp … … 301 277 zdudz_pyc(:,:,:) = 0._wp ; zdvdz_pyc(:,:,:) = 0._wp 302 278 ! 303 zdtdz_ext(:,:) = 0._wp ; zdsdz_ext(:,:) = 0._wp ; zdbdz_ext(:,:) = 0._wp304 279 ! Flux-Gradient arrays. 305 280 zdifml_sc(:,:) = 0._wp ; zvisml_sc(:,:) = 0._wp ; zdifpyc_sc(:,:) = 0._wp … … 312 287 ghams(:,:,:) = 0._wp ; ghamu(:,:,:) = 0._wp ; ghamv(:,:,:) = 0._wp 313 288 314 zdhdt_2(:,:) = 0._wp315 289 ! hbl = MAX(hbl,epsln) 316 290 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 367 341 zus_y = zsin_wind(ji,jj) * zustar(ji,jj) / 0.3**2 368 342 zustke(ji,jj) = MAX ( SQRT( zus_x*zus_x + zus_y*zus_y), 1.0e-8 ) 369 dstokes(ji,jj) = rn_osm_dstokes370 343 ! dstokes(ji,jj) set to constant value rn_osm_dstokes from namelist in zdf_osm_init 371 344 END DO … … 377 350 ! Use wind speed wndm included in sbc_oce module 378 351 zustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8 ) 379 dstokes(ji,jj) = MAX( 0.12 * wndm(ji,jj)**2 / grav, 5.e-1)352 dstokes(ji,jj) = 0.12 * wndm(ji,jj)**2 / grav 380 353 END DO 381 354 END DO … … 383 356 CASE(2) 384 357 zfac = 2.0_wp * rpi / 16.0_wp 385 386 358 DO jj = 2, jpjm1 387 359 DO ji = 2, jpim1 … … 390 362 ! It could represent the effects of the spread of wave directions 391 363 ! around the mean wind. The effect of this adjustment needs to be tested. 392 z absstke = SQRT(ut0sd(ji,jj)**2 + vt0sd(ji,jj)**2)393 zustke(ji,jj) = MAX (0.8 * ( zcos_wind(ji,jj) * ut0sd(ji,jj) + zsin_wind(ji,jj) * vt0sd(ji,jj) ), 1.0e-8)394 dstokes(ji,jj) = MAX(zfac * hsw(ji,jj)*hsw(ji,jj) / ( MAX(z absstke*wmp(ji,jj), 1.0e-7 ) ), 5.0e-1) !rn_osm_dstokes !364 zustke(ji,jj) = MAX ( 1.0 * ( zcos_wind(ji,jj) * ut0sd(ji,jj ) + zsin_wind(ji,jj) * vt0sd(ji,jj) ), & 365 & zustar(ji,jj) / ( 0.45 * 0.45 ) ) 366 dstokes(ji,jj) = MAX(zfac * hsw(ji,jj)*hsw(ji,jj) / ( MAX(zustke(ji,jj)*wmp(ji,jj), 1.0e-7 ) ), 5.0e-1) !rn_osm_dstokes ! 395 367 END DO 396 368 END DO … … 403 375 ! Langmuir velocity scale (zwstrl), at T-point 404 376 zwstrl(ji,jj) = ( zustar(ji,jj) * zustar(ji,jj) * zustke(ji,jj) )**pthird 405 zla(ji,jj) = MAX(MIN(SQRT ( zustar(ji,jj) / ( zwstrl(ji,jj) + epsln ) )**3, 4.0), 0.2) 406 IF(zla(ji,jj) > 0.45) dstokes(ji,jj) = MIN(dstokes(ji,jj), 0.5_wp*hbl(ji,jj)) 377 ! Modify zwstrl to allow for small and large values of dstokes/hbl. 378 ! Intended as a possible test. Doesn't affect LES results for entrainment, 379 ! but hasn't been shown to be correct as dstokes/h becomes large or small. 380 zwstrl(ji,jj) = zwstrl(ji,jj) * & 381 & (1.12 * ( 1.0 - ( 1.0 - EXP( -hbl(ji,jj) / dstokes(ji,jj) ) ) * dstokes(ji,jj) / hbl(ji,jj) ))**pthird * & 382 & ( 1.0 - EXP( -15.0 * dstokes(ji,jj) / hbl(ji,jj) )) 383 ! define La this way so effects of Stokes penetration depth on velocity scale are included 384 zla(ji,jj) = SQRT ( zustar(ji,jj) / zwstrl(ji,jj) )**3 407 385 ! Velocity scale that tends to zustar for large Langmuir numbers 408 386 zvstr(ji,jj) = ( zwstrl(ji,jj)**3 + & … … 411 389 ! limit maximum value of Langmuir number as approximate treatment for shear turbulence. 412 390 ! Note zustke and zwstrl are not amended. 391 IF ( zla(ji,jj) >= 0.45 ) zla(ji,jj) = 0.45 413 392 ! 414 393 ! get convective velocity (zwstrc), stabilty scale (zhol) and logical conection flag lconv … … 427 406 ! Mixed-layer model - calculate averages over the boundary layer, and the change in the boundary layer depth 428 407 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 429 ! BL must be always 4levels deep.430 hbl(:,:) = MAX(hbl(:,:), gdepw_n(:,:, 4) )431 ibld(:,:) = 4432 DO jk = 5, jpkm1408 ! BL must be always 2 levels deep. 409 hbl(:,:) = MAX(hbl(:,:), gdepw_n(:,:,3) ) 410 ibld(:,:) = 3 411 DO jk = 4, jpkm1 433 412 DO jj = 2, jpjm1 434 413 DO ji = 2, jpim1 … … 440 419 END DO 441 420 442 DO jj = 2, jpjm1 421 DO jj = 2, jpjm1 ! Vertical slab 443 422 DO ji = 2, jpim1 444 zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj)) 445 imld(ji,jj) = MAX(3,ibld(ji,jj) - MAX( INT( dh(ji,jj) / e3t_n(ji, jj, ibld(ji,jj) )) , 1 )) 446 zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 447 END DO 423 zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 424 zbeta = rab_n(ji,jj,1,jp_sal) 425 zt = 0._wp 426 zs = 0._wp 427 zu = 0._wp 428 zv = 0._wp 429 ! average over depth of boundary layer 430 zthick=0._wp 431 DO jm = 2, ibld(ji,jj) 432 zthick=zthick+e3t_n(ji,jj,jm) 433 zt = zt + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem) 434 zs = zs + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal) 435 zu = zu + e3t_n(ji,jj,jm) & 436 & * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) & 437 & / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 438 zv = zv + e3t_n(ji,jj,jm) & 439 & * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) & 440 & / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 441 END DO 442 zt_bl(ji,jj) = zt / zthick 443 zs_bl(ji,jj) = zs / zthick 444 zu_bl(ji,jj) = zu / zthick 445 zv_bl(ji,jj) = zv / zthick 446 zdt_bl(ji,jj) = zt_bl(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_tem) 447 zds_bl(ji,jj) = zs_bl(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal) 448 zdu_bl(ji,jj) = zu_bl(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) & 449 & / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 450 zdv_bl(ji,jj) = zv_bl(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) & 451 & / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 452 zdb_bl(ji,jj) = grav * zthermal * zdt_bl(ji,jj) - grav * zbeta * zds_bl(ji,jj) 453 IF ( lconv(ji,jj) ) THEN ! Convective 454 zwb_ent(ji,jj) = -( 2.0 * 0.2 * zwbav(ji,jj) & 455 & + 0.135 * zla(ji,jj) * zwstrl(ji,jj)**3/hbl(ji,jj) ) 456 457 zvel_max = - ( 1.0 + 1.0 * ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / hbl(ji,jj) ) & 458 & * zwb_ent(ji,jj) / ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 459 ! Entrainment including component due to shear turbulence. Modified Langmuir component, but gives same result for La=0.3 For testing uncomment. 460 ! zwb_ent(ji,jj) = -( 2.0 * 0.2 * zwbav(ji,jj) & 461 ! & + ( 0.15 * ( 1.0 - EXP( -0.5 * zla(ji,jj) ) ) + 0.03 / zla(ji,jj)**2 ) * zustar(ji,jj)**3/hbl(ji,jj) ) 462 463 ! zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / zhbl(ji,jj) ) * zwb_ent(ji,jj) / & 464 ! & ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 465 zzdhdt = - zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj),0.0) ) 466 ELSE ! Stable 467 zzdhdt = 0.32 * ( hbli(ji,jj) / hbl(ji,jj) -1.0 ) * zwstrl(ji,jj)**3 / hbli(ji,jj) & 468 & + ( ( 0.32 / 3.0 ) * exp ( -2.5 * ( hbli(ji,jj) / hbl(ji,jj) - 1.0 ) ) & 469 & - ( 0.32 / 3.0 - 0.135 * zla(ji,jj) ) * exp ( -12.5 * ( hbli(ji,jj) / hbl(ji,jj) ) ) ) & 470 & * zwstrl(ji,jj)**3 / hbli(ji,jj) 471 zzdhdt = zzdhdt + zwbav(ji,jj) 472 IF ( zzdhdt < 0._wp ) THEN 473 ! For long timsteps factor in brackets slows the rapid collapse of the OSBL 474 zpert = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_rdt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj) 475 ELSE 476 zpert = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_rdt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj) & 477 & + MAX( zdb_bl(ji,jj), 0.0 ) 478 ENDIF 479 zzdhdt = 2.0 * zzdhdt / zpert 480 ENDIF 481 zdhdt(ji,jj) = zzdhdt 482 END DO 448 483 END DO 449 ! Averages over well-mixed and boundary layer 450 ! Alan: do we need zb_nl?, zb_ml? 451 CALL zdf_osm_vertical_average(ibld, zt_bl, zs_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl) 452 CALL zdf_osm_vertical_average(imld, zt_ml, zs_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml) 453 ! External gradient 454 CALL zdf_osm_external_gradients( zdtdz_ext, zdsdz_ext, zdbdz_ext ) 455 456 ! Rate of change of hbl 457 CALL zdf_osm_calculate_dhdt( zdhdt, zdhdt_2 ) 484 458 485 ! Calculate averages over depth of boundary layer 459 486 imld = ibld ! use imld to hold previous blayer index 460 ibld(:,:) = 4 461 462 DO jj = 2, jpjm1 463 DO ji = 2, jpim1 464 zhbl_t(ji,jj) = hbl(ji,jj) + (zdhdt(ji,jj) - wn(ji,jj,ibld(ji,jj)))* rn_rdt ! certainly need wn here, so subtract it 465 zhbl_t(ji,jj) = MIN(zhbl_t(ji,jj), gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol)! ht_n(:,:)) 466 !zdhdt(ji,jj) = MIN(zdhdt(ji,jj), (zhbl_t(ji,jj) - hbl(ji,jj))/rn_rdt + wn(ji,jj,ibld(ji,jj))) 467 END DO 468 END DO 469 ! adjustment to represent limiting by ocean bottom 487 ibld(:,:) = 3 488 489 zhbl_t(:,:) = hbl(:,:) + (zdhdt(:,:) - wn(ji,jj,ibld(ji,jj)))* rn_rdt ! certainly need wb here, so subtract it 490 zhbl_t(:,:) = MIN(zhbl_t(:,:), ht_n(:,:)) 491 zdhdt(:,:) = MIN(zdhdt(:,:), (zhbl_t(:,:) - hbl(:,:))/rn_rdt + wn(ji,jj,ibld(ji,jj))) ! adjustment to represent limiting by ocean bottom 470 492 471 493 DO jk = 4, jpkm1 … … 473 495 DO ji = 2, jpim1 474 496 IF ( zhbl_t(ji,jj) >= gdepw_n(ji,jj,jk) ) THEN 475 ibld(ji,jj) = jk497 ibld(ji,jj) = MIN(mbkt(ji,jj), jk) 476 498 ENDIF 477 499 END DO … … 482 504 ! Step through model levels taking account of buoyancy change to determine the effect on dhdt 483 505 ! 484 CALL zdf_osm_timestep_hbl( zdhdt, zdhdt_2 )485 ! Alan: do we need zb_ml?486 CALL zdf_osm_vertical_average( ibld, zt_bl, zs_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl )506 DO jj = 2, jpjm1 507 DO ji = 2, jpim1 508 IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN 487 509 ! 510 ! If boundary layer changes by more than one level, need to check for stable layers between initial and final depths. 488 511 ! 489 CALL zdf_osm_pycnocline_thickness( dh, zdh ) 512 zhbl_s = hbl(ji,jj) 513 jm = imld(ji,jj) 514 zthermal = rab_n(ji,jj,1,jp_tem) 515 zbeta = rab_n(ji,jj,1,jp_sal) 516 IF ( lconv(ji,jj) ) THEN 517 !unstable 518 zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / hbl(ji,jj) ) & 519 & * zwb_ent(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 520 521 DO jk = imld(ji,jj), ibld(ji,jj) 522 zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) ) & 523 & - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ), 0.0 ) + zvel_max 524 525 zhbl_s = zhbl_s + MIN( - zwb_ent(ji,jj) / zdb * rn_rdt / FLOAT(ibld(ji,jj)-imld(ji,jj) ), e3w_n(ji,jj,jk) ) 526 zhbl_s = MIN(zhbl_s, ht_n(ji,jj)) 527 528 IF ( zhbl_s >= gdepw_n(ji,jj,jm+1) ) jm = jm + 1 529 END DO 530 hbl(ji,jj) = zhbl_s 531 ibld(ji,jj) = jm 532 hbli(ji,jj) = hbl(ji,jj) 533 ELSE 534 ! stable 535 DO jk = imld(ji,jj), ibld(ji,jj) 536 zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) ) & 537 & - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ), 0.0 ) & 538 & + 2.0 * zwstrl(ji,jj)**2 / zhbl_s 539 540 zhbl_s = zhbl_s + ( & 541 & 0.32 * ( hbli(ji,jj) / zhbl_s -1.0 ) & 542 & * zwstrl(ji,jj)**3 / hbli(ji,jj) & 543 & + ( ( 0.32 / 3.0 ) * EXP( - 2.5 * ( hbli(ji,jj) / zhbl_s -1.0 ) ) & 544 & - ( 0.32 / 3.0 - 0.0485 ) * EXP( - 12.5 * ( hbli(ji,jj) / zhbl_s ) ) ) & 545 & * zwstrl(ji,jj)**3 / hbli(ji,jj) ) / zdb * e3w_n(ji,jj,jk) / zdhdt(ji,jj) ! ALMG to investigate whether need to include wn here 546 547 zhbl_s = MIN(zhbl_s, ht_n(ji,jj)) 548 IF ( zhbl_s >= gdepw_n(ji,jj,jm) ) jm = jm + 1 549 END DO 550 hbl(ji,jj) = MAX(zhbl_s, gdepw_n(ji,jj,3) ) 551 ibld(ji,jj) = MAX(jm, 3 ) 552 IF ( hbl(ji,jj) > hbli(ji,jj) ) hbli(ji,jj) = hbl(ji,jj) 553 ENDIF ! IF ( lconv ) 554 ELSE 555 ! change zero or one model level. 556 hbl(ji,jj) = zhbl_t(ji,jj) 557 IF ( lconv(ji,jj) ) THEN 558 hbli(ji,jj) = hbl(ji,jj) 559 ELSE 560 hbl(ji,jj) = MAX(hbl(ji,jj), gdepw_n(ji,jj,3) ) 561 IF ( hbl(ji,jj) > hbli(ji,jj) ) hbli(ji,jj) = hbl(ji,jj) 562 ENDIF 563 ENDIF 564 zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj)) 565 END DO 566 END DO 490 567 dstokes(:,:) = MIN ( dstokes(:,:), hbl(:,:)/3. ) ! Limit delta for shallow boundary layers for calculating flux-gradient terms. 491 ! 492 ! Average over the depth of the mixed layer in the convective boundary layer 493 ! Alan: do we need zb_ml? 494 CALL zdf_osm_vertical_average( imld, zt_ml, zs_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml ) 568 569 ! Recalculate averages over boundary layer after depth updated 570 ! Consider later combining this into the loop above and looking for columns 571 ! where the index for base of the boundary layer have changed 572 DO jj = 2, jpjm1 ! Vertical slab 573 DO ji = 2, jpim1 574 zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 575 zbeta = rab_n(ji,jj,1,jp_sal) 576 zt = 0._wp 577 zs = 0._wp 578 zu = 0._wp 579 zv = 0._wp 580 ! average over depth of boundary layer 581 zthick=0._wp 582 DO jm = 2, ibld(ji,jj) 583 zthick=zthick+e3t_n(ji,jj,jm) 584 zt = zt + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem) 585 zs = zs + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal) 586 zu = zu + e3t_n(ji,jj,jm) & 587 & * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) & 588 & / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 589 zv = zv + e3t_n(ji,jj,jm) & 590 & * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) & 591 & / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 592 END DO 593 zt_bl(ji,jj) = zt / zthick 594 zs_bl(ji,jj) = zs / zthick 595 zu_bl(ji,jj) = zu / zthick 596 zv_bl(ji,jj) = zv / zthick 597 zdt_bl(ji,jj) = zt_bl(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_tem) 598 zds_bl(ji,jj) = zs_bl(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal) 599 zdu_bl(ji,jj) = zu_bl(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) & 600 & / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 601 zdv_bl(ji,jj) = zv_bl(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) & 602 & / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 603 zdb_bl(ji,jj) = grav * zthermal * zdt_bl(ji,jj) - grav * zbeta * zds_bl(ji,jj) 604 zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj)) 605 IF ( lconv(ji,jj) ) THEN 606 IF ( zdb_bl(ji,jj) > 0._wp )THEN 607 IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN ! near neutral stability 608 zari = 4.5 * ( zvstr(ji,jj)**2 ) & 609 & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01 610 ELSE ! unstable 611 zari = 4.5 * ( zwstrc(ji,jj)**2 ) & 612 & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01 613 ENDIF 614 IF ( zari > 0.2 ) THEN ! This test checks for weakly stratified pycnocline 615 zari = 0.2 616 zwb_ent(ji,jj) = 0._wp 617 ENDIF 618 inhml = MAX( INT( zari * zhbl(ji,jj) / e3t_n(ji,jj,ibld(ji,jj)) ) , 1 ) 619 imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 1) 620 zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 621 zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 622 ELSE ! IF (zdb_bl) 623 imld(ji,jj) = ibld(ji,jj) - 1 624 zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 625 zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 626 ENDIF 627 ELSE ! IF (lconv) 628 IF ( zdhdt(ji,jj) >= 0.0 ) THEN ! probably shouldn't include wm here 629 ! boundary layer deepening 630 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 631 ! pycnocline thickness set by stratification - use same relationship as for neutral conditions. 632 zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) & 633 & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01 , 0.2 ) 634 inhml = MAX( INT( zari * zhbl(ji,jj) / e3t_n(ji,jj,ibld(ji,jj)) ) , 1 ) 635 imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 1) 636 zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 637 zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 638 ELSE 639 imld(ji,jj) = ibld(ji,jj) - 1 640 zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 641 zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 642 ENDIF ! IF (zdb_bl > 0.0) 643 ELSE ! IF(dhdt >= 0) 644 ! boundary layer collapsing. 645 imld(ji,jj) = ibld(ji,jj) 646 zhml(ji,jj) = zhbl(ji,jj) 647 zdh(ji,jj) = 0._wp 648 ENDIF ! IF (dhdt >= 0) 649 ENDIF ! IF (lconv) 650 END DO 651 END DO 652 653 ! Average over the depth of the mixed layer in the convective boundary layer 654 ! Also calculate entrainment fluxes for temperature and salinity 655 DO jj = 2, jpjm1 ! Vertical slab 656 DO ji = 2, jpim1 657 zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 658 zbeta = rab_n(ji,jj,1,jp_sal) 659 IF ( lconv(ji,jj) ) THEN 660 zt = 0._wp 661 zs = 0._wp 662 zu = 0._wp 663 zv = 0._wp 664 ! average over depth of boundary layer 665 zthick=0._wp 666 DO jm = 2, imld(ji,jj) 667 zthick=zthick+e3t_n(ji,jj,jm) 668 zt = zt + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem) 669 zs = zs + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal) 670 zu = zu + e3t_n(ji,jj,jm) & 671 & * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) & 672 & / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 673 zv = zv + e3t_n(ji,jj,jm) & 674 & * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) & 675 & / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 676 END DO 677 zt_ml(ji,jj) = zt / zthick 678 zs_ml(ji,jj) = zs / zthick 679 zu_ml(ji,jj) = zu / zthick 680 zv_ml(ji,jj) = zv / zthick 681 zdt_ml(ji,jj) = zt_ml(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_tem) 682 zds_ml(ji,jj) = zs_ml(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal) 683 zdu_ml(ji,jj) = zu_ml(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) & 684 & / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 685 zdv_ml(ji,jj) = zv_ml(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) & 686 & / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 687 zdb_ml(ji,jj) = grav * zthermal * zdt_ml(ji,jj) - grav * zbeta * zds_ml(ji,jj) 688 ELSE 689 ! stable, if entraining calulate average below interface layer. 690 IF ( zdhdt(ji,jj) >= 0._wp ) THEN 691 zt = 0._wp 692 zs = 0._wp 693 zu = 0._wp 694 zv = 0._wp 695 ! average over depth of boundary layer 696 zthick=0._wp 697 DO jm = 2, imld(ji,jj) 698 zthick=zthick+e3t_n(ji,jj,jm) 699 zt = zt + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem) 700 zs = zs + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal) 701 zu = zu + e3t_n(ji,jj,jm) & 702 & * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) & 703 & / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 704 zv = zv + e3t_n(ji,jj,jm) & 705 & * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) & 706 & / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 707 END DO 708 zt_ml(ji,jj) = zt / zthick 709 zs_ml(ji,jj) = zs / zthick 710 zu_ml(ji,jj) = zu / zthick 711 zv_ml(ji,jj) = zv / zthick 712 zdt_ml(ji,jj) = zt_ml(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_tem) 713 zds_ml(ji,jj) = zs_ml(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal) 714 zdu_ml(ji,jj) = zu_ml(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) & 715 & / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 716 zdv_ml(ji,jj) = zv_ml(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) & 717 & / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 718 zdb_ml(ji,jj) = grav * zthermal * zdt_ml(ji,jj) - grav * zbeta * zds_ml(ji,jj) 719 ENDIF 720 ENDIF 721 END DO 722 END DO 723 ! 495 724 ! rotate mean currents and changes onto wind align co-ordinates 496 725 ! 497 CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_ml, zv_ml, zdu_ml, zdv_ml ) 498 CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_bl, zv_bl, zdu_bl, zdv_bl ) 726 727 DO jj = 2, jpjm1 728 DO ji = 2, jpim1 729 ztemp = zu_ml(ji,jj) 730 zu_ml(ji,jj) = zu_ml(ji,jj) * zcos_wind(ji,jj) + zv_ml(ji,jj) * zsin_wind(ji,jj) 731 zv_ml(ji,jj) = zv_ml(ji,jj) * zcos_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 732 ztemp = zdu_ml(ji,jj) 733 zdu_ml(ji,jj) = zdu_ml(ji,jj) * zcos_wind(ji,jj) + zdv_ml(ji,jj) * zsin_wind(ji,jj) 734 zdv_ml(ji,jj) = zdv_ml(ji,jj) * zsin_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 735 ! 736 ztemp = zu_bl(ji,jj) 737 zu_bl = zu_bl(ji,jj) * zcos_wind(ji,jj) + zv_bl(ji,jj) * zsin_wind(ji,jj) 738 zv_bl(ji,jj) = zv_bl(ji,jj) * zcos_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 739 ztemp = zdu_bl(ji,jj) 740 zdu_bl(ji,jj) = zdu_bl(ji,jj) * zcos_wind(ji,jj) + zdv_bl(ji,jj) * zsin_wind(ji,jj) 741 zdv_bl(ji,jj) = zdv_bl(ji,jj) * zsin_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 742 END DO 743 END DO 744 499 745 zuw_bse = 0._wp 500 746 zvw_bse = 0._wp 501 zwth_ent = 0._wp502 zws_ent = 0._wp503 747 DO jj = 2, jpjm1 504 748 DO ji = 2, jpim1 505 IF ( ibld(ji,jj) < mbkt(ji,jj) ) THEN 506 IF ( lconv(ji,jj) ) THEN 507 zuw_bse(ji,jj) = -0.0075*((zvstr(ji,jj)**3+0.5*zwstrc(ji,jj)**3)**pthird*zdu_ml(ji,jj) + & 508 & 1.5*zustar(ji,jj)**2*(zhbl(ji,jj)-zhml(ji,jj)) )/ & 509 & ( zhml(ji,jj)*MIN(zla(ji,jj)**(8./3.),1.) + epsln) 510 zvw_bse(ji,jj) = 0.01*(-(zvstr(ji,jj)**3+0.5*zwstrc(ji,jj)**3)**pthird*zdv_ml(ji,jj)+ & 511 & 2.0*ff_t(ji,jj)*zustke(ji,jj)*dstokes(ji,jj)*zla(ji,jj)) 512 IF ( zdb_ml(ji,jj) > 0._wp ) THEN 513 zwth_ent(ji,jj) = zwb_ent(ji,jj) * zdt_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 514 zws_ent(ji,jj) = zwb_ent(ji,jj) * zds_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 515 ENDIF 516 ELSE 517 zwth_ent(ji,jj) = -2.0 * zwthav(ji,jj) * ( (1.0 - 0.8) - ( 1.0 - 0.8)**(3.0/2.0) ) 518 zws_ent(ji,jj) = -2.0 * zwsav(ji,jj) * ( (1.0 - 0.8 ) - ( 1.0 - 0.8 )**(3.0/2.0) ) 749 750 IF ( lconv(ji,jj) ) THEN 751 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 752 zwth_ent(ji,jj) = zwb_ent(ji,jj) * zdt_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 753 zws_ent(ji,jj) = zwb_ent(ji,jj) * zds_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 519 754 ENDIF 755 ELSE 756 zwth_ent(ji,jj) = -2.0 * zwthav(ji,jj) * ( (1.0 - 0.8) - ( 1.0 - 0.8)**(3.0/2.0) ) 757 zws_ent(ji,jj) = -2.0 * zwsav(ji,jj) * ( (1.0 - 0.8 ) - ( 1.0 - 0.8 )**(3.0/2.0) ) 520 758 ENDIF 521 759 END DO … … 526 764 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 527 765 528 CALL zdf_osm_external_gradients( zdtdz_ext, zdsdz_ext, zdbdz_ext ) 529 CALL zdf_osm_pycnocline_scalar_profiles( zdtdz_pyc, zdsdz_pyc, zdbdz_pyc ) 530 CALL zdf_osm_pycnocline_shear_profiles( zdudz_pyc, zdvdz_pyc ) 766 DO jj = 2, jpjm1 767 DO ji = 2, jpim1 768 ! 769 IF ( lconv (ji,jj) ) THEN 770 ! Unstable conditions 771 IF( zdb_bl(ji,jj) > 0._wp ) THEN 772 ! calculate pycnocline profiles, no need if zdb_bl <= 0. since profile is zero and arrays have been initialized to zero 773 ztgrad = ( zdt_ml(ji,jj) / zdh(ji,jj) ) 774 zsgrad = ( zds_ml(ji,jj) / zdh(ji,jj) ) 775 zbgrad = ( zdb_ml(ji,jj) / zdh(ji,jj) ) 776 DO jk = 2 , ibld(ji,jj) 777 znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 778 zdtdz_pyc(ji,jj,jk) = ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 779 zdbdz_pyc(ji,jj,jk) = zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 780 zdsdz_pyc(ji,jj,jk) = zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 781 END DO 782 ENDIF 783 ELSE 784 ! stable conditions 785 ! if pycnocline profile only defined when depth steady of increasing. 786 IF ( zdhdt(ji,jj) >= 0.0 ) THEN ! Depth increasing, or steady. 787 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 788 IF ( zhol(ji,jj) >= 0.5 ) THEN ! Very stable - 'thick' pycnocline 789 ztgrad = zdt_bl(ji,jj) / zhbl(ji,jj) 790 zsgrad = zds_bl(ji,jj) / zhbl(ji,jj) 791 zbgrad = zdb_bl(ji,jj) / zhbl(ji,jj) 792 DO jk = 2, ibld(ji,jj) 793 znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 794 zdtdz_pyc(ji,jj,jk) = ztgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 795 zdbdz_pyc(ji,jj,jk) = zbgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 796 zdsdz_pyc(ji,jj,jk) = zsgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 797 END DO 798 ELSE ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form. 799 ztgrad = zdt_bl(ji,jj) / zdh(ji,jj) 800 zsgrad = zds_bl(ji,jj) / zdh(ji,jj) 801 zbgrad = zdb_bl(ji,jj) / zdh(ji,jj) 802 DO jk = 2, ibld(ji,jj) 803 znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 804 zdtdz_pyc(ji,jj,jk) = ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 805 zdbdz_pyc(ji,jj,jk) = zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 806 zdsdz_pyc(ji,jj,jk) = zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 807 END DO 808 ENDIF ! IF (zhol >=0.5) 809 ENDIF ! IF (zdb_bl> 0.) 810 ENDIF ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero, profile arrays are intialized to zero 811 ENDIF ! IF (lconv) 812 ! 813 END DO 814 END DO 815 ! 816 DO jj = 2, jpjm1 817 DO ji = 2, jpim1 818 ! 819 IF ( lconv (ji,jj) ) THEN 820 ! Unstable conditions 821 zugrad = ( zdu_ml(ji,jj) / zdh(ji,jj) ) + 0.275 * zustar(ji,jj)*zustar(ji,jj) / & 822 & (( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) / zla(ji,jj)**(8.0/3.0) 823 zvgrad = ( zdv_ml(ji,jj) / zdh(ji,jj) ) + 3.5 * ff_t(ji,jj) * zustke(ji,jj) / & 824 & ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 825 DO jk = 2 , ibld(ji,jj)-1 826 znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 827 zdudz_pyc(ji,jj,jk) = zugrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 828 zdvdz_pyc(ji,jj,jk) = zvgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 829 END DO 830 ELSE 831 ! stable conditions 832 zugrad = 3.25 * zdu_bl(ji,jj) / zhbl(ji,jj) 833 zvgrad = 2.75 * zdv_bl(ji,jj) / zhbl(ji,jj) 834 DO jk = 2, ibld(ji,jj) 835 znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 836 IF ( znd < 1.0 ) THEN 837 zdudz_pyc(ji,jj,jk) = zugrad * EXP( -40.0 * ( znd - 1.0 )**2 ) 838 ELSE 839 zdudz_pyc(ji,jj,jk) = zugrad * EXP( -20.0 * ( znd - 1.0 )**2 ) 840 ENDIF 841 zdvdz_pyc(ji,jj,jk) = zvgrad * EXP( -20.0 * ( znd - 0.85 )**2 ) 842 END DO 843 ENDIF 844 ! 845 END DO 846 END DO 531 847 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 532 848 ! Eddy viscosity/diffusivity and non-gradient terms in the flux-gradient relationship 533 849 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 534 850 851 ! WHERE ( lconv ) 852 ! zdifml_sc = zhml * ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird 853 ! zvisml_sc = zdifml_sc 854 ! zdifpyc_sc = 0.165 * ( zwstrl**3 + zwstrc**3 )**pthird * ( zhbl - zhml ) 855 ! zvispyc_sc = 0.142 * ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * ( zhbl - zhml ) 856 ! zbeta_d_sc = 1.0 - (0.165 / 0.8 * ( zhbl - zhml ) / zhbl )**p2third 857 ! zbeta_v_sc = 1.0 - 2.0 * (0.142 /0.375) * (zhbl - zhml ) / zhml 858 ! ELSEWHERE 859 ! zdifml_sc = zwstrl * zhbl * EXP ( -( zhol / 0.183_wp )**2 ) 860 ! zvisml_sc = zwstrl * zhbl * EXP ( -( zhol / 0.183_wp )**2 ) 861 ! ENDWHERE 535 862 DO jj = 2, jpjm1 536 863 DO ji = 2, jpim1 … … 541 868 zvispyc_sc(ji,jj) = 0.142 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zdh(ji,jj) 542 869 zbeta_d_sc(ji,jj) = 1.0 - (0.165 / 0.8 * zdh(ji,jj) / zhbl(ji,jj) )**p2third 543 zbeta_v_sc(ji,jj) = 1.0 - 2.0 * (0.142 /0.375) * zdh(ji,jj) / ( zhml(ji,jj) + epsln)870 zbeta_v_sc(ji,jj) = 1.0 - 2.0 * (0.142 /0.375) * zdh(ji,jj) / zhml(ji,jj) 544 871 ELSE 545 872 zdifml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ) 546 873 zvisml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ) 547 548 549 874 END IF 875 END DO 876 END DO 550 877 ! 551 878 DO jj = 2, jpjm1 … … 562 889 ! pycnocline - if present linear profile 563 890 IF ( zdh(ji,jj) > 0._wp ) THEN 564 zgamma_b = 6.0565 891 DO jk = imld(ji,jj)+1 , ibld(ji,jj) 566 892 zznd_pyc = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 567 893 ! 568 zdiffut(ji,jj,jk) = zdifpyc_sc(ji,jj) * EXP( zgamma_b *zznd_pyc )894 zdiffut(ji,jj,jk) = zdifpyc_sc(ji,jj) * ( 1.0 + zznd_pyc ) 569 895 ! 570 zviscos(ji,jj,jk) = zvispyc_sc(ji,jj) * EXP( zgamma_b *zznd_pyc )896 zviscos(ji,jj,jk) = zvispyc_sc(ji,jj) * ( 1.0 + zznd_pyc ) 571 897 END DO 572 IF ( ibld_ext == 0 ) THEN573 zdiffut(ji,jj,ibld(ji,jj)) = 0._wp574 zviscos(ji,jj,ibld(ji,jj)) = 0._wp575 ELSE576 zdiffut(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj) * ( hbl(ji,jj) - gdepw_n(ji, jj, ibld(ji,jj)-1) )577 zviscos(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj) * ( hbl(ji,jj) - gdepw_n(ji, jj, ibld(ji,jj)-1) )578 ENDIF579 898 ENDIF 580 ! Temporary fix to ensure zdiffut is +ve; won't be necessary with wn taken out 581 zdiffut(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj) * e3t_n(ji,jj,ibld(ji,jj)) 582 zviscos(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj) * e3t_n(ji,jj,ibld(ji,jj)) 899 ! Temporay fix to ensure zdiffut is +ve; won't be necessary with wn taken out 900 zdiffut(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj)* e3t_n(ji,jj,ibld(ji,jj)) 583 901 ! could be taken out, take account of entrainment represents as a diffusivity 584 902 ! should remove w from here, represents entrainment … … 590 908 zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * (1.0 - zznd_ml) * ( 1.0 - zznd_ml**2 ) 591 909 END DO 592 593 IF ( ibld_ext == 0 ) THEN594 zdiffut(ji,jj,ibld(ji,jj)) = 0._wp595 zviscos(ji,jj,ibld(ji,jj)) = 0._wp596 ELSE597 zdiffut(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 0._wp) * e3w_n(ji, jj, ibld(ji,jj))598 zviscos(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 0._wp) * e3w_n(ji, jj, ibld(ji,jj))599 ENDIF600 910 ENDIF ! end if ( lconv ) 601 911 ! 602 912 END DO ! end of ji loop 603 913 END DO ! end of jj loop … … 642 952 END DO ! end of jj loop 643 953 954 644 955 ! Stokes term in flux-gradient relationship (note in zsc_uw_n don't use zvstr since term needs to go to zero as zwstrl goes to zero) 645 956 WHERE ( lconv ) 646 zsc_uw_1 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke / MAX( ( 1.0 - 1.0 * 6.5 * zla**(8.0/3.0) ), 0.2)647 zsc_uw_2 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke / MIN( zla**(8.0/3.0) + epsln, 0.12)648 zsc_vw_1 = ff_t * zhml * zustke**3 * MIN( zla**(8.0/3.0), 0.12) / ( ( zvstr**3 + 0.5 * zwstrc**3 )**(2.0/3.0) + epsln )957 zsc_uw_1 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke /( 1.0 - 1.0 * 6.5 * zla**(8.0/3.0) ) 958 zsc_uw_2 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke / ( zla**(8.0/3.0) + epsln ) 959 zsc_vw_1 = ff_t * zhml * zustke**3 * zla**(8.0/3.0) / ( ( zvstr**3 + 0.5 * zwstrc**3 )**(2.0/3.0) + epsln ) 649 960 ELSEWHERE 650 961 zsc_uw_1 = zustar**2 651 zsc_vw_1 = ff_t * zhbl * zustke**3 * MIN( zla**(8.0/3.0), 0.12) / (zvstr**2 + epsln)962 zsc_vw_1 = ff_t * zhbl * zustke**3 * zla**(8.0/3.0) / (zvstr**2 + epsln) 652 963 ENDWHERE 653 IF(ln_dia_osm) THEN 654 IF ( iom_use("ghamu_00") ) CALL iom_put( "ghamu_00", wmask*ghamu ) 655 IF ( iom_use("ghamv_00") ) CALL iom_put( "ghamv_00", wmask*ghamv ) 656 END IF 964 657 965 DO jj = 2, jpjm1 658 966 DO ji = 2, jpim1 … … 699 1007 zl_l = 2.0 * ( 1.0 - EXP ( - 2.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) ) & 700 1008 & * ( 1.0 - EXP ( - 5.0 * ( 1.0 - zznd_ml ) ) ) * ( 1.0 + dstokes(ji,jj) / zhml (ji,jj) ) 701 zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0 + EXP ( -3.0 * LOG10 ( - zhol(ji,jj) ) ) ) ** (3.0/2.0)1009 zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0 + EXP ( 3.0 * LOG10 ( - zhol(ji,jj) ) ) ) ** (3.0/2.0) 702 1010 ! non-gradient buoyancy terms 703 1011 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * 0.5 * zsc_wth_1(ji,jj) * zl_eps * zhml(ji,jj) / ( 0.15 + zznd_ml ) … … 712 1020 END DO ! ji loop 713 1021 END DO ! jj loop 1022 714 1023 715 1024 WHERE ( lconv ) … … 742 1051 END DO ! jj loop 743 1052 744 IF(ln_dia_osm) THEN745 IF ( iom_use("ghamu_0") ) CALL iom_put( "ghamu_0", wmask*ghamu )746 IF ( iom_use("zsc_uw_1_0") ) CALL iom_put( "zsc_uw_1_0", tmask(:,:,1)*zsc_uw_1 )747 END IF748 1053 ! Transport term in flux-gradient relationship [note : includes ROI ratio (X0.3) ] 749 1054 … … 784 1089 END DO ! jj loop 785 1090 1091 786 1092 WHERE ( lconv ) 787 1093 zsc_uw_1 = zustar**2 … … 828 1134 END DO 829 1135 END DO 830 831 IF(ln_dia_osm) THEN832 IF ( iom_use("ghamu_f") ) CALL iom_put( "ghamu_f", wmask*ghamu )833 IF ( iom_use("ghamv_f") ) CALL iom_put( "ghamv_f", wmask*ghamv )834 IF ( iom_use("zsc_uw_1_f") ) CALL iom_put( "zsc_uw_1_f", tmask(:,:,1)*zsc_uw_1 )835 IF ( iom_use("zsc_vw_1_f") ) CALL iom_put( "zsc_vw_1_f", tmask(:,:,1)*zsc_vw_1 )836 IF ( iom_use("zsc_uw_2_f") ) CALL iom_put( "zsc_uw_2_f", tmask(:,:,1)*zsc_uw_2 )837 IF ( iom_use("zsc_vw_2_f") ) CALL iom_put( "zsc_vw_2_f", tmask(:,:,1)*zsc_vw_2 )838 END IF839 1136 ! 840 1137 ! Make surface forced velocity non-gradient terms go to zero at the base of the mixed layer. … … 868 1165 END DO 869 1166 870 IF(ln_dia_osm) THEN871 IF ( iom_use("ghamu_b") ) CALL iom_put( "ghamu_b", wmask*ghamu )872 IF ( iom_use("ghamv_b") ) CALL iom_put( "ghamv_b", wmask*ghamv )873 END IF874 1167 ! pynocline contributions 875 1168 ! Temporary fix to avoid instabilities when zdb_bl becomes very very small … … 877 1170 DO jj = 2, jpjm1 878 1171 DO ji = 2, jpim1 879 IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN 880 DO jk= 2, ibld(ji,jj) 881 znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 882 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk) 883 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk) 884 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc(ji,jj,jk) 885 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zsc_uw_1(ji,jj) * ( 1.0 - znd )**(7.0/4.0) * zdbdz_pyc(ji,jj,jk) 886 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc(ji,jj,jk) 887 END DO 888 END IF 889 END DO 1172 DO jk= 2, ibld(ji,jj) 1173 znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 1174 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk) 1175 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk) 1176 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc(ji,jj,jk) 1177 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zsc_uw_1(ji,jj) * ( 1.0 - znd )**(7.0/4.0) * zdbdz_pyc(ji,jj,jk) 1178 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc(ji,jj,jk) 1179 END DO 1180 END DO 890 1181 END DO 891 1182 … … 894 1185 DO jj=2, jpjm1 895 1186 DO ji = 2, jpim1 896 IF ( lconv(ji,jj) .AND. ibld(ji,jj) + ibld_ext < mbkt(ji,jj)) THEN1187 IF ( lconv(ji,jj) ) THEN 897 1188 DO jk = 1, imld(ji,jj) - 1 898 1189 znd=gdepw_n(ji,jj,jk) / zhml(ji,jj) 899 !ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * znd900 !ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * znd1190 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * znd 1191 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * znd 901 1192 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * znd 902 1193 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * znd … … 904 1195 DO jk = imld(ji,jj), ibld(ji,jj) 905 1196 znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 906 !ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * ( 1.0 + znd )907 !ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * ( 1.0 + znd )1197 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * ( 1.0 + znd ) 1198 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * ( 1.0 + znd ) 908 1199 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * ( 1.0 + znd ) 909 1200 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * ( 1.0 + znd ) 910 1201 END DO 911 1202 ENDIF 912 913 ghamt(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 914 ghams(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 915 ghamu(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 916 ghamv(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 1203 ghamt(ji,jj,ibld(ji,jj)) = 0._wp 1204 ghams(ji,jj,ibld(ji,jj)) = 0._wp 1205 ghamu(ji,jj,ibld(ji,jj)) = 0._wp 1206 ghamv(ji,jj,ibld(ji,jj)) = 0._wp 917 1207 END DO ! ji loop 918 1208 END DO ! jj loop 919 1209 920 IF(ln_dia_osm) THEN 921 IF ( iom_use("ghamu_1") ) CALL iom_put( "ghamu_1", wmask*ghamu ) 922 IF ( iom_use("ghamv_1") ) CALL iom_put( "ghamv_1", wmask*ghamv ) 923 IF ( iom_use("zuw_bse") ) CALL iom_put( "zuw_bse", tmask(:,:,1)*zuw_bse ) 924 IF ( iom_use("zvw_bse") ) CALL iom_put( "zvw_bse", tmask(:,:,1)*zvw_bse ) 925 IF ( iom_use("zdudz_pyc") ) CALL iom_put( "zdudz_pyc", wmask*zdudz_pyc ) 926 IF ( iom_use("zdvdz_pyc") ) CALL iom_put( "zdvdz_pyc", wmask*zdvdz_pyc ) 927 IF ( iom_use("zviscos") ) CALL iom_put( "zviscos", wmask*zviscos ) 928 END IF 1210 929 1211 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 930 1212 ! Need to put in code for contributions that are applied explicitly to … … 950 1232 IF(ln_dia_osm) THEN 951 1233 IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask*zdtdz_pyc ) 952 IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask*zdsdz_pyc )953 IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask*zdbdz_pyc )954 1234 END IF 955 1235 … … 1007 1287 1008 1288 ! Lateral boundary conditions on zvicos (sign unchanged), needed to caclulate viscosities on u and v grids 1009 !CALL lbc_lnk(zviscos(:,:,:), 'W', 1. )1289 CALL lbc_lnk( 'zdfosm', zviscos(:,:,:), 'W', 1. ) 1010 1290 1011 1291 ! GN 25/8: need to change tmask --> wmask … … 1039 1319 ! Lateral boundary conditions on final outputs for gham[uv], on [UV]-grid (sign unchanged) 1040 1320 CALL lbc_lnk_multi( 'zdfosm', ghamt, 'W', 1. , ghams, 'W', 1., & 1041 & ghamu, 'U', -1. , ghamv, 'V', -1. )1042 1043 IF(ln_dia_osm) THEN1321 & ghamu, 'U', 1. , ghamv, 'V', 1. ) 1322 1323 IF(ln_dia_osm) THEN 1044 1324 SELECT CASE (nn_osm_wave) 1045 1325 ! Stokes drift set by assumimg onstant La#=0.3(=0) or Pierson-Moskovitz spectrum (=1). … … 1050 1330 ! Stokes drift read in from sbcwave (=2). 1051 1331 CASE(2) 1052 IF ( iom_use("us_x") ) CALL iom_put( "us_x", ut0sd*umask(:,:,1) ) ! x surface Stokes drift 1053 IF ( iom_use("us_y") ) CALL iom_put( "us_y", vt0sd*vmask(:,:,1) ) ! y surface Stokes drift 1054 IF ( iom_use("wmp") ) CALL iom_put( "wmp", wmp*tmask(:,:,1) ) ! wave mean period 1055 IF ( iom_use("hsw") ) CALL iom_put( "hsw", hsw*tmask(:,:,1) ) ! significant wave height 1056 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 1057 IF ( iom_use("hsw_NP") ) CALL iom_put( "hsw_NP", (0.22/grav)*wndm**2*tmask(:,:,1) ) ! significant wave height from NP spectrum 1058 IF ( iom_use("wndm") ) CALL iom_put( "wndm", wndm*tmask(:,:,1) ) ! U_10 1332 IF ( iom_use("us_x") ) CALL iom_put( "us_x", ut0sd ) ! x surface Stokes drift 1333 IF ( iom_use("us_y") ) CALL iom_put( "us_y", vt0sd ) ! y surface Stokes drift 1059 1334 IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rau0*tmask(:,:,1)*zustar**2* & 1060 1335 & SQRT(ut0sd**2 + vt0sd**2 ) ) … … 1067 1342 IF ( iom_use("zws0") ) CALL iom_put( "zws0", tmask(:,:,1)*zws0 ) ! <Sw_0> 1068 1343 IF ( iom_use("hbl") ) CALL iom_put( "hbl", tmask(:,:,1)*hbl ) ! boundary-layer depth 1069 IF ( iom_use("ibld") ) CALL iom_put( "ibld", tmask(:,:,1)*ibld ) ! boundary-layer max k 1070 IF ( iom_use("zdt_bl") ) CALL iom_put( "zdt_bl", tmask(:,:,1)*zdt_bl ) ! dt at ml base 1071 IF ( iom_use("zds_bl") ) CALL iom_put( "zds_bl", tmask(:,:,1)*zds_bl ) ! ds at ml base 1072 IF ( iom_use("zdb_bl") ) CALL iom_put( "zdb_bl", tmask(:,:,1)*zdb_bl ) ! db at ml base 1073 IF ( iom_use("zdu_bl") ) CALL iom_put( "zdu_bl", tmask(:,:,1)*zdu_bl ) ! du at ml base 1074 IF ( iom_use("zdv_bl") ) CALL iom_put( "zdv_bl", tmask(:,:,1)*zdv_bl ) ! dv at ml base 1075 IF ( iom_use("dh") ) CALL iom_put( "dh", tmask(:,:,1)*dh ) ! Initial boundary-layer depth 1076 IF ( iom_use("hml") ) CALL iom_put( "hml", tmask(:,:,1)*hml ) ! Initial boundary-layer depth 1344 IF ( iom_use("hbli") ) CALL iom_put( "hbli", tmask(:,:,1)*hbli ) ! Initial boundary-layer depth 1077 1345 IF ( iom_use("dstokes") ) CALL iom_put( "dstokes", tmask(:,:,1)*dstokes ) ! Stokes drift penetration depth 1078 1346 IF ( iom_use("zustke") ) CALL iom_put( "zustke", tmask(:,:,1)*zustke ) ! Stokes drift magnitude at T-points … … 1080 1348 IF ( iom_use("zwstrl") ) CALL iom_put( "zwstrl", tmask(:,:,1)*zwstrl ) ! Langmuir velocity scale 1081 1349 IF ( iom_use("zustar") ) CALL iom_put( "zustar", tmask(:,:,1)*zustar ) ! friction velocity scale 1082 IF ( iom_use("zvstr") ) CALL iom_put( "zvstr", tmask(:,:,1)*zvstr ) ! mixed velocity scale1083 IF ( iom_use("zla") ) CALL iom_put( "zla", tmask(:,:,1)*zla ) ! langmuir #1084 1350 IF ( iom_use("wind_power") ) CALL iom_put( "wind_power", 1000.*rau0*tmask(:,:,1)*zustar**3 ) ! BL depth internal to zdf_osm routine 1085 1351 IF ( iom_use("wind_wave_power") ) CALL iom_put( "wind_wave_power", 1000.*rau0*tmask(:,:,1)*zustar**2*zustke ) 1086 1352 IF ( iom_use("zhbl") ) CALL iom_put( "zhbl", tmask(:,:,1)*zhbl ) ! BL depth internal to zdf_osm routine 1087 1353 IF ( iom_use("zhml") ) CALL iom_put( "zhml", tmask(:,:,1)*zhml ) ! ML depth internal to zdf_osm routine 1088 IF ( iom_use("imld") ) CALL iom_put( "imld", tmask(:,:,1)*imld ) ! index for ML depth internal to zdf_osm routine 1089 IF ( iom_use("zdh") ) CALL iom_put( "zdh", tmask(:,:,1)*zdh ) ! pyc thicknessh internal to zdf_osm routine 1354 IF ( iom_use("zdh") ) CALL iom_put( "zdh", tmask(:,:,1)*zdh ) ! ML depth internal to zdf_osm routine 1090 1355 IF ( iom_use("zhol") ) CALL iom_put( "zhol", tmask(:,:,1)*zhol ) ! ML depth internal to zdf_osm routine 1091 IF ( iom_use("zwthav") ) CALL iom_put( "zwthav", tmask(:,:,1)*zwthav ) ! upward BL-avged turb temp flux 1092 IF ( iom_use("zwth_ent") ) CALL iom_put( "zwth_ent", tmask(:,:,1)*zwth_ent ) ! upward turb temp entrainment flux 1093 IF ( iom_use("zwb_ent") ) CALL iom_put( "zwb_ent", tmask(:,:,1)*zwb_ent ) ! upward turb buoyancy entrainment flux 1094 IF ( iom_use("zws_ent") ) CALL iom_put( "zws_ent", tmask(:,:,1)*zws_ent ) ! upward turb salinity entrainment flux 1095 IF ( iom_use("zt_ml") ) CALL iom_put( "zt_ml", tmask(:,:,1)*zt_ml ) ! average T in ML 1356 IF ( iom_use("zwthav") ) CALL iom_put( "zwthav", tmask(:,:,1)*zwthav ) ! ML depth internal to zdf_osm routine 1357 IF ( iom_use("zwth_ent") ) CALL iom_put( "zwth_ent", tmask(:,:,1)*zwth_ent ) ! ML depth internal to zdf_osm routine 1358 IF ( iom_use("zt_ml") ) CALL iom_put( "zt_ml", tmask(:,:,1)*zt_ml ) ! average T in ML 1096 1359 END IF 1097 1098 CONTAINS 1099 1100 1101 ! Alan: do we need zb? 1102 SUBROUTINE zdf_osm_vertical_average( jnlev_av, zt, zs, zu, zv, zdt, zds, zdb, zdu, zdv ) 1103 !!--------------------------------------------------------------------- 1104 !! *** ROUTINE zdf_vertical_average *** 1105 !! 1106 !! ** Purpose : Determines vertical averages from surface to jnlev. 1107 !! 1108 !! ** Method : Averages are calculated from the surface to jnlev. 1109 !! The external level used to calculate differences is ibld+ibld_ext 1110 !! 1111 !!---------------------------------------------------------------------- 1112 1113 INTEGER, DIMENSION(jpi,jpj) :: jnlev_av ! Number of levels to average over. 1114 1115 ! Alan: do we need zb? 1116 REAL(wp), DIMENSION(jpi,jpj) :: zt, zs ! Average temperature and salinity 1117 REAL(wp), DIMENSION(jpi,jpj) :: zu,zv ! Average current components 1118 REAL(wp), DIMENSION(jpi,jpj) :: zdt, zds, zdb ! Difference between average and value at base of OSBL 1119 REAL(wp), DIMENSION(jpi,jpj) :: zdu, zdv ! Difference for velocity components. 1120 1121 INTEGER :: jk, ji, jj 1122 REAL(wp) :: zthick, zthermal, zbeta 1123 1124 1125 zt = 0._wp 1126 zs = 0._wp 1127 zu = 0._wp 1128 zv = 0._wp 1129 DO jj = 2, jpjm1 ! Vertical slab 1130 DO ji = 2, jpim1 1131 zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 1132 zbeta = rab_n(ji,jj,1,jp_sal) 1133 ! average over depth of boundary layer 1134 zthick = epsln 1135 DO jk = 2, jnlev_av(ji,jj) 1136 zthick = zthick + e3t_n(ji,jj,jk) 1137 zt(ji,jj) = zt(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem) 1138 zs(ji,jj) = zs(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) 1139 zu(ji,jj) = zu(ji,jj) + e3t_n(ji,jj,jk) & 1140 & * ( ub(ji,jj,jk) + ub(ji - 1,jj,jk) ) & 1141 & / MAX( 1. , umask(ji,jj,jk) + umask(ji - 1,jj,jk) ) 1142 zv(ji,jj) = zv(ji,jj) + e3t_n(ji,jj,jk) & 1143 & * ( vb(ji,jj,jk) + vb(ji,jj - 1,jk) ) & 1144 & / MAX( 1. , vmask(ji,jj,jk) + vmask(ji,jj - 1,jk) ) 1145 END DO 1146 zt(ji,jj) = zt(ji,jj) / zthick 1147 zs(ji,jj) = zs(ji,jj) / zthick 1148 zu(ji,jj) = zu(ji,jj) / zthick 1149 zv(ji,jj) = zv(ji,jj) / zthick 1150 ! Alan: do we need zb? 1151 zdt(ji,jj) = zt(ji,jj) - tsn(ji,jj,ibld(ji,jj)+ibld_ext,jp_tem) 1152 zds(ji,jj) = zs(ji,jj) - tsn(ji,jj,ibld(ji,jj)+ibld_ext,jp_sal) 1153 zdu(ji,jj) = zu(ji,jj) - ( ub(ji,jj,ibld(ji,jj)+ibld_ext) + ub(ji-1,jj,ibld(ji,jj)+ibld_ext ) ) & 1154 & / MAX(1. , umask(ji,jj,ibld(ji,jj)+ibld_ext ) + umask(ji-1,jj,ibld(ji,jj)+ibld_ext ) ) 1155 zdv(ji,jj) = zv(ji,jj) - ( vb(ji,jj,ibld(ji,jj)+ibld_ext) + vb(ji,jj-1,ibld(ji,jj)+ibld_ext ) ) & 1156 & / MAX(1. , vmask(ji,jj,ibld(ji,jj)+ibld_ext ) + vmask(ji,jj-1,ibld(ji,jj)+ibld_ext ) ) 1157 zdb(ji,jj) = grav * zthermal * zdt(ji,jj) - grav * zbeta * zds(ji,jj) 1158 END DO 1159 END DO 1160 END SUBROUTINE zdf_osm_vertical_average 1161 1162 SUBROUTINE zdf_osm_velocity_rotation( zcos_w, zsin_w, zu, zv, zdu, zdv ) 1163 !!--------------------------------------------------------------------- 1164 !! *** ROUTINE zdf_velocity_rotation *** 1165 !! 1166 !! ** Purpose : Rotates frame of reference of averaged velocity components. 1167 !! 1168 !! ** Method : The velocity components are rotated into frame specified by zcos_w and zsin_w 1169 !! 1170 !!---------------------------------------------------------------------- 1171 1172 REAL(wp), DIMENSION(jpi,jpj) :: zcos_w, zsin_w ! Cos and Sin of rotation angle 1173 REAL(wp), DIMENSION(jpi,jpj) :: zu, zv ! Components of current 1174 REAL(wp), DIMENSION(jpi,jpj) :: zdu, zdv ! Change in velocity components across pycnocline 1175 1176 INTEGER :: ji, jj 1177 REAL(wp) :: ztemp 1178 1179 DO jj = 2, jpjm1 1180 DO ji = 2, jpim1 1181 ztemp = zu(ji,jj) 1182 zu(ji,jj) = zu(ji,jj) * zcos_w(ji,jj) + zv(ji,jj) * zsin_w(ji,jj) 1183 zv(ji,jj) = zv(ji,jj) * zcos_w(ji,jj) - ztemp * zsin_w(ji,jj) 1184 ztemp = zdu(ji,jj) 1185 zdu(ji,jj) = zdu(ji,jj) * zcos_w(ji,jj) + zdv(ji,jj) * zsin_w(ji,jj) 1186 zdv(ji,jj) = zdv(ji,jj) * zsin_w(ji,jj) - ztemp * zsin_w(ji,jj) 1187 END DO 1188 END DO 1189 END SUBROUTINE zdf_osm_velocity_rotation 1190 1191 SUBROUTINE zdf_osm_external_gradients( zdtdz, zdsdz, zdbdz ) 1192 !!--------------------------------------------------------------------- 1193 !! *** ROUTINE zdf_osm_external_gradients *** 1194 !! 1195 !! ** Purpose : Calculates the gradients below the OSBL 1196 !! 1197 !! ** Method : Uses ibld and ibld_ext to determine levels to calculate the gradient. 1198 !! 1199 !!---------------------------------------------------------------------- 1200 1201 REAL(wp), DIMENSION(jpi,jpj) :: zdtdz, zdsdz, zdbdz ! External gradients of temperature, salinity and buoyancy. 1202 1203 INTEGER :: jj, ji, jkb, jkb1 1204 REAL(wp) :: zthermal, zbeta 1205 1206 1207 DO jj = 2, jpjm1 1208 DO ji = 2, jpim1 1209 IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN 1210 zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 1211 zbeta = rab_n(ji,jj,1,jp_sal) 1212 jkb = ibld(ji,jj) + ibld_ext 1213 jkb1 = MIN(jkb + 1, mbkt(ji,jj)) 1214 zdtdz(ji,jj) = - ( tsn(ji,jj,jkb1,jp_tem) - tsn(ji,jj,jkb,jp_tem ) ) & 1215 & / e3t_n(ji,jj,ibld(ji,jj)) 1216 zdsdz(ji,jj) = - ( tsn(ji,jj,jkb1,jp_sal) - tsn(ji,jj,jkb,jp_sal ) ) & 1217 & / e3t_n(ji,jj,ibld(ji,jj)) 1218 zdbdz(ji,jj) = grav * zthermal * zdtdz(ji,jj) - grav * zbeta * zdsdz(ji,jj) 1219 END IF 1220 END DO 1221 END DO 1222 END SUBROUTINE zdf_osm_external_gradients 1223 1224 SUBROUTINE zdf_osm_pycnocline_scalar_profiles( zdtdz, zdsdz, zdbdz ) 1225 1226 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdtdz, zdsdz, zdbdz ! gradients in the pycnocline 1227 1228 INTEGER :: jk, jj, ji 1229 REAL(wp) :: ztgrad, zsgrad, zbgrad 1230 REAL(wp) :: zgamma_b_nd, zgamma_c, znd 1231 REAL(wp) :: zzeta_s=0.3 1232 1233 DO jj = 2, jpjm1 1234 DO ji = 2, jpim1 1235 IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN 1236 IF ( lconv(ji,jj) ) THEN ! convective conditions 1237 IF ( zdb_bl(ji,jj) > 0._wp .AND. zdbdz_ext(ji,jj) > 0._wp ) THEN 1238 ztgrad = 0.5 * zdt_ml(ji,jj) / zdh(ji,jj) + zdtdz_ext(ji,jj) 1239 zsgrad = 0.5 * zds_ml(ji,jj) / zdh(ji,jj) + zdsdz_ext(ji,jj) 1240 zbgrad = 0.5 * zdb_ml(ji,jj) / zdh(ji,jj) + zdbdz_ext(ji,jj) 1241 zgamma_b_nd = zdbdz_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj) 1242 zgamma_c = ( 3.14159 / 4.0 ) * ( 0.5 + zgamma_b_nd ) /& 1243 & ( 1.0 - 0.25 * SQRT( 3.14159 / 6.0 ) - 2.0 * zgamma_b_nd * zzeta_s )**2 ! check 1244 DO jk = 2, ibld(ji,jj)+ibld_ext 1245 znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zdh(ji,jj) 1246 IF ( znd <= zzeta_s ) THEN 1247 zdtdz(ji,jj,jk) = zdtdz_ext(ji,jj) + 0.5 * zdt_ml(ji,jj) / zdh(ji,jj) * & 1248 & EXP( -6.0 * ( znd -zzeta_s )**2 ) 1249 zdsdz(ji,jj,jk) = zdsdz_ext(ji,jj) + 0.5 * zds_ml(ji,jj) / zdh(ji,jj) * & 1250 & EXP( -6.0 * ( znd -zzeta_s )**2 ) 1251 zdbdz(ji,jj,jk) = zdbdz_ext(ji,jj) + 0.5 * zdb_ml(ji,jj) / zdh(ji,jj) * & 1252 & EXP( -6.0 * ( znd -zzeta_s )**2 ) 1253 ELSE 1254 zdtdz(ji,jj,jk) = ztgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 ) 1255 zdsdz(ji,jj,jk) = zsgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 ) 1256 zdbdz(ji,jj,jk) = zbgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 ) 1257 ENDIF 1258 END DO 1259 ENDIF ! If condition not satisfied, no pycnocline present. Gradients have been initialised to zero, so do nothing 1260 ELSE 1261 ! stable conditions 1262 ! if pycnocline profile only defined when depth steady of increasing. 1263 IF ( zdhdt(ji,jj) >= 0.0 ) THEN ! Depth increasing, or steady. 1264 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 1265 IF ( zhol(ji,jj) >= 0.5 ) THEN ! Very stable - 'thick' pycnocline 1266 ztgrad = zdt_bl(ji,jj) / zhbl(ji,jj) 1267 zsgrad = zds_bl(ji,jj) / zhbl(ji,jj) 1268 zbgrad = zdb_bl(ji,jj) / zhbl(ji,jj) 1269 DO jk = 2, ibld(ji,jj) 1270 znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 1271 zdtdz(ji,jj,jk) = ztgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 1272 zdbdz(ji,jj,jk) = zbgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 1273 zdsdz(ji,jj,jk) = zsgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 1274 END DO 1275 ELSE ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form. 1276 ztgrad = zdt_bl(ji,jj) / zdh(ji,jj) 1277 zsgrad = zds_bl(ji,jj) / zdh(ji,jj) 1278 zbgrad = zdb_bl(ji,jj) / zdh(ji,jj) 1279 DO jk = 2, ibld(ji,jj) 1280 znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 1281 zdtdz(ji,jj,jk) = ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 1282 zdbdz(ji,jj,jk) = zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 1283 zdsdz(ji,jj,jk) = zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 1284 END DO 1285 ENDIF ! IF (zhol >=0.5) 1286 ENDIF ! IF (zdb_bl> 0.) 1287 ENDIF ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero 1288 ENDIF ! IF (lconv) 1289 END IF ! IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) 1290 END DO 1291 END DO 1292 1293 END SUBROUTINE zdf_osm_pycnocline_scalar_profiles 1294 1295 SUBROUTINE zdf_osm_pycnocline_shear_profiles( zdudz, zdvdz ) 1296 !!--------------------------------------------------------------------- 1297 !! *** ROUTINE zdf_osm_pycnocline_shear_profiles *** 1298 !! 1299 !! ** Purpose : Calculates velocity shear in the pycnocline 1300 !! 1301 !! ** Method : 1302 !! 1303 !!---------------------------------------------------------------------- 1304 1305 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdudz, zdvdz 1306 1307 INTEGER :: jk, jj, ji 1308 REAL(wp) :: zugrad, zvgrad, znd 1309 REAL(wp) :: zzeta_v = 0.45 1360 ! Lateral boundary conditions on p_avt (sign unchanged) 1361 CALL lbc_lnk( 'zdfosm', p_avt(:,:,:), 'W', 1. ) 1310 1362 ! 1311 DO jj = 2, jpjm1 1312 DO ji = 2, jpim1 1313 ! 1314 IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN 1315 IF ( lconv (ji,jj) ) THEN 1316 ! Unstable conditions 1317 zugrad = 0.7 * zdu_ml(ji,jj) / zdh(ji,jj) + 0.3 * zustar(ji,jj)*zustar(ji,jj) / & 1318 & ( ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) * & 1319 & MIN(zla(ji,jj)**(8.0/3.0) + epsln, 0.12 )) 1320 !Alan is this right? 1321 zvgrad = ( 0.7 * zdv_ml(ji,jj) + & 1322 & 2.0 * ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) / & 1323 & ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird + epsln ) & 1324 & )/ (zdh(ji,jj) + epsln ) 1325 DO jk = 2, ibld(ji,jj) - 1 + ibld_ext 1326 znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / (zdh(ji,jj) + epsln ) - zzeta_v 1327 IF ( znd <= 0.0 ) THEN 1328 zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( 3.0 * znd ) 1329 zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( 3.0 * znd ) 1330 ELSE 1331 zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( -2.0 * znd ) 1332 zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( -2.0 * znd ) 1333 ENDIF 1334 END DO 1335 ELSE 1336 ! stable conditions 1337 zugrad = 3.25 * zdu_bl(ji,jj) / zhbl(ji,jj) 1338 zvgrad = 2.75 * zdv_bl(ji,jj) / zhbl(ji,jj) 1339 DO jk = 2, ibld(ji,jj) 1340 znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 1341 IF ( znd < 1.0 ) THEN 1342 zdudz(ji,jj,jk) = zugrad * EXP( -40.0 * ( znd - 1.0 )**2 ) 1343 ELSE 1344 zdudz(ji,jj,jk) = zugrad * EXP( -20.0 * ( znd - 1.0 )**2 ) 1345 ENDIF 1346 zdvdz(ji,jj,jk) = zvgrad * EXP( -20.0 * ( znd - 0.85 )**2 ) 1347 END DO 1348 ENDIF 1349 ! 1350 END IF ! IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) 1351 END DO 1352 END DO 1353 END SUBROUTINE zdf_osm_pycnocline_shear_profiles 1354 1355 SUBROUTINE zdf_osm_calculate_dhdt( zdhdt, zdhdt_2 ) 1356 !!--------------------------------------------------------------------- 1357 !! *** ROUTINE zdf_osm_calculate_dhdt *** 1358 !! 1359 !! ** Purpose : Calculates the rate at which hbl changes. 1360 !! 1361 !! ** Method : 1362 !! 1363 !!---------------------------------------------------------------------- 1364 1365 REAL(wp), DIMENSION(jpi,jpj) :: zdhdt, zdhdt_2 ! Rate of change of hbl 1366 1367 INTEGER :: jj, ji 1368 REAL(wp) :: zgamma_b_nd, zgamma_dh_nd, zpert 1369 REAL(wp) :: zvel_max, zwb_min 1370 REAL(wp) :: zwcor, zrf_conv, zrf_shear, zrf_langmuir, zr_stokes 1371 REAL(wp) :: zzeta_m = 0.3 1372 REAL(wp) :: zgamma_c = 2.0 1373 REAL(wp) :: zdhoh = 0.1 1374 REAL(wp) :: alpha_bc = 0.5 1375 1376 DO jj = 2, jpjm1 1377 DO ji = 2, jpim1 1378 IF ( lconv(ji,jj) ) THEN ! Convective 1379 ! Alan is this right? Yes, it's a bit different from the previous relationship 1380 ! zwb_ent(ji,jj) = - 2.0 * 0.2 * zwbav(ji,jj) & 1381 ! & - ( 0.15 * ( 1.0 - EXP( -1.5 * zla(ji,jj) ) ) * zustar(ji,jj)**3 + 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj) 1382 zwcor = ABS(ff_t(ji,jj)) * zhbl(ji,jj) + epsln 1383 zrf_conv = TANH( ( zwstrc(ji,jj) / zwcor )**0.69 ) 1384 zrf_shear = TANH( ( zustar(ji,jj) / zwcor )**0.69 ) 1385 zrf_langmuir = TANH( ( zwstrl(ji,jj) / zwcor )**0.69 ) 1386 zr_stokes = 1.0 - EXP( -25.0 * dstokes(ji,jj) / hbl(ji,jj) & 1387 & * ( 1.0 + 4.0 * dstokes(ji,jj) / hbl(ji,jj) ) ) 1388 1389 zwb_ent(ji,jj) = - 2.0 * 0.2 * zrf_conv * zwbav(ji,jj) & 1390 & - 0.15 * zrf_shear * zustar(ji,jj)**3 /zhml(ji,jj) & 1391 & + zr_stokes * ( 0.15 * EXP( -1.5 * zla(ji,jj) ) * zrf_shear * zustar(ji,jj)**3 & 1392 & - zrf_langmuir * 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj) 1393 ! 1394 zwb_min = dh(ji,jj) * zwb0(ji,jj) / zhml(ji,jj) + zwb_ent(ji,jj) 1395 zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / hbl(ji,jj) ) * zwb_ent(ji,jj) / & 1396 & ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 1397 zdhdt(ji,jj) = -zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) 1398 ! added ajgn 23 July as temporay fix 1399 zdhdt_2(ji,jj) = 0._wp 1400 1401 ! commented out ajgn 23 July as temporay fix 1402 ! IF ( zdb_ml(ji,jj) > 0.0 .and. zdbdz_ext(ji,jj) > 0.0 ) THEN 1403 ! !additional term to account for thickness of pycnocline on dhdt. Small correction, so could get rid of this if necessary. 1404 ! zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 1405 ! zgamma_b_nd = zdbdz_ext(ji,jj) * zhml(ji,jj) / zdb_ml(ji,jj) 1406 ! zgamma_dh_nd = zdbdz_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj) 1407 ! zdhdt_2(ji,jj) = ( 1.0 - SQRT( 3.1415 / ( 4.0 * zgamma_c) ) * zdhoh ) * zdh(ji,jj) / zhml(ji,jj) 1408 ! zdhdt_2(ji,jj) = zdhdt_2(ji,jj) * ( zwb0(ji,jj) - (1.0 + zgamma_b_nd / alpha_bc ) * zwb_min ) 1409 ! ! Alan no idea what this should be? 1410 ! zdhdt_2(ji,jj) = alpha_bc / ( 4.0 * zgamma_c ) * zdhdt_2(ji,jj) & 1411 ! & + (alpha_bc + zgamma_dh_nd ) * ( 1.0 + SQRT( 3.1414 / ( 4.0 * zgamma_c ) ) * zdh(ji,jj) / zhbl(ji,jj) ) & 1412 ! & * (1.0 / ( 4.0 * zgamma_c * alpha_bc ) ) * zwb_min * zdh(ji,jj) / zhbl(ji,jj) 1413 ! zdhdt_2(ji,jj) = zdhdt_2(ji,jj) / ( zvel_max + MAX( zdb_bl(ji,jj), 1.0e-15 ) ) 1414 ! IF ( zdhdt_2(ji,jj) <= 0.2 * zdhdt(ji,jj) ) THEN 1415 ! zdhdt(ji,jj) = zdhdt(ji,jj) + zdhdt_2(ji,jj) 1416 ! ENDIF 1417 ELSE ! Stable 1418 zdhdt(ji,jj) = ( 0.06 + 0.52 * zhol(ji,jj) / 2.0 ) * zvstr(ji,jj)**3 / hbl(ji,jj) + zwbav(ji,jj) 1419 zdhdt_2(ji,jj) = 0._wp 1420 IF ( zdhdt(ji,jj) < 0._wp ) THEN 1421 ! For long timsteps factor in brackets slows the rapid collapse of the OSBL 1422 zpert = 2.0 * ( 1.0 + 0.0 * 2.0 * zvstr(ji,jj) * rn_rdt / hbl(ji,jj) ) * zvstr(ji,jj)**2 / hbl(ji,jj) 1423 ELSE 1424 zpert = MAX( 2.0 * ( 1.0 + 0.0 * 2.0 * zvstr(ji,jj) * rn_rdt / hbl(ji,jj) ) * zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) ) 1425 ENDIF 1426 zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / zpert 1427 ENDIF 1428 END DO 1429 END DO 1430 END SUBROUTINE zdf_osm_calculate_dhdt 1431 1432 SUBROUTINE zdf_osm_timestep_hbl( zdhdt, zdhdt_2 ) 1433 !!--------------------------------------------------------------------- 1434 !! *** ROUTINE zdf_osm_timestep_hbl *** 1435 !! 1436 !! ** Purpose : Increments hbl. 1437 !! 1438 !! ** Method : If thechange in hbl exceeds one model level the change is 1439 !! is calculated by moving down the grid, changing the buoyancy 1440 !! jump. This is to ensure that the change in hbl does not 1441 !! overshoot a stable layer. 1442 !! 1443 !!---------------------------------------------------------------------- 1444 1445 1446 REAL(wp), DIMENSION(jpi,jpj) :: zdhdt, zdhdt_2 ! rates of change of hbl. 1447 1448 INTEGER :: jk, jj, ji, jm 1449 REAL(wp) :: zhbl_s, zvel_max, zdb 1450 REAL(wp) :: zthermal, zbeta 1451 1452 DO jj = 2, jpjm1 1453 DO ji = 2, jpim1 1454 IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN 1455 ! 1456 ! If boundary layer changes by more than one level, need to check for stable layers between initial and final depths. 1457 ! 1458 zhbl_s = hbl(ji,jj) 1459 jm = imld(ji,jj) 1460 zthermal = rab_n(ji,jj,1,jp_tem) 1461 zbeta = rab_n(ji,jj,1,jp_sal) 1462 1463 1464 IF ( lconv(ji,jj) ) THEN 1465 !unstable 1466 zvel_max = -( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / hbl(ji,jj) ) * zwb_ent(ji,jj) / & 1467 & ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 1468 DO jk = imld(ji,jj), ibld(ji,jj) 1469 zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) ) & 1470 & - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ), & 1471 & 0.0 ) + zvel_max 1472 1473 zhbl_s = zhbl_s + MIN( & 1474 & rn_rdt * ( -zwb_ent(ji,jj) / zdb + zdhdt_2(ji,jj) ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), & 1475 & e3w_n(ji,jj,jm) ) 1476 zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol) 1477 1478 IF ( zhbl_s >= gdepw_n(ji,jj,jm+1) ) jm = jm + 1 1479 END DO 1480 hbl(ji,jj) = zhbl_s 1481 ibld(ji,jj) = jm 1482 ELSE 1483 ! stable 1484 DO jk = imld(ji,jj), ibld(ji,jj) 1485 zdb = MAX( & 1486 & grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) )& 1487 & - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ),& 1488 & 0.0 ) + & 1489 & 2.0 * zvstr(ji,jj)**2 / zhbl_s 1490 1491 ! Alan is thuis right? I have simply changed hbli to hbl 1492 zhol(ji,jj) = -zhbl_s / ( ( zvstr(ji,jj)**3 + epsln )/ zwbav(ji,jj) ) 1493 zdhdt(ji,jj) = -( zwbav(ji,jj) - 0.04 / 2.0 * zwstrl(ji,jj)**3 / zhbl_s - 0.15 / 2.0 * ( 1.0 - EXP( -1.5 * zla(ji,jj) ) ) * & 1494 & zustar(ji,jj)**3 / zhbl_s ) * ( 0.725 + 0.225 * EXP( -7.5 * zhol(ji,jj) ) ) 1495 zdhdt(ji,jj) = zdhdt(ji,jj) + zwbav(ji,jj) 1496 zhbl_s = zhbl_s + MIN( zdhdt(ji,jj) / zdb * rn_rdt / FLOAT( ibld(ji,jj) - imld(ji,jj) ), e3w_n(ji,jj,jm) ) 1497 1498 zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol) 1499 IF ( zhbl_s >= gdepw_n(ji,jj,jm) ) jm = jm + 1 1500 END DO 1501 ENDIF ! IF ( lconv ) 1502 hbl(ji,jj) = MAX(zhbl_s, gdepw_n(ji,jj,4) ) 1503 ibld(ji,jj) = MAX(jm, 4 ) 1504 ELSE 1505 ! change zero or one model level. 1506 hbl(ji,jj) = MAX(zhbl_t(ji,jj), gdepw_n(ji,jj,4) ) 1507 ENDIF 1508 zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj)) 1509 END DO 1510 END DO 1511 1512 END SUBROUTINE zdf_osm_timestep_hbl 1513 1514 SUBROUTINE zdf_osm_pycnocline_thickness( dh, zdh ) 1515 !!--------------------------------------------------------------------- 1516 !! *** ROUTINE zdf_osm_pycnocline_thickness *** 1517 !! 1518 !! ** Purpose : Calculates thickness of the pycnocline 1519 !! 1520 !! ** Method : The thickness is calculated from a prognostic equation 1521 !! that relaxes the pycnocine thickness to a diagnostic 1522 !! value. The time change is calculated assuming the 1523 !! thickness relaxes exponentially. This is done to deal 1524 !! with large timesteps. 1525 !! 1526 !!---------------------------------------------------------------------- 1527 1528 REAL(wp), DIMENSION(jpi,jpj) :: dh, zdh ! pycnocline thickness. 1529 ! 1530 INTEGER :: jj, ji 1531 INTEGER :: inhml 1532 REAL(wp) :: zari, ztau, zddhdt 1533 1534 1535 DO jj = 2, jpjm1 1536 DO ji = 2, jpim1 1537 IF ( lconv(ji,jj) ) THEN 1538 IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_ext(ji,jj) > 0._wp)THEN 1539 IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN ! near neutral stability 1540 zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / & 1541 (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zvstr(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 ) 1542 ELSE ! unstable 1543 zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / & 1544 (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zwstrc(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 ) 1545 ENDIF 1546 ztau = hbl(ji,jj) / (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird 1547 zddhdt = zari * hbl(ji,jj) 1548 ELSE 1549 ztau = hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 1550 zddhdt = 0.2 * hbl(ji,jj) 1551 ENDIF 1552 dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau ) + zddhdt * ( 1.0 - EXP( -rn_rdt / ztau ) ) 1553 ! Alan: this hml is never defined or used 1554 ELSE ! IF (lconv) 1555 ztau = hbl(ji,jj) / zvstr(ji,jj) 1556 IF ( zdhdt(ji,jj) >= 0.0 ) THEN ! probably shouldn't include wm here 1557 ! boundary layer deepening 1558 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 1559 ! pycnocline thickness set by stratification - use same relationship as for neutral conditions. 1560 zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) & 1561 & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01 , 0.2 ) 1562 zddhdt = MIN( zari, 0.2 ) * hbl(ji,jj) 1563 ELSE 1564 zddhdt = 0.2 * hbl(ji,jj) 1565 ENDIF 1566 ELSE ! IF(dhdt < 0) 1567 zddhdt = 0.2 * hbl(ji,jj) 1568 ENDIF ! IF (dhdt >= 0) 1569 dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau )+ zddhdt * ( 1.0 - EXP( -rn_rdt / ztau ) ) 1570 IF ( zdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zddhdt ! can be a problem with dh>hbl for rapid collapse 1571 ! Alan: this hml is never defined or used -- do we need it? 1572 ENDIF ! IF (lconv) 1573 1574 hml(ji,jj) = hbl(ji,jj) - dh(ji,jj) 1575 inhml = MAX( INT( dh(ji,jj) / e3t_n(ji,jj,ibld(ji,jj)) ) , 1 ) 1576 imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 3) 1577 zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 1578 zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 1579 END DO 1580 END DO 1581 1582 END SUBROUTINE zdf_osm_pycnocline_thickness 1583 1584 END SUBROUTINE zdf_osm 1363 END SUBROUTINE zdf_osm 1585 1364 1586 1365 … … 1599 1378 INTEGER :: ios ! local integer 1600 1379 INTEGER :: ji, jj, jk ! dummy loop indices 1601 REAL z1_t21602 1380 !! 1603 1381 NAMELIST/namzdf_osm/ ln_use_osm_la, rn_osm_la, rn_osm_dstokes, nn_ave & 1604 1382 & ,nn_osm_wave, ln_dia_osm, rn_osm_hbl0 & 1605 & ,ln_kpprimix, rn_riinfty, rn_difri, ln_convmix, rn_difconv , nn_osm_wave1383 & ,ln_kpprimix, rn_riinfty, rn_difri, ln_convmix, rn_difconv 1606 1384 !!---------------------------------------------------------------------- 1607 1385 ! 1608 1386 REWIND( numnam_ref ) ! Namelist namzdf_osm in reference namelist : Osmosis ML model 1609 1387 READ ( numnam_ref, namzdf_osm, IOSTAT = ios, ERR = 901) 1610 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_osm in reference namelist' , lwp)1388 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_osm in reference namelist' ) 1611 1389 1612 1390 REWIND( numnam_cfg ) ! Namelist namzdf_tke in configuration namelist : Turbulent Kinetic Energy 1613 1391 READ ( numnam_cfg, namzdf_osm, IOSTAT = ios, ERR = 902 ) 1614 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namzdf_osm in configuration namelist' , lwp)1392 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namzdf_osm in configuration namelist' ) 1615 1393 IF(lwm) WRITE ( numond, namzdf_osm ) 1616 1394 … … 1619 1397 WRITE(numout,*) 'zdf_osm_init : OSMOSIS Parameterisation' 1620 1398 WRITE(numout,*) '~~~~~~~~~~~~' 1621 WRITE(numout,*) ' Namelist namzdf_osm : set osmmixing parameters'1622 WRITE(numout,*) ' Use rn_osm_laln_use_osm_la = ', ln_use_osm_la1399 WRITE(numout,*) ' Namelist namzdf_osm : set tke mixing parameters' 1400 WRITE(numout,*) ' Use namelist rn_osm_la ln_use_osm_la = ', ln_use_osm_la 1623 1401 WRITE(numout,*) ' Turbulent Langmuir number rn_osm_la = ', rn_osm_la 1624 1402 WRITE(numout,*) ' Initial hbl for 1D runs rn_osm_hbl0 = ', rn_osm_hbl0 1625 WRITE(numout,*) ' Depth scale of Stokes drift 1403 WRITE(numout,*) ' Depth scale of Stokes drift rn_osm_dstokes = ', rn_osm_dstokes 1626 1404 WRITE(numout,*) ' horizontal average flag nn_ave = ', nn_ave 1627 1405 WRITE(numout,*) ' Stokes drift nn_osm_wave = ', nn_osm_wave … … 1645 1423 IF( zdf_osm_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_osm_init : unable to allocate arrays' ) 1646 1424 1647 1648 call osm_rst( nit000, 'READ' ) !* read or initialize hbl, dh, hmle 1649 1425 call osm_rst( nit000, 'READ' ) !* read or initialize hbl 1650 1426 1651 1427 IF( ln_zdfddm) THEN … … 1760 1536 REAL(wp) :: zN2_c ! local scalar 1761 1537 REAL(wp) :: rho_c = 0.01_wp !: density criterion for mixed layer depth 1762 INTEGER, DIMENSION( jpi,jpj):: imld_rst ! level of mixed-layer depth (pycnocline top)1538 INTEGER, DIMENSION(:,:), ALLOCATABLE :: imld_rst ! level of mixed-layer depth (pycnocline top) 1763 1539 !!---------------------------------------------------------------------- 1764 1540 ! … … 1775 1551 WRITE(numout,*) ' ===>>>> : wn not in restart file, set to zero initially' 1776 1552 END IF 1777 1778 1553 id1 = iom_varid( numror, 'hbl' , ldstop = .FALSE. ) 1779 id2 = iom_varid( numror, ' dh' , ldstop = .FALSE. )1554 id2 = iom_varid( numror, 'hbli' , ldstop = .FALSE. ) 1780 1555 IF( id1 > 0 .AND. id2 > 0) THEN ! 'hbl' exists; read and return 1781 1556 CALL iom_get( numror, jpdom_autoglo, 'hbl' , hbl , ldxios = lrxios ) 1782 CALL iom_get( numror, jpdom_autoglo, ' dh', dh, ldxios = lrxios )1783 WRITE(numout,*) ' ===>>>> : hbl & dhread from restart file'1557 CALL iom_get( numror, jpdom_autoglo, 'hbli', hbli, ldxios = lrxios ) 1558 WRITE(numout,*) ' ===>>>> : hbl & hbli read from restart file' 1784 1559 RETURN 1785 ELSE ! 'hbl' & ' dh' not in restart file, recalculate1560 ELSE ! 'hbl' & 'hbli' not in restart file, recalculate 1786 1561 WRITE(numout,*) ' ===>>>> : previous run without osmosis scheme, hbl computed from stratification' 1787 1562 END IF … … 1795 1570 CALL iom_rstput( kt, nitrst, numrow, 'wn' , wn , ldxios = lwxios ) 1796 1571 CALL iom_rstput( kt, nitrst, numrow, 'hbl' , hbl , ldxios = lwxios ) 1797 CALL iom_rstput( kt, nitrst, numrow, ' dh' , dh, ldxios = lwxios )1572 CALL iom_rstput( kt, nitrst, numrow, 'hbli' , hbli, ldxios = lwxios ) 1798 1573 RETURN 1799 1574 END IF … … 1803 1578 !!----------------------------------------------------------------------------- 1804 1579 IF( lwp ) WRITE(numout,*) ' ===>>>> : calculating hbl computed from stratification' 1580 ALLOCATE( imld_rst(jpi,jpj) ) 1805 1581 ! w-level of the mixing and mixed layers 1806 1582 CALL eos_rab( tsn, rab_n ) … … 1823 1599 DO jj = 1, jpj 1824 1600 DO ji = 1, jpi 1825 iiki = MAX(4,imld_rst(ji,jj)) 1826 hbl (ji,jj) = gdepw_n(ji,jj,iiki ) ! Turbocline depth 1827 dh (ji,jj) = e3t_n(ji,jj,iiki-1 ) ! Turbocline depth 1601 iiki = imld_rst(ji,jj) 1602 hbl (ji,jj) = gdepw_n(ji,jj,iiki ) * ssmask(ji,jj) ! Turbocline depth 1828 1603 END DO 1829 1604 END DO 1605 hbl = MAX(hbl,epsln) 1606 hbli(:,:) = hbl(:,:) 1607 DEALLOCATE( imld_rst ) 1830 1608 WRITE(numout,*) ' ===>>>> : hbl computed from stratification' 1831 wn(:,:,:) = 0._wp1832 WRITE(numout,*) ' ===>>>> : wn not in restart file, set to zero initially'1833 1609 END SUBROUTINE osm_rst 1834 1610 … … 1858 1634 ENDIF 1859 1635 1636 ! add non-local temperature and salinity flux 1860 1637 DO jk = 1, jpkm1 1861 1638 DO jj = 2, jpjm1 … … 1871 1648 END DO 1872 1649 1873 ! save the non-local tracer flux trends for diagnostics 1650 1651 ! save the non-local tracer flux trends for diagnostic 1874 1652 IF( l_trdtra ) THEN 1875 1653 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 1876 1654 ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 1877 1878 CALL trd_tra( kt, 'TRA', jp_tem, jptra_ osm, ztrdt )1879 CALL trd_tra( kt, 'TRA', jp_sal, jptra_ osm, ztrds )1655 !!bug gm jpttdzdf ==> jpttosm 1656 CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdf, ztrdt ) 1657 CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdf, ztrds ) 1880 1658 DEALLOCATE( ztrdt ) ; DEALLOCATE( ztrds ) 1881 1659 ENDIF … … 1945 1723 1946 1724 !!====================================================================== 1947 1948 1725 END MODULE zdfosm
Note: See TracChangeset
for help on using the changeset viewer.
