Changeset 12216 for NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
- Timestamp:
- 2019-12-12T17:01:18+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
r12178 r12216 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 ouyancy flux due to entrainment changed to include contribution from shear turbulence (for testing commented out).27 !! (14) Buoyancy flux due to entrainment changed to include contribution from shear turbulence. 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 layers 33 !! (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_osmldpth1 31 36 !!---------------------------------------------------------------------- 32 37 … … 40 45 !! trc_osm : compute and add to the passive tracer trend the non-local flux (TBD) 41 46 !! dyn_osm : compute and add to u & v trensd the non-local flux 47 !! 48 !! Subroutines in revised code. 42 49 !!---------------------------------------------------------------------- 43 50 USE oce ! ocean dynamics and active tracers … … 97 104 98 105 ! !!! ** General constants ** 99 REAL(wp) :: epsln = 1.0e-20_wp ! a small positive number 106 REAL(wp) :: epsln = 1.0e-20_wp ! a small positive number to ensure no div by zero 107 REAL(wp) :: depth_tol = 1.0e-6_wp ! a small-ish positive number to give a hbl slightly shallower than gdepw 100 108 REAL(wp) :: pthird = 1._wp/3._wp ! 1/3 101 109 REAL(wp) :: p2third = 2._wp/3._wp ! 2/3 … … 166 174 REAL(wp) :: zbeta, zthermal ! 167 175 REAL(wp) :: zehat, zeta, zhrib, zsig, zscale, zwst, zws, zwm ! Velocity scales 168 REAL(wp) :: zwsun, zwmun, zcons, zconm, zwcons, zwconm ! 176 REAL(wp) :: zwsun, zwmun, zcons, zconm, zwcons, zwconm ! 177 169 178 REAL(wp) :: zsr, zbw, ze, zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zcomp , zrhd,zrhdr,zbvzed ! In situ density 170 179 INTEGER :: jm ! dummy loop indices … … 196 205 REAL(wp), DIMENSION(jpi,jpj) :: zsin_wind ! Sin angle of surface stress 197 206 REAL(wp), DIMENSION(jpi,jpj) :: zhol ! Stability parameter for boundary layer 198 LOGICAL, DIMENSION( :,:), ALLOCATABLE :: lconv! unstable/stable bl207 LOGICAL, DIMENSION(jpi,jpj) :: lconv ! unstable/stable bl 199 208 200 209 ! mixed-layer variables … … 238 247 ! Temporary variables 239 248 INTEGER :: inhml 240 INTEGER :: i_lconv_alloc241 249 REAL(wp) :: znd,znd_d,zznd_ml,zznd_pyc,zznd_d ! temporary non-dimensional depths used in various routines 242 250 REAL(wp) :: ztemp, zari, zpert, zzdhdt, zdb ! temporary variables … … 248 256 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdiffut ! t-diffusivity 249 257 258 INTEGER :: ibld_ext=0 ! does not have to be zero for modified scheme 259 REAL(wp) :: zwb_min, zgamma_b_nd, zgamma_b, zdhoh, ztau, zddhdt 260 REAL(wp) :: zzeta_s = 0._wp 261 REAL(wp) :: zzeta_v = 0.46 262 REAL(wp) :: zabsstke 263 250 264 ! For debugging 251 265 INTEGER :: ikt 252 266 !!-------------------------------------------------------------------- 253 267 ! 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 257 268 ibld(:,:) = 0 ; imld(:,:) = 0 258 269 zrad0(:,:) = 0._wp ; zradh(:,:) = 0._wp ; zradav(:,:) = 0._wp ; zustar(:,:) = 0._wp … … 268 279 zt_bl(:,:) = 0._wp ; zs_bl(:,:) = 0._wp ; zu_bl(:,:) = 0._wp ; zv_bl(:,:) = 0._wp 269 280 zrh_bl(:,:) = 0._wp ; zt_ml(:,:) = 0._wp ; zs_ml(:,:) = 0._wp ; zu_ml(:,:) = 0._wp 281 270 282 zv_ml(:,:) = 0._wp ; zrh_ml(:,:) = 0._wp ; zdt_bl(:,:) = 0._wp ; zds_bl(:,:) = 0._wp 271 283 zdu_bl(:,:) = 0._wp ; zdv_bl(:,:) = 0._wp ; zdrh_bl(:,:) = 0._wp ; zdb_bl(:,:) = 0._wp … … 287 299 ghams(:,:,:) = 0._wp ; ghamu(:,:,:) = 0._wp ; ghamv(:,:,:) = 0._wp 288 300 301 zdhdt_2(:,:) = 0._wp 289 302 ! hbl = MAX(hbl,epsln) 290 303 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 350 363 ! Use wind speed wndm included in sbc_oce module 351 364 zustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8 ) 352 dstokes(ji,jj) = 0.12 * wndm(ji,jj)**2 / grav365 dstokes(ji,jj) = MAX( 0.12 * wndm(ji,jj)**2 / grav, 5.e-1) 353 366 END DO 354 367 END DO … … 362 375 ! It could represent the effects of the spread of wave directions 363 376 ! around the mean wind. The effect of this adjustment needs to be tested. 364 z ustke(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(z ustke(ji,jj)*wmp(ji,jj), 1.0e-7 ) ), 5.0e-1) !rn_osm_dstokes !377 zabsstke = SQRT(ut0sd(ji,jj)**2 + vt0sd(ji,jj)**2) 378 zustke(ji,jj) = MAX (0.8 * ( zcos_wind(ji,jj) * ut0sd(ji,jj) + zsin_wind(ji,jj) * vt0sd(ji,jj) ), 1.0e-8) 379 dstokes(ji,jj) = MAX(zfac * hsw(ji,jj)*hsw(ji,jj) / ( MAX(zabsstke*wmp(ji,jj), 1.0e-7 ) ), 5.0e-1) !rn_osm_dstokes ! 367 380 END DO 368 381 END DO … … 375 388 ! Langmuir velocity scale (zwstrl), at T-point 376 389 zwstrl(ji,jj) = ( zustar(ji,jj) * zustar(ji,jj) * zustke(ji,jj) )**pthird 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 390 zla(ji,jj) = MAX(MIN(SQRT ( zustar(ji,jj) / ( zwstrl(ji,jj) + epsln ) )**3, 4.0), 0.2) 391 IF(zla(ji,jj) > 0.45) dstokes(ji,jj) = MIN(dstokes(ji,jj), 0.5_wp*hbl(ji,jj)) 385 392 ! Velocity scale that tends to zustar for large Langmuir numbers 386 393 zvstr(ji,jj) = ( zwstrl(ji,jj)**3 + & … … 389 396 ! limit maximum value of Langmuir number as approximate treatment for shear turbulence. 390 397 ! Note zustke and zwstrl are not amended. 391 IF ( zla(ji,jj) >= 0.45 ) zla(ji,jj) = 0.45392 398 ! 393 399 ! get convective velocity (zwstrc), stabilty scale (zhol) and logical conection flag lconv … … 406 412 ! Mixed-layer model - calculate averages over the boundary layer, and the change in the boundary layer depth 407 413 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 408 ! BL must be always 2levels deep.409 hbl(:,:) = MAX(hbl(:,:), gdepw_n(:,:, 3) )410 ibld(:,:) = 3411 DO jk = 4, jpkm1414 ! BL must be always 4 levels deep. 415 hbl(:,:) = MAX(hbl(:,:), gdepw_n(:,:,4) ) 416 ibld(:,:) = 4 417 DO jk = 5, jpkm1 412 418 DO jj = 2, jpjm1 413 419 DO ji = 2, jpim1 … … 419 425 END DO 420 426 421 DO jj = 2, jpjm1 ! Vertical slab427 DO jj = 2, jpjm1 422 428 DO ji = 2, jpim1 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 429 zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj)) 430 imld(ji,jj) = MAX(3,ibld(ji,jj) - MAX( INT( dh(ji,jj) / e3t_n(ji, jj, ibld(ji,jj) )) , 1 )) 431 zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 432 END DO 483 433 END DO 484 434 … … 495 445 DO ji = 2, jpim1 496 446 IF ( zhbl_t(ji,jj) >= gdepw_n(ji,jj,jk) ) THEN 497 ibld(ji,jj) = MIN(mbkt(ji,jj), jk)447 ibld(ji,jj) = jk 498 448 ENDIF 499 449 END DO … … 504 454 ! Step through model levels taking account of buoyancy change to determine the effect on dhdt 505 455 ! 506 DO jj = 2, jpjm1507 DO ji = 2, jpim1508 IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN456 CALL zdf_osm_timestep_hbl( zdhdt, zdhdt_2 ) 457 ! Alan: do we need zb_ml? 458 CALL zdf_osm_vertical_average( ibld, zt_bl, zs_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl ) 509 459 ! 510 ! If boundary layer changes by more than one level, need to check for stable layers between initial and final depths.511 460 ! 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 461 CALL zdf_osm_pycnocline_thickness( dh, zdh ) 567 462 dstokes(:,:) = MIN ( dstokes(:,:), hbl(:,:)/3. ) ! Limit delta for shallow boundary layers for calculating flux-gradient terms. 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 ! 463 ! 464 ! Average over the depth of the mixed layer in the convective boundary layer 465 ! Alan: do we need zb_ml? 466 CALL zdf_osm_vertical_average( imld, zt_ml, zs_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml ) 724 467 ! rotate mean currents and changes onto wind align co-ordinates 725 468 ! 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 469 CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_ml, zv_ml, zdu_ml, zdv_ml ) 470 CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_bl, zv_bl, zdu_bl, zdv_bl ) 745 471 zuw_bse = 0._wp 746 472 zvw_bse = 0._wp 473 zwth_ent = 0._wp 474 zws_ent = 0._wp 747 475 DO jj = 2, jpjm1 748 476 DO ji = 2, jpim1 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) 477 IF ( ibld(ji,jj) < mbkt(ji,jj) ) THEN 478 IF ( lconv(ji,jj) ) THEN 479 zuw_bse(ji,jj) = -0.0075*((zvstr(ji,jj)**3+0.5*zwstrc(ji,jj)**3)**pthird*zdu_ml(ji,jj) + & 480 & 1.5*zustar(ji,jj)**2*(zhbl(ji,jj)-zhml(ji,jj)) )/ & 481 & ( zhml(ji,jj)*MIN(zla(ji,jj)**(8./3.),1.) + epsln) 482 zvw_bse(ji,jj) = 0.01*(-(zvstr(ji,jj)**3+0.5*zwstrc(ji,jj)**3)**pthird*zdv_ml(ji,jj)+ & 483 & 2.0*ff_t(ji,jj)*zustke(ji,jj)*dstokes(ji,jj)*zla(ji,jj)) 484 IF ( zdb_ml(ji,jj) > 0._wp ) THEN 485 zwth_ent(ji,jj) = zwb_ent(ji,jj) * zdt_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 486 zws_ent(ji,jj) = zwb_ent(ji,jj) * zds_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 487 ENDIF 488 ELSE 489 zwth_ent(ji,jj) = -2.0 * zwthav(ji,jj) * ( (1.0 - 0.8) - ( 1.0 - 0.8)**(3.0/2.0) ) 490 zws_ent(ji,jj) = -2.0 * zwsav(ji,jj) * ( (1.0 - 0.8 ) - ( 1.0 - 0.8 )**(3.0/2.0) ) 754 491 ENDIF 755 ELSE756 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) )758 492 ENDIF 759 493 END DO … … 764 498 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 765 499 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 500 CALL zdf_osm_external_gradients( zdtdz_ext, zdsdz_ext, zdbdz_ext ) 501 CALL zdf_osm_pycnocline_scalar_profiles( zdtdz_pyc, zdsdz_pyc, zdbdz_pyc ) 502 CALL zdf_osm_pycnocline_shear_profiles( zdudz_pyc, zdvdz_pyc ) 847 503 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 848 504 ! Eddy viscosity/diffusivity and non-gradient terms in the flux-gradient relationship 849 505 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 850 506 851 ! WHERE ( lconv )852 ! zdifml_sc = zhml * ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird853 ! zvisml_sc = zdifml_sc854 ! 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 )**p2third857 ! zbeta_v_sc = 1.0 - 2.0 * (0.142 /0.375) * (zhbl - zhml ) / zhml858 ! ELSEWHERE859 ! zdifml_sc = zwstrl * zhbl * EXP ( -( zhol / 0.183_wp )**2 )860 ! zvisml_sc = zwstrl * zhbl * EXP ( -( zhol / 0.183_wp )**2 )861 ! ENDWHERE862 507 DO jj = 2, jpjm1 863 508 DO ji = 2, jpim1 … … 868 513 zvispyc_sc(ji,jj) = 0.142 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zdh(ji,jj) 869 514 zbeta_d_sc(ji,jj) = 1.0 - (0.165 / 0.8 * zdh(ji,jj) / zhbl(ji,jj) )**p2third 870 zbeta_v_sc(ji,jj) = 1.0 - 2.0 * (0.142 /0.375) * zdh(ji,jj) / zhml(ji,jj)515 zbeta_v_sc(ji,jj) = 1.0 - 2.0 * (0.142 /0.375) * zdh(ji,jj) / ( zhml(ji,jj) + epsln ) 871 516 ELSE 872 517 zdifml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ) 873 518 zvisml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ) 874 END IF875 END DO876 END DO519 END IF 520 END DO 521 END DO 877 522 ! 878 523 DO jj = 2, jpjm1 … … 889 534 ! pycnocline - if present linear profile 890 535 IF ( zdh(ji,jj) > 0._wp ) THEN 536 zgamma_b = 6.0 891 537 DO jk = imld(ji,jj)+1 , ibld(ji,jj) 892 538 zznd_pyc = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 893 539 ! 894 zdiffut(ji,jj,jk) = zdifpyc_sc(ji,jj) * ( 1.0 +zznd_pyc )540 zdiffut(ji,jj,jk) = zdifpyc_sc(ji,jj) * EXP( zgamma_b * zznd_pyc ) 895 541 ! 896 zviscos(ji,jj,jk) = zvispyc_sc(ji,jj) * ( 1.0 +zznd_pyc )542 zviscos(ji,jj,jk) = zvispyc_sc(ji,jj) * EXP( zgamma_b * zznd_pyc ) 897 543 END DO 544 IF ( ibld_ext == 0 ) THEN 545 zdiffut(ji,jj,ibld(ji,jj)) = 0._wp 546 zviscos(ji,jj,ibld(ji,jj)) = 0._wp 547 ELSE 548 zdiffut(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj) * ( hbl(ji,jj) - gdepw_n(ji, jj, ibld(ji,jj)-1) ) 549 zviscos(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj) * ( hbl(ji,jj) - gdepw_n(ji, jj, ibld(ji,jj)-1) ) 550 ENDIF 898 551 ENDIF 899 552 ! Temporay fix to ensure zdiffut is +ve; won't be necessary with wn taken out … … 908 561 zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * (1.0 - zznd_ml) * ( 1.0 - zznd_ml**2 ) 909 562 END DO 563 564 IF ( ibld_ext == 0 ) THEN 565 zdiffut(ji,jj,ibld(ji,jj)) = 0._wp 566 zviscos(ji,jj,ibld(ji,jj)) = 0._wp 567 ELSE 568 zdiffut(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 0._wp) * e3w_n(ji, jj, ibld(ji,jj)) 569 zviscos(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 0._wp) * e3w_n(ji, jj, ibld(ji,jj)) 570 ENDIF 910 571 ENDIF ! end if ( lconv ) 911 !572 ! 912 573 END DO ! end of ji loop 913 574 END DO ! end of jj loop … … 952 613 END DO ! end of jj loop 953 614 954 955 615 ! 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) 956 616 WHERE ( lconv ) 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 )617 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 ) 618 zsc_uw_2 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke / MIN( zla**(8.0/3.0) + epsln, 0.12 ) 619 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 ) 960 620 ELSEWHERE 961 621 zsc_uw_1 = zustar**2 962 zsc_vw_1 = ff_t * zhbl * zustke**3 * zla**(8.0/3.0) / (zvstr**2 + epsln)622 zsc_vw_1 = ff_t * zhbl * zustke**3 * MIN( zla**(8.0/3.0), 0.12 ) / (zvstr**2 + epsln) 963 623 ENDWHERE 964 624 IF(ln_dia_osm) THEN 625 IF ( iom_use("ghamu_00") ) CALL iom_put( "ghamu_00", wmask*ghamu ) 626 IF ( iom_use("ghamv_00") ) CALL iom_put( "ghamv_00", wmask*ghamv ) 627 END IF 965 628 DO jj = 2, jpjm1 966 629 DO ji = 2, jpim1 … … 1007 670 zl_l = 2.0 * ( 1.0 - EXP ( - 2.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) ) & 1008 671 & * ( 1.0 - EXP ( - 5.0 * ( 1.0 - zznd_ml ) ) ) * ( 1.0 + dstokes(ji,jj) / zhml (ji,jj) ) 1009 zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0 + EXP ( 3.0 * LOG10 ( - zhol(ji,jj) ) ) ) ** (3.0/2.0)672 zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0 + EXP ( -3.0 * LOG10 ( - zhol(ji,jj) ) ) ) ** (3.0/2.0) 1010 673 ! non-gradient buoyancy terms 1011 674 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 ) … … 1020 683 END DO ! ji loop 1021 684 END DO ! jj loop 1022 1023 685 1024 686 WHERE ( lconv ) … … 1051 713 END DO ! jj loop 1052 714 715 IF(ln_dia_osm) THEN 716 IF ( iom_use("ghamu_0") ) CALL iom_put( "ghamu_0", wmask*ghamu ) 717 IF ( iom_use("zsc_uw_1_0") ) CALL iom_put( "zsc_uw_1_0", tmask(:,:,1)*zsc_uw_1 ) 718 END IF 1053 719 ! Transport term in flux-gradient relationship [note : includes ROI ratio (X0.3) ] 1054 720 … … 1089 755 END DO ! jj loop 1090 756 1091 1092 757 WHERE ( lconv ) 1093 758 zsc_uw_1 = zustar**2 … … 1134 799 END DO 1135 800 END DO 801 802 IF(ln_dia_osm) THEN 803 IF ( iom_use("ghamu_f") ) CALL iom_put( "ghamu_f", wmask*ghamu ) 804 IF ( iom_use("ghamv_f") ) CALL iom_put( "ghamv_f", wmask*ghamv ) 805 IF ( iom_use("zsc_uw_1_f") ) CALL iom_put( "zsc_uw_1_f", tmask(:,:,1)*zsc_uw_1 ) 806 IF ( iom_use("zsc_vw_1_f") ) CALL iom_put( "zsc_vw_1_f", tmask(:,:,1)*zsc_vw_1 ) 807 IF ( iom_use("zsc_uw_2_f") ) CALL iom_put( "zsc_uw_2_f", tmask(:,:,1)*zsc_uw_2 ) 808 IF ( iom_use("zsc_vw_2_f") ) CALL iom_put( "zsc_vw_2_f", tmask(:,:,1)*zsc_vw_2 ) 809 END IF 1136 810 ! 1137 811 ! Make surface forced velocity non-gradient terms go to zero at the base of the mixed layer. … … 1165 839 END DO 1166 840 841 IF(ln_dia_osm) THEN 842 IF ( iom_use("ghamu_b") ) CALL iom_put( "ghamu_b", wmask*ghamu ) 843 IF ( iom_use("ghamv_b") ) CALL iom_put( "ghamv_b", wmask*ghamv ) 844 END IF 1167 845 ! pynocline contributions 1168 846 ! Temporary fix to avoid instabilities when zdb_bl becomes very very small … … 1170 848 DO jj = 2, jpjm1 1171 849 DO ji = 2, jpim1 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 850 IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN 851 DO jk= 2, ibld(ji,jj) 852 znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 853 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk) 854 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk) 855 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc(ji,jj,jk) 856 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) 857 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc(ji,jj,jk) 858 END DO 859 END IF 860 END DO 1181 861 END DO 1182 862 … … 1185 865 DO jj=2, jpjm1 1186 866 DO ji = 2, jpim1 1187 IF ( lconv(ji,jj)) THEN867 IF ( lconv(ji,jj) .AND. ibld(ji,jj) + ibld_ext < mbkt(ji,jj)) THEN 1188 868 DO jk = 1, imld(ji,jj) - 1 1189 869 znd=gdepw_n(ji,jj,jk) / zhml(ji,jj) 1190 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * znd1191 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * znd870 ! ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * znd 871 ! ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * znd 1192 872 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * znd 1193 873 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * znd … … 1195 875 DO jk = imld(ji,jj), ibld(ji,jj) 1196 876 znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 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 )877 ! ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * ( 1.0 + znd ) 878 ! ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * ( 1.0 + znd ) 1199 879 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * ( 1.0 + znd ) 1200 880 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * ( 1.0 + znd ) 1201 881 END DO 1202 882 ENDIF 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 883 884 ghamt(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 885 ghams(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 886 ghamu(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 887 ghamv(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 1207 888 END DO ! ji loop 1208 889 END DO ! jj loop 1209 890 1210 891 IF(ln_dia_osm) THEN 892 IF ( iom_use("ghamu_1") ) CALL iom_put( "ghamu_1", wmask*ghamu ) 893 IF ( iom_use("ghamv_1") ) CALL iom_put( "ghamv_1", wmask*ghamv ) 894 IF ( iom_use("zuw_bse") ) CALL iom_put( "zuw_bse", tmask(:,:,1)*zuw_bse ) 895 IF ( iom_use("zvw_bse") ) CALL iom_put( "zvw_bse", tmask(:,:,1)*zvw_bse ) 896 IF ( iom_use("zdudz_pyc") ) CALL iom_put( "zdudz_pyc", wmask*zdudz_pyc ) 897 IF ( iom_use("zdvdz_pyc") ) CALL iom_put( "zdvdz_pyc", wmask*zdvdz_pyc ) 898 IF ( iom_use("zviscos") ) CALL iom_put( "zviscos", wmask*zviscos ) 899 END IF 1211 900 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 1212 901 ! Need to put in code for contributions that are applied explicitly to … … 1287 976 1288 977 ! Lateral boundary conditions on zvicos (sign unchanged), needed to caclulate viscosities on u and v grids 1289 CALL lbc_lnk( 'zdfosm',zviscos(:,:,:), 'W', 1. )978 !CALL lbc_lnk( zviscos(:,:,:), 'W', 1. ) 1290 979 1291 980 ! GN 25/8: need to change tmask --> wmask … … 1319 1008 ! Lateral boundary conditions on final outputs for gham[uv], on [UV]-grid (sign unchanged) 1320 1009 CALL lbc_lnk_multi( 'zdfosm', ghamt, 'W', 1. , ghams, 'W', 1., & 1321 & ghamu, 'U', 1. , ghamv, 'V',1. )1322 1323 1010 & ghamu, 'U', -1. , ghamv, 'V', -1. ) 1011 1012 IF(ln_dia_osm) THEN 1324 1013 SELECT CASE (nn_osm_wave) 1325 1014 ! Stokes drift set by assumimg onstant La#=0.3(=0) or Pierson-Moskovitz spectrum (=1). … … 1330 1019 ! Stokes drift read in from sbcwave (=2). 1331 1020 CASE(2) 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 1021 IF ( iom_use("us_x") ) CALL iom_put( "us_x", ut0sd*umask(:,:,1) ) ! x surface Stokes drift 1022 IF ( iom_use("us_y") ) CALL iom_put( "us_y", vt0sd*vmask(:,:,1) ) ! y surface Stokes drift 1023 IF ( iom_use("wmp") ) CALL iom_put( "wmp", wmp*tmask(:,:,1) ) ! wave mean period 1024 IF ( iom_use("hsw") ) CALL iom_put( "hsw", hsw*tmask(:,:,1) ) ! significant wave height 1025 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 1026 IF ( iom_use("hsw_NP") ) CALL iom_put( "hsw_NP", (0.22/grav)*wndm**2*tmask(:,:,1) ) ! significant wave height from NP spectrum 1027 IF ( iom_use("wndm") ) CALL iom_put( "wndm", wndm*tmask(:,:,1) ) ! U_10 1334 1028 IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rau0*tmask(:,:,1)*zustar**2* & 1335 1029 & SQRT(ut0sd**2 + vt0sd**2 ) ) … … 1342 1036 IF ( iom_use("zws0") ) CALL iom_put( "zws0", tmask(:,:,1)*zws0 ) ! <Sw_0> 1343 1037 IF ( iom_use("hbl") ) CALL iom_put( "hbl", tmask(:,:,1)*hbl ) ! boundary-layer depth 1344 IF ( iom_use("hbli") ) CALL iom_put( "hbli", tmask(:,:,1)*hbli ) ! Initial boundary-layer depth 1038 IF ( iom_use("ibld") ) CALL iom_put( "ibld", tmask(:,:,1)*ibld ) ! boundary-layer max k 1039 IF ( iom_use("zdt_bl") ) CALL iom_put( "zdt_bl", tmask(:,:,1)*zdt_bl ) ! dt at ml base 1040 IF ( iom_use("zds_bl") ) CALL iom_put( "zds_bl", tmask(:,:,1)*zds_bl ) ! ds at ml base 1041 IF ( iom_use("zdb_bl") ) CALL iom_put( "zdb_bl", tmask(:,:,1)*zdb_bl ) ! db at ml base 1042 IF ( iom_use("zdu_bl") ) CALL iom_put( "zdu_bl", tmask(:,:,1)*zdu_bl ) ! du at ml base 1043 IF ( iom_use("zdv_bl") ) CALL iom_put( "zdv_bl", tmask(:,:,1)*zdv_bl ) ! dv at ml base 1044 IF ( iom_use("dh") ) CALL iom_put( "dh", tmask(:,:,1)*dh ) ! Initial boundary-layer depth 1045 IF ( iom_use("hml") ) CALL iom_put( "hml", tmask(:,:,1)*hml ) ! Initial boundary-layer depth 1345 1046 IF ( iom_use("dstokes") ) CALL iom_put( "dstokes", tmask(:,:,1)*dstokes ) ! Stokes drift penetration depth 1346 1047 IF ( iom_use("zustke") ) CALL iom_put( "zustke", tmask(:,:,1)*zustke ) ! Stokes drift magnitude at T-points … … 1348 1049 IF ( iom_use("zwstrl") ) CALL iom_put( "zwstrl", tmask(:,:,1)*zwstrl ) ! Langmuir velocity scale 1349 1050 IF ( iom_use("zustar") ) CALL iom_put( "zustar", tmask(:,:,1)*zustar ) ! friction velocity scale 1051 IF ( iom_use("zvstr") ) CALL iom_put( "zvstr", tmask(:,:,1)*zvstr ) ! mixed velocity scale 1052 IF ( iom_use("zla") ) CALL iom_put( "zla", tmask(:,:,1)*zla ) ! langmuir # 1350 1053 IF ( iom_use("wind_power") ) CALL iom_put( "wind_power", 1000.*rau0*tmask(:,:,1)*zustar**3 ) ! BL depth internal to zdf_osm routine 1351 1054 IF ( iom_use("wind_wave_power") ) CALL iom_put( "wind_wave_power", 1000.*rau0*tmask(:,:,1)*zustar**2*zustke ) 1352 1055 IF ( iom_use("zhbl") ) CALL iom_put( "zhbl", tmask(:,:,1)*zhbl ) ! BL depth internal to zdf_osm routine 1353 1056 IF ( iom_use("zhml") ) CALL iom_put( "zhml", tmask(:,:,1)*zhml ) ! ML depth internal to zdf_osm routine 1354 IF ( iom_use("zdh") ) CALL iom_put( "zdh", tmask(:,:,1)*zdh ) ! ML depth internal to zdf_osm routine 1057 IF ( iom_use("imld") ) CALL iom_put( "imld", tmask(:,:,1)*imld ) ! index for ML depth internal to zdf_osm routine 1058 IF ( iom_use("zdh") ) CALL iom_put( "zdh", tmask(:,:,1)*zdh ) ! pyc thicknessh internal to zdf_osm routine 1355 1059 IF ( iom_use("zhol") ) CALL iom_put( "zhol", tmask(:,:,1)*zhol ) ! ML depth internal to zdf_osm routine 1356 1060 IF ( iom_use("zwthav") ) CALL iom_put( "zwthav", tmask(:,:,1)*zwthav ) ! ML depth internal to zdf_osm routine … … 1378 1082 INTEGER :: ios ! local integer 1379 1083 INTEGER :: ji, jj, jk ! dummy loop indices 1084 REAL z1_t2 1380 1085 !! 1381 1086 NAMELIST/namzdf_osm/ ln_use_osm_la, rn_osm_la, rn_osm_dstokes, nn_ave & … … 1401 1106 WRITE(numout,*) ' Turbulent Langmuir number rn_osm_la = ', rn_osm_la 1402 1107 WRITE(numout,*) ' Initial hbl for 1D runs rn_osm_hbl0 = ', rn_osm_hbl0 1403 WRITE(numout,*) ' Depth scale of Stokes drift rn_osm_dstokes = ', rn_osm_dstokes1108 WRITE(numout,*) ' Depth scale of Stokes drift rn_osm_dstokes = ', rn_osm_dstokes 1404 1109 WRITE(numout,*) ' horizontal average flag nn_ave = ', nn_ave 1405 1110 WRITE(numout,*) ' Stokes drift nn_osm_wave = ', nn_osm_wave … … 1536 1241 REAL(wp) :: zN2_c ! local scalar 1537 1242 REAL(wp) :: rho_c = 0.01_wp !: density criterion for mixed layer depth 1538 INTEGER, DIMENSION( :,:), ALLOCATABLE:: imld_rst ! level of mixed-layer depth (pycnocline top)1243 INTEGER, DIMENSION(jpi,jpj) :: imld_rst ! level of mixed-layer depth (pycnocline top) 1539 1244 !!---------------------------------------------------------------------- 1540 1245 ! … … 1551 1256 WRITE(numout,*) ' ===>>>> : wn not in restart file, set to zero initially' 1552 1257 END IF 1258 1553 1259 id1 = iom_varid( numror, 'hbl' , ldstop = .FALSE. ) 1554 id2 = iom_varid( numror, ' hbli' , ldstop = .FALSE. )1260 id2 = iom_varid( numror, 'dh' , ldstop = .FALSE. ) 1555 1261 IF( id1 > 0 .AND. id2 > 0) THEN ! 'hbl' exists; read and return 1556 1262 CALL iom_get( numror, jpdom_autoglo, 'hbl' , hbl , ldxios = lrxios ) 1557 CALL iom_get( numror, jpdom_autoglo, ' hbli', hbli, ldxios = lrxios )1558 WRITE(numout,*) ' ===>>>> : hbl & hbliread from restart file'1263 CALL iom_get( numror, jpdom_autoglo, 'dh', dh, ldxios = lrxios ) 1264 WRITE(numout,*) ' ===>>>> : hbl & dh read from restart file' 1559 1265 RETURN 1560 ELSE ! 'hbl' & ' hbli' not in restart file, recalculate1266 ELSE ! 'hbl' & 'dh' not in restart file, recalculate 1561 1267 WRITE(numout,*) ' ===>>>> : previous run without osmosis scheme, hbl computed from stratification' 1562 1268 END IF … … 1570 1276 CALL iom_rstput( kt, nitrst, numrow, 'wn' , wn , ldxios = lwxios ) 1571 1277 CALL iom_rstput( kt, nitrst, numrow, 'hbl' , hbl , ldxios = lwxios ) 1572 CALL iom_rstput( kt, nitrst, numrow, ' hbli' , hbli, ldxios = lwxios )1278 CALL iom_rstput( kt, nitrst, numrow, 'dh' , dh, ldxios = lwxios ) 1573 1279 RETURN 1574 1280 END IF … … 1578 1284 !!----------------------------------------------------------------------------- 1579 1285 IF( lwp ) WRITE(numout,*) ' ===>>>> : calculating hbl computed from stratification' 1580 ALLOCATE( imld_rst(jpi,jpj) )1581 1286 ! w-level of the mixing and mixed layers 1582 1287 CALL eos_rab( tsn, rab_n ) … … 1599 1304 DO jj = 1, jpj 1600 1305 DO ji = 1, jpi 1601 iiki = imld_rst(ji,jj) 1602 hbl (ji,jj) = gdepw_n(ji,jj,iiki ) * ssmask(ji,jj) ! Turbocline depth 1306 iiki = MAX(4,imld_rst(ji,jj)) 1307 hbl (ji,jj) = gdepw_n(ji,jj,iiki ) ! Turbocline depth 1308 dh (ji,jj) = e3t_n(ji,jj,iiki-1 ) ! Turbocline depth 1603 1309 END DO 1604 1310 END DO … … 1607 1313 DEALLOCATE( imld_rst ) 1608 1314 WRITE(numout,*) ' ===>>>> : hbl computed from stratification' 1315 wn(:,:,:) = 0._wp 1316 WRITE(numout,*) ' ===>>>> : wn not in restart file, set to zero initially' 1609 1317 END SUBROUTINE osm_rst 1610 1318 … … 1634 1342 ENDIF 1635 1343 1636 ! add non-local temperature and salinity flux1637 1344 DO jk = 1, jpkm1 1638 1345 DO jj = 2, jpjm1 … … 1648 1355 END DO 1649 1356 1650 1651 ! save the non-local tracer flux trends for diagnostic 1357 ! save the non-local tracer flux trends for diagnostics 1652 1358 IF( l_trdtra ) THEN 1653 1359 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 1654 1360 ztrds(:,:,:) = tsa(:,:,:,jp_sal) - 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 )1361 1362 CALL trd_tra( kt, 'TRA', jp_tem, jptra_osm, ztrdt ) 1363 CALL trd_tra( kt, 'TRA', jp_sal, jptra_osm, ztrds ) 1658 1364 DEALLOCATE( ztrdt ) ; DEALLOCATE( ztrds ) 1659 1365 ENDIF … … 1723 1429 1724 1430 !!====================================================================== 1431 1725 1432 END MODULE zdfosm
Note: See TracChangeset
for help on using the changeset viewer.
