- Timestamp:
- 2020-05-14T21:46:00+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/OCE/ZDF/zdfosm.F90 (modified) (65 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
- Property svn:externals
-
old new 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 9 # SETTE 10 ^/utils/CI/sette@HEAD sette
-
- Property svn:externals
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/ZDF/zdfosm.F90
r12384 r12928 49 49 !!---------------------------------------------------------------------- 50 50 USE oce ! ocean dynamics and active tracers 51 ! uses w nfrom previous time step (which is now wb) to calculate hbl51 ! uses ww from previous time step (which is now wb) to calculate hbl 52 52 USE dom_oce ! ocean space and time domain 53 53 USE zdf_oce ! ocean vertical physics … … 139 139 INTEGER :: idebug = 236 140 140 INTEGER :: jdebug = 228 141 !! * Substitutions 142 # include "do_loop_substitute.h90" 141 143 !!---------------------------------------------------------------------- 142 144 !! NEMO/OCE 4.0 , NEMO Consortium (2018) 143 !! $Id : zdfosm.F90 12317 2020-01-14 12:40:47Z agn$145 !! $Id$ 144 146 !! Software governed by the CeCILL license (see ./LICENSE) 145 147 !!---------------------------------------------------------------------- … … 171 173 172 174 173 SUBROUTINE zdf_osm( kt, p_avm, p_avt )175 SUBROUTINE zdf_osm( kt, Kbb, Kmm, Krhs, p_avm, p_avt ) 174 176 !!---------------------------------------------------------------------- 175 177 !! *** ROUTINE zdf_osm *** … … 206 208 !! the equation number. (LMD94, here after) 207 209 !!---------------------------------------------------------------------- 208 INTEGER , INTENT(in ) :: kt ! ocean time step 210 INTEGER , INTENT(in ) :: kt ! ocean time step 211 INTEGER , INTENT(in ) :: Kbb, Kmm, Krhs ! ocean time level indices 209 212 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: p_avm, p_avt ! momentum and tracer Kz (w-points) 210 213 !! … … 378 381 zz0 = rn_abs ! surface equi-partition in 2-bands 379 382 zz1 = 1. - rn_abs 380 DO jj = 2, jpjm1 381 DO ji = 2, jpim1 382 ! Surface downward irradiance (so always +ve) 383 zrad0(ji,jj) = qsr(ji,jj) * r1_rau0_rcp 384 ! Downwards irradiance at base of boundary layer 385 zradh(ji,jj) = zrad0(ji,jj) * ( zz0 * EXP( -hbl(ji,jj)/rn_si0 ) + zz1 * EXP( -hbl(ji,jj)/rn_si1) ) 386 ! Downwards irradiance averaged over depth of the OSBL 387 zradav(ji,jj) = zrad0(ji,jj) * ( zz0 * ( 1.0 - EXP( -hbl(ji,jj)/rn_si0 ) )*rn_si0 & 388 & + zz1 * ( 1.0 - EXP( -hbl(ji,jj)/rn_si1 ) )*rn_si1 ) / hbl(ji,jj) 389 END DO 390 END DO 383 DO_2D_00_00 384 ! Surface downward irradiance (so always +ve) 385 zrad0(ji,jj) = qsr(ji,jj) * r1_rho0_rcp 386 ! Downwards irradiance at base of boundary layer 387 zradh(ji,jj) = zrad0(ji,jj) * ( zz0 * EXP( -hbl(ji,jj)/rn_si0 ) + zz1 * EXP( -hbl(ji,jj)/rn_si1) ) 388 ! Downwards irradiance averaged over depth of the OSBL 389 zradav(ji,jj) = zrad0(ji,jj) * ( zz0 * ( 1.0 - EXP( -hbl(ji,jj)/rn_si0 ) )*rn_si0 & 390 & + zz1 * ( 1.0 - EXP( -hbl(ji,jj)/rn_si1 ) )*rn_si1 ) / hbl(ji,jj) 391 END_2D 391 392 ! Turbulent surface fluxes and fluxes averaged over depth of the OSBL 392 DO jj = 2, jpjm1 393 DO ji = 2, jpim1 394 zthermal = rab_n(ji,jj,1,jp_tem) 395 zbeta = rab_n(ji,jj,1,jp_sal) 396 ! Upwards surface Temperature flux for non-local term 397 zwth0(ji,jj) = - qns(ji,jj) * r1_rau0_rcp * tmask(ji,jj,1) 398 ! Upwards surface salinity flux for non-local term 399 zws0(ji,jj) = - ( ( emp(ji,jj)-rnf(ji,jj) ) * tsn(ji,jj,1,jp_sal) + sfx(ji,jj) ) * r1_rau0 * tmask(ji,jj,1) 400 ! Non radiative upwards surface buoyancy flux 401 zwb0(ji,jj) = grav * zthermal * zwth0(ji,jj) - grav * zbeta * zws0(ji,jj) 402 ! turbulent heat flux averaged over depth of OSBL 403 zwthav(ji,jj) = 0.5 * zwth0(ji,jj) - ( 0.5*( zrad0(ji,jj) + zradh(ji,jj) ) - zradav(ji,jj) ) 404 ! turbulent salinity flux averaged over depth of the OBSL 405 zwsav(ji,jj) = 0.5 * zws0(ji,jj) 406 ! turbulent buoyancy flux averaged over the depth of the OBSBL 407 zwbav(ji,jj) = grav * zthermal * zwthav(ji,jj) - grav * zbeta * zwsav(ji,jj) 408 ! Surface upward velocity fluxes 409 zuw0(ji,jj) = -utau(ji,jj) * r1_rau0 * tmask(ji,jj,1) 410 zvw0(ji,jj) = -vtau(ji,jj) * r1_rau0 * tmask(ji,jj,1) 411 ! Friction velocity (zustar), at T-point : LMD94 eq. 2 412 zustar(ji,jj) = MAX( SQRT( SQRT( zuw0(ji,jj) * zuw0(ji,jj) + zvw0(ji,jj) * zvw0(ji,jj) ) ), 1.0e-8 ) 413 zcos_wind(ji,jj) = -zuw0(ji,jj) / ( zustar(ji,jj) * zustar(ji,jj) ) 414 zsin_wind(ji,jj) = -zvw0(ji,jj) / ( zustar(ji,jj) * zustar(ji,jj) ) 415 END DO 416 END DO 393 DO_2D_00_00 394 zthermal = rab_n(ji,jj,1,jp_tem) 395 zbeta = rab_n(ji,jj,1,jp_sal) 396 ! Upwards surface Temperature flux for non-local term 397 zwth0(ji,jj) = - qns(ji,jj) * r1_rho0_rcp * tmask(ji,jj,1) 398 ! Upwards surface salinity flux for non-local term 399 zws0(ji,jj) = - ( ( emp(ji,jj)-rnf(ji,jj) ) * ts(ji,jj,1,jp_sal,Kmm) + sfx(ji,jj) ) * r1_rho0 * tmask(ji,jj,1) 400 ! Non radiative upwards surface buoyancy flux 401 zwb0(ji,jj) = grav * zthermal * zwth0(ji,jj) - grav * zbeta * zws0(ji,jj) 402 ! turbulent heat flux averaged over depth of OSBL 403 zwthav(ji,jj) = 0.5 * zwth0(ji,jj) - ( 0.5*( zrad0(ji,jj) + zradh(ji,jj) ) - zradav(ji,jj) ) 404 ! turbulent salinity flux averaged over depth of the OBSL 405 zwsav(ji,jj) = 0.5 * zws0(ji,jj) 406 ! turbulent buoyancy flux averaged over the depth of the OBSBL 407 zwbav(ji,jj) = grav * zthermal * zwthav(ji,jj) - grav * zbeta * zwsav(ji,jj) 408 ! Surface upward velocity fluxes 409 zuw0(ji,jj) = -utau(ji,jj) * r1_rho0 * tmask(ji,jj,1) 410 zvw0(ji,jj) = -vtau(ji,jj) * r1_rho0 * tmask(ji,jj,1) 411 ! Friction velocity (zustar), at T-point : LMD94 eq. 2 412 zustar(ji,jj) = MAX( SQRT( SQRT( zuw0(ji,jj) * zuw0(ji,jj) + zvw0(ji,jj) * zvw0(ji,jj) ) ), 1.0e-8 ) 413 zcos_wind(ji,jj) = -zuw0(ji,jj) / ( zustar(ji,jj) * zustar(ji,jj) ) 414 zsin_wind(ji,jj) = -zvw0(ji,jj) / ( zustar(ji,jj) * zustar(ji,jj) ) 415 END_2D 417 416 ! Calculate Stokes drift in direction of wind (zustke) and Stokes penetration depth (dstokes) 418 417 SELECT CASE (nn_osm_wave) 419 418 ! Assume constant La#=0.3 420 419 CASE(0) 421 DO jj = 2, jpjm1 422 DO ji = 2, jpim1 423 zus_x = zcos_wind(ji,jj) * zustar(ji,jj) / 0.3**2 424 zus_y = zsin_wind(ji,jj) * zustar(ji,jj) / 0.3**2 425 zustke(ji,jj) = MAX ( SQRT( zus_x*zus_x + zus_y*zus_y), 1.0e-8 ) 426 ! dstokes(ji,jj) set to constant value rn_osm_dstokes from namelist in zdf_osm_init 427 END DO 428 END DO 420 DO_2D_00_00 421 zus_x = zcos_wind(ji,jj) * zustar(ji,jj) / 0.3**2 422 zus_y = zsin_wind(ji,jj) * zustar(ji,jj) / 0.3**2 423 zustke(ji,jj) = MAX ( SQRT( zus_x*zus_x + zus_y*zus_y), 1.0e-8 ) 424 ! dstokes(ji,jj) set to constant value rn_osm_dstokes from namelist in zdf_osm_init 425 END_2D 429 426 ! Assume Pierson-Moskovitz wind-wave spectrum 430 427 CASE(1) 431 DO jj = 2, jpjm1 432 DO ji = 2, jpim1 433 ! Use wind speed wndm included in sbc_oce module 434 zustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8 ) 435 dstokes(ji,jj) = MAX( 0.12 * wndm(ji,jj)**2 / grav, 5.e-1) 436 END DO 437 END DO 428 DO_2D_00_00 429 ! Use wind speed wndm included in sbc_oce module 430 zustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8 ) 431 dstokes(ji,jj) = MAX( 0.12 * wndm(ji,jj)**2 / grav, 5.e-1) 432 END_2D 438 433 ! Use ECMWF wave fields as output from SBCWAVE 439 434 CASE(2) 440 435 zfac = 2.0_wp * rpi / 16.0_wp 441 DO jj = 2, jpjm1 442 DO ji = 2, jpim1 443 ! The Langmur number from the ECMWF model appears to give La<0.3 for wind-driven seas. 444 ! The coefficient 0.8 gives La=0.3 in this situation. 445 ! It could represent the effects of the spread of wave directions 446 ! around the mean wind. The effect of this adjustment needs to be tested. 447 zabsstke = SQRT(ut0sd(ji,jj)**2 + vt0sd(ji,jj)**2) 448 zustke(ji,jj) = MAX (0.8 * ( zcos_wind(ji,jj) * ut0sd(ji,jj) + zsin_wind(ji,jj) * vt0sd(ji,jj) ), 1.0e-8) 449 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 ! 450 END DO 451 END DO 436 DO_2D_00_00 437 ! The Langmur number from the ECMWF model appears to give La<0.3 for wind-driven seas. 438 ! The coefficient 0.8 gives La=0.3 in this situation. 439 ! It could represent the effects of the spread of wave directions 440 ! around the mean wind. The effect of this adjustment needs to be tested. 441 zabsstke = SQRT(ut0sd(ji,jj)**2 + vt0sd(ji,jj)**2) 442 zustke(ji,jj) = MAX (0.8 * ( zcos_wind(ji,jj) * ut0sd(ji,jj) + zsin_wind(ji,jj) * vt0sd(ji,jj) ), 1.0e-8) 443 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 ! 444 END_2D 452 445 END SELECT 453 446 454 447 ! Langmuir velocity scale (zwstrl), La # (zla) 455 448 ! mixed scale (zvstr), convective velocity scale (zwstrc) 456 DO jj = 2, jpjm1 457 DO ji = 2, jpim1 458 ! Langmuir velocity scale (zwstrl), at T-point 459 zwstrl(ji,jj) = ( zustar(ji,jj) * zustar(ji,jj) * zustke(ji,jj) )**pthird 460 zla(ji,jj) = MAX(MIN(SQRT ( zustar(ji,jj) / ( zwstrl(ji,jj) + epsln ) )**3, 4.0), 0.2) 461 IF(zla(ji,jj) > 0.45) dstokes(ji,jj) = MIN(dstokes(ji,jj), 0.5_wp*hbl(ji,jj)) 462 ! Velocity scale that tends to zustar for large Langmuir numbers 463 zvstr(ji,jj) = ( zwstrl(ji,jj)**3 + & 464 & ( 1.0 - EXP( -0.5 * zla(ji,jj)**2 ) ) * zustar(ji,jj) * zustar(ji,jj) * zustar(ji,jj) )**pthird 465 466 ! limit maximum value of Langmuir number as approximate treatment for shear turbulence. 467 ! Note zustke and zwstrl are not amended. 468 ! 469 ! get convective velocity (zwstrc), stabilty scale (zhol) and logical conection flag lconv 470 IF ( zwbav(ji,jj) > 0.0) THEN 471 zwstrc(ji,jj) = ( 2.0 * zwbav(ji,jj) * 0.9 * hbl(ji,jj) )**pthird 472 zhol(ji,jj) = -0.9 * hbl(ji,jj) * 2.0 * zwbav(ji,jj) / (zvstr(ji,jj)**3 + epsln ) 473 lconv(ji,jj) = .TRUE. 474 ELSE 475 zhol(ji,jj) = -hbl(ji,jj) * 2.0 * zwbav(ji,jj)/ (zvstr(ji,jj)**3 + epsln ) 476 lconv(ji,jj) = .FALSE. 477 ENDIF 478 END DO 479 END DO 449 DO_2D_00_00 450 ! Langmuir velocity scale (zwstrl), at T-point 451 zwstrl(ji,jj) = ( zustar(ji,jj) * zustar(ji,jj) * zustke(ji,jj) )**pthird 452 zla(ji,jj) = MAX(MIN(SQRT ( zustar(ji,jj) / ( zwstrl(ji,jj) + epsln ) )**3, 4.0), 0.2) 453 IF(zla(ji,jj) > 0.45) dstokes(ji,jj) = MIN(dstokes(ji,jj), 0.5_wp*hbl(ji,jj)) 454 ! Velocity scale that tends to zustar for large Langmuir numbers 455 zvstr(ji,jj) = ( zwstrl(ji,jj)**3 + & 456 & ( 1.0 - EXP( -0.5 * zla(ji,jj)**2 ) ) * zustar(ji,jj) * zustar(ji,jj) * zustar(ji,jj) )**pthird 457 458 ! limit maximum value of Langmuir number as approximate treatment for shear turbulence. 459 ! Note zustke and zwstrl are not amended. 460 ! 461 ! get convective velocity (zwstrc), stabilty scale (zhol) and logical conection flag lconv 462 IF ( zwbav(ji,jj) > 0.0) THEN 463 zwstrc(ji,jj) = ( 2.0 * zwbav(ji,jj) * 0.9 * hbl(ji,jj) )**pthird 464 zhol(ji,jj) = -0.9 * hbl(ji,jj) * 2.0 * zwbav(ji,jj) / (zvstr(ji,jj)**3 + epsln ) 465 lconv(ji,jj) = .TRUE. 466 ELSE 467 zhol(ji,jj) = -hbl(ji,jj) * 2.0 * zwbav(ji,jj)/ (zvstr(ji,jj)**3 + epsln ) 468 lconv(ji,jj) = .FALSE. 469 ENDIF 470 END_2D 480 471 481 472 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 486 477 ! zdf_osm_zmld_horizontal_gradients need halo values for ibld, so must 487 478 ! previously exist for hbl also. 488 hbl(:,:) = MAX(hbl(:,:), gdepw _n(:,:,4) )479 hbl(:,:) = MAX(hbl(:,:), gdepw(:,:,4,Kmm) ) 489 480 ibld(:,:) = 4 490 DO jk = 5, jpkm1 491 DO jj = 1, jpj 492 DO ji = 1, jpi 493 IF ( hbl(ji,jj) >= gdepw_n(ji,jj,jk) ) THEN 494 ibld(ji,jj) = MIN(mbkt(ji,jj), jk) 495 ENDIF 496 END DO 497 END DO 498 END DO 499 500 DO jj = 2, jpjm1 501 DO ji = 2, jpim1 502 zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj)) 503 imld(ji,jj) = MAX(3,ibld(ji,jj) - MAX( INT( dh(ji,jj) / e3t_n(ji, jj, ibld(ji,jj) )) , 1 )) 504 zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 505 END DO 506 END DO 481 DO_3D_11_11( 5, jpkm1 ) 482 IF ( hbl(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN 483 ibld(ji,jj) = MIN(mbkt(ji,jj), jk) 484 ENDIF 485 END_3D 486 487 DO_2D_00_00 488 zhbl(ji,jj) = gdepw(ji,jj,ibld(ji,jj),Kmm) 489 imld(ji,jj) = MAX(3,ibld(ji,jj) - MAX( INT( dh(ji,jj) / e3t(ji, jj, ibld(ji,jj), Kmm )) , 1 )) 490 zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm) 491 END_2D 507 492 ! Averages over well-mixed and boundary layer 508 493 ! Alan: do we need zb_nl?, zb_ml? 509 CALL zdf_osm_vertical_average(ibld, zt_bl, zs_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl )510 CALL zdf_osm_vertical_average(imld, zt_ml, zs_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml )494 CALL zdf_osm_vertical_average(ibld, zt_bl, zs_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl, Kmm, Kbb ) 495 CALL zdf_osm_vertical_average(imld, zt_ml, zs_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml, Kmm, Kbb ) 511 496 ! External gradient 512 CALL zdf_osm_external_gradients( zdtdz_ext, zdsdz_ext, zdbdz_ext )497 CALL zdf_osm_external_gradients( zdtdz_ext, zdsdz_ext, zdbdz_ext, Kmm ) 513 498 514 499 515 500 IF ( ln_osm_mle ) THEN 516 CALL zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle )517 CALL zdf_osm_mle_parameters( hmle, zwb_fk, zvel_mle, zdiff_mle )501 CALL zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle, Kmm ) 502 CALL zdf_osm_mle_parameters( hmle, zwb_fk, zvel_mle, zdiff_mle, Kmm ) 518 503 ENDIF 519 504 … … 521 506 CALL zdf_osm_calculate_dhdt( zdhdt, zdhdt_2 ) 522 507 ! Calculate averages over depth of boundary layer 523 DO jj = 2, jpjm1 524 DO ji = 2, jpim1 525 zhbl_t(ji,jj) = hbl(ji,jj) + (zdhdt(ji,jj) - wn(ji,jj,ibld(ji,jj)))* rn_rdt ! certainly need wn here, so subtract it 526 ! adjustment to represent limiting by ocean bottom 527 zhbl_t(ji,jj) = MIN(zhbl_t(ji,jj), gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol)! ht_n(:,:)) 528 END DO 529 END DO 508 DO_2D_00_00 509 zhbl_t(ji,jj) = hbl(ji,jj) + (zdhdt(ji,jj) - ww(ji,jj,ibld(ji,jj)))* rn_Dt ! certainly need wn here, so subtract it 510 ! adjustment to represent limiting by ocean bottom 511 zhbl_t(ji,jj) = MIN(zhbl_t(ji,jj), gdepw(ji,jj, mbkt(ji,jj) + 1,Kmm) - depth_tol)! ht_n(:,:)) 512 END_2D 530 513 531 514 imld(:,:) = ibld(:,:) ! use imld to hold previous blayer index 532 515 ibld(:,:) = 4 533 516 534 DO jk = 4, jpkm1 535 DO jj = 2, jpjm1 536 DO ji = 2, jpim1 537 IF ( zhbl_t(ji,jj) >= gdepw_n(ji,jj,jk) ) THEN 538 ibld(ji,jj) = jk 539 ENDIF 540 END DO 541 END DO 542 END DO 517 518 DO_3D_00_00( 4, jpkm1 ) 519 IF ( zhbl_t(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN 520 ibld(ji,jj) = jk 521 ENDIF 522 END_3D 543 523 544 524 ! 545 525 ! Step through model levels taking account of buoyancy change to determine the effect on dhdt 546 526 ! 547 CALL zdf_osm_timestep_hbl( zdhdt, zdhdt_2 )527 CALL zdf_osm_timestep_hbl( zdhdt, zdhdt_2, Kmm ) 548 528 ! Alan: do we need zb_ml? 549 CALL zdf_osm_vertical_average( ibld, zt_bl, zs_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl )529 CALL zdf_osm_vertical_average( ibld, zt_bl, zs_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl, Kmm, Kbb ) 550 530 ! 551 531 ! 552 CALL zdf_osm_pycnocline_thickness( dh, zdh )532 CALL zdf_osm_pycnocline_thickness( dh, zdh, Kmm ) 553 533 dstokes(:,:) = MIN ( dstokes(:,:), hbl(:,:)/3. ) ! Limit delta for shallow boundary layers for calculating flux-gradient terms. 554 534 ! 555 535 ! Average over the depth of the mixed layer in the convective boundary layer 556 536 ! Alan: do we need zb_ml? 557 CALL zdf_osm_vertical_average( imld, zt_ml, zs_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml )537 CALL zdf_osm_vertical_average( imld, zt_ml, zs_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml, Kmm, Kbb ) 558 538 ! rotate mean currents and changes onto wind align co-ordinates 559 539 ! … … 564 544 zwth_ent = 0._wp 565 545 zws_ent = 0._wp 566 DO jj = 2, jpjm1 567 DO ji = 2, jpim1 568 IF ( ibld(ji,jj) < mbkt(ji,jj) ) THEN 569 IF ( lconv(ji,jj) ) THEN 570 zuw_bse(ji,jj) = -0.0075*((zvstr(ji,jj)**3+0.5*zwstrc(ji,jj)**3)**pthird*zdu_ml(ji,jj) + & 571 & 1.5*zustar(ji,jj)**2*(zhbl(ji,jj)-zhml(ji,jj)) )/ & 572 & ( zhml(ji,jj)*MIN(zla(ji,jj)**(8./3.),1.) + epsln) 573 zvw_bse(ji,jj) = 0.01*(-(zvstr(ji,jj)**3+0.5*zwstrc(ji,jj)**3)**pthird*zdv_ml(ji,jj)+ & 574 & 2.0*ff_t(ji,jj)*zustke(ji,jj)*dstokes(ji,jj)*zla(ji,jj)) 575 IF ( zdb_ml(ji,jj) > 0._wp ) THEN 576 zwth_ent(ji,jj) = zwb_ent(ji,jj) * zdt_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 577 zws_ent(ji,jj) = zwb_ent(ji,jj) * zds_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 578 ENDIF 579 ELSE 580 zwth_ent(ji,jj) = -2.0 * zwthav(ji,jj) * ( (1.0 - 0.8) - ( 1.0 - 0.8)**(3.0/2.0) ) 581 zws_ent(ji,jj) = -2.0 * zwsav(ji,jj) * ( (1.0 - 0.8 ) - ( 1.0 - 0.8 )**(3.0/2.0) ) 546 DO_2D_00_00 547 IF ( ibld(ji,jj) < mbkt(ji,jj) ) THEN 548 IF ( lconv(ji,jj) ) THEN 549 zuw_bse(ji,jj) = -0.0075*((zvstr(ji,jj)**3+0.5*zwstrc(ji,jj)**3)**pthird*zdu_ml(ji,jj) + & 550 & 1.5*zustar(ji,jj)**2*(zhbl(ji,jj)-zhml(ji,jj)) )/ & 551 & ( zhml(ji,jj)*MIN(zla(ji,jj)**(8./3.),1.) + epsln) 552 zvw_bse(ji,jj) = 0.01*(-(zvstr(ji,jj)**3+0.5*zwstrc(ji,jj)**3)**pthird*zdv_ml(ji,jj)+ & 553 & 2.0*ff_t(ji,jj)*zustke(ji,jj)*dstokes(ji,jj)*zla(ji,jj)) 554 IF ( zdb_ml(ji,jj) > 0._wp ) THEN 555 zwth_ent(ji,jj) = zwb_ent(ji,jj) * zdt_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 556 zws_ent(ji,jj) = zwb_ent(ji,jj) * zds_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 582 557 ENDIF 558 ELSE 559 zwth_ent(ji,jj) = -2.0 * zwthav(ji,jj) * ( (1.0 - 0.8) - ( 1.0 - 0.8)**(3.0/2.0) ) 560 zws_ent(ji,jj) = -2.0 * zwsav(ji,jj) * ( (1.0 - 0.8 ) - ( 1.0 - 0.8 )**(3.0/2.0) ) 583 561 ENDIF 584 END DO585 END DO562 ENDIF 563 END_2D 586 564 587 565 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 589 567 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 590 568 591 CALL zdf_osm_external_gradients( zdtdz_ext, zdsdz_ext, zdbdz_ext )592 CALL zdf_osm_pycnocline_scalar_profiles( zdtdz_pyc, zdsdz_pyc, zdbdz_pyc )593 CALL zdf_osm_pycnocline_shear_profiles( zdudz_pyc, zdvdz_pyc )569 CALL zdf_osm_external_gradients( zdtdz_ext, zdsdz_ext, zdbdz_ext, Kmm ) 570 CALL zdf_osm_pycnocline_scalar_profiles( zdtdz_pyc, zdsdz_pyc, zdbdz_pyc, Kmm ) 571 CALL zdf_osm_pycnocline_shear_profiles( zdudz_pyc, zdvdz_pyc, Kmm ) 594 572 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 595 573 ! Eddy viscosity/diffusivity and non-gradient terms in the flux-gradient relationship 596 574 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 597 575 598 DO jj = 2, jpjm1 599 DO ji = 2, jpim1 600 IF ( lconv(ji,jj) ) THEN 601 zdifml_sc(ji,jj) = zhml(ji,jj) * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 602 zvisml_sc(ji,jj) = zdifml_sc(ji,jj) 603 zdifpyc_sc(ji,jj) = 0.165 * ( zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3 )**pthird * zdh(ji,jj) 604 zvispyc_sc(ji,jj) = 0.142 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zdh(ji,jj) 605 zbeta_d_sc(ji,jj) = 1.0 - (0.165 / 0.8 * zdh(ji,jj) / zhbl(ji,jj) )**p2third 606 zbeta_v_sc(ji,jj) = 1.0 - 2.0 * (0.142 /0.375) * zdh(ji,jj) / ( zhml(ji,jj) + epsln ) 607 ELSE 608 zdifml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ) 609 zvisml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ) 610 END IF 611 END DO 612 END DO 576 DO_2D_00_00 577 IF ( lconv(ji,jj) ) THEN 578 zdifml_sc(ji,jj) = zhml(ji,jj) * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 579 zvisml_sc(ji,jj) = zdifml_sc(ji,jj) 580 zdifpyc_sc(ji,jj) = 0.165 * ( zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3 )**pthird * zdh(ji,jj) 581 zvispyc_sc(ji,jj) = 0.142 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zdh(ji,jj) 582 zbeta_d_sc(ji,jj) = 1.0 - (0.165 / 0.8 * zdh(ji,jj) / zhbl(ji,jj) )**p2third 583 zbeta_v_sc(ji,jj) = 1.0 - 2.0 * (0.142 /0.375) * zdh(ji,jj) / ( zhml(ji,jj) + epsln ) 584 ELSE 585 zdifml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ) 586 zvisml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ) 587 END IF 588 END_2D 613 589 ! 614 DO jj = 2, jpjm1 615 DO ji = 2, jpim1 616 IF ( lconv(ji,jj) ) THEN 617 DO jk = 2, imld(ji,jj) ! mixed layer diffusivity 618 zznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj) 590 DO_2D_00_00 591 IF ( lconv(ji,jj) ) THEN 592 DO jk = 2, imld(ji,jj) ! mixed layer diffusivity 593 zznd_ml = gdepw(ji,jj,jk,Kmm) / zhml(ji,jj) 594 ! 595 zdiffut(ji,jj,jk) = 0.8 * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_d_sc(ji,jj) * zznd_ml )**1.5 596 ! 597 zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_v_sc(ji,jj) * zznd_ml ) & 598 & * ( 1.0 - 0.5 * zznd_ml**2 ) 599 END DO 600 ! pycnocline - if present linear profile 601 IF ( zdh(ji,jj) > 0._wp ) THEN 602 zgamma_b = 6.0 603 DO jk = imld(ji,jj)+1 , ibld(ji,jj) 604 zznd_pyc = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj) 619 605 ! 620 zdiffut(ji,jj,jk) = 0.8 * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_d_sc(ji,jj) * zznd_ml )**1.5606 zdiffut(ji,jj,jk) = zdifpyc_sc(ji,jj) * EXP( zgamma_b * zznd_pyc ) 621 607 ! 622 zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_v_sc(ji,jj) * zznd_ml ) & 623 & * ( 1.0 - 0.5 * zznd_ml**2 ) 608 zviscos(ji,jj,jk) = zvispyc_sc(ji,jj) * EXP( zgamma_b * zznd_pyc ) 624 609 END DO 625 ! pycnocline - if present linear profile626 IF ( zdh(ji,jj) > 0._wp ) THEN627 zgamma_b = 6.0628 DO jk = imld(ji,jj)+1 , ibld(ji,jj)629 zznd_pyc = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj)630 !631 zdiffut(ji,jj,jk) = zdifpyc_sc(ji,jj) * EXP( zgamma_b * zznd_pyc )632 !633 zviscos(ji,jj,jk) = zvispyc_sc(ji,jj) * EXP( zgamma_b * zznd_pyc )634 END DO635 IF ( ibld_ext == 0 ) THEN636 zdiffut(ji,jj,ibld(ji,jj)) = 0._wp637 zviscos(ji,jj,ibld(ji,jj)) = 0._wp638 ELSE639 zdiffut(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj) * ( hbl(ji,jj) - gdepw_n(ji, jj, ibld(ji,jj)-1) )640 zviscos(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj) * ( hbl(ji,jj) - gdepw_n(ji, jj, ibld(ji,jj)-1) )641 ENDIF642 ENDIF643 ! Temporary fix to ensure zdiffut is +ve; won't be necessary with wn taken out644 zdiffut(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj) * e3t_n(ji,jj,ibld(ji,jj)), 1.e-6)645 zviscos(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj) * e3t_n(ji,jj,ibld(ji,jj)), 1.e-5)646 ! could be taken out, take account of entrainment represents as a diffusivity647 ! should remove w from here, represents entrainment648 ELSE649 ! stable conditions650 DO jk = 2, ibld(ji,jj)651 zznd_ml = gdepw_n(ji,jj,jk) / zhbl(ji,jj)652 zdiffut(ji,jj,jk) = 0.75 * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zznd_ml )**1.5653 zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * (1.0 - zznd_ml) * ( 1.0 - zznd_ml**2 )654 END DO655 656 610 IF ( ibld_ext == 0 ) THEN 657 611 zdiffut(ji,jj,ibld(ji,jj)) = 0._wp 658 612 zviscos(ji,jj,ibld(ji,jj)) = 0._wp 659 613 ELSE 660 zdiffut(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 0._wp) * e3w_n(ji, jj, ibld(ji,jj))661 zviscos(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 0._wp) * e3w_n(ji, jj, ibld(ji,jj))614 zdiffut(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj) * ( hbl(ji,jj) - gdepw(ji, jj, ibld(ji,jj)-1, Kmm) ) 615 zviscos(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj) * ( hbl(ji,jj) - gdepw(ji, jj, ibld(ji,jj)-1, Kmm) ) 662 616 ENDIF 663 ENDIF ! end if ( lconv ) 664 ! 665 END DO ! end of ji loop 666 END DO ! end of jj loop 617 ENDIF 618 ! Temporary fix to ensure zdiffut is +ve; won't be necessary with ww taken out 619 zdiffut(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj) * e3t(ji,jj,ibld(ji,jj), Kmm), 1.e-6) 620 zviscos(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj) * e3t(ji,jj,ibld(ji,jj), Kmm), 1.e-5) 621 ! could be taken out, take account of entrainment represents as a diffusivity 622 ! should remove w from here, represents entrainment 623 ELSE 624 ! stable conditions 625 DO jk = 2, ibld(ji,jj) 626 zznd_ml = gdepw(ji,jj,jk, Kmm) / zhbl(ji,jj) 627 zdiffut(ji,jj,jk) = 0.75 * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zznd_ml )**1.5 628 zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * (1.0 - zznd_ml) * ( 1.0 - zznd_ml**2 ) 629 END DO 630 631 IF ( ibld_ext == 0 ) THEN 632 zdiffut(ji,jj,ibld(ji,jj)) = 0._wp 633 zviscos(ji,jj,ibld(ji,jj)) = 0._wp 634 ELSE 635 zdiffut(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 0._wp) * e3w(ji, jj, ibld(ji,jj), Kmm) 636 zviscos(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 0._wp) * e3w(ji, jj, ibld(ji,jj), Kmm) 637 ENDIF 638 ENDIF ! end if ( lconv ) 639 ! 640 END_2D 667 641 668 642 ! … … 681 655 682 656 683 DO jj = 2, jpjm1 684 DO ji = 2, jpim1 685 IF ( lconv(ji,jj) ) THEN 686 DO jk = 2, imld(ji,jj) 687 zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj) 688 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.35 * EXP ( -zznd_d ) * ( 1.0 - EXP ( -2.0 * zznd_d ) ) * zsc_wth_1(ji,jj) 689 ! 690 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 1.35 * EXP ( -zznd_d ) * ( 1.0 - EXP ( -2.0 * zznd_d ) ) * zsc_ws_1(ji,jj) 691 END DO ! end jk loop 692 ELSE ! else for if (lconv) 657 DO_2D_00_00 658 IF ( lconv(ji,jj) ) THEN 659 DO jk = 2, imld(ji,jj) 660 zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 661 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.35 * EXP ( -zznd_d ) * ( 1.0 - EXP ( -2.0 * zznd_d ) ) * zsc_wth_1(ji,jj) 662 ! 663 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 1.35 * EXP ( -zznd_d ) * ( 1.0 - EXP ( -2.0 * zznd_d ) ) * zsc_ws_1(ji,jj) 664 END DO ! end jk loop 665 ELSE ! else for if (lconv) 693 666 ! Stable conditions 694 DO jk = 2, ibld(ji,jj) 695 zznd_d=gdepw_n(ji,jj,jk) / dstokes(ji,jj) 696 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.5 * EXP ( -0.9 * zznd_d ) & 697 & * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_wth_1(ji,jj) 698 ! 699 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 1.5 * EXP ( -0.9 * zznd_d ) & 700 & * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_ws_1(ji,jj) 701 END DO 702 ENDIF ! endif for check on lconv 703 704 END DO ! end of ji loop 705 END DO ! end of jj loop 667 DO jk = 2, ibld(ji,jj) 668 zznd_d=gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 669 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.5 * EXP ( -0.9 * zznd_d ) & 670 & * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_wth_1(ji,jj) 671 ! 672 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 1.5 * EXP ( -0.9 * zznd_d ) & 673 & * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_ws_1(ji,jj) 674 END DO 675 ENDIF ! endif for check on lconv 676 677 END_2D 706 678 707 679 ! 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) … … 718 690 IF ( iom_use("ghamv_00") ) CALL iom_put( "ghamv_00", wmask*ghamv ) 719 691 END IF 720 DO jj = 2, jpjm1 721 DO ji = 2, jpim1 722 IF ( lconv(ji,jj) ) THEN 723 DO jk = 2, imld(ji,jj) 724 zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj) 725 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + ( -0.05 * EXP ( -0.4 * zznd_d ) * zsc_uw_1(ji,jj) & 726 & + 0.00125 * EXP ( - zznd_d ) * zsc_uw_2(ji,jj) ) & 727 & * ( 1.0 - EXP ( -2.0 * zznd_d ) ) 692 DO_2D_00_00 693 IF ( lconv(ji,jj) ) THEN 694 DO jk = 2, imld(ji,jj) 695 zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 696 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + ( -0.05 * EXP ( -0.4 * zznd_d ) * zsc_uw_1(ji,jj) & 697 & + 0.00125 * EXP ( - zznd_d ) * zsc_uw_2(ji,jj) ) & 698 & * ( 1.0 - EXP ( -2.0 * zznd_d ) ) 728 699 ! 729 730 731 732 700 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.65 * 0.15 * EXP ( - zznd_d ) & 701 & * ( 1.0 - EXP ( -2.0 * zznd_d ) ) * zsc_vw_1(ji,jj) 702 END DO ! end jk loop 703 ELSE 733 704 ! Stable conditions 734 DO jk = 2, ibld(ji,jj) ! corrected to ibld 735 zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj) 736 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.75 * 1.3 * EXP ( -0.5 * zznd_d ) & 737 & * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_uw_1(ji,jj) 738 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + 0._wp 739 END DO ! end jk loop 740 ENDIF 741 END DO ! ji loop 742 END DO ! jj loo 705 DO jk = 2, ibld(ji,jj) ! corrected to ibld 706 zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 707 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.75 * 1.3 * EXP ( -0.5 * zznd_d ) & 708 & * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_uw_1(ji,jj) 709 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + 0._wp 710 END DO ! end jk loop 711 ENDIF 712 END_2D 743 713 744 714 ! Buoyancy term in flux-gradient relationship [note : includes ROI ratio (X0.3) and pressure (X0.5)] … … 752 722 ENDWHERE 753 723 754 DO jj = 2, jpjm1 755 DO ji = 2, jpim1 756 IF (lconv(ji,jj) ) THEN 757 DO jk = 2, imld(ji,jj) 758 zznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj) 759 ! calculate turbulent length scale 760 zl_c = 0.9 * ( 1.0 - EXP ( - 7.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) ) & 761 & * ( 1.0 - EXP ( -15.0 * ( 1.1 - zznd_ml ) ) ) 762 zl_l = 2.0 * ( 1.0 - EXP ( - 2.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) ) & 763 & * ( 1.0 - EXP ( - 5.0 * ( 1.0 - zznd_ml ) ) ) * ( 1.0 + dstokes(ji,jj) / zhml (ji,jj) ) 764 zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0 + EXP ( -3.0 * LOG10 ( - zhol(ji,jj) ) ) ) ** (3.0/2.0) 765 ! non-gradient buoyancy terms 766 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 ) 767 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * 0.5 * zsc_ws_1(ji,jj) * zl_eps * zhml(ji,jj) / ( 0.15 + zznd_ml ) 768 END DO 769 ELSE 770 DO jk = 2, ibld(ji,jj) 771 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zsc_wth_1(ji,jj) 772 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zsc_ws_1(ji,jj) 773 END DO 774 ENDIF 775 END DO ! ji loop 776 END DO ! jj loop 724 DO_2D_00_00 725 IF (lconv(ji,jj) ) THEN 726 DO jk = 2, imld(ji,jj) 727 zznd_ml = gdepw(ji,jj,jk,Kmm) / zhml(ji,jj) 728 ! calculate turbulent length scale 729 zl_c = 0.9 * ( 1.0 - EXP ( - 7.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) ) & 730 & * ( 1.0 - EXP ( -15.0 * ( 1.1 - zznd_ml ) ) ) 731 zl_l = 2.0 * ( 1.0 - EXP ( - 2.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) ) & 732 & * ( 1.0 - EXP ( - 5.0 * ( 1.0 - zznd_ml ) ) ) * ( 1.0 + dstokes(ji,jj) / zhml (ji,jj) ) 733 zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0 + EXP ( -3.0 * LOG10 ( - zhol(ji,jj) ) ) ) ** (3.0/2.0) 734 ! non-gradient buoyancy terms 735 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 ) 736 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * 0.5 * zsc_ws_1(ji,jj) * zl_eps * zhml(ji,jj) / ( 0.15 + zznd_ml ) 737 END DO 738 ELSE 739 DO jk = 2, ibld(ji,jj) 740 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zsc_wth_1(ji,jj) 741 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zsc_ws_1(ji,jj) 742 END DO 743 ENDIF 744 END_2D 777 745 778 746 WHERE ( lconv ) … … 785 753 ENDWHERE 786 754 787 DO jj = 2, jpjm1 788 DO ji = 2, jpim1 789 IF ( lconv(ji,jj) ) THEN 790 DO jk = 2 , imld(ji,jj) 791 zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj) 792 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.3 * 0.5 * ( zsc_uw_1(ji,jj) + 0.125 * EXP( -0.5 * zznd_d ) & 793 & * ( 1.0 - EXP( -0.5 * zznd_d ) ) & 794 & * zsc_uw_2(ji,jj) ) 795 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj) 796 END DO ! jk loop 797 ELSE 798 ! stable conditions 799 DO jk = 2, ibld(ji,jj) 800 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zsc_uw_1(ji,jj) 801 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj) 802 END DO 803 ENDIF 804 END DO ! ji loop 805 END DO ! jj loop 755 DO_2D_00_00 756 IF ( lconv(ji,jj) ) THEN 757 DO jk = 2 , imld(ji,jj) 758 zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 759 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.3 * 0.5 * ( zsc_uw_1(ji,jj) + 0.125 * EXP( -0.5 * zznd_d ) & 760 & * ( 1.0 - EXP( -0.5 * zznd_d ) ) & 761 & * zsc_uw_2(ji,jj) ) 762 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj) 763 END DO ! jk loop 764 ELSE 765 ! stable conditions 766 DO jk = 2, ibld(ji,jj) 767 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zsc_uw_1(ji,jj) 768 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj) 769 END DO 770 ENDIF 771 END_2D 806 772 807 773 IF(ln_dia_osm) THEN … … 819 785 ENDWHERE 820 786 821 DO jj = 2, jpjm1 822 DO ji = 2, jpim1 823 IF ( lconv(ji,jj) ) THEN 824 DO jk = 2, imld(ji,jj) 825 zznd_ml=gdepw_n(ji,jj,jk) / zhml(ji,jj) 826 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * zsc_wth_1(ji,jj) & 827 & * ( -2.0 + 2.75 * ( ( 1.0 + 0.6 * zznd_ml**4 ) & 828 & - EXP( - 6.0 * zznd_ml ) ) ) & 829 & * ( 1.0 - EXP( - 15.0 * ( 1.0 - zznd_ml ) ) ) 830 ! 831 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * zsc_ws_1(ji,jj) & 832 & * ( -2.0 + 2.75 * ( ( 1.0 + 0.6 * zznd_ml**4 ) & 833 & - EXP( - 6.0 * zznd_ml ) ) ) & 834 & * ( 1.0 - EXP ( -15.0 * ( 1.0 - zznd_ml ) ) ) 835 END DO 836 ELSE 837 DO jk = 2, ibld(ji,jj) 838 zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj) 839 znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 840 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + & 841 & 7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_wth_1(ji,jj) 842 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + & 843 & 7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_ws_1(ji,jj) 844 END DO 845 ENDIF 846 ENDDO ! ji loop 847 END DO ! jj loop 787 DO_2D_00_00 788 IF ( lconv(ji,jj) ) THEN 789 DO jk = 2, imld(ji,jj) 790 zznd_ml=gdepw(ji,jj,jk,Kmm) / zhml(ji,jj) 791 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * zsc_wth_1(ji,jj) & 792 & * ( -2.0 + 2.75 * ( ( 1.0 + 0.6 * zznd_ml**4 ) & 793 & - EXP( - 6.0 * zznd_ml ) ) ) & 794 & * ( 1.0 - EXP( - 15.0 * ( 1.0 - zznd_ml ) ) ) 795 ! 796 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * zsc_ws_1(ji,jj) & 797 & * ( -2.0 + 2.75 * ( ( 1.0 + 0.6 * zznd_ml**4 ) & 798 & - EXP( - 6.0 * zznd_ml ) ) ) & 799 & * ( 1.0 - EXP ( -15.0 * ( 1.0 - zznd_ml ) ) ) 800 END DO 801 ELSE 802 DO jk = 2, ibld(ji,jj) 803 zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 804 znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 805 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + & 806 & 7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_wth_1(ji,jj) 807 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + & 808 & 7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_ws_1(ji,jj) 809 END DO 810 ENDIF 811 END_2D 848 812 849 813 WHERE ( lconv ) … … 857 821 ENDWHERE 858 822 859 DO jj = 2, jpjm1 860 DO ji = 2, jpim1 861 IF ( lconv(ji,jj) ) THEN 862 DO jk = 2, imld(ji,jj) 863 zznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj) 864 zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj) 823 DO_2D_00_00 824 IF ( lconv(ji,jj) ) THEN 825 DO jk = 2, imld(ji,jj) 826 zznd_ml = gdepw(ji,jj,jk,Kmm) / zhml(ji,jj) 827 zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 828 ghamu(ji,jj,jk) = ghamu(ji,jj,jk)& 829 & + 0.3 * ( -2.0 + 2.5 * ( 1.0 + 0.1 * zznd_ml**4 ) - EXP ( -8.0 * zznd_ml ) ) * zsc_uw_1(ji,jj) 830 ! 831 ghamv(ji,jj,jk) = ghamv(ji,jj,jk)& 832 & + 0.3 * 0.1 * ( EXP( -zznd_d ) + EXP( -5.0 * ( 1.0 - zznd_ml ) ) ) * zsc_vw_1(ji,jj) 833 END DO 834 ELSE 835 DO jk = 2, ibld(ji,jj) 836 znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 837 zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 838 IF ( zznd_d <= 2.0 ) THEN 839 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.5 * 0.3 & 840 &* ( 2.25 - 3.0 * ( 1.0 - EXP( - 1.25 * zznd_d ) ) * ( 1.0 - EXP( -2.0 * zznd_d ) ) ) * zsc_uw_1(ji,jj) 841 ! 842 ELSE 865 843 ghamu(ji,jj,jk) = ghamu(ji,jj,jk)& 866 & + 0. 3 * ( -2.0 + 2.5 * ( 1.0 + 0.1 * zznd_ml**4 ) - EXP ( -8.0 * zznd_ml ) ) * zsc_uw_1(ji,jj)844 & + 0.5 * 0.3 * ( 1.0 - EXP( -5.0 * ( 1.0 - znd ) ) ) * zsc_uw_2(ji,jj) 867 845 ! 868 ghamv(ji,jj,jk) = ghamv(ji,jj,jk)& 869 & + 0.3 * 0.1 * ( EXP( -zznd_d ) + EXP( -5.0 * ( 1.0 - zznd_ml ) ) ) * zsc_vw_1(ji,jj) 870 END DO 871 ELSE 872 DO jk = 2, ibld(ji,jj) 873 znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 874 zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj) 875 IF ( zznd_d <= 2.0 ) THEN 876 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.5 * 0.3 & 877 &* ( 2.25 - 3.0 * ( 1.0 - EXP( - 1.25 * zznd_d ) ) * ( 1.0 - EXP( -2.0 * zznd_d ) ) ) * zsc_uw_1(ji,jj) 878 ! 879 ELSE 880 ghamu(ji,jj,jk) = ghamu(ji,jj,jk)& 881 & + 0.5 * 0.3 * ( 1.0 - EXP( -5.0 * ( 1.0 - znd ) ) ) * zsc_uw_2(ji,jj) 882 ! 883 ENDIF 884 885 ghamv(ji,jj,jk) = ghamv(ji,jj,jk)& 886 & + 0.3 * 0.15 * SIN( 3.14159 * ( 0.65 * zznd_d ) ) * EXP( -0.25 * zznd_d**2 ) * zsc_vw_1(ji,jj) 887 ghamv(ji,jj,jk) = ghamv(ji,jj,jk)& 888 & + 0.3 * 0.15 * EXP( -5.0 * ( 1.0 - znd ) ) * ( 1.0 - EXP( -20.0 * ( 1.0 - znd ) ) ) * zsc_vw_2(ji,jj) 889 END DO 890 ENDIF 891 END DO 892 END DO 846 ENDIF 847 848 ghamv(ji,jj,jk) = ghamv(ji,jj,jk)& 849 & + 0.3 * 0.15 * SIN( 3.14159 * ( 0.65 * zznd_d ) ) * EXP( -0.25 * zznd_d**2 ) * zsc_vw_1(ji,jj) 850 ghamv(ji,jj,jk) = ghamv(ji,jj,jk)& 851 & + 0.3 * 0.15 * EXP( -5.0 * ( 1.0 - znd ) ) * ( 1.0 - EXP( -20.0 * ( 1.0 - znd ) ) ) * zsc_vw_2(ji,jj) 852 END DO 853 ENDIF 854 END_2D 893 855 894 856 IF(ln_dia_osm) THEN … … 903 865 ! Make surface forced velocity non-gradient terms go to zero at the base of the mixed layer. 904 866 905 DO jj = 2, jpjm1 906 DO ji = 2, jpim1 907 IF ( lconv(ji,jj) ) THEN 908 DO jk = 2, ibld(ji,jj) 909 znd = ( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zhml(ji,jj) !ALMG to think about 910 IF ( znd >= 0.0 ) THEN 911 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -30.0 * znd**2 ) ) 912 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * ( 1.0 - EXP( -30.0 * znd**2 ) ) 913 ELSE 914 ghamu(ji,jj,jk) = 0._wp 915 ghamv(ji,jj,jk) = 0._wp 916 ENDIF 917 END DO 918 ELSE 919 DO jk = 2, ibld(ji,jj) 920 znd = ( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zhml(ji,jj) !ALMG to think about 921 IF ( znd >= 0.0 ) THEN 922 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) ) 923 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) ) 924 ELSE 925 ghamu(ji,jj,jk) = 0._wp 926 ghamv(ji,jj,jk) = 0._wp 927 ENDIF 928 END DO 929 ENDIF 930 END DO 931 END DO 867 DO_2D_00_00 868 IF ( lconv(ji,jj) ) THEN 869 DO jk = 2, ibld(ji,jj) 870 znd = ( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zhml(ji,jj) !ALMG to think about 871 IF ( znd >= 0.0 ) THEN 872 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -30.0 * znd**2 ) ) 873 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * ( 1.0 - EXP( -30.0 * znd**2 ) ) 874 ELSE 875 ghamu(ji,jj,jk) = 0._wp 876 ghamv(ji,jj,jk) = 0._wp 877 ENDIF 878 END DO 879 ELSE 880 DO jk = 2, ibld(ji,jj) 881 znd = ( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zhml(ji,jj) !ALMG to think about 882 IF ( znd >= 0.0 ) THEN 883 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) ) 884 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) ) 885 ELSE 886 ghamu(ji,jj,jk) = 0._wp 887 ghamv(ji,jj,jk) = 0._wp 888 ENDIF 889 END DO 890 ENDIF 891 END_2D 932 892 933 893 IF(ln_dia_osm) THEN … … 938 898 ! Temporary fix to avoid instabilities when zdb_bl becomes very very small 939 899 zsc_uw_1 = 0._wp ! 50.0 * zla**(8.0/3.0) * zustar**2 * zhbl / ( zdb_bl + epsln ) 940 DO jj = 2, jpjm1 941 DO ji = 2, jpim1 942 IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN 943 DO jk= 2, ibld(ji,jj) 944 znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 945 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk) 946 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk) 947 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc(ji,jj,jk) 948 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) 949 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc(ji,jj,jk) 950 END DO 951 END IF 952 END DO 953 END DO 900 DO_2D_00_00 901 IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN 902 DO jk= 2, ibld(ji,jj) 903 znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 904 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk) 905 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk) 906 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc(ji,jj,jk) 907 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) 908 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc(ji,jj,jk) 909 END DO 910 END IF 911 END_2D 954 912 955 913 ! Entrainment contribution. 956 914 957 DO jj=2, jpjm1 958 DO ji = 2, jpim1 959 IF ( lconv(ji,jj) .AND. ibld(ji,jj) + ibld_ext < mbkt(ji,jj)) THEN 960 DO jk = 1, imld(ji,jj) - 1 961 znd=gdepw_n(ji,jj,jk) / zhml(ji,jj) 962 ! ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * znd 963 ! ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * znd 964 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * znd 965 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * znd 966 END DO 967 DO jk = imld(ji,jj), ibld(ji,jj) 968 znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 969 ! ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * ( 1.0 + znd ) 970 ! ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * ( 1.0 + znd ) 971 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * ( 1.0 + znd ) 972 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * ( 1.0 + znd ) 973 END DO 974 ENDIF 975 976 ghamt(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 977 ghams(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 978 ghamu(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 979 ghamv(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 980 END DO ! ji loop 981 END DO ! jj loop 915 DO_2D_00_00 916 IF ( lconv(ji,jj) .AND. ibld(ji,jj) + ibld_ext < mbkt(ji,jj)) THEN 917 DO jk = 1, imld(ji,jj) - 1 918 znd=gdepw(ji,jj,jk,kmm) / zhml(ji,jj) 919 ! ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * znd 920 ! ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * znd 921 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * znd 922 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * znd 923 END DO 924 DO jk = imld(ji,jj), ibld(ji,jj) 925 znd = -( gdepw(ji,jj,jk,kmm) - zhml(ji,jj) ) / zdh(ji,jj) 926 ! ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * ( 1.0 + znd ) 927 ! ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * ( 1.0 + znd ) 928 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * ( 1.0 + znd ) 929 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * ( 1.0 + znd ) 930 END DO 931 ENDIF 932 933 ghamt(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 934 ghams(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 935 ghamu(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 936 ghamv(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 937 END_2D 982 938 983 939 IF(ln_dia_osm) THEN … … 1001 957 ! rotate non-gradient velocity terms back to model reference frame 1002 958 1003 DO jj = 2, jpjm1 1004 DO ji = 2, jpim1 1005 DO jk = 2, ibld(ji,jj) 1006 ztemp = ghamu(ji,jj,jk) 1007 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * zcos_wind(ji,jj) - ghamv(ji,jj,jk) * zsin_wind(ji,jj) 1008 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * zcos_wind(ji,jj) + ztemp * zsin_wind(ji,jj) 1009 END DO 959 DO_2D_00_00 960 DO jk = 2, ibld(ji,jj) 961 ztemp = ghamu(ji,jj,jk) 962 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * zcos_wind(ji,jj) - ghamv(ji,jj,jk) * zsin_wind(ji,jj) 963 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * zcos_wind(ji,jj) + ztemp * zsin_wind(ji,jj) 1010 964 END DO 1011 END DO965 END_2D 1012 966 1013 967 IF(ln_dia_osm) THEN … … 1019 973 ! KPP-style Ri# mixing 1020 974 IF( ln_kpprimix) THEN 1021 DO jk = 2, jpkm1 !* Shear production at uw- and vw-points (energy conserving form) 1022 DO jj = 1, jpjm1 1023 DO ji = 1, jpim1 ! vector opt. 1024 z3du(ji,jj,jk) = 0.5 * ( un(ji,jj,jk-1) - un(ji ,jj,jk) ) & 1025 & * ( ub(ji,jj,jk-1) - ub(ji ,jj,jk) ) * wumask(ji,jj,jk) & 1026 & / ( e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) ) 1027 z3dv(ji,jj,jk) = 0.5 * ( vn(ji,jj,jk-1) - vn(ji,jj ,jk) ) & 1028 & * ( vb(ji,jj,jk-1) - vb(ji,jj ,jk) ) * wvmask(ji,jj,jk) & 1029 & / ( e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) ) 1030 END DO 975 DO_3D_10_10( 2, jpkm1 ) 976 z3du(ji,jj,jk) = 0.5 * ( uu(ji,jj,jk-1,Kmm) - uu(ji ,jj,jk,Kmm) ) & 977 & * ( uu(ji,jj,jk-1,Kbb) - uu(ji ,jj,jk,Kbb) ) * wumask(ji,jj,jk) & 978 & / ( e3uw(ji,jj,jk,Kmm) * e3uw(ji,jj,jk,Kbb) ) 979 z3dv(ji,jj,jk) = 0.5 * ( vv(ji,jj,jk-1,Kmm) - vv(ji,jj ,jk,Kmm) ) & 980 & * ( vv(ji,jj,jk-1,Kbb) - vv(ji,jj ,jk,Kbb) ) * wvmask(ji,jj,jk) & 981 & / ( e3vw(ji,jj,jk,Kmm) * e3vw(ji,jj,jk,Kbb) ) 982 END_3D 983 ! 984 DO_3D_00_00( 2, jpkm1 ) 985 ! ! shear prod. at w-point weightened by mask 986 zesh2 = ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) & 987 & + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) ) 988 ! ! local Richardson number 989 zri = MAX( rn2b(ji,jj,jk), 0._wp ) / MAX(zesh2, epsln) 990 zfri = MIN( zri / rn_riinfty , 1.0_wp ) 991 zfri = ( 1.0_wp - zfri * zfri ) 992 zrimix(ji,jj,jk) = zfri * zfri * zfri * wmask(ji, jj, jk) 993 END_3D 994 995 DO_2D_00_00 996 DO jk = ibld(ji,jj) + 1, jpkm1 997 zdiffut(ji,jj,jk) = zrimix(ji,jj,jk)*rn_difri 998 zviscos(ji,jj,jk) = zrimix(ji,jj,jk)*rn_difri 1031 999 END DO 1032 END DO 1033 ! 1034 DO jk = 2, jpkm1 1035 DO jj = 2, jpjm1 1036 DO ji = 2, jpim1 ! vector opt. 1037 ! ! shear prod. at w-point weightened by mask 1038 zesh2 = ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) & 1039 & + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) ) 1040 ! ! local Richardson number 1041 zri = MAX( rn2b(ji,jj,jk), 0._wp ) / MAX(zesh2, epsln) 1042 zfri = MIN( zri / rn_riinfty , 1.0_wp ) 1043 zfri = ( 1.0_wp - zfri * zfri ) 1044 zrimix(ji,jj,jk) = zfri * zfri * zfri * wmask(ji, jj, jk) 1045 END DO 1046 END DO 1047 END DO 1048 1049 DO jj = 2, jpjm1 1050 DO ji = 2, jpim1 1051 DO jk = ibld(ji,jj) + 1, jpkm1 1052 zdiffut(ji,jj,jk) = zrimix(ji,jj,jk)*rn_difri 1053 zviscos(ji,jj,jk) = zrimix(ji,jj,jk)*rn_difri 1054 END DO 1055 END DO 1056 END DO 1000 END_2D 1057 1001 1058 1002 END IF ! ln_kpprimix = .true. … … 1060 1004 ! KPP-style set diffusivity large if unstable below BL 1061 1005 IF( ln_convmix) THEN 1062 DO jj = 2, jpjm1 1063 DO ji = 2, jpim1 1064 DO jk = ibld(ji,jj) + 1, jpkm1 1065 IF( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) zdiffut(ji,jj,jk) = rn_difconv 1066 END DO 1006 DO_2D_00_00 1007 DO jk = ibld(ji,jj) + 1, jpkm1 1008 IF( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) zdiffut(ji,jj,jk) = rn_difconv 1067 1009 END DO 1068 END DO1010 END_2D 1069 1011 END IF ! ln_convmix = .true. 1070 1012 … … 1072 1014 1073 1015 IF ( ln_osm_mle ) THEN ! set up diffusivity and non-gradient mixing 1074 DO jj = 2 , jpjm1 1075 DO ji = 2, jpim1 1076 IF ( lconv(ji,jj) .AND. mld_prof(ji,jj) - ibld(ji,jj) > 1 ) THEN ! MLE mmixing extends below the OSBL. 1077 ! Calculate MLE flux profiles 1078 ! DO jk = 1, mld_prof(ji,jj) 1079 ! znd = - gdepw_n(ji,jj,jk) / MAX(zhmle(ji,jj),epsln) 1080 ! ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + & 1081 ! & zwt_fk(ji,jj) * ( 1.0 - ( 2.0 * znd + 1.0 )**2 ) * ( 1.0 + r5_21 * ( 2.0 * znd + 1.0 )**2 ) 1082 ! ghams(ji,jj,jk) = ghams(ji,jj,jk) + & 1083 ! & zws_fk(ji,jj) * ( 1.0 - ( 2.0 * znd + 1.0 )**2 ) * ( 1.0 + r5_21 * ( 2.0 * znd + 1.0 )**2 ) 1084 ! END DO 1085 ! Calculate MLE flux contribution from surface fluxes 1086 DO jk = 1, ibld(ji,jj) 1087 znd = gdepw_n(ji,jj,jk) / MAX(zhbl(ji,jj),epsln) 1088 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) - zwth0(ji,jj) * ( 1.0 - znd ) 1089 ghams(ji,jj,jk) = ghams(ji,jj,jk) - zws0(ji,jj) * ( 1.0 - znd ) 1090 END DO 1091 DO jk = 1, mld_prof(ji,jj) 1092 znd = gdepw_n(ji,jj,jk) / MAX(zhmle(ji,jj),epsln) 1093 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth0(ji,jj) * ( 1.0 - znd ) 1094 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws0(ji,jj) * ( 1.0 -znd ) 1095 END DO 1096 ! Viscosity for MLEs 1097 DO jk = ibld(ji,jj), mld_prof(ji,jj) 1098 zdiffut(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), zdiff_mle(ji,jj) ) 1099 END DO 1100 ! Iterate to find approx vertical index for depth 1.1*zhmle(ji,jj) 1101 jl = MIN(mld_prof(ji,jj) + 2, mbkt(ji,jj)) 1102 jl = MIN( MAX(INT( 0.1 * zhmle(ji,jj) / e3t_n(ji,jj,jl)), 2 ) + mld_prof(ji,jj), mbkt(ji,jj)) 1103 DO jk = mld_prof(ji,jj), jl 1104 zdiffut(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), zdiff_mle(ji,jj) * & 1105 & ( gdepw_n(ji,jj,jk) - gdepw_n(ji,jj,jl) ) / & 1106 & ( gdepw_n(ji,jj,mld_prof(ji,jj)) - gdepw_n(ji,jj,jl) - epsln)) 1107 END DO 1108 ENDIF 1109 END DO 1110 END DO 1016 DO_2D_00_00 1017 IF ( lconv(ji,jj) .AND. mld_prof(ji,jj) - ibld(ji,jj) > 1 ) THEN ! MLE mmixing extends below the OSBL. 1018 ! Calculate MLE flux profiles 1019 ! DO jk = 1, mld_prof(ji,jj) 1020 ! znd = - gdepw(ji,jj,jk, Kmm ) / MAX(zhmle(ji,jj),epsln) 1021 ! ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + & 1022 ! & zwt_fk(ji,jj) * ( 1.0 - ( 2.0 * znd + 1.0 )**2 ) * ( 1.0 + r5_21 * ( 2.0 * znd + 1.0 )**2 ) 1023 ! ghams(ji,jj,jk) = ghams(ji,jj,jk) + & 1024 ! & zws_fk(ji,jj) * ( 1.0 - ( 2.0 * znd + 1.0 )**2 ) * ( 1.0 + r5_21 * ( 2.0 * znd + 1.0 )**2 ) 1025 ! END DO 1026 ! Calculate MLE flux contribution from surface fluxes 1027 DO jk = 1, ibld(ji,jj) 1028 znd = gdepw(ji,jj,jk, Kmm ) / MAX(zhbl(ji,jj),epsln) 1029 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) - zwth0(ji,jj) * ( 1.0 - znd ) 1030 ghams(ji,jj,jk) = ghams(ji,jj,jk) - zws0(ji,jj) * ( 1.0 - znd ) 1031 END DO 1032 DO jk = 1, mld_prof(ji,jj) 1033 znd = gdepw(ji,jj,jk, Kmm ) / MAX(zhmle(ji,jj),epsln) 1034 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth0(ji,jj) * ( 1.0 - znd ) 1035 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws0(ji,jj) * ( 1.0 -znd ) 1036 END DO 1037 ! Viscosity for MLEs 1038 DO jk = ibld(ji,jj), mld_prof(ji,jj) 1039 zdiffut(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), zdiff_mle(ji,jj) ) 1040 END DO 1041 ! Iterate to find approx vertical index for depth 1.1*zhmle(ji,jj) 1042 jl = MIN(mld_prof(ji,jj) + 2, mbkt(ji,jj)) 1043 jl = MIN( MAX(INT( 0.1 * zhmle(ji,jj) / e3t(ji,jj,jl, Kmm ) ), 2 ) + mld_prof(ji,jj), mbkt(ji,jj)) 1044 DO jk = mld_prof(ji,jj), jl 1045 zdiffut(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), zdiff_mle(ji,jj) * & 1046 & ( gdepw(ji,jj,jk, Kmm ) - gdepw(ji,jj,jl, Kmm ) ) / & 1047 & ( gdepw(ji,jj,mld_prof(ji,jj), Kmm ) - gdepw(ji,jj,jl, Kmm ) - epsln)) 1048 END DO 1049 ENDIF 1050 END_2D 1111 1051 ENDIF 1112 1052 … … 1123 1063 ! GN 25/8: need to change tmask --> wmask 1124 1064 1125 DO jk = 2, jpkm1 1126 DO jj = 2, jpjm1 1127 DO ji = 2, jpim1 1128 p_avt(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), avtb(jk) ) * tmask(ji,jj,jk) 1129 p_avm(ji,jj,jk) = MAX( zviscos(ji,jj,jk), avmb(jk) ) * tmask(ji,jj,jk) 1130 END DO 1131 END DO 1132 END DO 1065 DO_3D_00_00( 2, jpkm1 ) 1066 p_avt(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), avtb(jk) ) * tmask(ji,jj,jk) 1067 p_avm(ji,jj,jk) = MAX( zviscos(ji,jj,jk), avmb(jk) ) * tmask(ji,jj,jk) 1068 END_3D 1133 1069 ! Lateral boundary conditions on ghamu and ghamv, currently on W-grid (sign unchanged), needed to caclulate gham[uv] on u and v grids 1134 1070 CALL lbc_lnk_multi( 'zdfosm', p_avt, 'W', 1. , p_avm, 'W', 1., & 1135 1071 & ghamu, 'W', 1. , ghamv, 'W', 1. ) 1136 DO jk = 2, jpkm1 1137 DO jj = 2, jpjm1 1138 DO ji = 2, jpim1 1139 ghamu(ji,jj,jk) = ( ghamu(ji,jj,jk) + ghamu(ji+1,jj,jk) ) & 1140 & / MAX( 1., tmask(ji,jj,jk) + tmask (ji + 1,jj,jk) ) * umask(ji,jj,jk) 1141 1142 ghamv(ji,jj,jk) = ( ghamv(ji,jj,jk) + ghamv(ji,jj+1,jk) ) & 1143 & / MAX( 1., tmask(ji,jj,jk) + tmask (ji,jj+1,jk) ) * vmask(ji,jj,jk) 1144 1145 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) * tmask(ji,jj,jk) 1146 ghams(ji,jj,jk) = ghams(ji,jj,jk) * tmask(ji,jj,jk) 1147 END DO 1148 END DO 1149 END DO 1072 DO_3D_00_00( 2, jpkm1 ) 1073 ghamu(ji,jj,jk) = ( ghamu(ji,jj,jk) + ghamu(ji+1,jj,jk) ) & 1074 & / MAX( 1., tmask(ji,jj,jk) + tmask (ji + 1,jj,jk) ) * umask(ji,jj,jk) 1075 1076 ghamv(ji,jj,jk) = ( ghamv(ji,jj,jk) + ghamv(ji,jj+1,jk) ) & 1077 & / MAX( 1., tmask(ji,jj,jk) + tmask (ji,jj+1,jk) ) * vmask(ji,jj,jk) 1078 1079 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) * tmask(ji,jj,jk) 1080 ghams(ji,jj,jk) = ghams(ji,jj,jk) * tmask(ji,jj,jk) 1081 END_3D 1150 1082 ! Lateral boundary conditions on final outputs for hbl, on T-grid (sign unchanged) 1151 1083 CALL lbc_lnk_multi( 'zdfosm', hbl, 'T', 1., dh, 'T', 1., hmle, 'T', 1. ) … … 1161 1093 IF ( iom_use("us_x") ) CALL iom_put( "us_x", tmask(:,:,1)*zustke*zcos_wind ) ! x surface Stokes drift 1162 1094 IF ( iom_use("us_y") ) CALL iom_put( "us_y", tmask(:,:,1)*zustke*zsin_wind ) ! y surface Stokes drift 1163 IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*r au0*tmask(:,:,1)*zustar**2*zustke )1095 IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rho0*tmask(:,:,1)*zustar**2*zustke ) 1164 1096 ! Stokes drift read in from sbcwave (=2). 1165 1097 CASE(2) … … 1171 1103 IF ( iom_use("hsw_NP") ) CALL iom_put( "hsw_NP", (0.22/grav)*wndm**2*tmask(:,:,1) ) ! significant wave height from NP spectrum 1172 1104 IF ( iom_use("wndm") ) CALL iom_put( "wndm", wndm*tmask(:,:,1) ) ! U_10 1173 IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*r au0*tmask(:,:,1)*zustar**2* &1105 IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rho0*tmask(:,:,1)*zustar**2* & 1174 1106 & SQRT(ut0sd**2 + vt0sd**2 ) ) 1175 1107 END SELECT … … 1196 1128 IF ( iom_use("zvstr") ) CALL iom_put( "zvstr", tmask(:,:,1)*zvstr ) ! mixed velocity scale 1197 1129 IF ( iom_use("zla") ) CALL iom_put( "zla", tmask(:,:,1)*zla ) ! langmuir # 1198 IF ( iom_use("wind_power") ) CALL iom_put( "wind_power", 1000.*r au0*tmask(:,:,1)*zustar**3 ) ! BL depth internal to zdf_osm routine1199 IF ( iom_use("wind_wave_power") ) CALL iom_put( "wind_wave_power", 1000.*r au0*tmask(:,:,1)*zustar**2*zustke )1130 IF ( iom_use("wind_power") ) CALL iom_put( "wind_power", 1000.*rho0*tmask(:,:,1)*zustar**3 ) ! BL depth internal to zdf_osm routine 1131 IF ( iom_use("wind_wave_power") ) CALL iom_put( "wind_wave_power", 1000.*rho0*tmask(:,:,1)*zustar**2*zustke ) 1200 1132 IF ( iom_use("zhbl") ) CALL iom_put( "zhbl", tmask(:,:,1)*zhbl ) ! BL depth internal to zdf_osm routine 1201 1133 IF ( iom_use("zhml") ) CALL iom_put( "zhml", tmask(:,:,1)*zhml ) ! ML depth internal to zdf_osm routine … … 1229 1161 1230 1162 ! Alan: do we need zb? 1231 SUBROUTINE zdf_osm_vertical_average( jnlev_av, zt, zs, zu, zv, zdt, zds, zdb, zdu, zdv )1163 SUBROUTINE zdf_osm_vertical_average( jnlev_av, zt, zs, zu, zv, zdt, zds, zdb, zdu, zdv, Kmm, Kbb ) 1232 1164 !!--------------------------------------------------------------------- 1233 1165 !! *** ROUTINE zdf_vertical_average *** … … 1240 1172 !!---------------------------------------------------------------------- 1241 1173 1242 INTEGER, DIMENSION(jpi,jpj) :: jnlev_av ! Number of levels to average over. 1174 INTEGER, INTENT(in) :: Kmm, Kbb ! time-level indices 1175 INTEGER, DIMENSION(jpi,jpj) :: jnlev_av ! Number of levels to average over. 1243 1176 1244 1177 ! Alan: do we need zb? … … 1256 1189 zu = 0._wp 1257 1190 zv = 0._wp 1258 DO jj = 2, jpjm1 ! Vertical slab 1259 DO ji = 2, jpim1 1260 zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 1261 zbeta = rab_n(ji,jj,1,jp_sal) 1262 ! average over depth of boundary layer 1263 zthick = epsln 1264 DO jk = 2, jnlev_av(ji,jj) 1265 zthick = zthick + e3t_n(ji,jj,jk) 1266 zt(ji,jj) = zt(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem) 1267 zs(ji,jj) = zs(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) 1268 zu(ji,jj) = zu(ji,jj) + e3t_n(ji,jj,jk) & 1269 & * ( ub(ji,jj,jk) + ub(ji - 1,jj,jk) ) & 1270 & / MAX( 1. , umask(ji,jj,jk) + umask(ji - 1,jj,jk) ) 1271 zv(ji,jj) = zv(ji,jj) + e3t_n(ji,jj,jk) & 1272 & * ( vb(ji,jj,jk) + vb(ji,jj - 1,jk) ) & 1273 & / MAX( 1. , vmask(ji,jj,jk) + vmask(ji,jj - 1,jk) ) 1274 END DO 1275 zt(ji,jj) = zt(ji,jj) / zthick 1276 zs(ji,jj) = zs(ji,jj) / zthick 1277 zu(ji,jj) = zu(ji,jj) / zthick 1278 zv(ji,jj) = zv(ji,jj) / zthick 1279 ! Alan: do we need zb? 1280 zdt(ji,jj) = zt(ji,jj) - tsn(ji,jj,ibld(ji,jj)+ibld_ext,jp_tem) 1281 zds(ji,jj) = zs(ji,jj) - tsn(ji,jj,ibld(ji,jj)+ibld_ext,jp_sal) 1282 zdu(ji,jj) = zu(ji,jj) - ( ub(ji,jj,ibld(ji,jj)+ibld_ext) + ub(ji-1,jj,ibld(ji,jj)+ibld_ext ) ) & 1283 & / MAX(1. , umask(ji,jj,ibld(ji,jj)+ibld_ext ) + umask(ji-1,jj,ibld(ji,jj)+ibld_ext ) ) 1284 zdv(ji,jj) = zv(ji,jj) - ( vb(ji,jj,ibld(ji,jj)+ibld_ext) + vb(ji,jj-1,ibld(ji,jj)+ibld_ext ) ) & 1285 & / MAX(1. , vmask(ji,jj,ibld(ji,jj)+ibld_ext ) + vmask(ji,jj-1,ibld(ji,jj)+ibld_ext ) ) 1286 zdb(ji,jj) = grav * zthermal * zdt(ji,jj) - grav * zbeta * zds(ji,jj) 1191 DO_2D_00_00 1192 zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 1193 zbeta = rab_n(ji,jj,1,jp_sal) 1194 ! average over depth of boundary layer 1195 zthick = epsln 1196 DO jk = 2, jnlev_av(ji,jj) 1197 zthick = zthick + e3t(ji,jj,jk, Kmm ) 1198 zt(ji,jj) = zt(ji,jj) + e3t(ji,jj,jk, Kmm ) * ts(ji,jj,jk,jp_tem, Kmm ) 1199 zs(ji,jj) = zs(ji,jj) + e3t(ji,jj,jk, Kmm ) * ts(ji,jj,jk,jp_sal, Kmm ) 1200 zu(ji,jj) = zu(ji,jj) + e3t(ji,jj,jk, Kmm ) & 1201 & * ( uu(ji,jj,jk, Kbb ) + uu(ji - 1,jj,jk, Kbb ) ) & 1202 & / MAX( 1. , umask(ji,jj,jk) + umask(ji - 1,jj,jk) ) 1203 zv(ji,jj) = zv(ji,jj) + e3t(ji,jj,jk, Kmm ) & 1204 & * ( vv(ji,jj,jk, Kbb ) + vv(ji,jj - 1,jk, Kbb ) ) & 1205 & / MAX( 1. , vmask(ji,jj,jk) + vmask(ji,jj - 1,jk) ) 1287 1206 END DO 1288 END DO 1207 zt(ji,jj) = zt(ji,jj) / zthick 1208 zs(ji,jj) = zs(ji,jj) / zthick 1209 zu(ji,jj) = zu(ji,jj) / zthick 1210 zv(ji,jj) = zv(ji,jj) / zthick 1211 ! Alan: do we need zb? 1212 zdt(ji,jj) = zt(ji,jj) - ts(ji,jj,ibld(ji,jj)+ibld_ext,jp_tem, Kmm ) 1213 zds(ji,jj) = zs(ji,jj) - ts(ji,jj,ibld(ji,jj)+ibld_ext,jp_sal, Kmm ) 1214 zdu(ji,jj) = zu(ji,jj) - ( uu(ji,jj,ibld(ji,jj)+ibld_ext, Kbb ) + uu(ji-1,jj,ibld(ji,jj)+ibld_ext, Kbb ) ) & 1215 & / MAX(1. , umask(ji,jj,ibld(ji,jj)+ibld_ext ) + umask(ji-1,jj,ibld(ji,jj)+ibld_ext ) ) 1216 zdv(ji,jj) = zv(ji,jj) - ( vv(ji,jj,ibld(ji,jj)+ibld_ext, Kbb ) + vv(ji,jj-1,ibld(ji,jj)+ibld_ext, Kbb ) ) & 1217 & / MAX(1. , vmask(ji,jj,ibld(ji,jj)+ibld_ext ) + vmask(ji,jj-1,ibld(ji,jj)+ibld_ext ) ) 1218 zdb(ji,jj) = grav * zthermal * zdt(ji,jj) - grav * zbeta * zds(ji,jj) 1219 END_2D 1289 1220 END SUBROUTINE zdf_osm_vertical_average 1290 1221 … … 1306 1237 REAL(wp) :: ztemp 1307 1238 1308 DO jj = 2, jpjm1 1309 DO ji = 2, jpim1 1310 ztemp = zu(ji,jj) 1311 zu(ji,jj) = zu(ji,jj) * zcos_w(ji,jj) + zv(ji,jj) * zsin_w(ji,jj) 1312 zv(ji,jj) = zv(ji,jj) * zcos_w(ji,jj) - ztemp * zsin_w(ji,jj) 1313 ztemp = zdu(ji,jj) 1314 zdu(ji,jj) = zdu(ji,jj) * zcos_w(ji,jj) + zdv(ji,jj) * zsin_w(ji,jj) 1315 zdv(ji,jj) = zdv(ji,jj) * zsin_w(ji,jj) - ztemp * zsin_w(ji,jj) 1316 END DO 1317 END DO 1239 DO_2D_00_00 1240 ztemp = zu(ji,jj) 1241 zu(ji,jj) = zu(ji,jj) * zcos_w(ji,jj) + zv(ji,jj) * zsin_w(ji,jj) 1242 zv(ji,jj) = zv(ji,jj) * zcos_w(ji,jj) - ztemp * zsin_w(ji,jj) 1243 ztemp = zdu(ji,jj) 1244 zdu(ji,jj) = zdu(ji,jj) * zcos_w(ji,jj) + zdv(ji,jj) * zsin_w(ji,jj) 1245 zdv(ji,jj) = zdv(ji,jj) * zsin_w(ji,jj) - ztemp * zsin_w(ji,jj) 1246 END_2D 1318 1247 END SUBROUTINE zdf_osm_velocity_rotation 1319 1248 1320 SUBROUTINE zdf_osm_external_gradients( zdtdz, zdsdz, zdbdz )1249 SUBROUTINE zdf_osm_external_gradients( zdtdz, zdsdz, zdbdz, Kmm ) 1321 1250 !!--------------------------------------------------------------------- 1322 1251 !! *** ROUTINE zdf_osm_external_gradients *** … … 1328 1257 !!---------------------------------------------------------------------- 1329 1258 1259 INTEGER, INTENT(in) :: Kmm ! time-level index 1330 1260 REAL(wp), DIMENSION(jpi,jpj) :: zdtdz, zdsdz, zdbdz ! External gradients of temperature, salinity and buoyancy. 1331 1261 … … 1334 1264 1335 1265 1336 DO jj = 2, jpjm1 1337 DO ji = 2, jpim1 1338 IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN 1339 zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 1340 zbeta = rab_n(ji,jj,1,jp_sal) 1341 jkb = ibld(ji,jj) + ibld_ext 1342 jkb1 = MIN(jkb + 1, mbkt(ji,jj)) 1343 zdtdz(ji,jj) = - ( tsn(ji,jj,jkb1,jp_tem) - tsn(ji,jj,jkb,jp_tem ) ) & 1344 & / e3t_n(ji,jj,ibld(ji,jj)) 1345 zdsdz(ji,jj) = - ( tsn(ji,jj,jkb1,jp_sal) - tsn(ji,jj,jkb,jp_sal ) ) & 1346 & / e3t_n(ji,jj,ibld(ji,jj)) 1347 zdbdz(ji,jj) = grav * zthermal * zdtdz(ji,jj) - grav * zbeta * zdsdz(ji,jj) 1348 END IF 1349 END DO 1350 END DO 1266 DO_2D_00_00 1267 IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN 1268 zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 1269 zbeta = rab_n(ji,jj,1,jp_sal) 1270 jkb = ibld(ji,jj) + ibld_ext 1271 jkb1 = MIN(jkb + 1, mbkt(ji,jj)) 1272 zdtdz(ji,jj) = - ( ts(ji,jj,jkb1,jp_tem, Kmm ) - ts(ji,jj,jkb,jp_tem, Kmm ) ) & 1273 & / e3t(ji,jj,ibld(ji,jj), Kmm ) 1274 zdsdz(ji,jj) = - ( ts(ji,jj,jkb1,jp_sal, Kmm ) - ts(ji,jj,jkb,jp_sal, Kmm ) ) & 1275 & / e3t(ji,jj,ibld(ji,jj), Kmm ) 1276 zdbdz(ji,jj) = grav * zthermal * zdtdz(ji,jj) - grav * zbeta * zdsdz(ji,jj) 1277 END IF 1278 END_2D 1351 1279 END SUBROUTINE zdf_osm_external_gradients 1352 1280 1353 SUBROUTINE zdf_osm_pycnocline_scalar_profiles( zdtdz, zdsdz, zdbdz ) 1354 1281 SUBROUTINE zdf_osm_pycnocline_scalar_profiles( zdtdz, zdsdz, zdbdz, Kmm ) 1282 1283 INTEGER, INTENT(in) :: Kmm ! time-level index 1355 1284 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdtdz, zdsdz, zdbdz ! gradients in the pycnocline 1356 1285 … … 1360 1289 REAL(wp) :: zzeta_s=0.3, ztmp 1361 1290 1362 DO jj = 2, jpjm1 1363 DO ji = 2, jpim1 1364 IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN 1365 IF ( lconv(ji,jj) ) THEN ! convective conditions 1366 IF ( zdbdz_ext(ji,jj) > 0._wp .AND. & 1367 & (zdhdt(ji,jj) > 0._wp .AND. ln_osm_mle .AND. zdb_bl(ji,jj) > rn_osm_mle_thresh & 1368 & .OR. zdb_bl(ji,jj) > 0._wp)) THEN ! zdhdt could be <0 due to FK, hence check zdhdt>0 1291 DO_2D_00_00 1292 IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN 1293 IF ( lconv(ji,jj) ) THEN ! convective conditions 1294 IF ( zdbdz_ext(ji,jj) > 0._wp .AND. & 1295 & (zdhdt(ji,jj) > 0._wp .AND. ln_osm_mle .AND. zdb_bl(ji,jj) > rn_osm_mle_thresh & 1296 & .OR. zdb_bl(ji,jj) > 0._wp)) THEN ! zdhdt could be <0 due to FK, hence check zdhdt>0 1297 ztmp = 1._wp/MAX(zdh(ji,jj), epsln) 1298 ztgrad = 0.5 * zdt_ml(ji,jj) * ztmp + zdtdz_ext(ji,jj) 1299 zsgrad = 0.5 * zds_ml(ji,jj) * ztmp + zdsdz_ext(ji,jj) 1300 zbgrad = 0.5 * zdb_ml(ji,jj) * ztmp + zdbdz_ext(ji,jj) 1301 zgamma_b_nd = zdbdz_ext(ji,jj) * zdh(ji,jj) / MAX(zdb_ml(ji,jj), epsln) 1302 zgamma_c = ( 3.14159 / 4.0 ) * ( 0.5 + zgamma_b_nd ) /& 1303 & ( 1.0 - 0.25 * SQRT( 3.14159 / 6.0 ) - 2.0 * zgamma_b_nd * zzeta_s )**2 ! check 1304 DO jk = 2, ibld(ji,jj)+ibld_ext 1305 znd = -( gdepw(ji,jj,jk, Kmm ) - zhbl(ji,jj) ) * ztmp 1306 IF ( znd <= zzeta_s ) THEN 1307 zdtdz(ji,jj,jk) = zdtdz_ext(ji,jj) + 0.5 * zdt_ml(ji,jj) * ztmp * & 1308 & EXP( -6.0 * ( znd -zzeta_s )**2 ) 1309 zdsdz(ji,jj,jk) = zdsdz_ext(ji,jj) + 0.5 * zds_ml(ji,jj) * ztmp * & 1310 & EXP( -6.0 * ( znd -zzeta_s )**2 ) 1311 zdbdz(ji,jj,jk) = zdbdz_ext(ji,jj) + 0.5 * zdb_ml(ji,jj) * ztmp * & 1312 & EXP( -6.0 * ( znd -zzeta_s )**2 ) 1313 ELSE 1314 zdtdz(ji,jj,jk) = ztgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 ) 1315 zdsdz(ji,jj,jk) = zsgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 ) 1316 zdbdz(ji,jj,jk) = zbgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 ) 1317 ENDIF 1318 END DO 1319 ENDIF ! If condition not satisfied, no pycnocline present. Gradients have been initialised to zero, so do nothing 1320 ELSE 1321 ! stable conditions 1322 ! if pycnocline profile only defined when depth steady of increasing. 1323 IF ( zdhdt(ji,jj) >= 0.0 ) THEN ! Depth increasing, or steady. 1324 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 1325 IF ( zhol(ji,jj) >= 0.5 ) THEN ! Very stable - 'thick' pycnocline 1326 ztmp = 1._wp/MAX(zhbl(ji,jj), epsln) 1327 ztgrad = zdt_bl(ji,jj) * ztmp 1328 zsgrad = zds_bl(ji,jj) * ztmp 1329 zbgrad = zdb_bl(ji,jj) * ztmp 1330 DO jk = 2, ibld(ji,jj) 1331 znd = gdepw(ji,jj,jk, Kmm ) * ztmp 1332 zdtdz(ji,jj,jk) = ztgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 1333 zdbdz(ji,jj,jk) = zbgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 1334 zdsdz(ji,jj,jk) = zsgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 1335 END DO 1336 ELSE ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form. 1369 1337 ztmp = 1._wp/MAX(zdh(ji,jj), epsln) 1370 ztgrad = 0.5 * zdt_ml(ji,jj) * ztmp + zdtdz_ext(ji,jj) 1371 zsgrad = 0.5 * zds_ml(ji,jj) * ztmp + zdsdz_ext(ji,jj) 1372 zbgrad = 0.5 * zdb_ml(ji,jj) * ztmp + zdbdz_ext(ji,jj) 1373 zgamma_b_nd = zdbdz_ext(ji,jj) * zdh(ji,jj) / MAX(zdb_ml(ji,jj), epsln) 1374 zgamma_c = ( 3.14159 / 4.0 ) * ( 0.5 + zgamma_b_nd ) /& 1375 & ( 1.0 - 0.25 * SQRT( 3.14159 / 6.0 ) - 2.0 * zgamma_b_nd * zzeta_s )**2 ! check 1376 DO jk = 2, ibld(ji,jj)+ibld_ext 1377 znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) * ztmp 1378 IF ( znd <= zzeta_s ) THEN 1379 zdtdz(ji,jj,jk) = zdtdz_ext(ji,jj) + 0.5 * zdt_ml(ji,jj) * ztmp * & 1380 & EXP( -6.0 * ( znd -zzeta_s )**2 ) 1381 zdsdz(ji,jj,jk) = zdsdz_ext(ji,jj) + 0.5 * zds_ml(ji,jj) * ztmp * & 1382 & EXP( -6.0 * ( znd -zzeta_s )**2 ) 1383 zdbdz(ji,jj,jk) = zdbdz_ext(ji,jj) + 0.5 * zdb_ml(ji,jj) * ztmp * & 1384 & EXP( -6.0 * ( znd -zzeta_s )**2 ) 1385 ELSE 1386 zdtdz(ji,jj,jk) = ztgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 ) 1387 zdsdz(ji,jj,jk) = zsgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 ) 1388 zdbdz(ji,jj,jk) = zbgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 ) 1389 ENDIF 1338 ztgrad = zdt_bl(ji,jj) * ztmp 1339 zsgrad = zds_bl(ji,jj) * ztmp 1340 zbgrad = zdb_bl(ji,jj) * ztmp 1341 DO jk = 2, ibld(ji,jj) 1342 znd = -( gdepw(ji,jj,jk, Kmm ) - zhml(ji,jj) ) * ztmp 1343 zdtdz(ji,jj,jk) = ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 1344 zdbdz(ji,jj,jk) = zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 1345 zdsdz(ji,jj,jk) = zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 1390 1346 END DO 1391 ENDIF ! If condition not satisfied, no pycnocline present. Gradients have been initialised to zero, so do nothing 1392 ELSE 1393 ! stable conditions 1394 ! if pycnocline profile only defined when depth steady of increasing. 1395 IF ( zdhdt(ji,jj) >= 0.0 ) THEN ! Depth increasing, or steady. 1396 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 1397 IF ( zhol(ji,jj) >= 0.5 ) THEN ! Very stable - 'thick' pycnocline 1398 ztmp = 1._wp/MAX(zhbl(ji,jj), epsln) 1399 ztgrad = zdt_bl(ji,jj) * ztmp 1400 zsgrad = zds_bl(ji,jj) * ztmp 1401 zbgrad = zdb_bl(ji,jj) * ztmp 1402 DO jk = 2, ibld(ji,jj) 1403 znd = gdepw_n(ji,jj,jk) * ztmp 1404 zdtdz(ji,jj,jk) = ztgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 1405 zdbdz(ji,jj,jk) = zbgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 1406 zdsdz(ji,jj,jk) = zsgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 1407 END DO 1408 ELSE ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form. 1409 ztmp = 1._wp/MAX(zdh(ji,jj), epsln) 1410 ztgrad = zdt_bl(ji,jj) * ztmp 1411 zsgrad = zds_bl(ji,jj) * ztmp 1412 zbgrad = zdb_bl(ji,jj) * ztmp 1413 DO jk = 2, ibld(ji,jj) 1414 znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) * ztmp 1415 zdtdz(ji,jj,jk) = ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 1416 zdbdz(ji,jj,jk) = zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 1417 zdsdz(ji,jj,jk) = zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 1418 END DO 1419 ENDIF ! IF (zhol >=0.5) 1420 ENDIF ! IF (zdb_bl> 0.) 1421 ENDIF ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero 1422 ENDIF ! IF (lconv) 1423 END IF ! IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) 1424 END DO 1425 END DO 1347 ENDIF ! IF (zhol >=0.5) 1348 ENDIF ! IF (zdb_bl> 0.) 1349 ENDIF ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero 1350 ENDIF ! IF (lconv) 1351 END IF ! IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) 1352 END_2D 1426 1353 1427 1354 END SUBROUTINE zdf_osm_pycnocline_scalar_profiles 1428 1355 1429 SUBROUTINE zdf_osm_pycnocline_shear_profiles( zdudz, zdvdz )1356 SUBROUTINE zdf_osm_pycnocline_shear_profiles( zdudz, zdvdz, Kmm ) 1430 1357 !!--------------------------------------------------------------------- 1431 1358 !! *** ROUTINE zdf_osm_pycnocline_shear_profiles *** … … 1437 1364 !!---------------------------------------------------------------------- 1438 1365 1366 INTEGER, INTENT(in) :: Kmm ! time-level index 1439 1367 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdudz, zdvdz 1440 1368 … … 1443 1371 REAL(wp) :: zzeta_v = 0.45 1444 1372 ! 1445 DO jj = 2, jpjm1 1446 DO ji = 2, jpim1 1373 DO_2D_00_00 1374 ! 1375 IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN 1376 IF ( lconv (ji,jj) ) THEN 1377 ! Unstable conditions 1378 zugrad = 0.7 * zdu_ml(ji,jj) / zdh(ji,jj) + 0.3 * zustar(ji,jj)*zustar(ji,jj) / & 1379 & ( ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) * & 1380 & MIN(zla(ji,jj)**(8.0/3.0) + epsln, 0.12 )) 1381 !Alan is this right? 1382 zvgrad = ( 0.7 * zdv_ml(ji,jj) + & 1383 & 2.0 * ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) / & 1384 & ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird + epsln ) & 1385 & )/ (zdh(ji,jj) + epsln ) 1386 DO jk = 2, ibld(ji,jj) - 1 + ibld_ext 1387 znd = -( gdepw(ji,jj,jk, Kmm ) - zhbl(ji,jj) ) / (zdh(ji,jj) + epsln ) - zzeta_v 1388 IF ( znd <= 0.0 ) THEN 1389 zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( 3.0 * znd ) 1390 zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( 3.0 * znd ) 1391 ELSE 1392 zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( -2.0 * znd ) 1393 zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( -2.0 * znd ) 1394 ENDIF 1395 END DO 1396 ELSE 1397 ! stable conditions 1398 zugrad = 3.25 * zdu_bl(ji,jj) / zhbl(ji,jj) 1399 zvgrad = 2.75 * zdv_bl(ji,jj) / zhbl(ji,jj) 1400 DO jk = 2, ibld(ji,jj) 1401 znd = gdepw(ji,jj,jk, Kmm ) / zhbl(ji,jj) 1402 IF ( znd < 1.0 ) THEN 1403 zdudz(ji,jj,jk) = zugrad * EXP( -40.0 * ( znd - 1.0 )**2 ) 1404 ELSE 1405 zdudz(ji,jj,jk) = zugrad * EXP( -20.0 * ( znd - 1.0 )**2 ) 1406 ENDIF 1407 zdvdz(ji,jj,jk) = zvgrad * EXP( -20.0 * ( znd - 0.85 )**2 ) 1408 END DO 1409 ENDIF 1447 1410 ! 1448 IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN 1449 IF ( lconv (ji,jj) ) THEN 1450 ! Unstable conditions 1451 zugrad = 0.7 * zdu_ml(ji,jj) / zdh(ji,jj) + 0.3 * zustar(ji,jj)*zustar(ji,jj) / & 1452 & ( ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) * & 1453 & MIN(zla(ji,jj)**(8.0/3.0) + epsln, 0.12 )) 1454 !Alan is this right? 1455 zvgrad = ( 0.7 * zdv_ml(ji,jj) + & 1456 & 2.0 * ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) / & 1457 & ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird + epsln ) & 1458 & )/ (zdh(ji,jj) + epsln ) 1459 DO jk = 2, ibld(ji,jj) - 1 + ibld_ext 1460 znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / (zdh(ji,jj) + epsln ) - zzeta_v 1461 IF ( znd <= 0.0 ) THEN 1462 zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( 3.0 * znd ) 1463 zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( 3.0 * znd ) 1464 ELSE 1465 zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( -2.0 * znd ) 1466 zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( -2.0 * znd ) 1467 ENDIF 1468 END DO 1469 ELSE 1470 ! stable conditions 1471 zugrad = 3.25 * zdu_bl(ji,jj) / zhbl(ji,jj) 1472 zvgrad = 2.75 * zdv_bl(ji,jj) / zhbl(ji,jj) 1473 DO jk = 2, ibld(ji,jj) 1474 znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 1475 IF ( znd < 1.0 ) THEN 1476 zdudz(ji,jj,jk) = zugrad * EXP( -40.0 * ( znd - 1.0 )**2 ) 1477 ELSE 1478 zdudz(ji,jj,jk) = zugrad * EXP( -20.0 * ( znd - 1.0 )**2 ) 1479 ENDIF 1480 zdvdz(ji,jj,jk) = zvgrad * EXP( -20.0 * ( znd - 0.85 )**2 ) 1481 END DO 1482 ENDIF 1483 ! 1484 END IF ! IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) 1485 END DO 1486 END DO 1411 END IF ! IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) 1412 END_2D 1487 1413 END SUBROUTINE zdf_osm_pycnocline_shear_profiles 1488 1414 … … 1508 1434 REAL(wp) :: alpha_bc = 0.5 1509 1435 1510 DO jj = 2, jpjm1 1511 DO ji = 2, jpim1 1512 IF ( lconv(ji,jj) ) THEN ! Convective 1513 ! Alan is this right? Yes, it's a bit different from the previous relationship 1514 ! zwb_ent(ji,jj) = - 2.0 * 0.2 * zwbav(ji,jj) & 1515 ! & - ( 0.15 * ( 1.0 - EXP( -1.5 * zla(ji,jj) ) ) * zustar(ji,jj)**3 + 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj) 1516 zwcor = ABS(ff_t(ji,jj)) * zhbl(ji,jj) + epsln 1517 zrf_conv = TANH( ( zwstrc(ji,jj) / zwcor )**0.69 ) 1518 zrf_shear = TANH( ( zustar(ji,jj) / zwcor )**0.69 ) 1519 zrf_langmuir = TANH( ( zwstrl(ji,jj) / zwcor )**0.69 ) 1520 zr_stokes = 1.0 - EXP( -25.0 * dstokes(ji,jj) / hbl(ji,jj) & 1521 & * ( 1.0 + 4.0 * dstokes(ji,jj) / hbl(ji,jj) ) ) 1522 1523 zwb_ent(ji,jj) = - 2.0 * 0.2 * zrf_conv * zwbav(ji,jj) & 1524 & - 0.15 * zrf_shear * zustar(ji,jj)**3 /zhml(ji,jj) & 1525 & + zr_stokes * ( 0.15 * EXP( -1.5 * zla(ji,jj) ) * zrf_shear * zustar(ji,jj)**3 & 1526 & - zrf_langmuir * 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj) 1527 ! 1528 zwb_min = dh(ji,jj) * zwb0(ji,jj) / zhml(ji,jj) + zwb_ent(ji,jj) 1529 1530 IF ( ln_osm_mle ) THEN 1531 ! zwb_fk_b(ji,jj) = zwb_fk(ji,jj) * hmle(ji,jj) / ( 6.0 * hbl(ji,jj) ) * ( 6.0 * hbl(ji,jj) / hmle(ji,jj) - 1.0 + & 1532 ! & ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3 ) ! Fox-Kemper buoyancy flux average over OSBL 1533 IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN 1534 zwb_fk_b(ji,jj) = zwb_fk(ji,jj) * & 1535 (1.0 + hmle(ji,jj) / ( 6.0 * hbl(ji,jj) ) * (-1.0 + ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3) ) 1436 DO_2D_00_00 1437 IF ( lconv(ji,jj) ) THEN ! Convective 1438 ! Alan is this right? Yes, it's a bit different from the previous relationship 1439 ! zwb_ent(ji,jj) = - 2.0 * 0.2 * zwbav(ji,jj) & 1440 ! & - ( 0.15 * ( 1.0 - EXP( -1.5 * zla(ji,jj) ) ) * zustar(ji,jj)**3 + 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj) 1441 zwcor = ABS(ff_t(ji,jj)) * zhbl(ji,jj) + epsln 1442 zrf_conv = TANH( ( zwstrc(ji,jj) / zwcor )**0.69 ) 1443 zrf_shear = TANH( ( zustar(ji,jj) / zwcor )**0.69 ) 1444 zrf_langmuir = TANH( ( zwstrl(ji,jj) / zwcor )**0.69 ) 1445 zr_stokes = 1.0 - EXP( -25.0 * dstokes(ji,jj) / hbl(ji,jj) & 1446 & * ( 1.0 + 4.0 * dstokes(ji,jj) / hbl(ji,jj) ) ) 1447 1448 zwb_ent(ji,jj) = - 2.0 * 0.2 * zrf_conv * zwbav(ji,jj) & 1449 & - 0.15 * zrf_shear * zustar(ji,jj)**3 /zhml(ji,jj) & 1450 & + zr_stokes * ( 0.15 * EXP( -1.5 * zla(ji,jj) ) * zrf_shear * zustar(ji,jj)**3 & 1451 & - zrf_langmuir * 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj) 1452 ! 1453 zwb_min = dh(ji,jj) * zwb0(ji,jj) / zhml(ji,jj) + zwb_ent(ji,jj) 1454 1455 IF ( ln_osm_mle ) THEN 1456 ! zwb_fk_b(ji,jj) = zwb_fk(ji,jj) * hmle(ji,jj) / ( 6.0 * hbl(ji,jj) ) * ( 6.0 * hbl(ji,jj) / hmle(ji,jj) - 1.0 + & 1457 ! & ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3 ) ! Fox-Kemper buoyancy flux average over OSBL 1458 IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN 1459 zwb_fk_b(ji,jj) = zwb_fk(ji,jj) * & 1460 (1.0 + hmle(ji,jj) / ( 6.0 * hbl(ji,jj) ) * (-1.0 + ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3) ) 1461 ELSE 1462 zwb_fk_b(ji,jj) = 0.5 * zwb_fk(ji,jj) * hmle(ji,jj) / hbl(ji,jj) 1463 ENDIF 1464 zvel_max = ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**p2third / hbl(ji,jj) 1465 IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0.0 ) THEN 1466 ! OSBL is deepening, entrainment > restratification 1467 IF ( zdb_bl(ji,jj) > 0.0 .and. zdbdz_ext(ji,jj) > 0.0 ) THEN 1468 zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) 1536 1469 ELSE 1537 z wb_fk_b(ji,jj) = 0.5 * zwb_fk(ji,jj) * hmle(ji,jj) / hbl(ji,jj)1470 zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / MAX( zvel_max, 1.0e-15) 1538 1471 ENDIF 1539 zvel_max = ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**p2third / hbl(ji,jj)1540 IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0.0 ) THEN1541 ! OSBL is deepening, entrainment > restratification1542 IF ( zdb_bl(ji,jj) > 0.0 .and. zdbdz_ext(ji,jj) > 0.0 ) THEN1543 zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) )1544 ELSE1545 zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / MAX( zvel_max, 1.0e-15)1546 ENDIF1547 ELSE1548 ! OSBL shoaling due to restratification flux. This is the velocity defined in Fox-Kemper et al (2008)1549 zdhdt(ji,jj) = - zvel_mle(ji,jj)1550 1551 1552 ENDIF1553 1554 1472 ELSE 1555 ! Fox-Kemper not used. 1556 1557 zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / hbl(ji,jj) ) * zwb_ent(ji,jj) / & 1558 & ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 1559 zdhdt(ji,jj) = -zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) 1560 ! added ajgn 23 July as temporay fix 1473 ! OSBL shoaling due to restratification flux. This is the velocity defined in Fox-Kemper et al (2008) 1474 zdhdt(ji,jj) = - zvel_mle(ji,jj) 1475 1561 1476 1562 1477 ENDIF 1563 1478 1564 zdhdt_2(ji,jj) = 0._wp 1565 1566 ! commented out ajgn 23 July as temporay fix 1479 ELSE 1480 ! Fox-Kemper not used. 1481 1482 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) / & 1483 & ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 1484 zdhdt(ji,jj) = -zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) 1485 ! added ajgn 23 July as temporay fix 1486 1487 ENDIF 1488 1489 zdhdt_2(ji,jj) = 0._wp 1490 1491 ! commented out ajgn 23 July as temporay fix 1567 1492 ! IF ( zdb_ml(ji,jj) > 0.0 .and. zdbdz_ext(ji,jj) > 0.0 ) THEN 1568 1493 ! !additional term to account for thickness of pycnocline on dhdt. Small correction, so could get rid of this if necessary. … … 1580 1505 ! zdhdt(ji,jj) = zdhdt(ji,jj) + zdhdt_2(ji,jj) 1581 1506 ! ENDIF 1582 ELSE ! Stable 1583 zdhdt(ji,jj) = ( 0.06 + 0.52 * zhol(ji,jj) / 2.0 ) * zvstr(ji,jj)**3 / hbl(ji,jj) + zwbav(ji,jj) 1584 zdhdt_2(ji,jj) = 0._wp 1585 IF ( zdhdt(ji,jj) < 0._wp ) THEN 1586 ! For long timsteps factor in brackets slows the rapid collapse of the OSBL 1587 zpert = 2.0 * ( 1.0 + 0.0 * 2.0 * zvstr(ji,jj) * rn_rdt / hbl(ji,jj) ) * zvstr(ji,jj)**2 / hbl(ji,jj) 1588 ELSE 1589 zpert = MAX( 2.0 * ( 1.0 + 0.0 * 2.0 * zvstr(ji,jj) * rn_rdt / hbl(ji,jj) ) * zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) ) 1590 ENDIF 1591 zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / zpert 1592 ENDIF 1593 END DO 1594 END DO 1507 ELSE ! Stable 1508 zdhdt(ji,jj) = ( 0.06 + 0.52 * zhol(ji,jj) / 2.0 ) * zvstr(ji,jj)**3 / hbl(ji,jj) + zwbav(ji,jj) 1509 zdhdt_2(ji,jj) = 0._wp 1510 IF ( zdhdt(ji,jj) < 0._wp ) THEN 1511 ! For long timsteps factor in brackets slows the rapid collapse of the OSBL 1512 zpert = 2.0 * ( 1.0 + 0.0 * 2.0 * zvstr(ji,jj) * rn_Dt / hbl(ji,jj) ) * zvstr(ji,jj)**2 / hbl(ji,jj) 1513 ELSE 1514 zpert = MAX( 2.0 * ( 1.0 + 0.0 * 2.0 * zvstr(ji,jj) * rn_Dt / hbl(ji,jj) ) * zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) ) 1515 ENDIF 1516 zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / zpert 1517 ENDIF 1518 END_2D 1595 1519 END SUBROUTINE zdf_osm_calculate_dhdt 1596 1520 1597 SUBROUTINE zdf_osm_timestep_hbl( zdhdt, zdhdt_2 )1521 SUBROUTINE zdf_osm_timestep_hbl( zdhdt, zdhdt_2, Kmm ) 1598 1522 !!--------------------------------------------------------------------- 1599 1523 !! *** ROUTINE zdf_osm_timestep_hbl *** … … 1609 1533 1610 1534 1535 INTEGER, INTENT(in) :: Kmm ! time-level index 1611 1536 REAL(wp), DIMENSION(jpi,jpj) :: zdhdt, zdhdt_2 ! rates of change of hbl. 1612 1537 … … 1615 1540 REAL(wp) :: zthermal, zbeta 1616 1541 1617 DO jj = 2, jpjm1 1618 DO ji = 2, jpim1 1619 IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN 1542 DO_2D_00_00 1543 IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN 1620 1544 ! 1621 1545 ! If boundary layer changes by more than one level, need to check for stable layers between initial and final depths. 1622 1546 ! 1623 1624 1625 1626 1627 1628 1629 1547 zhbl_s = hbl(ji,jj) 1548 jm = imld(ji,jj) 1549 zthermal = rab_n(ji,jj,1,jp_tem) 1550 zbeta = rab_n(ji,jj,1,jp_sal) 1551 1552 1553 IF ( lconv(ji,jj) ) THEN 1630 1554 !unstable 1631 1555 1632 IF( ln_osm_mle ) THEN 1633 zvel_max = ( zwstrl(ji,jj)**3 + zwstrc(ji,jj)**3 )**p2third / hbl(ji,jj) 1556 IF( ln_osm_mle ) THEN 1557 zvel_max = ( zwstrl(ji,jj)**3 + zwstrc(ji,jj)**3 )**p2third / hbl(ji,jj) 1558 ELSE 1559 1560 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) / & 1561 & ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 1562 1563 ENDIF 1564 1565 DO jk = imld(ji,jj), ibld(ji,jj) 1566 zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - ts(ji,jj,jm,jp_tem,Kmm) ) & 1567 & - zbeta * ( zs_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ), & 1568 & 0.0 ) + zvel_max 1569 1570 1571 IF ( ln_osm_mle ) THEN 1572 zhbl_s = zhbl_s + MIN( & 1573 & rn_Dt * ( ( -zwb_ent(ji,jj) - 2.0 * zwb_fk_b(ji,jj) )/ zdb + zdhdt_2(ji,jj) ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), & 1574 & e3w(ji,jj,jm, Kmm ) ) 1634 1575 ELSE 1635 1636 zvel_max = -( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / hbl(ji,jj) ) * zwb_ent(ji,jj) / & 1637 & ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 1638 1576 zhbl_s = zhbl_s + MIN( & 1577 & rn_Dt * ( -zwb_ent(ji,jj) / zdb + zdhdt_2(ji,jj) ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), & 1578 & e3w(ji,jj,jm, Kmm ) ) 1639 1579 ENDIF 1640 1580 1641 DO jk = imld(ji,jj), ibld(ji,jj) 1642 zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) ) & 1643 & - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ), & 1644 & 0.0 ) + zvel_max 1645 1646 1647 IF ( ln_osm_mle ) THEN 1648 zhbl_s = zhbl_s + MIN( & 1649 & rn_rdt * ( ( -zwb_ent(ji,jj) - 2.0 * zwb_fk_b(ji,jj) )/ zdb + zdhdt_2(ji,jj) ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), & 1650 & e3w_n(ji,jj,jm) ) 1651 ELSE 1652 zhbl_s = zhbl_s + MIN( & 1653 & rn_rdt * ( -zwb_ent(ji,jj) / zdb + zdhdt_2(ji,jj) ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), & 1654 & e3w_n(ji,jj,jm) ) 1655 ENDIF 1656 1657 zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol) 1658 1659 IF ( zhbl_s >= gdepw_n(ji,jj,jm+1) ) jm = jm + 1 1660 END DO 1661 hbl(ji,jj) = zhbl_s 1662 ibld(ji,jj) = jm 1663 ELSE 1581 zhbl_s = MIN(zhbl_s, gdepw(ji,jj, mbkt(ji,jj) + 1, Kmm ) - depth_tol) 1582 1583 IF ( zhbl_s >= gdepw(ji,jj,jm+1, Kmm) ) jm = jm + 1 1584 END DO 1585 hbl(ji,jj) = zhbl_s 1586 ibld(ji,jj) = jm 1587 ELSE 1664 1588 ! stable 1665 1666 1667 & grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) )&1668 & - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ),&1669 1670 1671 1672 1673 1674 1675 1676 1677 zhbl_s = zhbl_s + MIN( zdhdt(ji,jj) / zdb * rn_rdt / FLOAT( ibld(ji,jj) - imld(ji,jj) ), e3w_n(ji,jj,jm) )1678 1679 zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol)1680 IF ( zhbl_s >= gdepw_n(ji,jj,jm) ) jm = jm + 11681 1682 1683 hbl(ji,jj) = MAX(zhbl_s, gdepw_n(ji,jj,4) )1684 1685 1589 DO jk = imld(ji,jj), ibld(ji,jj) 1590 zdb = MAX( & 1591 & grav * ( zthermal * ( zt_bl(ji,jj) - ts(ji,jj,jm,jp_tem, Kmm) )& 1592 & - zbeta * ( zs_bl(ji,jj) - ts(ji,jj,jm,jp_sal, Kmm) ) ),& 1593 & 0.0 ) + & 1594 & 2.0 * zvstr(ji,jj)**2 / zhbl_s 1595 1596 ! Alan is thuis right? I have simply changed hbli to hbl 1597 zhol(ji,jj) = -zhbl_s / ( ( zvstr(ji,jj)**3 + epsln )/ zwbav(ji,jj) ) 1598 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) ) ) * & 1599 & zustar(ji,jj)**3 / zhbl_s ) * ( 0.725 + 0.225 * EXP( -7.5 * zhol(ji,jj) ) ) 1600 zdhdt(ji,jj) = zdhdt(ji,jj) + zwbav(ji,jj) 1601 zhbl_s = zhbl_s + MIN( zdhdt(ji,jj) / zdb * rn_Dt / FLOAT( ibld(ji,jj) - imld(ji,jj) ), e3w(ji,jj,jm, Kmm ) ) 1602 1603 zhbl_s = MIN(zhbl_s, gdepw(ji,jj, mbkt(ji,jj) + 1, Kmm ) - depth_tol) 1604 IF ( zhbl_s >= gdepw(ji,jj,jm, Kmm ) ) jm = jm + 1 1605 END DO 1606 ENDIF ! IF ( lconv ) 1607 hbl(ji,jj) = MAX(zhbl_s, gdepw(ji,jj,4, Kmm ) ) 1608 ibld(ji,jj) = MAX(jm, 4 ) 1609 ELSE 1686 1610 ! change zero or one model level. 1687 hbl(ji,jj) = MAX(zhbl_t(ji,jj), gdepw_n(ji,jj,4) ) 1688 ENDIF 1689 zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj)) 1690 END DO 1691 END DO 1611 hbl(ji,jj) = MAX(zhbl_t(ji,jj), gdepw(ji,jj,4, Kmm ) ) 1612 ENDIF 1613 zhbl(ji,jj) = gdepw(ji,jj,ibld(ji,jj), Kmm ) 1614 END_2D 1692 1615 1693 1616 END SUBROUTINE zdf_osm_timestep_hbl 1694 1617 1695 SUBROUTINE zdf_osm_pycnocline_thickness( dh, zdh )1618 SUBROUTINE zdf_osm_pycnocline_thickness( dh, zdh, Kmm ) 1696 1619 !!--------------------------------------------------------------------- 1697 1620 !! *** ROUTINE zdf_osm_pycnocline_thickness *** … … 1707 1630 !!---------------------------------------------------------------------- 1708 1631 1632 INTEGER, INTENT(in) :: Kmm ! time-level index 1709 1633 REAL(wp), DIMENSION(jpi,jpj) :: dh, zdh ! pycnocline thickness. 1710 1634 ! … … 1714 1638 1715 1639 1716 DO jj = 2, jpjm1 1717 DO ji = 2, jpim1 1718 IF ( lconv(ji,jj) ) THEN 1719 1720 IF( ln_osm_mle ) THEN 1721 IF ( ( zwb_ent(ji,jj) + zwb_fk_b(ji,jj) ) < 0._wp ) THEN 1722 ! OSBL is deepening. Note wb_fk_b is zero if ln_osm_mle=F 1723 IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_ext(ji,jj) > 0._wp)THEN 1724 IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN ! near neutral stability 1725 zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / & 1726 (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zvstr(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 ) 1727 ELSE ! unstable 1728 zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / & 1729 (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zwstrc(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 ) 1730 ENDIF 1731 ztau = 0.2 * hbl(ji,jj) / (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird 1732 zddhdt = zari * hbl(ji,jj) 1733 ELSE 1734 ztau = 0.2 * hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 1735 zddhdt = 0.2 * hbl(ji,jj) 1736 ENDIF 1737 ELSE 1738 ztau = 0.2 * hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 1739 zddhdt = 0.2 * hbl(ji,jj) 1740 ENDIF 1741 ELSE ! ln_osm_mle 1640 DO_2D_00_00 1641 IF ( lconv(ji,jj) ) THEN 1642 1643 IF( ln_osm_mle ) THEN 1644 IF ( ( zwb_ent(ji,jj) + zwb_fk_b(ji,jj) ) < 0._wp ) THEN 1645 ! OSBL is deepening. Note wb_fk_b is zero if ln_osm_mle=F 1742 1646 IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_ext(ji,jj) > 0._wp)THEN 1743 1647 IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN ! near neutral stability … … 1748 1652 (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zwstrc(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 ) 1749 1653 ENDIF 1750 ztau = hbl(ji,jj) / (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird1654 ztau = 0.2 * hbl(ji,jj) / (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird 1751 1655 zddhdt = zari * hbl(ji,jj) 1752 1656 ELSE 1753 ztau = hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird1657 ztau = 0.2 * hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 1754 1658 zddhdt = 0.2 * hbl(ji,jj) 1755 1659 ENDIF 1756 1757 END IF ! ln_osm_mle 1758 1759 dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau ) + zddhdt * ( 1.0 - EXP( -rn_rdt / ztau ) ) 1760 ! Alan: this hml is never defined or used 1761 ELSE ! IF (lconv) 1762 ztau = hbl(ji,jj) / zvstr(ji,jj) 1763 IF ( zdhdt(ji,jj) >= 0.0 ) THEN ! probably shouldn't include wm here 1764 ! boundary layer deepening 1765 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 1766 ! pycnocline thickness set by stratification - use same relationship as for neutral conditions. 1767 zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) & 1768 & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01 , 0.2 ) 1769 zddhdt = MIN( zari, 0.2 ) * hbl(ji,jj) 1770 ELSE 1771 zddhdt = 0.2 * hbl(ji,jj) 1660 ELSE 1661 ztau = 0.2 * hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 1662 zddhdt = 0.2 * hbl(ji,jj) 1663 ENDIF 1664 ELSE ! ln_osm_mle 1665 IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_ext(ji,jj) > 0._wp)THEN 1666 IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN ! near neutral stability 1667 zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / & 1668 (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zvstr(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 ) 1669 ELSE ! unstable 1670 zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / & 1671 (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zwstrc(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 ) 1772 1672 ENDIF 1773 ELSE ! IF(dhdt < 0) 1673 ztau = hbl(ji,jj) / (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird 1674 zddhdt = zari * hbl(ji,jj) 1675 ELSE 1676 ztau = hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 1774 1677 zddhdt = 0.2 * hbl(ji,jj) 1775 ENDIF ! IF (dhdt >= 0) 1776 dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau )+ zddhdt * ( 1.0 - EXP( -rn_rdt / ztau ) ) 1777 IF ( zdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zddhdt ! can be a problem with dh>hbl for rapid collapse 1778 ! Alan: this hml is never defined or used -- do we need it? 1779 ENDIF ! IF (lconv) 1780 1781 hml(ji,jj) = hbl(ji,jj) - dh(ji,jj) 1782 inhml = MAX( INT( dh(ji,jj) / e3t_n(ji,jj,ibld(ji,jj)) ) , 1 ) 1783 imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 3) 1784 zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 1785 zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 1786 END DO 1787 END DO 1678 ENDIF 1679 1680 END IF ! ln_osm_mle 1681 1682 dh(ji,jj) = dh(ji,jj) * EXP( -rn_Dt / ztau ) + zddhdt * ( 1.0 - EXP( -rn_Dt / ztau ) ) 1683 ! Alan: this hml is never defined or used 1684 ELSE ! IF (lconv) 1685 ztau = hbl(ji,jj) / zvstr(ji,jj) 1686 IF ( zdhdt(ji,jj) >= 0.0 ) THEN ! probably shouldn't include wm here 1687 ! boundary layer deepening 1688 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 1689 ! pycnocline thickness set by stratification - use same relationship as for neutral conditions. 1690 zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) & 1691 & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01 , 0.2 ) 1692 zddhdt = MIN( zari, 0.2 ) * hbl(ji,jj) 1693 ELSE 1694 zddhdt = 0.2 * hbl(ji,jj) 1695 ENDIF 1696 ELSE ! IF(dhdt < 0) 1697 zddhdt = 0.2 * hbl(ji,jj) 1698 ENDIF ! IF (dhdt >= 0) 1699 dh(ji,jj) = dh(ji,jj) * EXP( -rn_Dt / ztau )+ zddhdt * ( 1.0 - EXP( -rn_Dt / ztau ) ) 1700 IF ( zdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zddhdt ! can be a problem with dh>hbl for rapid collapse 1701 ! Alan: this hml is never defined or used -- do we need it? 1702 ENDIF ! IF (lconv) 1703 1704 hml(ji,jj) = hbl(ji,jj) - dh(ji,jj) 1705 inhml = MAX( INT( dh(ji,jj) / e3t(ji,jj,ibld(ji,jj), Kmm) ) , 1 ) 1706 imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 3) 1707 zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj), Kmm) 1708 zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 1709 END_2D 1788 1710 1789 1711 END SUBROUTINE zdf_osm_pycnocline_thickness 1790 1712 1791 1713 1792 SUBROUTINE zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle )1714 SUBROUTINE zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle, Kmm ) 1793 1715 !!---------------------------------------------------------------------- 1794 1716 !! *** ROUTINE zdf_osm_horizontal_gradients *** … … 1801 1723 !! Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008 1802 1724 1803 1725 INTEGER, INTENT(in) :: Kmm ! time-level index 1804 1726 REAL(wp), DIMENSION(jpi,jpj) :: dbdx_mle, dbdy_mle ! MLE horiz gradients at u & v points 1805 1727 REAL(wp), DIMENSION(jpi,jpj) :: zmld ! == estimated FK BLD used for MLE horiz gradients == ! … … 1819 1741 mld_prof(:,:) = nlb10 ! Initialization to the number of w ocean point 1820 1742 zmld(:,:) = 0._wp ! here hmlp used as a dummy variable, integrating vertically N^2 1821 zN2_c = grav * rn_osm_mle_rho_c * r1_rau0 ! convert density criteria into N^2 criteria 1822 DO jk = nlb10, jpkm1 1823 DO jj = 1, jpj ! Mixed layer level: w-level 1824 DO ji = 1, jpi 1825 ikt = mbkt(ji,jj) 1826 zmld(ji,jj) = zmld(ji,jj) + MAX( rn2b(ji,jj,jk) , 0._wp ) * e3w_n(ji,jj,jk) 1827 IF( zmld(ji,jj) < zN2_c ) mld_prof(ji,jj) = MIN( jk , ikt ) + 1 ! Mixed layer level 1828 END DO 1829 END DO 1830 END DO 1831 DO jj = 1, jpj 1832 DO ji = 1, jpi 1833 mld_prof(ji,jj) = MAX(mld_prof(ji,jj),ibld(ji,jj)) 1834 zmld(ji,jj) = gdepw_n(ji,jj,mld_prof(ji,jj)) 1835 END DO 1836 END DO 1743 zN2_c = grav * rn_osm_mle_rho_c * r1_rho0 ! convert density criteria into N^2 criteria 1744 DO_3D_11_11( nlb10, jpkm1 ) ! Mixed layer level: w-level 1745 ikt = mbkt(ji,jj) 1746 zmld(ji,jj) = zmld(ji,jj) + MAX( rn2b(ji,jj,jk) , 0._wp ) * e3w(ji,jj,jk, Kmm ) 1747 IF( zmld(ji,jj) < zN2_c ) mld_prof(ji,jj) = MIN( jk , ikt ) + 1 ! Mixed layer level 1748 END_3D 1749 DO_2D_11_11 1750 mld_prof(ji,jj) = MAX(mld_prof(ji,jj),ibld(ji,jj)) 1751 zmld(ji,jj) = gdepw(ji,jj,mld_prof(ji,jj), Kmm ) 1752 END_2D 1837 1753 ! ensure mld_prof .ge. ibld 1838 1754 ! … … 1841 1757 ztm(:,:) = 0._wp 1842 1758 zsm(:,:) = 0._wp 1843 DO jk = 1, ikmax ! MLD and mean buoyancy and N2 over the mixed layer 1844 DO jj = 1, jpj 1845 DO ji = 1, jpi 1846 zc = e3t_n(ji,jj,jk) * REAL( MIN( MAX( 0, mld_prof(ji,jj)-jk ) , 1 ) ) ! zc being 0 outside the ML t-points 1847 ztm(ji,jj) = ztm(ji,jj) + zc * tsn(ji,jj,jk,jp_tem) 1848 zsm(ji,jj) = zsm(ji,jj) + zc * tsn(ji,jj,jk,jp_sal) 1849 END DO 1850 END DO 1851 END DO 1759 DO_3D_11_11( 1, ikmax ) ! MLD and mean buoyancy and N2 over the mixed layer 1760 zc = e3t(ji,jj,jk, Kmm ) * REAL( MIN( MAX( 0, mld_prof(ji,jj)-jk ) , 1 ) ) ! zc being 0 outside the ML t-points 1761 ztm(ji,jj) = ztm(ji,jj) + zc * ts(ji,jj,jk,jp_tem, Kmm ) 1762 zsm(ji,jj) = zsm(ji,jj) + zc * ts(ji,jj,jk,jp_sal, Kmm ) 1763 END_3D 1852 1764 ! average temperature and salinity. 1853 ztm(:,:) = ztm(:,:) / MAX( e3t _n(:,:,1), zmld(:,:) )1854 zsm(:,:) = zsm(:,:) / MAX( e3t _n(:,:,1), zmld(:,:) )1765 ztm(:,:) = ztm(:,:) / MAX( e3t(:,:,1, Kmm), zmld(:,:) ) 1766 zsm(:,:) = zsm(:,:) / MAX( e3t(:,:,1, Kmm), zmld(:,:) ) 1855 1767 ! calculate horizontal gradients at u & v points 1856 1768 1857 DO jj = 2, jpjm1 1858 DO ji = 1, jpim1 1859 zdtdx(ji,jj) = ( ztm(ji+1,jj) - ztm( ji,jj) ) * umask(ji,jj,1) / e1u(ji,jj) 1860 zdsdx(ji,jj) = ( zsm(ji+1,jj) - zsm( ji,jj) ) * umask(ji,jj,1) / e1u(ji,jj) 1861 zmld_midu(ji,jj) = 0.25_wp * (zmld(ji+1,jj) + zmld( ji,jj)) 1862 ztsm_midu(ji,jj,jp_tem) = 0.5_wp * ( ztm(ji+1,jj) + ztm( ji,jj) ) 1863 ztsm_midu(ji,jj,jp_sal) = 0.5_wp * ( zsm(ji+1,jj) + zsm( ji,jj) ) 1864 END DO 1865 END DO 1866 1867 DO jj = 1, jpjm1 1868 DO ji = 2, jpim1 1869 zdtdy(ji,jj) = ( ztm(ji,jj+1) - ztm( ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj) 1870 zdsdy(ji,jj) = ( zsm(ji,jj+1) - zsm( ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj) 1871 zmld_midv(ji,jj) = 0.25_wp * (zmld(ji,jj+1) + zmld( ji,jj)) 1872 ztsm_midv(ji,jj,jp_tem) = 0.5_wp * ( ztm(ji,jj+1) + ztm( ji,jj) ) 1873 ztsm_midv(ji,jj,jp_sal) = 0.5_wp * ( zsm(ji,jj+1) + zsm( ji,jj) ) 1874 END DO 1875 END DO 1876 1877 CALL eos_rab(ztsm_midu, zmld_midu, zabu) 1878 CALL eos_rab(ztsm_midv, zmld_midv, zabv) 1879 1880 DO jj = 2, jpjm1 1881 DO ji = 1, jpim1 1882 dbdx_mle(ji,jj) = grav*(zdtdx(ji,jj)*zabu(ji,jj,jp_tem) - zdsdx(ji,jj)*zabu(ji,jj,jp_sal)) 1883 END DO 1884 END DO 1885 DO jj = 1, jpjm1 1886 DO ji = 2, jpim1 1887 dbdy_mle(ji,jj) = grav*(zdtdy(ji,jj)*zabv(ji,jj,jp_tem) - zdsdy(ji,jj)*zabv(ji,jj,jp_sal)) 1888 END DO 1889 END DO 1769 DO_2D_00_10 1770 zdtdx(ji,jj) = ( ztm(ji+1,jj) - ztm( ji,jj) ) * umask(ji,jj,1) / e1u(ji,jj) 1771 zdsdx(ji,jj) = ( zsm(ji+1,jj) - zsm( ji,jj) ) * umask(ji,jj,1) / e1u(ji,jj) 1772 zmld_midu(ji,jj) = 0.25_wp * (zmld(ji+1,jj) + zmld( ji,jj)) 1773 ztsm_midu(ji,jj,jp_tem) = 0.5_wp * ( ztm(ji+1,jj) + ztm( ji,jj) ) 1774 ztsm_midu(ji,jj,jp_sal) = 0.5_wp * ( zsm(ji+1,jj) + zsm( ji,jj) ) 1775 END_2D 1776 1777 DO_2D_10_00 1778 zdtdy(ji,jj) = ( ztm(ji,jj+1) - ztm( ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj) 1779 zdsdy(ji,jj) = ( zsm(ji,jj+1) - zsm( ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj) 1780 zmld_midv(ji,jj) = 0.25_wp * (zmld(ji,jj+1) + zmld( ji,jj)) 1781 ztsm_midv(ji,jj,jp_tem) = 0.5_wp * ( ztm(ji,jj+1) + ztm( ji,jj) ) 1782 ztsm_midv(ji,jj,jp_sal) = 0.5_wp * ( zsm(ji,jj+1) + zsm( ji,jj) ) 1783 END_2D 1784 1785 CALL eos_rab(ztsm_midu, zmld_midu, zabu, Kmm ) 1786 CALL eos_rab(ztsm_midv, zmld_midv, zabv, Kmm ) 1787 1788 DO_2D_00_10 1789 dbdx_mle(ji,jj) = grav*(zdtdx(ji,jj)*zabu(ji,jj,jp_tem) - zdsdx(ji,jj)*zabu(ji,jj,jp_sal)) 1790 END_2D 1791 DO_2D_10_00 1792 dbdy_mle(ji,jj) = grav*(zdtdy(ji,jj)*zabv(ji,jj,jp_tem) - zdsdy(ji,jj)*zabv(ji,jj,jp_sal)) 1793 END_2D 1890 1794 1891 1795 END SUBROUTINE zdf_osm_zmld_horizontal_gradients 1892 SUBROUTINE zdf_osm_mle_parameters( hmle, zwb_fk, zvel_mle, zdiff_mle )1796 SUBROUTINE zdf_osm_mle_parameters( hmle, zwb_fk, zvel_mle, zdiff_mle, Kmm ) 1893 1797 !!---------------------------------------------------------------------- 1894 1798 !! *** ROUTINE zdf_osm_mle_parameters *** … … 1901 1805 !! Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008 1902 1806 1807 INTEGER, INTENT(in) :: Kmm ! time-level index 1903 1808 REAL(wp), DIMENSION(jpi,jpj) :: hmle, zwb_fk, zvel_mle, zdiff_mle 1904 1809 INTEGER :: ji, jj, jk ! dummy loop indices … … 1908 1813 1909 1814 mld_prof(:,:) = 4 1910 DO jk = 5, jpkm1 1911 DO jj = 2, jpjm1 1912 DO ji = 2, jpim1 1913 IF ( hmle(ji,jj) >= gdepw_n(ji,jj,jk) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk) 1914 END DO 1915 END DO 1916 END DO 1917 ! DO jj = 2, jpjm1 1918 ! DO ji = 1, jpim1 1919 ! zhmle(ji,jj) = gdepw_n(ji,jj,mld_prof(ji,jj)) 1920 ! END DO 1921 ! END DO 1815 DO_3D_00_00( 5, jpkm1 ) 1816 IF ( hmle(ji,jj) >= gdepw(ji,jj,jk, Kmm ) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk) 1817 END_3D 1818 ! DO_2D_00_10 1819 ! zhmle(ji,jj) = gdepw(ji,jj,mld_prof(ji,jj), Kmm ) 1820 ! END_2D 1922 1821 ! Timestep mixed layer eddy depth. 1923 DO jj = 2, jpjm1 1924 DO ji = 2, jpim1 1925 zdb_mle = grav * (rhop(ji,jj,mld_prof(ji,jj)) - rhop(ji,jj,ibld(ji,jj) )) * r1_rau0 ! check ALMG 1926 IF ( lconv(ji,jj) .and. ( zdb_bl(ji,jj) < rn_osm_mle_thresh .and. mld_prof(ji,jj) > ibld(ji,jj) .and. zdb_mle > 0.0 ) ) THEN 1927 hmle(ji,jj) = hmle(ji,jj) + zwb0(ji,jj) * rn_rdt / MAX( zdb_mle, rn_osm_mle_thresh ) ! MLE layer deepening through encroachment. Don't have a good maximum value for deepening, so use threshold buoyancy. 1822 DO_2D_00_00 1823 zdb_mle = grav * (rhop(ji,jj,mld_prof(ji,jj)) - rhop(ji,jj,ibld(ji,jj) )) * r1_rho0 ! check ALMG 1824 IF ( lconv(ji,jj) .and. ( zdb_bl(ji,jj) < rn_osm_mle_thresh .and. mld_prof(ji,jj) > ibld(ji,jj) .and. zdb_mle > 0.0 ) ) THEN 1825 hmle(ji,jj) = hmle(ji,jj) + zwb0(ji,jj) * rn_Dt / MAX( zdb_mle, rn_osm_mle_thresh ) ! MLE layer deepening through encroachment. Don't have a good maximum value for deepening, so use threshold buoyancy. 1826 ELSE 1827 ! MLE layer relaxes towards mixed layer depth on timescale tau_mle, or tau_mle/10 1828 ! IF ( hmle(ji,jj) > zmld(ji,jj) ) THEN 1829 ! hmle(ji,jj) = hmle(ji,jj) - ( hmle(ji,jj) - zmld(ji,jj) ) * rn_Dt / rn_osm_mle_tau 1830 ! ELSE 1831 ! hmle(ji,jj) = hmle(ji,jj) - 10.0 * ( hmle(ji,jj) - zmld(ji,jj) ) * rn_Dt / rn_osm_mle_tau ! fast relaxation if MLE layer shallower than MLD 1832 ! ENDIF 1833 IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN 1834 hmle(ji,jj) = hmle(ji,jj) - ( hmle(ji,jj) - hbl(ji,jj) ) * rn_Dt / rn_osm_mle_tau 1928 1835 ELSE 1929 ! MLE layer relaxes towards mixed layer depth on timescale tau_mle, or tau_mle/10 1930 ! IF ( hmle(ji,jj) > zmld(ji,jj) ) THEN 1931 ! hmle(ji,jj) = hmle(ji,jj) - ( hmle(ji,jj) - zmld(ji,jj) ) * rn_rdt / rn_osm_mle_tau 1932 ! ELSE 1933 ! hmle(ji,jj) = hmle(ji,jj) - 10.0 * ( hmle(ji,jj) - zmld(ji,jj) ) * rn_rdt / rn_osm_mle_tau ! fast relaxation if MLE layer shallower than MLD 1934 ! ENDIF 1935 IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN 1936 hmle(ji,jj) = hmle(ji,jj) - ( hmle(ji,jj) - hbl(ji,jj) ) * rn_rdt / rn_osm_mle_tau 1937 ELSE 1938 hmle(ji,jj) = hmle(ji,jj) - 10.0 * ( hmle(ji,jj) - hbl(ji,jj) ) * rn_rdt / rn_osm_mle_tau ! fast relaxation if MLE layer shallower than MLD 1939 ENDIF 1836 hmle(ji,jj) = hmle(ji,jj) - 10.0 * ( hmle(ji,jj) - hbl(ji,jj) ) * rn_Dt / rn_osm_mle_tau ! fast relaxation if MLE layer shallower than MLD 1940 1837 ENDIF 1941 hmle(ji,jj) = MIN(hmle(ji,jj), ht_n(ji,jj))1942 END DO1943 END DO1838 ENDIF 1839 hmle(ji,jj) = MIN(hmle(ji,jj), ht(ji,jj)) 1840 END_2D 1944 1841 1945 1842 mld_prof = 4 1946 DO jk = 5, jpkm1 1947 DO jj = 2, jpjm1 1948 DO ji = 2, jpim1 1949 IF ( hmle(ji,jj) >= gdepw_n(ji,jj,jk) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk) 1950 END DO 1951 END DO 1952 END DO 1953 DO jj = 2, jpjm1 1954 DO ji = 2, jpim1 1955 zhmle(ji,jj) = gdepw_n(ji,jj, mld_prof(ji,jj)) 1956 END DO 1957 END DO 1843 DO_3D_00_00( 5, jpkm1 ) 1844 IF ( hmle(ji,jj) >= gdepw(ji,jj,jk, Kmm ) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk) 1845 END_3D 1846 DO_2D_00_00 1847 zhmle(ji,jj) = gdepw(ji,jj, mld_prof(ji,jj), Kmm ) 1848 END_2D 1958 1849 ! Calculate vertical buoyancy, heat and salinity fluxes due to MLE. 1959 1850 1960 DO jj = 2, jpjm1 1961 DO ji = 2, jpim1 1962 IF ( lconv(ji,jj) ) THEN 1963 ztmp = r1_ft(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf 1964 ! zwt_fk(ji,jj) = 0.5_wp * ztmp * ( dbdx_mle(ji,jj) * zdtdx(ji,jj) + dbdy_mle(ji,jj) * zdtdy(ji,jj) & 1965 ! & + dbdx_mle(ji-1,jj) * zdtdx(ji-1,jj) + dbdy_mle(ji,jj-1) * zdtdy(ji,jj-1) ) 1966 ! zws_fk(ji,jj) = 0.5_wp * ztmp * ( dbdx_mle(ji,jj) * zdsdx(ji,jj) + dbdy_mle(ji,jj) * zdsdy(ji,jj) & 1967 ! & + dbdx_mle(ji-1,jj) * zdsdx(ji-1,jj) + dbdy_mle(ji,jj-1) * zdsdy(ji,jj-1) ) 1968 zdbds_mle = SQRT( 0.5_wp * ( dbdx_mle(ji,jj) * dbdx_mle(ji,jj) + dbdy_mle(ji,jj) * dbdy_mle(ji,jj) & 1969 & + dbdx_mle(ji-1,jj) * dbdx_mle(ji-1,jj) + dbdy_mle(ji,jj-1) * dbdy_mle(ji,jj-1) ) ) 1970 zwb_fk(ji,jj) = rn_osm_mle_ce * hmle(ji,jj) * hmle(ji,jj) *ztmp * zdbds_mle * zdbds_mle 1971 ! This vbelocity scale, defined in Fox-Kemper et al (2008), is needed for calculating dhdt. 1972 zvel_mle(ji,jj) = zdbds_mle * ztmp * hmle(ji,jj) * tmask(ji,jj,1) 1973 zdiff_mle(ji,jj) = 1.e-4_wp * ztmp * zdbds_mle * zhmle(ji,jj)**3 / rn_osm_mle_lf 1974 ENDIF 1975 END DO 1976 END DO 1851 DO_2D_00_00 1852 IF ( lconv(ji,jj) ) THEN 1853 ztmp = r1_ft(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf 1854 ! zwt_fk(ji,jj) = 0.5_wp * ztmp * ( dbdx_mle(ji,jj) * zdtdx(ji,jj) + dbdy_mle(ji,jj) * zdtdy(ji,jj) & 1855 ! & + dbdx_mle(ji-1,jj) * zdtdx(ji-1,jj) + dbdy_mle(ji,jj-1) * zdtdy(ji,jj-1) ) 1856 ! zws_fk(ji,jj) = 0.5_wp * ztmp * ( dbdx_mle(ji,jj) * zdsdx(ji,jj) + dbdy_mle(ji,jj) * zdsdy(ji,jj) & 1857 ! & + dbdx_mle(ji-1,jj) * zdsdx(ji-1,jj) + dbdy_mle(ji,jj-1) * zdsdy(ji,jj-1) ) 1858 zdbds_mle = SQRT( 0.5_wp * ( dbdx_mle(ji,jj) * dbdx_mle(ji,jj) + dbdy_mle(ji,jj) * dbdy_mle(ji,jj) & 1859 & + dbdx_mle(ji-1,jj) * dbdx_mle(ji-1,jj) + dbdy_mle(ji,jj-1) * dbdy_mle(ji,jj-1) ) ) 1860 zwb_fk(ji,jj) = rn_osm_mle_ce * hmle(ji,jj) * hmle(ji,jj) *ztmp * zdbds_mle * zdbds_mle 1861 ! This vbelocity scale, defined in Fox-Kemper et al (2008), is needed for calculating dhdt. 1862 zvel_mle(ji,jj) = zdbds_mle * ztmp * hmle(ji,jj) * tmask(ji,jj,1) 1863 zdiff_mle(ji,jj) = 1.e-4_wp * ztmp * zdbds_mle * zhmle(ji,jj)**3 / rn_osm_mle_lf 1864 ENDIF 1865 END_2D 1977 1866 END SUBROUTINE zdf_osm_mle_parameters 1978 1867 … … 1980 1869 1981 1870 1982 SUBROUTINE zdf_osm_init 1871 SUBROUTINE zdf_osm_init( Kmm ) 1983 1872 !!---------------------------------------------------------------------- 1984 1873 !! *** ROUTINE zdf_osm_init *** … … 1992 1881 !! ** input : Namlist namosm 1993 1882 !!---------------------------------------------------------------------- 1883 INTEGER, INTENT(in) :: Kmm ! time level index (middle) 1884 ! 1994 1885 INTEGER :: ios ! local integer 1995 1886 INTEGER :: ji, jj, jk ! dummy loop indices … … 2006 1897 !!---------------------------------------------------------------------- 2007 1898 ! 2008 REWIND( numnam_ref ) ! Namelist namzdf_osm in reference namelist : Osmosis ML model2009 1899 READ ( numnam_ref, namzdf_osm, IOSTAT = ios, ERR = 901) 2010 1900 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_osm in reference namelist' ) 2011 1901 2012 REWIND( numnam_cfg ) ! Namelist namzdf_tke in configuration namelist : Turbulent Kinetic Energy2013 1902 READ ( numnam_cfg, namzdf_osm, IOSTAT = ios, ERR = 902 ) 2014 1903 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namzdf_osm in configuration namelist' ) … … 2049 1938 IF( ln_osm_mle ) THEN 2050 1939 ! Initialise Fox-Kemper parametrization 2051 REWIND( numnam_ref )! Namelist namosm_mle in reference namelist : Tracer advection scheme1940 ! Namelist namosm_mle in reference namelist : Tracer advection scheme 2052 1941 READ ( numnam_ref, namosm_mle, IOSTAT = ios, ERR = 903) 2053 1942 903 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namosm_mle in reference namelist') 2054 1943 2055 REWIND( numnam_cfg )! Namelist namosm_mle in configuration namelist : Tracer advection scheme1944 ! Namelist namosm_mle in configuration namelist : Tracer advection scheme 2056 1945 READ ( numnam_cfg, namosm_mle, IOSTAT = ios, ERR = 904 ) 2057 1946 904 IF( ios > 0 ) CALL ctl_nam ( ios , 'namosm_mle in configuration namelist') … … 2087 1976 IF( ln_osm_mle ) THEN ! MLE initialisation 2088 1977 ! 2089 rb_c = grav * rn_osm_mle_rho_c /r au0 ! Mixed Layer buoyancy criteria1978 rb_c = grav * rn_osm_mle_rho_c /rho0 ! Mixed Layer buoyancy criteria 2090 1979 IF(lwp) WRITE(numout,*) 2091 1980 IF(lwp) WRITE(numout,*) ' ML buoyancy criteria = ', rb_c, ' m/s2 ' … … 2110 1999 ENDIF 2111 2000 2112 call osm_rst( nit000, 'READ' ) !* read or initialize hbl, dh, hmle2001 call osm_rst( nit000, Kmm, 'READ' ) !* read or initialize hbl, dh, hmle 2113 2002 2114 2003 … … 2147 2036 etmean(:,:,:) = 0.e0 2148 2037 2149 DO jk = 1, jpkm1 2150 DO jj = 2, jpjm1 2151 DO ji = 2, jpim1 ! vector opt. 2152 etmean(ji,jj,jk) = tmask(ji,jj,jk) & 2153 & / MAX( 1., umask(ji-1,jj ,jk) + umask(ji,jj,jk) & 2154 & + vmask(ji ,jj-1,jk) + vmask(ji,jj,jk) ) 2155 END DO 2156 END DO 2157 END DO 2038 DO_3D_00_00( 1, jpkm1 ) 2039 etmean(ji,jj,jk) = tmask(ji,jj,jk) & 2040 & / MAX( 1., umask(ji-1,jj ,jk) + umask(ji,jj,jk) & 2041 & + vmask(ji ,jj-1,jk) + vmask(ji,jj,jk) ) 2042 END_3D 2158 2043 2159 2044 CASE ( 1 ) ! horizontal average … … 2165 2050 etmean(:,:,:) = 0.e0 2166 2051 2167 DO jk = 1, jpkm1 2168 DO jj = 2, jpjm1 2169 DO ji = 2, jpim1 ! vector opt. 2170 etmean(ji,jj,jk) = tmask(ji, jj,jk) & 2171 & / MAX( 1., 2.* tmask(ji,jj,jk) & 2172 & +.5 * ( tmask(ji-1,jj+1,jk) + tmask(ji-1,jj-1,jk) & 2173 & +tmask(ji+1,jj+1,jk) + tmask(ji+1,jj-1,jk) ) & 2174 & +1. * ( tmask(ji-1,jj ,jk) + tmask(ji ,jj+1,jk) & 2175 & +tmask(ji ,jj-1,jk) + tmask(ji+1,jj ,jk) ) ) 2176 END DO 2177 END DO 2178 END DO 2052 DO_3D_00_00( 1, jpkm1 ) 2053 etmean(ji,jj,jk) = tmask(ji, jj,jk) & 2054 & / MAX( 1., 2.* tmask(ji,jj,jk) & 2055 & +.5 * ( tmask(ji-1,jj+1,jk) + tmask(ji-1,jj-1,jk) & 2056 & +tmask(ji+1,jj+1,jk) + tmask(ji+1,jj-1,jk) ) & 2057 & +1. * ( tmask(ji-1,jj ,jk) + tmask(ji ,jj+1,jk) & 2058 & +tmask(ji ,jj-1,jk) + tmask(ji+1,jj ,jk) ) ) 2059 END_3D 2179 2060 2180 2061 CASE DEFAULT … … 2208 2089 2209 2090 2210 SUBROUTINE osm_rst( kt, cdrw )2091 SUBROUTINE osm_rst( kt, Kmm, cdrw ) 2211 2092 !!--------------------------------------------------------------------- 2212 2093 !! *** ROUTINE osm_rst *** … … 2218 2099 !!---------------------------------------------------------------------- 2219 2100 2220 INTEGER, INTENT(in) :: kt 2101 INTEGER , INTENT(in) :: kt ! ocean time step index 2102 INTEGER , INTENT(in) :: Kmm ! ocean time level index (middle) 2221 2103 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 2222 2104 … … 2236 2118 id1 = iom_varid( numror, 'wn' , ldstop = .FALSE. ) 2237 2119 IF( id1 > 0 ) THEN ! 'wn' exists; read 2238 CALL iom_get( numror, jpdom_autoglo, 'wn', w n, ldxios = lrxios )2239 WRITE(numout,*) ' ===>>>> : w nread from restart file'2120 CALL iom_get( numror, jpdom_autoglo, 'wn', ww, ldxios = lrxios ) 2121 WRITE(numout,*) ' ===>>>> : ww read from restart file' 2240 2122 ELSE 2241 w n(:,:,:) = 0._wp2242 WRITE(numout,*) ' ===>>>> : w nnot in restart file, set to zero initially'2123 ww(:,:,:) = 0._wp 2124 WRITE(numout,*) ' ===>>>> : ww not in restart file, set to zero initially' 2243 2125 END IF 2244 2126 … … 2270 2152 IF( TRIM(cdrw) == 'WRITE') THEN !* Write hbli into the restart file, then return 2271 2153 IF(lwp) WRITE(numout,*) '---- osm-rst ----' 2272 CALL iom_rstput( kt, nitrst, numrow, 'wn' , w n,ldxios = lwxios )2154 CALL iom_rstput( kt, nitrst, numrow, 'wn' , ww , ldxios = lwxios ) 2273 2155 CALL iom_rstput( kt, nitrst, numrow, 'hbl' , hbl, ldxios = lwxios ) 2274 2156 CALL iom_rstput( kt, nitrst, numrow, 'dh' , dh, ldxios = lwxios ) … … 2284 2166 IF( lwp ) WRITE(numout,*) ' ===>>>> : calculating hbl computed from stratification' 2285 2167 ! w-level of the mixing and mixed layers 2286 CALL eos_rab( ts n, rab_n)2287 CALL bn2(ts n, rab_n, rn2)2168 CALL eos_rab( ts(:,:,:,:,Kmm), rab_n, Kmm ) 2169 CALL bn2(ts(:,:,:,:,Kmm), rab_n, rn2, Kmm) 2288 2170 imld_rst(:,:) = nlb10 ! Initialization to the number of w ocean point 2289 2171 hbl(:,:) = 0._wp ! here hbl used as a dummy variable, integrating vertically N^2 2290 zN2_c = grav * rho_c * r1_r au0 ! convert density criteria into N^2 criteria2172 zN2_c = grav * rho_c * r1_rho0 ! convert density criteria into N^2 criteria 2291 2173 ! 2292 2174 hbl(:,:) = 0._wp ! here hbl used as a dummy variable, integrating vertically N^2 2293 DO jk = 1, jpkm1 2294 DO jj = 1, jpj ! Mixed layer level: w-level 2295 DO ji = 1, jpi 2296 ikt = mbkt(ji,jj) 2297 hbl(ji,jj) = hbl(ji,jj) + MAX( rn2(ji,jj,jk) , 0._wp ) * e3w_n(ji,jj,jk) 2298 IF( hbl(ji,jj) < zN2_c ) imld_rst(ji,jj) = MIN( jk , ikt ) + 1 ! Mixed layer level 2299 END DO 2300 END DO 2301 END DO 2175 DO_3D_11_11( 1, jpkm1 ) 2176 ikt = mbkt(ji,jj) 2177 hbl(ji,jj) = hbl(ji,jj) + MAX( rn2(ji,jj,jk) , 0._wp ) * e3w(ji,jj,jk,Kmm) 2178 IF( hbl(ji,jj) < zN2_c ) imld_rst(ji,jj) = MIN( jk , ikt ) + 1 ! Mixed layer level 2179 END_3D 2302 2180 ! 2303 DO jj = 1, jpj 2304 DO ji = 1, jpi 2305 iiki = MAX(4,imld_rst(ji,jj)) 2306 hbl (ji,jj) = gdepw_n(ji,jj,iiki ) ! Turbocline depth 2307 dh (ji,jj) = e3t_n(ji,jj,iiki-1 ) ! Turbocline depth 2308 END DO 2309 END DO 2181 DO_2D_11_11 2182 iiki = MAX(4,imld_rst(ji,jj)) 2183 hbl (ji,jj) = gdepw(ji,jj,iiki, Kmm ) ! Turbocline depth 2184 dh (ji,jj) = e3t(ji,jj,iiki-1, Kmm ) ! Turbocline depth 2185 END_2D 2310 2186 2311 2187 WRITE(numout,*) ' ===>>>> : hbl computed from stratification' … … 2316 2192 END IF 2317 2193 2318 w n(:,:,:) = 0._wp2194 ww(:,:,:) = 0._wp 2319 2195 WRITE(numout,*) ' ===>>>> : wn not in restart file, set to zero initially' 2320 2196 END SUBROUTINE osm_rst 2321 2197 2322 2198 2323 SUBROUTINE tra_osm( kt )2199 SUBROUTINE tra_osm( kt, Kmm, pts, Krhs ) 2324 2200 !!---------------------------------------------------------------------- 2325 2201 !! *** ROUTINE tra_osm *** … … 2331 2207 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdt, ztrds ! 3D workspace 2332 2208 !!---------------------------------------------------------------------- 2333 INTEGER, INTENT(in) :: kt 2209 INTEGER , INTENT(in) :: kt ! time step index 2210 INTEGER , INTENT(in) :: Kmm, Krhs ! time level indices 2211 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers and RHS of tracer equation 2212 ! 2334 2213 INTEGER :: ji, jj, jk 2335 2214 ! … … 2341 2220 2342 2221 IF( l_trdtra ) THEN !* Save ta and sa trends 2343 ALLOCATE( ztrdt(jpi,jpj,jpk) ) ; ztrdt(:,:,:) = tsa(:,:,:,jp_tem)2344 ALLOCATE( ztrds(jpi,jpj,jpk) ) ; ztrds(:,:,:) = tsa(:,:,:,jp_sal)2222 ALLOCATE( ztrdt(jpi,jpj,jpk) ) ; ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) 2223 ALLOCATE( ztrds(jpi,jpj,jpk) ) ; ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) 2345 2224 ENDIF 2346 2225 2347 DO jk = 1, jpkm1 2348 DO jj = 2, jpjm1 2349 DO ji = 2, jpim1 2350 tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) & 2351 & - ( ghamt(ji,jj,jk ) & 2352 & - ghamt(ji,jj,jk+1) ) /e3t_n(ji,jj,jk) 2353 tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) & 2354 & - ( ghams(ji,jj,jk ) & 2355 & - ghams(ji,jj,jk+1) ) / e3t_n(ji,jj,jk) 2356 END DO 2357 END DO 2358 END DO 2226 DO_3D_00_00( 1, jpkm1 ) 2227 pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) & 2228 & - ( ghamt(ji,jj,jk ) & 2229 & - ghamt(ji,jj,jk+1) ) /e3t(ji,jj,jk,Kmm) 2230 pts(ji,jj,jk,jp_sal,Krhs) = pts(ji,jj,jk,jp_sal,Krhs) & 2231 & - ( ghams(ji,jj,jk ) & 2232 & - ghams(ji,jj,jk+1) ) / e3t(ji,jj,jk,Kmm) 2233 END_3D 2359 2234 2360 2235 ! save the non-local tracer flux trends for diagnostics 2361 2236 IF( l_trdtra ) THEN 2362 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 2363 ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 2364 2365 CALL trd_tra( kt, 'TRA', jp_tem, jptra_osm, ztrdt ) 2366 CALL trd_tra( kt, 'TRA', jp_sal, jptra_osm, ztrds ) 2237 ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:) 2238 ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) - ztrds(:,:,:) 2239 CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_osm, ztrdt ) 2240 CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_osm, ztrds ) 2367 2241 DEALLOCATE( ztrdt ) ; DEALLOCATE( ztrds ) 2368 2242 ENDIF 2369 2243 2370 IF( ln_ctl) THEN2371 CALL prt_ctl( tab3d_1= tsa(:,:,:,jp_tem), clinfo1=' osm - Ta: ', mask1=tmask, &2372 & tab3d_2= tsa(:,:,:,jp_sal), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' )2244 IF(sn_cfctl%l_prtctl) THEN 2245 CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' osm - Ta: ', mask1=tmask, & 2246 & tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) 2373 2247 ENDIF 2374 2248 ! … … 2393 2267 2394 2268 2395 SUBROUTINE dyn_osm( kt )2269 SUBROUTINE dyn_osm( kt, Kmm, puu, pvv, Krhs ) 2396 2270 !!---------------------------------------------------------------------- 2397 2271 !! *** ROUTINE dyn_osm *** … … 2402 2276 !! ** Method : ??? 2403 2277 !!---------------------------------------------------------------------- 2404 INTEGER, INTENT(in) :: kt ! 2278 INTEGER , INTENT( in ) :: kt ! ocean time step index 2279 INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices 2280 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 2405 2281 ! 2406 2282 INTEGER :: ji, jj, jk ! dummy loop indices … … 2414 2290 !code saving tracer trends removed, replace with trdmxl_oce 2415 2291 2416 DO jk = 1, jpkm1 ! add non-local u and v fluxes 2417 DO jj = 2, jpjm1 2418 DO ji = 2, jpim1 2419 ua(ji,jj,jk) = ua(ji,jj,jk) & 2420 & - ( ghamu(ji,jj,jk ) & 2421 & - ghamu(ji,jj,jk+1) ) / e3u_n(ji,jj,jk) 2422 va(ji,jj,jk) = va(ji,jj,jk) & 2423 & - ( ghamv(ji,jj,jk ) & 2424 & - ghamv(ji,jj,jk+1) ) / e3v_n(ji,jj,jk) 2425 END DO 2426 END DO 2427 END DO 2292 DO_3D_00_00( 1, jpkm1 ) 2293 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) & 2294 & - ( ghamu(ji,jj,jk ) & 2295 & - ghamu(ji,jj,jk+1) ) / e3u(ji,jj,jk,Kmm) 2296 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) & 2297 & - ( ghamv(ji,jj,jk ) & 2298 & - ghamv(ji,jj,jk+1) ) / e3v(ji,jj,jk,Kmm) 2299 END_3D 2428 2300 ! 2429 2301 ! code for saving tracer trends removed
Note: See TracChangeset
for help on using the changeset viewer.
