Changeset 12912 for NEMO/branches/UKMO/NEMO_4.0.2_FKOSM/src
- Timestamp:
- 2020-05-12T14:51:37+02:00 (4 years ago)
- Location:
- NEMO/branches/UKMO/NEMO_4.0.2_FKOSM/src/OCE
- Files:
-
- 7 edited
-
DIA/diawri.F90 (modified) (2 diffs)
-
SBC/sbcblk.F90 (modified) (11 diffs)
-
TRA/tramle.F90 (modified) (2 diffs)
-
TRD/trd_oce.F90 (modified) (2 diffs)
-
TRD/trdtra.F90 (modified) (1 diff)
-
ZDF/zdfosm.F90 (modified) (58 diffs)
-
ZDF/zdfphy.F90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/UKMO/NEMO_4.0.2_FKOSM/src/OCE/DIA/diawri.F90
r12658 r12912 43 43 USE zdfdrg ! ocean vertical physics: top/bottom friction 44 44 USE zdfmxl ! mixed layer 45 USE zdfosm ! mixed layer 45 46 ! 46 47 USE lbclnk ! ocean lateral boundary conditions (or mpp link) … … 923 924 CALL iom_rstput( 0, 0, inum, 'sdvecrtz', wsd ) ! now StokesDrift k-velocity 924 925 ENDIF 926 927 IF( ln_zdfosm ) THEN 928 CALL iom_rstput( 0, 0, inum, 'hbl', hbl*tmask(:,:,1) ) ! now boundary-layer depth 929 CALL iom_rstput( 0, 0, inum, 'hml', hml*tmask(:,:,1) ) ! now mixed-layer depth 930 CALL iom_rstput( 0, 0, inum, 'avt_k', avt_k*wmask ) ! w-level diffusion 931 CALL iom_rstput( 0, 0, inum, 'avm_k', avm_k*wmask ) ! now w-level viscosity 932 CALL iom_rstput( 0, 0, inum, 'ghamt', ghamt*wmask ) ! non-local t forcing 933 CALL iom_rstput( 0, 0, inum, 'ghams', ghams*wmask ) ! non-local s forcing 934 CALL iom_rstput( 0, 0, inum, 'ghamu', ghamu*wmask ) ! non-local u forcing 935 CALL iom_rstput( 0, 0, inum, 'ghamv', ghamu*wmask ) ! non-local v forcing 936 IF( ln_osm_mle ) THEN 937 CALL iom_rstput( 0, 0, inum, 'hmle', hmle*tmask(:,:,1) ) ! now transition-layer depth 938 END IF 939 ENDIF 925 940 926 941 #if defined key_si3 -
NEMO/branches/UKMO/NEMO_4.0.2_FKOSM/src/OCE/SBC/sbcblk.F90
r12658 r12912 124 124 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: cdn_oce, chn_oce, cen_oce ! needed by Lupkes 2015 bulk scheme 125 125 126 LOGICAL :: ln_humi_dpt = .FALSE. ! calculate specific hunidity from dewpoint 127 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: qair ! specific humidity of air at input height 128 126 129 INTEGER :: nblk ! choice of the bulk algorithm 127 130 ! ! associated indices: … … 145 148 !!------------------------------------------------------------------- 146 149 ALLOCATE( Cd_atm (jpi,jpj), Ch_atm (jpi,jpj), Ce_atm (jpi,jpj), t_zu(jpi,jpj), q_zu(jpi,jpj), & 147 & cdn_oce(jpi,jpj), chn_oce(jpi,jpj), cen_oce(jpi,jpj), STAT=sbc_blk_alloc )150 & cdn_oce(jpi,jpj), chn_oce(jpi,jpj), cen_oce(jpi,jpj), qair(jpi,jpj), STAT=sbc_blk_alloc ) 148 151 ! 149 152 CALL mpp_sum ( 'sbcblk', sbc_blk_alloc ) … … 171 174 NAMELIST/namsbc_blk/ sn_wndi, sn_wndj, sn_humi, sn_qsr, sn_qlw , & ! input fields 172 175 & sn_tair, sn_prec, sn_snow, sn_slp, sn_tdif, & 173 & ln_NCAR, ln_COARE_3p0, ln_COARE_3p5, ln_ECMWF, 176 & ln_NCAR, ln_COARE_3p0, ln_COARE_3p5, ln_ECMWF, ln_humi_dpt,& ! bulk algorithm 174 177 & cn_dir , ln_taudif, rn_zqt, rn_zu, & 175 178 & rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15 … … 323 326 ! 324 327 ! ! compute the surface ocean fluxes using bulk formulea 328 ! ..... if dewpoint supplied instead of specific humidaity, calculate specific humidity 329 IF(ln_humi_dpt) THEN 330 qair(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 331 ELSE 332 qair(:,:) = sf(jp_humi)%fnow(:,:,1) 333 END IF 334 325 335 IF( MOD( kt - 1, nn_fsbc ) == 0 ) CALL blk_oce( kt, sf, sst_m, ssu_m, ssv_m ) 326 336 … … 332 342 ENDIF 333 343 tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1) 334 qatm_ice(:,:) = sf(jp_humi)%fnow(:,:,1)344 qatm_ice(:,:) = qair(:,:) 335 345 tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac 336 346 sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac … … 434 444 !! (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2 435 445 !! (since reanalysis products provide T at z, not theta !) 436 ztpot = sf(jp_tair)%fnow(:,:,1) + gamma_moist( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1) ) * rn_zqt446 ztpot(:,:) = sf(jp_tair)%fnow(:,:,1) + gamma_moist( sf(jp_tair)%fnow(:,:,1), qair(:,:) ) * rn_zqt 437 447 438 448 SELECT CASE( nblk ) !== transfer coefficients ==! Cd, Ch, Ce at T-point 439 449 ! 440 CASE( np_NCAR ) ; CALL turb_ncar ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! NCAR-COREv2450 CASE( np_NCAR ) ; CALL turb_ncar ( rn_zqt, rn_zu, zst, ztpot, zsq, qair, wndm, & ! NCAR-COREv2 441 451 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 442 CASE( np_COARE_3p0 ) ; CALL turb_coare ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! COARE v3.0452 CASE( np_COARE_3p0 ) ; CALL turb_coare ( rn_zqt, rn_zu, zst, ztpot, zsq, qair, wndm, & ! COARE v3.0 443 453 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 444 CASE( np_COARE_3p5 ) ; CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! COARE v3.5454 CASE( np_COARE_3p5 ) ; CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, qair, wndm, & ! COARE v3.5 445 455 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 446 CASE( np_ECMWF ) ; CALL turb_ecmwf ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! ECMWF456 CASE( np_ECMWF ) ; CALL turb_ecmwf ( rn_zqt, rn_zu, zst, ztpot, zsq, qair, wndm, & ! ECMWF 447 457 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 448 458 CASE DEFAULT … … 454 464 zrhoa(:,:) = rho_air( t_zu(:,:) , q_zu(:,:) , sf(jp_slp)%fnow(:,:,1) ) 455 465 ELSE ! At zt: 456 zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )466 zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), qair(:,:), sf(jp_slp)%fnow(:,:,1) ) 457 467 END IF 458 468 … … 495 505 IF( ABS( rn_zu - rn_zqt) < 0.01_wp ) THEN 496 506 !! q_air and t_air are given at 10m (wind reference height) 497 zevap(:,:) = rn_efac*MAX( 0._wp, zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - sf(jp_humi)%fnow(:,:,1)) ) ! Evaporation, using bulk wind speed498 zqsb (:,:) = cp_air( sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - ztpot(:,:) ) ! Sensible Heat, using bulk wind speed507 zevap(:,:) = rn_efac*MAX( 0._wp, zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - qair(:,:)) ) ! Evaporation, using bulk wind speed 508 zqsb (:,:) = cp_air(qair(:,:))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - ztpot(:,:) ) ! Sensible Heat, using bulk wind speed 499 509 ELSE 500 510 !! q_air and t_air are not given at 10m (wind reference height) 501 511 ! Values of temp. and hum. adjusted to height of wind during bulk algorithm iteration must be used!!! 502 512 zevap(:,:) = rn_efac*MAX( 0._wp, zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - q_zu(:,:) ) ) ! Evaporation, using bulk wind speed 503 zqsb (:,:) = cp_air( sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - t_zu(:,:) ) ! Sensible Heat, using bulk wind speed513 zqsb (:,:) = cp_air(qair(:,:))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - t_zu(:,:) ) ! Sensible Heat, using bulk wind speed 504 514 ENDIF 505 515 … … 742 752 ! local scalars ( place there for vector optimisation purposes) 743 753 ! Computing density of air! Way denser that 1.2 over sea-ice !!! 744 zrhoa (:,:) = rho_air(sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1))754 zrhoa (:,:) = rho_air(sf(jp_tair)%fnow(:,:,1), qair(:,:), sf(jp_slp)%fnow(:,:,1)) 745 755 746 756 !!gm brutal.... … … 807 817 zcoef_dqla = -Ls * 11637800. * (-5897.8) 808 818 ! 809 zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )819 zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), qair(:,:), sf(jp_slp)%fnow(:,:,1) ) 810 820 ! 811 821 zztmp = 1. / ( 1. - albo ) … … 838 848 ! Latent Heat 839 849 qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * Ls * Ch_atm(ji,jj) * wndm_ice(ji,jj) * & 840 & ( 11637800. * EXP( -5897.8 * z1_st(ji,jj,jl) ) / zrhoa(ji,jj) - sf(jp_humi)%fnow(ji,jj,1) ) )850 & ( 11637800. * EXP( -5897.8 * z1_st(ji,jj,jl) ) / zrhoa(ji,jj) - qair(ji,jj) ) ) 841 851 ! Latent heat sensitivity for ice (Dqla/Dt) 842 852 IF( qla_ice(ji,jj,jl) > 0._wp ) THEN -
NEMO/branches/UKMO/NEMO_4.0.2_FKOSM/src/OCE/TRA/tramle.F90
r12658 r12912 21 21 USE lbclnk ! lateral boundary condition / mpp link 22 22 23 ! where OSMOSIS_OBL is used with integrated FK 24 USE zdf_oce, ONLY : ln_zdfosm 25 USE zdfosm, ONLY : ln_osm_mle, hmle, dbdx_mle, dbdy_mle, mld_prof 26 23 27 IMPLICIT NONE 24 28 PRIVATE … … 56 60 CONTAINS 57 61 58 SUBROUTINE tra_mle_trp( kt, kit000, pu, pv, pw, cdtype ) 59 !!---------------------------------------------------------------------- 60 !! *** ROUTINE tra_mle_trp *** 61 !! 62 !! ** Purpose : Add to the transport the Mixed Layer Eddy induced transport 63 !! 64 !! ** Method : The 3 components of the Mixed Layer Eddy (MLE) induced 65 !! transport are computed as follows : 66 !! zu_mle = dk[ zpsi_uw ] 67 !! zv_mle = dk[ zpsi_vw ] 68 !! zw_mle = - di[ zpsi_uw ] - dj[ zpsi_vw ] 69 !! where zpsi is the MLE streamfunction at uw and vw points (see the doc) 70 !! and added to the input velocity : 71 !! p.n = p.n + z._mle 72 !! 73 !! ** Action : - (pun,pvn,pwn) increased by the mle transport 74 !! CAUTION, the transport is not updated at the last line/raw 75 !! this may be a problem for some advection schemes 76 !! 77 !! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008 78 !! Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008 79 !!---------------------------------------------------------------------- 80 INTEGER , INTENT(in ) :: kt ! ocean time-step index 81 INTEGER , INTENT(in ) :: kit000 ! first time step index 82 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 83 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu ! in : 3 ocean transport components 84 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pv ! out: same 3 transport components 85 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pw ! increased by the MLE induced transport 86 ! 87 INTEGER :: ji, jj, jk ! dummy loop indices 88 INTEGER :: ii, ij, ik, ikmax ! local integers 89 REAL(wp) :: zcuw, zmuw, zc ! local scalar 90 REAL(wp) :: zcvw, zmvw ! - - 91 INTEGER , DIMENSION(jpi,jpj) :: inml_mle 92 REAL(wp), DIMENSION(jpi,jpj) :: zpsim_u, zpsim_v, zmld, zbm, zhu, zhv, zn2, zLf_NH, zLf_MH 93 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpsi_uw, zpsi_vw 94 !!---------------------------------------------------------------------- 95 ! 96 ! !== MLD used for MLE ==! 97 ! ! compute from the 10m density to deal with the diurnal cycle 98 inml_mle(:,:) = mbkt(:,:) + 1 ! init. to number of ocean w-level (T-level + 1) 99 IF ( nla10 > 0 ) THEN ! avoid case where first level is thicker than 10m 100 DO jk = jpkm1, nlb10, -1 ! from the bottom to nlb10 (10m) 101 DO jj = 1, jpj 102 DO ji = 1, jpi ! index of the w-level at the ML based 103 IF( rhop(ji,jj,jk) > rhop(ji,jj,nla10) + rn_rho_c_mle ) inml_mle(ji,jj) = jk ! Mixed layer 104 END DO 105 END DO 106 END DO 107 ENDIF 108 ikmax = MIN( MAXVAL( inml_mle(:,:) ), jpkm1 ) ! max level of the computation 109 ! 110 ! 111 zmld(:,:) = 0._wp !== Horizontal shape of the MLE ==! 112 zbm (:,:) = 0._wp 113 zn2 (:,:) = 0._wp 114 DO jk = 1, ikmax ! MLD and mean buoyancy and N2 over the mixed layer 115 DO jj = 1, jpj 116 DO ji = 1, jpi 117 zc = e3t_n(ji,jj,jk) * REAL( MIN( MAX( 0, inml_mle(ji,jj)-jk ) , 1 ) ) ! zc being 0 outside the ML t-points 118 zmld(ji,jj) = zmld(ji,jj) + zc 119 zbm (ji,jj) = zbm (ji,jj) + zc * (rau0 - rhop(ji,jj,jk) ) * r1_rau0 120 zn2 (ji,jj) = zn2 (ji,jj) + zc * (rn2(ji,jj,jk)+rn2(ji,jj,jk+1))*0.5_wp 121 END DO 122 END DO 123 END DO 124 125 SELECT CASE( nn_mld_uv ) ! MLD at u- & v-pts 126 CASE ( 0 ) != min of the 2 neighbour MLDs 127 DO jj = 1, jpjm1 128 DO ji = 1, fs_jpim1 ! vector opt. 129 zhu(ji,jj) = MIN( zmld(ji+1,jj), zmld(ji,jj) ) 130 zhv(ji,jj) = MIN( zmld(ji,jj+1), zmld(ji,jj) ) 131 END DO 132 END DO 133 CASE ( 1 ) != average of the 2 neighbour MLDs 134 DO jj = 1, jpjm1 135 DO ji = 1, fs_jpim1 ! vector opt. 136 zhu(ji,jj) = ( zmld(ji+1,jj) + zmld(ji,jj) ) * 0.5_wp 137 zhv(ji,jj) = ( zmld(ji,jj+1) + zmld(ji,jj) ) * 0.5_wp 138 END DO 139 END DO 140 CASE ( 2 ) != max of the 2 neighbour MLDs 141 DO jj = 1, jpjm1 142 DO ji = 1, fs_jpim1 ! vector opt. 143 zhu(ji,jj) = MAX( zmld(ji+1,jj), zmld(ji,jj) ) 144 zhv(ji,jj) = MAX( zmld(ji,jj+1), zmld(ji,jj) ) 145 END DO 146 END DO 147 END SELECT 148 ! ! convert density into buoyancy 149 zbm(:,:) = + grav * zbm(:,:) / MAX( e3t_n(:,:,1), zmld(:,:) ) 150 ! 151 ! 152 ! !== Magnitude of the MLE stream function ==! 153 ! 154 ! di[bm] Ds 155 ! Psi = Ce H^2 ---------------- e2u mu(z) where fu Lf = MAX( fu*rn_fl , (Db H)^1/2 ) 156 ! e1u Lf fu and the e2u for the "transport" 157 ! (not *e3u as divided by e3u at the end) 158 ! 159 IF( nn_mle == 0 ) THEN ! Fox-Kemper et al. 2010 formulation 160 DO jj = 1, jpjm1 161 DO ji = 1, fs_jpim1 ! vector opt. 162 zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj) * e2_e1u(ji,jj) & 163 & * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) ) & 164 & / ( MAX( rn_lf * rfu(ji,jj) , SQRT( rb_c * zhu(ji,jj) ) ) ) 165 ! 166 zpsim_v(ji,jj) = rn_ce * zhv(ji,jj) * zhv(ji,jj) * e1_e2v(ji,jj) & 167 & * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) ) & 168 & / ( MAX( rn_lf * rfv(ji,jj) , SQRT( rb_c * zhv(ji,jj) ) ) ) 169 END DO 170 END DO 171 ! 172 ELSEIF( nn_mle == 1 ) THEN ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat) 173 DO jj = 1, jpjm1 174 DO ji = 1, fs_jpim1 ! vector opt. 175 zpsim_u(ji,jj) = rc_f * zhu(ji,jj) * zhu(ji,jj) * e2_e1u(ji,jj) & 176 & * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) ) 177 ! 178 zpsim_v(ji,jj) = rc_f * zhv(ji,jj) * zhv(ji,jj) * e1_e2v(ji,jj) & 179 & * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) ) 180 END DO 181 END DO 182 ENDIF 183 ! 184 IF( nn_conv == 1 ) THEN ! No MLE in case of convection 185 DO jj = 1, jpjm1 186 DO ji = 1, fs_jpim1 ! vector opt. 187 IF( MIN( zn2(ji,jj) , zn2(ji+1,jj) ) < 0._wp ) zpsim_u(ji,jj) = 0._wp 188 IF( MIN( zn2(ji,jj) , zn2(ji,jj+1) ) < 0._wp ) zpsim_v(ji,jj) = 0._wp 189 END DO 190 END DO 191 ENDIF 192 ! 193 ! !== structure function value at uw- and vw-points ==! 194 DO jj = 1, jpjm1 195 DO ji = 1, fs_jpim1 ! vector opt. 196 zhu(ji,jj) = 1._wp / zhu(ji,jj) ! hu --> 1/hu 197 zhv(ji,jj) = 1._wp / zhv(ji,jj) 198 END DO 199 END DO 200 ! 201 zpsi_uw(:,:,:) = 0._wp 202 zpsi_vw(:,:,:) = 0._wp 203 ! 204 DO jk = 2, ikmax ! start from 2 : surface value = 0 205 DO jj = 1, jpjm1 206 DO ji = 1, fs_jpim1 ! vector opt. 207 zcuw = 1._wp - ( gdepw_n(ji+1,jj,jk) + gdepw_n(ji,jj,jk) ) * zhu(ji,jj) 208 zcvw = 1._wp - ( gdepw_n(ji,jj+1,jk) + gdepw_n(ji,jj,jk) ) * zhv(ji,jj) 209 zcuw = zcuw * zcuw 210 zcvw = zcvw * zcvw 211 zmuw = MAX( 0._wp , ( 1._wp - zcuw ) * ( 1._wp + r5_21 * zcuw ) ) 212 zmvw = MAX( 0._wp , ( 1._wp - zcvw ) * ( 1._wp + r5_21 * zcvw ) ) 213 ! 214 zpsi_uw(ji,jj,jk) = zpsim_u(ji,jj) * zmuw * umask(ji,jj,jk) 215 zpsi_vw(ji,jj,jk) = zpsim_v(ji,jj) * zmvw * vmask(ji,jj,jk) 216 END DO 217 END DO 218 END DO 219 ! 220 ! !== transport increased by the MLE induced transport ==! 221 DO jk = 1, ikmax 222 DO jj = 1, jpjm1 ! CAUTION pu,pv must be defined at row/column i=1 / j=1 223 DO ji = 1, fs_jpim1 ! vector opt. 224 pu(ji,jj,jk) = pu(ji,jj,jk) + ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji,jj,jk+1) ) 225 pv(ji,jj,jk) = pv(ji,jj,jk) + ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj,jk+1) ) 226 END DO 227 END DO 228 DO jj = 2, jpjm1 229 DO ji = fs_2, fs_jpim1 ! vector opt. 230 pw(ji,jj,jk) = pw(ji,jj,jk) - ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji-1,jj,jk) & 62 SUBROUTINE tra_mle_trp( kt, kit000, pu, pv, pw, cdtype ) 63 !!---------------------------------------------------------------------- 64 !! *** ROUTINE tra_mle_trp *** 65 !! 66 !! ** Purpose : Add to the transport the Mixed Layer Eddy induced transport 67 !! 68 !! ** Method : The 3 components of the Mixed Layer Eddy (MLE) induced 69 !! transport are computed as follows : 70 !! zu_mle = dk[ zpsi_uw ] 71 !! zv_mle = dk[ zpsi_vw ] 72 !! zw_mle = - di[ zpsi_uw ] - dj[ zpsi_vw ] 73 !! where zpsi is the MLE streamfunction at uw and vw points (see the doc) 74 !! and added to the input velocity : 75 !! p.n = p.n + z._mle 76 !! 77 !! ** Action : - (pun,pvn,pwn) increased by the mle transport 78 !! CAUTION, the transport is not updated at the last line/raw 79 !! this may be a problem for some advection schemes 80 !! 81 !! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008 82 !! Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008 83 !!---------------------------------------------------------------------- 84 INTEGER , INTENT(in ) :: kt ! ocean time-step index 85 INTEGER , INTENT(in ) :: kit000 ! first time step index 86 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 87 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu ! in : 3 ocean transport components 88 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pv ! out: same 3 transport components 89 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pw ! increased by the MLE induced transport 90 ! 91 INTEGER :: ji, jj, jk ! dummy loop indices 92 INTEGER :: ii, ij, ik, ikmax ! local integers 93 REAL(wp) :: zcuw, zmuw, zc ! local scalar 94 REAL(wp) :: zcvw, zmvw ! - - 95 INTEGER , DIMENSION(jpi,jpj) :: inml_mle 96 REAL(wp), DIMENSION(jpi,jpj) :: zpsim_u, zpsim_v, zmld, zbm, zhu, zhv, zn2, zLf_NH, zLf_MH 97 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpsi_uw, zpsi_vw 98 !!---------------------------------------------------------------------- 99 ! 100 ! 101 IF(ln_osm_mle.and.ln_zdfosm) THEN 102 ikmax = MIN( MAXVAL( mld_prof(:,:) ), jpkm1 ) ! max level of the computation 103 ! 104 ! 105 SELECT CASE( nn_mld_uv ) ! MLD at u- & v-pts 106 CASE ( 0 ) != min of the 2 neighbour MLDs 107 DO jj = 1, jpjm1 108 DO ji = 1, fs_jpim1 ! vector opt. 109 zhu(ji,jj) = MIN( hmle(ji+1,jj), hmle(ji,jj) ) 110 zhv(ji,jj) = MIN( hmle(ji,jj+1), hmle(ji,jj) ) 111 END DO 112 END DO 113 CASE ( 1 ) != average of the 2 neighbour MLDs 114 DO jj = 1, jpjm1 115 DO ji = 1, fs_jpim1 ! vector opt. 116 zhu(ji,jj) = MAX( hmle(ji+1,jj), hmle(ji,jj) ) 117 zhv(ji,jj) = MAX( hmle(ji,jj+1), hmle(ji,jj) ) 118 END DO 119 END DO 120 CASE ( 2 ) != max of the 2 neighbour MLDs 121 DO jj = 1, jpjm1 122 DO ji = 1, fs_jpim1 ! vector opt. 123 zhu(ji,jj) = MAX( hmle(ji+1,jj), hmle(ji,jj) ) 124 zhv(ji,jj) = MAX( hmle(ji,jj+1), hmle(ji,jj) ) 125 END DO 126 END DO 127 END SELECT 128 IF( nn_mle == 0 ) THEN ! Fox-Kemper et al. 2010 formulation 129 DO jj = 1, jpjm1 130 DO ji = 1, fs_jpim1 ! vector opt. 131 zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj) * e2u(ji,jj) & 132 & * dbdx_mle(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) & 133 & / ( MAX( rn_lf * rfu(ji,jj) , SQRT( rb_c * zhu(ji,jj) ) ) ) 134 ! 135 zpsim_v(ji,jj) = rn_ce * zhv(ji,jj) * zhv(ji,jj) * e1v(ji,jj) & 136 & * dbdy_mle(ji,jj) * MIN( 111.e3_wp , e2v(ji,jj) ) & 137 & / ( MAX( rn_lf * rfv(ji,jj) , SQRT( rb_c * zhv(ji,jj) ) ) ) 138 END DO 139 END DO 140 ! 141 ELSEIF( nn_mle == 1 ) THEN ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat) 142 DO jj = 1, jpjm1 143 DO ji = 1, fs_jpim1 ! vector opt. 144 zpsim_u(ji,jj) = rc_f * zhu(ji,jj) * zhu(ji,jj) * e2u(ji,jj) & 145 & * dbdx_mle(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) 146 ! 147 zpsim_v(ji,jj) = rc_f * zhv(ji,jj) * zhv(ji,jj) * e1v(ji,jj) & 148 & * dbdy_mle(ji,jj) * MIN( 111.e3_wp , e2v(ji,jj) ) 149 END DO 150 END DO 151 ENDIF 152 153 ELSE !do not use osn_mle 154 ! !== MLD used for MLE ==! 155 ! ! compute from the 10m density to deal with the diurnal cycle 156 inml_mle(:,:) = mbkt(:,:) + 1 ! init. to number of ocean w-level (T-level + 1) 157 IF ( nla10 > 0 ) THEN ! avoid case where first level is thicker than 10m 158 DO jk = jpkm1, nlb10, -1 ! from the bottom to nlb10 (10m) 159 DO jj = 1, jpj 160 DO ji = 1, jpi ! index of the w-level at the ML based 161 IF( rhop(ji,jj,jk) > rhop(ji,jj,nla10) + rn_rho_c_mle ) inml_mle(ji,jj) = jk ! Mixed layer 162 END DO 163 END DO 164 END DO 165 ENDIF 166 ikmax = MIN( MAXVAL( inml_mle(:,:) ), jpkm1 ) ! max level of the computation 167 168 ! 169 ! 170 zmld(:,:) = 0._wp !== Horizontal shape of the MLE ==! 171 zbm (:,:) = 0._wp 172 zn2 (:,:) = 0._wp 173 DO jk = 1, ikmax ! MLD and mean buoyancy and N2 over the mixed layer 174 DO jj = 1, jpj 175 DO ji = 1, jpi 176 zc = e3t_n(ji,jj,jk) * REAL( MIN( MAX( 0, inml_mle(ji,jj)-jk ) , 1 ) ) ! zc being 0 outside the ML t-points 177 zmld(ji,jj) = zmld(ji,jj) + zc 178 zbm (ji,jj) = zbm (ji,jj) + zc * (rau0 - rhop(ji,jj,jk) ) * r1_rau0 179 zn2 (ji,jj) = zn2 (ji,jj) + zc * (rn2(ji,jj,jk)+rn2(ji,jj,jk+1))*0.5_wp 180 END DO 181 END DO 182 END DO 183 184 SELECT CASE( nn_mld_uv ) ! MLD at u- & v-pts 185 CASE ( 0 ) != min of the 2 neighbour MLDs 186 DO jj = 1, jpjm1 187 DO ji = 1, fs_jpim1 ! vector opt. 188 zhu(ji,jj) = MIN( zmld(ji+1,jj), zmld(ji,jj) ) 189 zhv(ji,jj) = MIN( zmld(ji,jj+1), zmld(ji,jj) ) 190 END DO 191 END DO 192 CASE ( 1 ) != average of the 2 neighbour MLDs 193 DO jj = 1, jpjm1 194 DO ji = 1, fs_jpim1 ! vector opt. 195 zhu(ji,jj) = ( zmld(ji+1,jj) + zmld(ji,jj) ) * 0.5_wp 196 zhv(ji,jj) = ( zmld(ji,jj+1) + zmld(ji,jj) ) * 0.5_wp 197 END DO 198 END DO 199 CASE ( 2 ) != max of the 2 neighbour MLDs 200 DO jj = 1, jpjm1 201 DO ji = 1, fs_jpim1 ! vector opt. 202 zhu(ji,jj) = MAX( zmld(ji+1,jj), zmld(ji,jj) ) 203 zhv(ji,jj) = MAX( zmld(ji,jj+1), zmld(ji,jj) ) 204 END DO 205 END DO 206 END SELECT 207 ! ! convert density into buoyancy 208 zbm(:,:) = + grav * zbm(:,:) / MAX( e3t_n(:,:,1), zmld(:,:) ) 209 ! 210 ! 211 ! !== Magnitude of the MLE stream function ==! 212 ! 213 ! di[bm] Ds 214 ! Psi = Ce H^2 ---------------- e2u mu(z) where fu Lf = MAX( fu*rn_fl , (Db H)^1/2 ) 215 ! e1u Lf fu and the e2u for the "transport" 216 ! (not *e3u as divided by e3u at the end) 217 ! 218 IF( nn_mle == 0 ) THEN ! Fox-Kemper et al. 2010 formulation 219 DO jj = 1, jpjm1 220 DO ji = 1, fs_jpim1 ! vector opt. 221 zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj) * e2_e1u(ji,jj) & 222 & * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) ) & 223 & / ( MAX( rn_lf * rfu(ji,jj) , SQRT( rb_c * zhu(ji,jj) ) ) ) 224 ! 225 zpsim_v(ji,jj) = rn_ce * zhv(ji,jj) * zhv(ji,jj) * e1_e2v(ji,jj) & 226 & * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) ) & 227 & / ( MAX( rn_lf * rfv(ji,jj) , SQRT( rb_c * zhv(ji,jj) ) ) ) 228 END DO 229 END DO 230 ! 231 ELSEIF( nn_mle == 1 ) THEN ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat) 232 DO jj = 1, jpjm1 233 DO ji = 1, fs_jpim1 ! vector opt. 234 zpsim_u(ji,jj) = rc_f * zhu(ji,jj) * zhu(ji,jj) * e2_e1u(ji,jj) & 235 & * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) ) 236 ! 237 zpsim_v(ji,jj) = rc_f * zhv(ji,jj) * zhv(ji,jj) * e1_e2v(ji,jj) & 238 & * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) ) 239 END DO 240 END DO 241 ENDIF 242 ! 243 IF( nn_conv == 1 ) THEN ! No MLE in case of convection 244 DO jj = 1, jpjm1 245 DO ji = 1, fs_jpim1 ! vector opt. 246 IF( MIN( zn2(ji,jj) , zn2(ji+1,jj) ) < 0._wp ) zpsim_u(ji,jj) = 0._wp 247 IF( MIN( zn2(ji,jj) , zn2(ji,jj+1) ) < 0._wp ) zpsim_v(ji,jj) = 0._wp 248 END DO 249 END DO 250 ENDIF 251 ! 252 ENDIF ! end of lm_osm_mle loop 253 ! !== structure function value at uw- and vw-points ==! 254 DO jj = 1, jpjm1 255 DO ji = 1, fs_jpim1 ! vector opt. 256 zhu(ji,jj) = 1._wp / MAX(zhu(ji,jj), rsmall) ! hu --> 1/hu 257 zhv(ji,jj) = 1._wp / MAX(zhv(ji,jj), rsmall) 258 END DO 259 END DO 260 ! 261 zpsi_uw(:,:,:) = 0._wp 262 zpsi_vw(:,:,:) = 0._wp 263 ! 264 DO jk = 2, ikmax ! start from 2 : surface value = 0 265 DO jj = 1, jpjm1 266 DO ji = 1, fs_jpim1 ! vector opt. 267 zcuw = 1._wp - ( gdepw_n(ji+1,jj,jk) + gdepw_n(ji,jj,jk) ) * zhu(ji,jj) 268 zcvw = 1._wp - ( gdepw_n(ji,jj+1,jk) + gdepw_n(ji,jj,jk) ) * zhv(ji,jj) 269 zcuw = zcuw * zcuw 270 zcvw = zcvw * zcvw 271 zmuw = MAX( 0._wp , ( 1._wp - zcuw ) * ( 1._wp + r5_21 * zcuw ) ) 272 zmvw = MAX( 0._wp , ( 1._wp - zcvw ) * ( 1._wp + r5_21 * zcvw ) ) 273 ! 274 zpsi_uw(ji,jj,jk) = zpsim_u(ji,jj) * zmuw * umask(ji,jj,jk) 275 zpsi_vw(ji,jj,jk) = zpsim_v(ji,jj) * zmvw * vmask(ji,jj,jk) 276 END DO 277 END DO 278 END DO 279 ! 280 ! !== transport increased by the MLE induced transport ==! 281 DO jk = 1, ikmax 282 DO jj = 1, jpjm1 ! CAUTION pu,pv must be defined at row/column i=1 / j=1 283 DO ji = 1, fs_jpim1 ! vector opt. 284 pu(ji,jj,jk) = pu(ji,jj,jk) + ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji,jj,jk+1) ) 285 pv(ji,jj,jk) = pv(ji,jj,jk) + ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj,jk+1) ) 286 END DO 287 END DO 288 DO jj = 2, jpjm1 289 DO ji = fs_2, fs_jpim1 ! vector opt. 290 pw(ji,jj,jk) = pw(ji,jj,jk) - ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji-1,jj,jk) & 231 291 & + zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj-1,jk) ) 232 END DO 233 END DO 234 END DO 235 236 IF( cdtype == 'TRA') THEN !== outputs ==! 237 ! 238 zLf_NH(:,:) = SQRT( rb_c * zmld(:,:) ) * r1_ft(:,:) ! Lf = N H / f 239 CALL iom_put( "Lf_NHpf" , zLf_NH ) ! Lf = N H / f 240 ! 241 ! divide by cross distance to give streamfunction with dimensions m^2/s 242 DO jk = 1, ikmax+1 243 zpsi_uw(:,:,jk) = zpsi_uw(:,:,jk) * r1_e2u(:,:) 244 zpsi_vw(:,:,jk) = zpsi_vw(:,:,jk) * r1_e1v(:,:) 245 END DO 246 CALL iom_put( "psiu_mle", zpsi_uw ) ! i-mle streamfunction 247 CALL iom_put( "psiv_mle", zpsi_vw ) ! j-mle streamfunction 248 ENDIF 249 ! 250 END SUBROUTINE tra_mle_trp 292 END DO 293 END DO 294 END DO 295 296 IF( cdtype == 'TRA') THEN !== outputs ==! 297 ! 298 IF (ln_osm_mle.and.ln_zdfosm) THEN 299 zLf_NH(:,:) = SQRT( rb_c * hmle(:,:) ) * r1_ft(:,:) ! Lf = N H / f 300 ELSE 301 zLf_NH(:,:) = SQRT( rb_c * zmld(:,:) ) * r1_ft(:,:) ! Lf = N H / f 302 END IF 303 CALL iom_put( "Lf_NHpf" , zLf_NH ) ! Lf = N H / f 304 ! 305 ! divide by cross distance to give streamfunction with dimensions m^2/s 306 DO jk = 1, ikmax+1 307 zpsi_uw(:,:,jk) = zpsi_uw(:,:,jk) * r1_e2u(:,:) 308 zpsi_vw(:,:,jk) = zpsi_vw(:,:,jk) * r1_e1v(:,:) 309 END DO 310 CALL iom_put( "psiu_mle", zpsi_uw ) ! i-mle streamfunction 311 CALL iom_put( "psiv_mle", zpsi_vw ) ! j-mle streamfunction 312 ENDIF 313 ! 314 END SUBROUTINE tra_mle_trp 251 315 252 316 -
NEMO/branches/UKMO/NEMO_4.0.2_FKOSM/src/OCE/TRD/trd_oce.F90
r12658 r12912 33 33 # endif 34 34 ! !!!* Active tracers trends indexes 35 INTEGER, PUBLIC, PARAMETER :: jptot_tra = 2 0!: Total trend nb: change it when adding/removing one indice below35 INTEGER, PUBLIC, PARAMETER :: jptot_tra = 21 !: Total trend nb: change it when adding/removing one indice below 36 36 ! =============== ! 37 37 INTEGER, PUBLIC, PARAMETER :: jptra_xad = 1 !: x- horizontal advection … … 46 46 INTEGER, PUBLIC, PARAMETER :: jptra_bbc = 10 !: Bottom Boundary Condition (geoth. heating) 47 47 INTEGER, PUBLIC, PARAMETER :: jptra_bbl = 11 !: Bottom Boundary Layer (diffusive and/or advective) 48 INTEGER, PUBLIC, PARAMETER :: jptra_osm = 21 !: Non-local terms from OSMOSIS OBL model 48 49 INTEGER, PUBLIC, PARAMETER :: jptra_npc = 12 !: non-penetrative convection treatment 49 50 INTEGER, PUBLIC, PARAMETER :: jptra_dmp = 13 !: internal restoring (damping) -
NEMO/branches/UKMO/NEMO_4.0.2_FKOSM/src/OCE/TRD/trdtra.F90
r12658 r12912 346 346 CASE( jptra_bbl ) ; CALL iom_put( "ttrd_bbl" , ptrdx ) ! bottom boundary layer 347 347 CALL iom_put( "strd_bbl" , ptrdy ) 348 CASE( jptra_osm ) ; CALL iom_put( "ttrd_osm" , ptrdx ) ! OSMOSIS non-local forcing 349 CALL iom_put( "strd_osm" , ptrdy ) 348 350 CASE( jptra_npc ) ; CALL iom_put( "ttrd_npc" , ptrdx ) ! static instability mixing 349 351 CALL iom_put( "strd_npc" , ptrdy ) -
NEMO/branches/UKMO/NEMO_4.0.2_FKOSM/src/OCE/ZDF/zdfosm.F90
r12658 r12912 25 25 !! (12) Replace zwstrl with zvstr in calculation of eddy viscosity. 26 26 !! 27/09/2017 (13) Calculate Stokes drift and Stokes penetration depth from wave information 27 !! (14) B ouyancy flux due to entrainment changed to include contribution from shear turbulence (for testing commented out).27 !! (14) Buoyancy flux due to entrainment changed to include contribution from shear turbulence. 28 28 !! 28/09/2017 (15) Calculation of Stokes drift moved into separate do-loops to allow for different options for the determining the Stokes drift to be added. 29 29 !! (16) Calculation of Stokes drift from windspeed for PM spectrum (for testing, commented out) 30 30 !! (17) Modification to Langmuir velocity scale to include effects due to the Stokes penetration depth (for testing, commented out) 31 !! ??/??/2018 (18) Revision to code structure, selected using key_osmldpth1. Inline code moved into subroutines. Changes to physics made, 32 !! (a) Pycnocline temperature and salinity profies changed for unstable layers 33 !! (b) The stable OSBL depth parametrization changed. 34 !! 16/05/2019 (19) Fox-Kemper parametrization of restratification through mixed layer eddies added to revised code. 35 !! 23/05/19 (20) Old code where key_osmldpth1` is *not* set removed, together with the key key_osmldpth1 31 36 !!---------------------------------------------------------------------- 32 37 … … 40 45 !! trc_osm : compute and add to the passive tracer trend the non-local flux (TBD) 41 46 !! dyn_osm : compute and add to u & v trensd the non-local flux 47 !! 48 !! Subroutines in revised code. 42 49 !!---------------------------------------------------------------------- 43 50 USE oce ! ocean dynamics and active tracers … … 69 76 PUBLIC tra_osm ! routine called by step.F90 70 77 PUBLIC trc_osm ! routine called by trcstp.F90 71 PUBLIC dyn_osm ! routine called by 'step.F90' 78 PUBLIC dyn_osm ! routine called by step.F90 79 80 PUBLIC ln_osm_mle ! logical needed by tra_mle_init in tramle.F90 72 81 73 82 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ghamu !: non-local u-momentum flux … … 77 86 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: etmean !: averaging operator for avt 78 87 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hbl !: boundary layer depth 79 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hbli !: intial boundary layer depth for stable blayer 88 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: dh ! depth of pycnocline 89 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hml ! ML depth 80 90 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: dstokes !: penetration depth of the Stokes drift. 91 92 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_ft ! inverse of the modified Coriolis parameter at t-pts 93 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hmle ! Depth of layer affexted by mixed layer eddies in Fox-Kemper parametrization 94 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: dbdx_mle ! zonal buoyancy gradient in ML 95 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: dbdy_mle ! meridional buoyancy gradient in ML 96 INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mld_prof ! level of base of MLE layer. 81 97 82 98 ! !!** Namelist namzdf_osm ** 83 99 LOGICAL :: ln_use_osm_la ! Use namelist rn_osm_la 100 101 LOGICAL :: ln_osm_mle !: flag to activate the Mixed Layer Eddy (MLE) parameterisation 102 84 103 REAL(wp) :: rn_osm_la ! Turbulent Langmuir number 85 104 REAL(wp) :: rn_osm_dstokes ! Depth scale of Stokes drift … … 96 115 REAL(wp) :: rn_difconv = 1._wp ! diffusivity when unstable below BL (m2/s) 97 116 117 ! OSMOSIS mixed layer eddy parametrization constants 118 INTEGER :: nn_osm_mle ! = 0/1 flag for horizontal average on avt 119 REAL(wp) :: rn_osm_mle_ce ! MLE coefficient 120 ! ! parameters used in nn_osm_mle = 0 case 121 REAL(wp) :: rn_osm_mle_lf ! typical scale of mixed layer front 122 REAL(wp) :: rn_osm_mle_time ! time scale for mixing momentum across the mixed layer 123 ! ! parameters used in nn_osm_mle = 1 case 124 REAL(wp) :: rn_osm_mle_lat ! reference latitude for a 5 km scale of ML front 125 REAL(wp) :: rn_osm_mle_rho_c ! Density criterion for definition of MLD used by FK 126 REAL(wp) :: r5_21 = 5.e0 / 21.e0 ! factor used in mle streamfunction computation 127 REAL(wp) :: rb_c ! ML buoyancy criteria = g rho_c /rau0 where rho_c is defined in zdfmld 128 REAL(wp) :: rc_f ! MLE coefficient (= rn_ce / (5 km * fo) ) in nn_osm_mle=1 case 129 REAL(wp) :: rn_osm_mle_thresh ! Threshold buoyancy for deepening of MLE layer below OSBL base. 130 REAL(wp) :: rn_osm_mle_tau ! Adjustment timescale for MLE. 131 132 98 133 ! !!! ** General constants ** 99 REAL(wp) :: epsln = 1.0e-20_wp ! a small positive number 134 REAL(wp) :: epsln = 1.0e-20_wp ! a small positive number to ensure no div by zero 135 REAL(wp) :: depth_tol = 1.0e-6_wp ! a small-ish positive number to give a hbl slightly shallower than gdepw 100 136 REAL(wp) :: pthird = 1._wp/3._wp ! 1/3 101 137 REAL(wp) :: p2third = 2._wp/3._wp ! 2/3 … … 105 141 !!---------------------------------------------------------------------- 106 142 !! NEMO/OCE 4.0 , NEMO Consortium (2018) 107 !! $Id $143 !! $Id: zdfosm.F90 12317 2020-01-14 12:40:47Z agn $ 108 144 !! Software governed by the CeCILL license (see ./LICENSE) 109 145 !!---------------------------------------------------------------------- … … 114 150 !! *** FUNCTION zdf_osm_alloc *** 115 151 !!---------------------------------------------------------------------- 116 ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk), ghams(jpi,jpj,jpk), & 117 & hbl(jpi,jpj) , hbli(jpi,jpj) , dstokes(jpi, jpj) , & 118 & etmean(jpi,jpj,jpk), STAT= zdf_osm_alloc ) 152 ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk),ghams(jpi,jpj,jpk), & 153 & hbl(jpi,jpj), dh(jpi,jpj), hml(jpi,jpj), dstokes(jpi, jpj), & 154 & etmean(jpi,jpj,jpk), STAT= zdf_osm_alloc ) 155 156 ALLOCATE( hmle(jpi,jpj), r1_ft(jpi,jpj), dbdx_mle(jpi,jpj), dbdy_mle(jpi,jpj), & 157 & mld_prof(jpi,jpj), STAT= zdf_osm_alloc ) 158 159 ! ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk),ghams(jpi,jpj,jpk), & ! would ths be better ? 160 ! & hbl(jpi,jpj), dh(jpi,jpj), hml(jpi,jpj), dstokes(jpi, jpj), & 161 ! & etmean(jpi,jpj,jpk), STAT= zdf_osm_alloc ) 162 ! IF( zdf_osm_alloc /= 0 ) CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm arrays') 163 ! IF ( ln_osm_mle ) THEN 164 ! Allocate( hmle(jpi,jpj), r1_ft(jpi,jpj), STAT= zdf_osm_alloc ) 165 ! IF( zdf_osm_alloc /= 0 ) CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm mle arrays') 166 ! ENDIF 167 119 168 IF( zdf_osm_alloc /= 0 ) CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm arrays') 120 169 CALL mpp_sum ( 'zdfosm', zdf_osm_alloc ) … … 161 210 !! 162 211 INTEGER :: ji, jj, jk ! dummy loop indices 212 213 INTEGER :: jl ! dummy loop indices 214 163 215 INTEGER :: ikbot, jkmax, jkm1, jkp2 ! 164 216 … … 166 218 REAL(wp) :: zbeta, zthermal ! 167 219 REAL(wp) :: zehat, zeta, zhrib, zsig, zscale, zwst, zws, zwm ! Velocity scales 168 REAL(wp) :: zwsun, zwmun, zcons, zconm, zwcons, zwconm ! 220 REAL(wp) :: zwsun, zwmun, zcons, zconm, zwcons, zwconm ! 221 169 222 REAL(wp) :: zsr, zbw, ze, zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zcomp , zrhd,zrhdr,zbvzed ! In situ density 170 223 INTEGER :: jm ! dummy loop indices … … 191 244 REAL(wp), DIMENSION(jpi,jpj) :: zwbav ! Buoyancy flux - bl average 192 245 REAL(wp), DIMENSION(jpi,jpj) :: zwb_ent ! Buoyancy entrainment flux 246 247 248 REAL(wp), DIMENSION(jpi,jpj) :: zwb_fk_b ! MLE buoyancy flux averaged over OSBL 249 REAL(wp), DIMENSION(jpi,jpj) :: zwb_fk ! max MLE buoyancy flux 250 REAL(wp), DIMENSION(jpi,jpj) :: zdiff_mle ! extra MLE vertical diff 251 REAL(wp), DIMENSION(jpi,jpj) :: zvel_mle ! velocity scale for dhdt with stable ML and FK 252 193 253 REAL(wp), DIMENSION(jpi,jpj) :: zustke ! Surface Stokes drift 194 254 REAL(wp), DIMENSION(jpi,jpj) :: zla ! Trubulent Langmuir number … … 196 256 REAL(wp), DIMENSION(jpi,jpj) :: zsin_wind ! Sin angle of surface stress 197 257 REAL(wp), DIMENSION(jpi,jpj) :: zhol ! Stability parameter for boundary layer 198 LOGICAL, DIMENSION( :,:), ALLOCATABLE :: lconv! unstable/stable bl258 LOGICAL, DIMENSION(jpi,jpj) :: lconv ! unstable/stable bl 199 259 200 260 ! mixed-layer variables … … 208 268 REAL(wp), DIMENSION(jpi,jpj) :: zhbl ! bl depth - grid 209 269 REAL(wp), DIMENSION(jpi,jpj) :: zhml ! ml depth - grid 270 271 REAL(wp), DIMENSION(jpi,jpj) :: zhmle ! MLE depth - grid 272 REAL(wp), DIMENSION(jpi,jpj) :: zmld ! ML depth on grid 273 210 274 REAL(wp), DIMENSION(jpi,jpj) :: zdh ! pycnocline depth - grid 211 275 REAL(wp), DIMENSION(jpi,jpj) :: zdhdt ! BL depth tendency 276 REAL(wp), DIMENSION(jpi,jpj) :: zdhdt_2 ! correction to dhdt due to internal structure. 277 REAL(wp), DIMENSION(jpi,jpj) :: zdtdz_ext,zdsdz_ext,zdbdz_ext ! external temperature/salinity and buoyancy gradients 278 279 REAL(wp), DIMENSION(jpi,jpj) :: zdtdx, zdtdy, zdsdx, zdsdy ! horizontal gradients for Fox-Kemper parametrization. 280 212 281 REAL(wp), DIMENSION(jpi,jpj) :: zt_bl,zs_bl,zu_bl,zv_bl,zrh_bl ! averages over the depth of the blayer 213 282 REAL(wp), DIMENSION(jpi,jpj) :: zt_ml,zs_ml,zu_ml,zv_ml,zrh_ml ! averages over the depth of the mixed layer … … 238 307 ! Temporary variables 239 308 INTEGER :: inhml 240 INTEGER :: i_lconv_alloc241 309 REAL(wp) :: znd,znd_d,zznd_ml,zznd_pyc,zznd_d ! temporary non-dimensional depths used in various routines 242 310 REAL(wp) :: ztemp, zari, zpert, zzdhdt, zdb ! temporary variables … … 248 316 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdiffut ! t-diffusivity 249 317 318 INTEGER :: ibld_ext=0 ! does not have to be zero for modified scheme 319 REAL(wp) :: zwb_min, zgamma_b_nd, zgamma_b, zdhoh, ztau, zddhdt 320 REAL(wp) :: zzeta_s = 0._wp 321 REAL(wp) :: zzeta_v = 0.46 322 REAL(wp) :: zabsstke 323 250 324 ! For debugging 251 325 INTEGER :: ikt 252 326 !!-------------------------------------------------------------------- 253 327 ! 254 ALLOCATE( lconv(jpi,jpj), STAT= i_lconv_alloc )255 IF( i_lconv_alloc /= 0 ) CALL ctl_warn('zdf_osm: failed to allocate lconv')256 257 328 ibld(:,:) = 0 ; imld(:,:) = 0 258 329 zrad0(:,:) = 0._wp ; zradh(:,:) = 0._wp ; zradav(:,:) = 0._wp ; zustar(:,:) = 0._wp … … 268 339 zt_bl(:,:) = 0._wp ; zs_bl(:,:) = 0._wp ; zu_bl(:,:) = 0._wp ; zv_bl(:,:) = 0._wp 269 340 zrh_bl(:,:) = 0._wp ; zt_ml(:,:) = 0._wp ; zs_ml(:,:) = 0._wp ; zu_ml(:,:) = 0._wp 341 270 342 zv_ml(:,:) = 0._wp ; zrh_ml(:,:) = 0._wp ; zdt_bl(:,:) = 0._wp ; zds_bl(:,:) = 0._wp 271 343 zdu_bl(:,:) = 0._wp ; zdv_bl(:,:) = 0._wp ; zdrh_bl(:,:) = 0._wp ; zdb_bl(:,:) = 0._wp … … 277 349 zdudz_pyc(:,:,:) = 0._wp ; zdvdz_pyc(:,:,:) = 0._wp 278 350 ! 351 zdtdz_ext(:,:) = 0._wp ; zdsdz_ext(:,:) = 0._wp ; zdbdz_ext(:,:) = 0._wp 352 353 IF ( ln_osm_mle ) THEN ! only initialise arrays if needed 354 zdtdx(:,:) = 0._wp ; zdtdy(:,:) = 0._wp ; zdsdx(:,:) = 0._wp 355 zdsdy(:,:) = 0._wp ; dbdx_mle(:,:) = 0._wp ; dbdy_mle(:,:) = 0._wp 356 zwb_fk(:,:) = 0._wp ; zvel_mle(:,:) = 0._wp; zdiff_mle(:,:) = 0._wp 357 zhmle(:,:) = 0._wp ; zmld(:,:) = 0._wp 358 ENDIF 359 zwb_fk_b(:,:) = 0._wp ! must be initialised even with ln_osm_mle=F as used in zdf_osm_calculate_dhdt 360 279 361 ! Flux-Gradient arrays. 280 362 zdifml_sc(:,:) = 0._wp ; zvisml_sc(:,:) = 0._wp ; zdifpyc_sc(:,:) = 0._wp … … 287 369 ghams(:,:,:) = 0._wp ; ghamu(:,:,:) = 0._wp ; ghamv(:,:,:) = 0._wp 288 370 371 zdhdt_2(:,:) = 0._wp 289 372 ! hbl = MAX(hbl,epsln) 290 373 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 350 433 ! Use wind speed wndm included in sbc_oce module 351 434 zustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8 ) 352 dstokes(ji,jj) = 0.12 * wndm(ji,jj)**2 / grav435 dstokes(ji,jj) = MAX( 0.12 * wndm(ji,jj)**2 / grav, 5.e-1) 353 436 END DO 354 437 END DO … … 362 445 ! It could represent the effects of the spread of wave directions 363 446 ! around the mean wind. The effect of this adjustment needs to be tested. 364 z ustke(ji,jj) = MAX ( 1.0 * ( zcos_wind(ji,jj) * ut0sd(ji,jj ) + zsin_wind(ji,jj) * vt0sd(ji,jj) ), &365 & zustar(ji,jj) / ( 0.45 * 0.45 ))366 dstokes(ji,jj) = MAX(zfac * hsw(ji,jj)*hsw(ji,jj) / ( MAX(z ustke(ji,jj)*wmp(ji,jj), 1.0e-7 ) ), 5.0e-1) !rn_osm_dstokes !447 zabsstke = SQRT(ut0sd(ji,jj)**2 + vt0sd(ji,jj)**2) 448 zustke(ji,jj) = MAX (0.8 * ( zcos_wind(ji,jj) * ut0sd(ji,jj) + zsin_wind(ji,jj) * vt0sd(ji,jj) ), 1.0e-8) 449 dstokes(ji,jj) = MAX(zfac * hsw(ji,jj)*hsw(ji,jj) / ( MAX(zabsstke*wmp(ji,jj), 1.0e-7 ) ), 5.0e-1) !rn_osm_dstokes ! 367 450 END DO 368 451 END DO … … 375 458 ! Langmuir velocity scale (zwstrl), at T-point 376 459 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 460 zla(ji,jj) = MAX(MIN(SQRT ( zustar(ji,jj) / ( zwstrl(ji,jj) + epsln ) )**3, 4.0), 0.2) 461 IF(zla(ji,jj) > 0.45) dstokes(ji,jj) = MIN(dstokes(ji,jj), 0.5_wp*hbl(ji,jj)) 385 462 ! Velocity scale that tends to zustar for large Langmuir numbers 386 463 zvstr(ji,jj) = ( zwstrl(ji,jj)**3 + & … … 389 466 ! limit maximum value of Langmuir number as approximate treatment for shear turbulence. 390 467 ! Note zustke and zwstrl are not amended. 391 IF ( zla(ji,jj) >= 0.45 ) zla(ji,jj) = 0.45392 468 ! 393 469 ! get convective velocity (zwstrc), stabilty scale (zhol) and logical conection flag lconv … … 406 482 ! Mixed-layer model - calculate averages over the boundary layer, and the change in the boundary layer depth 407 483 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 408 ! BL must be always 2 levels deep. 409 hbl(:,:) = MAX(hbl(:,:), gdepw_n(:,:,3) ) 410 ibld(:,:) = 3 411 DO jk = 4, jpkm1 412 DO jj = 2, jpjm1 413 DO ji = 2, jpim1 484 ! BL must be always 4 levels deep. 485 ! For calculation of lateral buoyancy gradients for FK in 486 ! zdf_osm_zmld_horizontal_gradients need halo values for ibld, so must 487 ! previously exist for hbl also. 488 hbl(:,:) = MAX(hbl(:,:), gdepw_n(:,:,4) ) 489 ibld(:,:) = 4 490 DO jk = 5, jpkm1 491 DO jj = 1, jpj 492 DO ji = 1, jpi 414 493 IF ( hbl(ji,jj) >= gdepw_n(ji,jj,jk) ) THEN 415 494 ibld(ji,jj) = MIN(mbkt(ji,jj), jk) … … 419 498 END DO 420 499 421 DO jj = 2, jpjm1 ! Vertical slab500 DO jj = 2, jpjm1 422 501 DO ji = 2, jpim1 423 zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 424 zbeta = rab_n(ji,jj,1,jp_sal) 425 zt = 0._wp 426 zs = 0._wp 427 zu = 0._wp 428 zv = 0._wp 429 ! average over depth of boundary layer 430 zthick=0._wp 431 DO jm = 2, ibld(ji,jj) 432 zthick=zthick+e3t_n(ji,jj,jm) 433 zt = zt + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem) 434 zs = zs + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal) 435 zu = zu + e3t_n(ji,jj,jm) & 436 & * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) & 437 & / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 438 zv = zv + e3t_n(ji,jj,jm) & 439 & * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) & 440 & / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 441 END DO 442 zt_bl(ji,jj) = zt / zthick 443 zs_bl(ji,jj) = zs / zthick 444 zu_bl(ji,jj) = zu / zthick 445 zv_bl(ji,jj) = zv / zthick 446 zdt_bl(ji,jj) = zt_bl(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_tem) 447 zds_bl(ji,jj) = zs_bl(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal) 448 zdu_bl(ji,jj) = zu_bl(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) & 449 & / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 450 zdv_bl(ji,jj) = zv_bl(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) & 451 & / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 452 zdb_bl(ji,jj) = grav * zthermal * zdt_bl(ji,jj) - grav * zbeta * zds_bl(ji,jj) 453 IF ( lconv(ji,jj) ) THEN ! Convective 454 zwb_ent(ji,jj) = -( 2.0 * 0.2 * zwbav(ji,jj) & 455 & + 0.135 * zla(ji,jj) * zwstrl(ji,jj)**3/hbl(ji,jj) ) 456 457 zvel_max = - ( 1.0 + 1.0 * ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / hbl(ji,jj) ) & 458 & * zwb_ent(ji,jj) / ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 459 ! Entrainment including component due to shear turbulence. Modified Langmuir component, but gives same result for La=0.3 For testing uncomment. 460 ! zwb_ent(ji,jj) = -( 2.0 * 0.2 * zwbav(ji,jj) & 461 ! & + ( 0.15 * ( 1.0 - EXP( -0.5 * zla(ji,jj) ) ) + 0.03 / zla(ji,jj)**2 ) * zustar(ji,jj)**3/hbl(ji,jj) ) 462 463 ! zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / zhbl(ji,jj) ) * zwb_ent(ji,jj) / & 464 ! & ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 465 zzdhdt = - zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj),0.0) ) 466 ELSE ! Stable 467 zzdhdt = 0.32 * ( hbli(ji,jj) / hbl(ji,jj) -1.0 ) * zwstrl(ji,jj)**3 / hbli(ji,jj) & 468 & + ( ( 0.32 / 3.0 ) * exp ( -2.5 * ( hbli(ji,jj) / hbl(ji,jj) - 1.0 ) ) & 469 & - ( 0.32 / 3.0 - 0.135 * zla(ji,jj) ) * exp ( -12.5 * ( hbli(ji,jj) / hbl(ji,jj) ) ) ) & 470 & * zwstrl(ji,jj)**3 / hbli(ji,jj) 471 zzdhdt = zzdhdt + zwbav(ji,jj) 472 IF ( zzdhdt < 0._wp ) THEN 473 ! For long timsteps factor in brackets slows the rapid collapse of the OSBL 474 zpert = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_rdt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj) 475 ELSE 476 zpert = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_rdt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj) & 477 & + MAX( zdb_bl(ji,jj), 0.0 ) 478 ENDIF 479 zzdhdt = 2.0 * zzdhdt / zpert 480 ENDIF 481 zdhdt(ji,jj) = zzdhdt 482 END DO 502 zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj)) 503 imld(ji,jj) = MAX(3,ibld(ji,jj) - MAX( INT( dh(ji,jj) / e3t_n(ji, jj, ibld(ji,jj) )) , 1 )) 504 zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 505 END DO 483 506 END DO 484 507 ! Averages over well-mixed and boundary layer 508 ! Alan: do we need zb_nl?, zb_ml? 509 CALL zdf_osm_vertical_average(ibld, zt_bl, zs_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl) 510 CALL zdf_osm_vertical_average(imld, zt_ml, zs_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml) 511 ! External gradient 512 CALL zdf_osm_external_gradients( zdtdz_ext, zdsdz_ext, zdbdz_ext ) 513 514 515 IF ( ln_osm_mle ) THEN 516 CALL zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle ) 517 CALL zdf_osm_mle_parameters( hmle, zwb_fk, zvel_mle, zdiff_mle ) 518 ENDIF 519 520 ! Rate of change of hbl 521 CALL zdf_osm_calculate_dhdt( zdhdt, zdhdt_2 ) 485 522 ! Calculate averages over depth of boundary layer 486 imld = ibld ! use imld to hold previous blayer index 487 ibld(:,:) = 3 488 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 523 DO jj = 2, jpjm1 524 DO ji = 2, jpim1 525 zhbl_t(ji,jj) = hbl(ji,jj) + (zdhdt(ji,jj) - wn(ji,jj,ibld(ji,jj)))* rn_rdt ! certainly need wn here, so subtract it 526 ! adjustment to represent limiting by ocean bottom 527 zhbl_t(ji,jj) = MIN(zhbl_t(ji,jj), gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol)! ht_n(:,:)) 528 END DO 529 END DO 530 531 imld(:,:) = ibld(:,:) ! use imld to hold previous blayer index 532 ibld(:,:) = 4 492 533 493 534 DO jk = 4, jpkm1 … … 495 536 DO ji = 2, jpim1 496 537 IF ( zhbl_t(ji,jj) >= gdepw_n(ji,jj,jk) ) THEN 497 ibld(ji,jj) = MIN(mbkt(ji,jj), jk)538 ibld(ji,jj) = jk 498 539 ENDIF 499 540 END DO … … 504 545 ! Step through model levels taking account of buoyancy change to determine the effect on dhdt 505 546 ! 506 DO jj = 2, jpjm1507 DO ji = 2, jpim1508 IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN547 CALL zdf_osm_timestep_hbl( zdhdt, zdhdt_2 ) 548 ! Alan: do we need zb_ml? 549 CALL zdf_osm_vertical_average( ibld, zt_bl, zs_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl ) 509 550 ! 510 ! If boundary layer changes by more than one level, need to check for stable layers between initial and final depths.511 551 ! 512 zhbl_s = hbl(ji,jj) 513 jm = imld(ji,jj) 514 zthermal = rab_n(ji,jj,1,jp_tem) 515 zbeta = rab_n(ji,jj,1,jp_sal) 516 IF ( lconv(ji,jj) ) THEN 517 !unstable 518 zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / hbl(ji,jj) ) & 519 & * zwb_ent(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 520 521 DO jk = imld(ji,jj), ibld(ji,jj) 522 zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) ) & 523 & - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ), 0.0 ) + zvel_max 524 525 zhbl_s = zhbl_s + MIN( - zwb_ent(ji,jj) / zdb * rn_rdt / FLOAT(ibld(ji,jj)-imld(ji,jj) ), e3w_n(ji,jj,jk) ) 526 zhbl_s = MIN(zhbl_s, ht_n(ji,jj)) 527 528 IF ( zhbl_s >= gdepw_n(ji,jj,jm+1) ) jm = jm + 1 529 END DO 530 hbl(ji,jj) = zhbl_s 531 ibld(ji,jj) = jm 532 hbli(ji,jj) = hbl(ji,jj) 533 ELSE 534 ! stable 535 DO jk = imld(ji,jj), ibld(ji,jj) 536 zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) ) & 537 & - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ), 0.0 ) & 538 & + 2.0 * zwstrl(ji,jj)**2 / zhbl_s 539 540 zhbl_s = zhbl_s + ( & 541 & 0.32 * ( hbli(ji,jj) / zhbl_s -1.0 ) & 542 & * zwstrl(ji,jj)**3 / hbli(ji,jj) & 543 & + ( ( 0.32 / 3.0 ) * EXP( - 2.5 * ( hbli(ji,jj) / zhbl_s -1.0 ) ) & 544 & - ( 0.32 / 3.0 - 0.0485 ) * EXP( - 12.5 * ( hbli(ji,jj) / zhbl_s ) ) ) & 545 & * zwstrl(ji,jj)**3 / hbli(ji,jj) ) / zdb * e3w_n(ji,jj,jk) / zdhdt(ji,jj) ! ALMG to investigate whether need to include wn here 546 547 zhbl_s = MIN(zhbl_s, ht_n(ji,jj)) 548 IF ( zhbl_s >= gdepw_n(ji,jj,jm) ) jm = jm + 1 549 END DO 550 hbl(ji,jj) = MAX(zhbl_s, gdepw_n(ji,jj,3) ) 551 ibld(ji,jj) = MAX(jm, 3 ) 552 IF ( hbl(ji,jj) > hbli(ji,jj) ) hbli(ji,jj) = hbl(ji,jj) 553 ENDIF ! IF ( lconv ) 554 ELSE 555 ! change zero or one model level. 556 hbl(ji,jj) = zhbl_t(ji,jj) 557 IF ( lconv(ji,jj) ) THEN 558 hbli(ji,jj) = hbl(ji,jj) 559 ELSE 560 hbl(ji,jj) = MAX(hbl(ji,jj), gdepw_n(ji,jj,3) ) 561 IF ( hbl(ji,jj) > hbli(ji,jj) ) hbli(ji,jj) = hbl(ji,jj) 562 ENDIF 563 ENDIF 564 zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj)) 565 END DO 566 END DO 552 CALL zdf_osm_pycnocline_thickness( dh, zdh ) 567 553 dstokes(:,:) = MIN ( dstokes(:,:), hbl(:,:)/3. ) ! Limit delta for shallow boundary layers for calculating flux-gradient terms. 568 569 ! Recalculate averages over boundary layer after depth updated 570 ! Consider later combining this into the loop above and looking for columns 571 ! where the index for base of the boundary layer have changed 572 DO jj = 2, jpjm1 ! Vertical slab 573 DO ji = 2, jpim1 574 zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 575 zbeta = rab_n(ji,jj,1,jp_sal) 576 zt = 0._wp 577 zs = 0._wp 578 zu = 0._wp 579 zv = 0._wp 580 ! average over depth of boundary layer 581 zthick=0._wp 582 DO jm = 2, ibld(ji,jj) 583 zthick=zthick+e3t_n(ji,jj,jm) 584 zt = zt + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem) 585 zs = zs + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal) 586 zu = zu + e3t_n(ji,jj,jm) & 587 & * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) & 588 & / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 589 zv = zv + e3t_n(ji,jj,jm) & 590 & * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) & 591 & / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 592 END DO 593 zt_bl(ji,jj) = zt / zthick 594 zs_bl(ji,jj) = zs / zthick 595 zu_bl(ji,jj) = zu / zthick 596 zv_bl(ji,jj) = zv / zthick 597 zdt_bl(ji,jj) = zt_bl(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_tem) 598 zds_bl(ji,jj) = zs_bl(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal) 599 zdu_bl(ji,jj) = zu_bl(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) & 600 & / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 601 zdv_bl(ji,jj) = zv_bl(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) & 602 & / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 603 zdb_bl(ji,jj) = grav * zthermal * zdt_bl(ji,jj) - grav * zbeta * zds_bl(ji,jj) 604 zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj)) 605 IF ( lconv(ji,jj) ) THEN 606 IF ( zdb_bl(ji,jj) > 0._wp )THEN 607 IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN ! near neutral stability 608 zari = 4.5 * ( zvstr(ji,jj)**2 ) & 609 & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01 610 ELSE ! unstable 611 zari = 4.5 * ( zwstrc(ji,jj)**2 ) & 612 & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01 613 ENDIF 614 IF ( zari > 0.2 ) THEN ! This test checks for weakly stratified pycnocline 615 zari = 0.2 616 zwb_ent(ji,jj) = 0._wp 617 ENDIF 618 inhml = MAX( INT( zari * zhbl(ji,jj) / e3t_n(ji,jj,ibld(ji,jj)) ) , 1 ) 619 imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 1) 620 zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 621 zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 622 ELSE ! IF (zdb_bl) 623 imld(ji,jj) = ibld(ji,jj) - 1 624 zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 625 zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 626 ENDIF 627 ELSE ! IF (lconv) 628 IF ( zdhdt(ji,jj) >= 0.0 ) THEN ! probably shouldn't include wm here 629 ! boundary layer deepening 630 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 631 ! pycnocline thickness set by stratification - use same relationship as for neutral conditions. 632 zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) & 633 & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01 , 0.2 ) 634 inhml = MAX( INT( zari * zhbl(ji,jj) / e3t_n(ji,jj,ibld(ji,jj)) ) , 1 ) 635 imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 1) 636 zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 637 zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 638 ELSE 639 imld(ji,jj) = ibld(ji,jj) - 1 640 zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 641 zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 642 ENDIF ! IF (zdb_bl > 0.0) 643 ELSE ! IF(dhdt >= 0) 644 ! boundary layer collapsing. 645 imld(ji,jj) = ibld(ji,jj) 646 zhml(ji,jj) = zhbl(ji,jj) 647 zdh(ji,jj) = 0._wp 648 ENDIF ! IF (dhdt >= 0) 649 ENDIF ! IF (lconv) 650 END DO 651 END DO 652 653 ! Average over the depth of the mixed layer in the convective boundary layer 654 ! Also calculate entrainment fluxes for temperature and salinity 655 DO jj = 2, jpjm1 ! Vertical slab 656 DO ji = 2, jpim1 657 zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 658 zbeta = rab_n(ji,jj,1,jp_sal) 659 IF ( lconv(ji,jj) ) THEN 660 zt = 0._wp 661 zs = 0._wp 662 zu = 0._wp 663 zv = 0._wp 664 ! average over depth of boundary layer 665 zthick=0._wp 666 DO jm = 2, imld(ji,jj) 667 zthick=zthick+e3t_n(ji,jj,jm) 668 zt = zt + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem) 669 zs = zs + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal) 670 zu = zu + e3t_n(ji,jj,jm) & 671 & * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) & 672 & / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 673 zv = zv + e3t_n(ji,jj,jm) & 674 & * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) & 675 & / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 676 END DO 677 zt_ml(ji,jj) = zt / zthick 678 zs_ml(ji,jj) = zs / zthick 679 zu_ml(ji,jj) = zu / zthick 680 zv_ml(ji,jj) = zv / zthick 681 zdt_ml(ji,jj) = zt_ml(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_tem) 682 zds_ml(ji,jj) = zs_ml(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal) 683 zdu_ml(ji,jj) = zu_ml(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) & 684 & / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 685 zdv_ml(ji,jj) = zv_ml(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) & 686 & / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 687 zdb_ml(ji,jj) = grav * zthermal * zdt_ml(ji,jj) - grav * zbeta * zds_ml(ji,jj) 688 ELSE 689 ! stable, if entraining calulate average below interface layer. 690 IF ( zdhdt(ji,jj) >= 0._wp ) THEN 691 zt = 0._wp 692 zs = 0._wp 693 zu = 0._wp 694 zv = 0._wp 695 ! average over depth of boundary layer 696 zthick=0._wp 697 DO jm = 2, imld(ji,jj) 698 zthick=zthick+e3t_n(ji,jj,jm) 699 zt = zt + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem) 700 zs = zs + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal) 701 zu = zu + e3t_n(ji,jj,jm) & 702 & * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) & 703 & / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 704 zv = zv + e3t_n(ji,jj,jm) & 705 & * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) & 706 & / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 707 END DO 708 zt_ml(ji,jj) = zt / zthick 709 zs_ml(ji,jj) = zs / zthick 710 zu_ml(ji,jj) = zu / zthick 711 zv_ml(ji,jj) = zv / zthick 712 zdt_ml(ji,jj) = zt_ml(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_tem) 713 zds_ml(ji,jj) = zs_ml(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal) 714 zdu_ml(ji,jj) = zu_ml(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) & 715 & / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 716 zdv_ml(ji,jj) = zv_ml(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) & 717 & / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 718 zdb_ml(ji,jj) = grav * zthermal * zdt_ml(ji,jj) - grav * zbeta * zds_ml(ji,jj) 719 ENDIF 720 ENDIF 721 END DO 722 END DO 723 ! 554 ! 555 ! Average over the depth of the mixed layer in the convective boundary layer 556 ! Alan: do we need zb_ml? 557 CALL zdf_osm_vertical_average( imld, zt_ml, zs_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml ) 724 558 ! rotate mean currents and changes onto wind align co-ordinates 725 559 ! 726 727 DO jj = 2, jpjm1 728 DO ji = 2, jpim1 729 ztemp = zu_ml(ji,jj) 730 zu_ml(ji,jj) = zu_ml(ji,jj) * zcos_wind(ji,jj) + zv_ml(ji,jj) * zsin_wind(ji,jj) 731 zv_ml(ji,jj) = zv_ml(ji,jj) * zcos_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 732 ztemp = zdu_ml(ji,jj) 733 zdu_ml(ji,jj) = zdu_ml(ji,jj) * zcos_wind(ji,jj) + zdv_ml(ji,jj) * zsin_wind(ji,jj) 734 zdv_ml(ji,jj) = zdv_ml(ji,jj) * zsin_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 735 ! 736 ztemp = zu_bl(ji,jj) 737 zu_bl = zu_bl(ji,jj) * zcos_wind(ji,jj) + zv_bl(ji,jj) * zsin_wind(ji,jj) 738 zv_bl(ji,jj) = zv_bl(ji,jj) * zcos_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 739 ztemp = zdu_bl(ji,jj) 740 zdu_bl(ji,jj) = zdu_bl(ji,jj) * zcos_wind(ji,jj) + zdv_bl(ji,jj) * zsin_wind(ji,jj) 741 zdv_bl(ji,jj) = zdv_bl(ji,jj) * zsin_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 742 END DO 743 END DO 744 560 CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_ml, zv_ml, zdu_ml, zdv_ml ) 561 CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_bl, zv_bl, zdu_bl, zdv_bl ) 745 562 zuw_bse = 0._wp 746 563 zvw_bse = 0._wp 564 zwth_ent = 0._wp 565 zws_ent = 0._wp 747 566 DO jj = 2, jpjm1 748 567 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) 568 IF ( ibld(ji,jj) < mbkt(ji,jj) ) THEN 569 IF ( lconv(ji,jj) ) THEN 570 zuw_bse(ji,jj) = -0.0075*((zvstr(ji,jj)**3+0.5*zwstrc(ji,jj)**3)**pthird*zdu_ml(ji,jj) + & 571 & 1.5*zustar(ji,jj)**2*(zhbl(ji,jj)-zhml(ji,jj)) )/ & 572 & ( zhml(ji,jj)*MIN(zla(ji,jj)**(8./3.),1.) + epsln) 573 zvw_bse(ji,jj) = 0.01*(-(zvstr(ji,jj)**3+0.5*zwstrc(ji,jj)**3)**pthird*zdv_ml(ji,jj)+ & 574 & 2.0*ff_t(ji,jj)*zustke(ji,jj)*dstokes(ji,jj)*zla(ji,jj)) 575 IF ( zdb_ml(ji,jj) > 0._wp ) THEN 576 zwth_ent(ji,jj) = zwb_ent(ji,jj) * zdt_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 577 zws_ent(ji,jj) = zwb_ent(ji,jj) * zds_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 578 ENDIF 579 ELSE 580 zwth_ent(ji,jj) = -2.0 * zwthav(ji,jj) * ( (1.0 - 0.8) - ( 1.0 - 0.8)**(3.0/2.0) ) 581 zws_ent(ji,jj) = -2.0 * zwsav(ji,jj) * ( (1.0 - 0.8 ) - ( 1.0 - 0.8 )**(3.0/2.0) ) 754 582 ENDIF 755 ELSE756 zwth_ent(ji,jj) = -2.0 * zwthav(ji,jj) * ( (1.0 - 0.8) - ( 1.0 - 0.8)**(3.0/2.0) )757 zws_ent(ji,jj) = -2.0 * zwsav(ji,jj) * ( (1.0 - 0.8 ) - ( 1.0 - 0.8 )**(3.0/2.0) )758 583 ENDIF 759 584 END DO … … 764 589 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 765 590 766 DO jj = 2, jpjm1 767 DO ji = 2, jpim1 768 ! 769 IF ( lconv (ji,jj) ) THEN 770 ! Unstable conditions 771 IF( zdb_bl(ji,jj) > 0._wp ) THEN 772 ! calculate pycnocline profiles, no need if zdb_bl <= 0. since profile is zero and arrays have been initialized to zero 773 ztgrad = ( zdt_ml(ji,jj) / zdh(ji,jj) ) 774 zsgrad = ( zds_ml(ji,jj) / zdh(ji,jj) ) 775 zbgrad = ( zdb_ml(ji,jj) / zdh(ji,jj) ) 776 DO jk = 2 , ibld(ji,jj) 777 znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 778 zdtdz_pyc(ji,jj,jk) = ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 779 zdbdz_pyc(ji,jj,jk) = zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 780 zdsdz_pyc(ji,jj,jk) = zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 781 END DO 782 ENDIF 783 ELSE 784 ! stable conditions 785 ! if pycnocline profile only defined when depth steady of increasing. 786 IF ( zdhdt(ji,jj) >= 0.0 ) THEN ! Depth increasing, or steady. 787 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 788 IF ( zhol(ji,jj) >= 0.5 ) THEN ! Very stable - 'thick' pycnocline 789 ztgrad = zdt_bl(ji,jj) / zhbl(ji,jj) 790 zsgrad = zds_bl(ji,jj) / zhbl(ji,jj) 791 zbgrad = zdb_bl(ji,jj) / zhbl(ji,jj) 792 DO jk = 2, ibld(ji,jj) 793 znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 794 zdtdz_pyc(ji,jj,jk) = ztgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 795 zdbdz_pyc(ji,jj,jk) = zbgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 796 zdsdz_pyc(ji,jj,jk) = zsgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 797 END DO 798 ELSE ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form. 799 ztgrad = zdt_bl(ji,jj) / zdh(ji,jj) 800 zsgrad = zds_bl(ji,jj) / zdh(ji,jj) 801 zbgrad = zdb_bl(ji,jj) / zdh(ji,jj) 802 DO jk = 2, ibld(ji,jj) 803 znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 804 zdtdz_pyc(ji,jj,jk) = ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 805 zdbdz_pyc(ji,jj,jk) = zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 806 zdsdz_pyc(ji,jj,jk) = zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 807 END DO 808 ENDIF ! IF (zhol >=0.5) 809 ENDIF ! IF (zdb_bl> 0.) 810 ENDIF ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero, profile arrays are intialized to zero 811 ENDIF ! IF (lconv) 812 ! 813 END DO 814 END DO 815 ! 816 DO jj = 2, jpjm1 817 DO ji = 2, jpim1 818 ! 819 IF ( lconv (ji,jj) ) THEN 820 ! Unstable conditions 821 zugrad = ( zdu_ml(ji,jj) / zdh(ji,jj) ) + 0.275 * zustar(ji,jj)*zustar(ji,jj) / & 822 & (( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) / zla(ji,jj)**(8.0/3.0) 823 zvgrad = ( zdv_ml(ji,jj) / zdh(ji,jj) ) + 3.5 * ff_t(ji,jj) * zustke(ji,jj) / & 824 & ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 825 DO jk = 2 , ibld(ji,jj)-1 826 znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 827 zdudz_pyc(ji,jj,jk) = zugrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 828 zdvdz_pyc(ji,jj,jk) = zvgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 829 END DO 830 ELSE 831 ! stable conditions 832 zugrad = 3.25 * zdu_bl(ji,jj) / zhbl(ji,jj) 833 zvgrad = 2.75 * zdv_bl(ji,jj) / zhbl(ji,jj) 834 DO jk = 2, ibld(ji,jj) 835 znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 836 IF ( znd < 1.0 ) THEN 837 zdudz_pyc(ji,jj,jk) = zugrad * EXP( -40.0 * ( znd - 1.0 )**2 ) 838 ELSE 839 zdudz_pyc(ji,jj,jk) = zugrad * EXP( -20.0 * ( znd - 1.0 )**2 ) 840 ENDIF 841 zdvdz_pyc(ji,jj,jk) = zvgrad * EXP( -20.0 * ( znd - 0.85 )**2 ) 842 END DO 843 ENDIF 844 ! 845 END DO 846 END DO 591 CALL zdf_osm_external_gradients( zdtdz_ext, zdsdz_ext, zdbdz_ext ) 592 CALL zdf_osm_pycnocline_scalar_profiles( zdtdz_pyc, zdsdz_pyc, zdbdz_pyc ) 593 CALL zdf_osm_pycnocline_shear_profiles( zdudz_pyc, zdvdz_pyc ) 847 594 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 848 595 ! Eddy viscosity/diffusivity and non-gradient terms in the flux-gradient relationship 849 596 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 850 597 851 ! WHERE ( lconv )852 ! zdifml_sc = zhml * ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird853 ! zvisml_sc = zdifml_sc854 ! zdifpyc_sc = 0.165 * ( zwstrl**3 + zwstrc**3 )**pthird * ( zhbl - zhml )855 ! zvispyc_sc = 0.142 * ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * ( zhbl - zhml )856 ! zbeta_d_sc = 1.0 - (0.165 / 0.8 * ( zhbl - zhml ) / zhbl )**p2third857 ! zbeta_v_sc = 1.0 - 2.0 * (0.142 /0.375) * (zhbl - zhml ) / zhml858 ! ELSEWHERE859 ! zdifml_sc = zwstrl * zhbl * EXP ( -( zhol / 0.183_wp )**2 )860 ! zvisml_sc = zwstrl * zhbl * EXP ( -( zhol / 0.183_wp )**2 )861 ! ENDWHERE862 598 DO jj = 2, jpjm1 863 599 DO ji = 2, jpim1 … … 868 604 zvispyc_sc(ji,jj) = 0.142 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zdh(ji,jj) 869 605 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)606 zbeta_v_sc(ji,jj) = 1.0 - 2.0 * (0.142 /0.375) * zdh(ji,jj) / ( zhml(ji,jj) + epsln ) 871 607 ELSE 872 608 zdifml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ) 873 609 zvisml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ) 874 END IF875 END DO876 END DO610 END IF 611 END DO 612 END DO 877 613 ! 878 614 DO jj = 2, jpjm1 … … 889 625 ! pycnocline - if present linear profile 890 626 IF ( zdh(ji,jj) > 0._wp ) THEN 627 zgamma_b = 6.0 891 628 DO jk = imld(ji,jj)+1 , ibld(ji,jj) 892 629 zznd_pyc = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 893 630 ! 894 zdiffut(ji,jj,jk) = zdifpyc_sc(ji,jj) * ( 1.0 +zznd_pyc )631 zdiffut(ji,jj,jk) = zdifpyc_sc(ji,jj) * EXP( zgamma_b * zznd_pyc ) 895 632 ! 896 zviscos(ji,jj,jk) = zvispyc_sc(ji,jj) * ( 1.0 +zznd_pyc )633 zviscos(ji,jj,jk) = zvispyc_sc(ji,jj) * EXP( zgamma_b * zznd_pyc ) 897 634 END DO 635 IF ( ibld_ext == 0 ) THEN 636 zdiffut(ji,jj,ibld(ji,jj)) = 0._wp 637 zviscos(ji,jj,ibld(ji,jj)) = 0._wp 638 ELSE 639 zdiffut(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj) * ( hbl(ji,jj) - gdepw_n(ji, jj, ibld(ji,jj)-1) ) 640 zviscos(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj) * ( hbl(ji,jj) - gdepw_n(ji, jj, ibld(ji,jj)-1) ) 641 ENDIF 898 642 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)) 643 ! Temporary fix to ensure zdiffut is +ve; won't be necessary with wn taken out 644 zdiffut(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj) * e3t_n(ji,jj,ibld(ji,jj)), 1.e-6) 645 zviscos(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj) * e3t_n(ji,jj,ibld(ji,jj)), 1.e-5) 901 646 ! could be taken out, take account of entrainment represents as a diffusivity 902 647 ! should remove w from here, represents entrainment … … 908 653 zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * (1.0 - zznd_ml) * ( 1.0 - zznd_ml**2 ) 909 654 END DO 655 656 IF ( ibld_ext == 0 ) THEN 657 zdiffut(ji,jj,ibld(ji,jj)) = 0._wp 658 zviscos(ji,jj,ibld(ji,jj)) = 0._wp 659 ELSE 660 zdiffut(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 0._wp) * e3w_n(ji, jj, ibld(ji,jj)) 661 zviscos(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 0._wp) * e3w_n(ji, jj, ibld(ji,jj)) 662 ENDIF 910 663 ENDIF ! end if ( lconv ) 911 !664 ! 912 665 END DO ! end of ji loop 913 666 END DO ! end of jj loop … … 952 705 END DO ! end of jj loop 953 706 954 955 707 ! Stokes term in flux-gradient relationship (note in zsc_uw_n don't use zvstr since term needs to go to zero as zwstrl goes to zero) 956 708 WHERE ( lconv ) 957 zsc_uw_1 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke / ( 1.0 - 1.0 * 6.5 * zla**(8.0/3.0))958 zsc_uw_2 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke / ( zla**(8.0/3.0) + epsln)959 zsc_vw_1 = ff_t * zhml * zustke**3 * zla**(8.0/3.0) / ( ( zvstr**3 + 0.5 * zwstrc**3 )**(2.0/3.0) + epsln )709 zsc_uw_1 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke / MAX( ( 1.0 - 1.0 * 6.5 * zla**(8.0/3.0) ), 0.2 ) 710 zsc_uw_2 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke / MIN( zla**(8.0/3.0) + epsln, 0.12 ) 711 zsc_vw_1 = ff_t * zhml * zustke**3 * MIN( zla**(8.0/3.0), 0.12 ) / ( ( zvstr**3 + 0.5 * zwstrc**3 )**(2.0/3.0) + epsln ) 960 712 ELSEWHERE 961 713 zsc_uw_1 = zustar**2 962 zsc_vw_1 = ff_t * zhbl * zustke**3 * zla**(8.0/3.0) / (zvstr**2 + epsln)714 zsc_vw_1 = ff_t * zhbl * zustke**3 * MIN( zla**(8.0/3.0), 0.12 ) / (zvstr**2 + epsln) 963 715 ENDWHERE 964 716 IF(ln_dia_osm) THEN 717 IF ( iom_use("ghamu_00") ) CALL iom_put( "ghamu_00", wmask*ghamu ) 718 IF ( iom_use("ghamv_00") ) CALL iom_put( "ghamv_00", wmask*ghamv ) 719 END IF 965 720 DO jj = 2, jpjm1 966 721 DO ji = 2, jpim1 … … 1007 762 zl_l = 2.0 * ( 1.0 - EXP ( - 2.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) ) & 1008 763 & * ( 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)764 zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0 + EXP ( -3.0 * LOG10 ( - zhol(ji,jj) ) ) ) ** (3.0/2.0) 1010 765 ! non-gradient buoyancy terms 1011 766 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * 0.5 * zsc_wth_1(ji,jj) * zl_eps * zhml(ji,jj) / ( 0.15 + zznd_ml ) … … 1020 775 END DO ! ji loop 1021 776 END DO ! jj loop 1022 1023 777 1024 778 WHERE ( lconv ) … … 1051 805 END DO ! jj loop 1052 806 807 IF(ln_dia_osm) THEN 808 IF ( iom_use("ghamu_0") ) CALL iom_put( "ghamu_0", wmask*ghamu ) 809 IF ( iom_use("zsc_uw_1_0") ) CALL iom_put( "zsc_uw_1_0", tmask(:,:,1)*zsc_uw_1 ) 810 END IF 1053 811 ! Transport term in flux-gradient relationship [note : includes ROI ratio (X0.3) ] 1054 812 … … 1089 847 END DO ! jj loop 1090 848 1091 1092 849 WHERE ( lconv ) 1093 850 zsc_uw_1 = zustar**2 … … 1134 891 END DO 1135 892 END DO 893 894 IF(ln_dia_osm) THEN 895 IF ( iom_use("ghamu_f") ) CALL iom_put( "ghamu_f", wmask*ghamu ) 896 IF ( iom_use("ghamv_f") ) CALL iom_put( "ghamv_f", wmask*ghamv ) 897 IF ( iom_use("zsc_uw_1_f") ) CALL iom_put( "zsc_uw_1_f", tmask(:,:,1)*zsc_uw_1 ) 898 IF ( iom_use("zsc_vw_1_f") ) CALL iom_put( "zsc_vw_1_f", tmask(:,:,1)*zsc_vw_1 ) 899 IF ( iom_use("zsc_uw_2_f") ) CALL iom_put( "zsc_uw_2_f", tmask(:,:,1)*zsc_uw_2 ) 900 IF ( iom_use("zsc_vw_2_f") ) CALL iom_put( "zsc_vw_2_f", tmask(:,:,1)*zsc_vw_2 ) 901 END IF 1136 902 ! 1137 903 ! Make surface forced velocity non-gradient terms go to zero at the base of the mixed layer. … … 1165 931 END DO 1166 932 933 IF(ln_dia_osm) THEN 934 IF ( iom_use("ghamu_b") ) CALL iom_put( "ghamu_b", wmask*ghamu ) 935 IF ( iom_use("ghamv_b") ) CALL iom_put( "ghamv_b", wmask*ghamv ) 936 END IF 1167 937 ! pynocline contributions 1168 938 ! Temporary fix to avoid instabilities when zdb_bl becomes very very small … … 1170 940 DO jj = 2, jpjm1 1171 941 DO ji = 2, jpim1 1172 DO jk= 2, ibld(ji,jj) 1173 znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 1174 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk) 1175 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk) 1176 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc(ji,jj,jk) 1177 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zsc_uw_1(ji,jj) * ( 1.0 - znd )**(7.0/4.0) * zdbdz_pyc(ji,jj,jk) 1178 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc(ji,jj,jk) 1179 END DO 1180 END DO 942 IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN 943 DO jk= 2, ibld(ji,jj) 944 znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 945 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk) 946 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk) 947 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc(ji,jj,jk) 948 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zsc_uw_1(ji,jj) * ( 1.0 - znd )**(7.0/4.0) * zdbdz_pyc(ji,jj,jk) 949 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc(ji,jj,jk) 950 END DO 951 END IF 952 END DO 1181 953 END DO 1182 954 … … 1185 957 DO jj=2, jpjm1 1186 958 DO ji = 2, jpim1 1187 IF ( lconv(ji,jj)) THEN959 IF ( lconv(ji,jj) .AND. ibld(ji,jj) + ibld_ext < mbkt(ji,jj)) THEN 1188 960 DO jk = 1, imld(ji,jj) - 1 1189 961 znd=gdepw_n(ji,jj,jk) / zhml(ji,jj) 1190 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * znd1191 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * znd962 ! ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * znd 963 ! ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * znd 1192 964 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * znd 1193 965 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * znd … … 1195 967 DO jk = imld(ji,jj), ibld(ji,jj) 1196 968 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 )969 ! ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * ( 1.0 + znd ) 970 ! ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * ( 1.0 + znd ) 1199 971 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * ( 1.0 + znd ) 1200 972 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * ( 1.0 + znd ) 1201 973 END DO 1202 974 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 975 976 ghamt(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 977 ghams(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 978 ghamu(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 979 ghamv(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 1207 980 END DO ! ji loop 1208 981 END DO ! jj loop 1209 982 1210 983 IF(ln_dia_osm) THEN 984 IF ( iom_use("ghamu_1") ) CALL iom_put( "ghamu_1", wmask*ghamu ) 985 IF ( iom_use("ghamv_1") ) CALL iom_put( "ghamv_1", wmask*ghamv ) 986 IF ( iom_use("zuw_bse") ) CALL iom_put( "zuw_bse", tmask(:,:,1)*zuw_bse ) 987 IF ( iom_use("zvw_bse") ) CALL iom_put( "zvw_bse", tmask(:,:,1)*zvw_bse ) 988 IF ( iom_use("zdudz_pyc") ) CALL iom_put( "zdudz_pyc", wmask*zdudz_pyc ) 989 IF ( iom_use("zdvdz_pyc") ) CALL iom_put( "zdvdz_pyc", wmask*zdvdz_pyc ) 990 IF ( iom_use("zviscos") ) CALL iom_put( "zviscos", wmask*zviscos ) 991 END IF 1211 992 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 1212 993 ! Need to put in code for contributions that are applied explicitly to … … 1232 1013 IF(ln_dia_osm) THEN 1233 1014 IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask*zdtdz_pyc ) 1015 IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask*zdsdz_pyc ) 1016 IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask*zdbdz_pyc ) 1234 1017 END IF 1235 1018 … … 1286 1069 END IF ! ln_convmix = .true. 1287 1070 1071 1072 1073 IF ( ln_osm_mle ) THEN ! set up diffusivity and non-gradient mixing 1074 DO jj = 2 , jpjm1 1075 DO ji = 2, jpim1 1076 IF ( lconv(ji,jj) .AND. mld_prof(ji,jj) - ibld(ji,jj) > 1 ) THEN ! MLE mmixing extends below the OSBL. 1077 ! Calculate MLE flux profiles 1078 ! DO jk = 1, mld_prof(ji,jj) 1079 ! znd = - gdepw_n(ji,jj,jk) / MAX(zhmle(ji,jj),epsln) 1080 ! ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + & 1081 ! & zwt_fk(ji,jj) * ( 1.0 - ( 2.0 * znd + 1.0 )**2 ) * ( 1.0 + r5_21 * ( 2.0 * znd + 1.0 )**2 ) 1082 ! ghams(ji,jj,jk) = ghams(ji,jj,jk) + & 1083 ! & zws_fk(ji,jj) * ( 1.0 - ( 2.0 * znd + 1.0 )**2 ) * ( 1.0 + r5_21 * ( 2.0 * znd + 1.0 )**2 ) 1084 ! END DO 1085 ! Calculate MLE flux contribution from surface fluxes 1086 DO jk = 1, ibld(ji,jj) 1087 znd = gdepw_n(ji,jj,jk) / MAX(zhbl(ji,jj),epsln) 1088 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) - zwth0(ji,jj) * ( 1.0 - znd ) 1089 ghams(ji,jj,jk) = ghams(ji,jj,jk) - zws0(ji,jj) * ( 1.0 - znd ) 1090 END DO 1091 DO jk = 1, mld_prof(ji,jj) 1092 znd = gdepw_n(ji,jj,jk) / MAX(zhmle(ji,jj),epsln) 1093 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth0(ji,jj) * ( 1.0 - znd ) 1094 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws0(ji,jj) * ( 1.0 -znd ) 1095 END DO 1096 ! Viscosity for MLEs 1097 DO jk = ibld(ji,jj), mld_prof(ji,jj) 1098 zdiffut(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), zdiff_mle(ji,jj) ) 1099 END DO 1100 ! Iterate to find approx vertical index for depth 1.1*zhmle(ji,jj) 1101 jl = MIN(mld_prof(ji,jj) + 2, mbkt(ji,jj)) 1102 jl = MIN( MAX(INT( 0.1 * zhmle(ji,jj) / e3t_n(ji,jj,jl)), 2 ) + mld_prof(ji,jj), mbkt(ji,jj)) 1103 DO jk = mld_prof(ji,jj), jl 1104 zdiffut(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), zdiff_mle(ji,jj) * & 1105 & ( gdepw_n(ji,jj,jk) - gdepw_n(ji,jj,jl) ) / & 1106 & ( gdepw_n(ji,jj,mld_prof(ji,jj)) - gdepw_n(ji,jj,jl) - epsln)) 1107 END DO 1108 ENDIF 1109 END DO 1110 END DO 1111 ENDIF 1112 1113 IF(ln_dia_osm) THEN 1114 IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask*zdtdz_pyc ) 1115 IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask*zdsdz_pyc ) 1116 IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask*zdbdz_pyc ) 1117 END IF 1118 1119 1288 1120 ! Lateral boundary conditions on zvicos (sign unchanged), needed to caclulate viscosities on u and v grids 1289 CALL lbc_lnk( 'zdfosm',zviscos(:,:,:), 'W', 1. )1121 !CALL lbc_lnk( zviscos(:,:,:), 'W', 1. ) 1290 1122 1291 1123 ! GN 25/8: need to change tmask --> wmask … … 1316 1148 END DO 1317 1149 END DO 1150 ! Lateral boundary conditions on final outputs for hbl, on T-grid (sign unchanged) 1151 CALL lbc_lnk_multi( 'zdfosm', hbl, 'T', 1., dh, 'T', 1., hmle, 'T', 1. ) 1318 1152 ! Lateral boundary conditions on final outputs for gham[ts], on W-grid (sign unchanged) 1319 1153 ! Lateral boundary conditions on final outputs for gham[uv], on [UV]-grid (sign unchanged) 1320 1154 CALL lbc_lnk_multi( 'zdfosm', ghamt, 'W', 1. , ghams, 'W', 1., & 1321 & ghamu, 'U', 1. , ghamv, 'V',1. )1322 1323 1155 & ghamu, 'U', -1. , ghamv, 'V', -1. ) 1156 1157 IF(ln_dia_osm) THEN 1324 1158 SELECT CASE (nn_osm_wave) 1325 1159 ! Stokes drift set by assumimg onstant La#=0.3(=0) or Pierson-Moskovitz spectrum (=1). … … 1330 1164 ! Stokes drift read in from sbcwave (=2). 1331 1165 CASE(2) 1332 IF ( iom_use("us_x") ) CALL iom_put( "us_x", ut0sd ) ! x surface Stokes drift 1333 IF ( iom_use("us_y") ) CALL iom_put( "us_y", vt0sd ) ! y surface Stokes drift 1166 IF ( iom_use("us_x") ) CALL iom_put( "us_x", ut0sd*umask(:,:,1) ) ! x surface Stokes drift 1167 IF ( iom_use("us_y") ) CALL iom_put( "us_y", vt0sd*vmask(:,:,1) ) ! y surface Stokes drift 1168 IF ( iom_use("wmp") ) CALL iom_put( "wmp", wmp*tmask(:,:,1) ) ! wave mean period 1169 IF ( iom_use("hsw") ) CALL iom_put( "hsw", hsw*tmask(:,:,1) ) ! significant wave height 1170 IF ( iom_use("wmp_NP") ) CALL iom_put( "wmp_NP", (2.*rpi*1.026/(0.877*grav) )*wndm*tmask(:,:,1) ) ! wave mean period from NP spectrum 1171 IF ( iom_use("hsw_NP") ) CALL iom_put( "hsw_NP", (0.22/grav)*wndm**2*tmask(:,:,1) ) ! significant wave height from NP spectrum 1172 IF ( iom_use("wndm") ) CALL iom_put( "wndm", wndm*tmask(:,:,1) ) ! U_10 1334 1173 IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rau0*tmask(:,:,1)*zustar**2* & 1335 1174 & SQRT(ut0sd**2 + vt0sd**2 ) ) … … 1342 1181 IF ( iom_use("zws0") ) CALL iom_put( "zws0", tmask(:,:,1)*zws0 ) ! <Sw_0> 1343 1182 IF ( iom_use("hbl") ) CALL iom_put( "hbl", tmask(:,:,1)*hbl ) ! boundary-layer depth 1344 IF ( iom_use("hbli") ) CALL iom_put( "hbli", tmask(:,:,1)*hbli ) ! Initial boundary-layer depth 1183 IF ( iom_use("ibld") ) CALL iom_put( "ibld", tmask(:,:,1)*ibld ) ! boundary-layer max k 1184 IF ( iom_use("zdt_bl") ) CALL iom_put( "zdt_bl", tmask(:,:,1)*zdt_bl ) ! dt at ml base 1185 IF ( iom_use("zds_bl") ) CALL iom_put( "zds_bl", tmask(:,:,1)*zds_bl ) ! ds at ml base 1186 IF ( iom_use("zdb_bl") ) CALL iom_put( "zdb_bl", tmask(:,:,1)*zdb_bl ) ! db at ml base 1187 IF ( iom_use("zdu_bl") ) CALL iom_put( "zdu_bl", tmask(:,:,1)*zdu_bl ) ! du at ml base 1188 IF ( iom_use("zdv_bl") ) CALL iom_put( "zdv_bl", tmask(:,:,1)*zdv_bl ) ! dv at ml base 1189 IF ( iom_use("dh") ) CALL iom_put( "dh", tmask(:,:,1)*dh ) ! Initial boundary-layer depth 1190 IF ( iom_use("hml") ) CALL iom_put( "hml", tmask(:,:,1)*hml ) ! Initial boundary-layer depth 1345 1191 IF ( iom_use("dstokes") ) CALL iom_put( "dstokes", tmask(:,:,1)*dstokes ) ! Stokes drift penetration depth 1346 1192 IF ( iom_use("zustke") ) CALL iom_put( "zustke", tmask(:,:,1)*zustke ) ! Stokes drift magnitude at T-points … … 1348 1194 IF ( iom_use("zwstrl") ) CALL iom_put( "zwstrl", tmask(:,:,1)*zwstrl ) ! Langmuir velocity scale 1349 1195 IF ( iom_use("zustar") ) CALL iom_put( "zustar", tmask(:,:,1)*zustar ) ! friction velocity scale 1196 IF ( iom_use("zvstr") ) CALL iom_put( "zvstr", tmask(:,:,1)*zvstr ) ! mixed velocity scale 1197 IF ( iom_use("zla") ) CALL iom_put( "zla", tmask(:,:,1)*zla ) ! langmuir # 1350 1198 IF ( iom_use("wind_power") ) CALL iom_put( "wind_power", 1000.*rau0*tmask(:,:,1)*zustar**3 ) ! BL depth internal to zdf_osm routine 1351 1199 IF ( iom_use("wind_wave_power") ) CALL iom_put( "wind_wave_power", 1000.*rau0*tmask(:,:,1)*zustar**2*zustke ) 1352 1200 IF ( iom_use("zhbl") ) CALL iom_put( "zhbl", tmask(:,:,1)*zhbl ) ! BL depth internal to zdf_osm routine 1353 1201 IF ( iom_use("zhml") ) CALL iom_put( "zhml", tmask(:,:,1)*zhml ) ! ML depth internal to zdf_osm routine 1354 IF ( iom_use("zdh") ) CALL iom_put( "zdh", tmask(:,:,1)*zdh ) ! ML depth internal to zdf_osm routine 1202 IF ( iom_use("imld") ) CALL iom_put( "imld", tmask(:,:,1)*imld ) ! index for ML depth internal to zdf_osm routine 1203 IF ( iom_use("zdh") ) CALL iom_put( "zdh", tmask(:,:,1)*zdh ) ! pyc thicknessh internal to zdf_osm routine 1355 1204 IF ( iom_use("zhol") ) CALL iom_put( "zhol", tmask(:,:,1)*zhol ) ! ML depth internal to zdf_osm routine 1356 IF ( iom_use("zwthav") ) CALL iom_put( "zwthav", tmask(:,:,1)*zwthav ) ! ML depth internal to zdf_osm routine 1357 IF ( iom_use("zwth_ent") ) CALL iom_put( "zwth_ent", tmask(:,:,1)*zwth_ent ) ! ML depth internal to zdf_osm routine 1358 IF ( iom_use("zt_ml") ) CALL iom_put( "zt_ml", tmask(:,:,1)*zt_ml ) ! average T in ML 1205 IF ( iom_use("zwthav") ) CALL iom_put( "zwthav", tmask(:,:,1)*zwthav ) ! upward BL-avged turb temp flux 1206 IF ( iom_use("zwth_ent") ) CALL iom_put( "zwth_ent", tmask(:,:,1)*zwth_ent ) ! upward turb temp entrainment flux 1207 IF ( iom_use("zwb_ent") ) CALL iom_put( "zwb_ent", tmask(:,:,1)*zwb_ent ) ! upward turb buoyancy entrainment flux 1208 IF ( iom_use("zws_ent") ) CALL iom_put( "zws_ent", tmask(:,:,1)*zws_ent ) ! upward turb salinity entrainment flux 1209 IF ( iom_use("zt_ml") ) CALL iom_put( "zt_ml", tmask(:,:,1)*zt_ml ) ! average T in ML 1210 1211 IF ( iom_use("hmle") ) CALL iom_put( "hmle", tmask(:,:,1)*hmle ) ! FK layer depth 1212 IF ( iom_use("zmld") ) CALL iom_put( "zmld", tmask(:,:,1)*zmld ) ! FK target layer depth 1213 IF ( iom_use("zwb_fk") ) CALL iom_put( "zwb_fk", tmask(:,:,1)*zwb_fk ) ! FK b flux 1214 IF ( iom_use("zwb_fk_b") ) CALL iom_put( "zwb_fk_b", tmask(:,:,1)*zwb_fk_b ) ! FK b flux averaged over ML 1215 IF ( iom_use("mld_prof") ) CALL iom_put( "mld_prof", tmask(:,:,1)*mld_prof )! FK layer max k 1216 IF ( iom_use("zdtdx") ) CALL iom_put( "zdtdx", umask(:,:,1)*zdtdx ) ! FK dtdx at u-pt 1217 IF ( iom_use("zdtdy") ) CALL iom_put( "zdtdy", vmask(:,:,1)*zdtdy ) ! FK dtdy at v-pt 1218 IF ( iom_use("zdsdx") ) CALL iom_put( "zdsdx", umask(:,:,1)*zdsdx ) ! FK dtdx at u-pt 1219 IF ( iom_use("zdsdy") ) CALL iom_put( "zdsdy", vmask(:,:,1)*zdsdy ) ! FK dsdy at v-pt 1220 IF ( iom_use("dbdx_mle") ) CALL iom_put( "dbdx_mle", umask(:,:,1)*dbdx_mle ) ! FK dbdx at u-pt 1221 IF ( iom_use("dbdy_mle") ) CALL iom_put( "dbdy_mle", vmask(:,:,1)*dbdy_mle ) ! FK dbdy at v-pt 1222 IF ( iom_use("zdiff_mle") ) CALL iom_put( "zdiff_mle", tmask(:,:,1)*zdiff_mle )! FK diff in MLE at t-pt 1223 IF ( iom_use("zvel_mle") ) CALL iom_put( "zvel_mle", tmask(:,:,1)*zdiff_mle )! FK diff in MLE at t-pt 1224 1359 1225 END IF 1360 ! Lateral boundary conditions on p_avt (sign unchanged) 1361 CALL lbc_lnk( 'zdfosm', p_avt(:,:,:), 'W', 1. ) 1226 1227 CONTAINS 1228 1229 1230 ! Alan: do we need zb? 1231 SUBROUTINE zdf_osm_vertical_average( jnlev_av, zt, zs, zu, zv, zdt, zds, zdb, zdu, zdv ) 1232 !!--------------------------------------------------------------------- 1233 !! *** ROUTINE zdf_vertical_average *** 1234 !! 1235 !! ** Purpose : Determines vertical averages from surface to jnlev. 1236 !! 1237 !! ** Method : Averages are calculated from the surface to jnlev. 1238 !! The external level used to calculate differences is ibld+ibld_ext 1239 !! 1240 !!---------------------------------------------------------------------- 1241 1242 INTEGER, DIMENSION(jpi,jpj) :: jnlev_av ! Number of levels to average over. 1243 1244 ! Alan: do we need zb? 1245 REAL(wp), DIMENSION(jpi,jpj) :: zt, zs ! Average temperature and salinity 1246 REAL(wp), DIMENSION(jpi,jpj) :: zu,zv ! Average current components 1247 REAL(wp), DIMENSION(jpi,jpj) :: zdt, zds, zdb ! Difference between average and value at base of OSBL 1248 REAL(wp), DIMENSION(jpi,jpj) :: zdu, zdv ! Difference for velocity components. 1249 1250 INTEGER :: jk, ji, jj 1251 REAL(wp) :: zthick, zthermal, zbeta 1252 1253 1254 zt = 0._wp 1255 zs = 0._wp 1256 zu = 0._wp 1257 zv = 0._wp 1258 DO jj = 2, jpjm1 ! Vertical slab 1259 DO ji = 2, jpim1 1260 zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 1261 zbeta = rab_n(ji,jj,1,jp_sal) 1262 ! average over depth of boundary layer 1263 zthick = epsln 1264 DO jk = 2, jnlev_av(ji,jj) 1265 zthick = zthick + e3t_n(ji,jj,jk) 1266 zt(ji,jj) = zt(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem) 1267 zs(ji,jj) = zs(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) 1268 zu(ji,jj) = zu(ji,jj) + e3t_n(ji,jj,jk) & 1269 & * ( ub(ji,jj,jk) + ub(ji - 1,jj,jk) ) & 1270 & / MAX( 1. , umask(ji,jj,jk) + umask(ji - 1,jj,jk) ) 1271 zv(ji,jj) = zv(ji,jj) + e3t_n(ji,jj,jk) & 1272 & * ( vb(ji,jj,jk) + vb(ji,jj - 1,jk) ) & 1273 & / MAX( 1. , vmask(ji,jj,jk) + vmask(ji,jj - 1,jk) ) 1274 END DO 1275 zt(ji,jj) = zt(ji,jj) / zthick 1276 zs(ji,jj) = zs(ji,jj) / zthick 1277 zu(ji,jj) = zu(ji,jj) / zthick 1278 zv(ji,jj) = zv(ji,jj) / zthick 1279 ! Alan: do we need zb? 1280 zdt(ji,jj) = zt(ji,jj) - tsn(ji,jj,ibld(ji,jj)+ibld_ext,jp_tem) 1281 zds(ji,jj) = zs(ji,jj) - tsn(ji,jj,ibld(ji,jj)+ibld_ext,jp_sal) 1282 zdu(ji,jj) = zu(ji,jj) - ( ub(ji,jj,ibld(ji,jj)+ibld_ext) + ub(ji-1,jj,ibld(ji,jj)+ibld_ext ) ) & 1283 & / MAX(1. , umask(ji,jj,ibld(ji,jj)+ibld_ext ) + umask(ji-1,jj,ibld(ji,jj)+ibld_ext ) ) 1284 zdv(ji,jj) = zv(ji,jj) - ( vb(ji,jj,ibld(ji,jj)+ibld_ext) + vb(ji,jj-1,ibld(ji,jj)+ibld_ext ) ) & 1285 & / MAX(1. , vmask(ji,jj,ibld(ji,jj)+ibld_ext ) + vmask(ji,jj-1,ibld(ji,jj)+ibld_ext ) ) 1286 zdb(ji,jj) = grav * zthermal * zdt(ji,jj) - grav * zbeta * zds(ji,jj) 1287 END DO 1288 END DO 1289 END SUBROUTINE zdf_osm_vertical_average 1290 1291 SUBROUTINE zdf_osm_velocity_rotation( zcos_w, zsin_w, zu, zv, zdu, zdv ) 1292 !!--------------------------------------------------------------------- 1293 !! *** ROUTINE zdf_velocity_rotation *** 1294 !! 1295 !! ** Purpose : Rotates frame of reference of averaged velocity components. 1296 !! 1297 !! ** Method : The velocity components are rotated into frame specified by zcos_w and zsin_w 1298 !! 1299 !!---------------------------------------------------------------------- 1300 1301 REAL(wp), DIMENSION(jpi,jpj) :: zcos_w, zsin_w ! Cos and Sin of rotation angle 1302 REAL(wp), DIMENSION(jpi,jpj) :: zu, zv ! Components of current 1303 REAL(wp), DIMENSION(jpi,jpj) :: zdu, zdv ! Change in velocity components across pycnocline 1304 1305 INTEGER :: ji, jj 1306 REAL(wp) :: ztemp 1307 1308 DO jj = 2, jpjm1 1309 DO ji = 2, jpim1 1310 ztemp = zu(ji,jj) 1311 zu(ji,jj) = zu(ji,jj) * zcos_w(ji,jj) + zv(ji,jj) * zsin_w(ji,jj) 1312 zv(ji,jj) = zv(ji,jj) * zcos_w(ji,jj) - ztemp * zsin_w(ji,jj) 1313 ztemp = zdu(ji,jj) 1314 zdu(ji,jj) = zdu(ji,jj) * zcos_w(ji,jj) + zdv(ji,jj) * zsin_w(ji,jj) 1315 zdv(ji,jj) = zdv(ji,jj) * zsin_w(ji,jj) - ztemp * zsin_w(ji,jj) 1316 END DO 1317 END DO 1318 END SUBROUTINE zdf_osm_velocity_rotation 1319 1320 SUBROUTINE zdf_osm_external_gradients( zdtdz, zdsdz, zdbdz ) 1321 !!--------------------------------------------------------------------- 1322 !! *** ROUTINE zdf_osm_external_gradients *** 1323 !! 1324 !! ** Purpose : Calculates the gradients below the OSBL 1325 !! 1326 !! ** Method : Uses ibld and ibld_ext to determine levels to calculate the gradient. 1327 !! 1328 !!---------------------------------------------------------------------- 1329 1330 REAL(wp), DIMENSION(jpi,jpj) :: zdtdz, zdsdz, zdbdz ! External gradients of temperature, salinity and buoyancy. 1331 1332 INTEGER :: jj, ji, jkb, jkb1 1333 REAL(wp) :: zthermal, zbeta 1334 1335 1336 DO jj = 2, jpjm1 1337 DO ji = 2, jpim1 1338 IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN 1339 zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 1340 zbeta = rab_n(ji,jj,1,jp_sal) 1341 jkb = ibld(ji,jj) + ibld_ext 1342 jkb1 = MIN(jkb + 1, mbkt(ji,jj)) 1343 zdtdz(ji,jj) = - ( tsn(ji,jj,jkb1,jp_tem) - tsn(ji,jj,jkb,jp_tem ) ) & 1344 & / e3t_n(ji,jj,ibld(ji,jj)) 1345 zdsdz(ji,jj) = - ( tsn(ji,jj,jkb1,jp_sal) - tsn(ji,jj,jkb,jp_sal ) ) & 1346 & / e3t_n(ji,jj,ibld(ji,jj)) 1347 zdbdz(ji,jj) = grav * zthermal * zdtdz(ji,jj) - grav * zbeta * zdsdz(ji,jj) 1348 END IF 1349 END DO 1350 END DO 1351 END SUBROUTINE zdf_osm_external_gradients 1352 1353 SUBROUTINE zdf_osm_pycnocline_scalar_profiles( zdtdz, zdsdz, zdbdz ) 1354 1355 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdtdz, zdsdz, zdbdz ! gradients in the pycnocline 1356 1357 INTEGER :: jk, jj, ji 1358 REAL(wp) :: ztgrad, zsgrad, zbgrad 1359 REAL(wp) :: zgamma_b_nd, zgamma_c, znd 1360 REAL(wp) :: zzeta_s=0.3, ztmp 1361 1362 DO jj = 2, jpjm1 1363 DO ji = 2, jpim1 1364 IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN 1365 IF ( lconv(ji,jj) ) THEN ! convective conditions 1366 IF ( zdbdz_ext(ji,jj) > 0._wp .AND. & 1367 & (zdhdt(ji,jj) > 0._wp .AND. ln_osm_mle .AND. zdb_bl(ji,jj) > rn_osm_mle_thresh & 1368 & .OR. zdb_bl(ji,jj) > 0._wp)) THEN ! zdhdt could be <0 due to FK, hence check zdhdt>0 1369 ztmp = 1._wp/MAX(zdh(ji,jj), epsln) 1370 ztgrad = 0.5 * zdt_ml(ji,jj) * ztmp + zdtdz_ext(ji,jj) 1371 zsgrad = 0.5 * zds_ml(ji,jj) * ztmp + zdsdz_ext(ji,jj) 1372 zbgrad = 0.5 * zdb_ml(ji,jj) * ztmp + zdbdz_ext(ji,jj) 1373 zgamma_b_nd = zdbdz_ext(ji,jj) * zdh(ji,jj) / MAX(zdb_ml(ji,jj), epsln) 1374 zgamma_c = ( 3.14159 / 4.0 ) * ( 0.5 + zgamma_b_nd ) /& 1375 & ( 1.0 - 0.25 * SQRT( 3.14159 / 6.0 ) - 2.0 * zgamma_b_nd * zzeta_s )**2 ! check 1376 DO jk = 2, ibld(ji,jj)+ibld_ext 1377 znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) * ztmp 1378 IF ( znd <= zzeta_s ) THEN 1379 zdtdz(ji,jj,jk) = zdtdz_ext(ji,jj) + 0.5 * zdt_ml(ji,jj) * ztmp * & 1380 & EXP( -6.0 * ( znd -zzeta_s )**2 ) 1381 zdsdz(ji,jj,jk) = zdsdz_ext(ji,jj) + 0.5 * zds_ml(ji,jj) * ztmp * & 1382 & EXP( -6.0 * ( znd -zzeta_s )**2 ) 1383 zdbdz(ji,jj,jk) = zdbdz_ext(ji,jj) + 0.5 * zdb_ml(ji,jj) * ztmp * & 1384 & EXP( -6.0 * ( znd -zzeta_s )**2 ) 1385 ELSE 1386 zdtdz(ji,jj,jk) = ztgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 ) 1387 zdsdz(ji,jj,jk) = zsgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 ) 1388 zdbdz(ji,jj,jk) = zbgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 ) 1389 ENDIF 1390 END DO 1391 ENDIF ! If condition not satisfied, no pycnocline present. Gradients have been initialised to zero, so do nothing 1392 ELSE 1393 ! stable conditions 1394 ! if pycnocline profile only defined when depth steady of increasing. 1395 IF ( zdhdt(ji,jj) >= 0.0 ) THEN ! Depth increasing, or steady. 1396 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 1397 IF ( zhol(ji,jj) >= 0.5 ) THEN ! Very stable - 'thick' pycnocline 1398 ztmp = 1._wp/MAX(zhbl(ji,jj), epsln) 1399 ztgrad = zdt_bl(ji,jj) * ztmp 1400 zsgrad = zds_bl(ji,jj) * ztmp 1401 zbgrad = zdb_bl(ji,jj) * ztmp 1402 DO jk = 2, ibld(ji,jj) 1403 znd = gdepw_n(ji,jj,jk) * ztmp 1404 zdtdz(ji,jj,jk) = ztgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 1405 zdbdz(ji,jj,jk) = zbgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 1406 zdsdz(ji,jj,jk) = zsgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 1407 END DO 1408 ELSE ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form. 1409 ztmp = 1._wp/MAX(zdh(ji,jj), epsln) 1410 ztgrad = zdt_bl(ji,jj) * ztmp 1411 zsgrad = zds_bl(ji,jj) * ztmp 1412 zbgrad = zdb_bl(ji,jj) * ztmp 1413 DO jk = 2, ibld(ji,jj) 1414 znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) * ztmp 1415 zdtdz(ji,jj,jk) = ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 1416 zdbdz(ji,jj,jk) = zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 1417 zdsdz(ji,jj,jk) = zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 1418 END DO 1419 ENDIF ! IF (zhol >=0.5) 1420 ENDIF ! IF (zdb_bl> 0.) 1421 ENDIF ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero 1422 ENDIF ! IF (lconv) 1423 END IF ! IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) 1424 END DO 1425 END DO 1426 1427 END SUBROUTINE zdf_osm_pycnocline_scalar_profiles 1428 1429 SUBROUTINE zdf_osm_pycnocline_shear_profiles( zdudz, zdvdz ) 1430 !!--------------------------------------------------------------------- 1431 !! *** ROUTINE zdf_osm_pycnocline_shear_profiles *** 1432 !! 1433 !! ** Purpose : Calculates velocity shear in the pycnocline 1434 !! 1435 !! ** Method : 1436 !! 1437 !!---------------------------------------------------------------------- 1438 1439 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdudz, zdvdz 1440 1441 INTEGER :: jk, jj, ji 1442 REAL(wp) :: zugrad, zvgrad, znd 1443 REAL(wp) :: zzeta_v = 0.45 1362 1444 ! 1363 END SUBROUTINE zdf_osm 1445 DO jj = 2, jpjm1 1446 DO ji = 2, jpim1 1447 ! 1448 IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN 1449 IF ( lconv (ji,jj) ) THEN 1450 ! Unstable conditions 1451 zugrad = 0.7 * zdu_ml(ji,jj) / zdh(ji,jj) + 0.3 * zustar(ji,jj)*zustar(ji,jj) / & 1452 & ( ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) * & 1453 & MIN(zla(ji,jj)**(8.0/3.0) + epsln, 0.12 )) 1454 !Alan is this right? 1455 zvgrad = ( 0.7 * zdv_ml(ji,jj) + & 1456 & 2.0 * ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) / & 1457 & ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird + epsln ) & 1458 & )/ (zdh(ji,jj) + epsln ) 1459 DO jk = 2, ibld(ji,jj) - 1 + ibld_ext 1460 znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / (zdh(ji,jj) + epsln ) - zzeta_v 1461 IF ( znd <= 0.0 ) THEN 1462 zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( 3.0 * znd ) 1463 zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( 3.0 * znd ) 1464 ELSE 1465 zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( -2.0 * znd ) 1466 zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( -2.0 * znd ) 1467 ENDIF 1468 END DO 1469 ELSE 1470 ! stable conditions 1471 zugrad = 3.25 * zdu_bl(ji,jj) / zhbl(ji,jj) 1472 zvgrad = 2.75 * zdv_bl(ji,jj) / zhbl(ji,jj) 1473 DO jk = 2, ibld(ji,jj) 1474 znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 1475 IF ( znd < 1.0 ) THEN 1476 zdudz(ji,jj,jk) = zugrad * EXP( -40.0 * ( znd - 1.0 )**2 ) 1477 ELSE 1478 zdudz(ji,jj,jk) = zugrad * EXP( -20.0 * ( znd - 1.0 )**2 ) 1479 ENDIF 1480 zdvdz(ji,jj,jk) = zvgrad * EXP( -20.0 * ( znd - 0.85 )**2 ) 1481 END DO 1482 ENDIF 1483 ! 1484 END IF ! IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) 1485 END DO 1486 END DO 1487 END SUBROUTINE zdf_osm_pycnocline_shear_profiles 1488 1489 SUBROUTINE zdf_osm_calculate_dhdt( zdhdt, zdhdt_2 ) 1490 !!--------------------------------------------------------------------- 1491 !! *** ROUTINE zdf_osm_calculate_dhdt *** 1492 !! 1493 !! ** Purpose : Calculates the rate at which hbl changes. 1494 !! 1495 !! ** Method : 1496 !! 1497 !!---------------------------------------------------------------------- 1498 1499 REAL(wp), DIMENSION(jpi,jpj) :: zdhdt, zdhdt_2 ! Rate of change of hbl 1500 1501 INTEGER :: jj, ji 1502 REAL(wp) :: zgamma_b_nd, zgamma_dh_nd, zpert 1503 REAL(wp) :: zvel_max, zwb_min 1504 REAL(wp) :: zwcor, zrf_conv, zrf_shear, zrf_langmuir, zr_stokes 1505 REAL(wp) :: zzeta_m = 0.3 1506 REAL(wp) :: zgamma_c = 2.0 1507 REAL(wp) :: zdhoh = 0.1 1508 REAL(wp) :: alpha_bc = 0.5 1509 1510 DO jj = 2, jpjm1 1511 DO ji = 2, jpim1 1512 IF ( lconv(ji,jj) ) THEN ! Convective 1513 ! Alan is this right? Yes, it's a bit different from the previous relationship 1514 ! zwb_ent(ji,jj) = - 2.0 * 0.2 * zwbav(ji,jj) & 1515 ! & - ( 0.15 * ( 1.0 - EXP( -1.5 * zla(ji,jj) ) ) * zustar(ji,jj)**3 + 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj) 1516 zwcor = ABS(ff_t(ji,jj)) * zhbl(ji,jj) + epsln 1517 zrf_conv = TANH( ( zwstrc(ji,jj) / zwcor )**0.69 ) 1518 zrf_shear = TANH( ( zustar(ji,jj) / zwcor )**0.69 ) 1519 zrf_langmuir = TANH( ( zwstrl(ji,jj) / zwcor )**0.69 ) 1520 zr_stokes = 1.0 - EXP( -25.0 * dstokes(ji,jj) / hbl(ji,jj) & 1521 & * ( 1.0 + 4.0 * dstokes(ji,jj) / hbl(ji,jj) ) ) 1522 1523 zwb_ent(ji,jj) = - 2.0 * 0.2 * zrf_conv * zwbav(ji,jj) & 1524 & - 0.15 * zrf_shear * zustar(ji,jj)**3 /zhml(ji,jj) & 1525 & + zr_stokes * ( 0.15 * EXP( -1.5 * zla(ji,jj) ) * zrf_shear * zustar(ji,jj)**3 & 1526 & - zrf_langmuir * 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj) 1527 ! 1528 zwb_min = dh(ji,jj) * zwb0(ji,jj) / zhml(ji,jj) + zwb_ent(ji,jj) 1529 1530 IF ( ln_osm_mle ) THEN 1531 ! zwb_fk_b(ji,jj) = zwb_fk(ji,jj) * hmle(ji,jj) / ( 6.0 * hbl(ji,jj) ) * ( 6.0 * hbl(ji,jj) / hmle(ji,jj) - 1.0 + & 1532 ! & ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3 ) ! Fox-Kemper buoyancy flux average over OSBL 1533 IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN 1534 zwb_fk_b(ji,jj) = zwb_fk(ji,jj) * & 1535 (1.0 + hmle(ji,jj) / ( 6.0 * hbl(ji,jj) ) * (-1.0 + ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3) ) 1536 ELSE 1537 zwb_fk_b(ji,jj) = 0.5 * zwb_fk(ji,jj) * hmle(ji,jj) / hbl(ji,jj) 1538 ENDIF 1539 zvel_max = ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**p2third / hbl(ji,jj) 1540 IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0.0 ) THEN 1541 ! OSBL is deepening, entrainment > restratification 1542 IF ( zdb_bl(ji,jj) > 0.0 .and. zdbdz_ext(ji,jj) > 0.0 ) THEN 1543 zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) 1544 ELSE 1545 zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / MAX( zvel_max, 1.0e-15) 1546 ENDIF 1547 ELSE 1548 ! OSBL shoaling due to restratification flux. This is the velocity defined in Fox-Kemper et al (2008) 1549 zdhdt(ji,jj) = - zvel_mle(ji,jj) 1550 1551 1552 ENDIF 1553 1554 ELSE 1555 ! Fox-Kemper not used. 1556 1557 zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / hbl(ji,jj) ) * zwb_ent(ji,jj) / & 1558 & ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 1559 zdhdt(ji,jj) = -zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) 1560 ! added ajgn 23 July as temporay fix 1561 1562 ENDIF 1563 1564 zdhdt_2(ji,jj) = 0._wp 1565 1566 ! commented out ajgn 23 July as temporay fix 1567 ! IF ( zdb_ml(ji,jj) > 0.0 .and. zdbdz_ext(ji,jj) > 0.0 ) THEN 1568 ! !additional term to account for thickness of pycnocline on dhdt. Small correction, so could get rid of this if necessary. 1569 ! zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 1570 ! zgamma_b_nd = zdbdz_ext(ji,jj) * zhml(ji,jj) / zdb_ml(ji,jj) 1571 ! zgamma_dh_nd = zdbdz_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj) 1572 ! zdhdt_2(ji,jj) = ( 1.0 - SQRT( 3.1415 / ( 4.0 * zgamma_c) ) * zdhoh ) * zdh(ji,jj) / zhml(ji,jj) 1573 ! zdhdt_2(ji,jj) = zdhdt_2(ji,jj) * ( zwb0(ji,jj) - (1.0 + zgamma_b_nd / alpha_bc ) * zwb_min ) 1574 ! ! Alan no idea what this should be? 1575 ! zdhdt_2(ji,jj) = alpha_bc / ( 4.0 * zgamma_c ) * zdhdt_2(ji,jj) & 1576 ! & + (alpha_bc + zgamma_dh_nd ) * ( 1.0 + SQRT( 3.1414 / ( 4.0 * zgamma_c ) ) * zdh(ji,jj) / zhbl(ji,jj) ) & 1577 ! & * (1.0 / ( 4.0 * zgamma_c * alpha_bc ) ) * zwb_min * zdh(ji,jj) / zhbl(ji,jj) 1578 ! zdhdt_2(ji,jj) = zdhdt_2(ji,jj) / ( zvel_max + MAX( zdb_bl(ji,jj), 1.0e-15 ) ) 1579 ! IF ( zdhdt_2(ji,jj) <= 0.2 * zdhdt(ji,jj) ) THEN 1580 ! zdhdt(ji,jj) = zdhdt(ji,jj) + zdhdt_2(ji,jj) 1581 ! ENDIF 1582 ELSE ! Stable 1583 zdhdt(ji,jj) = ( 0.06 + 0.52 * zhol(ji,jj) / 2.0 ) * zvstr(ji,jj)**3 / hbl(ji,jj) + zwbav(ji,jj) 1584 zdhdt_2(ji,jj) = 0._wp 1585 IF ( zdhdt(ji,jj) < 0._wp ) THEN 1586 ! For long timsteps factor in brackets slows the rapid collapse of the OSBL 1587 zpert = 2.0 * ( 1.0 + 0.0 * 2.0 * zvstr(ji,jj) * rn_rdt / hbl(ji,jj) ) * zvstr(ji,jj)**2 / hbl(ji,jj) 1588 ELSE 1589 zpert = MAX( 2.0 * ( 1.0 + 0.0 * 2.0 * zvstr(ji,jj) * rn_rdt / hbl(ji,jj) ) * zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) ) 1590 ENDIF 1591 zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / zpert 1592 ENDIF 1593 END DO 1594 END DO 1595 END SUBROUTINE zdf_osm_calculate_dhdt 1596 1597 SUBROUTINE zdf_osm_timestep_hbl( zdhdt, zdhdt_2 ) 1598 !!--------------------------------------------------------------------- 1599 !! *** ROUTINE zdf_osm_timestep_hbl *** 1600 !! 1601 !! ** Purpose : Increments hbl. 1602 !! 1603 !! ** Method : If thechange in hbl exceeds one model level the change is 1604 !! is calculated by moving down the grid, changing the buoyancy 1605 !! jump. This is to ensure that the change in hbl does not 1606 !! overshoot a stable layer. 1607 !! 1608 !!---------------------------------------------------------------------- 1609 1610 1611 REAL(wp), DIMENSION(jpi,jpj) :: zdhdt, zdhdt_2 ! rates of change of hbl. 1612 1613 INTEGER :: jk, jj, ji, jm 1614 REAL(wp) :: zhbl_s, zvel_max, zdb 1615 REAL(wp) :: zthermal, zbeta 1616 1617 DO jj = 2, jpjm1 1618 DO ji = 2, jpim1 1619 IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN 1620 ! 1621 ! If boundary layer changes by more than one level, need to check for stable layers between initial and final depths. 1622 ! 1623 zhbl_s = hbl(ji,jj) 1624 jm = imld(ji,jj) 1625 zthermal = rab_n(ji,jj,1,jp_tem) 1626 zbeta = rab_n(ji,jj,1,jp_sal) 1627 1628 1629 IF ( lconv(ji,jj) ) THEN 1630 !unstable 1631 1632 IF( ln_osm_mle ) THEN 1633 zvel_max = ( zwstrl(ji,jj)**3 + zwstrc(ji,jj)**3 )**p2third / hbl(ji,jj) 1634 ELSE 1635 1636 zvel_max = -( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / hbl(ji,jj) ) * zwb_ent(ji,jj) / & 1637 & ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 1638 1639 ENDIF 1640 1641 DO jk = imld(ji,jj), ibld(ji,jj) 1642 zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) ) & 1643 & - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ), & 1644 & 0.0 ) + zvel_max 1645 1646 1647 IF ( ln_osm_mle ) THEN 1648 zhbl_s = zhbl_s + MIN( & 1649 & rn_rdt * ( ( -zwb_ent(ji,jj) - 2.0 * zwb_fk_b(ji,jj) )/ zdb + zdhdt_2(ji,jj) ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), & 1650 & e3w_n(ji,jj,jm) ) 1651 ELSE 1652 zhbl_s = zhbl_s + MIN( & 1653 & rn_rdt * ( -zwb_ent(ji,jj) / zdb + zdhdt_2(ji,jj) ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), & 1654 & e3w_n(ji,jj,jm) ) 1655 ENDIF 1656 1657 zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol) 1658 1659 IF ( zhbl_s >= gdepw_n(ji,jj,jm+1) ) jm = jm + 1 1660 END DO 1661 hbl(ji,jj) = zhbl_s 1662 ibld(ji,jj) = jm 1663 ELSE 1664 ! stable 1665 DO jk = imld(ji,jj), ibld(ji,jj) 1666 zdb = MAX( & 1667 & grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) )& 1668 & - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ),& 1669 & 0.0 ) + & 1670 & 2.0 * zvstr(ji,jj)**2 / zhbl_s 1671 1672 ! Alan is thuis right? I have simply changed hbli to hbl 1673 zhol(ji,jj) = -zhbl_s / ( ( zvstr(ji,jj)**3 + epsln )/ zwbav(ji,jj) ) 1674 zdhdt(ji,jj) = -( zwbav(ji,jj) - 0.04 / 2.0 * zwstrl(ji,jj)**3 / zhbl_s - 0.15 / 2.0 * ( 1.0 - EXP( -1.5 * zla(ji,jj) ) ) * & 1675 & zustar(ji,jj)**3 / zhbl_s ) * ( 0.725 + 0.225 * EXP( -7.5 * zhol(ji,jj) ) ) 1676 zdhdt(ji,jj) = zdhdt(ji,jj) + zwbav(ji,jj) 1677 zhbl_s = zhbl_s + MIN( zdhdt(ji,jj) / zdb * rn_rdt / FLOAT( ibld(ji,jj) - imld(ji,jj) ), e3w_n(ji,jj,jm) ) 1678 1679 zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol) 1680 IF ( zhbl_s >= gdepw_n(ji,jj,jm) ) jm = jm + 1 1681 END DO 1682 ENDIF ! IF ( lconv ) 1683 hbl(ji,jj) = MAX(zhbl_s, gdepw_n(ji,jj,4) ) 1684 ibld(ji,jj) = MAX(jm, 4 ) 1685 ELSE 1686 ! change zero or one model level. 1687 hbl(ji,jj) = MAX(zhbl_t(ji,jj), gdepw_n(ji,jj,4) ) 1688 ENDIF 1689 zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj)) 1690 END DO 1691 END DO 1692 1693 END SUBROUTINE zdf_osm_timestep_hbl 1694 1695 SUBROUTINE zdf_osm_pycnocline_thickness( dh, zdh ) 1696 !!--------------------------------------------------------------------- 1697 !! *** ROUTINE zdf_osm_pycnocline_thickness *** 1698 !! 1699 !! ** Purpose : Calculates thickness of the pycnocline 1700 !! 1701 !! ** Method : The thickness is calculated from a prognostic equation 1702 !! that relaxes the pycnocine thickness to a diagnostic 1703 !! value. The time change is calculated assuming the 1704 !! thickness relaxes exponentially. This is done to deal 1705 !! with large timesteps. 1706 !! 1707 !!---------------------------------------------------------------------- 1708 1709 REAL(wp), DIMENSION(jpi,jpj) :: dh, zdh ! pycnocline thickness. 1710 ! 1711 INTEGER :: jj, ji 1712 INTEGER :: inhml 1713 REAL(wp) :: zari, ztau, zddhdt 1714 1715 1716 DO jj = 2, jpjm1 1717 DO ji = 2, jpim1 1718 IF ( lconv(ji,jj) ) THEN 1719 1720 IF( ln_osm_mle ) THEN 1721 IF ( ( zwb_ent(ji,jj) + zwb_fk_b(ji,jj) ) < 0._wp ) THEN 1722 ! OSBL is deepening. Note wb_fk_b is zero if ln_osm_mle=F 1723 IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_ext(ji,jj) > 0._wp)THEN 1724 IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN ! near neutral stability 1725 zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / & 1726 (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zvstr(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 ) 1727 ELSE ! unstable 1728 zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / & 1729 (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zwstrc(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 ) 1730 ENDIF 1731 ztau = 0.2 * hbl(ji,jj) / (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird 1732 zddhdt = zari * hbl(ji,jj) 1733 ELSE 1734 ztau = 0.2 * hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 1735 zddhdt = 0.2 * hbl(ji,jj) 1736 ENDIF 1737 ELSE 1738 ztau = 0.2 * hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 1739 zddhdt = 0.2 * hbl(ji,jj) 1740 ENDIF 1741 ELSE ! ln_osm_mle 1742 IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_ext(ji,jj) > 0._wp)THEN 1743 IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN ! near neutral stability 1744 zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / & 1745 (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zvstr(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 ) 1746 ELSE ! unstable 1747 zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / & 1748 (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zwstrc(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 ) 1749 ENDIF 1750 ztau = hbl(ji,jj) / (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird 1751 zddhdt = zari * hbl(ji,jj) 1752 ELSE 1753 ztau = hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 1754 zddhdt = 0.2 * hbl(ji,jj) 1755 ENDIF 1756 1757 END IF ! ln_osm_mle 1758 1759 dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau ) + zddhdt * ( 1.0 - EXP( -rn_rdt / ztau ) ) 1760 ! Alan: this hml is never defined or used 1761 ELSE ! IF (lconv) 1762 ztau = hbl(ji,jj) / zvstr(ji,jj) 1763 IF ( zdhdt(ji,jj) >= 0.0 ) THEN ! probably shouldn't include wm here 1764 ! boundary layer deepening 1765 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 1766 ! pycnocline thickness set by stratification - use same relationship as for neutral conditions. 1767 zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) & 1768 & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01 , 0.2 ) 1769 zddhdt = MIN( zari, 0.2 ) * hbl(ji,jj) 1770 ELSE 1771 zddhdt = 0.2 * hbl(ji,jj) 1772 ENDIF 1773 ELSE ! IF(dhdt < 0) 1774 zddhdt = 0.2 * hbl(ji,jj) 1775 ENDIF ! IF (dhdt >= 0) 1776 dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau )+ zddhdt * ( 1.0 - EXP( -rn_rdt / ztau ) ) 1777 IF ( zdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zddhdt ! can be a problem with dh>hbl for rapid collapse 1778 ! Alan: this hml is never defined or used -- do we need it? 1779 ENDIF ! IF (lconv) 1780 1781 hml(ji,jj) = hbl(ji,jj) - dh(ji,jj) 1782 inhml = MAX( INT( dh(ji,jj) / e3t_n(ji,jj,ibld(ji,jj)) ) , 1 ) 1783 imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 3) 1784 zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 1785 zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 1786 END DO 1787 END DO 1788 1789 END SUBROUTINE zdf_osm_pycnocline_thickness 1790 1791 1792 SUBROUTINE zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle ) 1793 !!---------------------------------------------------------------------- 1794 !! *** ROUTINE zdf_osm_horizontal_gradients *** 1795 !! 1796 !! ** Purpose : Calculates horizontal gradients of buoyancy for use with Fox-Kemper parametrization. 1797 !! 1798 !! ** Method : 1799 !! 1800 !! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008 1801 !! Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008 1802 1803 1804 REAL(wp), DIMENSION(jpi,jpj) :: dbdx_mle, dbdy_mle ! MLE horiz gradients at u & v points 1805 REAL(wp), DIMENSION(jpi,jpj) :: zmld ! == estimated FK BLD used for MLE horiz gradients == ! 1806 REAL(wp), DIMENSION(jpi,jpj) :: zdtdx, zdtdy, zdsdx, zdsdy 1807 1808 INTEGER :: ji, jj, jk ! dummy loop indices 1809 INTEGER :: ii, ij, ik, ikmax ! local integers 1810 REAL(wp) :: zc 1811 REAL(wp) :: zN2_c ! local buoyancy difference from 10m value 1812 REAL(wp), DIMENSION(jpi,jpj) :: ztm, zsm, zLf_NH, zLf_MH 1813 REAL(wp), DIMENSION(jpi,jpj,jpts):: ztsm_midu, ztsm_midv, zabu, zabv 1814 REAL(wp), DIMENSION(jpi,jpj) :: zmld_midu, zmld_midv 1815 !!---------------------------------------------------------------------- 1816 ! 1817 ! !== MLD used for MLE ==! 1818 1819 mld_prof(:,:) = nlb10 ! Initialization to the number of w ocean point 1820 zmld(:,:) = 0._wp ! here hmlp used as a dummy variable, integrating vertically N^2 1821 zN2_c = grav * rn_osm_mle_rho_c * r1_rau0 ! convert density criteria into N^2 criteria 1822 DO jk = nlb10, jpkm1 1823 DO jj = 1, jpj ! Mixed layer level: w-level 1824 DO ji = 1, jpi 1825 ikt = mbkt(ji,jj) 1826 zmld(ji,jj) = zmld(ji,jj) + MAX( rn2b(ji,jj,jk) , 0._wp ) * e3w_n(ji,jj,jk) 1827 IF( zmld(ji,jj) < zN2_c ) mld_prof(ji,jj) = MIN( jk , ikt ) + 1 ! Mixed layer level 1828 END DO 1829 END DO 1830 END DO 1831 DO jj = 1, jpj 1832 DO ji = 1, jpi 1833 mld_prof(ji,jj) = MAX(mld_prof(ji,jj),ibld(ji,jj)) 1834 zmld(ji,jj) = gdepw_n(ji,jj,mld_prof(ji,jj)) 1835 END DO 1836 END DO 1837 ! ensure mld_prof .ge. ibld 1838 ! 1839 ikmax = MIN( MAXVAL( mld_prof(:,:) ), jpkm1 ) ! max level of the computation 1840 ! 1841 ztm(:,:) = 0._wp 1842 zsm(:,:) = 0._wp 1843 DO jk = 1, ikmax ! MLD and mean buoyancy and N2 over the mixed layer 1844 DO jj = 1, jpj 1845 DO ji = 1, jpi 1846 zc = e3t_n(ji,jj,jk) * REAL( MIN( MAX( 0, mld_prof(ji,jj)-jk ) , 1 ) ) ! zc being 0 outside the ML t-points 1847 ztm(ji,jj) = ztm(ji,jj) + zc * tsn(ji,jj,jk,jp_tem) 1848 zsm(ji,jj) = zsm(ji,jj) + zc * tsn(ji,jj,jk,jp_sal) 1849 END DO 1850 END DO 1851 END DO 1852 ! average temperature and salinity. 1853 ztm(:,:) = ztm(:,:) / MAX( e3t_n(:,:,1), zmld(:,:) ) 1854 zsm(:,:) = zsm(:,:) / MAX( e3t_n(:,:,1), zmld(:,:) ) 1855 ! calculate horizontal gradients at u & v points 1856 1857 DO jj = 2, jpjm1 1858 DO ji = 1, jpim1 1859 zdtdx(ji,jj) = ( ztm(ji+1,jj) - ztm( ji,jj) ) * umask(ji,jj,1) / e1u(ji,jj) 1860 zdsdx(ji,jj) = ( zsm(ji+1,jj) - zsm( ji,jj) ) * umask(ji,jj,1) / e1u(ji,jj) 1861 zmld_midu(ji,jj) = 0.25_wp * (zmld(ji+1,jj) + zmld( ji,jj)) 1862 ztsm_midu(ji,jj,jp_tem) = 0.5_wp * ( ztm(ji+1,jj) + ztm( ji,jj) ) 1863 ztsm_midu(ji,jj,jp_sal) = 0.5_wp * ( zsm(ji+1,jj) + zsm( ji,jj) ) 1864 END DO 1865 END DO 1866 1867 DO jj = 1, jpjm1 1868 DO ji = 2, jpim1 1869 zdtdy(ji,jj) = ( ztm(ji,jj+1) - ztm( ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj) 1870 zdsdy(ji,jj) = ( zsm(ji,jj+1) - zsm( ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj) 1871 zmld_midv(ji,jj) = 0.25_wp * (zmld(ji,jj+1) + zmld( ji,jj)) 1872 ztsm_midv(ji,jj,jp_tem) = 0.5_wp * ( ztm(ji,jj+1) + ztm( ji,jj) ) 1873 ztsm_midv(ji,jj,jp_sal) = 0.5_wp * ( zsm(ji,jj+1) + zsm( ji,jj) ) 1874 END DO 1875 END DO 1876 1877 CALL eos_rab(ztsm_midu, zmld_midu, zabu) 1878 CALL eos_rab(ztsm_midv, zmld_midv, zabv) 1879 1880 DO jj = 2, jpjm1 1881 DO ji = 1, jpim1 1882 dbdx_mle(ji,jj) = grav*(zdtdx(ji,jj)*zabu(ji,jj,jp_tem) - zdsdx(ji,jj)*zabu(ji,jj,jp_sal)) 1883 END DO 1884 END DO 1885 DO jj = 1, jpjm1 1886 DO ji = 2, jpim1 1887 dbdy_mle(ji,jj) = grav*(zdtdy(ji,jj)*zabv(ji,jj,jp_tem) - zdsdy(ji,jj)*zabv(ji,jj,jp_sal)) 1888 END DO 1889 END DO 1890 1891 END SUBROUTINE zdf_osm_zmld_horizontal_gradients 1892 SUBROUTINE zdf_osm_mle_parameters( hmle, zwb_fk, zvel_mle, zdiff_mle ) 1893 !!---------------------------------------------------------------------- 1894 !! *** ROUTINE zdf_osm_mle_parameters *** 1895 !! 1896 !! ** Purpose : Timesteps the mixed layer eddy depth, hmle and calculates the mixed layer eddy fluxes for buoyancy, heat and salinity. 1897 !! 1898 !! ** Method : 1899 !! 1900 !! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008 1901 !! Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008 1902 1903 REAL(wp), DIMENSION(jpi,jpj) :: hmle, zwb_fk, zvel_mle, zdiff_mle 1904 INTEGER :: ji, jj, jk ! dummy loop indices 1905 INTEGER :: ii, ij, ik ! local integers 1906 INTEGER , DIMENSION(jpi,jpj) :: inml_mle 1907 REAL(wp) :: zdb_mle, ztmp, zdbds_mle 1908 1909 mld_prof(:,:) = 4 1910 DO jk = 5, jpkm1 1911 DO jj = 2, jpjm1 1912 DO ji = 2, jpim1 1913 IF ( hmle(ji,jj) >= gdepw_n(ji,jj,jk) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk) 1914 END DO 1915 END DO 1916 END DO 1917 ! DO jj = 2, jpjm1 1918 ! DO ji = 1, jpim1 1919 ! zhmle(ji,jj) = gdepw_n(ji,jj,mld_prof(ji,jj)) 1920 ! END DO 1921 ! END DO 1922 ! Timestep mixed layer eddy depth. 1923 DO jj = 2, jpjm1 1924 DO ji = 2, jpim1 1925 zdb_mle = grav * (rhop(ji,jj,mld_prof(ji,jj)) - rhop(ji,jj,ibld(ji,jj) )) * r1_rau0 ! check ALMG 1926 IF ( lconv(ji,jj) .and. ( zdb_bl(ji,jj) < rn_osm_mle_thresh .and. mld_prof(ji,jj) > ibld(ji,jj) .and. zdb_mle > 0.0 ) ) THEN 1927 hmle(ji,jj) = hmle(ji,jj) + zwb0(ji,jj) * rn_rdt / MAX( zdb_mle, rn_osm_mle_thresh ) ! MLE layer deepening through encroachment. Don't have a good maximum value for deepening, so use threshold buoyancy. 1928 ELSE 1929 ! MLE layer relaxes towards mixed layer depth on timescale tau_mle, or tau_mle/10 1930 ! IF ( hmle(ji,jj) > zmld(ji,jj) ) THEN 1931 ! hmle(ji,jj) = hmle(ji,jj) - ( hmle(ji,jj) - zmld(ji,jj) ) * rn_rdt / rn_osm_mle_tau 1932 ! ELSE 1933 ! hmle(ji,jj) = hmle(ji,jj) - 10.0 * ( hmle(ji,jj) - zmld(ji,jj) ) * rn_rdt / rn_osm_mle_tau ! fast relaxation if MLE layer shallower than MLD 1934 ! ENDIF 1935 IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN 1936 hmle(ji,jj) = hmle(ji,jj) - ( hmle(ji,jj) - hbl(ji,jj) ) * rn_rdt / rn_osm_mle_tau 1937 ELSE 1938 hmle(ji,jj) = hmle(ji,jj) - 10.0 * ( hmle(ji,jj) - hbl(ji,jj) ) * rn_rdt / rn_osm_mle_tau ! fast relaxation if MLE layer shallower than MLD 1939 ENDIF 1940 ENDIF 1941 hmle(ji,jj) = MIN(hmle(ji,jj), ht_n(ji,jj)) 1942 hmle(ji,jj) = MIN(hmle(ji,jj), 1.2*zmld(ji,jj)) 1943 END DO 1944 END DO 1945 1946 mld_prof = 4 1947 DO jk = 5, jpkm1 1948 DO jj = 2, jpjm1 1949 DO ji = 2, jpim1 1950 IF ( hmle(ji,jj) >= gdepw_n(ji,jj,jk) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk) 1951 END DO 1952 END DO 1953 END DO 1954 DO jj = 2, jpjm1 1955 DO ji = 2, jpim1 1956 zhmle(ji,jj) = gdepw_n(ji,jj, mld_prof(ji,jj)) 1957 END DO 1958 END DO 1959 ! Calculate vertical buoyancy, heat and salinity fluxes due to MLE. 1960 1961 DO jj = 2, jpjm1 1962 DO ji = 2, jpim1 1963 IF ( lconv(ji,jj) ) THEN 1964 ztmp = r1_ft(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf 1965 ! zwt_fk(ji,jj) = 0.5_wp * ztmp * ( dbdx_mle(ji,jj) * zdtdx(ji,jj) + dbdy_mle(ji,jj) * zdtdy(ji,jj) & 1966 ! & + dbdx_mle(ji-1,jj) * zdtdx(ji-1,jj) + dbdy_mle(ji,jj-1) * zdtdy(ji,jj-1) ) 1967 ! zws_fk(ji,jj) = 0.5_wp * ztmp * ( dbdx_mle(ji,jj) * zdsdx(ji,jj) + dbdy_mle(ji,jj) * zdsdy(ji,jj) & 1968 ! & + dbdx_mle(ji-1,jj) * zdsdx(ji-1,jj) + dbdy_mle(ji,jj-1) * zdsdy(ji,jj-1) ) 1969 zdbds_mle = SQRT( 0.5_wp * ( dbdx_mle(ji,jj) * dbdx_mle(ji,jj) + dbdy_mle(ji,jj) * dbdy_mle(ji,jj) & 1970 & + dbdx_mle(ji-1,jj) * dbdx_mle(ji-1,jj) + dbdy_mle(ji,jj-1) * dbdy_mle(ji,jj-1) ) ) 1971 zwb_fk(ji,jj) = rn_osm_mle_ce * hmle(ji,jj) * hmle(ji,jj) *ztmp * zdbds_mle * zdbds_mle 1972 ! This vbelocity scale, defined in Fox-Kemper et al (2008), is needed for calculating dhdt. 1973 zvel_mle(ji,jj) = zdbds_mle * ztmp * hmle(ji,jj) * tmask(ji,jj,1) 1974 zdiff_mle(ji,jj) = 1.e-4_wp * ztmp * zdbds_mle * zhmle(ji,jj)**3 / rn_osm_mle_lf 1975 ENDIF 1976 END DO 1977 END DO 1978 END SUBROUTINE zdf_osm_mle_parameters 1979 1980 END SUBROUTINE zdf_osm 1364 1981 1365 1982 … … 1378 1995 INTEGER :: ios ! local integer 1379 1996 INTEGER :: ji, jj, jk ! dummy loop indices 1997 REAL z1_t2 1380 1998 !! 1381 1999 NAMELIST/namzdf_osm/ ln_use_osm_la, rn_osm_la, rn_osm_dstokes, nn_ave & 1382 2000 & ,nn_osm_wave, ln_dia_osm, rn_osm_hbl0 & 1383 & ,ln_kpprimix, rn_riinfty, rn_difri, ln_convmix, rn_difconv 2001 & ,ln_kpprimix, rn_riinfty, rn_difri, ln_convmix, rn_difconv, nn_osm_wave & 2002 & ,ln_osm_mle 2003 ! Namelist for Fox-Kemper parametrization. 2004 NAMELIST/namosm_mle/ nn_osm_mle, rn_osm_mle_ce, rn_osm_mle_lf, rn_osm_mle_time, rn_osm_mle_lat,& 2005 & rn_osm_mle_rho_c,rn_osm_mle_thresh,rn_osm_mle_tau 2006 1384 2007 !!---------------------------------------------------------------------- 1385 2008 ! … … 1397 2020 WRITE(numout,*) 'zdf_osm_init : OSMOSIS Parameterisation' 1398 2021 WRITE(numout,*) '~~~~~~~~~~~~' 1399 WRITE(numout,*) ' Namelist namzdf_osm : set tke mixing parameters' 1400 WRITE(numout,*) ' Use namelist rn_osm_la ln_use_osm_la = ', ln_use_osm_la 2022 WRITE(numout,*) ' Namelist namzdf_osm : set osm mixing parameters' 2023 WRITE(numout,*) ' Use rn_osm_la ln_use_osm_la = ', ln_use_osm_la 2024 WRITE(numout,*) ' Use MLE in OBL, i.e. Fox-Kemper param ln_osm_mle = ', ln_osm_mle 1401 2025 WRITE(numout,*) ' Turbulent Langmuir number rn_osm_la = ', rn_osm_la 1402 2026 WRITE(numout,*) ' Initial hbl for 1D runs rn_osm_hbl0 = ', rn_osm_hbl0 1403 WRITE(numout,*) ' Depth scale of Stokes drift rn_osm_dstokes = ', rn_osm_dstokes2027 WRITE(numout,*) ' Depth scale of Stokes drift rn_osm_dstokes = ', rn_osm_dstokes 1404 2028 WRITE(numout,*) ' horizontal average flag nn_ave = ', nn_ave 1405 2029 WRITE(numout,*) ' Stokes drift nn_osm_wave = ', nn_osm_wave … … 1423 2047 IF( zdf_osm_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_osm_init : unable to allocate arrays' ) 1424 2048 1425 call osm_rst( nit000, 'READ' ) !* read or initialize hbl 2049 2050 IF( ln_osm_mle ) THEN 2051 ! Initialise Fox-Kemper parametrization 2052 REWIND( numnam_ref ) ! Namelist namosm_mle in reference namelist : Tracer advection scheme 2053 READ ( numnam_ref, namosm_mle, IOSTAT = ios, ERR = 903) 2054 903 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namosm_mle in reference namelist') 2055 2056 REWIND( numnam_cfg ) ! Namelist namosm_mle in configuration namelist : Tracer advection scheme 2057 READ ( numnam_cfg, namosm_mle, IOSTAT = ios, ERR = 904 ) 2058 904 IF( ios > 0 ) CALL ctl_nam ( ios , 'namosm_mle in configuration namelist') 2059 IF(lwm) WRITE ( numond, namosm_mle ) 2060 2061 IF(lwp) THEN ! Namelist print 2062 WRITE(numout,*) 2063 WRITE(numout,*) 'zdf_osm_init : initialise mixed layer eddy (MLE)' 2064 WRITE(numout,*) '~~~~~~~~~~~~~' 2065 WRITE(numout,*) ' Namelist namosm_mle : ' 2066 WRITE(numout,*) ' MLE type: =0 standard Fox-Kemper ; =1 new formulation nn_osm_mle = ', nn_osm_mle 2067 WRITE(numout,*) ' magnitude of the MLE (typical value: 0.06 to 0.08) rn_osm_mle_ce = ', rn_osm_mle_ce 2068 WRITE(numout,*) ' scale of ML front (ML radius of deformation) (nn_osm_mle=0) rn_osm_mle_lf = ', rn_osm_mle_lf, 'm' 2069 WRITE(numout,*) ' maximum time scale of MLE (nn_osm_mle=0) rn_osm_mle_time = ', rn_osm_mle_time, 's' 2070 WRITE(numout,*) ' reference latitude (degrees) of MLE coef. (nn_osm_mle=1) rn_osm_mle_lat = ', rn_osm_mle_lat, 'deg' 2071 WRITE(numout,*) ' Density difference used to define ML for FK rn_osm_mle_rho_c = ', rn_osm_mle_rho_c 2072 WRITE(numout,*) ' Threshold used to define ML for FK rn_osm_mle_thresh = ', rn_osm_mle_thresh, 'm^2/s' 2073 WRITE(numout,*) ' Timescale for OSM-FK rn_osm_mle_tau = ', rn_osm_mle_tau, 's' 2074 ENDIF ! 2075 ENDIF 2076 ! 2077 IF(lwp) THEN 2078 WRITE(numout,*) 2079 IF( ln_osm_mle ) THEN 2080 WRITE(numout,*) ' ==>>> Mixed Layer Eddy induced transport added to OSMOSIS BL calculation' 2081 IF( nn_osm_mle == 0 ) WRITE(numout,*) ' Fox-Kemper et al 2010 formulation' 2082 IF( nn_osm_mle == 1 ) WRITE(numout,*) ' New formulation' 2083 ELSE 2084 WRITE(numout,*) ' ==>>> Mixed Layer induced transport NOT added to OSMOSIS BL calculation' 2085 ENDIF 2086 ENDIF 2087 ! 2088 IF( ln_osm_mle ) THEN ! MLE initialisation 2089 ! 2090 rb_c = grav * rn_osm_mle_rho_c /rau0 ! Mixed Layer buoyancy criteria 2091 IF(lwp) WRITE(numout,*) 2092 IF(lwp) WRITE(numout,*) ' ML buoyancy criteria = ', rb_c, ' m/s2 ' 2093 IF(lwp) WRITE(numout,*) ' associated ML density criteria defined in zdfmxl = ', rn_osm_mle_rho_c, 'kg/m3' 2094 ! 2095 IF( nn_osm_mle == 0 ) THEN ! MLE array allocation & initialisation ! 2096 ! 2097 ELSEIF( nn_osm_mle == 1 ) THEN ! MLE array allocation & initialisation 2098 rc_f = rn_osm_mle_ce/ ( 5.e3_wp * 2._wp * omega * SIN( rad * rn_osm_mle_lat ) ) 2099 ! 2100 ENDIF 2101 ! ! 1/(f^2+tau^2)^1/2 at t-point (needed in both nn_osm_mle case) 2102 z1_t2 = 2.e-5 2103 do jj=1,jpj 2104 do ji = 1,jpi 2105 r1_ft(ji,jj) = MIN(1./( ABS(ff_t(ji,jj)) + epsln ), ABS(ff_t(ji,jj))/z1_t2**2) 2106 end do 2107 end do 2108 ! z1_t2 = 1._wp / ( rn_osm_mle_time * rn_osm_mle_timeji,jj ) 2109 ! r1_ft(:,:) = 1._wp / SQRT( ff_t(:,:) * ff_t(:,:) + z1_t2 ) 2110 ! 2111 ENDIF 2112 2113 call osm_rst( nit000, 'READ' ) !* read or initialize hbl, dh, hmle 2114 1426 2115 1427 2116 IF( ln_zdfddm) THEN … … 1512 2201 CALL iom_set_rstw_var_active('wn') 1513 2202 CALL iom_set_rstw_var_active('hbl') 1514 CALL iom_set_rstw_var_active('hbli') 2203 CALL iom_set_rstw_var_active('dh') 2204 IF( ln_osm_mle ) THEN 2205 CALL iom_set_rstw_var_active('hmle') 2206 END IF 1515 2207 ENDIF 1516 2208 END SUBROUTINE zdf_osm_init … … 1530 2222 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 1531 2223 1532 INTEGER :: id1, id2 ! iom enquiry index2224 INTEGER :: id1, id2, id3 ! iom enquiry index 1533 2225 INTEGER :: ji, jj, jk ! dummy loop indices 1534 2226 INTEGER :: iiki, ikt ! local integer … … 1536 2228 REAL(wp) :: zN2_c ! local scalar 1537 2229 REAL(wp) :: rho_c = 0.01_wp !: density criterion for mixed layer depth 1538 INTEGER, DIMENSION( :,:), ALLOCATABLE:: imld_rst ! level of mixed-layer depth (pycnocline top)2230 INTEGER, DIMENSION(jpi,jpj) :: imld_rst ! level of mixed-layer depth (pycnocline top) 1539 2231 !!---------------------------------------------------------------------- 1540 2232 ! … … 1551 2243 WRITE(numout,*) ' ===>>>> : wn not in restart file, set to zero initially' 1552 2244 END IF 2245 1553 2246 id1 = iom_varid( numror, 'hbl' , ldstop = .FALSE. ) 1554 id2 = iom_varid( numror, ' hbli' , ldstop = .FALSE. )2247 id2 = iom_varid( numror, 'dh' , ldstop = .FALSE. ) 1555 2248 IF( id1 > 0 .AND. id2 > 0) THEN ! 'hbl' exists; read and return 1556 2249 CALL iom_get( numror, jpdom_autoglo, 'hbl' , hbl , ldxios = lrxios ) 1557 CALL iom_get( numror, jpdom_autoglo, 'hbli', hbli, ldxios = lrxios ) 1558 WRITE(numout,*) ' ===>>>> : hbl & hbli read from restart file' 2250 CALL iom_get( numror, jpdom_autoglo, 'dh', dh, ldxios = lrxios ) 2251 WRITE(numout,*) ' ===>>>> : hbl & dh read from restart file' 2252 IF( ln_osm_mle ) THEN 2253 id3 = iom_varid( numror, 'hmle' , ldstop = .FALSE. ) 2254 IF( id3 > 0) THEN 2255 CALL iom_get( numror, jpdom_autoglo, 'hmle' , hmle , ldxios = lrxios ) 2256 WRITE(numout,*) ' ===>>>> : hmle read from restart file' 2257 ELSE 2258 WRITE(numout,*) ' ===>>>> : hmle not found, set to hbl' 2259 hmle(:,:) = hbl(:,:) ! Initialise MLE depth. 2260 END IF 2261 END IF 1559 2262 RETURN 1560 ELSE ! 'hbl' & ' hbli' not in restart file, recalculate2263 ELSE ! 'hbl' & 'dh' not in restart file, recalculate 1561 2264 WRITE(numout,*) ' ===>>>> : previous run without osmosis scheme, hbl computed from stratification' 1562 2265 END IF … … 1568 2271 IF( TRIM(cdrw) == 'WRITE') THEN !* Write hbli into the restart file, then return 1569 2272 IF(lwp) WRITE(numout,*) '---- osm-rst ----' 1570 CALL iom_rstput( kt, nitrst, numrow, 'wn' , wn , ldxios = lwxios ) 1571 CALL iom_rstput( kt, nitrst, numrow, 'hbl' , hbl , ldxios = lwxios ) 1572 CALL iom_rstput( kt, nitrst, numrow, 'hbli' , hbli, ldxios = lwxios ) 2273 CALL iom_rstput( kt, nitrst, numrow, 'wn' , wn, ldxios = lwxios ) 2274 CALL iom_rstput( kt, nitrst, numrow, 'hbl' , hbl, ldxios = lwxios ) 2275 CALL iom_rstput( kt, nitrst, numrow, 'dh' , dh, ldxios = lwxios ) 2276 IF( ln_osm_mle ) THEN 2277 CALL iom_rstput( kt, nitrst, numrow, 'hmle', hmle, ldxios = lwxios ) 2278 END IF 1573 2279 RETURN 1574 2280 END IF … … 1578 2284 !!----------------------------------------------------------------------------- 1579 2285 IF( lwp ) WRITE(numout,*) ' ===>>>> : calculating hbl computed from stratification' 1580 ALLOCATE( imld_rst(jpi,jpj) )1581 2286 ! w-level of the mixing and mixed layers 1582 2287 CALL eos_rab( tsn, rab_n ) … … 1599 2304 DO jj = 1, jpj 1600 2305 DO ji = 1, jpi 1601 iiki = imld_rst(ji,jj) 1602 hbl (ji,jj) = gdepw_n(ji,jj,iiki ) * ssmask(ji,jj) ! Turbocline depth 2306 iiki = MAX(4,imld_rst(ji,jj)) 2307 hbl (ji,jj) = gdepw_n(ji,jj,iiki ) ! Turbocline depth 2308 dh (ji,jj) = e3t_n(ji,jj,iiki-1 ) ! Turbocline depth 1603 2309 END DO 1604 2310 END DO 1605 hbl = MAX(hbl,epsln) 1606 hbli(:,:) = hbl(:,:) 1607 DEALLOCATE( imld_rst ) 2311 1608 2312 WRITE(numout,*) ' ===>>>> : hbl computed from stratification' 2313 2314 IF( ln_osm_mle ) THEN 2315 hmle(:,:) = hbl(:,:) ! Initialise MLE depth. 2316 WRITE(numout,*) ' ===>>>> : hmle set = to hbl' 2317 END IF 2318 2319 wn(:,:,:) = 0._wp 2320 WRITE(numout,*) ' ===>>>> : wn not in restart file, set to zero initially' 1609 2321 END SUBROUTINE osm_rst 1610 2322 … … 1634 2346 ENDIF 1635 2347 1636 ! add non-local temperature and salinity flux1637 2348 DO jk = 1, jpkm1 1638 2349 DO jj = 2, jpjm1 … … 1648 2359 END DO 1649 2360 1650 1651 ! save the non-local tracer flux trends for diagnostic 2361 ! save the non-local tracer flux trends for diagnostics 1652 2362 IF( l_trdtra ) THEN 1653 2363 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 1654 2364 ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 1655 !!bug gm jpttdzdf ==> jpttosm 1656 CALL trd_tra( kt, 'TRA', jp_tem, jptra_ zdf, ztrdt )1657 CALL trd_tra( kt, 'TRA', jp_sal, jptra_ zdf, ztrds )2365 2366 CALL trd_tra( kt, 'TRA', jp_tem, jptra_osm, ztrdt ) 2367 CALL trd_tra( kt, 'TRA', jp_sal, jptra_osm, ztrds ) 1658 2368 DEALLOCATE( ztrdt ) ; DEALLOCATE( ztrds ) 1659 2369 ENDIF … … 1723 2433 1724 2434 !!====================================================================== 2435 1725 2436 END MODULE zdfosm -
NEMO/branches/UKMO/NEMO_4.0.2_FKOSM/src/OCE/ZDF/zdfphy.F90
r12658 r12912 56 56 !!---------------------------------------------------------------------- 57 57 !! NEMO/OCE 4.0 , NEMO Consortium (2018) 58 !! $Id $58 !! $Id: zdfphy.F90 12178 2019-12-11 11:02:38Z agn $ 59 59 !! Software governed by the CeCILL license (see ./LICENSE) 60 60 !!---------------------------------------------------------------------- … … 172 172 IF( ln_zdfosm .AND. ln_zdfevd ) CALL ctl_stop( 'zdf_phy_init: chose between ln_zdfosm and ln_zdfevd' ) 173 173 IF( lk_top .AND. ln_zdfnpc ) CALL ctl_stop( 'zdf_phy_init: npc scheme is not working with key_top' ) 174 IF( lk_top .AND. ln_zdfosm ) CALL ctl_ stop( 'zdf_phy_init: osmosis scheme is not working with key_top' )174 IF( lk_top .AND. ln_zdfosm ) CALL ctl_warn( 'zdf_phy_init: osmosis gives no non-local fluxes for TOP tracers yet' ) 175 175 IF(lwp) THEN 176 176 WRITE(numout,*)
Note: See TracChangeset
for help on using the changeset viewer.
