- Timestamp:
- 2020-09-14T17:40:34+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11351_fldread_with_XIOS
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/OCE/ZDF/zdfosm.F90 (modified) (40 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11351_fldread_with_XIOS
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev @HEADext/AGRIF5 ^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 9 # SETTE 10 ^/utils/CI/sette@13382 sette
-
- Property svn:externals
-
NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/OCE/ZDF/zdfosm.F90
r11405 r13463 42 42 !!---------------------------------------------------------------------- 43 43 USE oce ! ocean dynamics and active tracers 44 ! uses w nfrom previous time step (which is now wb) to calculate hbl44 ! uses ww from previous time step (which is now wb) to calculate hbl 45 45 USE dom_oce ! ocean space and time domain 46 46 USE zdf_oce ! ocean vertical physics … … 103 103 INTEGER :: idebug = 236 104 104 INTEGER :: jdebug = 228 105 106 !! * Substitutions 107 # include "do_loop_substitute.h90" 108 # include "domzgr_substitute.h90" 105 109 !!---------------------------------------------------------------------- 106 110 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 122 126 123 127 124 SUBROUTINE zdf_osm( kt, p_avm, p_avt )128 SUBROUTINE zdf_osm( kt, Kbb, Kmm, Krhs, p_avm, p_avt ) 125 129 !!---------------------------------------------------------------------- 126 130 !! *** ROUTINE zdf_osm *** … … 157 161 !! the equation number. (LMD94, here after) 158 162 !!---------------------------------------------------------------------- 159 INTEGER , INTENT(in ) :: kt ! ocean time step 163 INTEGER , INTENT(in ) :: kt ! ocean time step 164 INTEGER , INTENT(in ) :: Kbb, Kmm, Krhs ! ocean time level indices 160 165 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: p_avm, p_avt ! momentum and tracer Kz (w-points) 161 166 !! … … 295 300 zz0 = rn_abs ! surface equi-partition in 2-bands 296 301 zz1 = 1. - rn_abs 297 DO jj = 2, jpjm1 298 DO ji = 2, jpim1 299 ! Surface downward irradiance (so always +ve) 300 zrad0(ji,jj) = qsr(ji,jj) * r1_rau0_rcp 301 ! Downwards irradiance at base of boundary layer 302 zradh(ji,jj) = zrad0(ji,jj) * ( zz0 * EXP( -hbl(ji,jj)/rn_si0 ) + zz1 * EXP( -hbl(ji,jj)/rn_si1) ) 303 ! Downwards irradiance averaged over depth of the OSBL 304 zradav(ji,jj) = zrad0(ji,jj) * ( zz0 * ( 1.0 - EXP( -hbl(ji,jj)/rn_si0 ) )*rn_si0 & 305 & + zz1 * ( 1.0 - EXP( -hbl(ji,jj)/rn_si1 ) )*rn_si1 ) / hbl(ji,jj) 306 END DO 307 END DO 302 DO_2D( 0, 0, 0, 0 ) 303 ! Surface downward irradiance (so always +ve) 304 zrad0(ji,jj) = qsr(ji,jj) * r1_rho0_rcp 305 ! Downwards irradiance at base of boundary layer 306 zradh(ji,jj) = zrad0(ji,jj) * ( zz0 * EXP( -hbl(ji,jj)/rn_si0 ) + zz1 * EXP( -hbl(ji,jj)/rn_si1) ) 307 ! Downwards irradiance averaged over depth of the OSBL 308 zradav(ji,jj) = zrad0(ji,jj) * ( zz0 * ( 1.0 - EXP( -hbl(ji,jj)/rn_si0 ) )*rn_si0 & 309 & + zz1 * ( 1.0 - EXP( -hbl(ji,jj)/rn_si1 ) )*rn_si1 ) / hbl(ji,jj) 310 END_2D 308 311 ! Turbulent surface fluxes and fluxes averaged over depth of the OSBL 309 DO jj = 2, jpjm1 310 DO ji = 2, jpim1 311 zthermal = rab_n(ji,jj,1,jp_tem) 312 zbeta = rab_n(ji,jj,1,jp_sal) 313 ! Upwards surface Temperature flux for non-local term 314 zwth0(ji,jj) = - qns(ji,jj) * r1_rau0_rcp * tmask(ji,jj,1) 315 ! Upwards surface salinity flux for non-local term 316 zws0(ji,jj) = - ( ( emp(ji,jj)-rnf(ji,jj) ) * tsn(ji,jj,1,jp_sal) + sfx(ji,jj) ) * r1_rau0 * tmask(ji,jj,1) 317 ! Non radiative upwards surface buoyancy flux 318 zwb0(ji,jj) = grav * zthermal * zwth0(ji,jj) - grav * zbeta * zws0(ji,jj) 319 ! turbulent heat flux averaged over depth of OSBL 320 zwthav(ji,jj) = 0.5 * zwth0(ji,jj) - ( 0.5*( zrad0(ji,jj) + zradh(ji,jj) ) - zradav(ji,jj) ) 321 ! turbulent salinity flux averaged over depth of the OBSL 322 zwsav(ji,jj) = 0.5 * zws0(ji,jj) 323 ! turbulent buoyancy flux averaged over the depth of the OBSBL 324 zwbav(ji,jj) = grav * zthermal * zwthav(ji,jj) - grav * zbeta * zwsav(ji,jj) 325 ! Surface upward velocity fluxes 326 zuw0(ji,jj) = -utau(ji,jj) * r1_rau0 * tmask(ji,jj,1) 327 zvw0(ji,jj) = -vtau(ji,jj) * r1_rau0 * tmask(ji,jj,1) 328 ! Friction velocity (zustar), at T-point : LMD94 eq. 2 329 zustar(ji,jj) = MAX( SQRT( SQRT( zuw0(ji,jj) * zuw0(ji,jj) + zvw0(ji,jj) * zvw0(ji,jj) ) ), 1.0e-8 ) 330 zcos_wind(ji,jj) = -zuw0(ji,jj) / ( zustar(ji,jj) * zustar(ji,jj) ) 331 zsin_wind(ji,jj) = -zvw0(ji,jj) / ( zustar(ji,jj) * zustar(ji,jj) ) 332 END DO 333 END DO 312 DO_2D( 0, 0, 0, 0 ) 313 zthermal = rab_n(ji,jj,1,jp_tem) 314 zbeta = rab_n(ji,jj,1,jp_sal) 315 ! Upwards surface Temperature flux for non-local term 316 zwth0(ji,jj) = - qns(ji,jj) * r1_rho0_rcp * tmask(ji,jj,1) 317 ! Upwards surface salinity flux for non-local term 318 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) 319 ! Non radiative upwards surface buoyancy flux 320 zwb0(ji,jj) = grav * zthermal * zwth0(ji,jj) - grav * zbeta * zws0(ji,jj) 321 ! turbulent heat flux averaged over depth of OSBL 322 zwthav(ji,jj) = 0.5 * zwth0(ji,jj) - ( 0.5*( zrad0(ji,jj) + zradh(ji,jj) ) - zradav(ji,jj) ) 323 ! turbulent salinity flux averaged over depth of the OBSL 324 zwsav(ji,jj) = 0.5 * zws0(ji,jj) 325 ! turbulent buoyancy flux averaged over the depth of the OBSBL 326 zwbav(ji,jj) = grav * zthermal * zwthav(ji,jj) - grav * zbeta * zwsav(ji,jj) 327 ! Surface upward velocity fluxes 328 zuw0(ji,jj) = -utau(ji,jj) * r1_rho0 * tmask(ji,jj,1) 329 zvw0(ji,jj) = -vtau(ji,jj) * r1_rho0 * tmask(ji,jj,1) 330 ! Friction velocity (zustar), at T-point : LMD94 eq. 2 331 zustar(ji,jj) = MAX( SQRT( SQRT( zuw0(ji,jj) * zuw0(ji,jj) + zvw0(ji,jj) * zvw0(ji,jj) ) ), 1.0e-8 ) 332 zcos_wind(ji,jj) = -zuw0(ji,jj) / ( zustar(ji,jj) * zustar(ji,jj) ) 333 zsin_wind(ji,jj) = -zvw0(ji,jj) / ( zustar(ji,jj) * zustar(ji,jj) ) 334 END_2D 334 335 ! Calculate Stokes drift in direction of wind (zustke) and Stokes penetration depth (dstokes) 335 336 SELECT CASE (nn_osm_wave) 336 337 ! Assume constant La#=0.3 337 338 CASE(0) 338 DO jj = 2, jpjm1 339 DO ji = 2, jpim1 340 zus_x = zcos_wind(ji,jj) * zustar(ji,jj) / 0.3**2 341 zus_y = zsin_wind(ji,jj) * zustar(ji,jj) / 0.3**2 342 zustke(ji,jj) = MAX ( SQRT( zus_x*zus_x + zus_y*zus_y), 1.0e-8 ) 343 ! dstokes(ji,jj) set to constant value rn_osm_dstokes from namelist in zdf_osm_init 344 END DO 345 END DO 339 DO_2D( 0, 0, 0, 0 ) 340 zus_x = zcos_wind(ji,jj) * zustar(ji,jj) / 0.3**2 341 zus_y = zsin_wind(ji,jj) * zustar(ji,jj) / 0.3**2 342 zustke(ji,jj) = MAX ( SQRT( zus_x*zus_x + zus_y*zus_y), 1.0e-8 ) 343 ! dstokes(ji,jj) set to constant value rn_osm_dstokes from namelist in zdf_osm_init 344 END_2D 346 345 ! Assume Pierson-Moskovitz wind-wave spectrum 347 346 CASE(1) 348 DO jj = 2, jpjm1 349 DO ji = 2, jpim1 350 ! Use wind speed wndm included in sbc_oce module 351 zustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8 ) 352 dstokes(ji,jj) = 0.12 * wndm(ji,jj)**2 / grav 353 END DO 354 END DO 347 DO_2D( 0, 0, 0, 0 ) 348 ! Use wind speed wndm included in sbc_oce module 349 zustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8 ) 350 dstokes(ji,jj) = 0.12 * wndm(ji,jj)**2 / grav 351 END_2D 355 352 ! Use ECMWF wave fields as output from SBCWAVE 356 353 CASE(2) 357 354 zfac = 2.0_wp * rpi / 16.0_wp 358 DO jj = 2, jpjm1 359 DO ji = 2, jpim1 360 ! The Langmur number from the ECMWF model appears to give La<0.3 for wind-driven seas. 361 ! The coefficient 0.8 gives La=0.3 in this situation. 362 ! It could represent the effects of the spread of wave directions 363 ! around the mean wind. The effect of this adjustment needs to be tested. 364 zustke(ji,jj) = MAX ( 1.0 * ( zcos_wind(ji,jj) * ut0sd(ji,jj ) + zsin_wind(ji,jj) * vt0sd(ji,jj) ), & 365 & zustar(ji,jj) / ( 0.45 * 0.45 ) ) 366 dstokes(ji,jj) = MAX(zfac * hsw(ji,jj)*hsw(ji,jj) / ( MAX(zustke(ji,jj)*wmp(ji,jj), 1.0e-7 ) ), 5.0e-1) !rn_osm_dstokes ! 367 END DO 368 END DO 355 DO_2D( 0, 0, 0, 0 ) 356 ! The Langmur number from the ECMWF model appears to give La<0.3 for wind-driven seas. 357 ! The coefficient 0.8 gives La=0.3 in this situation. 358 ! It could represent the effects of the spread of wave directions 359 ! around the mean wind. The effect of this adjustment needs to be tested. 360 zustke(ji,jj) = MAX ( 1.0 * ( zcos_wind(ji,jj) * ut0sd(ji,jj ) + zsin_wind(ji,jj) * vt0sd(ji,jj) ), & 361 & zustar(ji,jj) / ( 0.45 * 0.45 ) ) 362 dstokes(ji,jj) = MAX(zfac * hsw(ji,jj)*hsw(ji,jj) / ( MAX(zustke(ji,jj)*wmp(ji,jj), 1.0e-7 ) ), 5.0e-1) !rn_osm_dstokes ! 363 END_2D 369 364 END SELECT 370 365 371 366 ! Langmuir velocity scale (zwstrl), La # (zla) 372 367 ! mixed scale (zvstr), convective velocity scale (zwstrc) 373 DO jj = 2, jpjm1 374 DO ji = 2, jpim1 375 ! Langmuir velocity scale (zwstrl), at T-point 376 zwstrl(ji,jj) = ( zustar(ji,jj) * zustar(ji,jj) * zustke(ji,jj) )**pthird 377 ! Modify zwstrl to allow for small and large values of dstokes/hbl. 378 ! Intended as a possible test. Doesn't affect LES results for entrainment, 379 ! but hasn't been shown to be correct as dstokes/h becomes large or small. 380 zwstrl(ji,jj) = zwstrl(ji,jj) * & 381 & (1.12 * ( 1.0 - ( 1.0 - EXP( -hbl(ji,jj) / dstokes(ji,jj) ) ) * dstokes(ji,jj) / hbl(ji,jj) ))**pthird * & 382 & ( 1.0 - EXP( -15.0 * dstokes(ji,jj) / hbl(ji,jj) )) 383 ! define La this way so effects of Stokes penetration depth on velocity scale are included 384 zla(ji,jj) = SQRT ( zustar(ji,jj) / zwstrl(ji,jj) )**3 385 ! Velocity scale that tends to zustar for large Langmuir numbers 386 zvstr(ji,jj) = ( zwstrl(ji,jj)**3 + & 387 & ( 1.0 - EXP( -0.5 * zla(ji,jj)**2 ) ) * zustar(ji,jj) * zustar(ji,jj) * zustar(ji,jj) )**pthird 388 389 ! limit maximum value of Langmuir number as approximate treatment for shear turbulence. 390 ! Note zustke and zwstrl are not amended. 391 IF ( zla(ji,jj) >= 0.45 ) zla(ji,jj) = 0.45 392 ! 393 ! get convective velocity (zwstrc), stabilty scale (zhol) and logical conection flag lconv 394 IF ( zwbav(ji,jj) > 0.0) THEN 395 zwstrc(ji,jj) = ( 2.0 * zwbav(ji,jj) * 0.9 * hbl(ji,jj) )**pthird 396 zhol(ji,jj) = -0.9 * hbl(ji,jj) * 2.0 * zwbav(ji,jj) / (zvstr(ji,jj)**3 + epsln ) 397 lconv(ji,jj) = .TRUE. 398 ELSE 399 zhol(ji,jj) = -hbl(ji,jj) * 2.0 * zwbav(ji,jj)/ (zvstr(ji,jj)**3 + epsln ) 400 lconv(ji,jj) = .FALSE. 401 ENDIF 402 END DO 403 END DO 368 DO_2D( 0, 0, 0, 0 ) 369 ! Langmuir velocity scale (zwstrl), at T-point 370 zwstrl(ji,jj) = ( zustar(ji,jj) * zustar(ji,jj) * zustke(ji,jj) )**pthird 371 ! Modify zwstrl to allow for small and large values of dstokes/hbl. 372 ! Intended as a possible test. Doesn't affect LES results for entrainment, 373 ! but hasn't been shown to be correct as dstokes/h becomes large or small. 374 zwstrl(ji,jj) = zwstrl(ji,jj) * & 375 & (1.12 * ( 1.0 - ( 1.0 - EXP( -hbl(ji,jj) / dstokes(ji,jj) ) ) * dstokes(ji,jj) / hbl(ji,jj) ))**pthird * & 376 & ( 1.0 - EXP( -15.0 * dstokes(ji,jj) / hbl(ji,jj) )) 377 ! define La this way so effects of Stokes penetration depth on velocity scale are included 378 zla(ji,jj) = SQRT ( zustar(ji,jj) / zwstrl(ji,jj) )**3 379 ! Velocity scale that tends to zustar for large Langmuir numbers 380 zvstr(ji,jj) = ( zwstrl(ji,jj)**3 + & 381 & ( 1.0 - EXP( -0.5 * zla(ji,jj)**2 ) ) * zustar(ji,jj) * zustar(ji,jj) * zustar(ji,jj) )**pthird 382 383 ! limit maximum value of Langmuir number as approximate treatment for shear turbulence. 384 ! Note zustke and zwstrl are not amended. 385 IF ( zla(ji,jj) >= 0.45 ) zla(ji,jj) = 0.45 386 ! 387 ! get convective velocity (zwstrc), stabilty scale (zhol) and logical conection flag lconv 388 IF ( zwbav(ji,jj) > 0.0) THEN 389 zwstrc(ji,jj) = ( 2.0 * zwbav(ji,jj) * 0.9 * hbl(ji,jj) )**pthird 390 zhol(ji,jj) = -0.9 * hbl(ji,jj) * 2.0 * zwbav(ji,jj) / (zvstr(ji,jj)**3 + epsln ) 391 lconv(ji,jj) = .TRUE. 392 ELSE 393 zhol(ji,jj) = -hbl(ji,jj) * 2.0 * zwbav(ji,jj)/ (zvstr(ji,jj)**3 + epsln ) 394 lconv(ji,jj) = .FALSE. 395 ENDIF 396 END_2D 404 397 405 398 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 407 400 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 408 401 ! BL must be always 2 levels deep. 409 hbl(:,:) = MAX(hbl(:,:), gdepw _n(:,:,3) )402 hbl(:,:) = MAX(hbl(:,:), gdepw(:,:,3,Kmm) ) 410 403 ibld(:,:) = 3 411 DO jk = 4, jpkm1 412 DO jj = 2, jpjm1 413 DO ji = 2, jpim1 414 IF ( hbl(ji,jj) >= gdepw_n(ji,jj,jk) ) THEN 415 ibld(ji,jj) = MIN(mbkt(ji,jj), jk) 416 ENDIF 404 DO_3D( 0, 0, 0, 0, 4, jpkm1 ) 405 IF ( hbl(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN 406 ibld(ji,jj) = MIN(mbkt(ji,jj), jk) 407 ENDIF 408 END_3D 409 410 DO_2D( 0, 0, 0, 0 ) 411 zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 412 zbeta = rab_n(ji,jj,1,jp_sal) 413 zt = 0._wp 414 zs = 0._wp 415 zu = 0._wp 416 zv = 0._wp 417 ! average over depth of boundary layer 418 zthick=0._wp 419 DO jm = 2, ibld(ji,jj) 420 zthick=zthick+e3t(ji,jj,jm,Kmm) 421 zt = zt + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_tem,Kmm) 422 zs = zs + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_sal,Kmm) 423 zu = zu + e3t(ji,jj,jm,Kmm) & 424 & * ( uu(ji,jj,jm,Kbb) + uu(ji - 1,jj,jm,Kbb) ) & 425 & / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 426 zv = zv + e3t(ji,jj,jm,Kmm) & 427 & * ( vv(ji,jj,jm,Kbb) + vv(ji,jj - 1,jm,Kbb) ) & 428 & / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 417 429 END DO 418 END DO 419 END DO 420 421 DO jj = 2, jpjm1 ! Vertical slab 422 DO ji = 2, jpim1 423 zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 424 zbeta = rab_n(ji,jj,1,jp_sal) 425 zt = 0._wp 426 zs = 0._wp 427 zu = 0._wp 428 zv = 0._wp 429 ! average over depth of boundary layer 430 zthick=0._wp 431 DO jm = 2, ibld(ji,jj) 432 zthick=zthick+e3t_n(ji,jj,jm) 433 zt = zt + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem) 434 zs = zs + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal) 435 zu = zu + e3t_n(ji,jj,jm) & 436 & * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) & 437 & / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 438 zv = zv + e3t_n(ji,jj,jm) & 439 & * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) & 440 & / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 441 END DO 442 zt_bl(ji,jj) = zt / zthick 443 zs_bl(ji,jj) = zs / zthick 444 zu_bl(ji,jj) = zu / zthick 445 zv_bl(ji,jj) = zv / zthick 446 zdt_bl(ji,jj) = zt_bl(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_tem) 447 zds_bl(ji,jj) = zs_bl(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal) 448 zdu_bl(ji,jj) = zu_bl(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) & 449 & / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 450 zdv_bl(ji,jj) = zv_bl(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) & 451 & / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 452 zdb_bl(ji,jj) = grav * zthermal * zdt_bl(ji,jj) - grav * zbeta * zds_bl(ji,jj) 453 IF ( lconv(ji,jj) ) THEN ! Convective 454 zwb_ent(ji,jj) = -( 2.0 * 0.2 * zwbav(ji,jj) & 455 & + 0.135 * zla(ji,jj) * zwstrl(ji,jj)**3/hbl(ji,jj) ) 456 457 zvel_max = - ( 1.0 + 1.0 * ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / hbl(ji,jj) ) & 458 & * zwb_ent(ji,jj) / ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 430 zt_bl(ji,jj) = zt / zthick 431 zs_bl(ji,jj) = zs / zthick 432 zu_bl(ji,jj) = zu / zthick 433 zv_bl(ji,jj) = zv / zthick 434 zdt_bl(ji,jj) = zt_bl(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_tem,Kmm) 435 zds_bl(ji,jj) = zs_bl(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_sal,Kmm) 436 zdu_bl(ji,jj) = zu_bl(ji,jj) - ( uu(ji,jj,ibld(ji,jj),Kbb) + uu(ji-1,jj,ibld(ji,jj) ,Kbb) ) & 437 & / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 438 zdv_bl(ji,jj) = zv_bl(ji,jj) - ( vv(ji,jj,ibld(ji,jj),Kbb) + vv(ji,jj-1,ibld(ji,jj) ,Kbb) ) & 439 & / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 440 zdb_bl(ji,jj) = grav * zthermal * zdt_bl(ji,jj) - grav * zbeta * zds_bl(ji,jj) 441 IF ( lconv(ji,jj) ) THEN ! Convective 442 zwb_ent(ji,jj) = -( 2.0 * 0.2 * zwbav(ji,jj) & 443 & + 0.135 * zla(ji,jj) * zwstrl(ji,jj)**3/hbl(ji,jj) ) 444 445 zvel_max = - ( 1.0 + 1.0 * ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_Dt / hbl(ji,jj) ) & 446 & * zwb_ent(ji,jj) / ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 459 447 ! Entrainment including component due to shear turbulence. Modified Langmuir component, but gives same result for La=0.3 For testing uncomment. 460 448 ! zwb_ent(ji,jj) = -( 2.0 * 0.2 * zwbav(ji,jj) & 461 449 ! & + ( 0.15 * ( 1.0 - EXP( -0.5 * zla(ji,jj) ) ) + 0.03 / zla(ji,jj)**2 ) * zustar(ji,jj)**3/hbl(ji,jj) ) 462 450 463 ! zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_ rdt / zhbl(ji,jj) ) * zwb_ent(ji,jj) / &451 ! zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_Dt / zhbl(ji,jj) ) * zwb_ent(ji,jj) / & 464 452 ! & ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 465 zzdhdt = - zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj),0.0) ) 466 ELSE ! Stable 467 zzdhdt = 0.32 * ( hbli(ji,jj) / hbl(ji,jj) -1.0 ) * zwstrl(ji,jj)**3 / hbli(ji,jj) & 468 & + ( ( 0.32 / 3.0 ) * exp ( -2.5 * ( hbli(ji,jj) / hbl(ji,jj) - 1.0 ) ) & 469 & - ( 0.32 / 3.0 - 0.135 * zla(ji,jj) ) * exp ( -12.5 * ( hbli(ji,jj) / hbl(ji,jj) ) ) ) & 470 & * zwstrl(ji,jj)**3 / hbli(ji,jj) 471 zzdhdt = zzdhdt + zwbav(ji,jj) 472 IF ( zzdhdt < 0._wp ) THEN 473 ! For long timsteps factor in brackets slows the rapid collapse of the OSBL 474 zpert = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_rdt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj) 475 ELSE 476 zpert = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_rdt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj) & 477 & + MAX( zdb_bl(ji,jj), 0.0 ) 478 ENDIF 479 zzdhdt = 2.0 * zzdhdt / zpert 480 ENDIF 481 zdhdt(ji,jj) = zzdhdt 482 END DO 483 END DO 453 zzdhdt = - zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj),0.0) ) 454 ELSE ! Stable 455 zzdhdt = 0.32 * ( hbli(ji,jj) / hbl(ji,jj) -1.0 ) * zwstrl(ji,jj)**3 / hbli(ji,jj) & 456 & + ( ( 0.32 / 3.0 ) * exp ( -2.5 * ( hbli(ji,jj) / hbl(ji,jj) - 1.0 ) ) & 457 & - ( 0.32 / 3.0 - 0.135 * zla(ji,jj) ) * exp ( -12.5 * ( hbli(ji,jj) / hbl(ji,jj) ) ) ) & 458 & * zwstrl(ji,jj)**3 / hbli(ji,jj) 459 zzdhdt = zzdhdt + zwbav(ji,jj) 460 IF ( zzdhdt < 0._wp ) THEN 461 ! For long timsteps factor in brackets slows the rapid collapse of the OSBL 462 zpert = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_Dt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj) 463 ELSE 464 zpert = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_Dt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj) & 465 & + MAX( zdb_bl(ji,jj), 0.0 ) 466 ENDIF 467 zzdhdt = 2.0 * zzdhdt / zpert 468 ENDIF 469 zdhdt(ji,jj) = zzdhdt 470 END_2D 484 471 485 472 ! Calculate averages over depth of boundary layer … … 487 474 ibld(:,:) = 3 488 475 489 zhbl_t(:,:) = hbl(:,:) + (zdhdt(:,:) - wn(ji,jj,ibld(ji,jj)))* rn_rdt ! certainly need wb here, so subtract it 490 zhbl_t(:,:) = MIN(zhbl_t(:,:), ht_n(:,:)) 491 zdhdt(:,:) = MIN(zdhdt(:,:), (zhbl_t(:,:) - hbl(:,:))/rn_rdt + wn(ji,jj,ibld(ji,jj))) ! adjustment to represent limiting by ocean bottom 492 493 DO jk = 4, jpkm1 494 DO jj = 2, jpjm1 495 DO ji = 2, jpim1 496 IF ( zhbl_t(ji,jj) >= gdepw_n(ji,jj,jk) ) THEN 497 ibld(ji,jj) = MIN(mbkt(ji,jj), jk) 498 ENDIF 499 END DO 500 END DO 501 END DO 476 zhbl_t(:,:) = hbl(:,:) + (zdhdt(:,:) - ww(ji,jj,ibld(ji,jj)))* rn_Dt ! certainly need wb here, so subtract it 477 zhbl_t(:,:) = MIN(zhbl_t(:,:), ht(:,:)) 478 zdhdt(:,:) = MIN(zdhdt(:,:), (zhbl_t(:,:) - hbl(:,:))/rn_Dt + ww(ji,jj,ibld(ji,jj))) ! adjustment to represent limiting by ocean bottom 479 480 DO_3D( 0, 0, 0, 0, 4, jpkm1 ) 481 IF ( zhbl_t(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN 482 ibld(ji,jj) = MIN(mbkt(ji,jj), jk) 483 ENDIF 484 END_3D 502 485 503 486 ! 504 487 ! Step through model levels taking account of buoyancy change to determine the effect on dhdt 505 488 ! 506 DO jj = 2, jpjm1 507 DO ji = 2, jpim1 508 IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN 489 DO_2D( 0, 0, 0, 0 ) 490 IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN 509 491 ! 510 492 ! If boundary layer changes by more than one level, need to check for stable layers between initial and final depths. 511 493 ! 512 513 514 515 516 494 zhbl_s = hbl(ji,jj) 495 jm = imld(ji,jj) 496 zthermal = rab_n(ji,jj,1,jp_tem) 497 zbeta = rab_n(ji,jj,1,jp_sal) 498 IF ( lconv(ji,jj) ) THEN 517 499 !unstable 518 zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / hbl(ji,jj) ) & 519 & * zwb_ent(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 520 521 DO jk = imld(ji,jj), ibld(ji,jj) 522 zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) ) & 523 & - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ), 0.0 ) + zvel_max 524 525 zhbl_s = zhbl_s + MIN( - zwb_ent(ji,jj) / zdb * rn_rdt / FLOAT(ibld(ji,jj)-imld(ji,jj) ), e3w_n(ji,jj,jk) ) 526 zhbl_s = MIN(zhbl_s, ht_n(ji,jj)) 527 528 IF ( zhbl_s >= gdepw_n(ji,jj,jm+1) ) jm = jm + 1 529 END DO 530 hbl(ji,jj) = zhbl_s 531 ibld(ji,jj) = jm 532 hbli(ji,jj) = hbl(ji,jj) 533 ELSE 500 zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_Dt / hbl(ji,jj) ) & 501 & * zwb_ent(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 502 503 DO jk = imld(ji,jj), ibld(ji,jj) 504 zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - ts(ji,jj,jm,jp_tem,Kmm) ) & 505 & - zbeta * ( zs_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ), 0.0 ) + zvel_max 506 507 zhbl_s = zhbl_s + MIN( - zwb_ent(ji,jj) / zdb * rn_Dt / FLOAT(ibld(ji,jj)-imld(ji,jj) ), & 508 & e3w(ji,jj,jk,Kmm) ) 509 zhbl_s = MIN(zhbl_s, ht(ji,jj)) 510 511 IF ( zhbl_s >= gdepw(ji,jj,jm+1,Kmm) ) jm = jm + 1 512 END DO 513 hbl(ji,jj) = zhbl_s 514 ibld(ji,jj) = jm 515 hbli(ji,jj) = hbl(ji,jj) 516 ELSE 534 517 ! stable 535 DO jk = imld(ji,jj), ibld(ji,jj) 536 zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) ) & 537 & - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ), 0.0 ) & 538 & + 2.0 * zwstrl(ji,jj)**2 / zhbl_s 539 540 zhbl_s = zhbl_s + ( & 541 & 0.32 * ( hbli(ji,jj) / zhbl_s -1.0 ) & 542 & * zwstrl(ji,jj)**3 / hbli(ji,jj) & 543 & + ( ( 0.32 / 3.0 ) * EXP( - 2.5 * ( hbli(ji,jj) / zhbl_s -1.0 ) ) & 544 & - ( 0.32 / 3.0 - 0.0485 ) * EXP( - 12.5 * ( hbli(ji,jj) / zhbl_s ) ) ) & 545 & * zwstrl(ji,jj)**3 / hbli(ji,jj) ) / zdb * e3w_n(ji,jj,jk) / zdhdt(ji,jj) ! ALMG to investigate whether need to include wn here 546 547 zhbl_s = MIN(zhbl_s, ht_n(ji,jj)) 548 IF ( zhbl_s >= gdepw_n(ji,jj,jm) ) jm = jm + 1 549 END DO 550 hbl(ji,jj) = MAX(zhbl_s, gdepw_n(ji,jj,3) ) 551 ibld(ji,jj) = MAX(jm, 3 ) 552 IF ( hbl(ji,jj) > hbli(ji,jj) ) hbli(ji,jj) = hbl(ji,jj) 553 ENDIF ! IF ( lconv ) 518 DO jk = imld(ji,jj), ibld(ji,jj) 519 zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - ts(ji,jj,jm,jp_tem,Kmm) ) & 520 & - zbeta * ( zs_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ), 0.0 ) & 521 & + 2.0 * zwstrl(ji,jj)**2 / zhbl_s 522 523 zhbl_s = zhbl_s + ( & 524 & 0.32 * ( hbli(ji,jj) / zhbl_s -1.0 ) & 525 & * zwstrl(ji,jj)**3 / hbli(ji,jj) & 526 & + ( ( 0.32 / 3.0 ) * EXP( - 2.5 * ( hbli(ji,jj) / zhbl_s -1.0 ) ) & 527 & - ( 0.32 / 3.0 - 0.0485 ) * EXP( - 12.5 * ( hbli(ji,jj) / zhbl_s ) ) ) & 528 & * zwstrl(ji,jj)**3 / hbli(ji,jj) ) / zdb * e3w(ji,jj,jk,Kmm) / zdhdt(ji,jj) ! ALMG to investigate whether need to include ww here 529 530 zhbl_s = MIN(zhbl_s, ht(ji,jj)) 531 IF ( zhbl_s >= gdepw(ji,jj,jm,Kmm) ) jm = jm + 1 532 END DO 533 hbl(ji,jj) = MAX(zhbl_s, gdepw(ji,jj,3,Kmm) ) 534 ibld(ji,jj) = MAX(jm, 3 ) 535 IF ( hbl(ji,jj) > hbli(ji,jj) ) hbli(ji,jj) = hbl(ji,jj) 536 ENDIF ! IF ( lconv ) 537 ELSE 538 ! change zero or one model level. 539 hbl(ji,jj) = zhbl_t(ji,jj) 540 IF ( lconv(ji,jj) ) THEN 541 hbli(ji,jj) = hbl(ji,jj) 554 542 ELSE 555 ! change zero or one model level. 556 hbl(ji,jj) = zhbl_t(ji,jj) 557 IF ( lconv(ji,jj) ) THEN 558 hbli(ji,jj) = hbl(ji,jj) 559 ELSE 560 hbl(ji,jj) = MAX(hbl(ji,jj), gdepw_n(ji,jj,3) ) 561 IF ( hbl(ji,jj) > hbli(ji,jj) ) hbli(ji,jj) = hbl(ji,jj) 562 ENDIF 543 hbl(ji,jj) = MAX(hbl(ji,jj), gdepw(ji,jj,3,Kmm) ) 544 IF ( hbl(ji,jj) > hbli(ji,jj) ) hbli(ji,jj) = hbl(ji,jj) 563 545 ENDIF 564 zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj))565 END DO566 END DO546 ENDIF 547 zhbl(ji,jj) = gdepw(ji,jj,ibld(ji,jj),Kmm) 548 END_2D 567 549 dstokes(:,:) = MIN ( dstokes(:,:), hbl(:,:)/3. ) ! Limit delta for shallow boundary layers for calculating flux-gradient terms. 568 550 … … 570 552 ! Consider later combining this into the loop above and looking for columns 571 553 ! where the index for base of the boundary layer have changed 572 DO jj = 2, jpjm1 ! Vertical slab 573 DO ji = 2, jpim1 574 zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 575 zbeta = rab_n(ji,jj,1,jp_sal) 576 zt = 0._wp 577 zs = 0._wp 578 zu = 0._wp 579 zv = 0._wp 580 ! average over depth of boundary layer 581 zthick=0._wp 582 DO jm = 2, ibld(ji,jj) 583 zthick=zthick+e3t_n(ji,jj,jm) 584 zt = zt + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem) 585 zs = zs + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal) 586 zu = zu + e3t_n(ji,jj,jm) & 587 & * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) & 588 & / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 589 zv = zv + e3t_n(ji,jj,jm) & 590 & * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) & 591 & / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 592 END DO 593 zt_bl(ji,jj) = zt / zthick 594 zs_bl(ji,jj) = zs / zthick 595 zu_bl(ji,jj) = zu / zthick 596 zv_bl(ji,jj) = zv / zthick 597 zdt_bl(ji,jj) = zt_bl(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_tem) 598 zds_bl(ji,jj) = zs_bl(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal) 599 zdu_bl(ji,jj) = zu_bl(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) & 600 & / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 601 zdv_bl(ji,jj) = zv_bl(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) & 602 & / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 603 zdb_bl(ji,jj) = grav * zthermal * zdt_bl(ji,jj) - grav * zbeta * zds_bl(ji,jj) 604 zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj)) 605 IF ( lconv(ji,jj) ) THEN 606 IF ( zdb_bl(ji,jj) > 0._wp )THEN 607 IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN ! near neutral stability 608 zari = 4.5 * ( zvstr(ji,jj)**2 ) & 609 & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01 610 ELSE ! unstable 611 zari = 4.5 * ( zwstrc(ji,jj)**2 ) & 612 & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01 613 ENDIF 614 IF ( zari > 0.2 ) THEN ! This test checks for weakly stratified pycnocline 615 zari = 0.2 616 zwb_ent(ji,jj) = 0._wp 617 ENDIF 618 inhml = MAX( INT( zari * zhbl(ji,jj) / e3t_n(ji,jj,ibld(ji,jj)) ) , 1 ) 554 DO_2D( 0, 0, 0, 0 ) 555 zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 556 zbeta = rab_n(ji,jj,1,jp_sal) 557 zt = 0._wp 558 zs = 0._wp 559 zu = 0._wp 560 zv = 0._wp 561 ! average over depth of boundary layer 562 zthick=0._wp 563 DO jm = 2, ibld(ji,jj) 564 zthick=zthick+e3t(ji,jj,jm,Kmm) 565 zt = zt + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_tem,Kmm) 566 zs = zs + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_sal,Kmm) 567 zu = zu + e3t(ji,jj,jm,Kmm) & 568 & * ( uu(ji,jj,jm,Kbb) + uu(ji - 1,jj,jm,Kbb) ) & 569 & / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 570 zv = zv + e3t(ji,jj,jm,Kmm) & 571 & * ( vv(ji,jj,jm,Kbb) + vv(ji,jj - 1,jm,Kbb) ) & 572 & / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 573 END DO 574 zt_bl(ji,jj) = zt / zthick 575 zs_bl(ji,jj) = zs / zthick 576 zu_bl(ji,jj) = zu / zthick 577 zv_bl(ji,jj) = zv / zthick 578 zdt_bl(ji,jj) = zt_bl(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_tem,Kmm) 579 zds_bl(ji,jj) = zs_bl(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_sal,Kmm) 580 zdu_bl(ji,jj) = zu_bl(ji,jj) - ( uu(ji,jj,ibld(ji,jj),Kbb) + uu(ji-1,jj,ibld(ji,jj) ,Kbb) ) & 581 & / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 582 zdv_bl(ji,jj) = zv_bl(ji,jj) - ( vv(ji,jj,ibld(ji,jj),Kbb) + vv(ji,jj-1,ibld(ji,jj) ,Kbb) ) & 583 & / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 584 zdb_bl(ji,jj) = grav * zthermal * zdt_bl(ji,jj) - grav * zbeta * zds_bl(ji,jj) 585 zhbl(ji,jj) = gdepw(ji,jj,ibld(ji,jj),Kmm) 586 IF ( lconv(ji,jj) ) THEN 587 IF ( zdb_bl(ji,jj) > 0._wp )THEN 588 IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN ! near neutral stability 589 zari = 4.5 * ( zvstr(ji,jj)**2 ) & 590 & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01 591 ELSE ! unstable 592 zari = 4.5 * ( zwstrc(ji,jj)**2 ) & 593 & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01 594 ENDIF 595 IF ( zari > 0.2 ) THEN ! This test checks for weakly stratified pycnocline 596 zari = 0.2 597 zwb_ent(ji,jj) = 0._wp 598 ENDIF 599 inhml = MAX( INT( zari * zhbl(ji,jj) & 600 & / e3t(ji,jj,ibld(ji,jj),Kmm) ), 1 ) 601 imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 1) 602 zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm) 603 zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 604 ELSE ! IF (zdb_bl) 605 imld(ji,jj) = ibld(ji,jj) - 1 606 zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm) 607 zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 608 ENDIF 609 ELSE ! IF (lconv) 610 IF ( zdhdt(ji,jj) >= 0.0 ) THEN ! probably shouldn't include wm here 611 ! boundary layer deepening 612 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 613 ! pycnocline thickness set by stratification - use same relationship as for neutral conditions. 614 zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) & 615 & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01 , 0.2 ) 616 inhml = MAX( INT( zari * zhbl(ji,jj) & 617 & / e3t(ji,jj,ibld(ji,jj),Kmm) ), 1 ) 619 618 imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 1) 620 zhml(ji,jj) = gdepw _n(ji,jj,imld(ji,jj))619 zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm) 621 620 zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 622 ELSE ! IF (zdb_bl)621 ELSE 623 622 imld(ji,jj) = ibld(ji,jj) - 1 624 zhml(ji,jj) = gdepw _n(ji,jj,imld(ji,jj))623 zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm) 625 624 zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 626 ENDIF 627 ELSE ! IF (lconv) 628 IF ( zdhdt(ji,jj) >= 0.0 ) THEN ! probably shouldn't include wm here 629 ! boundary layer deepening 630 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 631 ! pycnocline thickness set by stratification - use same relationship as for neutral conditions. 632 zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) & 633 & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01 , 0.2 ) 634 inhml = MAX( INT( zari * zhbl(ji,jj) / e3t_n(ji,jj,ibld(ji,jj)) ) , 1 ) 635 imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 1) 636 zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 637 zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 638 ELSE 639 imld(ji,jj) = ibld(ji,jj) - 1 640 zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 641 zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 642 ENDIF ! IF (zdb_bl > 0.0) 643 ELSE ! IF(dhdt >= 0) 644 ! boundary layer collapsing. 645 imld(ji,jj) = ibld(ji,jj) 646 zhml(ji,jj) = zhbl(ji,jj) 647 zdh(ji,jj) = 0._wp 648 ENDIF ! IF (dhdt >= 0) 649 ENDIF ! IF (lconv) 650 END DO 651 END DO 625 ENDIF ! IF (zdb_bl > 0.0) 626 ELSE ! IF(dhdt >= 0) 627 ! boundary layer collapsing. 628 imld(ji,jj) = ibld(ji,jj) 629 zhml(ji,jj) = zhbl(ji,jj) 630 zdh(ji,jj) = 0._wp 631 ENDIF ! IF (dhdt >= 0) 632 ENDIF ! IF (lconv) 633 END_2D 652 634 653 635 ! Average over the depth of the mixed layer in the convective boundary layer 654 636 ! Also calculate entrainment fluxes for temperature and salinity 655 DO jj = 2, jpjm1 ! Vertical slab 656 DO ji = 2, jpim1 657 zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 658 zbeta = rab_n(ji,jj,1,jp_sal) 659 IF ( lconv(ji,jj) ) THEN 637 DO_2D( 0, 0, 0, 0 ) 638 zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 639 zbeta = rab_n(ji,jj,1,jp_sal) 640 IF ( lconv(ji,jj) ) THEN 641 zt = 0._wp 642 zs = 0._wp 643 zu = 0._wp 644 zv = 0._wp 645 ! average over depth of boundary layer 646 zthick=0._wp 647 DO jm = 2, imld(ji,jj) 648 zthick=zthick+e3t(ji,jj,jm,Kmm) 649 zt = zt + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_tem,Kmm) 650 zs = zs + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_sal,Kmm) 651 zu = zu + e3t(ji,jj,jm,Kmm) & 652 & * ( uu(ji,jj,jm,Kbb) + uu(ji - 1,jj,jm,Kbb) ) & 653 & / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 654 zv = zv + e3t(ji,jj,jm,Kmm) & 655 & * ( vv(ji,jj,jm,Kbb) + vv(ji,jj - 1,jm,Kbb) ) & 656 & / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 657 END DO 658 zt_ml(ji,jj) = zt / zthick 659 zs_ml(ji,jj) = zs / zthick 660 zu_ml(ji,jj) = zu / zthick 661 zv_ml(ji,jj) = zv / zthick 662 zdt_ml(ji,jj) = zt_ml(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_tem,Kmm) 663 zds_ml(ji,jj) = zs_ml(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_sal,Kmm) 664 zdu_ml(ji,jj) = zu_ml(ji,jj) - ( uu(ji,jj,ibld(ji,jj),Kbb) + uu(ji-1,jj,ibld(ji,jj) ,Kbb) ) & 665 & / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 666 zdv_ml(ji,jj) = zv_ml(ji,jj) - ( vv(ji,jj,ibld(ji,jj),Kbb) + vv(ji,jj-1,ibld(ji,jj) ,Kbb) ) & 667 & / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 668 zdb_ml(ji,jj) = grav * zthermal * zdt_ml(ji,jj) - grav * zbeta * zds_ml(ji,jj) 669 ELSE 670 ! stable, if entraining calulate average below interface layer. 671 IF ( zdhdt(ji,jj) >= 0._wp ) THEN 660 672 zt = 0._wp 661 673 zs = 0._wp … … 665 677 zthick=0._wp 666 678 DO jm = 2, imld(ji,jj) 667 zthick=zthick+e3t _n(ji,jj,jm)668 zt = zt + e3t _n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem)669 zs = zs + e3t _n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal)670 zu = zu + e3t _n(ji,jj,jm) &671 & * ( u b(ji,jj,jm) + ub(ji - 1,jj,jm) ) &679 zthick=zthick+e3t(ji,jj,jm,Kmm) 680 zt = zt + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_tem,Kmm) 681 zs = zs + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_sal,Kmm) 682 zu = zu + e3t(ji,jj,jm,Kmm) & 683 & * ( uu(ji,jj,jm,Kbb) + uu(ji - 1,jj,jm,Kbb) ) & 672 684 & / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 673 zv = zv + e3t _n(ji,jj,jm) &674 & * ( v b(ji,jj,jm) + vb(ji,jj - 1,jm) ) &685 zv = zv + e3t(ji,jj,jm,Kmm) & 686 & * ( vv(ji,jj,jm,Kbb) + vv(ji,jj - 1,jm,Kbb) ) & 675 687 & / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 676 688 END DO … … 679 691 zu_ml(ji,jj) = zu / zthick 680 692 zv_ml(ji,jj) = zv / zthick 681 zdt_ml(ji,jj) = zt_ml(ji,jj) - ts n(ji,jj,ibld(ji,jj),jp_tem)682 zds_ml(ji,jj) = zs_ml(ji,jj) - ts n(ji,jj,ibld(ji,jj),jp_sal)683 zdu_ml(ji,jj) = zu_ml(ji,jj) - ( u b(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj)) ) &693 zdt_ml(ji,jj) = zt_ml(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_tem,Kmm) 694 zds_ml(ji,jj) = zs_ml(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_sal,Kmm) 695 zdu_ml(ji,jj) = zu_ml(ji,jj) - ( uu(ji,jj,ibld(ji,jj),Kbb) + uu(ji-1,jj,ibld(ji,jj) ,Kbb) ) & 684 696 & / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 685 zdv_ml(ji,jj) = zv_ml(ji,jj) - ( v b(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj)) ) &697 zdv_ml(ji,jj) = zv_ml(ji,jj) - ( vv(ji,jj,ibld(ji,jj),Kbb) + vv(ji,jj-1,ibld(ji,jj) ,Kbb) ) & 686 698 & / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 687 699 zdb_ml(ji,jj) = grav * zthermal * zdt_ml(ji,jj) - grav * zbeta * zds_ml(ji,jj) 688 ELSE689 ! stable, if entraining calulate average below interface layer.690 IF ( zdhdt(ji,jj) >= 0._wp ) THEN691 zt = 0._wp692 zs = 0._wp693 zu = 0._wp694 zv = 0._wp695 ! average over depth of boundary layer696 zthick=0._wp697 DO jm = 2, imld(ji,jj)698 zthick=zthick+e3t_n(ji,jj,jm)699 zt = zt + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem)700 zs = zs + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal)701 zu = zu + e3t_n(ji,jj,jm) &702 & * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) &703 & / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) )704 zv = zv + e3t_n(ji,jj,jm) &705 & * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) &706 & / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) )707 END DO708 zt_ml(ji,jj) = zt / zthick709 zs_ml(ji,jj) = zs / zthick710 zu_ml(ji,jj) = zu / zthick711 zv_ml(ji,jj) = zv / zthick712 zdt_ml(ji,jj) = zt_ml(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_tem)713 zds_ml(ji,jj) = zs_ml(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal)714 zdu_ml(ji,jj) = zu_ml(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) &715 & / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) )716 zdv_ml(ji,jj) = zv_ml(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) &717 & / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) )718 zdb_ml(ji,jj) = grav * zthermal * zdt_ml(ji,jj) - grav * zbeta * zds_ml(ji,jj)719 ENDIF720 700 ENDIF 721 END DO722 END DO701 ENDIF 702 END_2D 723 703 ! 724 704 ! rotate mean currents and changes onto wind align co-ordinates 725 705 ! 726 706 727 DO jj = 2, jpjm1 728 DO ji = 2, jpim1 729 ztemp = zu_ml(ji,jj) 730 zu_ml(ji,jj) = zu_ml(ji,jj) * zcos_wind(ji,jj) + zv_ml(ji,jj) * zsin_wind(ji,jj) 731 zv_ml(ji,jj) = zv_ml(ji,jj) * zcos_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 732 ztemp = zdu_ml(ji,jj) 733 zdu_ml(ji,jj) = zdu_ml(ji,jj) * zcos_wind(ji,jj) + zdv_ml(ji,jj) * zsin_wind(ji,jj) 734 zdv_ml(ji,jj) = zdv_ml(ji,jj) * zsin_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 735 ! 736 ztemp = zu_bl(ji,jj) 737 zu_bl = zu_bl(ji,jj) * zcos_wind(ji,jj) + zv_bl(ji,jj) * zsin_wind(ji,jj) 738 zv_bl(ji,jj) = zv_bl(ji,jj) * zcos_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 739 ztemp = zdu_bl(ji,jj) 740 zdu_bl(ji,jj) = zdu_bl(ji,jj) * zcos_wind(ji,jj) + zdv_bl(ji,jj) * zsin_wind(ji,jj) 741 zdv_bl(ji,jj) = zdv_bl(ji,jj) * zsin_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 742 END DO 743 END DO 707 DO_2D( 0, 0, 0, 0 ) 708 ztemp = zu_ml(ji,jj) 709 zu_ml(ji,jj) = zu_ml(ji,jj) * zcos_wind(ji,jj) + zv_ml(ji,jj) * zsin_wind(ji,jj) 710 zv_ml(ji,jj) = zv_ml(ji,jj) * zcos_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 711 ztemp = zdu_ml(ji,jj) 712 zdu_ml(ji,jj) = zdu_ml(ji,jj) * zcos_wind(ji,jj) + zdv_ml(ji,jj) * zsin_wind(ji,jj) 713 zdv_ml(ji,jj) = zdv_ml(ji,jj) * zsin_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 714 ! 715 ztemp = zu_bl(ji,jj) 716 zu_bl = zu_bl(ji,jj) * zcos_wind(ji,jj) + zv_bl(ji,jj) * zsin_wind(ji,jj) 717 zv_bl(ji,jj) = zv_bl(ji,jj) * zcos_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 718 ztemp = zdu_bl(ji,jj) 719 zdu_bl(ji,jj) = zdu_bl(ji,jj) * zcos_wind(ji,jj) + zdv_bl(ji,jj) * zsin_wind(ji,jj) 720 zdv_bl(ji,jj) = zdv_bl(ji,jj) * zsin_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 721 END_2D 744 722 745 723 zuw_bse = 0._wp 746 724 zvw_bse = 0._wp 747 DO jj = 2, jpjm1 748 DO ji = 2, jpim1 749 750 IF ( lconv(ji,jj) ) THEN 751 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 752 zwth_ent(ji,jj) = zwb_ent(ji,jj) * zdt_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 753 zws_ent(ji,jj) = zwb_ent(ji,jj) * zds_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 754 ENDIF 755 ELSE 756 zwth_ent(ji,jj) = -2.0 * zwthav(ji,jj) * ( (1.0 - 0.8) - ( 1.0 - 0.8)**(3.0/2.0) ) 757 zws_ent(ji,jj) = -2.0 * zwsav(ji,jj) * ( (1.0 - 0.8 ) - ( 1.0 - 0.8 )**(3.0/2.0) ) 725 DO_2D( 0, 0, 0, 0 ) 726 727 IF ( lconv(ji,jj) ) THEN 728 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 729 zwth_ent(ji,jj) = zwb_ent(ji,jj) * zdt_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 730 zws_ent(ji,jj) = zwb_ent(ji,jj) * zds_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 758 731 ENDIF 759 END DO 760 END DO 732 ELSE 733 zwth_ent(ji,jj) = -2.0 * zwthav(ji,jj) * ( (1.0 - 0.8) - ( 1.0 - 0.8)**(3.0/2.0) ) 734 zws_ent(ji,jj) = -2.0 * zwsav(ji,jj) * ( (1.0 - 0.8 ) - ( 1.0 - 0.8 )**(3.0/2.0) ) 735 ENDIF 736 END_2D 761 737 762 738 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 764 740 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 765 741 766 DO jj = 2, jpjm1 767 DO ji = 2, jpim1 768 ! 769 IF ( lconv (ji,jj) ) THEN 770 ! Unstable conditions 771 IF( zdb_bl(ji,jj) > 0._wp ) THEN 772 ! calculate pycnocline profiles, no need if zdb_bl <= 0. since profile is zero and arrays have been initialized to zero 773 ztgrad = ( zdt_ml(ji,jj) / zdh(ji,jj) ) 774 zsgrad = ( zds_ml(ji,jj) / zdh(ji,jj) ) 775 zbgrad = ( zdb_ml(ji,jj) / zdh(ji,jj) ) 776 DO jk = 2 , ibld(ji,jj) 777 znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 778 zdtdz_pyc(ji,jj,jk) = ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 779 zdbdz_pyc(ji,jj,jk) = zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 780 zdsdz_pyc(ji,jj,jk) = zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 781 END DO 782 ENDIF 783 ELSE 784 ! stable conditions 785 ! if pycnocline profile only defined when depth steady of increasing. 786 IF ( zdhdt(ji,jj) >= 0.0 ) THEN ! Depth increasing, or steady. 787 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 788 IF ( zhol(ji,jj) >= 0.5 ) THEN ! Very stable - 'thick' pycnocline 789 ztgrad = zdt_bl(ji,jj) / zhbl(ji,jj) 790 zsgrad = zds_bl(ji,jj) / zhbl(ji,jj) 791 zbgrad = zdb_bl(ji,jj) / zhbl(ji,jj) 792 DO jk = 2, ibld(ji,jj) 793 znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 794 zdtdz_pyc(ji,jj,jk) = ztgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 795 zdbdz_pyc(ji,jj,jk) = zbgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 796 zdsdz_pyc(ji,jj,jk) = zsgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 797 END DO 798 ELSE ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form. 799 ztgrad = zdt_bl(ji,jj) / zdh(ji,jj) 800 zsgrad = zds_bl(ji,jj) / zdh(ji,jj) 801 zbgrad = zdb_bl(ji,jj) / zdh(ji,jj) 802 DO jk = 2, ibld(ji,jj) 803 znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 804 zdtdz_pyc(ji,jj,jk) = ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 805 zdbdz_pyc(ji,jj,jk) = zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 806 zdsdz_pyc(ji,jj,jk) = zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 807 END DO 808 ENDIF ! IF (zhol >=0.5) 809 ENDIF ! IF (zdb_bl> 0.) 810 ENDIF ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero, profile arrays are intialized to zero 811 ENDIF ! IF (lconv) 812 ! 813 END DO 814 END DO 815 ! 816 DO jj = 2, jpjm1 817 DO ji = 2, jpim1 818 ! 819 IF ( lconv (ji,jj) ) THEN 820 ! Unstable conditions 821 zugrad = ( zdu_ml(ji,jj) / zdh(ji,jj) ) + 0.275 * zustar(ji,jj)*zustar(ji,jj) / & 822 & (( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) / zla(ji,jj)**(8.0/3.0) 823 zvgrad = ( zdv_ml(ji,jj) / zdh(ji,jj) ) + 3.5 * ff_t(ji,jj) * zustke(ji,jj) / & 824 & ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 825 DO jk = 2 , ibld(ji,jj)-1 826 znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 827 zdudz_pyc(ji,jj,jk) = zugrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 828 zdvdz_pyc(ji,jj,jk) = zvgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 829 END DO 830 ELSE 831 ! stable conditions 832 zugrad = 3.25 * zdu_bl(ji,jj) / zhbl(ji,jj) 833 zvgrad = 2.75 * zdv_bl(ji,jj) / zhbl(ji,jj) 834 DO jk = 2, ibld(ji,jj) 835 znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 836 IF ( znd < 1.0 ) THEN 837 zdudz_pyc(ji,jj,jk) = zugrad * EXP( -40.0 * ( znd - 1.0 )**2 ) 838 ELSE 839 zdudz_pyc(ji,jj,jk) = zugrad * EXP( -20.0 * ( znd - 1.0 )**2 ) 840 ENDIF 841 zdvdz_pyc(ji,jj,jk) = zvgrad * EXP( -20.0 * ( znd - 0.85 )**2 ) 742 DO_2D( 0, 0, 0, 0 ) 743 ! 744 IF ( lconv (ji,jj) ) THEN 745 ! Unstable conditions 746 IF( zdb_bl(ji,jj) > 0._wp ) THEN 747 ! calculate pycnocline profiles, no need if zdb_bl <= 0. since profile is zero and arrays have been initialized to zero 748 ztgrad = ( zdt_ml(ji,jj) / zdh(ji,jj) ) 749 zsgrad = ( zds_ml(ji,jj) / zdh(ji,jj) ) 750 zbgrad = ( zdb_ml(ji,jj) / zdh(ji,jj) ) 751 DO jk = 2 , ibld(ji,jj) 752 znd = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj) 753 zdtdz_pyc(ji,jj,jk) = ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 754 zdbdz_pyc(ji,jj,jk) = zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 755 zdsdz_pyc(ji,jj,jk) = zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 842 756 END DO 843 757 ENDIF 844 ! 845 END DO 846 END DO 758 ELSE 759 ! stable conditions 760 ! if pycnocline profile only defined when depth steady of increasing. 761 IF ( zdhdt(ji,jj) >= 0.0 ) THEN ! Depth increasing, or steady. 762 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 763 IF ( zhol(ji,jj) >= 0.5 ) THEN ! Very stable - 'thick' pycnocline 764 ztgrad = zdt_bl(ji,jj) / zhbl(ji,jj) 765 zsgrad = zds_bl(ji,jj) / zhbl(ji,jj) 766 zbgrad = zdb_bl(ji,jj) / zhbl(ji,jj) 767 DO jk = 2, ibld(ji,jj) 768 znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 769 zdtdz_pyc(ji,jj,jk) = ztgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 770 zdbdz_pyc(ji,jj,jk) = zbgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 771 zdsdz_pyc(ji,jj,jk) = zsgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 772 END DO 773 ELSE ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form. 774 ztgrad = zdt_bl(ji,jj) / zdh(ji,jj) 775 zsgrad = zds_bl(ji,jj) / zdh(ji,jj) 776 zbgrad = zdb_bl(ji,jj) / zdh(ji,jj) 777 DO jk = 2, ibld(ji,jj) 778 znd = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj) 779 zdtdz_pyc(ji,jj,jk) = ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 780 zdbdz_pyc(ji,jj,jk) = zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 781 zdsdz_pyc(ji,jj,jk) = zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 782 END DO 783 ENDIF ! IF (zhol >=0.5) 784 ENDIF ! IF (zdb_bl> 0.) 785 ENDIF ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero, profile arrays are intialized to zero 786 ENDIF ! IF (lconv) 787 ! 788 END_2D 789 ! 790 DO_2D( 0, 0, 0, 0 ) 791 ! 792 IF ( lconv (ji,jj) ) THEN 793 ! Unstable conditions 794 zugrad = ( zdu_ml(ji,jj) / zdh(ji,jj) ) + 0.275 * zustar(ji,jj)*zustar(ji,jj) / & 795 & (( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) / zla(ji,jj)**(8.0/3.0) 796 zvgrad = ( zdv_ml(ji,jj) / zdh(ji,jj) ) + 3.5 * ff_t(ji,jj) * zustke(ji,jj) / & 797 & ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 798 DO jk = 2 , ibld(ji,jj)-1 799 znd = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj) 800 zdudz_pyc(ji,jj,jk) = zugrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 801 zdvdz_pyc(ji,jj,jk) = zvgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 802 END DO 803 ELSE 804 ! stable conditions 805 zugrad = 3.25 * zdu_bl(ji,jj) / zhbl(ji,jj) 806 zvgrad = 2.75 * zdv_bl(ji,jj) / zhbl(ji,jj) 807 DO jk = 2, ibld(ji,jj) 808 znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 809 IF ( znd < 1.0 ) THEN 810 zdudz_pyc(ji,jj,jk) = zugrad * EXP( -40.0 * ( znd - 1.0 )**2 ) 811 ELSE 812 zdudz_pyc(ji,jj,jk) = zugrad * EXP( -20.0 * ( znd - 1.0 )**2 ) 813 ENDIF 814 zdvdz_pyc(ji,jj,jk) = zvgrad * EXP( -20.0 * ( znd - 0.85 )**2 ) 815 END DO 816 ENDIF 817 ! 818 END_2D 847 819 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 848 820 ! Eddy viscosity/diffusivity and non-gradient terms in the flux-gradient relationship … … 860 832 ! zvisml_sc = zwstrl * zhbl * EXP ( -( zhol / 0.183_wp )**2 ) 861 833 ! ENDWHERE 862 DO jj = 2, jpjm1 863 DO ji = 2, jpim1 864 IF ( lconv(ji,jj) ) THEN 865 zdifml_sc(ji,jj) = zhml(ji,jj) * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 866 zvisml_sc(ji,jj) = zdifml_sc(ji,jj) 867 zdifpyc_sc(ji,jj) = 0.165 * ( zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3 )**pthird * zdh(ji,jj) 868 zvispyc_sc(ji,jj) = 0.142 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zdh(ji,jj) 869 zbeta_d_sc(ji,jj) = 1.0 - (0.165 / 0.8 * zdh(ji,jj) / zhbl(ji,jj) )**p2third 870 zbeta_v_sc(ji,jj) = 1.0 - 2.0 * (0.142 /0.375) * zdh(ji,jj) / zhml(ji,jj) 871 ELSE 872 zdifml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ) 873 zvisml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ) 874 END IF 875 END DO 876 END DO 834 DO_2D( 0, 0, 0, 0 ) 835 IF ( lconv(ji,jj) ) THEN 836 zdifml_sc(ji,jj) = zhml(ji,jj) * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 837 zvisml_sc(ji,jj) = zdifml_sc(ji,jj) 838 zdifpyc_sc(ji,jj) = 0.165 * ( zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3 )**pthird * zdh(ji,jj) 839 zvispyc_sc(ji,jj) = 0.142 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zdh(ji,jj) 840 zbeta_d_sc(ji,jj) = 1.0 - (0.165 / 0.8 * zdh(ji,jj) / zhbl(ji,jj) )**p2third 841 zbeta_v_sc(ji,jj) = 1.0 - 2.0 * (0.142 /0.375) * zdh(ji,jj) / zhml(ji,jj) 842 ELSE 843 zdifml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ) 844 zvisml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ) 845 END IF 846 END_2D 877 847 ! 878 DO jj = 2, jpjm1 879 DO ji = 2, jpim1 880 IF ( lconv(ji,jj) ) THEN 881 DO jk = 2, imld(ji,jj) ! mixed layer diffusivity 882 zznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj) 848 DO_2D( 0, 0, 0, 0 ) 849 IF ( lconv(ji,jj) ) THEN 850 DO jk = 2, imld(ji,jj) ! mixed layer diffusivity 851 zznd_ml = gdepw(ji,jj,jk,Kmm) / zhml(ji,jj) 852 ! 853 zdiffut(ji,jj,jk) = 0.8 * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_d_sc(ji,jj) * zznd_ml )**1.5 854 ! 855 zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_v_sc(ji,jj) * zznd_ml ) & 856 & * ( 1.0 - 0.5 * zznd_ml**2 ) 857 END DO 858 ! pycnocline - if present linear profile 859 IF ( zdh(ji,jj) > 0._wp ) THEN 860 DO jk = imld(ji,jj)+1 , ibld(ji,jj) 861 zznd_pyc = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj) 883 862 ! 884 zdiffut(ji,jj,jk) = 0.8 * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_d_sc(ji,jj) * zznd_ml )**1.5863 zdiffut(ji,jj,jk) = zdifpyc_sc(ji,jj) * ( 1.0 + zznd_pyc ) 885 864 ! 886 zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_v_sc(ji,jj) * zznd_ml ) & 887 & * ( 1.0 - 0.5 * zznd_ml**2 ) 865 zviscos(ji,jj,jk) = zvispyc_sc(ji,jj) * ( 1.0 + zznd_pyc ) 888 866 END DO 889 ! pycnocline - if present linear profile 890 IF ( zdh(ji,jj) > 0._wp ) THEN 891 DO jk = imld(ji,jj)+1 , ibld(ji,jj) 892 zznd_pyc = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 893 ! 894 zdiffut(ji,jj,jk) = zdifpyc_sc(ji,jj) * ( 1.0 + zznd_pyc ) 895 ! 896 zviscos(ji,jj,jk) = zvispyc_sc(ji,jj) * ( 1.0 + zznd_pyc ) 897 END DO 898 ENDIF 899 ! Temporay fix to ensure zdiffut is +ve; won't be necessary with wn taken out 900 zdiffut(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj)* e3t_n(ji,jj,ibld(ji,jj)) 901 ! could be taken out, take account of entrainment represents as a diffusivity 902 ! should remove w from here, represents entrainment 903 ELSE 904 ! stable conditions 905 DO jk = 2, ibld(ji,jj) 906 zznd_ml = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 907 zdiffut(ji,jj,jk) = 0.75 * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zznd_ml )**1.5 908 zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * (1.0 - zznd_ml) * ( 1.0 - zznd_ml**2 ) 909 END DO 910 ENDIF ! end if ( lconv ) 867 ENDIF 868 ! Temporay fix to ensure zdiffut is +ve; won't be necessary with ww taken out 869 zdiffut(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj)* e3t(ji,jj,ibld(ji,jj),Kmm) 870 ! could be taken out, take account of entrainment represents as a diffusivity 871 ! should remove w from here, represents entrainment 872 ELSE 873 ! stable conditions 874 DO jk = 2, ibld(ji,jj) 875 zznd_ml = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 876 zdiffut(ji,jj,jk) = 0.75 * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zznd_ml )**1.5 877 zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * (1.0 - zznd_ml) * ( 1.0 - zznd_ml**2 ) 878 END DO 879 ENDIF ! end if ( lconv ) 911 880 ! 912 END DO ! end of ji loop 913 END DO ! end of jj loop 881 END_2D 914 882 915 883 ! … … 928 896 929 897 930 DO jj = 2, jpjm1 931 DO ji = 2, jpim1 932 IF ( lconv(ji,jj) ) THEN 933 DO jk = 2, imld(ji,jj) 934 zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj) 935 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) 936 ! 937 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) 938 END DO ! end jk loop 939 ELSE ! else for if (lconv) 898 DO_2D( 0, 0, 0, 0 ) 899 IF ( lconv(ji,jj) ) THEN 900 DO jk = 2, imld(ji,jj) 901 zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 902 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) 903 ! 904 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) 905 END DO ! end jk loop 906 ELSE ! else for if (lconv) 940 907 ! Stable conditions 941 DO jk = 2, ibld(ji,jj) 942 zznd_d=gdepw_n(ji,jj,jk) / dstokes(ji,jj) 943 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.5 * EXP ( -0.9 * zznd_d ) & 944 & * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_wth_1(ji,jj) 945 ! 946 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 1.5 * EXP ( -0.9 * zznd_d ) & 947 & * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_ws_1(ji,jj) 948 END DO 949 ENDIF ! endif for check on lconv 950 951 END DO ! end of ji loop 952 END DO ! end of jj loop 908 DO jk = 2, ibld(ji,jj) 909 zznd_d=gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 910 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.5 * EXP ( -0.9 * zznd_d ) & 911 & * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_wth_1(ji,jj) 912 ! 913 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 1.5 * EXP ( -0.9 * zznd_d ) & 914 & * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_ws_1(ji,jj) 915 END DO 916 ENDIF ! endif for check on lconv 917 918 END_2D 953 919 954 920 … … 963 929 ENDWHERE 964 930 965 DO jj = 2, jpjm1 966 DO ji = 2, jpim1 967 IF ( lconv(ji,jj) ) THEN 968 DO jk = 2, imld(ji,jj) 969 zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj) 970 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + ( -0.05 * EXP ( -0.4 * zznd_d ) * zsc_uw_1(ji,jj) & 971 & + 0.00125 * EXP ( - zznd_d ) * zsc_uw_2(ji,jj) ) & 972 & * ( 1.0 - EXP ( -2.0 * zznd_d ) ) 931 DO_2D( 0, 0, 0, 0 ) 932 IF ( lconv(ji,jj) ) THEN 933 DO jk = 2, imld(ji,jj) 934 zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 935 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + ( -0.05 * EXP ( -0.4 * zznd_d ) * zsc_uw_1(ji,jj) & 936 & + 0.00125 * EXP ( - zznd_d ) * zsc_uw_2(ji,jj) ) & 937 & * ( 1.0 - EXP ( -2.0 * zznd_d ) ) 973 938 ! 974 975 976 977 939 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.65 * 0.15 * EXP ( - zznd_d ) & 940 & * ( 1.0 - EXP ( -2.0 * zznd_d ) ) * zsc_vw_1(ji,jj) 941 END DO ! end jk loop 942 ELSE 978 943 ! Stable conditions 979 DO jk = 2, ibld(ji,jj) ! corrected to ibld 980 zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj) 981 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.75 * 1.3 * EXP ( -0.5 * zznd_d ) & 982 & * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_uw_1(ji,jj) 983 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + 0._wp 984 END DO ! end jk loop 985 ENDIF 986 END DO ! ji loop 987 END DO ! jj loo 944 DO jk = 2, ibld(ji,jj) ! corrected to ibld 945 zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 946 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.75 * 1.3 * EXP ( -0.5 * zznd_d ) & 947 & * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_uw_1(ji,jj) 948 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + 0._wp 949 END DO ! end jk loop 950 ENDIF 951 END_2D 988 952 989 953 ! Buoyancy term in flux-gradient relationship [note : includes ROI ratio (X0.3) and pressure (X0.5)] … … 997 961 ENDWHERE 998 962 999 DO jj = 2, jpjm1 1000 DO ji = 2, jpim1 1001 IF (lconv(ji,jj) ) THEN 1002 DO jk = 2, imld(ji,jj) 1003 zznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj) 1004 ! calculate turbulent length scale 1005 zl_c = 0.9 * ( 1.0 - EXP ( - 7.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) ) & 1006 & * ( 1.0 - EXP ( -15.0 * ( 1.1 - zznd_ml ) ) ) 1007 zl_l = 2.0 * ( 1.0 - EXP ( - 2.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) ) & 1008 & * ( 1.0 - EXP ( - 5.0 * ( 1.0 - zznd_ml ) ) ) * ( 1.0 + dstokes(ji,jj) / zhml (ji,jj) ) 1009 zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0 + EXP ( 3.0 * LOG10 ( - zhol(ji,jj) ) ) ) ** (3.0/2.0) 1010 ! non-gradient buoyancy terms 1011 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * 0.5 * zsc_wth_1(ji,jj) * zl_eps * zhml(ji,jj) / ( 0.15 + zznd_ml ) 1012 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 ) 1013 END DO 1014 ELSE 1015 DO jk = 2, ibld(ji,jj) 1016 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zsc_wth_1(ji,jj) 1017 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zsc_ws_1(ji,jj) 1018 END DO 1019 ENDIF 1020 END DO ! ji loop 1021 END DO ! jj loop 963 DO_2D( 0, 0, 0, 0 ) 964 IF (lconv(ji,jj) ) THEN 965 DO jk = 2, imld(ji,jj) 966 zznd_ml = gdepw(ji,jj,jk,Kmm) / zhml(ji,jj) 967 ! calculate turbulent length scale 968 zl_c = 0.9 * ( 1.0 - EXP ( - 7.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) ) & 969 & * ( 1.0 - EXP ( -15.0 * ( 1.1 - zznd_ml ) ) ) 970 zl_l = 2.0 * ( 1.0 - EXP ( - 2.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) ) & 971 & * ( 1.0 - EXP ( - 5.0 * ( 1.0 - zznd_ml ) ) ) * ( 1.0 + dstokes(ji,jj) / zhml (ji,jj) ) 972 zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0 + EXP ( 3.0 * LOG10 ( - zhol(ji,jj) ) ) ) ** (3.0/2.0) 973 ! non-gradient buoyancy terms 974 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 ) 975 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 ) 976 END DO 977 ELSE 978 DO jk = 2, ibld(ji,jj) 979 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zsc_wth_1(ji,jj) 980 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zsc_ws_1(ji,jj) 981 END DO 982 ENDIF 983 END_2D 1022 984 1023 985 … … 1031 993 ENDWHERE 1032 994 1033 DO jj = 2, jpjm1 1034 DO ji = 2, jpim1 1035 IF ( lconv(ji,jj) ) THEN 1036 DO jk = 2 , imld(ji,jj) 1037 zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj) 1038 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.3 * 0.5 * ( zsc_uw_1(ji,jj) + 0.125 * EXP( -0.5 * zznd_d ) & 1039 & * ( 1.0 - EXP( -0.5 * zznd_d ) ) & 1040 & * zsc_uw_2(ji,jj) ) 1041 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj) 1042 END DO ! jk loop 1043 ELSE 1044 ! stable conditions 1045 DO jk = 2, ibld(ji,jj) 1046 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zsc_uw_1(ji,jj) 1047 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj) 1048 END DO 1049 ENDIF 1050 END DO ! ji loop 1051 END DO ! jj loop 995 DO_2D( 0, 0, 0, 0 ) 996 IF ( lconv(ji,jj) ) THEN 997 DO jk = 2 , imld(ji,jj) 998 zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 999 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.3 * 0.5 * ( zsc_uw_1(ji,jj) + 0.125 * EXP( -0.5 * zznd_d ) & 1000 & * ( 1.0 - EXP( -0.5 * zznd_d ) ) & 1001 & * zsc_uw_2(ji,jj) ) 1002 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj) 1003 END DO ! jk loop 1004 ELSE 1005 ! stable conditions 1006 DO jk = 2, ibld(ji,jj) 1007 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zsc_uw_1(ji,jj) 1008 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj) 1009 END DO 1010 ENDIF 1011 END_2D 1052 1012 1053 1013 ! Transport term in flux-gradient relationship [note : includes ROI ratio (X0.3) ] … … 1061 1021 ENDWHERE 1062 1022 1063 DO jj = 2, jpjm1 1064 DO ji = 2, jpim1 1065 IF ( lconv(ji,jj) ) THEN 1066 DO jk = 2, imld(ji,jj) 1067 zznd_ml=gdepw_n(ji,jj,jk) / zhml(ji,jj) 1068 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * zsc_wth_1(ji,jj) & 1069 & * ( -2.0 + 2.75 * ( ( 1.0 + 0.6 * zznd_ml**4 ) & 1070 & - EXP( - 6.0 * zznd_ml ) ) ) & 1071 & * ( 1.0 - EXP( - 15.0 * ( 1.0 - zznd_ml ) ) ) 1072 ! 1073 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * zsc_ws_1(ji,jj) & 1074 & * ( -2.0 + 2.75 * ( ( 1.0 + 0.6 * zznd_ml**4 ) & 1075 & - EXP( - 6.0 * zznd_ml ) ) ) & 1076 & * ( 1.0 - EXP ( -15.0 * ( 1.0 - zznd_ml ) ) ) 1077 END DO 1078 ELSE 1079 DO jk = 2, ibld(ji,jj) 1080 zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj) 1081 znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 1082 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + & 1083 & 7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_wth_1(ji,jj) 1084 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + & 1085 & 7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_ws_1(ji,jj) 1086 END DO 1087 ENDIF 1088 ENDDO ! ji loop 1089 END DO ! jj loop 1023 DO_2D( 0, 0, 0, 0 ) 1024 IF ( lconv(ji,jj) ) THEN 1025 DO jk = 2, imld(ji,jj) 1026 zznd_ml=gdepw(ji,jj,jk,Kmm) / zhml(ji,jj) 1027 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * zsc_wth_1(ji,jj) & 1028 & * ( -2.0 + 2.75 * ( ( 1.0 + 0.6 * zznd_ml**4 ) & 1029 & - EXP( - 6.0 * zznd_ml ) ) ) & 1030 & * ( 1.0 - EXP( - 15.0 * ( 1.0 - zznd_ml ) ) ) 1031 ! 1032 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * zsc_ws_1(ji,jj) & 1033 & * ( -2.0 + 2.75 * ( ( 1.0 + 0.6 * zznd_ml**4 ) & 1034 & - EXP( - 6.0 * zznd_ml ) ) ) & 1035 & * ( 1.0 - EXP ( -15.0 * ( 1.0 - zznd_ml ) ) ) 1036 END DO 1037 ELSE 1038 DO jk = 2, ibld(ji,jj) 1039 zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 1040 znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 1041 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + & 1042 & 7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_wth_1(ji,jj) 1043 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + & 1044 & 7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_ws_1(ji,jj) 1045 END DO 1046 ENDIF 1047 END_2D 1090 1048 1091 1049 … … 1100 1058 ENDWHERE 1101 1059 1102 DO jj = 2, jpjm1 1103 DO ji = 2, jpim1 1104 IF ( lconv(ji,jj) ) THEN 1105 DO jk = 2, imld(ji,jj) 1106 zznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj) 1107 zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj) 1060 DO_2D( 0, 0, 0, 0 ) 1061 IF ( lconv(ji,jj) ) THEN 1062 DO jk = 2, imld(ji,jj) 1063 zznd_ml = gdepw(ji,jj,jk,Kmm) / zhml(ji,jj) 1064 zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 1065 ghamu(ji,jj,jk) = ghamu(ji,jj,jk)& 1066 & + 0.3 * ( -2.0 + 2.5 * ( 1.0 + 0.1 * zznd_ml**4 ) - EXP ( -8.0 * zznd_ml ) ) * zsc_uw_1(ji,jj) 1067 ! 1068 ghamv(ji,jj,jk) = ghamv(ji,jj,jk)& 1069 & + 0.3 * 0.1 * ( EXP( -zznd_d ) + EXP( -5.0 * ( 1.0 - zznd_ml ) ) ) * zsc_vw_1(ji,jj) 1070 END DO 1071 ELSE 1072 DO jk = 2, ibld(ji,jj) 1073 znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 1074 zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 1075 IF ( zznd_d <= 2.0 ) THEN 1076 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.5 * 0.3 & 1077 &* ( 2.25 - 3.0 * ( 1.0 - EXP( - 1.25 * zznd_d ) ) * ( 1.0 - EXP( -2.0 * zznd_d ) ) ) * zsc_uw_1(ji,jj) 1078 ! 1079 ELSE 1108 1080 ghamu(ji,jj,jk) = ghamu(ji,jj,jk)& 1109 & + 0. 3 * ( -2.0 + 2.5 * ( 1.0 + 0.1 * zznd_ml**4 ) - EXP ( -8.0 * zznd_ml ) ) * zsc_uw_1(ji,jj)1081 & + 0.5 * 0.3 * ( 1.0 - EXP( -5.0 * ( 1.0 - znd ) ) ) * zsc_uw_2(ji,jj) 1110 1082 ! 1111 ghamv(ji,jj,jk) = ghamv(ji,jj,jk)& 1112 & + 0.3 * 0.1 * ( EXP( -zznd_d ) + EXP( -5.0 * ( 1.0 - zznd_ml ) ) ) * zsc_vw_1(ji,jj) 1113 END DO 1114 ELSE 1115 DO jk = 2, ibld(ji,jj) 1116 znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 1117 zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj) 1118 IF ( zznd_d <= 2.0 ) THEN 1119 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.5 * 0.3 & 1120 &* ( 2.25 - 3.0 * ( 1.0 - EXP( - 1.25 * zznd_d ) ) * ( 1.0 - EXP( -2.0 * zznd_d ) ) ) * zsc_uw_1(ji,jj) 1121 ! 1122 ELSE 1123 ghamu(ji,jj,jk) = ghamu(ji,jj,jk)& 1124 & + 0.5 * 0.3 * ( 1.0 - EXP( -5.0 * ( 1.0 - znd ) ) ) * zsc_uw_2(ji,jj) 1125 ! 1126 ENDIF 1127 1128 ghamv(ji,jj,jk) = ghamv(ji,jj,jk)& 1129 & + 0.3 * 0.15 * SIN( 3.14159 * ( 0.65 * zznd_d ) ) * EXP( -0.25 * zznd_d**2 ) * zsc_vw_1(ji,jj) 1130 ghamv(ji,jj,jk) = ghamv(ji,jj,jk)& 1131 & + 0.3 * 0.15 * EXP( -5.0 * ( 1.0 - znd ) ) * ( 1.0 - EXP( -20.0 * ( 1.0 - znd ) ) ) * zsc_vw_2(ji,jj) 1132 END DO 1133 ENDIF 1134 END DO 1135 END DO 1083 ENDIF 1084 1085 ghamv(ji,jj,jk) = ghamv(ji,jj,jk)& 1086 & + 0.3 * 0.15 * SIN( 3.14159 * ( 0.65 * zznd_d ) ) * EXP( -0.25 * zznd_d**2 ) * zsc_vw_1(ji,jj) 1087 ghamv(ji,jj,jk) = ghamv(ji,jj,jk)& 1088 & + 0.3 * 0.15 * EXP( -5.0 * ( 1.0 - znd ) ) * ( 1.0 - EXP( -20.0 * ( 1.0 - znd ) ) ) * zsc_vw_2(ji,jj) 1089 END DO 1090 ENDIF 1091 END_2D 1136 1092 ! 1137 1093 ! Make surface forced velocity non-gradient terms go to zero at the base of the mixed layer. 1138 1094 1139 DO jj = 2, jpjm1 1140 DO ji = 2, jpim1 1141 IF ( lconv(ji,jj) ) THEN 1142 DO jk = 2, ibld(ji,jj) 1143 znd = ( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zhml(ji,jj) !ALMG to think about 1144 IF ( znd >= 0.0 ) THEN 1145 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -30.0 * znd**2 ) ) 1146 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * ( 1.0 - EXP( -30.0 * znd**2 ) ) 1147 ELSE 1148 ghamu(ji,jj,jk) = 0._wp 1149 ghamv(ji,jj,jk) = 0._wp 1150 ENDIF 1151 END DO 1152 ELSE 1153 DO jk = 2, ibld(ji,jj) 1154 znd = ( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zhml(ji,jj) !ALMG to think about 1155 IF ( znd >= 0.0 ) THEN 1156 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) ) 1157 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) ) 1158 ELSE 1159 ghamu(ji,jj,jk) = 0._wp 1160 ghamv(ji,jj,jk) = 0._wp 1161 ENDIF 1162 END DO 1163 ENDIF 1164 END DO 1165 END DO 1095 DO_2D( 0, 0, 0, 0 ) 1096 IF ( lconv(ji,jj) ) THEN 1097 DO jk = 2, ibld(ji,jj) 1098 znd = ( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zhml(ji,jj) !ALMG to think about 1099 IF ( znd >= 0.0 ) THEN 1100 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -30.0 * znd**2 ) ) 1101 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * ( 1.0 - EXP( -30.0 * znd**2 ) ) 1102 ELSE 1103 ghamu(ji,jj,jk) = 0._wp 1104 ghamv(ji,jj,jk) = 0._wp 1105 ENDIF 1106 END DO 1107 ELSE 1108 DO jk = 2, ibld(ji,jj) 1109 znd = ( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zhml(ji,jj) !ALMG to think about 1110 IF ( znd >= 0.0 ) THEN 1111 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) ) 1112 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) ) 1113 ELSE 1114 ghamu(ji,jj,jk) = 0._wp 1115 ghamv(ji,jj,jk) = 0._wp 1116 ENDIF 1117 END DO 1118 ENDIF 1119 END_2D 1166 1120 1167 1121 ! pynocline contributions 1168 1122 ! Temporary fix to avoid instabilities when zdb_bl becomes very very small 1169 1123 zsc_uw_1 = 0._wp ! 50.0 * zla**(8.0/3.0) * zustar**2 * zhbl / ( zdb_bl + epsln ) 1170 DO jj = 2, jpjm1 1171 DO ji = 2, jpim1 1172 DO jk= 2, ibld(ji,jj) 1173 znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 1174 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk) 1175 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk) 1176 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc(ji,jj,jk) 1177 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zsc_uw_1(ji,jj) * ( 1.0 - znd )**(7.0/4.0) * zdbdz_pyc(ji,jj,jk) 1178 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc(ji,jj,jk) 1124 DO_2D( 0, 0, 0, 0 ) 1125 DO jk= 2, ibld(ji,jj) 1126 znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 1127 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk) 1128 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk) 1129 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc(ji,jj,jk) 1130 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) 1131 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc(ji,jj,jk) 1132 END DO 1133 END_2D 1134 1135 ! Entrainment contribution. 1136 1137 DO_2D( 0, 0, 0, 0 ) 1138 IF ( lconv(ji,jj) ) THEN 1139 DO jk = 1, imld(ji,jj) - 1 1140 znd=gdepw(ji,jj,jk,Kmm) / zhml(ji,jj) 1141 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * znd 1142 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * znd 1143 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * znd 1144 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * znd 1145 END DO 1146 DO jk = imld(ji,jj), ibld(ji,jj) 1147 znd = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj) 1148 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * ( 1.0 + znd ) 1149 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * ( 1.0 + znd ) 1150 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * ( 1.0 + znd ) 1151 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * ( 1.0 + znd ) 1179 1152 END DO 1180 END DO 1181 END DO 1182 1183 ! Entrainment contribution. 1184 1185 DO jj=2, jpjm1 1186 DO ji = 2, jpim1 1187 IF ( lconv(ji,jj) ) THEN 1188 DO jk = 1, imld(ji,jj) - 1 1189 znd=gdepw_n(ji,jj,jk) / zhml(ji,jj) 1190 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * znd 1191 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * znd 1192 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * znd 1193 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * znd 1194 END DO 1195 DO jk = imld(ji,jj), ibld(ji,jj) 1196 znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 1197 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * ( 1.0 + znd ) 1198 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * ( 1.0 + znd ) 1199 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * ( 1.0 + znd ) 1200 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * ( 1.0 + znd ) 1201 END DO 1202 ENDIF 1203 ghamt(ji,jj,ibld(ji,jj)) = 0._wp 1204 ghams(ji,jj,ibld(ji,jj)) = 0._wp 1205 ghamu(ji,jj,ibld(ji,jj)) = 0._wp 1206 ghamv(ji,jj,ibld(ji,jj)) = 0._wp 1207 END DO ! ji loop 1208 END DO ! jj loop 1153 ENDIF 1154 ghamt(ji,jj,ibld(ji,jj)) = 0._wp 1155 ghams(ji,jj,ibld(ji,jj)) = 0._wp 1156 ghamu(ji,jj,ibld(ji,jj)) = 0._wp 1157 ghamv(ji,jj,ibld(ji,jj)) = 0._wp 1158 END_2D 1209 1159 1210 1160 … … 1220 1170 ! rotate non-gradient velocity terms back to model reference frame 1221 1171 1222 DO jj = 2, jpjm1 1223 DO ji = 2, jpim1 1224 DO jk = 2, ibld(ji,jj) 1225 ztemp = ghamu(ji,jj,jk) 1226 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * zcos_wind(ji,jj) - ghamv(ji,jj,jk) * zsin_wind(ji,jj) 1227 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * zcos_wind(ji,jj) + ztemp * zsin_wind(ji,jj) 1228 END DO 1172 DO_2D( 0, 0, 0, 0 ) 1173 DO jk = 2, ibld(ji,jj) 1174 ztemp = ghamu(ji,jj,jk) 1175 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * zcos_wind(ji,jj) - ghamv(ji,jj,jk) * zsin_wind(ji,jj) 1176 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * zcos_wind(ji,jj) + ztemp * zsin_wind(ji,jj) 1229 1177 END DO 1230 END DO1178 END_2D 1231 1179 1232 1180 IF(ln_dia_osm) THEN … … 1236 1184 ! KPP-style Ri# mixing 1237 1185 IF( ln_kpprimix) THEN 1238 DO jk = 2, jpkm1 !* Shear production at uw- and vw-points (energy conserving form) 1239 DO jj = 1, jpjm1 1240 DO ji = 1, jpim1 ! vector opt. 1241 z3du(ji,jj,jk) = 0.5 * ( un(ji,jj,jk-1) - un(ji ,jj,jk) ) & 1242 & * ( ub(ji,jj,jk-1) - ub(ji ,jj,jk) ) * wumask(ji,jj,jk) & 1243 & / ( e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) ) 1244 z3dv(ji,jj,jk) = 0.5 * ( vn(ji,jj,jk-1) - vn(ji,jj ,jk) ) & 1245 & * ( vb(ji,jj,jk-1) - vb(ji,jj ,jk) ) * wvmask(ji,jj,jk) & 1246 & / ( e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) ) 1247 END DO 1186 DO_3D( 1, 0, 1, 0, 2, jpkm1 ) 1187 z3du(ji,jj,jk) = 0.5 * ( uu(ji,jj,jk-1,Kmm) - uu(ji ,jj,jk,Kmm) ) & 1188 & * ( uu(ji,jj,jk-1,Kbb) - uu(ji ,jj,jk,Kbb) ) * wumask(ji,jj,jk) & 1189 & / ( e3uw(ji,jj,jk,Kmm) * e3uw(ji,jj,jk,Kbb) ) 1190 z3dv(ji,jj,jk) = 0.5 * ( vv(ji,jj,jk-1,Kmm) - vv(ji,jj ,jk,Kmm) ) & 1191 & * ( vv(ji,jj,jk-1,Kbb) - vv(ji,jj ,jk,Kbb) ) * wvmask(ji,jj,jk) & 1192 & / ( e3vw(ji,jj,jk,Kmm) * e3vw(ji,jj,jk,Kbb) ) 1193 END_3D 1194 ! 1195 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 1196 ! ! shear prod. at w-point weightened by mask 1197 zesh2 = ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) & 1198 & + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) ) 1199 ! ! local Richardson number 1200 zri = MAX( rn2b(ji,jj,jk), 0._wp ) / MAX(zesh2, epsln) 1201 zfri = MIN( zri / rn_riinfty , 1.0_wp ) 1202 zfri = ( 1.0_wp - zfri * zfri ) 1203 zrimix(ji,jj,jk) = zfri * zfri * zfri * wmask(ji, jj, jk) 1204 END_3D 1205 1206 DO_2D( 0, 0, 0, 0 ) 1207 DO jk = ibld(ji,jj) + 1, jpkm1 1208 zdiffut(ji,jj,jk) = zrimix(ji,jj,jk)*rn_difri 1209 zviscos(ji,jj,jk) = zrimix(ji,jj,jk)*rn_difri 1248 1210 END DO 1249 END DO 1250 ! 1251 DO jk = 2, jpkm1 1252 DO jj = 2, jpjm1 1253 DO ji = 2, jpim1 ! vector opt. 1254 ! ! shear prod. at w-point weightened by mask 1255 zesh2 = ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) & 1256 & + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) ) 1257 ! ! local Richardson number 1258 zri = MAX( rn2b(ji,jj,jk), 0._wp ) / MAX(zesh2, epsln) 1259 zfri = MIN( zri / rn_riinfty , 1.0_wp ) 1260 zfri = ( 1.0_wp - zfri * zfri ) 1261 zrimix(ji,jj,jk) = zfri * zfri * zfri * wmask(ji, jj, jk) 1262 END DO 1263 END DO 1264 END DO 1265 1266 DO jj = 2, jpjm1 1267 DO ji = 2, jpim1 1268 DO jk = ibld(ji,jj) + 1, jpkm1 1269 zdiffut(ji,jj,jk) = zrimix(ji,jj,jk)*rn_difri 1270 zviscos(ji,jj,jk) = zrimix(ji,jj,jk)*rn_difri 1271 END DO 1272 END DO 1273 END DO 1211 END_2D 1274 1212 1275 1213 END IF ! ln_kpprimix = .true. … … 1277 1215 ! KPP-style set diffusivity large if unstable below BL 1278 1216 IF( ln_convmix) THEN 1279 DO jj = 2, jpjm1 1280 DO ji = 2, jpim1 1281 DO jk = ibld(ji,jj) + 1, jpkm1 1282 IF( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) zdiffut(ji,jj,jk) = rn_difconv 1283 END DO 1217 DO_2D( 0, 0, 0, 0 ) 1218 DO jk = ibld(ji,jj) + 1, jpkm1 1219 IF( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) zdiffut(ji,jj,jk) = rn_difconv 1284 1220 END DO 1285 END DO1221 END_2D 1286 1222 END IF ! ln_convmix = .true. 1287 1223 1288 1224 ! Lateral boundary conditions on zvicos (sign unchanged), needed to caclulate viscosities on u and v grids 1289 CALL lbc_lnk( 'zdfosm', zviscos(:,:,:), 'W', 1. )1225 CALL lbc_lnk( 'zdfosm', zviscos(:,:,:), 'W', 1.0_wp ) 1290 1226 1291 1227 ! GN 25/8: need to change tmask --> wmask 1292 1228 1293 DO jk = 2, jpkm1 1294 DO jj = 2, jpjm1 1295 DO ji = 2, jpim1 1296 p_avt(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), avtb(jk) ) * tmask(ji,jj,jk) 1297 p_avm(ji,jj,jk) = MAX( zviscos(ji,jj,jk), avmb(jk) ) * tmask(ji,jj,jk) 1298 END DO 1299 END DO 1300 END DO 1229 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 1230 p_avt(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), avtb(jk) ) * tmask(ji,jj,jk) 1231 p_avm(ji,jj,jk) = MAX( zviscos(ji,jj,jk), avmb(jk) ) * tmask(ji,jj,jk) 1232 END_3D 1301 1233 ! Lateral boundary conditions on ghamu and ghamv, currently on W-grid (sign unchanged), needed to caclulate gham[uv] on u and v grids 1302 CALL lbc_lnk_multi( 'zdfosm', p_avt, 'W', 1. , p_avm, 'W', 1., & 1303 & ghamu, 'W', 1. , ghamv, 'W', 1. ) 1304 DO jk = 2, jpkm1 1305 DO jj = 2, jpjm1 1306 DO ji = 2, jpim1 1307 ghamu(ji,jj,jk) = ( ghamu(ji,jj,jk) + ghamu(ji+1,jj,jk) ) & 1308 & / MAX( 1., tmask(ji,jj,jk) + tmask (ji + 1,jj,jk) ) * umask(ji,jj,jk) 1309 1310 ghamv(ji,jj,jk) = ( ghamv(ji,jj,jk) + ghamv(ji,jj+1,jk) ) & 1311 & / MAX( 1., tmask(ji,jj,jk) + tmask (ji,jj+1,jk) ) * vmask(ji,jj,jk) 1312 1313 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) * tmask(ji,jj,jk) 1314 ghams(ji,jj,jk) = ghams(ji,jj,jk) * tmask(ji,jj,jk) 1315 END DO 1316 END DO 1317 END DO 1234 CALL lbc_lnk_multi( 'zdfosm', p_avt, 'W', 1.0_wp , p_avm, 'W', 1.0_wp, & 1235 & ghamu, 'W', 1.0_wp , ghamv, 'W', 1.0_wp ) 1236 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 1237 ghamu(ji,jj,jk) = ( ghamu(ji,jj,jk) + ghamu(ji+1,jj,jk) ) & 1238 & / MAX( 1., tmask(ji,jj,jk) + tmask (ji + 1,jj,jk) ) * umask(ji,jj,jk) 1239 1240 ghamv(ji,jj,jk) = ( ghamv(ji,jj,jk) + ghamv(ji,jj+1,jk) ) & 1241 & / MAX( 1., tmask(ji,jj,jk) + tmask (ji,jj+1,jk) ) * vmask(ji,jj,jk) 1242 1243 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) * tmask(ji,jj,jk) 1244 ghams(ji,jj,jk) = ghams(ji,jj,jk) * tmask(ji,jj,jk) 1245 END_3D 1318 1246 ! Lateral boundary conditions on final outputs for gham[ts], on W-grid (sign unchanged) 1319 1247 ! Lateral boundary conditions on final outputs for gham[uv], on [UV]-grid (sign unchanged) 1320 CALL lbc_lnk_multi( 'zdfosm', ghamt, 'W', 1. , ghams, 'W', 1., &1321 & ghamu, 'U', 1. , ghamv, 'V', 1.)1248 CALL lbc_lnk_multi( 'zdfosm', ghamt, 'W', 1.0_wp , ghams, 'W', 1.0_wp, & 1249 & ghamu, 'U', 1.0_wp , ghamv, 'V', 1.0_wp ) 1322 1250 1323 1251 IF(ln_dia_osm) THEN … … 1327 1255 IF ( iom_use("us_x") ) CALL iom_put( "us_x", tmask(:,:,1)*zustke*zcos_wind ) ! x surface Stokes drift 1328 1256 IF ( iom_use("us_y") ) CALL iom_put( "us_y", tmask(:,:,1)*zustke*zsin_wind ) ! y surface Stokes drift 1329 IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*r au0*tmask(:,:,1)*zustar**2*zustke )1257 IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rho0*tmask(:,:,1)*zustar**2*zustke ) 1330 1258 ! Stokes drift read in from sbcwave (=2). 1331 1259 CASE(2) 1332 1260 IF ( iom_use("us_x") ) CALL iom_put( "us_x", ut0sd ) ! x surface Stokes drift 1333 1261 IF ( iom_use("us_y") ) CALL iom_put( "us_y", vt0sd ) ! y surface Stokes drift 1334 IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*r au0*tmask(:,:,1)*zustar**2* &1262 IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rho0*tmask(:,:,1)*zustar**2* & 1335 1263 & SQRT(ut0sd**2 + vt0sd**2 ) ) 1336 1264 END SELECT … … 1348 1276 IF ( iom_use("zwstrl") ) CALL iom_put( "zwstrl", tmask(:,:,1)*zwstrl ) ! Langmuir velocity scale 1349 1277 IF ( iom_use("zustar") ) CALL iom_put( "zustar", tmask(:,:,1)*zustar ) ! friction velocity scale 1350 IF ( iom_use("wind_power") ) CALL iom_put( "wind_power", 1000.*r au0*tmask(:,:,1)*zustar**3 ) ! BL depth internal to zdf_osm routine1351 IF ( iom_use("wind_wave_power") ) CALL iom_put( "wind_wave_power", 1000.*r au0*tmask(:,:,1)*zustar**2*zustke )1278 IF ( iom_use("wind_power") ) CALL iom_put( "wind_power", 1000.*rho0*tmask(:,:,1)*zustar**3 ) ! BL depth internal to zdf_osm routine 1279 IF ( iom_use("wind_wave_power") ) CALL iom_put( "wind_wave_power", 1000.*rho0*tmask(:,:,1)*zustar**2*zustke ) 1352 1280 IF ( iom_use("zhbl") ) CALL iom_put( "zhbl", tmask(:,:,1)*zhbl ) ! BL depth internal to zdf_osm routine 1353 1281 IF ( iom_use("zhml") ) CALL iom_put( "zhml", tmask(:,:,1)*zhml ) ! ML depth internal to zdf_osm routine … … 1359 1287 END IF 1360 1288 ! Lateral boundary conditions on p_avt (sign unchanged) 1361 CALL lbc_lnk( 'zdfosm', p_avt(:,:,:), 'W', 1. )1289 CALL lbc_lnk( 'zdfosm', p_avt(:,:,:), 'W', 1.0_wp ) 1362 1290 ! 1363 1291 END SUBROUTINE zdf_osm 1364 1292 1365 1293 1366 SUBROUTINE zdf_osm_init 1294 SUBROUTINE zdf_osm_init( Kmm ) 1367 1295 !!---------------------------------------------------------------------- 1368 1296 !! *** ROUTINE zdf_osm_init *** … … 1376 1304 !! ** input : Namlist namosm 1377 1305 !!---------------------------------------------------------------------- 1306 INTEGER, INTENT(in) :: Kmm ! time level index (middle) 1307 ! 1378 1308 INTEGER :: ios ! local integer 1379 1309 INTEGER :: ji, jj, jk ! dummy loop indices … … 1384 1314 !!---------------------------------------------------------------------- 1385 1315 ! 1386 REWIND( numnam_ref ) ! Namelist namzdf_osm in reference namelist : Osmosis ML model1387 1316 READ ( numnam_ref, namzdf_osm, IOSTAT = ios, ERR = 901) 1388 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_osm in reference namelist', lwp ) 1389 1390 REWIND( numnam_cfg ) ! Namelist namzdf_tke in configuration namelist : Turbulent Kinetic Energy 1317 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_osm in reference namelist' ) 1318 1391 1319 READ ( numnam_cfg, namzdf_osm, IOSTAT = ios, ERR = 902 ) 1392 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namzdf_osm in configuration namelist' , lwp)1320 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namzdf_osm in configuration namelist' ) 1393 1321 IF(lwm) WRITE ( numond, namzdf_osm ) 1394 1322 … … 1420 1348 ENDIF 1421 1349 1350 1351 ! ! Check wave coupling settings ! 1352 ! ! Further work needed - see ticket #2447 ! 1353 IF( nn_osm_wave == 2 ) THEN 1354 IF (.NOT. ( ln_wave .AND. ln_sdw )) & 1355 & CALL ctl_stop( 'zdf_osm_init : ln_zdfosm and nn_osm_wave=2, ln_wave and ln_sdw must be true' ) 1356 END IF 1357 1422 1358 ! ! allocate zdfosm arrays 1423 1359 IF( zdf_osm_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_osm_init : unable to allocate arrays' ) 1424 1360 1425 call osm_rst( nit000, 'READ' ) !* read or initialize hbl1361 call osm_rst( nit000, Kmm, 'READ' ) !* read or initialize hbl 1426 1362 1427 1363 IF( ln_zdfddm) THEN … … 1459 1395 etmean(:,:,:) = 0.e0 1460 1396 1461 DO jk = 1, jpkm1 1462 DO jj = 2, jpjm1 1463 DO ji = 2, jpim1 ! vector opt. 1464 etmean(ji,jj,jk) = tmask(ji,jj,jk) & 1465 & / MAX( 1., umask(ji-1,jj ,jk) + umask(ji,jj,jk) & 1466 & + vmask(ji ,jj-1,jk) + vmask(ji,jj,jk) ) 1467 END DO 1468 END DO 1469 END DO 1397 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 1398 etmean(ji,jj,jk) = tmask(ji,jj,jk) & 1399 & / MAX( 1., umask(ji-1,jj ,jk) + umask(ji,jj,jk) & 1400 & + vmask(ji ,jj-1,jk) + vmask(ji,jj,jk) ) 1401 END_3D 1470 1402 1471 1403 CASE ( 1 ) ! horizontal average … … 1477 1409 etmean(:,:,:) = 0.e0 1478 1410 1479 DO jk = 1, jpkm1 1480 DO jj = 2, jpjm1 1481 DO ji = 2, jpim1 ! vector opt. 1482 etmean(ji,jj,jk) = tmask(ji, jj,jk) & 1483 & / MAX( 1., 2.* tmask(ji,jj,jk) & 1484 & +.5 * ( tmask(ji-1,jj+1,jk) + tmask(ji-1,jj-1,jk) & 1485 & +tmask(ji+1,jj+1,jk) + tmask(ji+1,jj-1,jk) ) & 1486 & +1. * ( tmask(ji-1,jj ,jk) + tmask(ji ,jj+1,jk) & 1487 & +tmask(ji ,jj-1,jk) + tmask(ji+1,jj ,jk) ) ) 1488 END DO 1489 END DO 1490 END DO 1411 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 1412 etmean(ji,jj,jk) = tmask(ji, jj,jk) & 1413 & / MAX( 1., 2.* tmask(ji,jj,jk) & 1414 & +.5 * ( tmask(ji-1,jj+1,jk) + tmask(ji-1,jj-1,jk) & 1415 & +tmask(ji+1,jj+1,jk) + tmask(ji+1,jj-1,jk) ) & 1416 & +1. * ( tmask(ji-1,jj ,jk) + tmask(ji ,jj+1,jk) & 1417 & +tmask(ji ,jj-1,jk) + tmask(ji+1,jj ,jk) ) ) 1418 END_3D 1491 1419 1492 1420 CASE DEFAULT … … 1517 1445 1518 1446 1519 SUBROUTINE osm_rst( kt, cdrw )1447 SUBROUTINE osm_rst( kt, Kmm, cdrw ) 1520 1448 !!--------------------------------------------------------------------- 1521 1449 !! *** ROUTINE osm_rst *** … … 1527 1455 !!---------------------------------------------------------------------- 1528 1456 1529 INTEGER, INTENT(in) :: kt 1457 INTEGER , INTENT(in) :: kt ! ocean time step index 1458 INTEGER , INTENT(in) :: Kmm ! ocean time level index (middle) 1530 1459 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 1531 1460 … … 1546 1475 id1 = iom_varid( numror, 'wn' , ldstop = .FALSE. ) 1547 1476 IF( id1 > 0 ) THEN ! 'wn' exists; read 1548 CALL iom_get( numror, jpdom_auto glo, 'wn', wn, ldxios = lrxios )1549 WRITE(numout,*) ' ===>>>> : w nread from restart file'1477 CALL iom_get( numror, jpdom_auto, 'wn', ww, ldxios = lrxios ) 1478 WRITE(numout,*) ' ===>>>> : ww read from restart file' 1550 1479 ELSE 1551 w n(:,:,:) = 0._wp1552 WRITE(numout,*) ' ===>>>> : w nnot in restart file, set to zero initially'1480 ww(:,:,:) = 0._wp 1481 WRITE(numout,*) ' ===>>>> : ww not in restart file, set to zero initially' 1553 1482 END IF 1554 1483 id1 = iom_varid( numror, 'hbl' , ldstop = .FALSE. ) 1555 1484 id2 = iom_varid( numror, 'hbli' , ldstop = .FALSE. ) 1556 1485 IF( id1 > 0 .AND. id2 > 0) THEN ! 'hbl' exists; read and return 1557 CALL iom_get( numror, jpdom_auto glo, 'hbl' , hbl , ldxios = lrxios )1558 CALL iom_get( numror, jpdom_auto glo, 'hbli', hbli, ldxios = lrxios )1486 CALL iom_get( numror, jpdom_auto, 'hbl' , hbl , ldxios = lrxios ) 1487 CALL iom_get( numror, jpdom_auto, 'hbli', hbli, ldxios = lrxios ) 1559 1488 WRITE(numout,*) ' ===>>>> : hbl & hbli read from restart file' 1560 1489 RETURN … … 1570 1499 IF( TRIM(cdrw) == 'WRITE') THEN !* Write hbli into the restart file, then return 1571 1500 IF(lwp) WRITE(numout,*) '---- osm-rst ----' 1572 CALL iom_rstput( kt, nitrst, numrow, 'wn' , w n, ldxios = lwxios )1501 CALL iom_rstput( kt, nitrst, numrow, 'wn' , ww , ldxios = lwxios ) 1573 1502 CALL iom_rstput( kt, nitrst, numrow, 'hbl' , hbl , ldxios = lwxios ) 1574 1503 CALL iom_rstput( kt, nitrst, numrow, 'hbli' , hbli, ldxios = lwxios ) … … 1582 1511 ALLOCATE( imld_rst(jpi,jpj) ) 1583 1512 ! w-level of the mixing and mixed layers 1584 CALL eos_rab( ts n, rab_n)1585 CALL bn2(ts n, rab_n, rn2)1513 CALL eos_rab( ts(:,:,:,:,Kmm), rab_n, Kmm ) 1514 CALL bn2(ts(:,:,:,:,Kmm), rab_n, rn2, Kmm) 1586 1515 imld_rst(:,:) = nlb10 ! Initialization to the number of w ocean point 1587 1516 hbl(:,:) = 0._wp ! here hbl used as a dummy variable, integrating vertically N^2 1588 zN2_c = grav * rho_c * r1_r au0 ! convert density criteria into N^2 criteria1517 zN2_c = grav * rho_c * r1_rho0 ! convert density criteria into N^2 criteria 1589 1518 ! 1590 1519 hbl(:,:) = 0._wp ! here hbl used as a dummy variable, integrating vertically N^2 1591 DO jk = 1, jpkm1 1592 DO jj = 1, jpj ! Mixed layer level: w-level 1593 DO ji = 1, jpi 1594 ikt = mbkt(ji,jj) 1595 hbl(ji,jj) = hbl(ji,jj) + MAX( rn2(ji,jj,jk) , 0._wp ) * e3w_n(ji,jj,jk) 1596 IF( hbl(ji,jj) < zN2_c ) imld_rst(ji,jj) = MIN( jk , ikt ) + 1 ! Mixed layer level 1597 END DO 1598 END DO 1599 END DO 1520 DO_3D( 1, 1, 1, 1, 1, jpkm1 ) 1521 ikt = mbkt(ji,jj) 1522 hbl(ji,jj) = hbl(ji,jj) + MAX( rn2(ji,jj,jk) , 0._wp ) * e3w(ji,jj,jk,Kmm) 1523 IF( hbl(ji,jj) < zN2_c ) imld_rst(ji,jj) = MIN( jk , ikt ) + 1 ! Mixed layer level 1524 END_3D 1600 1525 ! 1601 DO jj = 1, jpj 1602 DO ji = 1, jpi 1603 iiki = imld_rst(ji,jj) 1604 hbl (ji,jj) = gdepw_n(ji,jj,iiki ) * ssmask(ji,jj) ! Turbocline depth 1605 END DO 1606 END DO 1526 DO_2D( 1, 1, 1, 1 ) 1527 iiki = imld_rst(ji,jj) 1528 hbl (ji,jj) = gdepw(ji,jj,iiki ,Kmm) * ssmask(ji,jj) ! Turbocline depth 1529 END_2D 1607 1530 hbl = MAX(hbl,epsln) 1608 1531 hbli(:,:) = hbl(:,:) … … 1612 1535 1613 1536 1614 SUBROUTINE tra_osm( kt )1537 SUBROUTINE tra_osm( kt, Kmm, pts, Krhs ) 1615 1538 !!---------------------------------------------------------------------- 1616 1539 !! *** ROUTINE tra_osm *** … … 1622 1545 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdt, ztrds ! 3D workspace 1623 1546 !!---------------------------------------------------------------------- 1624 INTEGER, INTENT(in) :: kt 1547 INTEGER , INTENT(in) :: kt ! time step index 1548 INTEGER , INTENT(in) :: Kmm, Krhs ! time level indices 1549 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers and RHS of tracer equation 1550 ! 1625 1551 INTEGER :: ji, jj, jk 1626 1552 ! … … 1632 1558 1633 1559 IF( l_trdtra ) THEN !* Save ta and sa trends 1634 ALLOCATE( ztrdt(jpi,jpj,jpk) ) ; ztrdt(:,:,:) = tsa(:,:,:,jp_tem)1635 ALLOCATE( ztrds(jpi,jpj,jpk) ) ; ztrds(:,:,:) = tsa(:,:,:,jp_sal)1560 ALLOCATE( ztrdt(jpi,jpj,jpk) ) ; ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) 1561 ALLOCATE( ztrds(jpi,jpj,jpk) ) ; ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) 1636 1562 ENDIF 1637 1563 1638 1564 ! add non-local temperature and salinity flux 1639 DO jk = 1, jpkm1 1640 DO jj = 2, jpjm1 1641 DO ji = 2, jpim1 1642 tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) & 1643 & - ( ghamt(ji,jj,jk ) & 1644 & - ghamt(ji,jj,jk+1) ) /e3t_n(ji,jj,jk) 1645 tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) & 1646 & - ( ghams(ji,jj,jk ) & 1647 & - ghams(ji,jj,jk+1) ) / e3t_n(ji,jj,jk) 1648 END DO 1649 END DO 1650 END DO 1565 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 1566 pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) & 1567 & - ( ghamt(ji,jj,jk ) & 1568 & - ghamt(ji,jj,jk+1) ) /e3t(ji,jj,jk,Kmm) 1569 pts(ji,jj,jk,jp_sal,Krhs) = pts(ji,jj,jk,jp_sal,Krhs) & 1570 & - ( ghams(ji,jj,jk ) & 1571 & - ghams(ji,jj,jk+1) ) / e3t(ji,jj,jk,Kmm) 1572 END_3D 1651 1573 1652 1574 1653 1575 ! save the non-local tracer flux trends for diagnostic 1654 1576 IF( l_trdtra ) THEN 1655 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)1656 ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)1577 ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:) 1578 ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) - ztrds(:,:,:) 1657 1579 !!bug gm jpttdzdf ==> jpttosm 1658 CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdf, ztrdt )1659 CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdf, ztrds )1580 CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_zdf, ztrdt ) 1581 CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_zdf, ztrds ) 1660 1582 DEALLOCATE( ztrdt ) ; DEALLOCATE( ztrds ) 1661 1583 ENDIF 1662 1584 1663 IF( ln_ctl) THEN1664 CALL prt_ctl( tab3d_1= tsa(:,:,:,jp_tem), clinfo1=' osm - Ta: ', mask1=tmask, &1665 & tab3d_2= tsa(:,:,:,jp_sal), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' )1585 IF(sn_cfctl%l_prtctl) THEN 1586 CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' osm - Ta: ', mask1=tmask, & 1587 & tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) 1666 1588 ENDIF 1667 1589 ! … … 1686 1608 1687 1609 1688 SUBROUTINE dyn_osm( kt )1610 SUBROUTINE dyn_osm( kt, Kmm, puu, pvv, Krhs ) 1689 1611 !!---------------------------------------------------------------------- 1690 1612 !! *** ROUTINE dyn_osm *** … … 1695 1617 !! ** Method : ??? 1696 1618 !!---------------------------------------------------------------------- 1697 INTEGER, INTENT(in) :: kt ! 1619 INTEGER , INTENT( in ) :: kt ! ocean time step index 1620 INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices 1621 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 1698 1622 ! 1699 1623 INTEGER :: ji, jj, jk ! dummy loop indices … … 1707 1631 !code saving tracer trends removed, replace with trdmxl_oce 1708 1632 1709 DO jk = 1, jpkm1 ! add non-local u and v fluxes 1710 DO jj = 2, jpjm1 1711 DO ji = 2, jpim1 1712 ua(ji,jj,jk) = ua(ji,jj,jk) & 1713 & - ( ghamu(ji,jj,jk ) & 1714 & - ghamu(ji,jj,jk+1) ) / e3u_n(ji,jj,jk) 1715 va(ji,jj,jk) = va(ji,jj,jk) & 1716 & - ( ghamv(ji,jj,jk ) & 1717 & - ghamv(ji,jj,jk+1) ) / e3v_n(ji,jj,jk) 1718 END DO 1719 END DO 1720 END DO 1633 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 1634 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) & 1635 & - ( ghamu(ji,jj,jk ) & 1636 & - ghamu(ji,jj,jk+1) ) / e3u(ji,jj,jk,Kmm) 1637 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) & 1638 & - ( ghamv(ji,jj,jk ) & 1639 & - ghamv(ji,jj,jk+1) ) / e3v(ji,jj,jk,Kmm) 1640 END_3D 1721 1641 ! 1722 1642 ! code for saving tracer trends removed
Note: See TracChangeset
for help on using the changeset viewer.
