Changeset 14758
- Timestamp:
- 2021-04-27T21:57:06+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
r14750 r14758 83 83 PUBLIC ln_osm_mle ! logical needed by tra_mle_init in tramle.F90 84 84 85 ! Scales 86 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: swth0 ! Surface heat flux (Kinematic) 87 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: sws0 ! Surface freshwater flux 88 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: swb0 ! Surface buoyancy flux 89 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: suw0 ! Surface u-momentum flux 90 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: sustar ! Friction velocity 91 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: scos_wind ! Cos angle of surface stress 92 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ssin_wind ! Sin angle of surface stress 93 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: swthav ! Heat flux - bl average 94 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: swsav ! Freshwater flux - bl average 95 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: swbav ! Buoyancy flux - bl average 96 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: sustke ! Surface Stokes drift 97 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: dstokes ! Penetration depth of the Stokes drift 98 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: swstrl ! Langmuir velocity scale 99 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: swstrc ! Convective velocity scale 100 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: sla ! Trubulent Langmuir number 101 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: svstr ! Velocity scale that tends to sustar for large Langmuir number 102 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: shol ! Stability parameter for boundary layer 103 85 104 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ghamu !: non-local u-momentum flux 86 105 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ghamv !: non-local v-momentum flux … … 91 110 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: dh ! depth of pycnocline 92 111 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hml ! ML depth 93 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: dstokes !: penetration depth of the Stokes drift.94 112 95 113 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_ft ! inverse of the modified Coriolis parameter at t-pts … … 179 197 !! *** FUNCTION zdf_osm_alloc *** 180 198 !!---------------------------------------------------------------------- 199 ! 200 INTEGER :: ierr 201 ! 202 ALLOCATE( swth0(jpi,jpj), sws0(jpi,jpj), swb0(jpi,jpj), suw0(jpi,jpj), & 203 & sustar(jpi,jpj), scos_wind(jpi,jpj), ssin_wind(jpi,jpj), swthav(jpi,jpj), & 204 & swsav(jpi,jpj), swbav(jpi,jpj), sustke(jpi,jpj), swstrl(jpi,jpj), & 205 & swstrc(jpi,jpj), sla(jpi,jpj), svstr(jpi,jpj), shol(jpi,jpj), & 206 & STAT=ierr ) 207 zdf_osm_alloc = ierr 208 ! 181 209 ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk),ghams(jpi,jpj,jpk), & 182 210 & hbl(jpi,jpj), dh(jpi,jpj), hml(jpi,jpj), dstokes(jpi, jpj), & 183 & etmean(jpi,jpj,jpk), STAT= zdf_osm_alloc ) 184 211 & etmean(jpi,jpj,jpk), STAT=ierr ) 212 zdf_osm_alloc = zdf_osm_alloc + ierr 213 ! 185 214 ALLOCATE( hmle(jpi,jpj), r1_ft(jpi,jpj), dbdx_mle(jpi,jpj), dbdy_mle(jpi,jpj), & 186 & mld_prof(jpi,jpj), STAT= zdf_osm_alloc ) 187 215 & mld_prof(jpi,jpj), STAT=ierr ) 216 zdf_osm_alloc = zdf_osm_alloc + ierr 217 ! 188 218 CALL mpp_sum ( 'zdfosm', zdf_osm_alloc ) 189 219 IF( zdf_osm_alloc /= 0 ) CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm arrays') 190 220 ! 191 221 END FUNCTION zdf_osm_alloc 192 222 … … 249 279 REAL(wp) :: zt,zs,zu,zv,zrh ! variables used in constructing averages 250 280 ! Scales 251 REAL(wp), DIMENSION(jpi,jpj) :: zrad0 ! Surface solar temperature flux (deg m/s) 252 REAL(wp), DIMENSION(jpi,jpj) :: zradh ! Radiative flux at bl base (Buoyancy units) 253 REAL(wp) :: zradav ! Radiative flux, bl average (Buoyancy Units) 254 REAL(wp), DIMENSION(jpi,jpj) :: zustar ! friction velocity 255 REAL(wp), DIMENSION(jpi,jpj) :: zwstrl ! Langmuir velocity scale 256 REAL(wp), DIMENSION(jpi,jpj) :: zvstr ! Velocity scale that ends to zustar for large Langmuir number. 257 REAL(wp), DIMENSION(jpi,jpj) :: zwstrc ! Convective velocity scale 258 REAL(wp), DIMENSION(jpi,jpj) :: zuw0 ! Surface u-momentum flux 259 REAL(wp) :: zvw0 ! Surface v-momentum flux 260 REAL(wp), DIMENSION(jpi,jpj) :: zwth0 ! Surface heat flux (Kinematic) 261 REAL(wp), DIMENSION(jpi,jpj) :: zws0 ! Surface freshwater flux 262 REAL(wp), DIMENSION(jpi,jpj) :: zwb0 ! Surface buoyancy flux 263 REAL(wp), DIMENSION(jpi,jpj) :: zwb0tot ! Total surface buoyancy flux including insolation 264 REAL(wp), DIMENSION(jpi,jpj) :: zwthav ! Heat flux - bl average 265 REAL(wp), DIMENSION(jpi,jpj) :: zwsav ! freshwater flux - bl average 266 REAL(wp), DIMENSION(jpi,jpj) :: zwbav ! Buoyancy flux - bl average 281 REAL(wp), DIMENSION(A2D(0)) :: zrad0 ! Surface solar temperature flux (deg m/s) 282 REAL(wp), DIMENSION(A2D(0)) :: zradh ! Radiative flux at bl base (Buoyancy units) 283 REAL(wp) :: zradav ! Radiative flux, bl average (Buoyancy Units) 284 REAL(wp) :: zvw0 ! Surface v-momentum flux 285 REAL(wp), DIMENSION(A2D(0)) :: zwb0tot ! Total surface buoyancy flux including insolation 267 286 REAL(wp), DIMENSION(jpi,jpj) :: zwb_ent ! Buoyancy entrainment flux 268 287 REAL(wp), DIMENSION(jpi,jpj) :: zwb_min … … 274 293 REAL(wp), DIMENSION(jpi,jpj) :: zvel_mle ! velocity scale for dhdt with stable ML and FK 275 294 276 REAL(wp), DIMENSION(jpi,jpj) :: zustke ! Surface Stokes drift277 REAL(wp), DIMENSION(jpi,jpj) :: zla ! Trubulent Langmuir number278 REAL(wp), DIMENSION(jpi,jpj) :: zcos_wind ! Cos angle of surface stress279 REAL(wp), DIMENSION(jpi,jpj) :: zsin_wind ! Sin angle of surface stress280 REAL(wp), DIMENSION(jpi,jpj) :: zhol ! Stability parameter for boundary layer281 295 LOGICAL, DIMENSION(jpi,jpj) :: lconv ! unstable/stable bl 282 296 LOGICAL, DIMENSION(jpi,jpj) :: lshear ! Shear layers … … 349 363 IF( ln_timing ) CALL timing_start('zdf_osm') 350 364 ibld(:,:) = 0 ; imld(:,:) = 0 351 zrad0(:,:) = zlarge ; zradh(:,:) = zlarge ; zustar(:,:)= zlarge352 zwstrl(:,:) = zlarge ; zvstr(:,:) = zlarge ; zwstrc(:,:) = zlarge ; zuw0(:,:) = zlarge353 zwth0(:,:) = zlarge ; zws0(:,:) = zlarge ; zwb0(:,:) = zlarge354 zwthav(:,:) = zlarge ; zwsav(:,:) = zlarge ; zwbav(:,:) = zlarge ; zwb_ent(:,:) = zlarge355 zustke(:,:) = zlarge ; zla(:,:) = zlarge ; zcos_wind(:,:) = zlarge ; zsin_wind(:,:) = zlarge356 zhol(:,:) = zlarge ; zwb0tot(:,:)= zlarge ; zalpha_pyc(:,:) = zlarge365 sustar(:,:) = zlarge 366 swstrl(:,:) = zlarge ; svstr(:,:) = zlarge ; swstrc(:,:) = zlarge ; suw0(:,:) = zlarge 367 swth0(:,:) = zlarge ; sws0(:,:) = zlarge ; swb0(:,:) = zlarge 368 swthav(:,:) = zlarge ; swsav(:,:) = zlarge ; swbav(:,:) = zlarge ; zwb_ent(:,:) = zlarge 369 sustke(:,:) = zlarge ; sla(:,:) = zlarge ; scos_wind(:,:) = zlarge ; ssin_wind(:,:) = zlarge 370 shol(:,:) = zlarge ; zalpha_pyc(:,:) = zlarge 357 371 lconv(:,:) = .FALSE.; lpyc(:,:) = .FALSE. ; lflux(:,:) = .FALSE. ; lmle(:,:) = .FALSE. 358 372 ! mixed layer … … 428 442 & ( zz0 * ( 1.0_wp - EXP( -hbl(ji,jj)/rn_si0 ) ) * rn_si0 + & 429 443 & zz1 * ( 1.0_wp - EXP( -hbl(ji,jj)/rn_si1 ) ) * rn_si1 ) / hbl(ji,jj) 430 zwth0(ji,jj) = - qns(ji,jj) * r1_rho0_rcp * tmask(ji,jj,1) ! Upwards surface Temperature flux for non-local term431 zwthav(ji,jj) = 0.5_wp * zwth0(ji,jj) - & ! Turbulent heat flux averaged over depth of OSBL444 swth0(ji,jj) = - qns(ji,jj) * r1_rho0_rcp * tmask(ji,jj,1) ! Upwards surface Temperature flux for non-local term 445 swthav(ji,jj) = 0.5_wp * swth0(ji,jj) - & ! Turbulent heat flux averaged over depth of OSBL 432 446 & ( 0.5_wp * ( zrad0(ji,jj) + zradh(ji,jj) ) - zradav ) 433 447 END_2D 434 448 DO_2D( 0, 0, 0, 0 ) 435 zws0(ji,jj) = -1.0_wp * & ! Upwards surface salinity flux for non-local term449 sws0(ji,jj) = -1.0_wp * & ! Upwards surface salinity flux for non-local term 436 450 & ( ( emp(ji,jj) - rnf(ji,jj) ) * ts(ji,jj,1,jp_sal,Kmm) + sfx(ji,jj) ) * r1_rho0 * tmask(ji,jj,1) 437 451 zthermal = rab_n(ji,jj,1,jp_tem) 438 452 zbeta = rab_n(ji,jj,1,jp_sal) 439 zwb0(ji,jj) = grav * zthermal * zwth0(ji,jj) - & ! Non radiative upwards surface buoyancy flux440 & grav * zbeta * zws0(ji,jj)441 zwb0tot(ji,jj) = zwb0(ji,jj) - grav * zthermal * & ! Total upwards surface buoyancy flux453 swb0(ji,jj) = grav * zthermal * swth0(ji,jj) - & ! Non radiative upwards surface buoyancy flux 454 & grav * zbeta * sws0(ji,jj) 455 zwb0tot(ji,jj) = swb0(ji,jj) - grav * zthermal * & ! Total upwards surface buoyancy flux 442 456 & ( zrad0(ji,jj) - zradh(ji,jj) ) 443 zwsav(ji,jj) = 0.5 * zws0(ji,jj) ! Turbulent salinity flux averaged over depth of the OBSL444 zwbav(ji,jj) = grav * zthermal * zwthav(ji,jj) - & ! Turbulent buoyancy flux averaged over the depth of the445 & grav * zbeta * zwsav(ji,jj) ! OBSBL457 swsav(ji,jj) = 0.5 * sws0(ji,jj) ! Turbulent salinity flux averaged over depth of the OBSL 458 swbav(ji,jj) = grav * zthermal * swthav(ji,jj) - & ! Turbulent buoyancy flux averaged over the depth of the 459 & grav * zbeta * swsav(ji,jj) ! OBSBL 446 460 END_2D 447 461 DO_2D( 0, 0, 0, 0 ) 448 zuw0(ji,jj) = - 0.5 * (utau(ji-1,jj) + utau(ji,jj)) * & ! Surface upward velocity fluxes462 suw0(ji,jj) = - 0.5 * (utau(ji-1,jj) + utau(ji,jj)) * & ! Surface upward velocity fluxes 449 463 & r1_rho0 * tmask(ji,jj,1) 450 464 zvw0 = - 0.5 * (vtau(ji,jj-1) + vtau(ji,jj)) * r1_rho0 * tmask(ji,jj,1) 451 zustar(ji,jj) = MAX( SQRT( SQRT( zuw0(ji,jj) * & ! Friction velocity (zustar), at T-point : LMD94 eq. 2452 & zuw0(ji,jj) + zvw0 * zvw0 ) ), 1.0e-8_wp )453 zcos_wind(ji,jj) = -zuw0(ji,jj) / ( zustar(ji,jj) * zustar(ji,jj) )454 zsin_wind(ji,jj) = -zvw0 / ( zustar(ji,jj) * zustar(ji,jj) )465 sustar(ji,jj) = MAX( SQRT( SQRT( suw0(ji,jj) * & ! Friction velocity (sustar), at T-point : LMD94 eq. 2 466 & suw0(ji,jj) + zvw0 * zvw0 ) ), 1.0e-8_wp ) 467 scos_wind(ji,jj) = -suw0(ji,jj) / ( sustar(ji,jj) * sustar(ji,jj) ) 468 ssin_wind(ji,jj) = -zvw0 / ( sustar(ji,jj) * sustar(ji,jj) ) 455 469 #ifdef key_osm_debug 456 470 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN … … 462 476 & 'after calculating fluxes: hbl=', hbl(ji,jj),' zthermal=',zthermal, ' zbeta=', zbeta,& 463 477 & ' zrad0=', zrad0(ji,jj),' zradh=', zradh(ji,jj), ' zradav=', zradav, & 464 & ' zwth0=', zwth0(ji,jj), ' zwthav=', zwthav(ji,jj), ' zws0=', zws0(ji,jj), &465 & ' zwb0=', zwb0(ji,jj), ' zwb0tot=', zwb0tot(ji,jj), ' zwb0tot_in hbl=', zwb0tot(ji,jj) + grav * zthermal * zradh(ji,jj),&466 & ' zwbav=', zwbav(ji,jj)478 & ' swth0=', swth0(ji,jj), ' swthav=', swthav(ji,jj), ' sws0=', sws0(ji,jj), & 479 & ' swb0=', swb0(ji,jj), ' zwb0tot=', zwb0tot(ji,jj), ' zwb0tot_in hbl=', zwb0tot(ji,jj) + grav * zthermal * zradh(ji,jj),& 480 & ' swbav=', swbav(ji,jj) 467 481 FLUSH(narea+100) 468 482 END IF 469 483 #endif 470 484 END_2D 471 ! Calculate Stokes drift in direction of wind ( zustke) and Stokes penetration depth (dstokes)485 ! Calculate Stokes drift in direction of wind (sustke) and Stokes penetration depth (dstokes) 472 486 SELECT CASE (nn_osm_wave) 473 487 ! Assume constant La#=0.3 474 488 CASE(0) 475 489 DO_2D( 0, 0, 0, 0 ) 476 zus_x = zcos_wind(ji,jj) * zustar(ji,jj) / 0.3**2477 zus_y = zsin_wind(ji,jj) * zustar(ji,jj) / 0.3**2490 zus_x = scos_wind(ji,jj) * sustar(ji,jj) / 0.3**2 491 zus_y = ssin_wind(ji,jj) * sustar(ji,jj) / 0.3**2 478 492 ! Linearly 479 zustke(ji,jj) = MAX ( SQRT( zus_x*zus_x + zus_y*zus_y), 1.0e-8 )493 sustke(ji,jj) = MAX ( SQRT( zus_x*zus_x + zus_y*zus_y), 1.0e-8 ) 480 494 dstokes(ji,jj) = rn_osm_dstokes 481 495 END_2D … … 484 498 DO_2D( 0, 0, 0, 0 ) 485 499 ! Use wind speed wndm included in sbc_oce module 486 zustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8 )500 sustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8 ) 487 501 dstokes(ji,jj) = MAX ( 0.12 * wndm(ji,jj)**2 / grav, 5.e-1) 488 502 END_2D … … 495 509 ! Use wave fields 496 510 zabsstke = SQRT(ut0sd(ji,jj)**2 + vt0sd(ji,jj)**2) 497 zustke(ji,jj) = MAX ( ( zcos_wind(ji,jj) * ut0sd(ji,jj) + zsin_wind(ji,jj) * vt0sd(ji,jj) ), 1.0e-8)511 sustke(ji,jj) = MAX ( ( scos_wind(ji,jj) * ut0sd(ji,jj) + ssin_wind(ji,jj) * vt0sd(ji,jj) ), 1.0e-8) 498 512 dstokes(ji,jj) = MAX (zfac * hsw(ji,jj)*hsw(ji,jj) / ( MAX(zabsstke * wmp(ji,jj), 1.0e-7 ) ), 5.0e-1) 499 513 ELSE 500 514 ! Assume masking issue (e.g. ice in ECMWF reanalysis but not in model run) 501 515 ! .. so default to Pierson-Moskowitz 502 zustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8 )516 sustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8 ) 503 517 dstokes(ji,jj) = MAX ( 0.12 * wndm(ji,jj)**2 / grav, 5.e-1) 504 518 END IF … … 508 522 IF(narea==nn_narea_db)THEN 509 523 WRITE(narea+100,'(2(a,g11.3))') & 510 & 'Before reduction: zustke=', zustke(iloc_db,jloc_db),' dstokes =',dstokes(iloc_db,jloc_db)524 & 'Before reduction: sustke=', sustke(iloc_db,jloc_db),' dstokes =',dstokes(iloc_db,jloc_db) 511 525 FLUSH(narea+100) 512 526 END IF … … 516 530 ! Reduce both Stokes drift and its depth scale by ocean fraction to represent sheltering by ice 517 531 DO_2D( 0, 0, 0, 0 ) 518 zustke(ji,jj) = zustke(ji,jj) * (1.0_wp - fr_i(ji,jj))532 sustke(ji,jj) = sustke(ji,jj) * (1.0_wp - fr_i(ji,jj)) 519 533 dstokes(ji,jj) = dstokes(ji,jj) * (1.0_wp - fr_i(ji,jj)) 520 534 END_2D … … 529 543 ! around the mean wind. The effect of this adjustment needs to be tested. 530 544 IF(nn_osm_wave > 0) THEN 531 zustke(2:jpim1,2:jpjm1) = rn_zdfosm_adjust_sd * zustke(2:jpim1,2:jpjm1)545 sustke(A2D(0)) = rn_zdfosm_adjust_sd * sustke(A2D(0)) 532 546 END IF 533 547 CASE(1) … … 542 556 zsqrt_depth = SQRT(z2k_times_thickness) 543 557 zexp_depth = EXP(-z2k_times_thickness) 544 zustke(ji,jj) = zustke(ji,jj) * (1.0_wp - zexp_depth &558 sustke(ji,jj) = sustke(ji,jj) * (1.0_wp - zexp_depth & 545 559 & - z_two_thirds * ( zsqrtpi*zsqrt_depth*z2k_times_thickness * ERFC(zsqrt_depth) & 546 560 & + 1.0_wp - (1.0_wp + z2k_times_thickness)*zexp_depth ) ) / z2k_times_thickness … … 567 581 zf = z2k_times_thickness*(1.0_wp/zexperfc - 1.0_wp) 568 582 dstokes(ji,jj) = 5.97 * zf * dstokes(ji,jj) 569 zustke(ji,jj) = zustke(ji,jj) * EXP(z2k_times_thickness * ( 1.0_wp / (2. * zf) - 1.0_wp )) * ( 1.0_wp - zexperfc)583 sustke(ji,jj) = sustke(ji,jj) * EXP(z2k_times_thickness * ( 1.0_wp / (2. * zf) - 1.0_wp )) * ( 1.0_wp - zexperfc) 570 584 END_2D 571 585 END SELECT 572 586 573 ! Langmuir velocity scale ( zwstrl), La # (zla)574 ! mixed scale ( zvstr), convective velocity scale (zwstrc)587 ! Langmuir velocity scale (swstrl), La # (sla) 588 ! mixed scale (svstr), convective velocity scale (swstrc) 575 589 DO_2D( 0, 0, 0, 0 ) 576 ! Langmuir velocity scale ( zwstrl), at T-point577 zwstrl(ji,jj) = ( zustar(ji,jj) * zustar(ji,jj) * zustke(ji,jj) )**pthird578 zla(ji,jj) = MAX(MIN(SQRT ( zustar(ji,jj) / ( zwstrl(ji,jj) + epsln ) )**3, 4.0), 0.2)579 IF( zla(ji,jj) > 0.45) dstokes(ji,jj) = MIN(dstokes(ji,jj), 0.5_wp*hbl(ji,jj))580 ! Velocity scale that tends to zustar for large Langmuir numbers581 zvstr(ji,jj) = ( zwstrl(ji,jj)**3 + &582 & ( 1.0 - EXP( -0.5 * zla(ji,jj)**2 ) ) * zustar(ji,jj) * zustar(ji,jj) * zustar(ji,jj) )**pthird590 ! Langmuir velocity scale (swstrl), at T-point 591 swstrl(ji,jj) = ( sustar(ji,jj) * sustar(ji,jj) * sustke(ji,jj) )**pthird 592 sla(ji,jj) = MAX(MIN(SQRT ( sustar(ji,jj) / ( swstrl(ji,jj) + epsln ) )**3, 4.0), 0.2) 593 IF(sla(ji,jj) > 0.45) dstokes(ji,jj) = MIN(dstokes(ji,jj), 0.5_wp*hbl(ji,jj)) 594 ! Velocity scale that tends to sustar for large Langmuir numbers 595 svstr(ji,jj) = ( swstrl(ji,jj)**3 + & 596 & ( 1.0 - EXP( -0.5 * sla(ji,jj)**2 ) ) * sustar(ji,jj) * sustar(ji,jj) * sustar(ji,jj) )**pthird 583 597 584 598 ! limit maximum value of Langmuir number as approximate treatment for shear turbulence. 585 ! Note zustke and zwstrl are not amended.599 ! Note sustke and swstrl are not amended. 586 600 ! 587 ! get convective velocity ( zwstrc), stabilty scale (zhol) and logical conection flag lconv588 IF ( zwbav(ji,jj) > 0.0) THEN589 zwstrc(ji,jj) = ( 2.0 * zwbav(ji,jj) * 0.9 * hbl(ji,jj) )**pthird590 zhol(ji,jj) = -0.9 * hbl(ji,jj) * 2.0 * zwbav(ji,jj) / (zvstr(ji,jj)**3 + epsln )601 ! get convective velocity (swstrc), stabilty scale (shol) and logical conection flag lconv 602 IF ( swbav(ji,jj) > 0.0) THEN 603 swstrc(ji,jj) = ( 2.0 * swbav(ji,jj) * 0.9 * hbl(ji,jj) )**pthird 604 shol(ji,jj) = -0.9 * hbl(ji,jj) * 2.0 * swbav(ji,jj) / (svstr(ji,jj)**3 + epsln ) 591 605 ELSE 592 zwstrc(ji,jj) = 0.0_wp593 zhol(ji,jj) = -hbl(ji,jj) * 2.0 * zwbav(ji,jj)/ (zvstr(ji,jj)**3 + epsln )606 swstrc(ji,jj) = 0.0_wp 607 shol(ji,jj) = -hbl(ji,jj) * 2.0 * swbav(ji,jj)/ (svstr(ji,jj)**3 + epsln ) 594 608 ENDIF 595 609 #ifdef key_osm_debug 596 610 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 597 611 WRITE(narea+100,'(2(a,g11.3),/,3(a,g11.3),/,3(a,g11.3),/)') & 598 & 'After reduction: zustke=', zustke(ji,jj), ' dstokes=', dstokes(ji,jj), &599 & ' zustar =', zustar(ji,jj), ' zwstrl=', zwstrl(ji,jj), ' zwstrc=', zwstrc(ji,jj),&600 & ' zhol=', zhol(ji,jj), ' zla=', zla(ji,jj), ' zvstr=', zvstr(ji,jj)612 & 'After reduction: sustke=', sustke(ji,jj), ' dstokes=', dstokes(ji,jj), & 613 & ' zustar =', sustar(ji,jj), ' swstrl=', swstrl(ji,jj), ' swstrc=', swstrc(ji,jj),& 614 & ' shol=', shol(ji,jj), ' sla=', sla(ji,jj), ' svstr=', svstr(ji,jj) 601 615 FLUSH(narea+100) 602 616 END IF … … 675 689 #endif 676 690 ! Velocity components in frame aligned with surface stress 677 CALL zdf_osm_velocity_rotation( zu_ml, zv_ml , zcos_wind, zsin_wind)678 CALL zdf_osm_velocity_rotation( zdu_ml, zdv_ml , zcos_wind, zsin_wind)679 CALL zdf_osm_velocity_rotation( zu_bl, zv_bl , zcos_wind, zsin_wind)680 CALL zdf_osm_velocity_rotation( zdu_bl, zdv_bl , zcos_wind, zsin_wind)691 CALL zdf_osm_velocity_rotation( zu_ml, zv_ml ) 692 CALL zdf_osm_velocity_rotation( zdu_ml, zdv_ml ) 693 CALL zdf_osm_velocity_rotation( zu_bl, zv_bl ) 694 CALL zdf_osm_velocity_rotation( zdu_bl, zdv_bl ) 681 695 #ifdef key_osm_debug 682 696 IF(narea==nn_narea_db) THEN … … 911 925 #endif 912 926 ! Rotate mean currents and changes onto wind aligned co-ordinates 913 CALL zdf_osm_velocity_rotation( zu_ml, zv_ml , zcos_wind, zsin_wind)914 CALL zdf_osm_velocity_rotation( zdu_ml, zdv_ml , zcos_wind, zsin_wind)915 CALL zdf_osm_velocity_rotation( zu_bl, zv_bl , zcos_wind, zsin_wind)916 CALL zdf_osm_velocity_rotation( zdu_bl, zdv_bl , zcos_wind, zsin_wind)927 CALL zdf_osm_velocity_rotation( zu_ml, zv_ml ) 928 CALL zdf_osm_velocity_rotation( zdu_ml, zdv_ml ) 929 CALL zdf_osm_velocity_rotation( zu_bl, zv_bl ) 930 CALL zdf_osm_velocity_rotation( zdu_bl, zdv_bl ) 917 931 #ifdef key_osm_debug 918 932 IF(narea==nn_narea_db) THEN … … 966 980 ! Calculate non-gradient components of the flux-gradient relationships 967 981 ! -------------------------------------------------------------------- 968 CALL zdf_osm_fgr_terms( Kmm, ibld, imld, jp_ext, lconv, lpyc, j_ddh, zhbl, zhml, zdh, zdhdt, zhol, zshear, & 969 & zustar, zwstrl, zvstr, zwstrc, zuw0, zwth0, zws0, zwb0, zwthav, zwsav, zwbav, zustke, zla, & 970 & zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml, & 982 CALL zdf_osm_fgr_terms( Kmm, ibld, imld, jp_ext, lconv, lpyc, j_ddh, zhbl, zhml, zdh, zdhdt, zshear, & 983 & zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml, & 971 984 & zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext, zdbdz_pyc, zalpha_pyc, zdiffut, zviscos ) 972 985 … … 979 992 980 993 ! Rotate non-gradient velocity terms back to model reference frame 981 CALL zdf_osm_velocity_rotation( ghamu, ghamv, zcos_wind, zsin_wind,.FALSE., 2, ibld )994 CALL zdf_osm_velocity_rotation( ghamu, ghamv, .FALSE., 2, ibld ) 982 995 983 996 ! KPP-style Ri# mixing … … 1042 1055 DO jk = 1, ibld(ji,jj) 1043 1056 znd = gdepw(ji,jj,jk,Kmm) / MAX(zhbl(ji,jj),epsln) 1044 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) - ( zwth0(ji,jj) - zrad0(ji,jj) + zradh(ji,jj) ) * ( 1.0 - znd )1045 ghams(ji,jj,jk) = ghams(ji,jj,jk) - zws0(ji,jj) * ( 1.0 - znd )1057 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) - ( swth0(ji,jj) - zrad0(ji,jj) + zradh(ji,jj) ) * ( 1.0 - znd ) 1058 ghams(ji,jj,jk) = ghams(ji,jj,jk) - sws0(ji,jj) * ( 1.0 - znd ) 1046 1059 END DO 1047 1060 DO jk = 1, mld_prof(ji,jj) 1048 1061 znd = gdepw(ji,jj,jk,Kmm) / MAX(zhmle(ji,jj),epsln) 1049 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + ( zwth0(ji,jj) - zrad0(ji,jj) + zradh(ji,jj) ) * ( 1.0 - znd )1050 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws0(ji,jj) * ( 1.0 -znd )1062 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + ( swth0(ji,jj) - zrad0(ji,jj) + zradh(ji,jj) ) * ( 1.0 - znd ) 1063 ghams(ji,jj,jk) = ghams(ji,jj,jk) + sws0(ji,jj) * ( 1.0 -znd ) 1051 1064 END DO 1052 1065 ! Viscosity for MLEs … … 1126 1139 ! Stokes drift set by assumimg onstant La#=0.3(=0) or Pierson-Moskovitz spectrum (=1). 1127 1140 CASE(0:1) 1128 IF ( iom_use("us_x") ) CALL iom_put( "us_x", tmask(:,:,1)* zustke*zcos_wind ) ! x surface Stokes drift1129 IF ( iom_use("us_y") ) CALL iom_put( "us_y", tmask(:,:,1)* zustke*zsin_wind ) ! y surface Stokes drift1130 IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rho0*tmask(:,:,1)* zustar**2*zustke )1141 IF ( iom_use("us_x") ) CALL iom_put( "us_x", tmask(:,:,1)*sustke*scos_wind ) ! x surface Stokes drift 1142 IF ( iom_use("us_y") ) CALL iom_put( "us_y", tmask(:,:,1)*sustke*ssin_wind ) ! y surface Stokes drift 1143 IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rho0*tmask(:,:,1)*sustar**2*sustke ) 1131 1144 ! Stokes drift read in from sbcwave (=2). 1132 1145 CASE(2:3) … … 1138 1151 IF ( iom_use("hsw_NP") ) CALL iom_put( "hsw_NP", (0.22/grav)*wndm**2*tmask(:,:,1) ) ! significant wave height from NP spectrum 1139 1152 IF ( iom_use("wndm") ) CALL iom_put( "wndm", wndm*tmask(:,:,1) ) ! U_10 1140 IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rho0*tmask(:,:,1)* zustar**2* &1153 IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rho0*tmask(:,:,1)*sustar**2* & 1141 1154 & SQRT(ut0sd**2 + vt0sd**2 ) ) 1142 1155 END SELECT … … 1145 1158 IF ( iom_use("ghamu") ) CALL iom_put( "ghamu", umask*ghamu ) ! <uw_NL> 1146 1159 IF ( iom_use("ghamv") ) CALL iom_put( "ghamv", vmask*ghamv ) ! <vw_NL> 1147 IF ( iom_use("zwth0") ) CALL iom_put( "zwth0", tmask(:,:,1)* zwth0 ) ! <Tw_0>1148 IF ( iom_use("zws0") ) CALL iom_put( "zws0", tmask(:,:,1)* zws0 ) ! <Sw_0>1149 IF ( iom_use("zwb0") ) CALL iom_put( "zwb0", tmask(:,:,1)* zwb0 ) ! <Sw_0>1150 IF ( iom_use("zwbav") ) CALL iom_put( "zwbav", tmask(:,:,1)* zwth0 ) ! upward BL-avged turb buoyancy flux1160 IF ( iom_use("zwth0") ) CALL iom_put( "zwth0", tmask(:,:,1)*swth0 ) ! <Tw_0> 1161 IF ( iom_use("zws0") ) CALL iom_put( "zws0", tmask(:,:,1)*sws0 ) ! <Sw_0> 1162 IF ( iom_use("zwb0") ) CALL iom_put( "zwb0", tmask(:,:,1)*swb0 ) ! <Sw_0> 1163 IF ( iom_use("zwbav") ) CALL iom_put( "zwbav", tmask(:,:,1)*swth0 ) ! upward BL-avged turb buoyancy flux 1151 1164 IF ( iom_use("hbl") ) CALL iom_put( "hbl", tmask(:,:,1)*hbl ) ! boundary-layer depth 1152 1165 IF ( iom_use("ibld") ) CALL iom_put( "ibld", tmask(:,:,1)*ibld ) ! boundary-layer max k … … 1162 1175 IF ( iom_use("zdb_ml") ) CALL iom_put( "zdb_ml", tmask(:,:,1)*zdb_ml ) ! db at ml base 1163 1176 IF ( iom_use("dstokes") ) CALL iom_put( "dstokes", tmask(:,:,1)*dstokes ) ! Stokes drift penetration depth 1164 IF ( iom_use("zustke") ) CALL iom_put( "zustke", tmask(:,:,1)* zustke ) ! Stokes drift magnitude at T-points1165 IF ( iom_use("zwstrc") ) CALL iom_put( "zwstrc", tmask(:,:,1)* zwstrc ) ! convective velocity scale1166 IF ( iom_use("zwstrl") ) CALL iom_put( "zwstrl", tmask(:,:,1)* zwstrl ) ! Langmuir velocity scale1167 IF ( iom_use("zustar") ) CALL iom_put( "zustar", tmask(:,:,1)* zustar ) ! friction velocity scale1168 IF ( iom_use("zvstr") ) CALL iom_put( "zvstr", tmask(:,:,1)* zvstr ) ! mixed velocity scale1169 IF ( iom_use("zla") ) CALL iom_put( "zla", tmask(:,:,1)* zla ) ! langmuir #1170 IF ( iom_use("wind_power") ) CALL iom_put( "wind_power", 1000.*rho0*tmask(:,:,1)* zustar**3 ) ! BL depth internal to zdf_osm routine1171 IF ( iom_use("wind_wave_power") ) CALL iom_put( "wind_wave_power", 1000.*rho0*tmask(:,:,1)* zustar**2*zustke )1177 IF ( iom_use("zustke") ) CALL iom_put( "zustke", tmask(:,:,1)*sustke ) ! Stokes drift magnitude at T-points 1178 IF ( iom_use("zwstrc") ) CALL iom_put( "zwstrc", tmask(:,:,1)*swstrc ) ! convective velocity scale 1179 IF ( iom_use("zwstrl") ) CALL iom_put( "zwstrl", tmask(:,:,1)*swstrl ) ! Langmuir velocity scale 1180 IF ( iom_use("zustar") ) CALL iom_put( "zustar", tmask(:,:,1)*sustar ) ! friction velocity scale 1181 IF ( iom_use("zvstr") ) CALL iom_put( "zvstr", tmask(:,:,1)*svstr ) ! mixed velocity scale 1182 IF ( iom_use("zla") ) CALL iom_put( "zla", tmask(:,:,1)*sla ) ! langmuir # 1183 IF ( iom_use("wind_power") ) CALL iom_put( "wind_power", 1000.*rho0*tmask(:,:,1)*sustar**3 ) ! BL depth internal to zdf_osm routine 1184 IF ( iom_use("wind_wave_power") ) CALL iom_put( "wind_wave_power", 1000.*rho0*tmask(:,:,1)*sustar**2*sustke ) 1172 1185 IF ( iom_use("zhbl") ) CALL iom_put( "zhbl", tmask(:,:,1)*zhbl ) ! BL depth internal to zdf_osm routine 1173 1186 IF ( iom_use("zhml") ) CALL iom_put( "zhml", tmask(:,:,1)*zhml ) ! ML depth internal to zdf_osm routine … … 1177 1190 IF ( iom_use("zshear") ) CALL iom_put( "zshear", tmask(:,:,1)*zshear ) ! shear production of TKE internal to zdf_osm routine 1178 1191 IF ( iom_use("zdh") ) CALL iom_put( "zdh", tmask(:,:,1)*zdh ) ! pyc thicknessh internal to zdf_osm routine 1179 IF ( iom_use("zhol") ) CALL iom_put( "zhol", tmask(:,:,1)* zhol ) ! ML depth internal to zdf_osm routine1192 IF ( iom_use("zhol") ) CALL iom_put( "zhol", tmask(:,:,1)*shol ) ! ML depth internal to zdf_osm routine 1180 1193 IF ( iom_use("zwb_ent") ) CALL iom_put( "zwb_ent", tmask(:,:,1)*zwb_ent ) ! upward turb buoyancy entrainment flux 1181 1194 IF ( iom_use("zt_ml") ) CALL iom_put( "zt_ml", tmask(:,:,1)*zt_ml ) ! average T in ML … … 1237 1250 IF ( lconv(ji,jj) ) THEN 1238 1251 1239 zvel_sc_pyc = ( 0.15 * zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 + 4.25 * zshear(ji,jj) * zhbl(ji,jj) )**pthird1240 zvel_sc_ml = ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird1241 zstab_fac = ( zhml(ji,jj) / zvel_sc_ml * ( 1.4 - 0.4 / ( 1.0 + EXP(-3.5 * LOG10(- zhol(ji,jj) ) ) )**1.25 ) )**21252 zvel_sc_pyc = ( 0.15 * svstr(ji,jj)**3 + swstrc(ji,jj)**3 + 4.25 * zshear(ji,jj) * zhbl(ji,jj) )**pthird 1253 zvel_sc_ml = ( svstr(ji,jj)**3 + 0.5 * swstrc(ji,jj)**3 )**pthird 1254 zstab_fac = ( zhml(ji,jj) / zvel_sc_ml * ( 1.4 - 0.4 / ( 1.0 + EXP(-3.5 * LOG10(-shol(ji,jj) ) ) )**1.25 ) )**2 1242 1255 1243 1256 zdifml_sc(ji,jj) = rn_dif_ml * zhml(ji,jj) * zvel_sc_ml … … 1362 1375 #endif 1363 1376 ELSE 1364 zdifml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * MAX( EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ), 0.2_wp)1365 zvisml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * MAX( EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ), 0.2_wp)1377 zdifml_sc(ji,jj) = svstr(ji,jj) * zhbl(ji,jj) * MAX( EXP ( -( shol(ji,jj) / 0.6_wp )**2 ), 0.2_wp) 1378 zvisml_sc(ji,jj) = svstr(ji,jj) * zhbl(ji,jj) * MAX( EXP ( -( shol(ji,jj) / 0.6_wp )**2 ), 0.2_wp) 1366 1379 #ifdef key_osm_debug 1367 1380 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN … … 1486 1499 ! Determins stability and set flag lconv 1487 1500 DO_2D( 0, 0, 0, 0 ) 1488 IF ( zhol(ji,jj) < 0._wp ) THEN1501 IF ( shol(ji,jj) < 0._wp ) THEN 1489 1502 lconv(ji,jj) = .TRUE. 1490 1503 ELSE … … 1493 1506 END_2D 1494 1507 1495 zekman(A2D(0)) = EXP( -1.0_wp * zek * ABS( ff_t(A2D(0)) ) * zhbl(A2D(0)) / MAX( zustar(A2D(0)), 1.e-8 ) )1508 zekman(A2D(0)) = EXP( -1.0_wp * zek * ABS( ff_t(A2D(0)) ) * zhbl(A2D(0)) / MAX(sustar(A2D(0)), 1.e-8 ) ) 1496 1509 1497 1510 zshear(A2D(0)) = 0._wp … … 1509 1522 IF ( lconv(ji,jj) ) THEN 1510 1523 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 1511 zri_p(ji,jj) = MAX ( SQRT( zdb_bl(ji,jj) * zdh(ji,jj) / MAX( zdu_bl(ji,jj)**2 + zdv_bl(ji,jj)**2, 1.e-8) ) * ( zhbl(ji,jj) / zdh(ji,jj) ) * ( zvstr(ji,jj) / MAX( zustar(ji,jj), 1.e-6 ) )**2 &1524 zri_p(ji,jj) = MAX ( SQRT( zdb_bl(ji,jj) * zdh(ji,jj) / MAX( zdu_bl(ji,jj)**2 + zdv_bl(ji,jj)**2, 1.e-8) ) * ( zhbl(ji,jj) / zdh(ji,jj) ) * ( svstr(ji,jj) / MAX( sustar(ji,jj), 1.e-6 ) )**2 & 1512 1525 & / MAX( zekman(ji,jj), 1.e-6 ) , 5._wp ) 1513 1526 … … 1519 1532 zri_b(ji,jj) = zdb_ml(ji,jj) * zdh(ji,jj) / ( MAX( zdu_ml(ji,jj), 1e-5_wp )**2 + MAX( zdv_ml(ji,jj), 1e-5_wp)**2 ) 1520 1533 END IF 1521 zshear(ji,jj) = za_shr * zekman(ji,jj) * ( MAX( zustar(ji,jj)**2 * zdu_ml(ji,jj) / zhbl(ji,jj), 0._wp ) + zb_shr * MAX( -ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) * zdv_ml(ji,jj) / zhbl(ji,jj), 0._wp ) )1534 zshear(ji,jj) = za_shr * zekman(ji,jj) * ( MAX( sustar(ji,jj)**2 * zdu_ml(ji,jj) / zhbl(ji,jj), 0._wp ) + zb_shr * MAX( -ff_t(ji,jj) * sustke(ji,jj) * dstokes(ji,jj) * zdv_ml(ji,jj) / zhbl(ji,jj), 0._wp ) ) 1522 1535 #ifdef key_osm_debug 1523 1536 ! IF(narea==nn_narea_db)THEN … … 1578 1591 IF ( lconv(ji,jj) ) THEN 1579 1592 zwcor = ABS(ff_t(ji,jj)) * zhbl(ji,jj) + epsln 1580 zrf_conv = TANH( ( zwstrc(ji,jj) / zwcor )**0.69 )1581 zrf_shear = TANH( ( zustar(ji,jj) / zwcor )**0.69 )1582 zrf_langmuir = TANH( ( zwstrl(ji,jj) / zwcor )**0.69 )1593 zrf_conv = TANH( ( swstrc(ji,jj) / zwcor )**0.69 ) 1594 zrf_shear = TANH( ( sustar(ji,jj) / zwcor )**0.69 ) 1595 zrf_langmuir = TANH( ( swstrl(ji,jj) / zwcor )**0.69 ) 1583 1596 IF (nn_osm_SD_reduce > 0 ) THEN 1584 1597 ! Effective Stokes drift already reduced from surface value … … 1590 1603 & * ( 1.0 + 4.0 * dstokes(ji,jj) / hbl(ji,jj) ) ) 1591 1604 END IF 1592 zwb_ent(ji,jj) = -2.0_wp * zalpha_c * zrf_conv * zwbav(ji,jj) &1593 & - zalpha_s * zrf_shear * zustar(ji,jj)**3 / zhml(ji,jj) &1594 & + zr_stokes * ( zalpha_s * EXP( -1.5_wp * zla(ji,jj) ) * zrf_shear * zustar(ji,jj)**3 &1595 & - zrf_langmuir * zalpha_lc * zwstrl(ji,jj)**3 ) / zhml(ji,jj)1605 zwb_ent(ji,jj) = -2.0_wp * zalpha_c * zrf_conv * swbav(ji,jj) & 1606 & - zalpha_s * zrf_shear * sustar(ji,jj)**3 / zhml(ji,jj) & 1607 & + zr_stokes * ( zalpha_s * EXP( -1.5_wp * sla(ji,jj) ) * zrf_shear * sustar(ji,jj)**3 & 1608 & - zrf_langmuir * zalpha_lc * swstrl(ji,jj)**3 ) / zhml(ji,jj) 1596 1609 ! 1597 1610 #ifdef key_osm_debug … … 1645 1658 IF ( lconv(ji,jj) ) THEN 1646 1659 ! Unstable OSBL 1647 zwb_min(ji,jj) = zwb_ent(ji,jj) + zdh(ji,jj) / zhbl(ji,jj) * 2.0_wp * zwbav(ji,jj)1660 zwb_min(ji,jj) = zwb_ent(ji,jj) + zdh(ji,jj) / zhbl(ji,jj) * 2.0_wp * swbav(ji,jj) 1648 1661 END IF ! lconv 1649 1662 #ifdef key_osm_debug 1650 1663 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 1651 1664 WRITE(narea+100,'(3(a,g11.3))')'end of zdf_osm_osbl_state:zwb_ent=',zwb_ent(ji,jj), & 1652 & ' zwb_min=',zwb_min(ji,jj), ' zwb0tot=', zwb0tot(ji,jj), ' zwbav= ', zwbav(ji,jj)1665 & ' zwb_min=',zwb_min(ji,jj), ' zwb0tot=', zwb0tot(ji,jj), ' swbav= ', swbav(ji,jj) 1653 1666 FLUSH(narea+100) 1654 1667 END IF … … 1832 1845 IF ( lconv(ji,jj) ) THEN ! convective conditions 1833 1846 IF ( lpyc(ji,jj) ) THEN 1834 zzeta_m = 0.1_wp + 0.3_wp / ( 1.0_wp + EXP( -3.5_wp * LOG10( -1.0_wp * zhol(ji,jj) ) ) )1847 zzeta_m = 0.1_wp + 0.3_wp / ( 1.0_wp + EXP( -3.5_wp * LOG10( -1.0_wp * shol(ji,jj) ) ) ) 1835 1848 palpha(ji,jj) = 2.0_wp * ( 1.0_wp - ( 0.80_wp * zzeta_m + 0.5_wp * SQRT( 3.14159_wp / ppgamma_b ) ) * & 1836 1849 & zdbdz_bl_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj) ) / & … … 1874 1887 IF ( zdhdt(ji,jj) > 0.0_wp ) THEN ! Depth increasing, or steady. 1875 1888 IF ( zdb_bl(ji,jj) > 0.0_wp ) THEN 1876 IF ( zhol(ji,jj) >= 0.5_wp ) THEN ! Very stable - 'thick' pycnocline1889 IF ( shol(ji,jj) >= 0.5_wp ) THEN ! Very stable - 'thick' pycnocline 1877 1890 ztmp = 1.0_wp / MAX( zhbl(ji,jj), epsln ) 1878 1891 zbgrad = zdb_bl(ji,jj) * ztmp … … 1888 1901 pdbdz(ji,jj,jk) = zbgrad * EXP( -1.75_wp * ( znd + 0.75_wp )**2 ) 1889 1902 END DO 1890 ENDIF ! IF ( zhol >=0.5)1903 ENDIF ! IF (shol >=0.5) 1891 1904 #ifdef key_osm_debug 1892 1905 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN … … 1947 1960 zwb_fk_b(ji,jj) = 0.5 * zwb_fk(ji,jj) * hmle(ji,jj) / hbl(ji,jj) 1948 1961 ENDIF 1949 zvel_max = ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**p2third / hbl(ji,jj)1962 zvel_max = ( svstr(ji,jj)**3 + 0.5 * swstrc(ji,jj)**3 )**p2third / hbl(ji,jj) 1950 1963 IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0.0 ) THEN 1951 1964 ! OSBL is deepening, entrainment > restratification … … 1953 1966 zgamma_b_nd = MAX( zdbdz_bl_ext(ji,jj), 0.0_wp ) * zdh(ji,jj) / ( zvel_max + MAX( zdb_bl(ji,jj), 1e-15 ) ) 1954 1967 zpsi = ( 1.0_wp - 0.5_wp * zdh(ji,jj) / zhbl(ji,jj) ) * & 1955 & ( zwb0(ji,jj) - MIN( ( zwb_min(ji,jj) + 2.0_wp * zwb_fk_b(ji,jj) ), 0.0_wp ) ) * zdh(ji,jj) / zhbl(ji,jj)1968 & ( swb0(ji,jj) - MIN( ( zwb_min(ji,jj) + 2.0_wp * zwb_fk_b(ji,jj) ), 0.0_wp ) ) * zdh(ji,jj) / zhbl(ji,jj) 1956 1969 zpsi = zpsi + 1.75_wp * ( 1.0_wp - 0.5_wp * zdh(ji,jj) / zhbl(ji,jj) ) * & 1957 1970 & ( zdh(ji,jj) / zhbl(ji,jj) + zgamma_b_nd ) * MIN( ( zwb_min(ji,jj) + 2.0_wp * zwb_fk_b(ji,jj) ), 0.0_wp ) … … 1968 1981 #endif 1969 1982 IF ( j_ddh(ji,jj) == 1 ) THEN 1970 IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN1971 zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zvstr(ji,jj)**2 , 1.e-12 )) ), 0.2d0 )1983 IF ( ( swstrc(ji,jj) / svstr(ji,jj) )**3 <= 0.5 ) THEN 1984 zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * svstr(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) 1972 1985 ELSE 1973 zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zwstrc(ji,jj)**2 , 1.e-12 )) ), 0.2d0 )1986 zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * swstrc(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) 1974 1987 ENDIF 1975 1988 ! Relaxation to dh_ref = zari * hbl … … 1987 2000 zddhdt = -1.0_wp * a_ddh * ( 1.0 - 1.6_wp * zdh(ji,jj) / zhbl(ji,jj) ) * zwb_ent(ji,jj) / & 1988 2001 & ( zvel_max + MAX( zdb_bl(ji,jj), 1e-15 ) ) 1989 zddhdt = EXP( -4.0_wp * ABS( ff_t(ji,jj) ) * zhbl(ji,jj) / MAX( zustar(ji,jj), 1e-8_wp ) ) * zddhdt2002 zddhdt = EXP( -4.0_wp * ABS( ff_t(ji,jj) ) * zhbl(ji,jj) / MAX(sustar(ji,jj), 1e-8_wp ) ) * zddhdt 1990 2003 ELSE 1991 2004 zddhdt = 0.0_wp … … 2006 2019 ! Fox-Kemper not used. 2007 2020 2008 zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_Dt / hbl(ji,jj) ) * zwb_ent(ji,jj) / &2009 & MAX(( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird, epsln)2021 zvel_max = - ( 1.0 + 1.0 * ( svstr(ji,jj)**3 + 0.5 * swstrc(ji,jj)**3 )**pthird * rn_Dt / hbl(ji,jj) ) * zwb_ent(ji,jj) / & 2022 & MAX((svstr(ji,jj)**3 + 0.5 * swstrc(ji,jj)**3 )**pthird, epsln) 2010 2023 zdhdt(ji,jj) = -zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) 2011 2024 ! added ajgn 23 July as temporay fix … … 2014 2027 2015 2028 ELSE ! lconv - Stable 2016 zdhdt(ji,jj) = ( 0.06 + 0.52 * zhol(ji,jj) / 2.0 ) * zvstr(ji,jj)**3 / hbl(ji,jj) + zwbav(ji,jj)2029 zdhdt(ji,jj) = ( 0.06 + 0.52 * shol(ji,jj) / 2.0 ) * svstr(ji,jj)**3 / hbl(ji,jj) + swbav(ji,jj) 2017 2030 IF ( zdhdt(ji,jj) < 0._wp ) THEN 2018 2031 ! For long timsteps factor in brackets slows the rapid collapse of the OSBL 2019 zpert = 2.0 * ( 1.0 + 0.0 * 2.0 * zvstr(ji,jj) * rn_Dt / hbl(ji,jj) ) * zvstr(ji,jj)**2 / hbl(ji,jj)2032 zpert = 2.0 * ( 1.0 + 0.0 * 2.0 * svstr(ji,jj) * rn_Dt / hbl(ji,jj) ) * svstr(ji,jj)**2 / hbl(ji,jj) 2020 2033 ELSE 2021 zpert = MAX( zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) )2034 zpert = MAX( svstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) ) 2022 2035 ENDIF 2023 2036 zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / MAX(zpert, epsln) … … 2036 2049 zwb_fk_b(ji,jj) = 0.5 * zwb_fk(ji,jj) * hmle(ji,jj) / hbl(ji,jj) 2037 2050 ENDIF 2038 zvel_max = ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**p2third / hbl(ji,jj)2051 zvel_max = ( swstrl(ji,jj)**3 + 0.5 * swstrc(ji,jj)**3 )**p2third / hbl(ji,jj) 2039 2052 IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0.0 ) THEN 2040 2053 ! OSBL is deepening, entrainment > restratification … … 2055 2068 2056 2069 zvel_max = -zwb_ent(ji,jj) / & 2057 & MAX(( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird, epsln)2070 & MAX((svstr(ji,jj)**3 + 0.5 * swstrc(ji,jj)**3 )**pthird, epsln) 2058 2071 zdhdt(ji,jj) = -zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) 2059 2072 ! added ajgn 23 July as temporay fix … … 2062 2075 2063 2076 ELSE ! Stable 2064 zdhdt(ji,jj) = ( 0.06 + 0.52 * zhol(ji,jj) / 2.0 ) * zvstr(ji,jj)**3 / hbl(ji,jj) + zwbav(ji,jj)2077 zdhdt(ji,jj) = ( 0.06 + 0.52 * shol(ji,jj) / 2.0 ) * svstr(ji,jj)**3 / hbl(ji,jj) + swbav(ji,jj) 2065 2078 IF ( zdhdt(ji,jj) < 0._wp ) THEN 2066 2079 ! For long timsteps factor in brackets slows the rapid collapse of the OSBL 2067 zpert = 2.0 * zvstr(ji,jj)**2 / hbl(ji,jj)2080 zpert = 2.0 * svstr(ji,jj)**2 / hbl(ji,jj) 2068 2081 ELSE 2069 zpert = MAX( zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) )2082 zpert = MAX( svstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) ) 2070 2083 ENDIF 2071 2084 zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / MAX(zpert, epsln) … … 2131 2144 2132 2145 IF( ln_osm_mle ) THEN 2133 zvel_max = ( zwstrl(ji,jj)**3 + zwstrc(ji,jj)**3 )**p2third / hbl(ji,jj)2146 zvel_max = ( swstrl(ji,jj)**3 + swstrc(ji,jj)**3 )**p2third / hbl(ji,jj) 2134 2147 ELSE 2135 2148 2136 zvel_max = -( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_Dt / hbl(ji,jj) ) * zwb_ent(ji,jj) / &2137 & ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird2149 zvel_max = -( 1.0 + 1.0 * ( svstr(ji,jj)**3 + 0.5 * swstrc(ji,jj)**3 )**pthird * rn_Dt / hbl(ji,jj) ) * zwb_ent(ji,jj) / & 2150 & ( svstr(ji,jj)**3 + 0.5 * swstrc(ji,jj)**3 )**pthird 2138 2151 2139 2152 ENDIF … … 2190 2203 & - zbeta * ( zs_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ),& 2191 2204 & 0.0 ) + & 2192 & 2.0 * zvstr(ji,jj)**2 / zhbl_s2205 & 2.0 * svstr(ji,jj)**2 / zhbl_s 2193 2206 2194 2207 ! Alan is thuis right? I have simply changed hbli to hbl 2195 zhol(ji,jj) = -zhbl_s / ( ( zvstr(ji,jj)**3 + epsln )/ zwbav(ji,jj) )2196 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) ) ) * &2197 & zustar(ji,jj)**3 / zhbl_s ) * ( 0.725 + 0.225 * EXP( -7.5 * zhol(ji,jj) ) )2198 zdhdt(ji,jj) = zdhdt(ji,jj) + zwbav(ji,jj)2208 shol(ji,jj) = -zhbl_s / ( ( svstr(ji,jj)**3 + epsln )/ swbav(ji,jj) ) 2209 zdhdt(ji,jj) = -( swbav(ji,jj) - 0.04 / 2.0 * swstrl(ji,jj)**3 / zhbl_s - 0.15 / 2.0 * ( 1.0 - EXP( -1.5 * sla(ji,jj) ) ) * & 2210 & sustar(ji,jj)**3 / zhbl_s ) * ( 0.725 + 0.225 * EXP( -7.5 * shol(ji,jj) ) ) 2211 zdhdt(ji,jj) = zdhdt(ji,jj) + swbav(ji,jj) 2199 2212 zhbl_s = zhbl_s + MIN( zdhdt(ji,jj) / zdb * rn_Dt / FLOAT( ibld(ji,jj) - imld(ji,jj) ), e3w(ji,jj,jm,Kmm) ) 2200 2213 … … 2208 2221 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 2209 2222 WRITE(narea+100,'(2(a,i7))')' jk=',jk,' jm=', jm 2210 WRITE(narea+100,'(4(a,g11.3),a,l7)')'zdb=',zdb,' zhol',zhol(ji,jj),' zdhdt',zdhdt(ji,jj),' zhbl_s=', zhbl_s,' lpyc=',lpyc(ji,jj)2223 WRITE(narea+100,'(4(a,g11.3),a,l7)')'zdb=',zdb,' shol',shol(ji,jj),' zdhdt',zdhdt(ji,jj),' zhbl_s=', zhbl_s,' lpyc=',lpyc(ji,jj) 2211 2224 FLUSH(narea+100) 2212 2225 END IF … … 2260 2273 IF ( zdb_bl(ji,jj) > 1e-15_wp ) THEN 2261 2274 IF ( j_ddh(ji,jj) == 0 ) THEN 2262 zvel_max = ( zvstr(ji,jj)**3 + 0.5_wp * zwstrc(ji,jj)**3 )**p2third / hbl(ji,jj)2275 zvel_max = ( svstr(ji,jj)**3 + 0.5_wp * swstrc(ji,jj)**3 )**p2third / hbl(ji,jj) 2263 2276 ! ddhdt for pycnocline determined in osm_calculate_dhdt 2264 2277 zddhdt = -a_ddh * ( 1.0 - 1.6 * zdh(ji,jj) / zhbl(ji,jj) ) * zwb_ent(ji,jj) / ( zvel_max + MAX( zdb_bl(ji,jj), 1e-15 ) ) 2265 zddhdt = EXP( -4.0_wp * ABS( ff_t(ji,jj) ) * zhbl(ji,jj) / MAX( zustar(ji,jj), 1e-8 ) ) * zddhdt2278 zddhdt = EXP( -4.0_wp * ABS( ff_t(ji,jj) ) * zhbl(ji,jj) / MAX( sustar(ji,jj), 1e-8 ) ) * zddhdt 2266 2279 ! maximum limit for how thick the shear layer can grow relative to the thickness of the boundary kayer 2267 2280 dh(ji,jj) = MIN( dh(ji,jj) + zddhdt * rn_Dt, 0.625_wp * hbl(ji,jj) ) 2268 2281 ELSE 2269 2282 ! Need to recalculate because hbl has been updated. 2270 IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN2271 zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zvstr(ji,jj)**2 , 1.e-12 )) ), 0.2d0 )2283 IF ( ( swstrc(ji,jj) / svstr(ji,jj) )**3 <= 0.5 ) THEN 2284 zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * svstr(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) 2272 2285 ELSE 2273 zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zwstrc(ji,jj)**2 , 1.e-12 )) ), 0.2d0 )2286 zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * swstrc(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) 2274 2287 ENDIF 2275 2288 ztau = MAX( zdb_bl(ji,jj) * ( zari * hbl(ji,jj) ) / ( a_ddh_2 * MAX(-zwb_ent(ji,jj), 1.e-12) ), 2.0 * rn_Dt ) … … 2278 2291 ENDIF 2279 2292 ELSE 2280 ztau = MAX( MAX( hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5_wp * zwstrc(ji,jj)**3 )**pthird, epsln), 2.0_wp * rn_Dt )2293 ztau = MAX( MAX( hbl(ji,jj) / ( svstr(ji,jj)**3 + 0.5_wp * swstrc(ji,jj)**3 )**pthird, epsln), 2.0_wp * rn_Dt ) 2281 2294 dh(ji,jj) = dh(ji,jj) * EXP( -1.0_wp * rn_Dt / ztau ) + 0.2_wp * zhbl(ji,jj) * ( 1.0_wp - EXP( -1.0_wp * rn_Dt / ztau ) ) 2282 2295 IF ( dh(ji,jj) > hbl(ji,jj) ) dh(ji,jj) = 0.2_wp * hbl(ji,jj) … … 2285 2298 ! Initially shear only for entraining OSBL. Stable code will be needed if extended to stable OSBL 2286 2299 2287 ztau = hbl(ji,jj) / MAX( zvstr(ji,jj), epsln)2300 ztau = hbl(ji,jj) / MAX(svstr(ji,jj), epsln) 2288 2301 IF ( zdhdt(ji,jj) >= 0.0 ) THEN ! probably shouldn't include wm here 2289 2302 ! boundary layer deepening 2290 2303 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 2291 2304 ! pycnocline thickness set by stratification - use same relationship as for neutral conditions. 2292 zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) &2305 zari = MIN( 4.5 * ( svstr(ji,jj)**2 ) & 2293 2306 & / MAX(zdb_bl(ji,jj) * zhbl(ji,jj), epsln ) + 0.01 , 0.2 ) 2294 2307 zdh_ref = MIN( zari, 0.2 ) * hbl(ji,jj) … … 2312 2325 ! OSBL is deepening. Note wb_fk_b is zero if ln_osm_mle=F 2313 2326 IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_bl_ext(ji,jj) > 0._wp)THEN 2314 IF ( ( zwstrc(ji,jj) / MAX(zvstr(ji,jj), epsln) )**3 <= 0.5 ) THEN ! near neutral stability2315 zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zvstr(ji,jj)**2 , 1.e-12 )) ), 0.2d0 )2327 IF ( ( swstrc(ji,jj) / MAX(svstr(ji,jj), epsln) )**3 <= 0.5 ) THEN ! near neutral stability 2328 zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * svstr(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) 2316 2329 ELSE ! unstable 2317 zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zwstrc(ji,jj)**2 , 1.e-12 )) ), 0.2d0 )2330 zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * swstrc(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) 2318 2331 ENDIF 2319 ztau = 0.2 * hbl(ji,jj) / MAX(epsln, ( zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird)2332 ztau = 0.2 * hbl(ji,jj) / MAX(epsln, (svstr(ji,jj)**3 + 0.5 *swstrc(ji,jj)**3)**pthird) 2320 2333 zdh_ref = zari * hbl(ji,jj) 2321 2334 ELSE 2322 ztau = 0.2 * hbl(ji,jj) / MAX(epsln, ( zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird)2335 ztau = 0.2 * hbl(ji,jj) / MAX(epsln, (svstr(ji,jj)**3 + 0.5 *swstrc(ji,jj)**3)**pthird) 2323 2336 zdh_ref = 0.2 * hbl(ji,jj) 2324 2337 ENDIF 2325 2338 ELSE 2326 ztau = 0.2 * hbl(ji,jj) / MAX(epsln, ( zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird)2339 ztau = 0.2 * hbl(ji,jj) / MAX(epsln, (svstr(ji,jj)**3 + 0.5 *swstrc(ji,jj)**3)**pthird) 2327 2340 zdh_ref = 0.2 * hbl(ji,jj) 2328 2341 ENDIF 2329 2342 ELSE ! ln_osm_mle 2330 2343 IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_bl_ext(ji,jj) > 0._wp)THEN 2331 IF ( ( zwstrc(ji,jj) / MAX(zvstr(ji,jj), epsln) )**3 <= 0.5 ) THEN ! near neutral stability2332 zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zvstr(ji,jj)**2 , 1.e-12 )) ), 0.2d0 )2344 IF ( ( swstrc(ji,jj) / MAX(svstr(ji,jj), epsln) )**3 <= 0.5 ) THEN ! near neutral stability 2345 zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * svstr(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) 2333 2346 ELSE ! unstable 2334 zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zwstrc(ji,jj)**2 , 1.e-12 )) ), 0.2d0 )2347 zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * swstrc(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) 2335 2348 ENDIF 2336 ztau = hbl(ji,jj) / MAX(epsln, ( zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird)2349 ztau = hbl(ji,jj) / MAX(epsln, (svstr(ji,jj)**3 + 0.5 *swstrc(ji,jj)**3)**pthird) 2337 2350 zdh_ref = zari * hbl(ji,jj) 2338 2351 ELSE 2339 ztau = hbl(ji,jj) / MAX(epsln, ( zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird)2352 ztau = hbl(ji,jj) / MAX(epsln, (svstr(ji,jj)**3 + 0.5 *swstrc(ji,jj)**3)**pthird) 2340 2353 zdh_ref = 0.2 * hbl(ji,jj) 2341 2354 ENDIF … … 2348 2361 ! Alan: this hml is never defined or used 2349 2362 ELSE ! IF (lconv) 2350 ztau = hbl(ji,jj) / MAX( zvstr(ji,jj), epsln)2363 ztau = hbl(ji,jj) / MAX(svstr(ji,jj), epsln) 2351 2364 IF ( zdhdt(ji,jj) >= 0.0 ) THEN ! probably shouldn't include wm here 2352 2365 ! boundary layer deepening 2353 2366 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 2354 2367 ! pycnocline thickness set by stratification - use same relationship as for neutral conditions. 2355 zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) &2368 zari = MIN( 4.5 * ( svstr(ji,jj)**2 ) & 2356 2369 & / MAX(zdb_bl(ji,jj) * zhbl(ji,jj), epsln ) + 0.01 , 0.2 ) 2357 2370 zdh_ref = MIN( zari, 0.2 ) * hbl(ji,jj) … … 2656 2669 END SUBROUTINE zdf_osm_vertical_average 2657 2670 2658 SUBROUTINE zdf_osm_velocity_rotation_2d( pu, pv, pcos_w, psin_w,fwd )2671 SUBROUTINE zdf_osm_velocity_rotation_2d( pu, pv, fwd ) 2659 2672 !!--------------------------------------------------------------------- 2660 2673 !! *** ROUTINE zdf_velocity_rotation_2d *** … … 2664 2677 !! 2665 2678 !! ** Method : The velocity components are rotated into (fwd=.TRUE.) or 2666 !! from (fwd=.FALSE.) the frame specified by pcos_w and psin_w 2679 !! from (fwd=.FALSE.) the frame specified by scos_wind and 2680 !! ssin_wind 2667 2681 !! 2668 2682 !!---------------------------------------------------------------------- 2669 2683 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) :: pu, pv ! Components of current 2670 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: pcos_w, psin_w ! Cos and sin of rotation angle2671 2684 LOGICAL, OPTIONAL, INTENT(in ) :: fwd ! Forward (default) or reverse rotation 2672 2685 ! … … 2680 2693 DO_2D( 0, 0, 0, 0 ) 2681 2694 ztmp = pu(ji,jj) 2682 pu(ji,jj) = pu(ji,jj) * pcos_w(ji,jj) + zfwd * pv(ji,jj) * psin_w(ji,jj)2683 pv(ji,jj) = pv(ji,jj) * pcos_w(ji,jj) - zfwd * ztmp * psin_w(ji,jj)2695 pu(ji,jj) = pu(ji,jj) * scos_wind(ji,jj) + zfwd * pv(ji,jj) * ssin_wind(ji,jj) 2696 pv(ji,jj) = pv(ji,jj) * scos_wind(ji,jj) - zfwd * ztmp * ssin_wind(ji,jj) 2684 2697 END_2D 2685 2698 ! … … 2688 2701 END SUBROUTINE zdf_osm_velocity_rotation_2d 2689 2702 2690 SUBROUTINE zdf_osm_velocity_rotation_3d( pu, pv, pcos_w, psin_w,fwd, ktop, knlev )2703 SUBROUTINE zdf_osm_velocity_rotation_3d( pu, pv, fwd, ktop, knlev ) 2691 2704 !!--------------------------------------------------------------------- 2692 2705 !! *** ROUTINE zdf_velocity_rotation_3d *** … … 2696 2709 !! 2697 2710 !! ** Method : The velocity components are rotated into (fwd=.TRUE.) or 2698 !! from (fwd=.FALSE.) the frame specified by pcos_wand2699 !! psin_w; optionally, the rotation can be restricted at each2700 !! water column to span from the a minimum index ktop to the2701 !! depth index specified in array knlev2711 !! from (fwd=.FALSE.) the frame specified by scos_wind and 2712 !! ssin_wind; optionally, the rotation can be restricted at 2713 !! each water column to span from the a minimum index ktop to 2714 !! the depth index specified in array knlev 2702 2715 !! 2703 2716 !!---------------------------------------------------------------------- 2704 2717 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pu, pv ! Components of current 2705 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: pcos_w, psin_w ! Cos and sin of rotation angle2706 2718 LOGICAL, OPTIONAL, INTENT(in ) :: fwd ! Forward (default) or reverse rotation 2707 2719 INTEGER, OPTIONAL, INTENT(in ) :: ktop ! Minimum depth index … … 2731 2743 IF ( llkbot .OR. knlev(ji,jj) >= jk ) THEN 2732 2744 ztmp = pu(ji,jj,jk) 2733 pu(ji,jj,jk) = pu(ji,jj,jk) * pcos_w(ji,jj) + zfwd * pv(ji,jj,jk) * psin_w(ji,jj)2734 pv(ji,jj,jk) = pv(ji,jj,jk) * pcos_w(ji,jj) - zfwd * ztmp * psin_w(ji,jj)2745 pu(ji,jj,jk) = pu(ji,jj,jk) * scos_wind(ji,jj) + zfwd * pv(ji,jj,jk) * ssin_wind(ji,jj) 2746 pv(ji,jj,jk) = pv(ji,jj,jk) * scos_wind(ji,jj) - zfwd * ztmp * ssin_wind(ji,jj) 2735 2747 END IF 2736 2748 END_3D … … 2740 2752 END SUBROUTINE zdf_osm_velocity_rotation_3d 2741 2753 2742 SUBROUTINE zdf_osm_fgr_terms( Kmm, kbld, kmld, kp_ext, ldconv, ldpyc, k_ddh, phbl, phml, pdh, pdhdt, phol, pshear, & 2743 & pustar, pwstrl, pvstr, pwstrc, puw0, pwth0, pws0, pwb0, pwthav, pwsav, pwbav, pustke, pla, & 2754 SUBROUTINE zdf_osm_fgr_terms( Kmm, kbld, kmld, kp_ext, ldconv, ldpyc, k_ddh, phbl, phml, pdh, pdhdt, pshear, & 2744 2755 & pdt_bl, pds_bl, pdb_bl, pdu_bl, pdv_bl, pdt_ml, pds_ml, pdb_ml, pdu_ml, pdv_ml, & 2745 2756 & pdtdz_bl_ext, pdsdz_bl_ext, pdbdz_bl_ext, pdbdz_pyc, palpha_pyc, pdiffut, pviscos ) … … 2763 2774 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdh ! Pycnocline depth 2764 2775 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdhdt ! BL depth tendency 2765 REAL(wp), DIMENSION(:,:), INTENT(in ) :: phol ! Stability parameter for boundary layer2766 2776 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pshear ! Shear production 2767 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pustar ! Friction velocity2768 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pwstrl ! Langmuir velocity scale2769 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pvstr ! Velocity scale (approaches zustar for large Langmuir number)2770 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pwstrc ! Convective velocity scale2771 REAL(wp), DIMENSION(:,:), INTENT(in ) :: puw0 ! Surface u-momentum flux2772 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pwth0 ! Surface heat flux2773 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pws0 ! Surface freshwater flux2774 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pwb0 ! Surface buoyancy flux2775 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pwthav ! BL average heat flux2776 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pwsav ! BL average freshwater flux2777 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pwbav ! BL average buoyancy flux2778 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pustke ! Surface Stokes drift2779 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pla ! Langmuir number2780 2777 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdt_bl ! Temperature diff. between BL average and basal value 2781 2778 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pds_bl ! Salinity diff. between BL average and basal value … … 2843 2840 ! ------------------------------------------------------ 2844 2841 WHERE ( ldconv(A2D(0)) ) 2845 zsc_wth_1(:,:) = pwstrl(A2D(0))**3 * pwth0(A2D(0)) / ( pvstr(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 + epsln )2846 zsc_ws_1(:,:) = pwstrl(A2D(0))**3 * pws0(A2D(0)) / ( pvstr(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 + epsln )2842 zsc_wth_1(:,:) = swstrl(A2D(0))**3 * swth0(A2D(0)) / ( svstr(A2D(0))**3 + 0.5_wp * swstrc(A2D(0))**3 + epsln ) 2843 zsc_ws_1(:,:) = swstrl(A2D(0))**3 * sws0(A2D(0)) / ( svstr(A2D(0))**3 + 0.5_wp * swstrc(A2D(0))**3 + epsln ) 2847 2844 ELSEWHERE 2848 zsc_wth_1(:,:) = 2.0_wp * pwthav(A2D(0))2849 zsc_ws_1(:,:) = 2.0_wp * pwsav(A2D(0))2845 zsc_wth_1(:,:) = 2.0_wp * swthav(A2D(0)) 2846 zsc_ws_1(:,:) = 2.0_wp * swsav(A2D(0)) 2850 2847 ENDWHERE 2851 2848 DO_3D( 0, 0, 0, 0, 2, MAX( jkm_mld, jkm_bld ) ) … … 2875 2872 ! 2876 2873 ! Stokes term in flux-gradient relationship (note in zsc_uw_n don't use 2877 ! zvstr since term needs to go to zero as zwstrl goes to zero)2874 ! svstr since term needs to go to zero as swstrl goes to zero) 2878 2875 ! --------------------------------------------------------------------- 2879 2876 WHERE ( ldconv(A2D(0)) ) 2880 zsc_uw_1(:,:) = ( pwstrl(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 )**pthird * pustke(A2D(0)) / &2881 & MAX( ( 1.0_wp - 1.0_wp * 6.5_wp * pla(A2D(0))**( 8.0_wp / 3.0_wp ) ), 0.2_wp )2882 zsc_uw_2(:,:) = ( pwstrl(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 )**pthird * pustke(A2D(0)) / &2883 & MIN( pla(A2D(0))**( 8.0_wp / 3.0_wp ) + epsln, 0.12_wp )2884 zsc_vw_1(:,:) = ff_t(A2D(0)) * phml(A2D(0)) * pustke(A2D(0))**3 * MIN( pla(A2D(0))**( 8.0_wp / 3.0_wp ), 0.12_wp ) / &2885 & ( ( pvstr(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 )**( 2.0_wp / 3.0_wp ) + epsln )2877 zsc_uw_1(:,:) = ( swstrl(A2D(0))**3 + 0.5_wp * swstrc(A2D(0))**3 )**pthird * sustke(A2D(0)) / & 2878 & MAX( ( 1.0_wp - 1.0_wp * 6.5_wp * sla(A2D(0))**( 8.0_wp / 3.0_wp ) ), 0.2_wp ) 2879 zsc_uw_2(:,:) = ( swstrl(A2D(0))**3 + 0.5_wp * swstrc(A2D(0))**3 )**pthird * sustke(A2D(0)) / & 2880 & MIN( sla(A2D(0))**( 8.0_wp / 3.0_wp ) + epsln, 0.12_wp ) 2881 zsc_vw_1(:,:) = ff_t(A2D(0)) * phml(A2D(0)) * sustke(A2D(0))**3 * MIN( sla(A2D(0))**( 8.0_wp / 3.0_wp ), 0.12_wp ) / & 2882 & ( ( svstr(A2D(0))**3 + 0.5_wp * swstrc(A2D(0))**3 )**( 2.0_wp / 3.0_wp ) + epsln ) 2886 2883 ELSEWHERE 2887 zsc_uw_1(:,:) = pustar(A2D(0))**22888 zsc_vw_1(:,:) = ff_t(A2D(0)) * phbl(A2D(0)) * pustke(A2D(0))**3 * MIN( pla(A2D(0))**( 8.0_wp / 3.0_wp ), 0.12_wp ) / &2889 & ( pvstr(A2D(0))**2 + epsln )2884 zsc_uw_1(:,:) = sustar(A2D(0))**2 2885 zsc_vw_1(:,:) = ff_t(A2D(0)) * phbl(A2D(0)) * sustke(A2D(0))**3 * MIN( sla(A2D(0))**( 8.0_wp / 3.0_wp ), 0.12_wp ) / & 2886 & ( svstr(A2D(0))**2 + epsln ) 2890 2887 ENDWHERE 2891 2888 DO_3D( 0, 0, 0, 0, 2, MAX( jkm_mld, jkm_bld ) ) … … 2931 2928 ! ---------------------------------------------------------------------- 2932 2929 WHERE ( ldconv(A2D(0)) ) 2933 zsc_wth_1(:,:) = pwbav(A2D(0)) * pwth0(A2D(0)) * ( 1.0_wp + EXP( 0.2_wp * phol(A2D(0)) ) ) * phml(A2D(0)) / &2934 & ( pvstr(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 + epsln )2935 zsc_ws_1(:,:) = pwbav(A2D(0)) * pws0(A2D(0)) * ( 1.0_wp + EXP( 0.2_wp * phol(A2D(0)) ) ) * phml(A2D(0)) / &2936 & ( pvstr(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 + epsln )2930 zsc_wth_1(:,:) = swbav(A2D(0)) * swth0(A2D(0)) * ( 1.0_wp + EXP( 0.2_wp * shol(A2D(0)) ) ) * phml(A2D(0)) / & 2931 & ( svstr(A2D(0))**3 + 0.5_wp * swstrc(A2D(0))**3 + epsln ) 2932 zsc_ws_1(:,:) = swbav(A2D(0)) * sws0(A2D(0)) * ( 1.0_wp + EXP( 0.2_wp * shol(A2D(0)) ) ) * phml(A2D(0)) / & 2933 & ( svstr(A2D(0))**3 + 0.5_wp * swstrc(A2D(0))**3 + epsln ) 2937 2934 ELSEWHERE 2938 2935 zsc_wth_1(:,:) = 0.0_wp … … 2948 2945 zl_l = 2.0_wp * ( 1.0_wp - EXP( -2.0_wp * ( zznd_ml + zznd_ml**3 / 3.0_wp ) ) ) * & 2949 2946 & ( 1.0_wp - EXP( -8.0_wp * ( 1.15_wp - zznd_ml ) ) ) * ( 1.0_wp + dstokes(ji,jj) / phml (ji,jj) ) 2950 zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0_wp + EXP( -3.0_wp * LOG10( -1.0_wp * phol(ji,jj) ) ) )**( 3.0_wp / 2.0_wp )2947 zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0_wp + EXP( -3.0_wp * LOG10( -1.0_wp * shol(ji,jj) ) ) )**( 3.0_wp / 2.0_wp ) 2951 2948 ! Non-gradient buoyancy terms 2952 2949 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3_wp * 0.4_wp * zsc_wth_1(ji,jj) * zl_eps / ( 0.15_wp + zznd_ml ) … … 2962 2959 DO_2D( 0, 0, 0, 0 ) 2963 2960 IF ( ldconv(ji,jj) .AND. ldpyc(ji,jj) ) THEN 2964 ztau_sc_u(ji,jj) = phml(ji,jj) / ( pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**pthird * &2965 & ( 1.4_wp - 0.4_wp / ( 1.0_wp + EXP( -3.5_wp * LOG10( -1.0_wp * phol(ji,jj) ) ) )**1.5_wp )2966 zwth_ent(ji,jj) = -0.003_wp * ( 0.15_wp * pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**pthird * &2961 ztau_sc_u(ji,jj) = phml(ji,jj) / ( svstr(ji,jj)**3 + swstrc(ji,jj)**3 )**pthird * & 2962 & ( 1.4_wp - 0.4_wp / ( 1.0_wp + EXP( -3.5_wp * LOG10( -1.0_wp * shol(ji,jj) ) ) )**1.5_wp ) 2963 zwth_ent(ji,jj) = -0.003_wp * ( 0.15_wp * svstr(ji,jj)**3 + swstrc(ji,jj)**3 )**pthird * & 2967 2964 & ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pdt_ml(ji,jj) 2968 zws_ent(ji,jj) = -0.003_wp * ( 0.15_wp * pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**pthird * &2965 zws_ent(ji,jj) = -0.003_wp * ( 0.15_wp * svstr(ji,jj)**3 + swstrc(ji,jj)**3 )**pthird * & 2969 2966 & ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pds_ml(ji,jj) 2970 2967 IF ( dh(ji,jj) < 0.2_wp * hbl(ji,jj) ) THEN 2971 2968 zbuoy_pyc_sc = 2.0_wp * MAX( pdb_ml(ji,jj), 0.0_wp ) / pdh(ji,jj) 2972 zdelta_pyc = ( pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**pthird / &2973 & SQRT( MAX( zbuoy_pyc_sc, ( pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**p2third / pdh(ji,jj)**2 ) )2969 zdelta_pyc = ( svstr(ji,jj)**3 + swstrc(ji,jj)**3 )**pthird / & 2970 & SQRT( MAX( zbuoy_pyc_sc, ( svstr(ji,jj)**3 + swstrc(ji,jj)**3 )**p2third / pdh(ji,jj)**2 ) ) 2974 2971 zwt_pyc_sc_1(ji,jj) = 0.325_wp * ( palpha_pyc(ji,jj) * pdt_ml(ji,jj) / pdh(ji,jj) + pdtdz_bl_ext(ji,jj) ) * & 2975 2972 & zdelta_pyc**2 / pdh(ji,jj) 2976 2973 zws_pyc_sc_1(ji,jj) = 0.325_wp * ( palpha_pyc(ji,jj) * pds_ml(ji,jj) / pdh(ji,jj) + pdsdz_bl_ext(ji,jj) ) * & 2977 2974 & zdelta_pyc**2 / pdh(ji,jj) 2978 zzeta_pyc(ji,jj) = 0.15_wp - 0.175_wp / ( 1.0_wp + EXP( -3.5_wp * LOG10( -1.0_wp * phol(ji,jj) ) ) )2975 zzeta_pyc(ji,jj) = 0.15_wp - 0.175_wp / ( 1.0_wp + EXP( -3.5_wp * LOG10( -1.0_wp * shol(ji,jj) ) ) ) 2979 2976 #ifdef key_osm_debug 2980 2977 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN … … 3016 3013 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.05_wp * zwt_pyc_sc_1(ji,jj) * & 3017 3014 & EXP( -0.25_wp * ( zznd_pyc / zzeta_pyc(ji,jj) )**2 ) * & 3018 & pdh(ji,jj) / ( pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**pthird3015 & pdh(ji,jj) / ( svstr(ji,jj)**3 + swstrc(ji,jj)**3 )**pthird 3019 3016 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.05_wp * zws_pyc_sc_1(ji,jj) * & 3020 3017 & EXP( -0.25_wp * ( zznd_pyc / zzeta_pyc(ji,jj) )**2 ) * & 3021 & pdh(ji,jj) / ( pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**pthird3018 & pdh(ji,jj) / ( svstr(ji,jj)**3 + swstrc(ji,jj)**3 )**pthird 3022 3019 END IF 3023 3020 END IF ! End of pycnocline … … 3031 3028 zsc_vw_1(:,:) = 0.0_wp 3032 3029 WHERE ( ldconv(A2D(0)) ) 3033 zsc_uw_1(:,:) = -1.0_wp * pwb0(A2D(0)) * pustar(A2D(0))**2 * phml(A2D(0)) / &3034 & ( pvstr(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 + epsln )3035 zsc_uw_2(:,:) = pwb0(A2D(0)) * pustke(A2D(0)) * phml(A2D(0)) / &3036 & ( pvstr(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 + epsln )**( 2.0_wp / 3.0_wp )3030 zsc_uw_1(:,:) = -1.0_wp * swb0(A2D(0)) * sustar(A2D(0))**2 * phml(A2D(0)) / & 3031 & ( svstr(A2D(0))**3 + 0.5_wp * swstrc(A2D(0))**3 + epsln ) 3032 zsc_uw_2(:,:) = swb0(A2D(0)) * sustke(A2D(0)) * phml(A2D(0)) / & 3033 & ( svstr(A2D(0))**3 + 0.5_wp * swstrc(A2D(0))**3 + epsln )**( 2.0_wp / 3.0_wp ) 3037 3034 ELSEWHERE 3038 3035 zsc_uw_1(:,:) = 0.0_wp … … 3059 3056 IF ( k_ddh(ji,jj) == 0 ) THEN 3060 3057 ! Place holding code. Parametrization needs checking for these conditions. 3061 zomega = ( 0.15_wp * pwstrl(ji,jj)**3 + pwstrc(ji,jj)**3 + 4.75_wp * ( pshear(ji,jj) * phbl(ji,jj) ) )**pthird3058 zomega = ( 0.15_wp * swstrl(ji,jj)**3 + swstrc(ji,jj)**3 + 4.75_wp * ( pshear(ji,jj) * phbl(ji,jj) ) )**pthird 3062 3059 zuw_bse(ji,jj) = -0.0035_wp * zomega * ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pdu_ml(ji,jj) 3063 3060 zvw_bse(ji,jj) = -0.0075_wp * zomega * ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pdv_ml(ji,jj) 3064 3061 ELSE 3065 zomega = ( 0.15_wp * pwstrl(ji,jj)**3 + pwstrc(ji,jj)**3 + 4.75_wp * ( pshear(ji,jj) * phbl(ji,jj) ) )**pthird3062 zomega = ( 0.15_wp * swstrl(ji,jj)**3 + swstrc(ji,jj)**3 + 4.75_wp * ( pshear(ji,jj) * phbl(ji,jj) ) )**pthird 3066 3063 zuw_bse(ji,jj) = -0.0035_wp * zomega * ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pdu_ml(ji,jj) 3067 3064 zvw_bse(ji,jj) = -0.0075_wp * zomega * ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pdv_ml(ji,jj) 3068 3065 ENDIF 3069 zb_cubic(ji,jj) = pdh(ji,jj) / phbl(ji,jj) * puw0(ji,jj) - ( 2.0 + pdh(ji,jj) / phml(ji,jj) ) * zuw_bse(ji,jj)3066 zb_cubic(ji,jj) = pdh(ji,jj) / phbl(ji,jj) * suw0(ji,jj) - ( 2.0 + pdh(ji,jj) / phml(ji,jj) ) * zuw_bse(ji,jj) 3070 3067 za_cubic(ji,jj) = zuw_bse(ji,jj) - zb_cubic(ji,jj) 3071 zvw_max = 0.7_wp * ff_t(ji,jj) * ( pustke(ji,jj) * dstokes(ji,jj) + 0.7_wp * pustar(ji,jj) * phml(ji,jj) )3068 zvw_max = 0.7_wp * ff_t(ji,jj) * ( sustke(ji,jj) * dstokes(ji,jj) + 0.7_wp * sustar(ji,jj) * phml(ji,jj) ) 3072 3069 zd_cubic(ji,jj) = zvw_max * pdh(ji,jj) / phml(ji,jj) - ( 2.0_wp + pdh(ji,jj) / phml(ji,jj) ) * zvw_bse(ji,jj) 3073 3070 zc_cubic(ji,jj) = zvw_bse(ji,jj) - zd_cubic(ji,jj) … … 3115 3112 ! ----------------------------------------------------------------------- 3116 3113 WHERE ( ldconv(A2D(0)) ) 3117 zsc_wth_1(:,:) = pwth0(A2D(0)) / ( 1.0_wp - 0.56_wp * EXP( phol(A2D(0)) ) )3118 zsc_ws_1(:,:) = pws0(A2D(0)) / ( 1.0_wp - 0.56_wp * EXP( phol(A2D(0)) ) )3114 zsc_wth_1(:,:) = swth0(A2D(0)) / ( 1.0_wp - 0.56_wp * EXP( shol(A2D(0)) ) ) 3115 zsc_ws_1(:,:) = sws0(A2D(0)) / ( 1.0_wp - 0.56_wp * EXP( shol(A2D(0)) ) ) 3119 3116 WHERE ( ldpyc(A2D(0)) ) ! Pycnocline scales 3120 zsc_wth_pyc(:,:) = -0.003_wp * pwstrc(A2D(0)) * ( 1.0_wp - pdh(A2D(0)) / phbl(A2D(0)) ) * pdt_ml(A2D(0))3121 zsc_ws_pyc(:,:) = -0.003_wp * pwstrc(A2D(0)) * ( 1.0_wp - pdh(A2D(0)) / phbl(A2D(0)) ) * pds_ml(A2D(0))3117 zsc_wth_pyc(:,:) = -0.003_wp * swstrc(A2D(0)) * ( 1.0_wp - pdh(A2D(0)) / phbl(A2D(0)) ) * pdt_ml(A2D(0)) 3118 zsc_ws_pyc(:,:) = -0.003_wp * swstrc(A2D(0)) * ( 1.0_wp - pdh(A2D(0)) / phbl(A2D(0)) ) * pds_ml(A2D(0)) 3122 3119 END WHERE 3123 3120 ELSEWHERE 3124 zsc_wth_1(:,:) = 2.0 * pwthav(A2D(0))3125 zsc_ws_1(:,:) = pws0(A2D(0))3121 zsc_wth_1(:,:) = 2.0 * swthav(A2D(0)) 3122 zsc_ws_1(:,:) = sws0(A2D(0)) 3126 3123 END WHERE 3127 3124 DO_3D( 0, 0, 0, 0, 1, MAX( jkm_mld, jkm_bld ) ) … … 3156 3153 ! 3157 3154 WHERE ( ldconv(A2D(0)) ) 3158 zsc_uw_1(:,:) = pustar(A2D(0))**23159 zsc_vw_1(:,:) = ff_t(A2D(0)) * pustke(A2D(0)) * phml(A2D(0))3155 zsc_uw_1(:,:) = sustar(A2D(0))**2 3156 zsc_vw_1(:,:) = ff_t(A2D(0)) * sustke(A2D(0)) * phml(A2D(0)) 3160 3157 ELSEWHERE 3161 zsc_uw_1(:,:) = pustar(A2D(0))**23158 zsc_uw_1(:,:) = sustar(A2D(0))**2 3162 3159 zsc_uw_2(:,:) = ( 2.25_wp - 3.0_wp * ( 1.0_wp - EXP( -1.25_wp * 2.0_wp ) ) ) * ( 1.0_wp - EXP( -4.0_wp * 2.0_wp ) ) * & 3163 3160 & zsc_uw_1(:,:) 3164 zsc_vw_1(:,:) = ff_t(A2D(0)) * pustke(A2D(0)) * phbl(A2D(0))3161 zsc_vw_1(:,:) = ff_t(A2D(0)) * sustke(A2D(0)) * phbl(A2D(0)) 3165 3162 zsc_vw_2(:,:) = -0.11_wp * SIN( 3.14159_wp * ( 2.0_wp + 0.4_wp ) ) * EXP( -1.0_wp * ( 1.5_wp + 2.0_wp )**2 ) * & 3166 3163 & zsc_vw_1(:,:) … … 3277 3274 IF ( pdhdt(ji,jj) > 0.0_wp ) THEN ! Depth increasing, or steady. 3278 3275 IF ( pdb_bl(ji,jj) > 0.0_wp ) THEN 3279 IF ( phol(ji,jj) >= 0.5_wp ) THEN ! Very stable - 'thick' pycnocline3276 IF ( shol(ji,jj) >= 0.5_wp ) THEN ! Very stable - 'thick' pycnocline 3280 3277 ztmp = 1.0_wp / MAX( phbl(ji,jj), epsln ) 3281 3278 ztgrad = pdt_bl(ji,jj) * ztmp … … 3309 3306 END IF 3310 3307 END IF 3311 ENDIF ! IF ( zhol >=0.5)3308 ENDIF ! IF (shol >=0.5) 3312 3309 ENDIF ! IF (zdb_bl> 0.) 3313 3310 ENDIF ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero
Note: See TracChangeset
for help on using the changeset viewer.
