Changeset 14305
- Timestamp:
- 2021-01-14T19:43:03+01:00 (3 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r14122_ENHANCE-14_smueller_OSMOSIS_streamlining/src/OCE/ZDF/zdfosm.F90
r14280 r14305 113 113 INTEGER :: nn_osm_SD_reduce ! = 0/1/2 flag for getting effective stokes drift from surface value 114 114 LOGICAL :: ln_dia_osm ! Use namelist rn_osm_la 115 LOGICAL :: ln_dia_pyc_scl = .FALSE. ! Output of pycnocline scalar-gradient profiles 116 LOGICAL :: ln_dia_pyc_shr = .FALSE. ! Output of pycnocline velocity-shear profiles 115 117 116 118 … … 215 217 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: p_avm, p_avt ! momentum and tracer Kz (w-points) 216 218 !! 217 INTEGER :: ji, jj, jk ! dummy loop indices219 INTEGER :: ji, jj, jk, jkflt ! dummy loop indices 218 220 219 221 INTEGER :: jl ! dummy loop indices 220 222 221 INTEGER :: ikbot, jkm ax, jkm1, jkp2 !223 INTEGER :: ikbot, jkm1, jkp2 ! 222 224 223 225 REAL(wp) :: ztx, zty, zflageos, zstabl, zbuofdep,zucube ! … … 233 235 REAL(wp) :: zt,zs,zu,zv,zrh ! variables used in constructing averages 234 236 ! Scales 235 REAL(wp) , DIMENSION(jpi,jpj) :: zrad0! Surface solar temperature flux (deg m/s)236 REAL(wp) , DIMENSION(jpi,jpj) :: zradh! Radiative flux at bl base (Buoyancy units)237 REAL(wp) , DIMENSION(jpi,jpj) :: zradav! Radiative flux, bl average (Buoyancy Units)237 REAL(wp) :: zrad0 ! Surface solar temperature flux (deg m/s) 238 REAL(wp) :: zradh ! Radiative flux at bl base (Buoyancy units) 239 REAL(wp) :: zradav ! Radiative flux, bl average (Buoyancy Units) 238 240 REAL(wp), DIMENSION(jpi,jpj) :: zustar ! friction velocity 239 241 REAL(wp), DIMENSION(jpi,jpj) :: zwstrl ! Langmuir velocity scale … … 241 243 REAL(wp), DIMENSION(jpi,jpj) :: zwstrc ! Convective velocity scale 242 244 REAL(wp), DIMENSION(jpi,jpj) :: zuw0 ! Surface u-momentum flux 243 REAL(wp) , DIMENSION(jpi,jpj) :: zvw0! Surface v-momentum flux245 REAL(wp) :: zvw0 ! Surface v-momentum flux 244 246 REAL(wp), DIMENSION(jpi,jpj) :: zwth0 ! Surface heat flux (Kinematic) 245 247 REAL(wp), DIMENSION(jpi,jpj) :: zws0 ! Surface freshwater flux … … 299 301 REAL(wp) :: zwth_ent,zws_ent ! heat and salinity fluxes at the top of the pycnocline 300 302 REAL(wp) :: zuw_bse,zvw_bse ! momentum fluxes at the top of the pycnocline 301 REAL(wp) , DIMENSION(jpi,jpj,jpk) :: zdtdz_pyc ! parametrized gradient of temperature in pycnocline302 REAL(wp) , DIMENSION(jpi,jpj,jpk) :: zdsdz_pyc ! parametrised gradient of salinity in pycnocline303 REAL(wp) :: zdtdz_pyc ! Parametrized gradient of temperature in pycnocline 304 REAL(wp) :: zdsdz_pyc ! Parametrised gradient of salinity in pycnocline 303 305 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdbdz_pyc ! parametrised gradient of buoyancy in the pycnocline 304 REAL(wp) , DIMENSION(jpi,jpj,jpk) ::zdudz_pyc ! u-shear across the pycnocline305 REAL(wp) , DIMENSION(jpi,jpj,jpk) ::zdvdz_pyc ! v-shear across the pycnocline306 REAL(wp) :: zdudz_pyc ! u-shear across the pycnocline 307 REAL(wp) :: zdvdz_pyc ! v-shear across the pycnocline 306 308 REAL(wp), DIMENSION(jpi,jpj) :: zdbds_mle ! Magnitude of horizontal buoyancy gradient. 307 309 ! Flux-gradient relationship variables … … 317 319 318 320 ! For calculating Ri#-dependent mixing 319 REAL(wp), DIMENSION(jpi,jpj ,jpk) :: z3du! u-shear^2320 REAL(wp), DIMENSION(jpi,jpj ,jpk) :: z3dv! v-shear^2321 REAL(wp) , DIMENSION(jpi,jpj,jpk) :: zrimix ! spatial form of ri#-induced diffusion321 REAL(wp), DIMENSION(jpi,jpj) :: z2du ! u-shear^2 322 REAL(wp), DIMENSION(jpi,jpj) :: z2dv ! v-shear^2 323 REAL(wp) :: zrimix ! Spatial form of ri#-induced diffusion 322 324 323 325 ! Temporary variables … … 342 344 REAL(wp) :: zsqrtpi, z_two_thirds, zproportion, ztransp, zthickness 343 345 REAL(wp) :: z2k_times_thickness, zsqrt_depth, zexp_depth, zdstokes0, zf, zexperfc 346 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: z3ddz_pyc_1, z3ddz_pyc_2 ! Pycnocline gradient/shear profiles 347 INTEGER :: istat ! Memory allocation status 344 348 345 349 ! For debugging … … 349 353 IF( ln_timing ) CALL timing_start('zdf_osm') 350 354 ibld(:,:) = 0 ; imld(:,:) = 0 351 z rad0(:,:) = 0._wp ; zradh(:,:) = 0._wp ; zradav(:,:) = 0._wp ; zustar(:,:) = 0._wp355 zustar(:,:) = 0.0_wp 352 356 zwstrl(:,:) = 0._wp ; zvstr(:,:) = 0._wp ; zwstrc(:,:) = 0._wp ; zuw0(:,:) = 0._wp 353 z vw0(:,:) = 0._wp ; zwth0(:,:) = 0._wp ; zws0(:,:) = 0._wp ; zwb0(:,:) = 0._wp357 zwth0(:,:) = 0.0_wp ; zws0(:,:) = 0.0_wp ; zwb0(:,:) = 0.0_wp 354 358 zwthav(:,:) = 0._wp ; zwsav(:,:) = 0._wp ; zwbav(:,:) = 0._wp ; zwb_ent(:,:) = 0._wp 355 359 zustke(:,:) = 0._wp ; zla(:,:) = 0._wp ; zcos_wind(:,:) = 0._wp ; zsin_wind(:,:) = 0._wp … … 370 374 zwth_ent = 0._wp ; zws_ent = 0._wp 371 375 ! 372 zdtdz_pyc(:,:,:) = 0._wp ; zdsdz_pyc(:,:,:) = 0._wp ; zdbdz_pyc(:,:,:) = 0._wp 373 zdudz_pyc(:,:,:) = 0._wp ; zdvdz_pyc(:,:,:) = 0._wp 376 zdbdz_pyc(:,:,:) = 0.0_wp 374 377 ! 375 378 zdtdz_bl_ext(:,:) = 0._wp ; zdsdz_bl_ext(:,:) = 0._wp ; zdbdz_bl_ext(:,:) = 0._wp … … 396 399 ! Calculate boundary layer scales 397 400 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 398 399 ! Assume two-band radiation model for depth of OSBL 400 zz0 = rn_abs ! surface equi-partition in 2-bands 401 zz1 = 1. - rn_abs 401 ! 402 ! Turbulent surface fluxes and fluxes averaged over depth of the OSBL 403 zz0 = rn_abs ! Assume two-band radiation model for depth of OSBL - surface equi-partition in 2-bands 404 zz1 = 1.0_wp - rn_abs 405 DO_2D( 0, 0, 0, 0 ) 406 zrad0 = qsr(ji,jj) * r1_rho0_rcp ! Surface downward irradiance (so always +ve) 407 zradh = zrad0 * & ! Downwards irradiance at base of boundary layer 408 & ( zz0 * EXP( -1.0_wp * hbl(ji,jj) / rn_si0 ) + zz1 * EXP( -1.0_wp * hbl(ji,jj) / rn_si1 ) ) 409 zradav = zrad0 * & ! Downwards irradiance averaged over depth of the OSBL 410 & ( zz0 * ( 1.0_wp - EXP( -hbl(ji,jj)/rn_si0 ) ) * rn_si0 + & 411 & zz1 * ( 1.0_wp - EXP( -hbl(ji,jj)/rn_si1 ) ) * rn_si1 ) / hbl(ji,jj) 412 zwth0(ji,jj) = - qns(ji,jj) * r1_rho0_rcp * tmask(ji,jj,1) ! Upwards surface Temperature flux for non-local term 413 zwthav(ji,jj) = 0.5_wp * zwth0(ji,jj) - & ! Turbulent heat flux averaged over depth of OSBL 414 & ( 0.5_wp * ( zrad0 + zradh ) - zradav ) 415 END_2D 402 416 DO_2D( 0, 0, 0, 0 ) 403 ! Surface downward irradiance (so always +ve) 404 zrad0(ji,jj) = qsr(ji,jj) * r1_rho0_rcp 405 ! Downwards irradiance at base of boundary layer 406 zradh(ji,jj) = zrad0(ji,jj) * ( zz0 * EXP( -hbl(ji,jj)/rn_si0 ) + zz1 * EXP( -hbl(ji,jj)/rn_si1) ) 407 ! Downwards irradiance averaged over depth of the OSBL 408 zradav(ji,jj) = zrad0(ji,jj) * ( zz0 * ( 1.0 - EXP( -hbl(ji,jj)/rn_si0 ) )*rn_si0 & 409 & + zz1 * ( 1.0 - EXP( -hbl(ji,jj)/rn_si1 ) )*rn_si1 ) / hbl(ji,jj) 417 zws0(ji,jj) = -1.0_wp * & ! Upwards surface salinity flux for non-local term 418 & ( ( emp(ji,jj) - rnf(ji,jj) ) * ts(ji,jj,1,jp_sal,Kmm) + sfx(ji,jj) ) * r1_rho0 * tmask(ji,jj,1) 419 zthermal = rab_n(ji,jj,1,jp_tem) 420 zbeta = rab_n(ji,jj,1,jp_sal) 421 zwb0(ji,jj) = grav * zthermal * zwth0(ji,jj) - & ! Non radiative upwards surface buoyancy flux 422 & grav * zbeta * zws0(ji,jj) 423 zwsav(ji,jj) = 0.5 * zws0(ji,jj) ! Turbulent salinity flux averaged over depth of the OBSL 424 zwbav(ji,jj) = grav * zthermal * zwthav(ji,jj) - & ! Turbulent buoyancy flux averaged over the depth of the 425 & grav * zbeta * zwsav(ji,jj) ! OBSBL 410 426 END_2D 411 ! Turbulent surface fluxes and fluxes averaged over depth of the OSBL412 427 DO_2D( 0, 0, 0, 0 ) 413 zthermal = rab_n(ji,jj,1,jp_tem) 414 zbeta = rab_n(ji,jj,1,jp_sal) 415 ! Upwards surface Temperature flux for non-local term 416 zwth0(ji,jj) = - qns(ji,jj) * r1_rho0_rcp * tmask(ji,jj,1) 417 ! Upwards surface salinity flux for non-local term 418 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) 419 ! Non radiative upwards surface buoyancy flux 420 zwb0(ji,jj) = grav * zthermal * zwth0(ji,jj) - grav * zbeta * zws0(ji,jj) 421 ! turbulent heat flux averaged over depth of OSBL 422 zwthav(ji,jj) = 0.5 * zwth0(ji,jj) - ( 0.5*( zrad0(ji,jj) + zradh(ji,jj) ) - zradav(ji,jj) ) 423 ! turbulent salinity flux averaged over depth of the OBSL 424 zwsav(ji,jj) = 0.5 * zws0(ji,jj) 425 ! turbulent buoyancy flux averaged over the depth of the OBSBL 426 zwbav(ji,jj) = grav * zthermal * zwthav(ji,jj) - grav * zbeta * zwsav(ji,jj) 427 ! Surface upward velocity fluxes 428 zuw0(ji,jj) = - 0.5 * (utau(ji-1,jj) + utau(ji,jj)) * r1_rho0 * tmask(ji,jj,1) 429 zvw0(ji,jj) = - 0.5 * (vtau(ji,jj-1) + vtau(ji,jj)) * r1_rho0 * tmask(ji,jj,1) 430 ! Friction velocity (zustar), at T-point : LMD94 eq. 2 431 zustar(ji,jj) = MAX( SQRT( SQRT( zuw0(ji,jj) * zuw0(ji,jj) + zvw0(ji,jj) * zvw0(ji,jj) ) ), 1.0e-8 ) 428 zuw0(ji,jj) = - 0.5 * (utau(ji-1,jj) + utau(ji,jj)) * & ! Surface upward velocity fluxes 429 & r1_rho0 * tmask(ji,jj,1) 430 zvw0 = - 0.5 * (vtau(ji,jj-1) + vtau(ji,jj)) * r1_rho0 * tmask(ji,jj,1) 431 zustar(ji,jj) = MAX( SQRT( SQRT( zuw0(ji,jj) * & ! Friction velocity (zustar), at T-point : LMD94 eq. 2 432 & zuw0(ji,jj) + zvw0 * zvw0 ) ), 1.0e-8_wp ) 432 433 zcos_wind(ji,jj) = -zuw0(ji,jj) / ( zustar(ji,jj) * zustar(ji,jj) ) 433 zsin_wind(ji,jj) = -zvw0 (ji,jj)/ ( zustar(ji,jj) * zustar(ji,jj) )434 zsin_wind(ji,jj) = -zvw0 / ( zustar(ji,jj) * zustar(ji,jj) ) 434 435 END_2D 435 436 ! Calculate Stokes drift in direction of wind (zustke) and Stokes penetration depth (dstokes) … … 699 700 700 701 CALL zdf_osm_external_gradients( ibld+2, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext ) 701 CALL zdf_osm_pycnocline_scalar_profiles( zdtdz_pyc, zdsdz_pyc, zdbdz_pyc, zalpha_pyc ) 702 CALL zdf_osm_pycnocline_shear_profiles( zdudz_pyc, zdvdz_pyc ) 702 CALL zdf_osm_pycnocline_buoyancy_profiles( zdbdz_pyc, zalpha_pyc ) 703 703 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 704 704 ! Eddy viscosity/diffusivity and non-gradient terms in the flux-gradient relationship … … 1028 1028 END_2D 1029 1029 1030 ! pynocline contributions 1031 DO_2D( 0, 0, 0, 0 ) 1032 IF ( .not. lconv(ji,jj) ) THEN 1033 IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN 1034 DO jk= 2, ibld(ji,jj) 1035 znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 1036 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk) 1037 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk) 1038 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc(ji,jj,jk) 1039 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc(ji,jj,jk) 1040 END DO 1041 END IF 1030 ! Pynocline contributions 1031 ! 1032 IF ( ln_dia_pyc_scl .OR. ln_dia_pyc_shr ) THEN ! Allocate arrays used for output of pycnocline gradient/shear profiles 1033 ALLOCATE( z3ddz_pyc_1(jpi,jpj,jpk), z3ddz_pyc_2(jpi,jpj,jpk), STAT=istat ) 1034 CALL mpp_sum( 'zdfosm', istat ) 1035 IF ( istat /= 0 ) CALL ctl_stop( 'zdf_osm: failed to allocate temporary arrays' ) 1036 z3ddz_pyc_1(:,:,:) = 0.0_wp 1037 z3ddz_pyc_2(:,:,:) = 0.0_wp 1038 END IF 1039 DO_2D( 0, 0, 0, 0 ) 1040 IF ( lconv (ji,jj) ) THEN 1041 ! Unstable conditions. Shouldn;t be needed with no pycnocline code. 1042 ! zugrad = 0.7 * zdu_ml(ji,jj) / zdh(ji,jj) + 0.3 * zustar(ji,jj)*zustar(ji,jj) / & 1043 ! & ( ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) * & 1044 ! & MIN(zla(ji,jj)**(8.0/3.0) + epsln, 0.12 )) 1045 !Alan is this right? 1046 ! zvgrad = ( 0.7 * zdv_ml(ji,jj) + & 1047 ! & 2.0 * ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) / & 1048 ! & ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird + epsln ) & 1049 ! & )/ (zdh(ji,jj) + epsln ) 1050 ! DO jk = 2, ibld(ji,jj) - 1 + ibld_ext 1051 ! znd = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / (zdh(ji,jj) + epsln ) - zzeta_v 1052 ! IF ( znd <= 0.0 ) THEN 1053 ! zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( 3.0 * znd ) 1054 ! zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( 3.0 * znd ) 1055 ! ELSE 1056 ! zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( -2.0 * znd ) 1057 ! zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( -2.0 * znd ) 1058 ! ENDIF 1059 ! END DO 1060 ELSE ! Stable conditions 1061 IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN 1062 ! If pycnocline profile only defined when depth steady of increasing. 1063 IF ( zdhdt(ji,jj) > 0.0_wp ) THEN ! Depth increasing, or steady. 1064 IF ( zdb_bl(ji,jj) > 0.0_wp ) THEN 1065 IF ( zhol(ji,jj) >= 0.5_wp ) THEN ! Very stable - 'thick' pycnocline 1066 ztmp = 1.0_wp / MAX( zhbl(ji,jj), epsln ) 1067 ztgrad = zdt_bl(ji,jj) * ztmp 1068 zsgrad = zds_bl(ji,jj) * ztmp 1069 zbgrad = zdb_bl(ji,jj) * ztmp 1070 DO jk = 2, ibld(ji,jj) 1071 znd = gdepw(ji,jj,jk,Kmm) * ztmp 1072 zdtdz_pyc = ztgrad * EXP( -15.0 * ( znd - 0.9_wp )**2 ) 1073 zdsdz_pyc = zsgrad * EXP( -15.0 * ( znd - 0.9_wp )**2 ) 1074 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc 1075 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc 1076 IF ( ln_dia_pyc_scl ) THEN 1077 z3ddz_pyc_1(ji,jj,jk) = zdtdz_pyc 1078 z3ddz_pyc_2(ji,jj,jk) = zdsdz_pyc 1079 END IF 1080 END DO 1081 ELSE ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form. 1082 ztmp = 1.0_wp / MAX( zdh(ji,jj), epsln ) 1083 ztgrad = zdt_bl(ji,jj) * ztmp 1084 zsgrad = zds_bl(ji,jj) * ztmp 1085 zbgrad = zdb_bl(ji,jj) * ztmp 1086 DO jk = 2, ibld(ji,jj) 1087 znd = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) * ztmp 1088 zdtdz_pyc = ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 1089 zdsdz_pyc = zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 1090 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc 1091 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc 1092 IF ( ln_dia_pyc_scl ) THEN 1093 z3ddz_pyc_1(ji,jj,jk) = zdtdz_pyc 1094 z3ddz_pyc_2(ji,jj,jk) = zdsdz_pyc 1095 END IF 1096 END DO 1097 ENDIF ! IF (zhol >=0.5) 1098 ENDIF ! IF (zdb_bl> 0.) 1099 ENDIF ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero 1100 END IF 1042 1101 END IF 1043 END_2D 1102 END_2D 1103 IF ( ln_dia_pyc_scl ) THEN ! Output of pycnocline gradient profiles 1104 IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask(:,:,:) * z3ddz_pyc_1(:,:,:) ) 1105 IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask(:,:,:) * z3ddz_pyc_2(:,:,:) ) 1106 IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask(:,:,:) * zdbdz_pyc(:,:,:) ) 1107 END IF 1108 DO_2D( 0, 0, 0, 0 ) 1109 IF ( .NOT. lconv (ji,jj) ) THEN 1110 IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN 1111 zugrad = 3.25_wp * zdu_bl(ji,jj) / zhbl(ji,jj) 1112 zvgrad = 2.75_wp * zdv_bl(ji,jj) / zhbl(ji,jj) 1113 DO jk = 2, ibld(ji,jj) 1114 znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 1115 IF ( znd < 1.0 ) THEN 1116 zdudz_pyc = zugrad * EXP( -40.0_wp * ( znd - 1.0_wp )**2 ) 1117 ELSE 1118 zdudz_pyc = zugrad * EXP( -20.0_wp * ( znd - 1.0_wp )**2 ) 1119 ENDIF 1120 zdvdz_pyc = zvgrad * EXP( -20.0_wp * ( znd - 0.85_wp )**2 ) 1121 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc 1122 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc 1123 IF ( ln_dia_pyc_shr ) THEN 1124 z3ddz_pyc_1(ji,jj,jk) = zdudz_pyc 1125 z3ddz_pyc_2(ji,jj,jk) = zdvdz_pyc 1126 END IF 1127 END DO 1128 END IF 1129 END IF 1130 END_2D 1131 IF ( ln_dia_pyc_shr ) THEN ! Output of pycnocline shear profiles 1132 IF ( iom_use("dudz_pyc") ) CALL iom_put( "zdudz_pyc", wmask(:,:,:) * z3ddz_pyc_1(:,:,:) ) 1133 IF ( iom_use("dvdz_pyc") ) CALL iom_put( "zdvdz_pyc", wmask(:,:,:) * z3ddz_pyc_2(:,:,:) ) 1134 END IF 1044 1135 IF(ln_dia_osm) THEN 1045 IF ( iom_use("ghamu_b") ) CALL iom_put( "ghamu_b", wmask*ghamu ) 1046 IF ( iom_use("ghamv_b") ) CALL iom_put( "ghamv_b", wmask*ghamv ) 1047 END IF 1136 IF ( iom_use("ghamu_b") ) CALL iom_put( "ghamu_b", wmask*ghamu ) 1137 IF ( iom_use("ghamv_b") ) CALL iom_put( "ghamv_b", wmask*ghamv ) 1138 END IF 1139 IF ( ln_dia_pyc_scl .OR. ln_dia_pyc_shr ) THEN ! Deallocate arrays used for output of pycnocline gradient/shear profiles 1140 DEALLOCATE( z3ddz_pyc_1, z3ddz_pyc_2 ) 1141 END IF 1048 1142 1049 1143 DO_2D( 0, 0, 0, 0 ) … … 1057 1151 IF ( iom_use("ghamu_1") ) CALL iom_put( "ghamu_1", wmask*ghamu ) 1058 1152 IF ( iom_use("ghamv_1") ) CALL iom_put( "ghamv_1", wmask*ghamv ) 1059 IF ( iom_use("zdudz_pyc") ) CALL iom_put( "zdudz_pyc", wmask*zdudz_pyc )1060 IF ( iom_use("zdvdz_pyc") ) CALL iom_put( "zdvdz_pyc", wmask*zdvdz_pyc )1061 1153 IF ( iom_use("zviscos") ) CALL iom_put( "zviscos", wmask*zviscos ) 1062 1154 END IF 1155 1063 1156 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 1064 1157 ! Need to put in code for contributions that are applied explicitly to … … 1080 1173 END_2D 1081 1174 1082 IF(ln_dia_osm) THEN 1083 IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask*zdtdz_pyc ) 1084 IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask*zdsdz_pyc ) 1085 IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask*zdbdz_pyc ) 1086 END IF 1087 1088 ! KPP-style Ri# mixing 1089 IF( ln_kpprimix) THEN 1090 DO_3D( 1, 0, 1, 0, 2, jpkm1 ) !* Shear production at uw- and vw-points (energy conserving form) 1091 z3du(ji,jj,jk) = 0.5 * ( uu(ji,jj,jk-1,Kmm) - uu(ji ,jj,jk,Kmm) ) & 1092 & * ( uu(ji,jj,jk-1,Kbb) - uu(ji ,jj,jk,Kbb) ) * wumask(ji,jj,jk) & 1093 & / ( e3uw(ji,jj,jk,Kmm) * e3uw(ji,jj,jk,Kbb) ) 1094 z3dv(ji,jj,jk) = 0.5 * ( vv(ji,jj,jk-1,Kmm) - vv(ji,jj ,jk,Kmm) ) & 1095 & * ( vv(ji,jj,jk-1,Kbb) - vv(ji,jj ,jk,Kbb) ) * wvmask(ji,jj,jk) & 1096 & / ( e3vw(ji,jj,jk,Kmm) * e3vw(ji,jj,jk,Kbb) ) 1097 END_3D 1098 ! 1099 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 1100 ! ! shear prod. at w-point weightened by mask 1101 zesh2 = ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) & 1102 & + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) ) 1103 ! ! local Richardson number 1104 zri = MAX( rn2b(ji,jj,jk), 0._wp ) / MAX(zesh2, epsln) 1105 zfri = MIN( zri / rn_riinfty , 1.0_wp ) 1106 zfri = ( 1.0_wp - zfri * zfri ) 1107 zrimix(ji,jj,jk) = zfri * zfri * zfri * wmask(ji, jj, jk) 1108 END_3D 1109 1175 ! KPP-style Ri# mixing 1176 IF ( ln_kpprimix ) THEN 1177 jkflt = jpk 1110 1178 DO_2D( 0, 0, 0, 0 ) 1111 DO jk = ibld(ji,jj) + 1, jpkm1 1112 zdiffut(ji,jj,jk) = zrimix(ji,jj,jk)*rn_difri 1113 zviscos(ji,jj,jk) = zrimix(ji,jj,jk)*rn_difri 1114 END DO 1179 IF ( ibld(ji,jj) < jkflt ) jkflt = ibld(ji,jj) 1115 1180 END_2D 1116 1181 DO jk = jkflt+1, jpkm1 1182 ! Shear production at uw- and vw-points (energy conserving form) 1183 DO_2D( 1, 0, 1, 0 ) 1184 IF ( jk > MIN( ibld(ji,jj), ibld(ji+1,jj) ) ) THEN 1185 z2du(ji,jj) = 0.5_wp * ( uu(ji,jj,jk-1,Kmm) - uu(ji,jj,jk,Kmm) ) * & 1186 & ( uu(ji,jj,jk-1,Kbb) - uu(ji,jj,jk,Kbb) ) * wumask(ji,jj,jk) / & 1187 & ( e3uw(ji,jj,jk,Kmm) * e3uw(ji,jj,jk,Kbb) ) 1188 END IF 1189 IF ( jk > MIN( ibld(ji,jj), ibld(ji,jj+1) ) ) THEN 1190 z2dv(ji,jj) = 0.5_wp * ( 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 IF 1194 END_2D 1195 DO_2D( 0, 0, 0, 0 ) 1196 IF ( jk > ibld(ji,jj) ) THEN 1197 ! Shear prod. at w-point weightened by mask 1198 zesh2 = ( z2du(ji-1,jj) + z2du(ji,jj) ) / MAX( 1.0_wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) & 1199 & + ( z2dv(ji,jj-1) + z2dv(ji,jj) ) / MAX( 1.0_wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) ) 1200 ! Local Richardson number 1201 zri = MAX( rn2b(ji,jj,jk), 0.0_wp ) / MAX(zesh2, epsln) 1202 zfri = MIN( zri / rn_riinfty , 1.0_wp ) 1203 zfri = ( 1.0_wp - zfri * zfri ) 1204 zrimix = zfri * zfri * zfri * wmask(ji, jj, jk) 1205 zdiffut(ji,jj,jk) = zrimix*rn_difri 1206 zviscos(ji,jj,jk) = zrimix*rn_difri 1207 END IF 1208 END_2D 1209 END DO 1117 1210 END IF ! ln_kpprimix = .true. 1118 1211 … … 1157 1250 END_2D 1158 1251 ENDIF 1159 1160 IF(ln_dia_osm) THEN1161 IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask*zdtdz_pyc )1162 IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask*zdsdz_pyc )1163 IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask*zdbdz_pyc )1164 END IF1165 1166 1252 1167 1253 ! Lateral boundary conditions on zvicos (sign unchanged), needed to caclulate viscosities on u and v grids … … 1727 1813 END SUBROUTINE zdf_osm_external_gradients 1728 1814 1729 SUBROUTINE zdf_osm_pycnocline_scalar_profiles( zdtdz, zdsdz, zdbdz, zalpha ) 1730 1731 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdtdz, zdsdz, zdbdz ! gradients in the pycnocline 1732 REAL(wp), DIMENSION(jpi,jpj) :: zalpha 1733 1734 INTEGER :: jk, jj, ji 1735 REAL(wp) :: ztgrad, zsgrad, zbgrad 1736 REAL(wp) :: zgamma_b_nd, znd 1737 REAL(wp) :: zzeta_m, zzeta_en, zbuoy_pyc_sc 1738 REAL(wp), PARAMETER :: zgamma_b = 2.25, zzeta_sh = 0.15 1739 1740 IF( ln_timing ) CALL timing_start('zdf_osm_pscp') 1741 DO_2D( 0, 0, 0, 0 ) 1742 IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN 1743 IF ( lconv(ji,jj) ) THEN ! convective conditions 1744 IF ( lpyc(ji,jj) ) THEN 1745 zzeta_m = 0.1 + 0.3 / ( 1.0 + EXP( -3.5 * LOG10( -zhol(ji,jj) ) ) ) 1746 zalpha(ji,jj) = 2.0 * ( 1.0 - ( 0.80 * zzeta_m + 0.5 * SQRT( 3.14159 / zgamma_b ) ) * zdbdz_bl_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj) ) / ( 0.723 + SQRT( 3.14159 / zgamma_b ) ) 1747 zalpha(ji,jj) = MAX( zalpha(ji,jj), 0._wp ) 1748 1749 ztmp = 1._wp/MAX(zdh(ji,jj), epsln) 1750 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1751 ! Commented lines in this section are not needed in new code, once tested ! 1752 ! can be removed ! 1753 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1754 ! ztgrad = zalpha * zdt_ml(ji,jj) * ztmp + zdtdz_bl_ext(ji,jj) 1755 ! zsgrad = zalpha * zds_ml(ji,jj) * ztmp + zdsdz_bl_ext(ji,jj) 1756 zbgrad = zalpha(ji,jj) * zdb_ml(ji,jj) * ztmp + zdbdz_bl_ext(ji,jj) 1757 zgamma_b_nd = zdbdz_bl_ext(ji,jj) * zdh(ji,jj) / MAX(zdb_ml(ji,jj), epsln) 1758 DO jk = 2, ibld(ji,jj)+ibld_ext 1759 znd = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) * ztmp 1760 IF ( znd <= zzeta_m ) THEN 1761 ! zdtdz(ji,jj,jk) = zdtdz_bl_ext(ji,jj) + zalpha * zdt_ml(ji,jj) * ztmp * & 1762 ! & EXP( -6.0 * ( znd -zzeta_m )**2 ) 1763 ! zdsdz(ji,jj,jk) = zdsdz_bl_ext(ji,jj) + zalpha * zds_ml(ji,jj) * ztmp * & 1764 ! & EXP( -6.0 * ( znd -zzeta_m )**2 ) 1765 zdbdz(ji,jj,jk) = zdbdz_bl_ext(ji,jj) + zalpha(ji,jj) * zdb_ml(ji,jj) * ztmp * & 1766 & EXP( -6.0 * ( znd -zzeta_m )**2 ) 1767 ELSE 1768 ! zdtdz(ji,jj,jk) = ztgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 ) 1769 ! zdsdz(ji,jj,jk) = zsgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 ) 1770 zdbdz(ji,jj,jk) = zbgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 ) 1771 ENDIF 1772 END DO 1773 ENDIF ! if no pycnocline pycnocline gradients set to zero 1774 ELSE 1775 ! stable conditions 1776 ! if pycnocline profile only defined when depth steady of increasing. 1777 IF ( zdhdt(ji,jj) > 0.0 ) THEN ! Depth increasing, or steady. 1778 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 1779 IF ( zhol(ji,jj) >= 0.5 ) THEN ! Very stable - 'thick' pycnocline 1780 ztmp = 1._wp/MAX(zhbl(ji,jj), epsln) 1781 ztgrad = zdt_bl(ji,jj) * ztmp 1782 zsgrad = zds_bl(ji,jj) * ztmp 1783 zbgrad = zdb_bl(ji,jj) * ztmp 1784 DO jk = 2, ibld(ji,jj) 1785 znd = gdepw(ji,jj,jk,Kmm) * ztmp 1786 zdtdz(ji,jj,jk) = ztgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 1787 zdbdz(ji,jj,jk) = zbgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 1788 zdsdz(ji,jj,jk) = zsgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 1789 END DO 1790 ELSE ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form. 1791 ztmp = 1._wp/MAX(zdh(ji,jj), epsln) 1792 ztgrad = zdt_bl(ji,jj) * ztmp 1793 zsgrad = zds_bl(ji,jj) * ztmp 1794 zbgrad = zdb_bl(ji,jj) * ztmp 1795 DO jk = 2, ibld(ji,jj) 1796 znd = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) * ztmp 1797 zdtdz(ji,jj,jk) = ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 1798 zdbdz(ji,jj,jk) = zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 1799 zdsdz(ji,jj,jk) = zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 1800 END DO 1801 ENDIF ! IF (zhol >=0.5) 1802 ENDIF ! IF (zdb_bl> 0.) 1803 ENDIF ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero 1804 ENDIF ! IF (lconv) 1805 ENDIF ! IF ( ibld(ji,jj) < mbkt(ji,jj) ) 1806 END_2D 1807 IF( ln_timing ) CALL timing_stop('zdf_osm_pscp') 1808 1809 END SUBROUTINE zdf_osm_pycnocline_scalar_profiles 1810 1811 SUBROUTINE zdf_osm_pycnocline_shear_profiles( zdudz, zdvdz ) 1812 !!--------------------------------------------------------------------- 1813 !! *** ROUTINE zdf_osm_pycnocline_shear_profiles *** 1814 !! 1815 !! ** Purpose : Calculates velocity shear in the pycnocline 1816 !! 1817 !! ** Method : 1818 !! 1819 !!---------------------------------------------------------------------- 1820 1821 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdudz, zdvdz 1822 1823 INTEGER :: jk, jj, ji 1824 REAL(wp) :: zugrad, zvgrad, znd 1825 REAL(wp) :: zzeta_v = 0.45 1815 SUBROUTINE zdf_osm_pycnocline_buoyancy_profiles( pdbdz, palpha ) 1816 REAL(wp), DIMENSION(:,:,:), INTENT( inout ) :: pdbdz ! Gradients in the pycnocline 1817 REAL(wp), DIMENSION(:,:), INTENT( inout ) :: palpha 1818 INTEGER :: jk, jj, ji 1819 REAL(wp) :: zbgrad 1820 REAL(wp) :: zgamma_b_nd, znd 1821 REAL(wp) :: zzeta_m 1822 REAL(wp), PARAMETER :: ppgamma_b = 2.25_wp 1826 1823 ! 1827 IF( ln_timing ) CALL timing_start('zdf_osm_pshp') 1824 IF( ln_timing ) CALL timing_start('zdf_osm_pscp') 1825 ! 1828 1826 DO_2D( 0, 0, 0, 0 ) 1829 !1830 1827 IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN 1831 IF ( lconv (ji,jj) ) THEN 1832 ! Unstable conditions. Shouldn;t be needed with no pycnocline code. 1833 ! zugrad = 0.7 * zdu_ml(ji,jj) / zdh(ji,jj) + 0.3 * zustar(ji,jj)*zustar(ji,jj) / & 1834 ! & ( ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) * & 1835 ! & MIN(zla(ji,jj)**(8.0/3.0) + epsln, 0.12 )) 1836 !Alan is this right? 1837 ! zvgrad = ( 0.7 * zdv_ml(ji,jj) + & 1838 ! & 2.0 * ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) / & 1839 ! & ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird + epsln ) & 1840 ! & )/ (zdh(ji,jj) + epsln ) 1841 ! DO jk = 2, ibld(ji,jj) - 1 + ibld_ext 1842 ! znd = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / (zdh(ji,jj) + epsln ) - zzeta_v 1843 ! IF ( znd <= 0.0 ) THEN 1844 ! zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( 3.0 * znd ) 1845 ! zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( 3.0 * znd ) 1846 ! ELSE 1847 ! zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( -2.0 * znd ) 1848 ! zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( -2.0 * znd ) 1849 ! ENDIF 1850 ! END DO 1851 ELSE 1852 ! stable conditions 1853 zugrad = 3.25 * zdu_bl(ji,jj) / zhbl(ji,jj) 1854 zvgrad = 2.75 * zdv_bl(ji,jj) / zhbl(ji,jj) 1855 DO jk = 2, ibld(ji,jj) 1856 znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 1857 IF ( znd < 1.0 ) THEN 1858 zdudz(ji,jj,jk) = zugrad * EXP( -40.0 * ( znd - 1.0 )**2 ) 1859 ELSE 1860 zdudz(ji,jj,jk) = zugrad * EXP( -20.0 * ( znd - 1.0 )**2 ) 1861 ENDIF 1862 zdvdz(ji,jj,jk) = zvgrad * EXP( -20.0 * ( znd - 0.85 )**2 ) 1863 END DO 1864 ENDIF 1865 ! 1866 END IF ! IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) 1828 IF ( lconv(ji,jj) ) THEN ! convective conditions 1829 IF ( lpyc(ji,jj) ) THEN 1830 zzeta_m = 0.1_wp + 0.3_wp / ( 1.0_wp + EXP( -3.5_wp * LOG10( -1.0_wp * zhol(ji,jj) ) ) ) 1831 palpha(ji,jj) = 2.0_wp * ( 1.0_wp - ( 0.80_wp * zzeta_m + 0.5_wp * SQRT( 3.14159_wp / ppgamma_b ) ) * & 1832 & zdbdz_bl_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj) ) / & 1833 & ( 0.723_wp + SQRT( 3.14159_wp / ppgamma_b ) ) 1834 palpha(ji,jj) = MAX( palpha(ji,jj), 0.0_wp ) 1835 ztmp = 1.0_wp / MAX( zdh(ji,jj), epsln ) 1836 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1837 ! Commented lines in this section are not needed in new code, once tested ! 1838 ! can be removed ! 1839 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1840 ! ztgrad = zalpha * zdt_ml(ji,jj) * ztmp + zdtdz_bl_ext(ji,jj) 1841 ! zsgrad = zalpha * zds_ml(ji,jj) * ztmp + zdsdz_bl_ext(ji,jj) 1842 zbgrad = palpha(ji,jj) * zdb_ml(ji,jj) * ztmp + zdbdz_bl_ext(ji,jj) 1843 zgamma_b_nd = zdbdz_bl_ext(ji,jj) * zdh(ji,jj) / MAX(zdb_ml(ji,jj), epsln) 1844 DO jk = 2, ibld(ji,jj)+ibld_ext 1845 znd = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) * ztmp 1846 IF ( znd <= zzeta_m ) THEN 1847 ! zdtdz(ji,jj,jk) = zdtdz_bl_ext(ji,jj) + zalpha * zdt_ml(ji,jj) * ztmp * & 1848 ! & EXP( -6.0 * ( znd -zzeta_m )**2 ) 1849 ! zdsdz(ji,jj,jk) = zdsdz_bl_ext(ji,jj) + zalpha * zds_ml(ji,jj) * ztmp * & 1850 ! & EXP( -6.0 * ( znd -zzeta_m )**2 ) 1851 pdbdz(ji,jj,jk) = zdbdz_bl_ext(ji,jj) + palpha(ji,jj) * zdb_ml(ji,jj) * ztmp * & 1852 & EXP( -6.0_wp * ( znd -zzeta_m )**2 ) 1853 ELSE 1854 ! zdtdz(ji,jj,jk) = ztgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 ) 1855 ! zdsdz(ji,jj,jk) = zsgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 ) 1856 pdbdz(ji,jj,jk) = zbgrad * EXP( -1.0_wp * ppgamma_b * ( znd - zzeta_m )**2 ) 1857 ENDIF 1858 END DO 1859 ENDIF ! If no pycnocline pycnocline gradients set to zero 1860 ELSE ! Stable conditions 1861 ! If pycnocline profile only defined when depth steady of increasing. 1862 IF ( zdhdt(ji,jj) > 0.0_wp ) THEN ! Depth increasing, or steady. 1863 IF ( zdb_bl(ji,jj) > 0.0_wp ) THEN 1864 IF ( zhol(ji,jj) >= 0.5_wp ) THEN ! Very stable - 'thick' pycnocline 1865 ztmp = 1.0_wp / MAX( zhbl(ji,jj), epsln ) 1866 zbgrad = zdb_bl(ji,jj) * ztmp 1867 DO jk = 2, ibld(ji,jj) 1868 znd = gdepw(ji,jj,jk,Kmm) * ztmp 1869 pdbdz(ji,jj,jk) = zbgrad * EXP( -15.0_wp * ( znd - 0.9_wp )**2 ) 1870 END DO 1871 ELSE ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form. 1872 ztmp = 1.0_wp / MAX( zdh(ji,jj), epsln ) 1873 zbgrad = zdb_bl(ji,jj) * ztmp 1874 DO jk = 2, ibld(ji,jj) 1875 znd = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) * ztmp 1876 pdbdz(ji,jj,jk) = zbgrad * EXP( -1.75_wp * ( znd + 0.75_wp )**2 ) 1877 END DO 1878 ENDIF ! IF (zhol >=0.5) 1879 ENDIF ! IF (zdb_bl> 0.) 1880 ENDIF ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero 1881 ENDIF ! IF (lconv) 1882 ENDIF ! IF ( ibld(ji,jj) < mbkt(ji,jj) ) 1867 1883 END_2D 1868 IF( ln_timing ) CALL timing_stop('zdf_osm_pshp') 1869 END SUBROUTINE zdf_osm_pycnocline_shear_profiles 1884 ! 1885 IF( ln_timing ) CALL timing_stop('zdf_osm_pscp') 1886 ! 1887 END SUBROUTINE zdf_osm_pycnocline_buoyancy_profiles 1870 1888 1871 1889 SUBROUTINE zdf_osm_calculate_dhdt( zdhdt, zddhdt ) … … 2609 2627 END IF 2610 2628 2629 ! Flags associated with diagnostic output 2630 IF ( ln_dia_osm .AND. ( iom_use("zdudz_pyc") .OR. iom_use("zdvdz_pyc") ) ) ln_dia_pyc_shr = .TRUE. 2631 IF ( ln_dia_osm .AND. ( iom_use("zdtdz_pyc") .OR. iom_use("zdsdz_pyc") .OR. iom_use("zdbdz_pyc" ) ) ) ln_dia_pyc_scl = .TRUE. 2632 2611 2633 ! ! allocate zdfosm arrays 2612 2634 IF( zdf_osm_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_osm_init : unable to allocate arrays' )
Note: See TracChangeset
for help on using the changeset viewer.