Changeset 13867
- Timestamp:
- 2020-11-24T21:32:46+01:00 (4 years ago)
- Location:
- NEMO/branches/2020/dev_r13787_OSMOSIS_IMMERSE/src/OCE
- Files:
-
- 5 edited
-
DIA/diawri.F90 (modified) (2 diffs)
-
TRA/traatf.F90 (modified) (1 diff)
-
TRA/tramle.F90 (modified) (3 diffs)
-
TRD/trd_oce.F90 (modified) (2 diffs)
-
ZDF/zdfosm.F90 (modified) (53 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r13787_OSMOSIS_IMMERSE/src/OCE/DIA/diawri.F90
r13497 r13867 46 46 USE zdfdrg ! ocean vertical physics: top/bottom friction 47 47 USE zdfmxl ! mixed layer 48 USE zdfosm ! mixed layer 48 49 ! 49 50 USE lbclnk ! ocean lateral boundary conditions (or mpp link) … … 1076 1077 CALL iom_rstput ( 0, 0, inum, "qz1_abl", tq_abl(:,:,2,nt_a,2) ) ! now first level humidity 1077 1078 ENDIF 1079 IF( ln_zdfosm ) THEN 1080 CALL iom_rstput( 0, 0, inum, 'hbl', hbl*tmask(:,:,1) ) ! now boundary-layer depth 1081 CALL iom_rstput( 0, 0, inum, 'hml', hml*tmask(:,:,1) ) ! now mixed-layer depth 1082 CALL iom_rstput( 0, 0, inum, 'avt_k', avt_k*wmask ) ! w-level diffusion 1083 CALL iom_rstput( 0, 0, inum, 'avm_k', avm_k*wmask ) ! now w-level viscosity 1084 CALL iom_rstput( 0, 0, inum, 'ghamt', ghamt*wmask ) ! non-local t forcing 1085 CALL iom_rstput( 0, 0, inum, 'ghams', ghams*wmask ) ! non-local s forcing 1086 CALL iom_rstput( 0, 0, inum, 'ghamu', ghamu*umask ) ! non-local u forcing 1087 CALL iom_rstput( 0, 0, inum, 'ghamv', ghamv*vmask ) ! non-local v forcing 1088 IF( ln_osm_mle ) THEN 1089 CALL iom_rstput( 0, 0, inum, 'hmle', hmle*tmask(:,:,1) ) ! now transition-layer depth 1090 END IF 1091 ENDIF 1078 1092 ! 1079 1093 CALL iom_close( inum ) -
NEMO/branches/2020/dev_r13787_OSMOSIS_IMMERSE/src/OCE/TRA/traatf.F90
r13295 r13867 117 117 IF( l_trdtra ) THEN 118 118 ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 119 ztrdt(:,:, jpk) = 0._wp120 ztrds(:,:, jpk) = 0._wp119 ztrdt(:,:,:) = 0._wp 120 ztrds(:,:,:) = 0._wp 121 121 IF( ln_traldf_iso ) THEN ! diagnose the "pure" Kz diffusive trend 122 122 CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_tem, jptra_zdfp, ztrdt ) -
NEMO/branches/2020/dev_r13787_OSMOSIS_IMMERSE/src/OCE/TRA/tramle.F90
r13497 r13867 20 20 USE lib_mpp ! MPP library 21 21 USE lbclnk ! lateral boundary condition / mpp link 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 22 26 23 27 IMPLICIT NONE … … 96 100 !!---------------------------------------------------------------------- 97 101 ! 98 ! !== MLD used for MLE ==! 99 ! ! compute from the 10m density to deal with the diurnal cycle 100 inml_mle(:,:) = mbkt(:,:) + 1 ! init. to number of ocean w-level (T-level + 1) 101 IF ( nla10 > 0 ) THEN ! avoid case where first level is thicker than 10m 102 DO_3DS( 1, 1, 1, 1, jpkm1, nlb10, -1 ) ! from the bottom to nlb10 (10m) 103 IF( rhop(ji,jj,jk) > rhop(ji,jj,nla10) + rn_rho_c_mle ) inml_mle(ji,jj) = jk ! Mixed layer 102 ! 103 IF(ln_osm_mle.and.ln_zdfosm) THEN 104 ikmax = MIN( MAXVAL( mld_prof(:,:) ), jpkm1 ) ! max level of the computation 105 ! 106 ! 107 SELECT CASE( nn_mld_uv ) ! MLD at u- & v-pts 108 CASE ( 0 ) != min of the 2 neighbour MLDs 109 DO_2D( 1, 0, 1, 0 ) 110 zhu(ji,jj) = MIN( hmle(ji+1,jj), hmle(ji,jj) ) 111 zhv(ji,jj) = MIN( hmle(ji,jj+1), hmle(ji,jj) ) 112 END_2D 113 CASE ( 1 ) != average of the 2 neighbour MLDs 114 DO_2D( 1, 0, 1, 0 ) 115 zhu(ji,jj) = MAX( hmle(ji+1,jj), hmle(ji,jj) ) 116 zhv(ji,jj) = MAX( hmle(ji,jj+1), hmle(ji,jj) ) 117 END_2D 118 CASE ( 2 ) != max of the 2 neighbour MLDs 119 DO_2D( 1, 0, 1, 0 ) 120 zhu(ji,jj) = MAX( hmle(ji+1,jj), hmle(ji,jj) ) 121 zhv(ji,jj) = MAX( hmle(ji,jj+1), hmle(ji,jj) ) 122 END_2D 123 END SELECT 124 IF( nn_mle == 0 ) THEN ! Fox-Kemper et al. 2010 formulation 125 DO_2D( 1, 0, 1, 0 ) 126 zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj) * e2u(ji,jj) & 127 & * dbdx_mle(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) & 128 & / ( MAX( rn_lf * rfu(ji,jj) , SQRT( rb_c * zhu(ji,jj) ) ) ) 129 ! 130 zpsim_v(ji,jj) = rn_ce * zhv(ji,jj) * zhv(ji,jj) * e1v(ji,jj) & 131 & * dbdy_mle(ji,jj) * MIN( 111.e3_wp , e2v(ji,jj) ) & 132 & / ( MAX( rn_lf * rfv(ji,jj) , SQRT( rb_c * zhv(ji,jj) ) ) ) 133 END_2D 134 ! 135 ELSEIF( nn_mle == 1 ) THEN ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat) 136 DO_2D( 1, 0, 1, 0 ) 137 zpsim_u(ji,jj) = rc_f * zhu(ji,jj) * zhu(ji,jj) * e2u(ji,jj) & 138 & * dbdx_mle(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) 139 ! 140 zpsim_v(ji,jj) = rc_f * zhv(ji,jj) * zhv(ji,jj) * e1v(ji,jj) & 141 & * dbdy_mle(ji,jj) * MIN( 111.e3_wp , e2v(ji,jj) ) 142 END_2D 143 ENDIF 144 145 ELSE !do not use osn_mle 146 ! !== MLD used for MLE ==! 147 ! ! compute from the 10m density to deal with the diurnal cycle 148 inml_mle(:,:) = mbkt(:,:) + 1 ! init. to number of ocean w-level (T-level + 1) 149 IF ( nla10 > 0 ) THEN ! avoid case where first level is thicker than 10m 150 DO_3DS( 1, 1, 1, 1, jpkm1, nlb10, -1 ) ! from the bottom to nlb10 (10m) 151 IF( rhop(ji,jj,jk) > rhop(ji,jj,nla10) + rn_rho_c_mle ) inml_mle(ji,jj) = jk ! Mixed layer 152 END_3D 153 ENDIF 154 ikmax = MIN( MAXVAL( inml_mle(:,:) ), jpkm1 ) ! max level of the computation 155 ! 156 ! 157 zmld(:,:) = 0._wp !== Horizontal shape of the MLE ==! 158 zbm (:,:) = 0._wp 159 zn2 (:,:) = 0._wp 160 DO_3D( 1, 1, 1, 1, 1, ikmax ) ! MLD and mean buoyancy and N2 over the mixed layer 161 zc = e3t(ji,jj,jk,Kmm) * REAL( MIN( MAX( 0, inml_mle(ji,jj)-jk ) , 1 ) ) ! zc being 0 outside the ML t-points 162 zmld(ji,jj) = zmld(ji,jj) + zc 163 zbm (ji,jj) = zbm (ji,jj) + zc * (rho0 - rhop(ji,jj,jk) ) * r1_rho0 164 zn2 (ji,jj) = zn2 (ji,jj) + zc * (rn2(ji,jj,jk)+rn2(ji,jj,jk+1))*0.5_wp 104 165 END_3D 105 ENDIF 106 ikmax = MIN( MAXVAL( inml_mle(:,:) ), jpkm1 ) ! max level of the computation 107 ! 108 ! 109 zmld(:,:) = 0._wp !== Horizontal shape of the MLE ==! 110 zbm (:,:) = 0._wp 111 zn2 (:,:) = 0._wp 112 DO_3D( 1, 1, 1, 1, 1, ikmax ) ! MLD and mean buoyancy and N2 over the mixed layer 113 zc = e3t(ji,jj,jk,Kmm) * REAL( MIN( MAX( 0, inml_mle(ji,jj)-jk ) , 1 ) ) ! zc being 0 outside the ML t-points 114 zmld(ji,jj) = zmld(ji,jj) + zc 115 zbm (ji,jj) = zbm (ji,jj) + zc * (rho0 - rhop(ji,jj,jk) ) * r1_rho0 116 zn2 (ji,jj) = zn2 (ji,jj) + zc * (rn2(ji,jj,jk)+rn2(ji,jj,jk+1))*0.5_wp 117 END_3D 118 119 SELECT CASE( nn_mld_uv ) ! MLD at u- & v-pts 120 CASE ( 0 ) != min of the 2 neighbour MLDs 121 DO_2D( 1, 0, 1, 0 ) 122 zhu(ji,jj) = MIN( zmld(ji+1,jj), zmld(ji,jj) ) 123 zhv(ji,jj) = MIN( zmld(ji,jj+1), zmld(ji,jj) ) 124 END_2D 125 CASE ( 1 ) != average of the 2 neighbour MLDs 126 DO_2D( 1, 0, 1, 0 ) 127 zhu(ji,jj) = ( zmld(ji+1,jj) + zmld(ji,jj) ) * 0.5_wp 128 zhv(ji,jj) = ( zmld(ji,jj+1) + zmld(ji,jj) ) * 0.5_wp 129 END_2D 130 CASE ( 2 ) != max of the 2 neighbour MLDs 131 DO_2D( 1, 0, 1, 0 ) 132 zhu(ji,jj) = MAX( zmld(ji+1,jj), zmld(ji,jj) ) 133 zhv(ji,jj) = MAX( zmld(ji,jj+1), zmld(ji,jj) ) 134 END_2D 135 END SELECT 136 ! ! convert density into buoyancy 137 zbm(:,:) = + grav * zbm(:,:) / MAX( e3t(:,:,1,Kmm), zmld(:,:) ) 138 ! 139 ! 140 ! !== Magnitude of the MLE stream function ==! 141 ! 142 ! di[bm] Ds 143 ! Psi = Ce H^2 ---------------- e2u mu(z) where fu Lf = MAX( fu*rn_fl , (Db H)^1/2 ) 144 ! e1u Lf fu and the e2u for the "transport" 145 ! (not *e3u as divided by e3u at the end) 146 ! 147 IF( nn_mle == 0 ) THEN ! Fox-Kemper et al. 2010 formulation 148 DO_2D( 1, 0, 1, 0 ) 149 zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj) * e2_e1u(ji,jj) & 150 & * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) ) & 151 & / ( MAX( rn_lf * rfu(ji,jj) , SQRT( rb_c * zhu(ji,jj) ) ) ) 166 167 SELECT CASE( nn_mld_uv ) ! MLD at u- & v-pts 168 CASE ( 0 ) != min of the 2 neighbour MLDs 169 DO_2D( 1, 0, 1, 0 ) 170 zhu(ji,jj) = MIN( zmld(ji+1,jj), zmld(ji,jj) ) 171 zhv(ji,jj) = MIN( zmld(ji,jj+1), zmld(ji,jj) ) 172 END_2D 173 CASE ( 1 ) != average of the 2 neighbour MLDs 174 DO_2D( 1, 0, 1, 0 ) 175 zhu(ji,jj) = ( zmld(ji+1,jj) + zmld(ji,jj) ) * 0.5_wp 176 zhv(ji,jj) = ( zmld(ji,jj+1) + zmld(ji,jj) ) * 0.5_wp 177 END_2D 178 CASE ( 2 ) != max of the 2 neighbour MLDs 179 DO_2D( 1, 0, 1, 0 ) 180 zhu(ji,jj) = MAX( zmld(ji+1,jj), zmld(ji,jj) ) 181 zhv(ji,jj) = MAX( zmld(ji,jj+1), zmld(ji,jj) ) 182 END_2D 183 END SELECT 184 ! ! convert density into buoyancy 185 zbm(:,:) = + grav * zbm(:,:) / MAX( e3t(:,:,1,Kmm), zmld(:,:) ) 186 ! 187 ! 188 ! !== Magnitude of the MLE stream function ==! 189 ! 190 ! di[bm] Ds 191 ! Psi = Ce H^2 ---------------- e2u mu(z) where fu Lf = MAX( fu*rn_fl , (Db H)^1/2 ) 192 ! e1u Lf fu and the e2u for the "transport" 193 ! (not *e3u as divided by e3u at the end) 194 ! 195 IF( nn_mle == 0 ) THEN ! Fox-Kemper et al. 2010 formulation 196 DO_2D( 1, 0, 1, 0 ) 197 zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj) * e2_e1u(ji,jj) & 198 & * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) ) & 199 & / ( MAX( rn_lf * rfu(ji,jj) , SQRT( rb_c * zhu(ji,jj) ) ) ) 152 200 ! 153 zpsim_v(ji,jj) = rn_ce * zhv(ji,jj) * zhv(ji,jj) * e1_e2v(ji,jj) &154 & * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) ) &155 & / ( MAX( rn_lf * rfv(ji,jj) , SQRT( rb_c * zhv(ji,jj) ) ) )156 END_2D157 !158 ELSEIF( nn_mle == 1 ) THEN ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat)159 DO_2D( 1, 0, 1, 0 )160 zpsim_u(ji,jj) = rc_f * zhu(ji,jj) * zhu(ji,jj) * e2_e1u(ji,jj) &161 & * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) )201 zpsim_v(ji,jj) = rn_ce * zhv(ji,jj) * zhv(ji,jj) * e1_e2v(ji,jj) & 202 & * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) ) & 203 & / ( MAX( rn_lf * rfv(ji,jj) , SQRT( rb_c * zhv(ji,jj) ) ) ) 204 END_2D 205 ! 206 ELSEIF( nn_mle == 1 ) THEN ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat) 207 DO_2D( 1, 0, 1, 0 ) 208 zpsim_u(ji,jj) = rc_f * zhu(ji,jj) * zhu(ji,jj) * e2_e1u(ji,jj) & 209 & * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) ) 162 210 ! 163 zpsim_v(ji,jj) = rc_f * zhv(ji,jj) * zhv(ji,jj) * e1_e2v(ji,jj) & 164 & * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) ) 165 END_2D 166 ENDIF 167 ! 168 IF( nn_conv == 1 ) THEN ! No MLE in case of convection 169 DO_2D( 1, 0, 1, 0 ) 170 IF( MIN( zn2(ji,jj) , zn2(ji+1,jj) ) < 0._wp ) zpsim_u(ji,jj) = 0._wp 171 IF( MIN( zn2(ji,jj) , zn2(ji,jj+1) ) < 0._wp ) zpsim_v(ji,jj) = 0._wp 172 END_2D 173 ENDIF 174 ! 175 ! !== structure function value at uw- and vw-points ==! 176 DO_2D( 1, 0, 1, 0 ) 177 zhu(ji,jj) = 1._wp / zhu(ji,jj) ! hu --> 1/hu 178 zhv(ji,jj) = 1._wp / zhv(ji,jj) 179 END_2D 180 ! 181 zpsi_uw(:,:,:) = 0._wp 182 zpsi_vw(:,:,:) = 0._wp 183 ! 211 zpsim_v(ji,jj) = rc_f * zhv(ji,jj) * zhv(ji,jj) * e1_e2v(ji,jj) & 212 & * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) ) 213 END_2D 214 ENDIF 215 ! 216 IF( nn_conv == 1 ) THEN ! No MLE in case of convection 217 DO_2D( 1, 0, 1, 0 ) 218 IF( MIN( zn2(ji,jj) , zn2(ji+1,jj) ) < 0._wp ) zpsim_u(ji,jj) = 0._wp 219 IF( MIN( zn2(ji,jj) , zn2(ji,jj+1) ) < 0._wp ) zpsim_v(ji,jj) = 0._wp 220 END_2D 221 ENDIF 222 ! 223 ENDIF ! end of ln_osm_mle conditional 224 ! !== structure function value at uw- and vw-points ==! 225 DO_2D( 1, 0, 1, 0 ) 226 zhu(ji,jj) = 1._wp / MAX(zhu(ji,jj), rsmall) ! hu --> 1/hu 227 zhv(ji,jj) = 1._wp / MAX(zhv(ji,jj), rsmall) 228 END_2D 229 ! 230 zpsi_uw(:,:,:) = 0._wp 231 zpsi_vw(:,:,:) = 0._wp 232 ! 184 233 DO_3D( 1, 0, 1, 0, 2, ikmax ) ! start from 2 : surface value = 0 185 234 zcuw = 1._wp - ( gdepw(ji+1,jj,jk,Kmm) + gdepw(ji,jj,jk,Kmm) ) * zhu(ji,jj) … … 206 255 END DO 207 256 208 IF( cdtype == 'TRA') THEN !== outputs ==! 209 ! 210 zLf_NH(:,:) = SQRT( rb_c * zmld(:,:) ) * r1_ft(:,:) ! Lf = N H / f 211 CALL iom_put( "Lf_NHpf" , zLf_NH ) ! Lf = N H / f 212 ! 213 ! divide by cross distance to give streamfunction with dimensions m^2/s 214 DO jk = 1, ikmax+1 215 zpsi_uw(:,:,jk) = zpsi_uw(:,:,jk) * r1_e2u(:,:) 216 zpsi_vw(:,:,jk) = zpsi_vw(:,:,jk) * r1_e1v(:,:) 217 END DO 218 CALL iom_put( "psiu_mle", zpsi_uw ) ! i-mle streamfunction 219 CALL iom_put( "psiv_mle", zpsi_vw ) ! j-mle streamfunction 220 ENDIF 221 ! 222 END SUBROUTINE tra_mle_trp 257 IF( cdtype == 'TRA') THEN !== outputs ==! 258 ! 259 IF (ln_osm_mle.and.ln_zdfosm) THEN 260 zLf_NH(:,:) = SQRT( rb_c * hmle(:,:) ) * r1_ft(:,:) ! Lf = N H / f 261 ELSE 262 zLf_NH(:,:) = SQRT( rb_c * zmld(:,:) ) * r1_ft(:,:) ! Lf = N H / f 263 END IF 264 CALL iom_put( "Lf_NHpf" , zLf_NH ) ! Lf = N H / f 265 ! 266 ! divide by cross distance to give streamfunction with dimensions m^2/s 267 DO jk = 1, ikmax+1 268 zpsi_uw(:,:,jk) = zpsi_uw(:,:,jk) * r1_e2u(:,:) 269 zpsi_vw(:,:,jk) = zpsi_vw(:,:,jk) * r1_e1v(:,:) 270 END DO 271 CALL iom_put( "psiu_mle", zpsi_uw ) ! i-mle streamfunction 272 CALL iom_put( "psiv_mle", zpsi_vw ) ! j-mle streamfunction 273 ENDIF 274 ! 275 END SUBROUTINE tra_mle_trp 223 276 224 277 -
NEMO/branches/2020/dev_r13787_OSMOSIS_IMMERSE/src/OCE/TRD/trd_oce.F90
r10068 r13867 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/2020/dev_r13787_OSMOSIS_IMMERSE/src/OCE/ZDF/zdfosm.F90
r13497 r13867 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 105 REAL(wp) :: rn_zdfosm_adjust_sd = 1.0 ! factor to reduce Stokes drift by 106 REAL(wp) :: rn_osm_hblfrac = 0.1! for nn_osm_wave = 3/4 specify fraction in top of hbl 107 LOGICAL :: ln_zdfosm_ice_shelter ! flag to activate ice sheltering 86 108 REAL(wp) :: rn_osm_hbl0 = 10._wp ! Initial value of hbl for 1D runs 87 109 INTEGER :: nn_ave ! = 0/1 flag for horizontal average on avt 88 110 INTEGER :: nn_osm_wave = 0 ! = 0/1/2 flag for getting stokes drift from La# / PM wind-waves/Inputs into sbcwave 111 INTEGER :: nn_osm_SD_reduce ! = 0/1/2 flag for getting effective stokes drift from surface value 89 112 LOGICAL :: ln_dia_osm ! Use namelist rn_osm_la 90 113 … … 96 119 REAL(wp) :: rn_difconv = 1._wp ! diffusivity when unstable below BL (m2/s) 97 120 121 ! OSMOSIS mixed layer eddy parametrization constants 122 INTEGER :: nn_osm_mle ! = 0/1 flag for horizontal average on avt 123 REAL(wp) :: rn_osm_mle_ce ! MLE coefficient 124 ! ! parameters used in nn_osm_mle = 0 case 125 REAL(wp) :: rn_osm_mle_lf ! typical scale of mixed layer front 126 REAL(wp) :: rn_osm_mle_time ! time scale for mixing momentum across the mixed layer 127 ! ! parameters used in nn_osm_mle = 1 case 128 REAL(wp) :: rn_osm_mle_lat ! reference latitude for a 5 km scale of ML front 129 LOGICAL :: ln_osm_hmle_limit ! If true arbitrarily restrict hmle to rn_osm_hmle_limit*zmld 130 REAL(wp) :: rn_osm_hmle_limit ! If ln_osm_hmle_limit true arbitrarily restrict hmle to rn_osm_hmle_limit*zmld 131 REAL(wp) :: rn_osm_mle_rho_c ! Density criterion for definition of MLD used by FK 132 REAL(wp) :: r5_21 = 5.e0 / 21.e0 ! factor used in mle streamfunction computation 133 REAL(wp) :: rb_c ! ML buoyancy criteria = g rho_c /rho0 where rho_c is defined in zdfmld 134 REAL(wp) :: rc_f ! MLE coefficient (= rn_ce / (5 km * fo) ) in nn_osm_mle=1 case 135 REAL(wp) :: rn_osm_mle_thresh ! Threshold buoyancy for deepening of MLE layer below OSBL base. 136 REAL(wp) :: rn_osm_bl_thresh ! Threshold buoyancy for deepening of OSBL base. 137 REAL(wp) :: rn_osm_mle_tau ! Adjustment timescale for MLE. 138 139 98 140 ! !!! ** General constants ** 99 REAL(wp) :: epsln = 1.0e-20_wp ! a small positive number 141 REAL(wp) :: epsln = 1.0e-20_wp ! a small positive number to ensure no div by zero 142 REAL(wp) :: depth_tol = 1.0e-6_wp ! a small-ish positive number to give a hbl slightly shallower than gdepw 100 143 REAL(wp) :: pthird = 1._wp/3._wp ! 1/3 101 144 REAL(wp) :: p2third = 2._wp/3._wp ! 2/3 … … 118 161 !! *** FUNCTION zdf_osm_alloc *** 119 162 !!---------------------------------------------------------------------- 120 ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk), ghams(jpi,jpj,jpk), & 121 & hbl(jpi,jpj) , hbli(jpi,jpj) , dstokes(jpi, jpj) , & 122 & etmean(jpi,jpj,jpk), STAT= zdf_osm_alloc ) 163 ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk),ghams(jpi,jpj,jpk), & 164 & hbl(jpi,jpj), dh(jpi,jpj), hml(jpi,jpj), dstokes(jpi, jpj), & 165 & etmean(jpi,jpj,jpk), STAT= zdf_osm_alloc ) 166 167 ALLOCATE( hmle(jpi,jpj), r1_ft(jpi,jpj), dbdx_mle(jpi,jpj), dbdy_mle(jpi,jpj), & 168 & mld_prof(jpi,jpj), STAT= zdf_osm_alloc ) 169 170 CALL mpp_sum ( 'zdfosm', zdf_osm_alloc ) 123 171 IF( zdf_osm_alloc /= 0 ) CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm arrays') 124 CALL mpp_sum ( 'zdfosm', zdf_osm_alloc ) 172 125 173 END FUNCTION zdf_osm_alloc 126 174 … … 166 214 !! 167 215 INTEGER :: ji, jj, jk ! dummy loop indices 216 217 INTEGER :: jl ! dummy loop indices 218 168 219 INTEGER :: ikbot, jkmax, jkm1, jkp2 ! 169 220 … … 196 247 REAL(wp), DIMENSION(jpi,jpj) :: zwbav ! Buoyancy flux - bl average 197 248 REAL(wp), DIMENSION(jpi,jpj) :: zwb_ent ! Buoyancy entrainment flux 249 REAL(wp), DIMENSION(jpi,jpj) :: zwb_min 250 251 252 REAL(wp), DIMENSION(jpi,jpj) :: zwb_fk_b ! MLE buoyancy flux averaged over OSBL 253 REAL(wp), DIMENSION(jpi,jpj) :: zwb_fk ! max MLE buoyancy flux 254 REAL(wp), DIMENSION(jpi,jpj) :: zdiff_mle ! extra MLE vertical diff 255 REAL(wp), DIMENSION(jpi,jpj) :: zvel_mle ! velocity scale for dhdt with stable ML and FK 256 198 257 REAL(wp), DIMENSION(jpi,jpj) :: zustke ! Surface Stokes drift 199 258 REAL(wp), DIMENSION(jpi,jpj) :: zla ! Trubulent Langmuir number … … 201 260 REAL(wp), DIMENSION(jpi,jpj) :: zsin_wind ! Sin angle of surface stress 202 261 REAL(wp), DIMENSION(jpi,jpj) :: zhol ! Stability parameter for boundary layer 203 LOGICAL, DIMENSION(:,:), ALLOCATABLE :: lconv ! unstable/stable bl 262 LOGICAL, DIMENSION(jpi,jpj) :: lconv ! unstable/stable bl 263 LOGICAL, DIMENSION(jpi,jpj) :: lshear ! Shear layers 264 LOGICAL, DIMENSION(jpi,jpj) :: lpyc ! OSBL pycnocline present 265 LOGICAL, DIMENSION(jpi,jpj) :: lflux ! surface flux extends below OSBL into MLE layer. 266 LOGICAL, DIMENSION(jpi,jpj) :: lmle ! MLE layer increases in hickness. 204 267 205 268 ! mixed-layer variables … … 207 270 INTEGER, DIMENSION(jpi,jpj) :: ibld ! level of boundary layer base 208 271 INTEGER, DIMENSION(jpi,jpj) :: imld ! level of mixed-layer depth (pycnocline top) 272 INTEGER, DIMENSION(jpi,jpj) :: jp_ext, jp_ext_mle ! offset for external level 273 INTEGER, DIMENSION(jpi, jpj) :: j_ddh ! Type of shear layer 209 274 210 275 REAL(wp) :: ztgrad,zsgrad,zbgrad ! Temporary variables used to calculate pycnocline gradients … … 213 278 REAL(wp), DIMENSION(jpi,jpj) :: zhbl ! bl depth - grid 214 279 REAL(wp), DIMENSION(jpi,jpj) :: zhml ! ml depth - grid 280 281 REAL(wp), DIMENSION(jpi,jpj) :: zhmle ! MLE depth - grid 282 REAL(wp), DIMENSION(jpi,jpj) :: zmld ! ML depth on grid 283 215 284 REAL(wp), DIMENSION(jpi,jpj) :: zdh ! pycnocline depth - grid 216 285 REAL(wp), DIMENSION(jpi,jpj) :: zdhdt ! BL depth tendency 217 REAL(wp), DIMENSION(jpi,jpj) :: zt_bl,zs_bl,zu_bl,zv_bl,zrh_bl ! averages over the depth of the blayer 218 REAL(wp), DIMENSION(jpi,jpj) :: zt_ml,zs_ml,zu_ml,zv_ml,zrh_ml ! averages over the depth of the mixed layer 219 REAL(wp), DIMENSION(jpi,jpj) :: zdt_bl,zds_bl,zdu_bl,zdv_bl,zdrh_bl,zdb_bl ! difference between blayer average and parameter at base of blayer 220 REAL(wp), DIMENSION(jpi,jpj) :: zdt_ml,zds_ml,zdu_ml,zdv_ml,zdrh_ml,zdb_ml ! difference between mixed layer average and parameter at base of blayer 221 REAL(wp), DIMENSION(jpi,jpj) :: zwth_ent,zws_ent ! heat and salinity fluxes at the top of the pycnocline 222 REAL(wp), DIMENSION(jpi,jpj) :: zuw_bse,zvw_bse ! momentum fluxes at the top of the pycnocline 286 REAL(wp), DIMENSION(jpi,jpj) :: zddhdt ! correction to dhdt due to internal structure. 287 REAL(wp), DIMENSION(jpi,jpj) :: zdtdz_bl_ext,zdsdz_bl_ext,zdbdz_bl_ext ! external temperature/salinity and buoyancy gradients 288 REAL(wp), DIMENSION(jpi,jpj) :: zdtdz_mle_ext,zdsdz_mle_ext,zdbdz_mle_ext ! external temperature/salinity and buoyancy gradients 289 REAL(wp), DIMENSION(jpi,jpj) :: zdtdx, zdtdy, zdsdx, zdsdy ! horizontal gradients for Fox-Kemper parametrization. 290 291 REAL(wp), DIMENSION(jpi,jpj) :: zt_bl,zs_bl,zu_bl,zv_bl,zb_bl ! averages over the depth of the blayer 292 REAL(wp), DIMENSION(jpi,jpj) :: zt_ml,zs_ml,zu_ml,zv_ml,zb_ml ! averages over the depth of the mixed layer 293 REAL(wp), DIMENSION(jpi,jpj) :: zt_mle,zs_mle,zu_mle,zv_mle,zb_mle ! averages over the depth of the MLE layer 294 REAL(wp), DIMENSION(jpi,jpj) :: zdt_bl,zds_bl,zdu_bl,zdv_bl,zdb_bl ! difference between blayer average and parameter at base of blayer 295 REAL(wp), DIMENSION(jpi,jpj) :: zdt_ml,zds_ml,zdu_ml,zdv_ml,zdb_ml ! difference between mixed layer average and parameter at base of blayer 296 REAL(wp), DIMENSION(jpi,jpj) :: zdt_mle,zds_mle,zdu_mle,zdv_mle,zdb_mle ! difference between MLE layer average and parameter at base of blayer 297 ! REAL(wp), DIMENSION(jpi,jpj) :: zwth_ent,zws_ent ! heat and salinity fluxes at the top of the pycnocline 298 REAL(wp) :: zwth_ent,zws_ent ! heat and salinity fluxes at the top of the pycnocline 299 REAL(wp) :: zuw_bse,zvw_bse ! momentum fluxes at the top of the pycnocline 223 300 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdtdz_pyc ! parametrized gradient of temperature in pycnocline 224 301 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdsdz_pyc ! parametrised gradient of salinity in pycnocline … … 226 303 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdudz_pyc ! u-shear across the pycnocline 227 304 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdvdz_pyc ! v-shear across the pycnocline 228 305 REAL(wp), DIMENSION(jpi,jpj) :: zdbds_mle ! Magnitude of horizontal buoyancy gradient. 229 306 ! Flux-gradient relationship variables 307 REAL(wp), DIMENSION(jpi, jpj) :: zshear, zri_i ! Shear production and interfacial richardon number. 230 308 231 309 REAL(wp) :: zl_c,zl_l,zl_eps ! Used to calculate turbulence length scale. 232 310 233 REAL(wp) , DIMENSION(jpi,jpj) :: zdifml_sc,zvisml_sc,zdifpyc_sc,zvispyc_sc,zbeta_d_sc,zbeta_v_sc ! Scales for eddy diffusivity/viscosity311 REAL(wp) :: za_cubic, zb_cubic, zc_cubic, zd_cubic ! coefficients in cubic polynomial specifying diffusivity in pycnocline. 234 312 REAL(wp), DIMENSION(jpi,jpj) :: zsc_wth_1,zsc_ws_1 ! Temporary scales used to calculate scalar non-gradient terms. 313 REAL(wp), DIMENSION(jpi,jpj) :: zsc_wth_pyc, zsc_ws_pyc ! Scales for pycnocline transport term/ 235 314 REAL(wp), DIMENSION(jpi,jpj) :: zsc_uw_1,zsc_uw_2,zsc_vw_1,zsc_vw_2 ! Temporary scales for non-gradient momentum flux terms. 236 315 REAL(wp), DIMENSION(jpi,jpj) :: zhbl_t ! holds boundary layer depth updated by full timestep … … 243 322 ! Temporary variables 244 323 INTEGER :: inhml 245 INTEGER :: i_lconv_alloc246 324 REAL(wp) :: znd,znd_d,zznd_ml,zznd_pyc,zznd_d ! temporary non-dimensional depths used in various routines 247 325 REAL(wp) :: ztemp, zari, zpert, zzdhdt, zdb ! temporary variables 248 326 REAL(wp) :: zthick, zz0, zz1 ! temporary variables 249 327 REAL(wp) :: zvel_max, zhbl_s ! temporary variables 250 REAL(wp) :: zfac 328 REAL(wp) :: zfac, ztmp ! temporary variable 251 329 REAL(wp) :: zus_x, zus_y ! temporary Stokes drift 252 330 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zviscos ! viscosity 253 331 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdiffut ! t-diffusivity 332 REAL(wp), DIMENSION(jpi,jpj) :: zalpha_pyc 333 REAL(wp), DIMENSION(jpi,jpj) :: ztau_sc_u ! dissipation timescale at baes of WML. 334 REAL(wp) :: zdelta_pyc, zwt_pyc_sc_1, zws_pyc_sc_1, zzeta_pyc 335 REAL(wp) :: zbuoy_pyc_sc, zomega, zvw_max 336 INTEGER :: ibld_ext=0 ! does not have to be zero for modified scheme 337 REAL(wp) :: zgamma_b_nd, zgamma_b, zdhoh, ztau 338 REAL(wp) :: zzeta_s = 0._wp 339 REAL(wp) :: zzeta_v = 0.46 340 REAL(wp) :: zabsstke 341 REAL(wp) :: zsqrtpi, z_two_thirds, zproportion, ztransp, zthickness 342 REAL(wp) :: z2k_times_thickness, zsqrt_depth, zexp_depth, zdstokes0, zf, zexperfc 254 343 255 344 ! For debugging … … 257 346 !!-------------------------------------------------------------------- 258 347 ! 259 ALLOCATE( lconv(jpi,jpj), STAT= i_lconv_alloc )260 IF( i_lconv_alloc /= 0 ) CALL ctl_warn('zdf_osm: failed to allocate lconv')261 262 348 ibld(:,:) = 0 ; imld(:,:) = 0 263 349 zrad0(:,:) = 0._wp ; zradh(:,:) = 0._wp ; zradav(:,:) = 0._wp ; zustar(:,:) = 0._wp … … 267 353 zustke(:,:) = 0._wp ; zla(:,:) = 0._wp ; zcos_wind(:,:) = 0._wp ; zsin_wind(:,:) = 0._wp 268 354 zhol(:,:) = 0._wp 269 lconv(:,:) = .FALSE. 355 lconv(:,:) = .FALSE.; lpyc(:,:) = .FALSE. ; lflux(:,:) = .FALSE. ; lmle(:,:) = .FALSE. 270 356 ! mixed layer 271 357 ! no initialization of zhbl or zhml (or zdh?) 272 358 zhbl(:,:) = 1._wp ; zhml(:,:) = 1._wp ; zdh(:,:) = 1._wp ; zdhdt(:,:) = 0._wp 273 zt_bl(:,:) = 0._wp ; zs_bl(:,:) = 0._wp ; zu_bl(:,:) = 0._wp ; zv_bl(:,:) = 0._wp 274 zrh_bl(:,:) = 0._wp ; zt_ml(:,:) = 0._wp ; zs_ml(:,:) = 0._wp ; zu_ml(:,:) = 0._wp 275 zv_ml(:,:) = 0._wp ; zrh_ml(:,:) = 0._wp ; zdt_bl(:,:) = 0._wp ; zds_bl(:,:) = 0._wp 276 zdu_bl(:,:) = 0._wp ; zdv_bl(:,:) = 0._wp ; zdrh_bl(:,:) = 0._wp ; zdb_bl(:,:) = 0._wp 359 zt_bl(:,:) = 0._wp ; zs_bl(:,:) = 0._wp ; zu_bl(:,:) = 0._wp 360 zv_bl(:,:) = 0._wp ; zb_bl(:,:) = 0._wp 361 zt_ml(:,:) = 0._wp ; zs_ml(:,:) = 0._wp ; zu_ml(:,:) = 0._wp 362 zt_mle(:,:) = 0._wp ; zs_mle(:,:) = 0._wp ; zu_mle(:,:) = 0._wp 363 zb_mle(:,:) = 0._wp 364 zv_ml(:,:) = 0._wp ; zdt_bl(:,:) = 0._wp ; zds_bl(:,:) = 0._wp 365 zdu_bl(:,:) = 0._wp ; zdv_bl(:,:) = 0._wp ; zdb_bl(:,:) = 0._wp 277 366 zdt_ml(:,:) = 0._wp ; zds_ml(:,:) = 0._wp ; zdu_ml(:,:) = 0._wp ; zdv_ml(:,:) = 0._wp 278 zdrh_ml(:,:) = 0._wp ; zdb_ml(:,:) = 0._wp ; zwth_ent(:,:) = 0._wp ; zws_ent(:,:) = 0._wp 279 zuw_bse(:,:) = 0._wp ; zvw_bse(:,:) = 0._wp 367 zdb_ml(:,:) = 0._wp 368 zdt_mle(:,:) = 0._wp ; zds_mle(:,:) = 0._wp ; zdu_mle(:,:) = 0._wp 369 zdv_mle(:,:) = 0._wp ; zdb_mle(:,:) = 0._wp 370 zwth_ent = 0._wp ; zws_ent = 0._wp 280 371 ! 281 372 zdtdz_pyc(:,:,:) = 0._wp ; zdsdz_pyc(:,:,:) = 0._wp ; zdbdz_pyc(:,:,:) = 0._wp 282 373 zdudz_pyc(:,:,:) = 0._wp ; zdvdz_pyc(:,:,:) = 0._wp 283 374 ! 375 zdtdz_bl_ext(:,:) = 0._wp ; zdsdz_bl_ext(:,:) = 0._wp ; zdbdz_bl_ext(:,:) = 0._wp 376 377 IF ( ln_osm_mle ) THEN ! only initialise arrays if needed 378 zdtdx(:,:) = 0._wp ; zdtdy(:,:) = 0._wp ; zdsdx(:,:) = 0._wp 379 zdsdy(:,:) = 0._wp ; dbdx_mle(:,:) = 0._wp ; dbdy_mle(:,:) = 0._wp 380 zwb_fk(:,:) = 0._wp ; zvel_mle(:,:) = 0._wp; zdiff_mle(:,:) = 0._wp 381 zhmle(:,:) = 0._wp ; zmld(:,:) = 0._wp 382 ENDIF 383 zwb_fk_b(:,:) = 0._wp ! must be initialised even with ln_osm_mle=F as used in zdf_osm_calculate_dhdt 384 284 385 ! Flux-Gradient arrays. 285 zdifml_sc(:,:) = 0._wp ; zvisml_sc(:,:) = 0._wp ; zdifpyc_sc(:,:) = 0._wp286 zvispyc_sc(:,:) = 0._wp ; zbeta_d_sc(:,:) = 0._wp ; zbeta_v_sc(:,:) = 0._wp287 386 zsc_wth_1(:,:) = 0._wp ; zsc_ws_1(:,:) = 0._wp ; zsc_uw_1(:,:) = 0._wp 288 387 zsc_uw_2(:,:) = 0._wp ; zsc_vw_1(:,:) = 0._wp ; zsc_vw_2(:,:) = 0._wp … … 292 391 ghams(:,:,:) = 0._wp ; ghamu(:,:,:) = 0._wp ; ghamv(:,:,:) = 0._wp 293 392 393 zddhdt(:,:) = 0._wp 294 394 ! hbl = MAX(hbl,epsln) 295 395 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 326 426 zwbav(ji,jj) = grav * zthermal * zwthav(ji,jj) - grav * zbeta * zwsav(ji,jj) 327 427 ! Surface upward velocity fluxes 328 zuw0(ji,jj) = - utau(ji,jj) * r1_rho0 * tmask(ji,jj,1)329 zvw0(ji,jj) = - vtau(ji,jj) * r1_rho0 * tmask(ji,jj,1)428 zuw0(ji,jj) = - 0.5 * (utau(ji-1,jj) + utau(ji,jj)) * r1_rho0 * tmask(ji,jj,1) 429 zvw0(ji,jj) = - 0.5 * (vtau(ji,jj-1) + vtau(ji,jj)) * r1_rho0 * tmask(ji,jj,1) 330 430 ! Friction velocity (zustar), at T-point : LMD94 eq. 2 331 431 zustar(ji,jj) = MAX( SQRT( SQRT( zuw0(ji,jj) * zuw0(ji,jj) + zvw0(ji,jj) * zvw0(ji,jj) ) ), 1.0e-8 ) … … 340 440 zus_x = zcos_wind(ji,jj) * zustar(ji,jj) / 0.3**2 341 441 zus_y = zsin_wind(ji,jj) * zustar(ji,jj) / 0.3**2 442 ! Linearly 342 443 zustke(ji,jj) = MAX ( SQRT( zus_x*zus_x + zus_y*zus_y), 1.0e-8 ) 343 ! dstokes(ji,jj) set to constant value rn_osm_dstokes from namelist in zdf_osm_init444 dstokes(ji,jj) = rn_osm_dstokes 344 445 END_2D 345 446 ! Assume Pierson-Moskovitz wind-wave spectrum … … 347 448 DO_2D( 0, 0, 0, 0 ) 348 449 ! Use wind speed wndm included in sbc_oce module 349 zustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8 )350 dstokes(ji,jj) = 0.12 * wndm(ji,jj)**2 / grav450 zustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8 ) 451 dstokes(ji,jj) = MAX ( 0.12 * wndm(ji,jj)**2 / grav, 5.e-1) 351 452 END_2D 352 453 ! Use ECMWF wave fields as output from SBCWAVE 353 454 CASE(2) 354 455 zfac = 2.0_wp * rpi / 16.0_wp 456 355 457 DO_2D( 0, 0, 0, 0 ) 356 ! The Langmur number from the ECMWF model appears to give La<0.3 for wind-driven seas. 357 ! The coefficient 0.8 gives La=0.3 in this situation. 358 ! It could represent the effects of the spread of wave directions 359 ! around the mean wind. The effect of this adjustment needs to be tested. 360 zustke(ji,jj) = MAX ( 1.0 * ( zcos_wind(ji,jj) * ut0sd(ji,jj ) + zsin_wind(ji,jj) * vt0sd(ji,jj) ), & 361 & zustar(ji,jj) / ( 0.45 * 0.45 ) ) 362 dstokes(ji,jj) = MAX(zfac * hsw(ji,jj)*hsw(ji,jj) / ( MAX(zustke(ji,jj)*wmp(ji,jj), 1.0e-7 ) ), 5.0e-1) !rn_osm_dstokes ! 458 IF (hsw(ji,jj) > 1.e-4) THEN 459 ! Use wave fields 460 zabsstke = SQRT(ut0sd(ji,jj)**2 + vt0sd(ji,jj)**2) 461 zustke(ji,jj) = MAX ( ( zcos_wind(ji,jj) * ut0sd(ji,jj) + zsin_wind(ji,jj) * vt0sd(ji,jj) ), 1.0e-8) 462 dstokes(ji,jj) = MAX (zfac * hsw(ji,jj)*hsw(ji,jj) / ( MAX(zabsstke * wmp(ji,jj), 1.0e-7 ) ), 5.0e-1) 463 ELSE 464 ! Assume masking issue (e.g. ice in ECMWF reanalysis but not in model run) 465 ! .. so default to Pierson-Moskowitz 466 zustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8 ) 467 dstokes(ji,jj) = MAX ( 0.12 * wndm(ji,jj)**2 / grav, 5.e-1) 468 END IF 469 END_2D 470 END SELECT 471 472 IF (ln_zdfosm_ice_shelter) THEN 473 ! Reduce both Stokes drift and its depth scale by ocean fraction to represent sheltering by ice 474 DO_2D( 0, 0, 0, 0 ) 475 zustke(ji,jj) = zustke(ji,jj) * (1.0_wp - fr_i(ji,jj)) 476 dstokes(ji,jj) = dstokes(ji,jj) * (1.0_wp - fr_i(ji,jj)) 477 END_2D 478 END IF 479 480 SELECT CASE (nn_osm_SD_reduce) 481 ! Reduce surface Stokes drift by a constant factor or following Breivik (2016) + van Roekel (2012) or Grant (2020). 482 CASE(0) 483 ! The Langmur number from the ECMWF model (or from PM) appears to give La<0.3 for wind-driven seas. 484 ! The coefficient rn_zdfosm_adjust_sd = 0.8 gives La=0.3 in this situation. 485 ! It could represent the effects of the spread of wave directions 486 ! around the mean wind. The effect of this adjustment needs to be tested. 487 IF(nn_osm_wave > 0) THEN 488 zustke(2:jpim1,2:jpjm1) = rn_zdfosm_adjust_sd * zustke(2:jpim1,2:jpjm1) 489 END IF 490 CASE(1) 491 ! van Roekel (2012): consider average SD over top 10% of boundary layer 492 ! assumes approximate depth profile of SD from Breivik (2016) 493 zsqrtpi = SQRT(rpi) 494 z_two_thirds = 2.0_wp / 3.0_wp 495 496 DO_2D( 0, 0, 0, 0 ) 497 zthickness = rn_osm_hblfrac*hbl(ji,jj) 498 z2k_times_thickness = zthickness * 2.0_wp / MAX( ABS( 5.97_wp * dstokes(ji,jj) ), 0.0000001_wp ) 499 zsqrt_depth = SQRT(z2k_times_thickness) 500 zexp_depth = EXP(-z2k_times_thickness) 501 zustke(ji,jj) = zustke(ji,jj) * (1.0_wp - zexp_depth & 502 & - z_two_thirds * ( zsqrtpi*zsqrt_depth*z2k_times_thickness * ERFC(zsqrt_depth) & 503 & + 1.0_wp - (1.0_wp + z2k_times_thickness)*zexp_depth ) ) / z2k_times_thickness 504 505 END_2D 506 CASE(2) 507 ! Grant (2020): Match to exponential with same SD and d/dz(Sd) at depth 10% of boundary layer 508 ! assumes approximate depth profile of SD from Breivik (2016) 509 zsqrtpi = SQRT(rpi) 510 511 DO_2D( 0, 0, 0, 0 ) 512 zthickness = rn_osm_hblfrac*hbl(ji,jj) 513 z2k_times_thickness = zthickness * 2.0_wp / MAX( ABS( 5.97_wp * dstokes(ji,jj) ), 0.0000001_wp ) 514 515 IF(z2k_times_thickness < 50._wp) THEN 516 zsqrt_depth = SQRT(z2k_times_thickness) 517 zexperfc = zsqrtpi * zsqrt_depth * ERFC(zsqrt_depth) * EXP(z2k_times_thickness) 518 ELSE 519 ! asymptotic expansion of sqrt(pi)*zsqrt_depth*EXP(z2k_times_thickness)*ERFC(zsqrt_depth) for large z2k_times_thickness 520 ! See Abramowitz and Stegun, Eq. 7.1.23 521 ! zexperfc = 1._wp - (1/2)/(z2k_times_thickness) + (3/4)/(z2k_times_thickness**2) - (15/8)/(z2k_times_thickness**3) 522 zexperfc = ((- 1.875_wp/z2k_times_thickness + 0.75_wp)/z2k_times_thickness - 0.5_wp)/z2k_times_thickness + 1.0_wp 523 END IF 524 zf = z2k_times_thickness*(1.0_wp/zexperfc - 1.0_wp) 525 dstokes(ji,jj) = 5.97 * zf * dstokes(ji,jj) 526 zustke(ji,jj) = zustke(ji,jj) * EXP(z2k_times_thickness * ( 1.0_wp / (2. * zf) - 1.0_wp )) * ( 1.0_wp - zexperfc) 363 527 END_2D 364 528 END SELECT … … 369 533 ! Langmuir velocity scale (zwstrl), at T-point 370 534 zwstrl(ji,jj) = ( zustar(ji,jj) * zustar(ji,jj) * zustke(ji,jj) )**pthird 371 ! Modify zwstrl to allow for small and large values of dstokes/hbl. 372 ! Intended as a possible test. Doesn't affect LES results for entrainment, 373 ! but hasn't been shown to be correct as dstokes/h becomes large or small. 374 zwstrl(ji,jj) = zwstrl(ji,jj) * & 375 & (1.12 * ( 1.0 - ( 1.0 - EXP( -hbl(ji,jj) / dstokes(ji,jj) ) ) * dstokes(ji,jj) / hbl(ji,jj) ))**pthird * & 376 & ( 1.0 - EXP( -15.0 * dstokes(ji,jj) / hbl(ji,jj) )) 377 ! define La this way so effects of Stokes penetration depth on velocity scale are included 378 zla(ji,jj) = SQRT ( zustar(ji,jj) / zwstrl(ji,jj) )**3 535 zla(ji,jj) = MAX(MIN(SQRT ( zustar(ji,jj) / ( zwstrl(ji,jj) + epsln ) )**3, 4.0), 0.2) 536 IF(zla(ji,jj) > 0.45) dstokes(ji,jj) = MIN(dstokes(ji,jj), 0.5_wp*hbl(ji,jj)) 379 537 ! Velocity scale that tends to zustar for large Langmuir numbers 380 538 zvstr(ji,jj) = ( zwstrl(ji,jj)**3 + & … … 383 541 ! limit maximum value of Langmuir number as approximate treatment for shear turbulence. 384 542 ! Note zustke and zwstrl are not amended. 385 IF ( zla(ji,jj) >= 0.45 ) zla(ji,jj) = 0.45386 543 ! 387 544 ! get convective velocity (zwstrc), stabilty scale (zhol) and logical conection flag lconv … … 389 546 zwstrc(ji,jj) = ( 2.0 * zwbav(ji,jj) * 0.9 * hbl(ji,jj) )**pthird 390 547 zhol(ji,jj) = -0.9 * hbl(ji,jj) * 2.0 * zwbav(ji,jj) / (zvstr(ji,jj)**3 + epsln ) 391 lconv(ji,jj) = .TRUE. 392 ELSE 548 ELSE 393 549 zhol(ji,jj) = -hbl(ji,jj) * 2.0 * zwbav(ji,jj)/ (zvstr(ji,jj)**3 + epsln ) 394 lconv(ji,jj) = .FALSE.395 550 ENDIF 396 551 END_2D … … 399 554 ! Mixed-layer model - calculate averages over the boundary layer, and the change in the boundary layer depth 400 555 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 401 ! BL must be always 2 levels deep. 402 hbl(:,:) = MAX(hbl(:,:), gdepw(:,:,3,Kmm) ) 403 ibld(:,:) = 3 404 DO_3D( 0, 0, 0, 0, 4, jpkm1 ) 556 ! BL must be always 4 levels deep. 557 ! For calculation of lateral buoyancy gradients for FK in 558 ! zdf_osm_zmld_horizontal_gradients need halo values for ibld, so must 559 ! previously exist for hbl also. 560 561 ! agn 23/6/20: not clear all this is needed, as hbl checked after it is re-calculated anyway 562 ! ########################################################################## 563 hbl(:,:) = MAX(hbl(:,:), gdepw(:,:,4,Kmm) ) 564 ibld(:,:) = 4 565 DO_3D( 1, 1, 1, 1, 5, jpkm1 ) 405 566 IF ( hbl(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN 406 567 ibld(ji,jj) = MIN(mbkt(ji,jj), jk) 407 568 ENDIF 408 569 END_3D 570 ! ########################################################################## 409 571 410 572 DO_2D( 0, 0, 0, 0 ) 411 zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 412 zbeta = rab_n(ji,jj,1,jp_sal) 413 zt = 0._wp 414 zs = 0._wp 415 zu = 0._wp 416 zv = 0._wp 417 ! average over depth of boundary layer 418 zthick=0._wp 419 DO jm = 2, ibld(ji,jj) 420 zthick=zthick+e3t(ji,jj,jm,Kmm) 421 zt = zt + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_tem,Kmm) 422 zs = zs + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_sal,Kmm) 423 zu = zu + e3t(ji,jj,jm,Kmm) & 424 & * ( uu(ji,jj,jm,Kbb) + uu(ji - 1,jj,jm,Kbb) ) & 425 & / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 426 zv = zv + e3t(ji,jj,jm,Kmm) & 427 & * ( vv(ji,jj,jm,Kbb) + vv(ji,jj - 1,jm,Kbb) ) & 428 & / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 429 END DO 430 zt_bl(ji,jj) = zt / zthick 431 zs_bl(ji,jj) = zs / zthick 432 zu_bl(ji,jj) = zu / zthick 433 zv_bl(ji,jj) = zv / zthick 434 zdt_bl(ji,jj) = zt_bl(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_tem,Kmm) 435 zds_bl(ji,jj) = zs_bl(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_sal,Kmm) 436 zdu_bl(ji,jj) = zu_bl(ji,jj) - ( uu(ji,jj,ibld(ji,jj),Kbb) + uu(ji-1,jj,ibld(ji,jj) ,Kbb) ) & 437 & / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 438 zdv_bl(ji,jj) = zv_bl(ji,jj) - ( vv(ji,jj,ibld(ji,jj),Kbb) + vv(ji,jj-1,ibld(ji,jj) ,Kbb) ) & 439 & / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 440 zdb_bl(ji,jj) = grav * zthermal * zdt_bl(ji,jj) - grav * zbeta * zds_bl(ji,jj) 441 IF ( lconv(ji,jj) ) THEN ! Convective 442 zwb_ent(ji,jj) = -( 2.0 * 0.2 * zwbav(ji,jj) & 443 & + 0.135 * zla(ji,jj) * zwstrl(ji,jj)**3/hbl(ji,jj) ) 444 445 zvel_max = - ( 1.0 + 1.0 * ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_Dt / hbl(ji,jj) ) & 446 & * zwb_ent(ji,jj) / ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 447 ! Entrainment including component due to shear turbulence. Modified Langmuir component, but gives same result for La=0.3 For testing uncomment. 448 ! zwb_ent(ji,jj) = -( 2.0 * 0.2 * zwbav(ji,jj) & 449 ! & + ( 0.15 * ( 1.0 - EXP( -0.5 * zla(ji,jj) ) ) + 0.03 / zla(ji,jj)**2 ) * zustar(ji,jj)**3/hbl(ji,jj) ) 450 451 ! zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_Dt / zhbl(ji,jj) ) * zwb_ent(ji,jj) / & 452 ! & ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 453 zzdhdt = - zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj),0.0) ) 454 ELSE ! Stable 455 zzdhdt = 0.32 * ( hbli(ji,jj) / hbl(ji,jj) -1.0 ) * zwstrl(ji,jj)**3 / hbli(ji,jj) & 456 & + ( ( 0.32 / 3.0 ) * exp ( -2.5 * ( hbli(ji,jj) / hbl(ji,jj) - 1.0 ) ) & 457 & - ( 0.32 / 3.0 - 0.135 * zla(ji,jj) ) * exp ( -12.5 * ( hbli(ji,jj) / hbl(ji,jj) ) ) ) & 458 & * zwstrl(ji,jj)**3 / hbli(ji,jj) 459 zzdhdt = zzdhdt + zwbav(ji,jj) 460 IF ( zzdhdt < 0._wp ) THEN 461 ! For long timsteps factor in brackets slows the rapid collapse of the OSBL 462 zpert = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_Dt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj) 463 ELSE 464 zpert = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_Dt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj) & 465 & + MAX( zdb_bl(ji,jj), 0.0 ) 466 ENDIF 467 zzdhdt = 2.0 * zzdhdt / zpert 468 ENDIF 469 zdhdt(ji,jj) = zzdhdt 573 zhbl(ji,jj) = gdepw(ji,jj,ibld(ji,jj),Kmm) 574 imld(ji,jj) = MAX(3,ibld(ji,jj) - MAX( INT( dh(ji,jj) / e3t(ji, jj, ibld(ji,jj), Kmm )) , 1 )) 575 zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm) 576 zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 470 577 END_2D 471 472 ! Calculate averages over depth of boundary layer 473 imld = ibld ! use imld to hold previous blayer index 474 ibld(:,:) = 3 475 476 zhbl_t(:,:) = hbl(:,:) + (zdhdt(:,:) - ww(ji,jj,ibld(ji,jj)))* rn_Dt ! certainly need wb here, so subtract it 477 zhbl_t(:,:) = MIN(zhbl_t(:,:), ht(:,:)) 478 zdhdt(:,:) = MIN(zdhdt(:,:), (zhbl_t(:,:) - hbl(:,:))/rn_Dt + ww(ji,jj,ibld(ji,jj))) ! adjustment to represent limiting by ocean bottom 578 ! Averages over well-mixed and boundary layer 579 jp_ext(:,:) = 2 580 CALL zdf_osm_vertical_average(ibld, jp_ext, zt_bl, zs_bl, zb_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl) 581 ! jp_ext(:,:) = ibld(:,:) - imld(:,:) + 1 582 CALL zdf_osm_vertical_average(ibld, ibld-imld+1, zt_ml, zs_ml, zb_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml) 583 ! Velocity components in frame aligned with surface stress. 584 CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_ml, zv_ml, zdu_ml, zdv_ml ) 585 CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_bl, zv_bl, zdu_bl, zdv_bl ) 586 ! Determine the state of the OSBL, stable/unstable, shear/no shear 587 CALL zdf_osm_osbl_state( lconv, lshear, j_ddh, zwb_ent, zwb_min, zshear, zri_i ) 588 589 IF ( ln_osm_mle ) THEN 590 ! Fox-Kemper Scheme 591 mld_prof = 4 592 DO_3D( 0, 0, 0, 0, 5, jpkm1 ) 593 IF ( hmle(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk) 594 END_3D 595 jp_ext_mle(:,:) = 2 596 CALL zdf_osm_vertical_average(mld_prof, jp_ext_mle, zt_mle, zs_mle, zb_mle, zu_mle, zv_mle, zdt_mle, zds_mle, zdb_mle, zdu_mle, zdv_mle) 597 598 DO_2D( 0, 0, 0, 0 ) 599 zhmle(ji,jj) = gdepw(ji,jj,mld_prof(ji,jj),Kmm) 600 END_2D 601 602 !! External gradient 603 CALL zdf_osm_external_gradients( ibld+2, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext ) 604 CALL zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle, zdbds_mle ) 605 CALL zdf_osm_external_gradients( mld_prof, zdtdz_mle_ext, zdsdz_mle_ext, zdbdz_mle_ext ) 606 CALL zdf_osm_osbl_state_fk( lpyc, lflux, lmle, zwb_fk ) 607 CALL zdf_osm_mle_parameters( mld_prof, hmle, zhmle, zvel_mle, zdiff_mle ) 608 ELSE ! ln_osm_mle 609 ! FK not selected, Boundary Layer only. 610 lpyc(:,:) = .TRUE. 611 lflux(:,:) = .FALSE. 612 lmle(:,:) = .FALSE. 613 DO_2D( 0, 0, 0, 0 ) 614 IF ( lconv(ji,jj) .AND. zdb_bl(ji,jj) < rn_osm_bl_thresh ) lpyc(ji,jj) = .FALSE. 615 END_2D 616 ENDIF ! ln_osm_mle 617 618 ! Test if pycnocline well resolved 619 DO_2D( 0, 0, 0, 0 ) 620 IF (lconv(ji,jj) ) THEN 621 ztmp = 0.2 * zhbl(ji,jj) / e3w(ji,jj,ibld(ji,jj),Kmm) 622 IF ( ztmp > 6 ) THEN 623 ! pycnocline well resolved 624 jp_ext(ji,jj) = 1 625 ELSE 626 ! pycnocline poorly resolved 627 jp_ext(ji,jj) = 0 628 ENDIF 629 ELSE 630 ! Stable conditions 631 jp_ext(ji,jj) = 0 632 ENDIF 633 END_2D 634 635 CALL zdf_osm_vertical_average(ibld, jp_ext, zt_bl, zs_bl, zb_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl ) 636 ! jp_ext = ibld-imld+1 637 CALL zdf_osm_vertical_average(imld-1, ibld-imld+1, zt_ml, zs_ml, zb_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml) 638 ! Rate of change of hbl 639 CALL zdf_osm_calculate_dhdt( zdhdt, zddhdt ) 640 DO_2D( 0, 0, 0, 0 ) 641 zhbl_t(ji,jj) = hbl(ji,jj) + (zdhdt(ji,jj) - ww(ji,jj,ibld(ji,jj)))* rn_Dt ! certainly need wn here, so subtract it 642 ! adjustment to represent limiting by ocean bottom 643 IF ( zhbl_t(ji,jj) >= gdepw(ji, jj, mbkt(ji,jj) + 1, Kmm ) ) THEN 644 zhbl_t(ji,jj) = MIN(zhbl_t(ji,jj), gdepw(ji,jj, mbkt(ji,jj) + 1, Kmm) - depth_tol)! ht(:,:)) 645 lpyc(ji,jj) = .FALSE. 646 ENDIF 647 END_2D 648 649 imld(:,:) = ibld(:,:) ! use imld to hold previous blayer index 650 ibld(:,:) = 4 479 651 480 652 DO_3D( 0, 0, 0, 0, 4, jpkm1 ) 481 653 IF ( zhbl_t(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN 482 ibld(ji,jj) = MIN(mbkt(ji,jj), jk)654 ibld(ji,jj) = jk 483 655 ENDIF 484 656 END_3D … … 487 659 ! Step through model levels taking account of buoyancy change to determine the effect on dhdt 488 660 ! 661 CALL zdf_osm_timestep_hbl( zdhdt ) 662 ! is external level in bounds? 663 664 CALL zdf_osm_vertical_average( ibld, jp_ext, zt_bl, zs_bl, zb_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl ) 665 ! 666 ! 667 ! Check to see if lpyc needs to be changed 668 669 CALL zdf_osm_pycnocline_thickness( dh, zdh ) 670 489 671 DO_2D( 0, 0, 0, 0 ) 490 IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN 672 IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh .or. ibld(ji,jj) + jp_ext(ji,jj) >= mbkt(ji,jj) .or. ibld(ji,jj)-imld(jj,ji) == 1 ) lpyc(ji,jj) = .FALSE. 673 END_2D 674 675 dstokes(:,:) = MIN ( dstokes(:,:), hbl(:,:)/3. ) ! Limit delta for shallow boundary layers for calculating flux-gradient terms. 491 676 ! 492 ! If boundary layer changes by more than one level, need to check for stable layers between initial and final depths. 493 ! 494 zhbl_s = hbl(ji,jj) 495 jm = imld(ji,jj) 496 zthermal = rab_n(ji,jj,1,jp_tem) 497 zbeta = rab_n(ji,jj,1,jp_sal) 498 IF ( lconv(ji,jj) ) THEN 499 !unstable 500 zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_Dt / hbl(ji,jj) ) & 501 & * zwb_ent(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 502 503 DO jk = imld(ji,jj), ibld(ji,jj) 504 zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - ts(ji,jj,jm,jp_tem,Kmm) ) & 505 & - zbeta * ( zs_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ), 0.0 ) + zvel_max 506 507 zhbl_s = zhbl_s + MIN( - zwb_ent(ji,jj) / zdb * rn_Dt / FLOAT(ibld(ji,jj)-imld(ji,jj) ), & 508 & e3w(ji,jj,jk,Kmm) ) 509 zhbl_s = MIN(zhbl_s, ht(ji,jj)) 510 511 IF ( zhbl_s >= gdepw(ji,jj,jm+1,Kmm) ) jm = jm + 1 512 END DO 513 hbl(ji,jj) = zhbl_s 514 ibld(ji,jj) = jm 515 hbli(ji,jj) = hbl(ji,jj) 516 ELSE 517 ! stable 518 DO jk = imld(ji,jj), ibld(ji,jj) 519 zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - ts(ji,jj,jm,jp_tem,Kmm) ) & 520 & - zbeta * ( zs_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ), 0.0 ) & 521 & + 2.0 * zwstrl(ji,jj)**2 / zhbl_s 522 523 zhbl_s = zhbl_s + ( & 524 & 0.32 * ( hbli(ji,jj) / zhbl_s -1.0 ) & 525 & * zwstrl(ji,jj)**3 / hbli(ji,jj) & 526 & + ( ( 0.32 / 3.0 ) * EXP( - 2.5 * ( hbli(ji,jj) / zhbl_s -1.0 ) ) & 527 & - ( 0.32 / 3.0 - 0.0485 ) * EXP( - 12.5 * ( hbli(ji,jj) / zhbl_s ) ) ) & 528 & * zwstrl(ji,jj)**3 / hbli(ji,jj) ) / zdb * e3w(ji,jj,jk,Kmm) / zdhdt(ji,jj) ! ALMG to investigate whether need to include ww here 529 530 zhbl_s = MIN(zhbl_s, ht(ji,jj)) 531 IF ( zhbl_s >= gdepw(ji,jj,jm,Kmm) ) jm = jm + 1 532 END DO 533 hbl(ji,jj) = MAX(zhbl_s, gdepw(ji,jj,3,Kmm) ) 534 ibld(ji,jj) = MAX(jm, 3 ) 535 IF ( hbl(ji,jj) > hbli(ji,jj) ) hbli(ji,jj) = hbl(ji,jj) 536 ENDIF ! IF ( lconv ) 537 ELSE 538 ! change zero or one model level. 539 hbl(ji,jj) = zhbl_t(ji,jj) 540 IF ( lconv(ji,jj) ) THEN 541 hbli(ji,jj) = hbl(ji,jj) 542 ELSE 543 hbl(ji,jj) = MAX(hbl(ji,jj), gdepw(ji,jj,3,Kmm) ) 544 IF ( hbl(ji,jj) > hbli(ji,jj) ) hbli(ji,jj) = hbl(ji,jj) 545 ENDIF 546 ENDIF 547 zhbl(ji,jj) = gdepw(ji,jj,ibld(ji,jj),Kmm) 548 END_2D 549 dstokes(:,:) = MIN ( dstokes(:,:), hbl(:,:)/3. ) ! Limit delta for shallow boundary layers for calculating flux-gradient terms. 550 551 ! Recalculate averages over boundary layer after depth updated 552 ! Consider later combining this into the loop above and looking for columns 553 ! where the index for base of the boundary layer have changed 554 DO_2D( 0, 0, 0, 0 ) 555 zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 556 zbeta = rab_n(ji,jj,1,jp_sal) 557 zt = 0._wp 558 zs = 0._wp 559 zu = 0._wp 560 zv = 0._wp 561 ! average over depth of boundary layer 562 zthick=0._wp 563 DO jm = 2, ibld(ji,jj) 564 zthick=zthick+e3t(ji,jj,jm,Kmm) 565 zt = zt + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_tem,Kmm) 566 zs = zs + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_sal,Kmm) 567 zu = zu + e3t(ji,jj,jm,Kmm) & 568 & * ( uu(ji,jj,jm,Kbb) + uu(ji - 1,jj,jm,Kbb) ) & 569 & / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 570 zv = zv + e3t(ji,jj,jm,Kmm) & 571 & * ( vv(ji,jj,jm,Kbb) + vv(ji,jj - 1,jm,Kbb) ) & 572 & / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 573 END DO 574 zt_bl(ji,jj) = zt / zthick 575 zs_bl(ji,jj) = zs / zthick 576 zu_bl(ji,jj) = zu / zthick 577 zv_bl(ji,jj) = zv / zthick 578 zdt_bl(ji,jj) = zt_bl(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_tem,Kmm) 579 zds_bl(ji,jj) = zs_bl(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_sal,Kmm) 580 zdu_bl(ji,jj) = zu_bl(ji,jj) - ( uu(ji,jj,ibld(ji,jj),Kbb) + uu(ji-1,jj,ibld(ji,jj) ,Kbb) ) & 581 & / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 582 zdv_bl(ji,jj) = zv_bl(ji,jj) - ( vv(ji,jj,ibld(ji,jj),Kbb) + vv(ji,jj-1,ibld(ji,jj) ,Kbb) ) & 583 & / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 584 zdb_bl(ji,jj) = grav * zthermal * zdt_bl(ji,jj) - grav * zbeta * zds_bl(ji,jj) 585 zhbl(ji,jj) = gdepw(ji,jj,ibld(ji,jj),Kmm) 586 IF ( lconv(ji,jj) ) THEN 587 IF ( zdb_bl(ji,jj) > 0._wp )THEN 588 IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN ! near neutral stability 589 zari = 4.5 * ( zvstr(ji,jj)**2 ) & 590 & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01 591 ELSE ! unstable 592 zari = 4.5 * ( zwstrc(ji,jj)**2 ) & 593 & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01 594 ENDIF 595 IF ( zari > 0.2 ) THEN ! This test checks for weakly stratified pycnocline 596 zari = 0.2 597 zwb_ent(ji,jj) = 0._wp 598 ENDIF 599 inhml = MAX( INT( zari * zhbl(ji,jj) & 600 & / e3t(ji,jj,ibld(ji,jj),Kmm) ), 1 ) 601 imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 1) 602 zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm) 603 zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 604 ELSE ! IF (zdb_bl) 605 imld(ji,jj) = ibld(ji,jj) - 1 606 zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm) 607 zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 608 ENDIF 609 ELSE ! IF (lconv) 610 IF ( zdhdt(ji,jj) >= 0.0 ) THEN ! probably shouldn't include wm here 611 ! boundary layer deepening 612 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 613 ! pycnocline thickness set by stratification - use same relationship as for neutral conditions. 614 zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) & 615 & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01 , 0.2 ) 616 inhml = MAX( INT( zari * zhbl(ji,jj) & 617 & / e3t(ji,jj,ibld(ji,jj),Kmm) ), 1 ) 618 imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 1) 619 zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm) 620 zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 621 ELSE 622 imld(ji,jj) = ibld(ji,jj) - 1 623 zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm) 624 zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 625 ENDIF ! IF (zdb_bl > 0.0) 626 ELSE ! IF(dhdt >= 0) 627 ! boundary layer collapsing. 628 imld(ji,jj) = ibld(ji,jj) 629 zhml(ji,jj) = zhbl(ji,jj) 630 zdh(ji,jj) = 0._wp 631 ENDIF ! IF (dhdt >= 0) 632 ENDIF ! IF (lconv) 633 END_2D 634 635 ! Average over the depth of the mixed layer in the convective boundary layer 636 ! Also calculate entrainment fluxes for temperature and salinity 637 DO_2D( 0, 0, 0, 0 ) 638 zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 639 zbeta = rab_n(ji,jj,1,jp_sal) 640 IF ( lconv(ji,jj) ) THEN 641 zt = 0._wp 642 zs = 0._wp 643 zu = 0._wp 644 zv = 0._wp 645 ! average over depth of boundary layer 646 zthick=0._wp 647 DO jm = 2, imld(ji,jj) 648 zthick=zthick+e3t(ji,jj,jm,Kmm) 649 zt = zt + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_tem,Kmm) 650 zs = zs + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_sal,Kmm) 651 zu = zu + e3t(ji,jj,jm,Kmm) & 652 & * ( uu(ji,jj,jm,Kbb) + uu(ji - 1,jj,jm,Kbb) ) & 653 & / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 654 zv = zv + e3t(ji,jj,jm,Kmm) & 655 & * ( vv(ji,jj,jm,Kbb) + vv(ji,jj - 1,jm,Kbb) ) & 656 & / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 657 END DO 658 zt_ml(ji,jj) = zt / zthick 659 zs_ml(ji,jj) = zs / zthick 660 zu_ml(ji,jj) = zu / zthick 661 zv_ml(ji,jj) = zv / zthick 662 zdt_ml(ji,jj) = zt_ml(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_tem,Kmm) 663 zds_ml(ji,jj) = zs_ml(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_sal,Kmm) 664 zdu_ml(ji,jj) = zu_ml(ji,jj) - ( uu(ji,jj,ibld(ji,jj),Kbb) + uu(ji-1,jj,ibld(ji,jj) ,Kbb) ) & 665 & / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 666 zdv_ml(ji,jj) = zv_ml(ji,jj) - ( vv(ji,jj,ibld(ji,jj),Kbb) + vv(ji,jj-1,ibld(ji,jj) ,Kbb) ) & 667 & / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 668 zdb_ml(ji,jj) = grav * zthermal * zdt_ml(ji,jj) - grav * zbeta * zds_ml(ji,jj) 669 ELSE 670 ! stable, if entraining calulate average below interface layer. 671 IF ( zdhdt(ji,jj) >= 0._wp ) THEN 672 zt = 0._wp 673 zs = 0._wp 674 zu = 0._wp 675 zv = 0._wp 676 ! average over depth of boundary layer 677 zthick=0._wp 678 DO jm = 2, imld(ji,jj) 679 zthick=zthick+e3t(ji,jj,jm,Kmm) 680 zt = zt + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_tem,Kmm) 681 zs = zs + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_sal,Kmm) 682 zu = zu + e3t(ji,jj,jm,Kmm) & 683 & * ( uu(ji,jj,jm,Kbb) + uu(ji - 1,jj,jm,Kbb) ) & 684 & / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 685 zv = zv + e3t(ji,jj,jm,Kmm) & 686 & * ( vv(ji,jj,jm,Kbb) + vv(ji,jj - 1,jm,Kbb) ) & 687 & / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 688 END DO 689 zt_ml(ji,jj) = zt / zthick 690 zs_ml(ji,jj) = zs / zthick 691 zu_ml(ji,jj) = zu / zthick 692 zv_ml(ji,jj) = zv / zthick 693 zdt_ml(ji,jj) = zt_ml(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_tem,Kmm) 694 zds_ml(ji,jj) = zs_ml(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_sal,Kmm) 695 zdu_ml(ji,jj) = zu_ml(ji,jj) - ( uu(ji,jj,ibld(ji,jj),Kbb) + uu(ji-1,jj,ibld(ji,jj) ,Kbb) ) & 696 & / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 697 zdv_ml(ji,jj) = zv_ml(ji,jj) - ( vv(ji,jj,ibld(ji,jj),Kbb) + vv(ji,jj-1,ibld(ji,jj) ,Kbb) ) & 698 & / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 699 zdb_ml(ji,jj) = grav * zthermal * zdt_ml(ji,jj) - grav * zbeta * zds_ml(ji,jj) 700 ENDIF 701 ENDIF 702 END_2D 703 ! 677 ! Average over the depth of the mixed layer in the convective boundary layer 678 ! jp_ext = ibld - imld +1 679 CALL zdf_osm_vertical_average( imld-1, ibld-imld+1, zt_ml, zs_ml, zb_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml ) 704 680 ! rotate mean currents and changes onto wind align co-ordinates 705 681 ! 706 707 DO_2D( 0, 0, 0, 0 ) 708 ztemp = zu_ml(ji,jj) 709 zu_ml(ji,jj) = zu_ml(ji,jj) * zcos_wind(ji,jj) + zv_ml(ji,jj) * zsin_wind(ji,jj) 710 zv_ml(ji,jj) = zv_ml(ji,jj) * zcos_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 711 ztemp = zdu_ml(ji,jj) 712 zdu_ml(ji,jj) = zdu_ml(ji,jj) * zcos_wind(ji,jj) + zdv_ml(ji,jj) * zsin_wind(ji,jj) 713 zdv_ml(ji,jj) = zdv_ml(ji,jj) * zsin_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 714 ! 715 ztemp = zu_bl(ji,jj) 716 zu_bl = zu_bl(ji,jj) * zcos_wind(ji,jj) + zv_bl(ji,jj) * zsin_wind(ji,jj) 717 zv_bl(ji,jj) = zv_bl(ji,jj) * zcos_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 718 ztemp = zdu_bl(ji,jj) 719 zdu_bl(ji,jj) = zdu_bl(ji,jj) * zcos_wind(ji,jj) + zdv_bl(ji,jj) * zsin_wind(ji,jj) 720 zdv_bl(ji,jj) = zdv_bl(ji,jj) * zsin_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 721 END_2D 722 723 zuw_bse = 0._wp 724 zvw_bse = 0._wp 725 DO_2D( 0, 0, 0, 0 ) 726 727 IF ( lconv(ji,jj) ) THEN 728 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 729 zwth_ent(ji,jj) = zwb_ent(ji,jj) * zdt_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 730 zws_ent(ji,jj) = zwb_ent(ji,jj) * zds_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 731 ENDIF 732 ELSE 733 zwth_ent(ji,jj) = -2.0 * zwthav(ji,jj) * ( (1.0 - 0.8) - ( 1.0 - 0.8)**(3.0/2.0) ) 734 zws_ent(ji,jj) = -2.0 * zwsav(ji,jj) * ( (1.0 - 0.8 ) - ( 1.0 - 0.8 )**(3.0/2.0) ) 735 ENDIF 736 END_2D 737 682 CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_ml, zv_ml, zdu_ml, zdv_ml ) 683 CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_bl, zv_bl, zdu_bl, zdv_bl ) 738 684 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 739 685 ! Pycnocline gradients for scalars and velocity 740 686 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 741 687 742 DO_2D( 0, 0, 0, 0 ) 743 ! 744 IF ( lconv (ji,jj) ) THEN 745 ! Unstable conditions 746 IF( zdb_bl(ji,jj) > 0._wp ) THEN 747 ! calculate pycnocline profiles, no need if zdb_bl <= 0. since profile is zero and arrays have been initialized to zero 748 ztgrad = ( zdt_ml(ji,jj) / zdh(ji,jj) ) 749 zsgrad = ( zds_ml(ji,jj) / zdh(ji,jj) ) 750 zbgrad = ( zdb_ml(ji,jj) / zdh(ji,jj) ) 751 DO jk = 2 , ibld(ji,jj) 752 znd = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj) 753 zdtdz_pyc(ji,jj,jk) = ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 754 zdbdz_pyc(ji,jj,jk) = zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 755 zdsdz_pyc(ji,jj,jk) = zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 756 END DO 757 ENDIF 758 ELSE 759 ! stable conditions 760 ! if pycnocline profile only defined when depth steady of increasing. 761 IF ( zdhdt(ji,jj) >= 0.0 ) THEN ! Depth increasing, or steady. 762 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 763 IF ( zhol(ji,jj) >= 0.5 ) THEN ! Very stable - 'thick' pycnocline 764 ztgrad = zdt_bl(ji,jj) / zhbl(ji,jj) 765 zsgrad = zds_bl(ji,jj) / zhbl(ji,jj) 766 zbgrad = zdb_bl(ji,jj) / zhbl(ji,jj) 767 DO jk = 2, ibld(ji,jj) 768 znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 769 zdtdz_pyc(ji,jj,jk) = ztgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 770 zdbdz_pyc(ji,jj,jk) = zbgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 771 zdsdz_pyc(ji,jj,jk) = zsgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 772 END DO 773 ELSE ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form. 774 ztgrad = zdt_bl(ji,jj) / zdh(ji,jj) 775 zsgrad = zds_bl(ji,jj) / zdh(ji,jj) 776 zbgrad = zdb_bl(ji,jj) / zdh(ji,jj) 777 DO jk = 2, ibld(ji,jj) 778 znd = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj) 779 zdtdz_pyc(ji,jj,jk) = ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 780 zdbdz_pyc(ji,jj,jk) = zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 781 zdsdz_pyc(ji,jj,jk) = zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 782 END DO 783 ENDIF ! IF (zhol >=0.5) 784 ENDIF ! IF (zdb_bl> 0.) 785 ENDIF ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero, profile arrays are intialized to zero 786 ENDIF ! IF (lconv) 787 ! 788 END_2D 789 ! 790 DO_2D( 0, 0, 0, 0 ) 791 ! 792 IF ( lconv (ji,jj) ) THEN 793 ! Unstable conditions 794 zugrad = ( zdu_ml(ji,jj) / zdh(ji,jj) ) + 0.275 * zustar(ji,jj)*zustar(ji,jj) / & 795 & (( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) / zla(ji,jj)**(8.0/3.0) 796 zvgrad = ( zdv_ml(ji,jj) / zdh(ji,jj) ) + 3.5 * ff_t(ji,jj) * zustke(ji,jj) / & 797 & ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 798 DO jk = 2 , ibld(ji,jj)-1 799 znd = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj) 800 zdudz_pyc(ji,jj,jk) = zugrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 801 zdvdz_pyc(ji,jj,jk) = zvgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 802 END DO 803 ELSE 804 ! stable conditions 805 zugrad = 3.25 * zdu_bl(ji,jj) / zhbl(ji,jj) 806 zvgrad = 2.75 * zdv_bl(ji,jj) / zhbl(ji,jj) 807 DO jk = 2, ibld(ji,jj) 808 znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 809 IF ( znd < 1.0 ) THEN 810 zdudz_pyc(ji,jj,jk) = zugrad * EXP( -40.0 * ( znd - 1.0 )**2 ) 811 ELSE 812 zdudz_pyc(ji,jj,jk) = zugrad * EXP( -20.0 * ( znd - 1.0 )**2 ) 813 ENDIF 814 zdvdz_pyc(ji,jj,jk) = zvgrad * EXP( -20.0 * ( znd - 0.85 )**2 ) 815 END DO 816 ENDIF 817 ! 818 END_2D 688 CALL zdf_osm_external_gradients( ibld+2, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext ) 689 CALL zdf_osm_pycnocline_scalar_profiles( zdtdz_pyc, zdsdz_pyc, zdbdz_pyc, zalpha_pyc ) 690 CALL zdf_osm_pycnocline_shear_profiles( zdudz_pyc, zdvdz_pyc ) 819 691 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 820 692 ! Eddy viscosity/diffusivity and non-gradient terms in the flux-gradient relationship 821 693 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 822 823 ! WHERE ( lconv ) 824 ! zdifml_sc = zhml * ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird 825 ! zvisml_sc = zdifml_sc 826 ! zdifpyc_sc = 0.165 * ( zwstrl**3 + zwstrc**3 )**pthird * ( zhbl - zhml ) 827 ! zvispyc_sc = 0.142 * ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * ( zhbl - zhml ) 828 ! zbeta_d_sc = 1.0 - (0.165 / 0.8 * ( zhbl - zhml ) / zhbl )**p2third 829 ! zbeta_v_sc = 1.0 - 2.0 * (0.142 /0.375) * (zhbl - zhml ) / zhml 830 ! ELSEWHERE 831 ! zdifml_sc = zwstrl * zhbl * EXP ( -( zhol / 0.183_wp )**2 ) 832 ! zvisml_sc = zwstrl * zhbl * EXP ( -( zhol / 0.183_wp )**2 ) 833 ! ENDWHERE 834 DO_2D( 0, 0, 0, 0 ) 835 IF ( lconv(ji,jj) ) THEN 836 zdifml_sc(ji,jj) = zhml(ji,jj) * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 837 zvisml_sc(ji,jj) = zdifml_sc(ji,jj) 838 zdifpyc_sc(ji,jj) = 0.165 * ( zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3 )**pthird * zdh(ji,jj) 839 zvispyc_sc(ji,jj) = 0.142 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zdh(ji,jj) 840 zbeta_d_sc(ji,jj) = 1.0 - (0.165 / 0.8 * zdh(ji,jj) / zhbl(ji,jj) )**p2third 841 zbeta_v_sc(ji,jj) = 1.0 - 2.0 * (0.142 /0.375) * zdh(ji,jj) / zhml(ji,jj) 842 ELSE 843 zdifml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ) 844 zvisml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ) 845 END IF 846 END_2D 847 ! 848 DO_2D( 0, 0, 0, 0 ) 849 IF ( lconv(ji,jj) ) THEN 850 DO jk = 2, imld(ji,jj) ! mixed layer diffusivity 851 zznd_ml = gdepw(ji,jj,jk,Kmm) / zhml(ji,jj) 852 ! 853 zdiffut(ji,jj,jk) = 0.8 * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_d_sc(ji,jj) * zznd_ml )**1.5 854 ! 855 zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_v_sc(ji,jj) * zznd_ml ) & 856 & * ( 1.0 - 0.5 * zznd_ml**2 ) 857 END DO 858 ! pycnocline - if present linear profile 859 IF ( zdh(ji,jj) > 0._wp ) THEN 860 DO jk = imld(ji,jj)+1 , ibld(ji,jj) 861 zznd_pyc = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj) 862 ! 863 zdiffut(ji,jj,jk) = zdifpyc_sc(ji,jj) * ( 1.0 + zznd_pyc ) 864 ! 865 zviscos(ji,jj,jk) = zvispyc_sc(ji,jj) * ( 1.0 + zznd_pyc ) 866 END DO 867 ENDIF 868 ! Temporay fix to ensure zdiffut is +ve; won't be necessary with ww taken out 869 zdiffut(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj)* e3t(ji,jj,ibld(ji,jj),Kmm) 870 ! could be taken out, take account of entrainment represents as a diffusivity 871 ! should remove w from here, represents entrainment 872 ELSE 873 ! stable conditions 874 DO jk = 2, ibld(ji,jj) 875 zznd_ml = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 876 zdiffut(ji,jj,jk) = 0.75 * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zznd_ml )**1.5 877 zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * (1.0 - zznd_ml) * ( 1.0 - zznd_ml**2 ) 878 END DO 879 ENDIF ! end if ( lconv ) 880 ! 881 END_2D 694 CALL zdf_osm_diffusivity_viscosity( zdiffut, zviscos ) 882 695 883 696 ! … … 918 731 END_2D 919 732 920 921 733 ! 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) 922 734 WHERE ( lconv ) 923 zsc_uw_1 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke / ( 1.0 - 1.0 * 6.5 * zla**(8.0/3.0))924 zsc_uw_2 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke / ( zla**(8.0/3.0) + epsln)925 zsc_vw_1 = ff_t * zhml * zustke**3 * zla**(8.0/3.0) / ( ( zvstr**3 + 0.5 * zwstrc**3 )**(2.0/3.0) + epsln )735 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 ) 736 zsc_uw_2 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke / MIN( zla**(8.0/3.0) + epsln, 0.12 ) 737 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 ) 926 738 ELSEWHERE 927 739 zsc_uw_1 = zustar**2 928 zsc_vw_1 = ff_t * zhbl * zustke**3 * zla**(8.0/3.0) / (zvstr**2 + epsln)740 zsc_vw_1 = ff_t * zhbl * zustke**3 * MIN( zla**(8.0/3.0), 0.12 ) / (zvstr**2 + epsln) 929 741 ENDWHERE 930 742 IF(ln_dia_osm) THEN 743 IF ( iom_use("ghamu_00") ) CALL iom_put( "ghamu_00", wmask*ghamu ) 744 IF ( iom_use("ghamv_00") ) CALL iom_put( "ghamv_00", wmask*ghamv ) 745 END IF 931 746 DO_2D( 0, 0, 0, 0 ) 932 747 IF ( lconv(ji,jj) ) THEN … … 970 785 zl_l = 2.0 * ( 1.0 - EXP ( - 2.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) ) & 971 786 & * ( 1.0 - EXP ( - 5.0 * ( 1.0 - zznd_ml ) ) ) * ( 1.0 + dstokes(ji,jj) / zhml (ji,jj) ) 972 zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0 + EXP ( 3.0 * LOG10 ( - zhol(ji,jj) ) ) ) ** (3.0/2.0)787 zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0 + EXP ( -3.0 * LOG10 ( - zhol(ji,jj) ) ) ) ** (3.0 / 2.0) 973 788 ! non-gradient buoyancy terms 974 789 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * 0.5 * zsc_wth_1(ji,jj) * zl_eps * zhml(ji,jj) / ( 0.15 + zznd_ml ) 975 790 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * 0.5 * zsc_ws_1(ji,jj) * zl_eps * zhml(ji,jj) / ( 0.15 + zznd_ml ) 976 791 END DO 977 ELSE 792 793 IF ( lpyc(ji,jj) ) THEN 794 ztau_sc_u(ji,jj) = zhml(ji,jj) / ( zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird 795 ztau_sc_u(ji,jj) = ztau_sc_u(ji,jj) * ( 1.4 -0.4 / ( 1.0 + EXP( -3.5 * LOG10( -zhol(ji,jj) ) ) )**1.5 ) 796 zwth_ent = -0.003 * ( 0.15 * zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird * ( 1.0 - zdh(ji,jj) /zhbl(ji,jj) ) * zdt_ml(ji,jj) 797 zws_ent = -0.003 * ( 0.15 * zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird * ( 1.0 - zdh(ji,jj) /zhbl(ji,jj) ) * zds_ml(ji,jj) 798 ! Cubic profile used for buoyancy term 799 za_cubic = 0.755 * ztau_sc_u(ji,jj) 800 zb_cubic = 0.25 * ztau_sc_u(ji,jj) 801 DO jk = 2, ibld(ji,jj) 802 zznd_pyc = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / zdh(ji,jj) 803 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) - 0.045 * ( ( zwth_ent * zdbdz_pyc(ji,jj,jk) ) * ztau_sc_u(ji,jj)**2 ) * MAX( ( 1.75 * zznd_pyc -0.15 * zznd_pyc**2 - 0.2 * zznd_pyc**3 ), 0.0 ) 804 805 ghams(ji,jj,jk) = ghams(ji,jj,jk) - 0.045 * ( ( zws_ent * zdbdz_pyc(ji,jj,jk) ) * ztau_sc_u(ji,jj)**2 ) * MAX( ( 1.75 * zznd_pyc -0.15 * zznd_pyc**2 - 0.2 * zznd_pyc**3 ), 0.0 ) 806 END DO 807 ! 808 zbuoy_pyc_sc = zalpha_pyc(ji,jj) * zdb_ml(ji,jj) / zdh(ji,jj) + zdbdz_bl_ext(ji,jj) 809 zdelta_pyc = ( zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird / SQRT( MAX( zbuoy_pyc_sc, ( zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**p2third / zdh(ji,jj)**2 ) ) 810 ! 811 zwt_pyc_sc_1 = 0.325 * ( zalpha_pyc(ji,jj) * zdt_ml(ji,jj) / zdh(ji,jj) + zdtdz_bl_ext(ji,jj) ) * zdelta_pyc**2 / zdh(ji,jj) 812 ! 813 zws_pyc_sc_1 = 0.325 * ( zalpha_pyc(ji,jj) * zds_ml(ji,jj) / zdh(ji,jj) + zdsdz_bl_ext(ji,jj) ) * zdelta_pyc**2 / zdh(ji,jj) 814 ! 815 zzeta_pyc = 0.15 - 0.175 / ( 1.0 + EXP( -3.5 * LOG10( -zhol(ji,jj) ) ) ) 816 DO jk = 2, ibld(ji,jj) 817 zznd_pyc = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / zdh(ji,jj) 818 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.05 * zwt_pyc_sc_1 * EXP( -0.25 * ( zznd_pyc / zzeta_pyc )**2 ) * zdh(ji,jj) / ( zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird 819 ! 820 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.05 * zws_pyc_sc_1 * EXP( -0.25 * ( zznd_pyc / zzeta_pyc )**2 ) * zdh(ji,jj) / ( zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird 821 END DO 822 ENDIF ! End of pycnocline 823 ELSE ! lconv test - stable conditions 978 824 DO jk = 2, ibld(ji,jj) 979 825 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zsc_wth_1(ji,jj) … … 982 828 ENDIF 983 829 END_2D 984 985 830 986 831 WHERE ( lconv ) … … 1011 856 END_2D 1012 857 858 DO_2D( 0, 0, 0, 0 ) 859 IF ( lpyc(ji,jj) ) THEN 860 IF ( j_ddh(ji,jj) == 0 ) THEN 861 ! Place holding code. Parametrization needs checking for these conditions. 862 zomega = ( 0.15 * zwstrl(ji,jj)**3 + zwstrc(ji,jj)**3 + 4.75 * ( zshear(ji,jj)* zhbl(ji,jj) )**pthird )**pthird 863 zuw_bse = -0.0035 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdu_ml(ji,jj) 864 zvw_bse = -0.0075 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdv_ml(ji,jj) 865 ELSE 866 zomega = ( 0.15 * zwstrl(ji,jj)**3 + zwstrc(ji,jj)**3 + 4.75 * ( zshear(ji,jj)* zhbl(ji,jj) )**pthird )**pthird 867 zuw_bse = -0.0035 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdu_ml(ji,jj) 868 zvw_bse = -0.0075 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdv_ml(ji,jj) 869 ENDIF 870 zd_cubic = zdh(ji,jj) / zhbl(ji,jj) * zuw0(ji,jj) - ( 2.0 + zdh(ji,jj) /zhml(ji,jj) ) * zuw_bse 871 zc_cubic = zuw_bse - zd_cubic 872 ! need ztau_sc_u to be available. Change to array. 873 DO jk = imld(ji,jj), ibld(ji,jj) 874 zznd_pyc = - ( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / zdh(ji,jj) 875 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.045 * ztau_sc_u(ji,jj)**2 * ( zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 ) * ( 0.75 + 0.25 * zznd_pyc )**2 * zdbdz_pyc(ji,jj,jk) 876 END DO 877 zvw_max = 0.7 * ff_t(ji,jj) * ( zustke(ji,jj) * dstokes(ji,jj) + 0.75 * zustar(ji,jj) * zhml(ji,jj) ) 878 zd_cubic = zvw_max * zdh(ji,jj) / zhml(ji,jj) - ( 2.0 + zdh(ji,jj) /zhml(ji,jj) ) * zvw_bse 879 zc_cubic = zvw_bse - zd_cubic 880 DO jk = imld(ji,jj), ibld(ji,jj) 881 zznd_pyc = -( gdepw(ji,jj,jk,Kmm) -zhbl(ji,jj) ) / zdh(ji,jj) 882 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.045 * ztau_sc_u(ji,jj)**2 * ( zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 ) * ( 0.75 + 0.25 * zznd_pyc )**2 * zdbdz_pyc(ji,jj,jk) 883 END DO 884 ENDIF ! lpyc 885 END_2D 886 887 IF(ln_dia_osm) THEN 888 IF ( iom_use("ghamu_0") ) CALL iom_put( "ghamu_0", wmask*ghamu ) 889 IF ( iom_use("zsc_uw_1_0") ) CALL iom_put( "zsc_uw_1_0", tmask(:,:,1)*zsc_uw_1 ) 890 END IF 1013 891 ! Transport term in flux-gradient relationship [note : includes ROI ratio (X0.3) ] 1014 892 1015 WHERE ( lconv ) 1016 zsc_wth_1 = zwth0 1017 zsc_ws_1 = zws0 1018 ELSEWHERE 1019 zsc_wth_1 = 2.0 * zwthav 1020 zsc_ws_1 = zws0 1021 ENDWHERE 893 DO_2D( 1, 0, 1, 0 ) 894 895 IF ( lconv(ji,jj) ) THEN 896 zsc_wth_1(ji,jj) = zwth0(ji,jj) / ( 1.0 - 0.56 * EXP( zhol(ji,jj) ) ) 897 zsc_ws_1(ji,jj) = zws0(ji,jj) / (1.0 - 0.56 *EXP( zhol(ji,jj) ) ) 898 IF ( lpyc(ji,jj) ) THEN 899 ! Pycnocline scales 900 zsc_wth_pyc(ji,jj) = -0.2 * zwb0(ji,jj) * zdt_bl(ji,jj) / zdb_bl(ji,jj) 901 zsc_ws_pyc(ji,jj) = -0.2 * zwb0(ji,jj) * zds_bl(ji,jj) / zdb_bl(ji,jj) 902 ENDIF 903 ELSE 904 zsc_wth_1(ji,jj) = 2.0 * zwthav(ji,jj) 905 zsc_ws_1(ji,jj) = zws0(ji,jj) 906 ENDIF 907 END_2D 1022 908 1023 909 DO_2D( 0, 0, 0, 0 ) … … 1035 921 & * ( 1.0 - EXP ( -15.0 * ( 1.0 - zznd_ml ) ) ) 1036 922 END DO 923 ! 924 IF ( lpyc(ji,jj) ) THEN 925 ! pycnocline 926 DO jk = imld(ji,jj), ibld(ji,jj) 927 zznd_pyc = - ( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / zdh(ji,jj) 928 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 4.0 * zsc_wth_pyc(ji,jj) * ( 0.48 - EXP( -1.5 * ( zznd_pyc -0.3)**2 ) ) 929 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 4.0 * zsc_ws_pyc(ji,jj) * ( 0.48 - EXP( -1.5 * ( zznd_pyc -0.3)**2 ) ) 930 END DO 931 ENDIF 1037 932 ELSE 1038 DO jk = 2, ibld(ji,jj) 1039 zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 1040 znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 1041 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + & 1042 & 7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_wth_1(ji,jj) 1043 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + & 1044 & 7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_ws_1(ji,jj) 1045 END DO 933 IF( zdhdt(ji,jj) > 0. ) THEN 934 DO jk = 2, ibld(ji,jj) 935 zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 936 znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 937 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + & 938 & 7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_wth_1(ji,jj) 939 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + & 940 & 7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_ws_1(ji,jj) 941 END DO 942 ENDIF 1046 943 ENDIF 1047 944 END_2D 1048 1049 945 1050 946 WHERE ( lconv ) … … 1090 986 ENDIF 1091 987 END_2D 988 989 IF(ln_dia_osm) THEN 990 IF ( iom_use("ghamu_f") ) CALL iom_put( "ghamu_f", wmask*ghamu ) 991 IF ( iom_use("ghamv_f") ) CALL iom_put( "ghamv_f", wmask*ghamv ) 992 IF ( iom_use("zsc_uw_1_f") ) CALL iom_put( "zsc_uw_1_f", tmask(:,:,1)*zsc_uw_1 ) 993 IF ( iom_use("zsc_vw_1_f") ) CALL iom_put( "zsc_vw_1_f", tmask(:,:,1)*zsc_vw_1 ) 994 IF ( iom_use("zsc_uw_2_f") ) CALL iom_put( "zsc_uw_2_f", tmask(:,:,1)*zsc_uw_2 ) 995 IF ( iom_use("zsc_vw_2_f") ) CALL iom_put( "zsc_vw_2_f", tmask(:,:,1)*zsc_vw_2 ) 996 END IF 1092 997 ! 1093 998 ! Make surface forced velocity non-gradient terms go to zero at the base of the mixed layer. 1094 999 1000 1001 ! Make surface forced velocity non-gradient terms go to zero at the base of the boundary layer. 1002 1095 1003 DO_2D( 0, 0, 0, 0 ) 1096 IF ( lconv(ji,jj) ) THEN1004 IF ( .not. lconv(ji,jj) ) THEN 1097 1005 DO jk = 2, ibld(ji,jj) 1098 znd = ( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zhml(ji,jj) !ALMG to think about 1099 IF ( znd >= 0.0 ) THEN 1100 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -30.0 * znd**2 ) ) 1101 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * ( 1.0 - EXP( -30.0 * znd**2 ) ) 1102 ELSE 1103 ghamu(ji,jj,jk) = 0._wp 1104 ghamv(ji,jj,jk) = 0._wp 1105 ENDIF 1106 END DO 1107 ELSE 1108 DO jk = 2, ibld(ji,jj) 1109 znd = ( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zhml(ji,jj) !ALMG to think about 1006 znd = ( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / zhbl(ji,jj) !ALMG to think about 1110 1007 IF ( znd >= 0.0 ) THEN 1111 1008 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) ) … … 1120 1017 1121 1018 ! pynocline contributions 1122 ! Temporary fix to avoid instabilities when zdb_bl becomes very very small1123 zsc_uw_1 = 0._wp ! 50.0 * zla**(8.0/3.0) * zustar**2 * zhbl / ( zdb_bl + epsln )1124 1019 DO_2D( 0, 0, 0, 0 ) 1125 DO jk= 2, ibld(ji,jj) 1126 znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 1127 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk) 1128 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk) 1129 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc(ji,jj,jk) 1130 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zsc_uw_1(ji,jj) * ( 1.0 - znd )**(7.0/4.0) * zdbdz_pyc(ji,jj,jk) 1131 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc(ji,jj,jk) 1132 END DO 1020 IF ( .not. lconv(ji,jj) ) THEN 1021 IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN 1022 DO jk= 2, ibld(ji,jj) 1023 znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 1024 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk) 1025 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk) 1026 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc(ji,jj,jk) 1027 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc(ji,jj,jk) 1028 END DO 1029 END IF 1030 END IF 1133 1031 END_2D 1134 1135 ! Entrainment contribution. 1032 IF(ln_dia_osm) THEN 1033 IF ( iom_use("ghamu_b") ) CALL iom_put( "ghamu_b", wmask*ghamu ) 1034 IF ( iom_use("ghamv_b") ) CALL iom_put( "ghamv_b", wmask*ghamv ) 1035 END IF 1136 1036 1137 1037 DO_2D( 0, 0, 0, 0 ) 1138 IF ( lconv(ji,jj) ) THEN 1139 DO jk = 1, imld(ji,jj) - 1 1140 znd=gdepw(ji,jj,jk,Kmm) / zhml(ji,jj) 1141 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * znd 1142 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * znd 1143 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * znd 1144 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * znd 1145 END DO 1146 DO jk = imld(ji,jj), ibld(ji,jj) 1147 znd = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj) 1148 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * ( 1.0 + znd ) 1149 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * ( 1.0 + znd ) 1150 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * ( 1.0 + znd ) 1151 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * ( 1.0 + znd ) 1152 END DO 1153 ENDIF 1154 ghamt(ji,jj,ibld(ji,jj)) = 0._wp 1155 ghams(ji,jj,ibld(ji,jj)) = 0._wp 1156 ghamu(ji,jj,ibld(ji,jj)) = 0._wp 1157 ghamv(ji,jj,ibld(ji,jj)) = 0._wp 1038 ghamt(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 1039 ghams(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 1040 ghamu(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 1041 ghamv(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 1158 1042 END_2D 1159 1043 1160 1044 IF(ln_dia_osm) THEN 1045 IF ( iom_use("ghamu_1") ) CALL iom_put( "ghamu_1", wmask*ghamu ) 1046 IF ( iom_use("ghamv_1") ) CALL iom_put( "ghamv_1", wmask*ghamv ) 1047 IF ( iom_use("zdudz_pyc") ) CALL iom_put( "zdudz_pyc", wmask*zdudz_pyc ) 1048 IF ( iom_use("zdvdz_pyc") ) CALL iom_put( "zdvdz_pyc", wmask*zdvdz_pyc ) 1049 IF ( iom_use("zviscos") ) CALL iom_put( "zviscos", wmask*zviscos ) 1050 END IF 1161 1051 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 1162 1052 ! Need to put in code for contributions that are applied explicitly to … … 1180 1070 IF(ln_dia_osm) THEN 1181 1071 IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask*zdtdz_pyc ) 1072 IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask*zdsdz_pyc ) 1073 IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask*zdbdz_pyc ) 1182 1074 END IF 1183 1075 … … 1222 1114 END IF ! ln_convmix = .true. 1223 1115 1116 1117 1118 IF ( ln_osm_mle ) THEN ! set up diffusivity and non-gradient mixing 1119 DO_2D( 0, 0, 0, 0 ) 1120 IF ( lflux(ji,jj) ) THEN ! MLE mixing extends below boundary layer 1121 ! Calculate MLE flux contribution from surface fluxes 1122 DO jk = 1, ibld(ji,jj) 1123 znd = gdepw(ji,jj,jk,Kmm) / MAX(zhbl(ji,jj),epsln) 1124 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) - zwth0(ji,jj) * ( 1.0 - znd ) 1125 ghams(ji,jj,jk) = ghams(ji,jj,jk) - zws0(ji,jj) * ( 1.0 - znd ) 1126 END DO 1127 DO jk = 1, mld_prof(ji,jj) 1128 znd = gdepw(ji,jj,jk,Kmm) / MAX(zhmle(ji,jj),epsln) 1129 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth0(ji,jj) * ( 1.0 - znd ) 1130 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws0(ji,jj) * ( 1.0 -znd ) 1131 END DO 1132 ! Viscosity for MLEs 1133 DO jk = 1, mld_prof(ji,jj) 1134 znd = -gdepw(ji,jj,jk,Kmm) / MAX(zhmle(ji,jj),epsln) 1135 zdiffut(ji,jj,jk) = zdiffut(ji,jj,jk) + zdiff_mle(ji,jj) * ( 1.0 - ( 2.0 * znd + 1.0 )**2 ) * ( 1.0 + 5.0 / 21.0 * ( 2.0 * znd + 1.0 )** 2 ) 1136 END DO 1137 ELSE 1138 ! Surface transports limited to OSBL. 1139 ! Viscosity for MLEs 1140 DO jk = 1, mld_prof(ji,jj) 1141 znd = -gdepw(ji,jj,jk,Kmm) / MAX(zhmle(ji,jj),epsln) 1142 zdiffut(ji,jj,jk) = zdiffut(ji,jj,jk) + zdiff_mle(ji,jj) * ( 1.0 - ( 2.0 * znd + 1.0 )**2 ) * ( 1.0 + 5.0 / 21.0 * ( 2.0 * znd + 1.0 )** 2 ) 1143 END DO 1144 ENDIF 1145 END_2D 1146 ENDIF 1147 1148 IF(ln_dia_osm) THEN 1149 IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask*zdtdz_pyc ) 1150 IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask*zdsdz_pyc ) 1151 IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask*zdbdz_pyc ) 1152 END IF 1153 1154 1224 1155 ! Lateral boundary conditions on zvicos (sign unchanged), needed to caclulate viscosities on u and v grids 1225 CALL lbc_lnk( 'zdfosm', zviscos(:,:,:), 'W', 1.0_wp)1156 !CALL lbc_lnk( zviscos(:,:,:), 'W', 1. ) 1226 1157 1227 1158 ! GN 25/8: need to change tmask --> wmask … … 1244 1175 ghams(ji,jj,jk) = ghams(ji,jj,jk) * tmask(ji,jj,jk) 1245 1176 END_3D 1177 ! Lateral boundary conditions on final outputs for hbl, on T-grid (sign unchanged) 1178 CALL lbc_lnk_multi( 'zdfosm', hbl, 'T', 1., dh, 'T', 1., hmle, 'T', 1. ) 1246 1179 ! Lateral boundary conditions on final outputs for gham[ts], on W-grid (sign unchanged) 1247 1180 ! Lateral boundary conditions on final outputs for gham[uv], on [UV]-grid (sign unchanged) … … 1249 1182 & ghamu, 'U', 1.0_wp , ghamv, 'V', 1.0_wp ) 1250 1183 1251 1184 IF(ln_dia_osm) THEN 1252 1185 SELECT CASE (nn_osm_wave) 1253 1186 ! Stokes drift set by assumimg onstant La#=0.3(=0) or Pierson-Moskovitz spectrum (=1). … … 1257 1190 IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rho0*tmask(:,:,1)*zustar**2*zustke ) 1258 1191 ! Stokes drift read in from sbcwave (=2). 1259 CASE(2) 1260 IF ( iom_use("us_x") ) CALL iom_put( "us_x", ut0sd ) ! x surface Stokes drift 1261 IF ( iom_use("us_y") ) CALL iom_put( "us_y", vt0sd ) ! y surface Stokes drift 1192 CASE(2:3) 1193 IF ( iom_use("us_x") ) CALL iom_put( "us_x", ut0sd*umask(:,:,1) ) ! x surface Stokes drift 1194 IF ( iom_use("us_y") ) CALL iom_put( "us_y", vt0sd*vmask(:,:,1) ) ! y surface Stokes drift 1195 IF ( iom_use("wmp") ) CALL iom_put( "wmp", wmp*tmask(:,:,1) ) ! wave mean period 1196 IF ( iom_use("hsw") ) CALL iom_put( "hsw", hsw*tmask(:,:,1) ) ! significant wave height 1197 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 1198 IF ( iom_use("hsw_NP") ) CALL iom_put( "hsw_NP", (0.22/grav)*wndm**2*tmask(:,:,1) ) ! significant wave height from NP spectrum 1199 IF ( iom_use("wndm") ) CALL iom_put( "wndm", wndm*tmask(:,:,1) ) ! U_10 1262 1200 IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rho0*tmask(:,:,1)*zustar**2* & 1263 1201 & SQRT(ut0sd**2 + vt0sd**2 ) ) … … 1270 1208 IF ( iom_use("zws0") ) CALL iom_put( "zws0", tmask(:,:,1)*zws0 ) ! <Sw_0> 1271 1209 IF ( iom_use("hbl") ) CALL iom_put( "hbl", tmask(:,:,1)*hbl ) ! boundary-layer depth 1272 IF ( iom_use("hbli") ) CALL iom_put( "hbli", tmask(:,:,1)*hbli ) ! Initial boundary-layer depth 1210 IF ( iom_use("ibld") ) CALL iom_put( "ibld", tmask(:,:,1)*ibld ) ! boundary-layer max k 1211 IF ( iom_use("zdt_bl") ) CALL iom_put( "zdt_bl", tmask(:,:,1)*zdt_bl ) ! dt at ml base 1212 IF ( iom_use("zds_bl") ) CALL iom_put( "zds_bl", tmask(:,:,1)*zds_bl ) ! ds at ml base 1213 IF ( iom_use("zdb_bl") ) CALL iom_put( "zdb_bl", tmask(:,:,1)*zdb_bl ) ! db at ml base 1214 IF ( iom_use("zdu_bl") ) CALL iom_put( "zdu_bl", tmask(:,:,1)*zdu_bl ) ! du at ml base 1215 IF ( iom_use("zdv_bl") ) CALL iom_put( "zdv_bl", tmask(:,:,1)*zdv_bl ) ! dv at ml base 1216 IF ( iom_use("dh") ) CALL iom_put( "dh", tmask(:,:,1)*dh ) ! Initial boundary-layer depth 1217 IF ( iom_use("hml") ) CALL iom_put( "hml", tmask(:,:,1)*hml ) ! Initial boundary-layer depth 1273 1218 IF ( iom_use("dstokes") ) CALL iom_put( "dstokes", tmask(:,:,1)*dstokes ) ! Stokes drift penetration depth 1274 1219 IF ( iom_use("zustke") ) CALL iom_put( "zustke", tmask(:,:,1)*zustke ) ! Stokes drift magnitude at T-points … … 1276 1221 IF ( iom_use("zwstrl") ) CALL iom_put( "zwstrl", tmask(:,:,1)*zwstrl ) ! Langmuir velocity scale 1277 1222 IF ( iom_use("zustar") ) CALL iom_put( "zustar", tmask(:,:,1)*zustar ) ! friction velocity scale 1223 IF ( iom_use("zvstr") ) CALL iom_put( "zvstr", tmask(:,:,1)*zvstr ) ! mixed velocity scale 1224 IF ( iom_use("zla") ) CALL iom_put( "zla", tmask(:,:,1)*zla ) ! langmuir # 1278 1225 IF ( iom_use("wind_power") ) CALL iom_put( "wind_power", 1000.*rho0*tmask(:,:,1)*zustar**3 ) ! BL depth internal to zdf_osm routine 1279 1226 IF ( iom_use("wind_wave_power") ) CALL iom_put( "wind_wave_power", 1000.*rho0*tmask(:,:,1)*zustar**2*zustke ) 1280 1227 IF ( iom_use("zhbl") ) CALL iom_put( "zhbl", tmask(:,:,1)*zhbl ) ! BL depth internal to zdf_osm routine 1281 1228 IF ( iom_use("zhml") ) CALL iom_put( "zhml", tmask(:,:,1)*zhml ) ! ML depth internal to zdf_osm routine 1282 IF ( iom_use("zdh") ) CALL iom_put( "zdh", tmask(:,:,1)*zdh ) ! ML depth internal to zdf_osm routine 1229 IF ( iom_use("imld") ) CALL iom_put( "imld", tmask(:,:,1)*imld ) ! index for ML depth internal to zdf_osm routine 1230 IF ( iom_use("zdh") ) CALL iom_put( "zdh", tmask(:,:,1)*zdh ) ! pyc thicknessh internal to zdf_osm routine 1283 1231 IF ( iom_use("zhol") ) CALL iom_put( "zhol", tmask(:,:,1)*zhol ) ! ML depth internal to zdf_osm routine 1284 IF ( iom_use("zwthav") ) CALL iom_put( "zwthav", tmask(:,:,1)*zwthav ) ! ML depth internal to zdf_osm routine 1285 IF ( iom_use("zwth_ent") ) CALL iom_put( "zwth_ent", tmask(:,:,1)*zwth_ent ) ! ML depth internal to zdf_osm routine 1286 IF ( iom_use("zt_ml") ) CALL iom_put( "zt_ml", tmask(:,:,1)*zt_ml ) ! average T in ML 1232 IF ( iom_use("zwthav") ) CALL iom_put( "zwthav", tmask(:,:,1)*zwthav ) ! upward BL-avged turb temp flux 1233 IF ( iom_use("zwth_ent") ) CALL iom_put( "zwth_ent", tmask(:,:,1)*zwth_ent ) ! upward turb temp entrainment flux 1234 IF ( iom_use("zwb_ent") ) CALL iom_put( "zwb_ent", tmask(:,:,1)*zwb_ent ) ! upward turb buoyancy entrainment flux 1235 IF ( iom_use("zws_ent") ) CALL iom_put( "zws_ent", tmask(:,:,1)*zws_ent ) ! upward turb salinity entrainment flux 1236 IF ( iom_use("zt_ml") ) CALL iom_put( "zt_ml", tmask(:,:,1)*zt_ml ) ! average T in ML 1237 1238 IF ( iom_use("hmle") ) CALL iom_put( "hmle", tmask(:,:,1)*hmle ) ! FK layer depth 1239 IF ( iom_use("zmld") ) CALL iom_put( "zmld", tmask(:,:,1)*zmld ) ! FK target layer depth 1240 IF ( iom_use("zwb_fk") ) CALL iom_put( "zwb_fk", tmask(:,:,1)*zwb_fk ) ! FK b flux 1241 IF ( iom_use("zwb_fk_b") ) CALL iom_put( "zwb_fk_b", tmask(:,:,1)*zwb_fk_b ) ! FK b flux averaged over ML 1242 IF ( iom_use("mld_prof") ) CALL iom_put( "mld_prof", tmask(:,:,1)*mld_prof )! FK layer max k 1243 IF ( iom_use("zdtdx") ) CALL iom_put( "zdtdx", umask(:,:,1)*zdtdx ) ! FK dtdx at u-pt 1244 IF ( iom_use("zdtdy") ) CALL iom_put( "zdtdy", vmask(:,:,1)*zdtdy ) ! FK dtdy at v-pt 1245 IF ( iom_use("zdsdx") ) CALL iom_put( "zdsdx", umask(:,:,1)*zdsdx ) ! FK dtdx at u-pt 1246 IF ( iom_use("zdsdy") ) CALL iom_put( "zdsdy", vmask(:,:,1)*zdsdy ) ! FK dsdy at v-pt 1247 IF ( iom_use("dbdx_mle") ) CALL iom_put( "dbdx_mle", umask(:,:,1)*dbdx_mle ) ! FK dbdx at u-pt 1248 IF ( iom_use("dbdy_mle") ) CALL iom_put( "dbdy_mle", vmask(:,:,1)*dbdy_mle ) ! FK dbdy at v-pt 1249 IF ( iom_use("zdiff_mle") ) CALL iom_put( "zdiff_mle", tmask(:,:,1)*zdiff_mle )! FK diff in MLE at t-pt 1250 IF ( iom_use("zvel_mle") ) CALL iom_put( "zvel_mle", tmask(:,:,1)*zdiff_mle )! FK diff in MLE at t-pt 1251 1287 1252 END IF 1288 ! Lateral boundary conditions on p_avt (sign unchanged) 1289 CALL lbc_lnk( 'zdfosm', p_avt(:,:,:), 'W', 1.0_wp ) 1253 1254 CONTAINS 1255 ! subroutine code changed, needs syntax checking. 1256 SUBROUTINE zdf_osm_diffusivity_viscosity( zdiffut, zviscos ) 1257 1258 !!--------------------------------------------------------------------- 1259 !! *** ROUTINE zdf_osm_diffusivity_viscosity *** 1260 !! 1261 !! ** Purpose : Determines the eddy diffusivity and eddy viscosity profiles in the mixed layer and the pycnocline. 1262 !! 1263 !! ** Method : 1264 !! 1265 !! !!---------------------------------------------------------------------- 1266 REAL(wp), DIMENSION(:,:,:) :: zdiffut 1267 REAL(wp), DIMENSION(:,:,:) :: zviscos 1268 ! local 1269 1270 ! Scales used to calculate eddy diffusivity and viscosity profiles 1271 REAL(wp), DIMENSION(jpi,jpj) :: zdifml_sc, zvisml_sc 1272 REAL(wp), DIMENSION(jpi,jpj) :: zdifpyc_n_sc, zdifpyc_s_sc, zdifpyc_shr 1273 REAL(wp), DIMENSION(jpi,jpj) :: zvispyc_n_sc, zvispyc_s_sc,zvispyc_shr 1274 REAL(wp), DIMENSION(jpi,jpj) :: zbeta_d_sc, zbeta_v_sc 1275 ! 1276 REAL(wp) :: zvel_sc_pyc, zvel_sc_ml, zstab_fac 1277 1278 REAL(wp), PARAMETER :: rn_dif_ml = 0.8, rn_vis_ml = 0.375 1279 REAL(wp), PARAMETER :: rn_dif_pyc = 0.15, rn_vis_pyc = 0.142 1280 REAL(wp), PARAMETER :: rn_vispyc_shr = 0.15 1281 1282 DO_2D( 0, 0, 0, 0 ) 1283 IF ( lconv(ji,jj) ) THEN 1284 1285 zvel_sc_pyc = ( 0.15 * zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 + 4.25 * zshear(ji,jj) * zhbl(ji,jj) )**pthird 1286 zvel_sc_ml = ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 1287 zstab_fac = ( zhml(ji,jj) / zvel_sc_ml * ( 1.4 - 0.4 / ( 1.0 + EXP(-3.5 * LOG10(-zhol(ji,jj) ) ) )**1.25 ) )**2 1288 1289 zdifml_sc(ji,jj) = rn_dif_ml * zhml(ji,jj) * zvel_sc_ml 1290 zvisml_sc(ji,jj) = rn_vis_ml * zdifml_sc(ji,jj) 1291 1292 IF ( lpyc(ji,jj) ) THEN 1293 zdifpyc_n_sc(ji,jj) = rn_dif_pyc * zvel_sc_ml * zdh(ji,jj) 1294 1295 IF ( lshear(ji,jj) .and. j_ddh(ji,jj) == 1 ) THEN 1296 zdifpyc_n_sc(ji,jj) = zdifpyc_n_sc(ji,jj) + rn_vispyc_shr * ( zshear(ji,jj) * zhbl(ji,jj) )**pthird * zhbl(ji,jj) 1297 ENDIF 1298 1299 zdifpyc_s_sc(ji,jj) = zwb_ent(ji,jj) + 0.0025 * zvel_sc_pyc * ( zhbl(ji,jj) / zdh(ji,jj) - 1.0 ) * ( zb_ml(ji,jj) - zb_bl(ji,jj) ) 1300 zdifpyc_s_sc(ji,jj) = 0.09 * zdifpyc_s_sc(ji,jj) * zstab_fac 1301 zdifpyc_s_sc(ji,jj) = MAX( zdifpyc_s_sc(ji,jj), -0.5 * zdifpyc_n_sc(ji,jj) ) 1302 1303 zvispyc_n_sc(ji,jj) = 0.09 * zvel_sc_pyc * ( 1.0 - zhbl(ji,jj) / zdh(ji,jj) )**2 * ( 0.005 * ( zu_ml(ji,jj)-zu_bl(ji,jj) )**2 + 0.0075 * ( zv_ml(ji,jj)-zv_bl(ji,jj) )**2 ) / zdh(ji,jj) 1304 zvispyc_n_sc(ji,jj) = rn_vis_pyc * zvel_sc_ml * zdh(ji,jj) + zvispyc_n_sc(ji,jj) * zstab_fac 1305 IF ( lshear(ji,jj) .and. j_ddh(ji,jj) == 1 ) THEN 1306 zvispyc_n_sc(ji,jj) = zvispyc_n_sc(ji,jj) + rn_vispyc_shr * ( zshear(ji,jj) * zhbl(ji,jj ) )**pthird * zhbl(ji,jj) 1307 ENDIF 1308 1309 zvispyc_s_sc(ji,jj) = 0.09 * ( zwb_min(ji,jj) + 0.0025 * zvel_sc_pyc * ( zhbl(ji,jj) / zdh(ji,jj) - 1.0 ) * ( zb_ml(ji,jj) - zb_bl(ji,jj) ) ) 1310 zvispyc_s_sc(ji,jj) = zvispyc_s_sc(ji,jj) * zstab_fac 1311 zvispyc_s_sc(ji,jj) = MAX( zvispyc_s_sc(ji,jj), -0.5 * zvispyc_s_sc(ji,jj) ) 1312 1313 zbeta_d_sc(ji,jj) = 1.0 - ( ( zdifpyc_n_sc(ji,jj) + 1.4 * zdifpyc_s_sc(ji,jj) ) / ( zdifml_sc(ji,jj) + epsln ) )**p2third 1314 zbeta_v_sc(ji,jj) = 1.0 - 2.0 * ( zvispyc_n_sc(ji,jj) + zvispyc_s_sc(ji,jj) ) / ( zvisml_sc(ji,jj) + epsln ) 1315 ELSE 1316 zbeta_d_sc = 1.0 1317 zbeta_v_sc = 1.0 1318 ENDIF 1319 ELSE 1320 zdifml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * MAX( EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ), 0.2_wp) 1321 zvisml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * MAX( EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ), 0.2_wp) 1322 END IF 1323 END_2D 1324 ! 1325 DO_2D( 0, 0, 0, 0 ) 1326 IF ( lconv(ji,jj) ) THEN 1327 DO jk = 2, imld(ji,jj) ! mixed layer diffusivity 1328 zznd_ml = gdepw(ji,jj,jk,Kmm) / zhml(ji,jj) 1329 ! 1330 zdiffut(ji,jj,jk) = zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_d_sc(ji,jj) * zznd_ml )**1.5 1331 ! 1332 zviscos(ji,jj,jk) = zvisml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_v_sc(ji,jj) * zznd_ml ) & 1333 & * ( 1.0 - 0.5 * zznd_ml**2 ) 1334 END DO 1335 ! pycnocline 1336 IF ( lpyc(ji,jj) ) THEN 1337 ! Diffusivity profile in the pycnocline given by cubic polynomial. 1338 za_cubic = 0.5 1339 zb_cubic = -1.75 * zdifpyc_s_sc(ji,jj) / zdifpyc_n_sc(ji,jj) 1340 zd_cubic = ( zdh(ji,jj) * zdifml_sc(ji,jj) / zhml(ji,jj) * SQRT( 1.0 - zbeta_d_sc(ji,jj) ) * ( 2.5 * zbeta_d_sc(ji,jj) - 1.0 ) & 1341 & - 0.85 * zdifpyc_s_sc(ji,jj) ) / MAX(zdifpyc_n_sc(ji,jj), 1.e-8) 1342 zd_cubic = zd_cubic - zb_cubic - 2.0 * ( 1.0 - za_cubic - zb_cubic ) 1343 zc_cubic = 1.0 - za_cubic - zb_cubic - zd_cubic 1344 DO jk = imld(ji,jj) , ibld(ji,jj) 1345 zznd_pyc = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / MAX(zdh(ji,jj), 1.e-6) 1346 ! 1347 zdiffut(ji,jj,jk) = zdifpyc_n_sc(ji,jj) * ( za_cubic + zb_cubic * zznd_pyc + zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 ) 1348 1349 zdiffut(ji,jj,jk) = zdiffut(ji,jj,jk) + zdifpyc_s_sc(ji,jj) * ( 1.75 * zznd_pyc - 0.15 * zznd_pyc**2 - 0.2 * zznd_pyc**3 ) 1350 END DO 1351 ! viscosity profiles. 1352 za_cubic = 0.5 1353 zb_cubic = -1.75 * zvispyc_s_sc(ji,jj) / zvispyc_n_sc(ji,jj) 1354 zd_cubic = ( 0.5 * zvisml_sc(ji,jj) * zdh(ji,jj) / zhml(ji,jj) - 0.85 * zvispyc_s_sc(ji,jj) ) / MAX(zvispyc_n_sc(ji,jj), 1.e-8) 1355 zd_cubic = zd_cubic - zb_cubic - 2.0 * ( 1.0 - za_cubic - zd_cubic ) 1356 zc_cubic = 1.0 - za_cubic - zb_cubic - zd_cubic 1357 DO jk = imld(ji,jj) , ibld(ji,jj) 1358 zznd_pyc = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / MAX(zdh(ji,jj), 1.e-6) 1359 zviscos(ji,jj,jk) = zvispyc_n_sc(ji,jj) * ( za_cubic + zb_cubic * zznd_pyc + zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 ) 1360 zviscos(ji,jj,jk) = zviscos(ji,jj,jk) + zvispyc_s_sc(ji,jj) * ( 1.75 * zznd_pyc - 0.15 * zznd_pyc**2 -0.2 * zznd_pyc**3 ) 1361 END DO 1362 IF ( zdhdt(ji,jj) > 0._wp ) THEN 1363 zdiffut(ji,jj,ibld(ji,jj)+1) = MAX( 0.5 * zdhdt(ji,jj) * e3w(ji,jj,ibld(ji,jj)+1,Kmm), 1.0e-6 ) 1364 zviscos(ji,jj,ibld(ji,jj)+1) = MAX( 0.5 * zdhdt(ji,jj) * e3w(ji,jj,ibld(ji,jj)+1,Kmm), 1.0e-6 ) 1365 ELSE 1366 zdiffut(ji,jj,ibld(ji,jj)) = 0._wp 1367 zviscos(ji,jj,ibld(ji,jj)) = 0._wp 1368 ENDIF 1369 ENDIF 1370 ELSE 1371 ! stable conditions 1372 DO jk = 2, ibld(ji,jj) 1373 zznd_ml = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 1374 zdiffut(ji,jj,jk) = 0.75 * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zznd_ml )**1.5 1375 zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * (1.0 - zznd_ml) * ( 1.0 - zznd_ml**2 ) 1376 END DO 1377 1378 IF ( zdhdt(ji,jj) > 0._wp ) THEN 1379 zdiffut(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 1.0e-6) * e3w(ji, jj, ibld(ji,jj), Kmm) 1380 zviscos(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 1.0e-6) * e3w(ji, jj, ibld(ji,jj), Kmm) 1381 ENDIF 1382 ENDIF ! end if ( lconv ) 1383 ! 1384 END_2D 1385 1386 END SUBROUTINE zdf_osm_diffusivity_viscosity 1387 1388 SUBROUTINE zdf_osm_osbl_state( lconv, lshear, j_ddh, zwb_ent, zwb_min, zshear, zri_i ) 1389 1390 !!--------------------------------------------------------------------- 1391 !! *** ROUTINE zdf_osm_osbl_state *** 1392 !! 1393 !! ** Purpose : Determines the state of the OSBL, stable/unstable, shear/ noshear. Also determines shear production, entrainment buoyancy flux and interfacial Richardson number 1394 !! 1395 !! ** Method : 1396 !! 1397 !! !!---------------------------------------------------------------------- 1398 1399 INTEGER, DIMENSION(jpi,jpj) :: j_ddh ! j_ddh = 0, active shear layer; j_ddh=1, shear layer not active; j_ddh=2 shear production low. 1400 1401 LOGICAL, DIMENSION(jpi,jpj) :: lconv, lshear 1402 1403 REAL(wp), DIMENSION(jpi,jpj) :: zwb_ent, zwb_min ! Buoyancy fluxes at base of well-mixed layer. 1404 REAL(wp), DIMENSION(jpi,jpj) :: zshear ! production of TKE due to shear across the pycnocline 1405 REAL(wp), DIMENSION(jpi,jpj) :: zri_i ! Interfacial Richardson Number 1406 1407 ! Local Variables 1408 1409 INTEGER :: jj, ji 1410 1411 REAL(wp), DIMENSION(jpi,jpj) :: zekman 1412 REAL(wp) :: zri_p, zri_b ! Richardson numbers 1413 REAL(wp) :: zshear_u, zshear_v, zwb_shr 1414 REAL(wp) :: zwcor, zrf_conv, zrf_shear, zrf_langmuir, zr_stokes 1415 1416 REAL, PARAMETER :: za_shr = 0.4, zb_shr = 6.5, za_wb_s = 0.1 1417 REAL, PARAMETER :: rn_ri_thres_a = 0.5, rn_ri_thresh_b = 0.59 1418 REAL, PARAMETER :: zalpha_c = 0.2, zalpha_lc = 0.04 1419 REAL, PARAMETER :: zalpha_ls = 0.06, zalpha_s = 0.15 1420 REAL, PARAMETER :: rn_ri_p_thresh = 27.0 1421 REAL, PARAMETER :: zrot=0._wp ! dummy rotation rate of surface stress. 1422 1423 ! Determins stability and set flag lconv 1424 DO_2D( 0, 0, 0, 0 ) 1425 IF ( zhol(ji,jj) < 0._wp ) THEN 1426 lconv(ji,jj) = .TRUE. 1427 ELSE 1428 lconv(ji,jj) = .FALSE. 1429 ENDIF 1430 END_2D 1431 1432 zekman(:,:) = EXP( - 4.0 * ABS( ff_t(:,:) ) * zhbl(:,:) / MAX(zustar(:,:), 1.e-8 ) ) 1433 1434 WHERE ( lconv ) 1435 zri_i = zdb_ml * zhml**2 / MAX( ( zvstr**3 + 0.5 * zwstrc**3 )**p2third * zdh, 1.e-12 ) 1436 END WHERE 1437 1438 zshear(:,:) = 0._wp 1439 j_ddh(:,:) = 1 1440 1441 DO_2D( 0, 0, 0, 0 ) 1442 IF ( lconv(ji,jj) ) THEN 1443 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 1444 zri_p = MAX ( SQRT( zdb_bl(ji,jj) * zdh(ji,jj) / MAX( zdu_bl(ji,jj)**2 + zdv_bl(ji,jj)**2, 1.e-8) ) * ( zhbl(ji,jj) / zdh(ji,jj) ) * ( zvstr(ji,jj) / MAX( zustar(ji,jj), 1.e-6 ) )**2 & 1445 & / MAX( zekman(ji,jj), 1.e-6 ) , 5._wp ) 1446 1447 zri_b = zdb_ml(ji,jj) * zdh(ji,jj) / MAX( zdu_ml(ji,jj)**2 + zdv_ml(ji,jj)**2, 1.e-8 ) 1448 1449 zshear(ji,jj) = za_shr * zekman(ji,jj) * ( MAX( zustar(ji,jj)**2 * zdu_ml(ji,jj) / zhbl(ji,jj), 0._wp ) + zb_shr * MAX( -ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) * zdv_ml(ji,jj) / zhbl(ji,jj), 0._wp ) ) 1450 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1451 ! Test ensures j_ddh=0 is not selected. Change to zri_p<27 when ! 1452 ! full code available ! 1453 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1454 IF ( zri_p < -rn_ri_p_thresh .and. zshear(ji,jj) > 0._wp ) THEN 1455 ! Growing shear layer 1456 j_ddh(ji,jj) = 0 1457 lshear(ji,jj) = .TRUE. 1458 ELSE 1459 j_ddh(ji,jj) = 1 1460 IF ( zri_b <= 1.5 .and. zshear(ji,jj) > 0._wp ) THEN 1461 ! shear production large enough to determine layer charcteristics, but can't maintain a shear layer. 1462 lshear(ji,jj) = .TRUE. 1463 ELSE 1464 ! Shear production may not be zero, but is small and doesn't determine characteristics of pycnocline. 1465 zshear(ji,jj) = 0.5 * zshear(ji,jj) 1466 lshear(ji,jj) = .FALSE. 1467 ENDIF 1468 ENDIF 1469 ELSE ! zdb_bl test, note zshear set to zero 1470 j_ddh(ji,jj) = 2 1471 lshear(ji,jj) = .FALSE. 1472 ENDIF 1473 ENDIF 1474 END_2D 1475 1476 ! Calculate entrainment buoyancy flux due to surface fluxes. 1477 1478 DO_2D( 0, 0, 0, 0 ) 1479 IF ( lconv(ji,jj) ) THEN 1480 zwcor = ABS(ff_t(ji,jj)) * zhbl(ji,jj) + epsln 1481 zrf_conv = TANH( ( zwstrc(ji,jj) / zwcor )**0.69 ) 1482 zrf_shear = TANH( ( zustar(ji,jj) / zwcor )**0.69 ) 1483 zrf_langmuir = TANH( ( zwstrl(ji,jj) / zwcor )**0.69 ) 1484 IF (nn_osm_SD_reduce > 0 ) THEN 1485 ! Effective Stokes drift already reduced from surface value 1486 zr_stokes = 1.0_wp 1487 ELSE 1488 ! Effective Stokes drift only reduced by factor rn_zdfodm_adjust_sd, 1489 ! requires further reduction where BL is deep 1490 zr_stokes = 1.0 - EXP( -25.0 * dstokes(ji,jj) / hbl(ji,jj) & 1491 & * ( 1.0 + 4.0 * dstokes(ji,jj) / hbl(ji,jj) ) ) 1492 END IF 1493 zwb_ent(ji,jj) = - 2.0 * 0.2 * zrf_conv * zwbav(ji,jj) & 1494 & - 0.15 * zrf_shear * zustar(ji,jj)**3 /zhml(ji,jj) & 1495 & + zr_stokes * ( 0.15 * EXP( -1.5 * zla(ji,jj) ) * zrf_shear * zustar(ji,jj)**3 & 1496 & - zrf_langmuir * 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj) 1497 ! 1498 ENDIF 1499 END_2D 1500 1501 zwb_min(:,:) = 0._wp 1502 1503 DO_2D( 0, 0, 0, 0 ) 1504 IF ( lshear(ji,jj) ) THEN 1505 IF ( lconv(ji,jj) ) THEN 1506 ! Unstable OSBL 1507 zwb_shr = -za_wb_s * zshear(ji,jj) 1508 IF ( j_ddh(ji,jj) == 0 ) THEN 1509 1510 ! Developing shear layer, additional shear production possible. 1511 1512 zshear_u = MAX( zustar(ji,jj)**2 * zdu_ml(ji,jj) / zhbl(ji,jj), 0._wp ) 1513 zshear(ji,jj) = zshear(ji,jj) + zshear_u * ( 1.0 - MIN( zri_p / rn_ri_p_thresh, 1.d0 ) ) 1514 zshear(ji,jj) = MIN( zshear(ji,jj), zshear_u ) 1515 1516 zwb_shr = -za_wb_s * zshear(ji,jj) 1517 1518 ENDIF 1519 zwb_ent(ji,jj) = zwb_ent(ji,jj) + zwb_shr 1520 zwb_min(ji,jj) = zwb_ent(ji,jj) + zdh(ji,jj) / zhbl(ji,jj) * zwb0(ji,jj) 1521 ELSE ! IF ( lconv ) THEN - ENDIF 1522 ! Stable OSBL - shear production not coded for first attempt. 1523 ENDIF ! lconv 1524 ELSE ! lshear 1525 IF ( lconv(ji,jj) ) THEN 1526 ! Unstable OSBL 1527 zwb_shr = -za_wb_s * zshear(ji,jj) 1528 zwb_ent(ji,jj) = zwb_ent(ji,jj) + zwb_shr 1529 zwb_min(ji,jj) = zwb_ent(ji,jj) + zdh(ji,jj) / zhbl(ji,jj) * zwb0(ji,jj) 1530 ENDIF ! lconv 1531 ENDIF ! lshear 1532 END_2D 1533 END SUBROUTINE zdf_osm_osbl_state 1534 1535 1536 SUBROUTINE zdf_osm_vertical_average( jnlev_av, jp_ext, zt, zs, zb, zu, zv, zdt, zds, zdb, zdu, zdv ) 1537 !!--------------------------------------------------------------------- 1538 !! *** ROUTINE zdf_vertical_average *** 1539 !! 1540 !! ** Purpose : Determines vertical averages from surface to jnlev. 1541 !! 1542 !! ** Method : Averages are calculated from the surface to jnlev. 1543 !! The external level used to calculate differences is ibld+ibld_ext 1544 !! 1545 !!---------------------------------------------------------------------- 1546 1547 INTEGER, DIMENSION(jpi,jpj) :: jnlev_av ! Number of levels to average over. 1548 INTEGER, DIMENSION(jpi,jpj) :: jp_ext 1549 1550 ! Alan: do we need zb? 1551 REAL(wp), DIMENSION(jpi,jpj) :: zt, zs, zb ! Average temperature and salinity 1552 REAL(wp), DIMENSION(jpi,jpj) :: zu,zv ! Average current components 1553 REAL(wp), DIMENSION(jpi,jpj) :: zdt, zds, zdb ! Difference between average and value at base of OSBL 1554 REAL(wp), DIMENSION(jpi,jpj) :: zdu, zdv ! Difference for velocity components. 1555 1556 INTEGER :: jk, ji, jj, ibld_ext 1557 REAL(wp) :: zthick, zthermal, zbeta 1558 1559 1560 zt = 0._wp 1561 zs = 0._wp 1562 zu = 0._wp 1563 zv = 0._wp 1564 DO_2D( 0, 0, 0, 0 ) 1565 zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 1566 zbeta = rab_n(ji,jj,1,jp_sal) 1567 ! average over depth of boundary layer 1568 zthick = epsln 1569 DO jk = 2, jnlev_av(ji,jj) 1570 zthick = zthick + e3t(ji,jj,jk,Kmm) 1571 zt(ji,jj) = zt(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm) 1572 zs(ji,jj) = zs(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) 1573 zu(ji,jj) = zu(ji,jj) + e3t(ji,jj,jk,Kmm) & 1574 & * ( uu(ji,jj,jk,Kbb) + uu(ji - 1,jj,jk,Kbb) ) & 1575 & / MAX( 1. , umask(ji,jj,jk) + umask(ji - 1,jj,jk) ) 1576 zv(ji,jj) = zv(ji,jj) + e3t(ji,jj,jk,Kmm) & 1577 & * ( vv(ji,jj,jk,Kbb) + vv(ji,jj - 1,jk,Kbb) ) & 1578 & / MAX( 1. , vmask(ji,jj,jk) + vmask(ji,jj - 1,jk) ) 1579 END DO 1580 zt(ji,jj) = zt(ji,jj) / zthick 1581 zs(ji,jj) = zs(ji,jj) / zthick 1582 zu(ji,jj) = zu(ji,jj) / zthick 1583 zv(ji,jj) = zv(ji,jj) / zthick 1584 zb(ji,jj) = grav * zthermal * zt(ji,jj) - grav * zbeta * zs(ji,jj) 1585 ibld_ext = jnlev_av(ji,jj) + jp_ext(ji,jj) 1586 IF ( ibld_ext < mbkt(ji,jj) ) THEN 1587 zdt(ji,jj) = zt(ji,jj) - ts(ji,jj,ibld_ext,jp_tem,Kmm) 1588 zds(ji,jj) = zs(ji,jj) - ts(ji,jj,ibld_ext,jp_sal,Kmm) 1589 zdu(ji,jj) = zu(ji,jj) - ( uu(ji,jj,ibld_ext,Kbb) + uu(ji-1,jj,ibld_ext,Kbb ) ) & 1590 & / MAX(1. , umask(ji,jj,ibld_ext ) + umask(ji-1,jj,ibld_ext ) ) 1591 zdv(ji,jj) = zv(ji,jj) - ( vv(ji,jj,ibld_ext,Kbb) + vv(ji,jj-1,ibld_ext,Kbb ) ) & 1592 & / MAX(1. , vmask(ji,jj,ibld_ext ) + vmask(ji,jj-1,ibld_ext ) ) 1593 zdb(ji,jj) = grav * zthermal * zdt(ji,jj) - grav * zbeta * zds(ji,jj) 1594 ELSE 1595 zdt(ji,jj) = 0._wp 1596 zds(ji,jj) = 0._wp 1597 zdu(ji,jj) = 0._wp 1598 zdv(ji,jj) = 0._wp 1599 zdb(ji,jj) = 0._wp 1600 ENDIF 1601 END_2D 1602 END SUBROUTINE zdf_osm_vertical_average 1603 1604 SUBROUTINE zdf_osm_velocity_rotation( zcos_w, zsin_w, zu, zv, zdu, zdv ) 1605 !!--------------------------------------------------------------------- 1606 !! *** ROUTINE zdf_velocity_rotation *** 1607 !! 1608 !! ** Purpose : Rotates frame of reference of averaged velocity components. 1609 !! 1610 !! ** Method : The velocity components are rotated into frame specified by zcos_w and zsin_w 1611 !! 1612 !!---------------------------------------------------------------------- 1613 1614 REAL(wp), DIMENSION(jpi,jpj) :: zcos_w, zsin_w ! Cos and Sin of rotation angle 1615 REAL(wp), DIMENSION(jpi,jpj) :: zu, zv ! Components of current 1616 REAL(wp), DIMENSION(jpi,jpj) :: zdu, zdv ! Change in velocity components across pycnocline 1617 1618 INTEGER :: ji, jj 1619 REAL(wp) :: ztemp 1620 1621 DO_2D( 0, 0, 0, 0 ) 1622 ztemp = zu(ji,jj) 1623 zu(ji,jj) = zu(ji,jj) * zcos_w(ji,jj) + zv(ji,jj) * zsin_w(ji,jj) 1624 zv(ji,jj) = zv(ji,jj) * zcos_w(ji,jj) - ztemp * zsin_w(ji,jj) 1625 ztemp = zdu(ji,jj) 1626 zdu(ji,jj) = zdu(ji,jj) * zcos_w(ji,jj) + zdv(ji,jj) * zsin_w(ji,jj) 1627 zdv(ji,jj) = zdv(ji,jj) * zsin_w(ji,jj) - ztemp * zsin_w(ji,jj) 1628 END_2D 1629 END SUBROUTINE zdf_osm_velocity_rotation 1630 1631 SUBROUTINE zdf_osm_osbl_state_fk( lpyc, lflux, lmle, zwb_fk ) 1632 !!--------------------------------------------------------------------- 1633 !! *** ROUTINE zdf_osm_osbl_state_fk *** 1634 !! 1635 !! ** Purpose : Determines the state of the OSBL and MLE layer. Info is returned in the logicals lpyc,lflux and lmle. Used with Fox-Kemper scheme. 1636 !! lpyc :: determines whether pycnocline flux-grad relationship needs to be determined 1637 !! lflux :: determines whether effects of surface flux extend below the base of the OSBL 1638 !! lmle :: determines whether the layer with MLE is increasing with time or if base is relaxing towards hbl. 1639 !! 1640 !! ** Method : 1641 !! 1642 !! 1643 !!---------------------------------------------------------------------- 1644 1645 ! Outputs 1646 LOGICAL, DIMENSION(jpi,jpj) :: lpyc, lflux, lmle 1647 REAL(wp), DIMENSION(jpi,jpj) :: zwb_fk 1648 ! 1649 REAL(wp), DIMENSION(jpi,jpj) :: znd_param 1650 REAL(wp) :: zbuoy, ztmp, zpe_mle_layer 1651 REAL(wp) :: zpe_mle_ref, zwb_ent, zdbdz_mle_int 1652 1653 znd_param(:,:) = 0._wp 1654 1655 DO_2D( 0, 0, 0, 0 ) 1656 ztmp = r1_ft(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf 1657 zwb_fk(ji,jj) = rn_osm_mle_ce * hmle(ji,jj) * hmle(ji,jj) * ztmp * zdbds_mle(ji,jj) * zdbds_mle(ji,jj) 1658 END_2D 1659 DO_2D( 0, 0, 0, 0 ) 1660 ! 1661 IF ( lconv(ji,jj) ) THEN 1662 IF ( zhmle(ji,jj) > 1.2 * zhbl(ji,jj) ) THEN 1663 zt_mle(ji,jj) = ( zt_mle(ji,jj) * zhmle(ji,jj) - zt_bl(ji,jj) * zhbl(ji,jj) ) / ( zhmle(ji,jj) - zhbl(ji,jj) ) 1664 zs_mle(ji,jj) = ( zs_mle(ji,jj) * zhmle(ji,jj) - zs_bl(ji,jj) * zhbl(ji,jj) ) / ( zhmle(ji,jj) - zhbl(ji,jj) ) 1665 zb_mle(ji,jj) = ( zb_mle(ji,jj) * zhmle(ji,jj) - zb_bl(ji,jj) * zhbl(ji,jj) ) / ( zhmle(ji,jj) - zhbl(ji,jj) ) 1666 zdbdz_mle_int = ( zb_bl(ji,jj) - ( 2.0 * zb_mle(ji,jj) -zb_bl(ji,jj) ) ) / ( zhmle(ji,jj) - zhbl(ji,jj) ) 1667 ! Calculate potential energies of actual profile and reference profile. 1668 zpe_mle_layer = 0._wp 1669 zpe_mle_ref = 0._wp 1670 DO jk = ibld(ji,jj), mld_prof(ji,jj) 1671 zbuoy = grav * ( zthermal * ts(ji,jj,jk,jp_tem,Kmm) - zbeta * ts(ji,jj,jk,jp_sal,Kmm) ) 1672 zpe_mle_layer = zpe_mle_layer + zbuoy * gdepw(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) 1673 zpe_mle_ref = zpe_mle_ref + ( zb_bl(ji,jj) - zdbdz_mle_int * ( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) ) * gdepw(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) 1674 END DO 1675 ! Non-dimensional parameter to diagnose the presence of thermocline 1676 1677 znd_param(ji,jj) = ( zpe_mle_layer - zpe_mle_ref ) * ABS( ff_t(ji,jj) ) / ( MAX( zwb_fk(ji,jj), 1.0e-10 ) * zhmle(ji,jj) ) 1678 ENDIF 1679 ENDIF 1680 END_2D 1681 1682 ! Diagnosis 1683 DO_2D( 0, 0, 0, 0 ) 1684 IF ( lconv(ji,jj) ) THEN 1685 zwb_ent = - 2.0 * 0.2 * zwbav(ji,jj) & 1686 & - 0.15 * zustar(ji,jj)**3 /zhml(ji,jj) & 1687 & + ( 0.15 * EXP( -1.5 * zla(ji,jj) ) * zustar(ji,jj)**3 & 1688 & - 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj) 1689 IF ( -2.0 * zwb_fk(ji,jj) / zwb_ent > 0.5 ) THEN 1690 IF ( zhmle(ji,jj) > 1.2 * zhbl(ji,jj) ) THEN 1691 ! MLE layer growing 1692 IF ( znd_param (ji,jj) > 100. ) THEN 1693 ! Thermocline present 1694 lflux(ji,jj) = .FALSE. 1695 lmle(ji,jj) =.FALSE. 1696 ELSE 1697 ! Thermocline not present 1698 lflux(ji,jj) = .TRUE. 1699 lmle(ji,jj) = .TRUE. 1700 ENDIF ! znd_param > 100 1701 ! 1702 IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh ) THEN 1703 lpyc(ji,jj) = .FALSE. 1704 ELSE 1705 lpyc = .TRUE. 1706 ENDIF 1707 ELSE 1708 ! MLE layer restricted to OSBL or just below. 1709 IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh ) THEN 1710 ! Weak stratification MLE layer can grow. 1711 lpyc(ji,jj) = .FALSE. 1712 lflux(ji,jj) = .TRUE. 1713 lmle(ji,jj) = .TRUE. 1714 ELSE 1715 ! Strong stratification 1716 lpyc(ji,jj) = .TRUE. 1717 lflux(ji,jj) = .FALSE. 1718 lmle(ji,jj) = .FALSE. 1719 ENDIF ! zdb_bl < rn_mle_thresh_bl and 1720 ENDIF ! zhmle > 1.2 zhbl 1721 ELSE 1722 lpyc(ji,jj) = .TRUE. 1723 lflux(ji,jj) = .FALSE. 1724 lmle(ji,jj) = .FALSE. 1725 IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh ) lpyc(ji,jj) = .FALSE. 1726 ENDIF ! -2.0 * zwb_fk(ji,jj) / zwb_ent > 0.5 1727 ELSE 1728 ! Stable Boundary Layer 1729 lpyc(ji,jj) = .FALSE. 1730 lflux(ji,jj) = .FALSE. 1731 lmle(ji,jj) = .FALSE. 1732 ENDIF ! lconv 1733 END_2D 1734 END SUBROUTINE zdf_osm_osbl_state_fk 1735 1736 SUBROUTINE zdf_osm_external_gradients(jbase, zdtdz, zdsdz, zdbdz ) 1737 !!--------------------------------------------------------------------- 1738 !! *** ROUTINE zdf_osm_external_gradients *** 1739 !! 1740 !! ** Purpose : Calculates the gradients below the OSBL 1741 !! 1742 !! ** Method : Uses ibld and ibld_ext to determine levels to calculate the gradient. 1743 !! 1744 !!---------------------------------------------------------------------- 1745 1746 INTEGER, DIMENSION(jpi,jpj) :: jbase 1747 REAL(wp), DIMENSION(jpi,jpj) :: zdtdz, zdsdz, zdbdz ! External gradients of temperature, salinity and buoyancy. 1748 1749 INTEGER :: jj, ji, jkb, jkb1 1750 REAL(wp) :: zthermal, zbeta 1751 1752 1753 DO_2D( 0, 0, 0, 0 ) 1754 IF ( jbase(ji,jj)+1 < mbkt(ji,jj) ) THEN 1755 zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 1756 zbeta = rab_n(ji,jj,1,jp_sal) 1757 jkb = jbase(ji,jj) 1758 jkb1 = MIN(jkb + 1, mbkt(ji,jj)) 1759 zdtdz(ji,jj) = - ( ts(ji,jj,jkb1,jp_tem,Kmm) - ts(ji,jj,jkb,jp_tem,Kmm ) ) & 1760 & / e3t(ji,jj,ibld(ji,jj),Kmm) 1761 zdsdz(ji,jj) = - ( ts(ji,jj,jkb1,jp_sal,Kmm) - ts(ji,jj,jkb,jp_sal,Kmm ) ) & 1762 & / e3t(ji,jj,ibld(ji,jj),Kmm) 1763 zdbdz(ji,jj) = grav * zthermal * zdtdz(ji,jj) - grav * zbeta * zdsdz(ji,jj) 1764 ELSE 1765 zdtdz(ji,jj) = 0._wp 1766 zdsdz(ji,jj) = 0._wp 1767 zdbdz(ji,jj) = 0._wp 1768 END IF 1769 END_2D 1770 END SUBROUTINE zdf_osm_external_gradients 1771 1772 SUBROUTINE zdf_osm_pycnocline_scalar_profiles( zdtdz, zdsdz, zdbdz, zalpha ) 1773 1774 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdtdz, zdsdz, zdbdz ! gradients in the pycnocline 1775 REAL(wp), DIMENSION(jpi,jpj) :: zalpha 1776 1777 INTEGER :: jk, jj, ji 1778 REAL(wp) :: ztgrad, zsgrad, zbgrad 1779 REAL(wp) :: zgamma_b_nd, znd 1780 REAL(wp) :: zzeta_m, zzeta_en, zbuoy_pyc_sc 1781 REAL(wp), PARAMETER :: zgamma_b = 2.25, zzeta_sh = 0.15 1782 1783 DO_2D( 0, 0, 0, 0 ) 1784 IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN 1785 IF ( lconv(ji,jj) ) THEN ! convective conditions 1786 IF ( lpyc(ji,jj) ) THEN 1787 zzeta_m = 0.1 + 0.3 / ( 1.0 + EXP( -3.5 * LOG10( -zhol(ji,jj) ) ) ) 1788 zalpha(ji,jj) = 2.0 * ( 1.0 - ( 0.80 * zzeta_m + 0.5 * SQRT( 3.14159 / zgamma_b ) ) * zdbdz_bl_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj) ) / ( 0.723 + SQRT( 3.14159 / zgamma_b ) ) 1789 zalpha(ji,jj) = MAX( zalpha(ji,jj), 0._wp ) 1790 1791 ztmp = 1._wp/MAX(zdh(ji,jj), epsln) 1792 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1793 ! Commented lines in this section are not needed in new code, once tested ! 1794 ! can be removed ! 1795 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1796 ! ztgrad = zalpha * zdt_ml(ji,jj) * ztmp + zdtdz_bl_ext(ji,jj) 1797 ! zsgrad = zalpha * zds_ml(ji,jj) * ztmp + zdsdz_bl_ext(ji,jj) 1798 zbgrad = zalpha(ji,jj) * zdb_ml(ji,jj) * ztmp + zdbdz_bl_ext(ji,jj) 1799 zgamma_b_nd = zdbdz_bl_ext(ji,jj) * zdh(ji,jj) / MAX(zdb_ml(ji,jj), epsln) 1800 DO jk = 2, ibld(ji,jj)+ibld_ext 1801 znd = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) * ztmp 1802 IF ( znd <= zzeta_m ) THEN 1803 ! zdtdz(ji,jj,jk) = zdtdz_bl_ext(ji,jj) + zalpha * zdt_ml(ji,jj) * ztmp * & 1804 ! & EXP( -6.0 * ( znd -zzeta_m )**2 ) 1805 ! zdsdz(ji,jj,jk) = zdsdz_bl_ext(ji,jj) + zalpha * zds_ml(ji,jj) * ztmp * & 1806 ! & EXP( -6.0 * ( znd -zzeta_m )**2 ) 1807 zdbdz(ji,jj,jk) = zdbdz_bl_ext(ji,jj) + zalpha(ji,jj) * zdb_ml(ji,jj) * ztmp * & 1808 & EXP( -6.0 * ( znd -zzeta_m )**2 ) 1809 ELSE 1810 ! zdtdz(ji,jj,jk) = ztgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 ) 1811 ! zdsdz(ji,jj,jk) = zsgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 ) 1812 zdbdz(ji,jj,jk) = zbgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 ) 1813 ENDIF 1814 END DO 1815 ENDIF ! if no pycnocline pycnocline gradients set to zero 1816 ELSE 1817 ! stable conditions 1818 ! if pycnocline profile only defined when depth steady of increasing. 1819 IF ( zdhdt(ji,jj) > 0.0 ) THEN ! Depth increasing, or steady. 1820 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 1821 IF ( zhol(ji,jj) >= 0.5 ) THEN ! Very stable - 'thick' pycnocline 1822 ztmp = 1._wp/MAX(zhbl(ji,jj), epsln) 1823 ztgrad = zdt_bl(ji,jj) * ztmp 1824 zsgrad = zds_bl(ji,jj) * ztmp 1825 zbgrad = zdb_bl(ji,jj) * ztmp 1826 DO jk = 2, ibld(ji,jj) 1827 znd = gdepw(ji,jj,jk,Kmm) * ztmp 1828 zdtdz(ji,jj,jk) = ztgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 1829 zdbdz(ji,jj,jk) = zbgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 1830 zdsdz(ji,jj,jk) = zsgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 1831 END DO 1832 ELSE ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form. 1833 ztmp = 1._wp/MAX(zdh(ji,jj), epsln) 1834 ztgrad = zdt_bl(ji,jj) * ztmp 1835 zsgrad = zds_bl(ji,jj) * ztmp 1836 zbgrad = zdb_bl(ji,jj) * ztmp 1837 DO jk = 2, ibld(ji,jj) 1838 znd = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) * ztmp 1839 zdtdz(ji,jj,jk) = ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 1840 zdbdz(ji,jj,jk) = zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 1841 zdsdz(ji,jj,jk) = zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 1842 END DO 1843 ENDIF ! IF (zhol >=0.5) 1844 ENDIF ! IF (zdb_bl> 0.) 1845 ENDIF ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero 1846 ENDIF ! IF (lconv) 1847 ENDIF ! IF ( ibld(ji,jj) < mbkt(ji,jj) ) 1848 END_2D 1849 1850 END SUBROUTINE zdf_osm_pycnocline_scalar_profiles 1851 1852 SUBROUTINE zdf_osm_pycnocline_shear_profiles( zdudz, zdvdz ) 1853 !!--------------------------------------------------------------------- 1854 !! *** ROUTINE zdf_osm_pycnocline_shear_profiles *** 1855 !! 1856 !! ** Purpose : Calculates velocity shear in the pycnocline 1857 !! 1858 !! ** Method : 1859 !! 1860 !!---------------------------------------------------------------------- 1861 1862 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdudz, zdvdz 1863 1864 INTEGER :: jk, jj, ji 1865 REAL(wp) :: zugrad, zvgrad, znd 1866 REAL(wp) :: zzeta_v = 0.45 1290 1867 ! 1291 END SUBROUTINE zdf_osm 1292 1293 1294 SUBROUTINE zdf_osm_init( Kmm ) 1868 DO_2D( 0, 0, 0, 0 ) 1869 ! 1870 IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN 1871 IF ( lconv (ji,jj) ) THEN 1872 ! Unstable conditions. Shouldn;t be needed with no pycnocline code. 1873 ! zugrad = 0.7 * zdu_ml(ji,jj) / zdh(ji,jj) + 0.3 * zustar(ji,jj)*zustar(ji,jj) / & 1874 ! & ( ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) * & 1875 ! & MIN(zla(ji,jj)**(8.0/3.0) + epsln, 0.12 )) 1876 !Alan is this right? 1877 ! zvgrad = ( 0.7 * zdv_ml(ji,jj) + & 1878 ! & 2.0 * ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) / & 1879 ! & ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird + epsln ) & 1880 ! & )/ (zdh(ji,jj) + epsln ) 1881 ! DO jk = 2, ibld(ji,jj) - 1 + ibld_ext 1882 ! znd = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / (zdh(ji,jj) + epsln ) - zzeta_v 1883 ! IF ( znd <= 0.0 ) THEN 1884 ! zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( 3.0 * znd ) 1885 ! zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( 3.0 * znd ) 1886 ! ELSE 1887 ! zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( -2.0 * znd ) 1888 ! zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( -2.0 * znd ) 1889 ! ENDIF 1890 ! END DO 1891 ELSE 1892 ! stable conditions 1893 zugrad = 3.25 * zdu_bl(ji,jj) / zhbl(ji,jj) 1894 zvgrad = 2.75 * zdv_bl(ji,jj) / zhbl(ji,jj) 1895 DO jk = 2, ibld(ji,jj) 1896 znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 1897 IF ( znd < 1.0 ) THEN 1898 zdudz(ji,jj,jk) = zugrad * EXP( -40.0 * ( znd - 1.0 )**2 ) 1899 ELSE 1900 zdudz(ji,jj,jk) = zugrad * EXP( -20.0 * ( znd - 1.0 )**2 ) 1901 ENDIF 1902 zdvdz(ji,jj,jk) = zvgrad * EXP( -20.0 * ( znd - 0.85 )**2 ) 1903 END DO 1904 ENDIF 1905 ! 1906 END IF ! IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) 1907 END_2D 1908 END SUBROUTINE zdf_osm_pycnocline_shear_profiles 1909 1910 SUBROUTINE zdf_osm_calculate_dhdt( zdhdt, zddhdt ) 1911 !!--------------------------------------------------------------------- 1912 !! *** ROUTINE zdf_osm_calculate_dhdt *** 1913 !! 1914 !! ** Purpose : Calculates the rate at which hbl changes. 1915 !! 1916 !! ** Method : 1917 !! 1918 !!---------------------------------------------------------------------- 1919 1920 REAL(wp), DIMENSION(jpi,jpj) :: zdhdt, zddhdt ! Rate of change of hbl 1921 1922 INTEGER :: jj, ji 1923 REAL(wp) :: zgamma_b_nd, zgamma_dh_nd, zpert, zpsi 1924 REAL(wp) :: zvel_max!, zwb_min 1925 REAL(wp) :: zzeta_m = 0.3 1926 REAL(wp) :: zgamma_c = 2.0 1927 REAL(wp) :: zdhoh = 0.1 1928 REAL(wp) :: alpha_bc = 0.5 1929 REAL, PARAMETER :: a_ddh = 2.5, a_ddh_2 = 3.5 ! also in pycnocline_depth 1930 1931 DO_2D( 0, 0, 0, 0 ) 1932 1933 IF ( lshear(ji,jj) ) THEN 1934 IF ( lconv(ji,jj) ) THEN ! Convective 1935 1936 IF ( ln_osm_mle ) THEN 1937 1938 IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN 1939 ! Fox-Kemper buoyancy flux average over OSBL 1940 zwb_fk_b(ji,jj) = zwb_fk(ji,jj) * & 1941 (1.0 + hmle(ji,jj) / ( 6.0 * hbl(ji,jj) ) * (-1.0 + ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3) ) 1942 ELSE 1943 zwb_fk_b(ji,jj) = 0.5 * zwb_fk(ji,jj) * hmle(ji,jj) / hbl(ji,jj) 1944 ENDIF 1945 zvel_max = ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**p2third / hbl(ji,jj) 1946 IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0.0 ) THEN 1947 ! OSBL is deepening, entrainment > restratification 1948 IF ( zdb_bl(ji,jj) > 0.0 .and. zdbdz_bl_ext(ji,jj) > 0.0 ) THEN 1949 ! *** Used for shear Needs to be changed to work stabily 1950 ! zgamma_b_nd = zdbdz_bl_ext * dh / zdb_ml 1951 ! zalpha_b = 6.7 * zgamma_b_nd / ( 1.0 + zgamma_b_nd ) 1952 ! zgamma_b = zgamma_b_nd / ( 0.12 * ( 1.25 + zgamma_b_nd ) ) 1953 ! za_1 = 1.0 / zgamma_b**2 - 0.017 1954 ! za_2 = 1.0 / zgamma_b**3 - 0.0025 1955 ! zpsi = zalpha_b * ( 1.0 + zgamma_b_nd ) * ( za_1 - 2.0 * za_2 * dh / hbl ) 1956 zpsi = 0._wp 1957 zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) 1958 zdhdt(ji,jj) = zdhdt(ji,jj)! - zpsi * ( -1.0 / zhml(ji,jj) + 2.4 * zdbdz_bl_ext(ji,jj) / zdb_ml(ji,jj) ) * zwb_min(ji,jj) * zdh(ji,jj) / zdb_bl(ji,jj) 1959 IF ( j_ddh(ji,jj) == 1 ) THEN 1960 IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN 1961 zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zvstr(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) 1962 ELSE 1963 zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zwstrc(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) 1964 ENDIF 1965 ! Relaxation to dh_ref = zari * hbl 1966 zddhdt(ji,jj) = -a_ddh_2 * ( 1.0 - zdh(ji,jj) / ( zari * zhbl(ji,jj) ) ) * zwb_ent(ji,jj) / zdb_bl(ji,jj) 1967 1968 ELSE ! j_ddh == 0 1969 ! Growing shear layer 1970 zddhdt(ji,jj) = -a_ddh * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zwb_ent(ji,jj) / zdb_bl(ji,jj) 1971 ENDIF ! j_ddh 1972 zdhdt(ji,jj) = zdhdt(ji,jj) ! + zpsi * zddhdt(ji,jj) 1973 ELSE ! zdb_bl >0 1974 zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / MAX( zvel_max, 1.0e-15) 1975 ENDIF 1976 ELSE ! zwb_min + 2*zwb_fk_b < 0 1977 ! OSBL shoaling due to restratification flux. This is the velocity defined in Fox-Kemper et al (2008) 1978 zdhdt(ji,jj) = - zvel_mle(ji,jj) 1979 1980 1981 ENDIF 1982 1983 ELSE 1984 ! Fox-Kemper not used. 1985 1986 zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_Dt / hbl(ji,jj) ) * zwb_ent(ji,jj) / & 1987 & MAX((zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird, epsln) 1988 zdhdt(ji,jj) = -zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) 1989 ! added ajgn 23 July as temporay fix 1990 1991 ENDIF ! ln_osm_mle 1992 1993 ELSE ! lconv - Stable 1994 zdhdt(ji,jj) = ( 0.06 + 0.52 * zhol(ji,jj) / 2.0 ) * zvstr(ji,jj)**3 / hbl(ji,jj) + zwbav(ji,jj) 1995 IF ( zdhdt(ji,jj) < 0._wp ) THEN 1996 ! For long timsteps factor in brackets slows the rapid collapse of the OSBL 1997 zpert = 2.0 * ( 1.0 + 0.0 * 2.0 * zvstr(ji,jj) * rn_Dt / hbl(ji,jj) ) * zvstr(ji,jj)**2 / hbl(ji,jj) 1998 ELSE 1999 zpert = MAX( 2.0 * zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) ) 2000 ENDIF 2001 zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / MAX(zpert, epsln) 2002 ENDIF ! lconv 2003 ELSE ! lshear 2004 IF ( lconv(ji,jj) ) THEN ! Convective 2005 2006 IF ( ln_osm_mle ) THEN 2007 2008 IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN 2009 ! Fox-Kemper buoyancy flux average over OSBL 2010 zwb_fk_b(ji,jj) = zwb_fk(ji,jj) * & 2011 (1.0 + hmle(ji,jj) / ( 6.0 * hbl(ji,jj) ) * (-1.0 + ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3) ) 2012 ELSE 2013 zwb_fk_b(ji,jj) = 0.5 * zwb_fk(ji,jj) * hmle(ji,jj) / hbl(ji,jj) 2014 ENDIF 2015 zvel_max = ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**p2third / hbl(ji,jj) 2016 IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0.0 ) THEN 2017 ! OSBL is deepening, entrainment > restratification 2018 IF ( zdb_bl(ji,jj) > 0.0 .and. zdbdz_bl_ext(ji,jj) > 0.0 ) THEN 2019 zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) 2020 ELSE 2021 zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / MAX( zvel_max, 1.0e-15) 2022 ENDIF 2023 ELSE 2024 ! OSBL shoaling due to restratification flux. This is the velocity defined in Fox-Kemper et al (2008) 2025 zdhdt(ji,jj) = - zvel_mle(ji,jj) 2026 2027 2028 ENDIF 2029 2030 ELSE 2031 ! Fox-Kemper not used. 2032 2033 zvel_max = -zwb_ent(ji,jj) / & 2034 & MAX((zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird, epsln) 2035 zdhdt(ji,jj) = -zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) 2036 ! added ajgn 23 July as temporay fix 2037 2038 ENDIF ! ln_osm_mle 2039 2040 ELSE ! Stable 2041 zdhdt(ji,jj) = ( 0.06 + 0.52 * zhol(ji,jj) / 2.0 ) * zvstr(ji,jj)**3 / hbl(ji,jj) + zwbav(ji,jj) 2042 IF ( zdhdt(ji,jj) < 0._wp ) THEN 2043 ! For long timsteps factor in brackets slows the rapid collapse of the OSBL 2044 zpert = 2.0 * zvstr(ji,jj)**2 / hbl(ji,jj) 2045 ELSE 2046 zpert = MAX( 2.0 * zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) ) 2047 ENDIF 2048 zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / MAX(zpert, epsln) 2049 ENDIF ! lconv 2050 ENDIF ! lshear 2051 END_2D 2052 END SUBROUTINE zdf_osm_calculate_dhdt 2053 2054 SUBROUTINE zdf_osm_timestep_hbl( zdhdt ) 2055 !!--------------------------------------------------------------------- 2056 !! *** ROUTINE zdf_osm_timestep_hbl *** 2057 !! 2058 !! ** Purpose : Increments hbl. 2059 !! 2060 !! ** Method : If thechange in hbl exceeds one model level the change is 2061 !! is calculated by moving down the grid, changing the buoyancy 2062 !! jump. This is to ensure that the change in hbl does not 2063 !! overshoot a stable layer. 2064 !! 2065 !!---------------------------------------------------------------------- 2066 2067 2068 REAL(wp), DIMENSION(jpi,jpj) :: zdhdt ! rates of change of hbl. 2069 2070 INTEGER :: jk, jj, ji, jm 2071 REAL(wp) :: zhbl_s, zvel_max, zdb 2072 REAL(wp) :: zthermal, zbeta 2073 2074 DO_2D( 0, 0, 0, 0 ) 2075 IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN 2076 ! 2077 ! If boundary layer changes by more than one level, need to check for stable layers between initial and final depths. 2078 ! 2079 zhbl_s = hbl(ji,jj) 2080 jm = imld(ji,jj) 2081 zthermal = rab_n(ji,jj,1,jp_tem) 2082 zbeta = rab_n(ji,jj,1,jp_sal) 2083 2084 2085 IF ( lconv(ji,jj) ) THEN 2086 !unstable 2087 2088 IF( ln_osm_mle ) THEN 2089 zvel_max = ( zwstrl(ji,jj)**3 + zwstrc(ji,jj)**3 )**p2third / hbl(ji,jj) 2090 ELSE 2091 2092 zvel_max = -( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_Dt / hbl(ji,jj) ) * zwb_ent(ji,jj) / & 2093 & ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 2094 2095 ENDIF 2096 2097 DO jk = imld(ji,jj), ibld(ji,jj) 2098 zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - ts(ji,jj,jm,jp_tem,Kmm) ) & 2099 & - zbeta * ( zs_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ), & 2100 & 0.0 ) + zvel_max 2101 2102 2103 IF ( ln_osm_mle ) THEN 2104 zhbl_s = zhbl_s + MIN( & 2105 & rn_Dt * ( ( -zwb_ent(ji,jj) - 2.0 * zwb_fk_b(ji,jj) )/ zdb ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), & 2106 & e3w(ji,jj,jm,Kmm) ) 2107 ELSE 2108 zhbl_s = zhbl_s + MIN( & 2109 & rn_Dt * ( -zwb_ent(ji,jj) / zdb ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), & 2110 & e3w(ji,jj,jm,Kmm) ) 2111 ENDIF 2112 2113 ! zhbl_s = MIN(zhbl_s, gdepw(ji,jj, mbkt(ji,jj) + 1,Kmm) - depth_tol) 2114 IF ( zhbl_s >= gdepw(ji,jj,mbkt(ji,jj) + 1,Kmm) ) THEN 2115 zhbl_s = MIN(zhbl_s, gdepw(ji,jj, mbkt(ji,jj) + 1,Kmm) - depth_tol) 2116 lpyc(ji,jj) = .FALSE. 2117 ENDIF 2118 IF ( zhbl_s >= gdepw(ji,jj,jm+1,Kmm) ) jm = jm + 1 2119 END DO 2120 hbl(ji,jj) = zhbl_s 2121 ibld(ji,jj) = jm 2122 ELSE 2123 ! stable 2124 DO jk = imld(ji,jj), ibld(ji,jj) 2125 zdb = MAX( & 2126 & grav * ( zthermal * ( zt_bl(ji,jj) - ts(ji,jj,jm,jp_tem,Kmm) )& 2127 & - zbeta * ( zs_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ),& 2128 & 0.0 ) + & 2129 & 2.0 * zvstr(ji,jj)**2 / zhbl_s 2130 2131 ! Alan is thuis right? I have simply changed hbli to hbl 2132 zhol(ji,jj) = -zhbl_s / ( ( zvstr(ji,jj)**3 + epsln )/ zwbav(ji,jj) ) 2133 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) ) ) * & 2134 & zustar(ji,jj)**3 / zhbl_s ) * ( 0.725 + 0.225 * EXP( -7.5 * zhol(ji,jj) ) ) 2135 zdhdt(ji,jj) = zdhdt(ji,jj) + zwbav(ji,jj) 2136 zhbl_s = zhbl_s + MIN( zdhdt(ji,jj) / zdb * rn_Dt / FLOAT( ibld(ji,jj) - imld(ji,jj) ), e3w(ji,jj,jm,Kmm) ) 2137 2138 ! zhbl_s = MIN(zhbl_s, gdepw(ji,jj, mbkt(ji,jj) + 1,Kmm) - depth_tol) 2139 IF ( zhbl_s >= mbkt(ji,jj) + 1 ) THEN 2140 zhbl_s = MIN(zhbl_s, gdepw(ji,jj, mbkt(ji,jj) + 1,Kmm) - depth_tol) 2141 lpyc(ji,jj) = .FALSE. 2142 ENDIF 2143 IF ( zhbl_s >= gdepw(ji,jj,jm,Kmm) ) jm = jm + 1 2144 END DO 2145 ENDIF ! IF ( lconv ) 2146 hbl(ji,jj) = MAX(zhbl_s, gdepw(ji,jj,4,Kmm) ) 2147 ibld(ji,jj) = MAX(jm, 4 ) 2148 ELSE 2149 ! change zero or one model level. 2150 hbl(ji,jj) = MAX(zhbl_t(ji,jj), gdepw(ji,jj,4,Kmm) ) 2151 ENDIF 2152 zhbl(ji,jj) = gdepw(ji,jj,ibld(ji,jj),Kmm) 2153 END_2D 2154 2155 END SUBROUTINE zdf_osm_timestep_hbl 2156 2157 SUBROUTINE zdf_osm_pycnocline_thickness( dh, zdh ) 2158 !!--------------------------------------------------------------------- 2159 !! *** ROUTINE zdf_osm_pycnocline_thickness *** 2160 !! 2161 !! ** Purpose : Calculates thickness of the pycnocline 2162 !! 2163 !! ** Method : The thickness is calculated from a prognostic equation 2164 !! that relaxes the pycnocine thickness to a diagnostic 2165 !! value. The time change is calculated assuming the 2166 !! thickness relaxes exponentially. This is done to deal 2167 !! with large timesteps. 2168 !! 2169 !!---------------------------------------------------------------------- 2170 2171 REAL(wp), DIMENSION(jpi,jpj) :: dh, zdh ! pycnocline thickness. 2172 ! 2173 INTEGER :: jj, ji 2174 INTEGER :: inhml 2175 REAL(wp) :: zari, ztau, zdh_ref 2176 REAL, PARAMETER :: a_ddh_2 = 3.5 ! also in pycnocline_depth 2177 2178 DO_2D( 0, 0, 0, 0 ) 2179 2180 IF ( lshear(ji,jj) ) THEN 2181 IF ( lconv(ji,jj) ) THEN 2182 IF ( j_ddh(ji,jj) == 0 ) THEN 2183 ! ddhdt for pycnocline determined in osm_calculate_dhdt 2184 dh(ji,jj) = dh(ji,jj) + zddhdt(ji,jj) * rn_Dt 2185 ELSE 2186 ! Temporary (probably) Recalculate dh_ref to ensure dh doesn't go negative. Can't do this using zddhdt from calculate_dhdt 2187 IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN 2188 zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zvstr(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) 2189 ELSE 2190 zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zwstrc(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) 2191 ENDIF 2192 ztau = MAX( zdb_bl(ji,jj) * ( zari * hbl(ji,jj) ) / ( a_ddh_2 * MAX(-zwb_ent(ji,jj), 1.e-12) ), 2.0 * rn_Dt ) 2193 dh(ji,jj) = dh(ji,jj) * EXP( -rn_Dt / ztau ) + zari * zhbl(ji,jj) * ( 1.0 - EXP( -rn_Dt / ztau ) ) 2194 IF ( dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zari * zhbl(ji,jj) 2195 ENDIF 2196 2197 ELSE ! lconv 2198 ! Initially shear only for entraining OSBL. Stable code will be needed if extended to stable OSBL 2199 2200 ztau = hbl(ji,jj) / MAX(zvstr(ji,jj), epsln) 2201 IF ( zdhdt(ji,jj) >= 0.0 ) THEN ! probably shouldn't include wm here 2202 ! boundary layer deepening 2203 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 2204 ! pycnocline thickness set by stratification - use same relationship as for neutral conditions. 2205 zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) & 2206 & / MAX(zdb_bl(ji,jj) * zhbl(ji,jj), epsln ) + 0.01 , 0.2 ) 2207 zdh_ref = MIN( zari, 0.2 ) * hbl(ji,jj) 2208 ELSE 2209 zdh_ref = 0.2 * hbl(ji,jj) 2210 ENDIF 2211 ELSE ! IF(dhdt < 0) 2212 zdh_ref = 0.2 * hbl(ji,jj) 2213 ENDIF ! IF (dhdt >= 0) 2214 dh(ji,jj) = dh(ji,jj) * EXP( -rn_Dt / ztau )+ zdh_ref * ( 1.0 - EXP( -rn_Dt / ztau ) ) 2215 IF ( zdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref ! can be a problem with dh>hbl for rapid collapse 2216 ! Alan: this hml is never defined or used -- do we need it? 2217 ENDIF 2218 2219 ELSE ! lshear 2220 ! for lshear = .FALSE. calculate ddhdt here 2221 2222 IF ( lconv(ji,jj) ) THEN 2223 2224 IF( ln_osm_mle ) THEN 2225 IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0._wp ) THEN 2226 ! OSBL is deepening. Note wb_fk_b is zero if ln_osm_mle=F 2227 IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_bl_ext(ji,jj) > 0._wp)THEN 2228 IF ( ( zwstrc(ji,jj) / MAX(zvstr(ji,jj), epsln) )**3 <= 0.5 ) THEN ! near neutral stability 2229 zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zvstr(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) 2230 ELSE ! unstable 2231 zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zwstrc(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) 2232 ENDIF 2233 ztau = 0.2 * hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird) 2234 zdh_ref = zari * hbl(ji,jj) 2235 ELSE 2236 ztau = 0.2 * hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird) 2237 zdh_ref = 0.2 * hbl(ji,jj) 2238 ENDIF 2239 ELSE 2240 ztau = 0.2 * hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird) 2241 zdh_ref = 0.2 * hbl(ji,jj) 2242 ENDIF 2243 ELSE ! ln_osm_mle 2244 IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_bl_ext(ji,jj) > 0._wp)THEN 2245 IF ( ( zwstrc(ji,jj) / MAX(zvstr(ji,jj), epsln) )**3 <= 0.5 ) THEN ! near neutral stability 2246 zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zvstr(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) 2247 ELSE ! unstable 2248 zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zwstrc(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) 2249 ENDIF 2250 ztau = hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird) 2251 zdh_ref = zari * hbl(ji,jj) 2252 ELSE 2253 ztau = hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird) 2254 zdh_ref = 0.2 * hbl(ji,jj) 2255 ENDIF 2256 2257 END IF ! ln_osm_mle 2258 2259 dh(ji,jj) = dh(ji,jj) * EXP( -rn_Dt / ztau ) + zdh_ref * ( 1.0 - EXP( -rn_Dt / ztau ) ) 2260 ! IF ( zdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref 2261 IF ( dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref 2262 ! Alan: this hml is never defined or used 2263 ELSE ! IF (lconv) 2264 ztau = hbl(ji,jj) / MAX(zvstr(ji,jj), epsln) 2265 IF ( zdhdt(ji,jj) >= 0.0 ) THEN ! probably shouldn't include wm here 2266 ! boundary layer deepening 2267 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 2268 ! pycnocline thickness set by stratification - use same relationship as for neutral conditions. 2269 zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) & 2270 & / MAX(zdb_bl(ji,jj) * zhbl(ji,jj), epsln ) + 0.01 , 0.2 ) 2271 zdh_ref = MIN( zari, 0.2 ) * hbl(ji,jj) 2272 ELSE 2273 zdh_ref = 0.2 * hbl(ji,jj) 2274 ENDIF 2275 ELSE ! IF(dhdt < 0) 2276 zdh_ref = 0.2 * hbl(ji,jj) 2277 ENDIF ! IF (dhdt >= 0) 2278 dh(ji,jj) = dh(ji,jj) * EXP( -rn_Dt / ztau )+ zdh_ref * ( 1.0 - EXP( -rn_Dt / ztau ) ) 2279 IF ( zdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref ! can be a problem with dh>hbl for rapid collapse 2280 ENDIF ! IF (lconv) 2281 ENDIF ! lshear 2282 2283 hml(ji,jj) = hbl(ji,jj) - dh(ji,jj) 2284 inhml = MAX( INT( dh(ji,jj) / MAX(e3t(ji,jj,ibld(ji,jj),Kmm), 1.e-3) ) , 1 ) 2285 imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 3) 2286 zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm) 2287 zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 2288 END_2D 2289 2290 END SUBROUTINE zdf_osm_pycnocline_thickness 2291 2292 2293 SUBROUTINE zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle, zdbds_mle ) 2294 !!---------------------------------------------------------------------- 2295 !! *** ROUTINE zdf_osm_horizontal_gradients *** 2296 !! 2297 !! ** Purpose : Calculates horizontal gradients of buoyancy for use with Fox-Kemper parametrization. 2298 !! 2299 !! ** Method : 2300 !! 2301 !! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008 2302 !! Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008 2303 2304 2305 REAL(wp), DIMENSION(jpi,jpj) :: dbdx_mle, dbdy_mle ! MLE horiz gradients at u & v points 2306 REAL(wp), DIMENSION(jpi,jpj) :: zdbds_mle ! Magnitude of horizontal buoyancy gradient. 2307 REAL(wp), DIMENSION(jpi,jpj) :: zmld ! == estimated FK BLD used for MLE horiz gradients == ! 2308 REAL(wp), DIMENSION(jpi,jpj) :: zdtdx, zdtdy, zdsdx, zdsdy 2309 2310 INTEGER :: ji, jj, jk ! dummy loop indices 2311 INTEGER :: ii, ij, ik, ikmax ! local integers 2312 REAL(wp) :: zc 2313 REAL(wp) :: zN2_c ! local buoyancy difference from 10m value 2314 REAL(wp), DIMENSION(jpi,jpj) :: ztm, zsm, zLf_NH, zLf_MH 2315 REAL(wp), DIMENSION(jpi,jpj,jpts):: ztsm_midu, ztsm_midv, zabu, zabv 2316 REAL(wp), DIMENSION(jpi,jpj) :: zmld_midu, zmld_midv 2317 !!---------------------------------------------------------------------- 2318 ! 2319 ! !== MLD used for MLE ==! 2320 2321 mld_prof(:,:) = nlb10 ! Initialization to the number of w ocean point 2322 zmld(:,:) = 0._wp ! here hmlp used as a dummy variable, integrating vertically N^2 2323 zN2_c = grav * rn_osm_mle_rho_c * r1_rho0 ! convert density criteria into N^2 criteria 2324 DO_3D( 1, 1, 1, 1, nlb10, jpkm1 ) 2325 ikt = mbkt(ji,jj) 2326 zmld(ji,jj) = zmld(ji,jj) + MAX( rn2b(ji,jj,jk) , 0._wp ) * e3w(ji,jj,jk,Kmm) 2327 IF( zmld(ji,jj) < zN2_c ) mld_prof(ji,jj) = MIN( jk , ikt ) + 1 ! Mixed layer level 2328 END_3D 2329 DO_2D( 1, 1, 1, 1 ) 2330 mld_prof(ji,jj) = MAX(mld_prof(ji,jj),ibld(ji,jj)) 2331 zmld(ji,jj) = gdepw(ji,jj,mld_prof(ji,jj),Kmm) 2332 END_2D 2333 ! ensure mld_prof .ge. ibld 2334 ! 2335 ikmax = MIN( MAXVAL( mld_prof(:,:) ), jpkm1 ) ! max level of the computation 2336 ! 2337 ztm(:,:) = 0._wp 2338 zsm(:,:) = 0._wp 2339 DO_3D( 1, 1, 1, 1, 1, ikmax ) 2340 zc = e3t(ji,jj,jk,Kmm) * REAL( MIN( MAX( 0, mld_prof(ji,jj)-jk ) , 1 ) ) ! zc being 0 outside the ML t-points 2341 ztm(ji,jj) = ztm(ji,jj) + zc * ts(ji,jj,jk,jp_tem,Kmm) 2342 zsm(ji,jj) = zsm(ji,jj) + zc * ts(ji,jj,jk,jp_sal,Kmm) 2343 END_3D 2344 ! average temperature and salinity. 2345 ztm(:,:) = ztm(:,:) / MAX( e3t(:,:,1,Kmm), zmld(:,:) ) 2346 zsm(:,:) = zsm(:,:) / MAX( e3t(:,:,1,Kmm), zmld(:,:) ) 2347 ! calculate horizontal gradients at u & v points 2348 2349 DO_2D( 0, 0, 1, 0 ) 2350 zdtdx(ji,jj) = ( ztm(ji+1,jj) - ztm( ji,jj) ) * umask(ji,jj,1) / e1u(ji,jj) 2351 zdsdx(ji,jj) = ( zsm(ji+1,jj) - zsm( ji,jj) ) * umask(ji,jj,1) / e1u(ji,jj) 2352 zmld_midu(ji,jj) = 0.25_wp * (zmld(ji+1,jj) + zmld( ji,jj)) 2353 ztsm_midu(ji,jj,jp_tem) = 0.5_wp * ( ztm(ji+1,jj) + ztm( ji,jj) ) 2354 ztsm_midu(ji,jj,jp_sal) = 0.5_wp * ( zsm(ji+1,jj) + zsm( ji,jj) ) 2355 END_2D 2356 2357 DO_2D( 1, 0, 0, 0 ) 2358 zdtdy(ji,jj) = ( ztm(ji,jj+1) - ztm( ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj) 2359 zdsdy(ji,jj) = ( zsm(ji,jj+1) - zsm( ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj) 2360 zmld_midv(ji,jj) = 0.25_wp * (zmld(ji,jj+1) + zmld( ji,jj)) 2361 ztsm_midv(ji,jj,jp_tem) = 0.5_wp * ( ztm(ji,jj+1) + ztm( ji,jj) ) 2362 ztsm_midv(ji,jj,jp_sal) = 0.5_wp * ( zsm(ji,jj+1) + zsm( ji,jj) ) 2363 END_2D 2364 2365 CALL eos_rab(ztsm_midu, zmld_midu, zabu, Kmm) 2366 CALL eos_rab(ztsm_midv, zmld_midv, zabv, Kmm) 2367 2368 DO_2D( 0, 0, 1, 0 ) 2369 dbdx_mle(ji,jj) = grav*(zdtdx(ji,jj)*zabu(ji,jj,jp_tem) - zdsdx(ji,jj)*zabu(ji,jj,jp_sal)) 2370 END_2D 2371 DO_2D( 1, 0, 0, 0 ) 2372 dbdy_mle(ji,jj) = grav*(zdtdy(ji,jj)*zabv(ji,jj,jp_tem) - zdsdy(ji,jj)*zabv(ji,jj,jp_sal)) 2373 END_2D 2374 2375 DO_2D( 0, 0, 0, 0 ) 2376 ztmp = r1_ft(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf 2377 zdbds_mle(ji,jj) = SQRT( 0.5_wp * ( dbdx_mle(ji,jj) * dbdx_mle(ji,jj) + dbdy_mle(ji,jj) * dbdy_mle(ji,jj) & 2378 & + dbdx_mle(ji-1,jj) * dbdx_mle(ji-1,jj) + dbdy_mle(ji,jj-1) * dbdy_mle(ji,jj-1) ) ) 2379 END_2D 2380 2381 END SUBROUTINE zdf_osm_zmld_horizontal_gradients 2382 SUBROUTINE zdf_osm_mle_parameters( mld_prof, hmle, zhmle, zvel_mle, zdiff_mle ) 2383 !!---------------------------------------------------------------------- 2384 !! *** ROUTINE zdf_osm_mle_parameters *** 2385 !! 2386 !! ** Purpose : Timesteps the mixed layer eddy depth, hmle and calculates the mixed layer eddy fluxes for buoyancy, heat and salinity. 2387 !! 2388 !! ** Method : 2389 !! 2390 !! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008 2391 !! Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008 2392 2393 INTEGER, DIMENSION(jpi,jpj) :: mld_prof 2394 REAL(wp), DIMENSION(jpi,jpj) :: hmle, zhmle, zwb_fk, zvel_mle, zdiff_mle 2395 INTEGER :: ji, jj, jk ! dummy loop indices 2396 INTEGER :: ii, ij, ik, jkb, jkb1 ! local integers 2397 INTEGER , DIMENSION(jpi,jpj) :: inml_mle 2398 REAL(wp) :: ztmp, zdbdz, zdtdz, zdsdz, zthermal,zbeta, zbuoy, zdb_mle 2399 2400 ! Calculate vertical buoyancy, heat and salinity fluxes due to MLE. 2401 2402 DO_2D( 0, 0, 0, 0 ) 2403 IF ( lconv(ji,jj) ) THEN 2404 ztmp = r1_ft(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf 2405 ! This velocity scale, defined in Fox-Kemper et al (2008), is needed for calculating dhdt. 2406 zvel_mle(ji,jj) = zdbds_mle(ji,jj) * ztmp * hmle(ji,jj) * tmask(ji,jj,1) 2407 zdiff_mle(ji,jj) = 5.e-4_wp * rn_osm_mle_ce * ztmp * zdbds_mle(ji,jj) * zhmle(ji,jj)**2 2408 ENDIF 2409 END_2D 2410 ! Timestep mixed layer eddy depth. 2411 DO_2D( 0, 0, 0, 0 ) 2412 IF ( lmle(ji,jj) ) THEN ! MLE layer growing. 2413 ! Buoyancy gradient at base of MLE layer. 2414 zthermal = rab_n(ji,jj,1,jp_tem) 2415 zbeta = rab_n(ji,jj,1,jp_sal) 2416 jkb = mld_prof(ji,jj) 2417 jkb1 = MIN(jkb + 1, mbkt(ji,jj)) 2418 ! 2419 zbuoy = grav * ( zthermal * ts(ji,jj,mld_prof(ji,jj)+2,jp_tem,Kmm) - zbeta * ts(ji,jj,mld_prof(ji,jj)+2,jp_sal,Kmm) ) 2420 zdb_mle = zb_bl(ji,jj) - zbuoy 2421 ! Timestep hmle. 2422 hmle(ji,jj) = hmle(ji,jj) + zwb0(ji,jj) * rn_Dt / zdb_mle 2423 ELSE 2424 IF ( zhmle(ji,jj) > zhbl(ji,jj) ) THEN 2425 hmle(ji,jj) = hmle(ji,jj) - ( hmle(ji,jj) - hbl(ji,jj) ) * rn_Dt / rn_osm_mle_tau 2426 ELSE 2427 hmle(ji,jj) = hmle(ji,jj) - 10.0 * ( hmle(ji,jj) - hbl(ji,jj) ) * rn_Dt /rn_osm_mle_tau 2428 ENDIF 2429 ENDIF 2430 hmle(ji,jj) = MIN(hmle(ji,jj), ht(ji,jj)) 2431 IF(ln_osm_hmle_limit) hmle(ji,jj) = MIN(hmle(ji,jj), MAX(rn_osm_hmle_limit,1.2*hbl(ji,jj)) ) 2432 END_2D 2433 2434 mld_prof = 4 2435 DO_3D( 0, 0, 0, 0, 5, jpkm1 ) 2436 IF ( hmle(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk) 2437 END_3D 2438 DO_2D( 0, 0, 0, 0 ) 2439 zhmle(ji,jj) = gdepw(ji,jj, mld_prof(ji,jj),Kmm) 2440 END_2D 2441 END SUBROUTINE zdf_osm_mle_parameters 2442 2443 END SUBROUTINE zdf_osm 2444 2445 2446 SUBROUTINE zdf_osm_init( Kmm ) 1295 2447 !!---------------------------------------------------------------------- 1296 2448 !! *** ROUTINE zdf_osm_init *** … … 1304 2456 !! ** input : Namlist namosm 1305 2457 !!---------------------------------------------------------------------- 1306 INTEGER, INTENT(in) :: Kmm ! time level index (middle) 1307 ! 2458 INTEGER, INTENT(in) :: Kmm ! time level 1308 2459 INTEGER :: ios ! local integer 1309 2460 INTEGER :: ji, jj, jk ! dummy loop indices 2461 REAL z1_t2 1310 2462 !! 1311 2463 NAMELIST/namzdf_osm/ ln_use_osm_la, rn_osm_la, rn_osm_dstokes, nn_ave & 1312 & ,nn_osm_wave, ln_dia_osm, rn_osm_hbl0 & 1313 & ,ln_kpprimix, rn_riinfty, rn_difri, ln_convmix, rn_difconv 2464 & ,nn_osm_wave, ln_dia_osm, rn_osm_hbl0, rn_zdfosm_adjust_sd & 2465 & ,ln_kpprimix, rn_riinfty, rn_difri, ln_convmix, rn_difconv, nn_osm_wave & 2466 & ,nn_osm_SD_reduce, ln_osm_mle, rn_osm_hblfrac, rn_osm_bl_thresh, ln_zdfosm_ice_shelter 2467 ! Namelist for Fox-Kemper parametrization. 2468 NAMELIST/namosm_mle/ nn_osm_mle, rn_osm_mle_ce, rn_osm_mle_lf, rn_osm_mle_time, rn_osm_mle_lat,& 2469 & rn_osm_mle_rho_c, rn_osm_mle_thresh, rn_osm_mle_tau, ln_osm_hmle_limit, rn_osm_hmle_limit 2470 1314 2471 !!---------------------------------------------------------------------- 1315 2472 ! … … 1325 2482 WRITE(numout,*) 'zdf_osm_init : OSMOSIS Parameterisation' 1326 2483 WRITE(numout,*) '~~~~~~~~~~~~' 1327 WRITE(numout,*) ' Namelist namzdf_osm : set tke mixing parameters' 1328 WRITE(numout,*) ' Use namelist rn_osm_la ln_use_osm_la = ', ln_use_osm_la 2484 WRITE(numout,*) ' Namelist namzdf_osm : set osm mixing parameters' 2485 WRITE(numout,*) ' Use rn_osm_la ln_use_osm_la = ', ln_use_osm_la 2486 WRITE(numout,*) ' Use MLE in OBL, i.e. Fox-Kemper param ln_osm_mle = ', ln_osm_mle 1329 2487 WRITE(numout,*) ' Turbulent Langmuir number rn_osm_la = ', rn_osm_la 2488 WRITE(numout,*) ' Stokes drift reduction factor rn_zdfosm_adjust_sd = ', rn_zdfosm_adjust_sd 1330 2489 WRITE(numout,*) ' Initial hbl for 1D runs rn_osm_hbl0 = ', rn_osm_hbl0 1331 WRITE(numout,*) ' Depth scale of Stokes drift rn_osm_dstokes = ', rn_osm_dstokes2490 WRITE(numout,*) ' Depth scale of Stokes drift rn_osm_dstokes = ', rn_osm_dstokes 1332 2491 WRITE(numout,*) ' horizontal average flag nn_ave = ', nn_ave 1333 2492 WRITE(numout,*) ' Stokes drift nn_osm_wave = ', nn_osm_wave … … 1339 2498 CASE(2) 1340 2499 WRITE(numout,*) ' calculated from ECMWF wave fields' 2500 END SELECT 2501 WRITE(numout,*) ' Stokes drift reduction nn_osm_SD_reduce', nn_osm_SD_reduce 2502 WRITE(numout,*) ' fraction of hbl to average SD over/fit' 2503 WRITE(numout,*) ' exponential with nn_osm_SD_reduce = 1 or 2 rn_osm_hblfrac = ', rn_osm_hblfrac 2504 SELECT CASE (nn_osm_SD_reduce) 2505 CASE(0) 2506 WRITE(numout,*) ' No reduction' 2507 CASE(1) 2508 WRITE(numout,*) ' Average SD over upper rn_osm_hblfrac of BL' 2509 CASE(2) 2510 WRITE(numout,*) ' Fit exponential to slope rn_osm_hblfrac of BL' 1341 2511 END SELECT 2512 WRITE(numout,*) ' reduce surface SD and depth scale under ice ln_zdfosm_ice_shelter=', ln_zdfosm_ice_shelter 1342 2513 WRITE(numout,*) ' Output osm diagnostics ln_dia_osm = ', ln_dia_osm 2514 WRITE(numout,*) ' Threshold used to define BL rn_osm_bl_thresh = ', rn_osm_bl_thresh, 'm^2/s' 1343 2515 WRITE(numout,*) ' Use KPP-style shear instability mixing ln_kpprimix = ', ln_kpprimix 1344 2516 WRITE(numout,*) ' local Richardson Number limit for shear instability rn_riinfty = ', rn_riinfty … … 1359 2531 IF( zdf_osm_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_osm_init : unable to allocate arrays' ) 1360 2532 1361 call osm_rst( nit000, Kmm, 'READ' ) !* read or initialize hbl 2533 2534 IF( ln_osm_mle ) THEN 2535 ! Initialise Fox-Kemper parametrization 2536 READ ( numnam_ref, namosm_mle, IOSTAT = ios, ERR = 903) 2537 903 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namosm_mle in reference namelist') 2538 2539 READ ( numnam_cfg, namosm_mle, IOSTAT = ios, ERR = 904 ) 2540 904 IF( ios > 0 ) CALL ctl_nam ( ios , 'namosm_mle in configuration namelist') 2541 IF(lwm) WRITE ( numond, namosm_mle ) 2542 2543 IF(lwp) THEN ! Namelist print 2544 WRITE(numout,*) 2545 WRITE(numout,*) 'zdf_osm_init : initialise mixed layer eddy (MLE)' 2546 WRITE(numout,*) '~~~~~~~~~~~~~' 2547 WRITE(numout,*) ' Namelist namosm_mle : ' 2548 WRITE(numout,*) ' MLE type: =0 standard Fox-Kemper ; =1 new formulation nn_osm_mle = ', nn_osm_mle 2549 WRITE(numout,*) ' magnitude of the MLE (typical value: 0.06 to 0.08) rn_osm_mle_ce = ', rn_osm_mle_ce 2550 WRITE(numout,*) ' scale of ML front (ML radius of deformation) (nn_osm_mle=0) rn_osm_mle_lf = ', rn_osm_mle_lf, 'm' 2551 WRITE(numout,*) ' maximum time scale of MLE (nn_osm_mle=0) rn_osm_mle_time = ', rn_osm_mle_time, 's' 2552 WRITE(numout,*) ' reference latitude (degrees) of MLE coef. (nn_osm_mle=1) rn_osm_mle_lat = ', rn_osm_mle_lat, 'deg' 2553 WRITE(numout,*) ' Density difference used to define ML for FK rn_osm_mle_rho_c = ', rn_osm_mle_rho_c 2554 WRITE(numout,*) ' Threshold used to define MLE for FK rn_osm_mle_thresh = ', rn_osm_mle_thresh, 'm^2/s' 2555 WRITE(numout,*) ' Timescale for OSM-FK rn_osm_mle_tau = ', rn_osm_mle_tau, 's' 2556 WRITE(numout,*) ' switch to limit hmle ln_osm_hmle_limit = ', ln_osm_hmle_limit 2557 WRITE(numout,*) ' fraction of zmld to limit hmle to if ln_osm_hmle_limit =.T. rn_osm_hmle_limit = ', rn_osm_hmle_limit 2558 ENDIF ! 2559 ENDIF 2560 ! 2561 IF(lwp) THEN 2562 WRITE(numout,*) 2563 IF( ln_osm_mle ) THEN 2564 WRITE(numout,*) ' ==>>> Mixed Layer Eddy induced transport added to OSMOSIS BL calculation' 2565 IF( nn_osm_mle == 0 ) WRITE(numout,*) ' Fox-Kemper et al 2010 formulation' 2566 IF( nn_osm_mle == 1 ) WRITE(numout,*) ' New formulation' 2567 ELSE 2568 WRITE(numout,*) ' ==>>> Mixed Layer induced transport NOT added to OSMOSIS BL calculation' 2569 ENDIF 2570 ENDIF 2571 ! 2572 IF( ln_osm_mle ) THEN ! MLE initialisation 2573 ! 2574 rb_c = grav * rn_osm_mle_rho_c /rho0 ! Mixed Layer buoyancy criteria 2575 IF(lwp) WRITE(numout,*) 2576 IF(lwp) WRITE(numout,*) ' ML buoyancy criteria = ', rb_c, ' m/s2 ' 2577 IF(lwp) WRITE(numout,*) ' associated ML density criteria defined in zdfmxl = ', rn_osm_mle_rho_c, 'kg/m3' 2578 ! 2579 IF( nn_osm_mle == 0 ) THEN ! MLE array allocation & initialisation ! 2580 ! 2581 ELSEIF( nn_osm_mle == 1 ) THEN ! MLE array allocation & initialisation 2582 rc_f = rn_osm_mle_ce/ ( 5.e3_wp * 2._wp * omega * SIN( rad * rn_osm_mle_lat ) ) 2583 ! 2584 ENDIF 2585 ! ! 1/(f^2+tau^2)^1/2 at t-point (needed in both nn_osm_mle case) 2586 z1_t2 = 2.e-5 2587 DO_2D( 1, 1, 1, 1 ) 2588 r1_ft(ji,jj) = MIN(1./( ABS(ff_t(ji,jj)) + epsln ), ABS(ff_t(ji,jj))/z1_t2**2) 2589 END_2D 2590 ! z1_t2 = 1._wp / ( rn_osm_mle_time * rn_osm_mle_timeji,jj ) 2591 ! r1_ft(:,:) = 1._wp / SQRT( ff_t(:,:) * ff_t(:,:) + z1_t2 ) 2592 ! 2593 ENDIF 2594 2595 call osm_rst( nit000, Kmm, 'READ' ) !* read or initialize hbl, dh, hmle 2596 1362 2597 1363 2598 IF( ln_zdfddm) THEN … … 1440 2675 CALL iom_set_rstw_var_active('wn') 1441 2676 CALL iom_set_rstw_var_active('hbl') 1442 CALL iom_set_rstw_var_active('hbli') 2677 CALL iom_set_rstw_var_active('dh') 2678 IF( ln_osm_mle ) THEN 2679 CALL iom_set_rstw_var_active('hmle') 2680 END IF 1443 2681 ENDIF 1444 2682 END SUBROUTINE zdf_osm_init … … 1459 2697 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 1460 2698 1461 INTEGER :: id1, id2 ! iom enquiry index2699 INTEGER :: id1, id2, id3 ! iom enquiry index 1462 2700 INTEGER :: ji, jj, jk ! dummy loop indices 1463 2701 INTEGER :: iiki, ikt ! local integer … … 1465 2703 REAL(wp) :: zN2_c ! local scalar 1466 2704 REAL(wp) :: rho_c = 0.01_wp !: density criterion for mixed layer depth 1467 INTEGER, DIMENSION( :,:), ALLOCATABLE:: imld_rst ! level of mixed-layer depth (pycnocline top)2705 INTEGER, DIMENSION(jpi,jpj) :: imld_rst ! level of mixed-layer depth (pycnocline top) 1468 2706 !!---------------------------------------------------------------------- 1469 2707 ! … … 1480 2718 WRITE(numout,*) ' ===>>>> : ww not in restart file, set to zero initially' 1481 2719 END IF 2720 1482 2721 id1 = iom_varid( numror, 'hbl' , ldstop = .FALSE. ) 1483 id2 = iom_varid( numror, ' hbli' , ldstop = .FALSE. )2722 id2 = iom_varid( numror, 'dh' , ldstop = .FALSE. ) 1484 2723 IF( id1 > 0 .AND. id2 > 0) THEN ! 'hbl' exists; read and return 1485 2724 CALL iom_get( numror, jpdom_auto, 'hbl' , hbl , ldxios = lrxios ) 1486 CALL iom_get( numror, jpdom_auto, 'hbli', hbli, ldxios = lrxios ) 1487 WRITE(numout,*) ' ===>>>> : hbl & hbli read from restart file' 2725 CALL iom_get( numror, jpdom_auto, 'dh', dh, ldxios = lrxios ) 2726 WRITE(numout,*) ' ===>>>> : hbl & dh read from restart file' 2727 IF( ln_osm_mle ) THEN 2728 id3 = iom_varid( numror, 'hmle' , ldstop = .FALSE. ) 2729 IF( id3 > 0) THEN 2730 CALL iom_get( numror, jpdom_auto, 'hmle' , hmle , ldxios = lrxios ) 2731 WRITE(numout,*) ' ===>>>> : hmle read from restart file' 2732 ELSE 2733 WRITE(numout,*) ' ===>>>> : hmle not found, set to hbl' 2734 hmle(:,:) = hbl(:,:) ! Initialise MLE depth. 2735 END IF 2736 END IF 1488 2737 RETURN 1489 ELSE ! 'hbl' & ' hbli' not in restart file, recalculate2738 ELSE ! 'hbl' & 'dh' not in restart file, recalculate 1490 2739 WRITE(numout,*) ' ===>>>> : previous run without osmosis scheme, hbl computed from stratification' 1491 2740 END IF … … 1495 2744 ! If READ/WRITE Flag is 'WRITE', write hbl into the restart file, then return 1496 2745 !!----------------------------------------------------------------------------- 1497 IF( TRIM(cdrw) == 'WRITE') THEN !* Write hbl iinto the restart file, then return2746 IF( TRIM(cdrw) == 'WRITE') THEN !* Write hbl into the restart file, then return 1498 2747 IF(lwp) WRITE(numout,*) '---- osm-rst ----' 1499 CALL iom_rstput( kt, nitrst, numrow, 'wn' , ww , ldxios = lwxios ) 1500 CALL iom_rstput( kt, nitrst, numrow, 'hbl' , hbl , ldxios = lwxios ) 1501 CALL iom_rstput( kt, nitrst, numrow, 'hbli' , hbli, ldxios = lwxios ) 2748 CALL iom_rstput( kt, nitrst, numrow, 'ww' , ww, ldxios = lwxios ) 2749 CALL iom_rstput( kt, nitrst, numrow, 'hbl' , hbl, ldxios = lwxios ) 2750 CALL iom_rstput( kt, nitrst, numrow, 'dh' , dh, ldxios = lwxios ) 2751 IF( ln_osm_mle ) THEN 2752 CALL iom_rstput( kt, nitrst, numrow, 'hmle', hmle, ldxios = lwxios ) 2753 END IF 1502 2754 RETURN 1503 2755 END IF … … 1507 2759 !!----------------------------------------------------------------------------- 1508 2760 IF( lwp ) WRITE(numout,*) ' ===>>>> : calculating hbl computed from stratification' 1509 ALLOCATE( imld_rst(jpi,jpj) )1510 2761 ! w-level of the mixing and mixed layers 1511 2762 CALL eos_rab( ts(:,:,:,:,Kmm), rab_n, Kmm ) … … 1516 2767 ! 1517 2768 hbl(:,:) = 0._wp ! here hbl used as a dummy variable, integrating vertically N^2 1518 DO_3D( 1, 1, 1, 1, 1, jpkm1 ) ! Mixed layer level: w-level2769 DO_3D( 1, 1, 1, 1, 1, jpkm1 ) 1519 2770 ikt = mbkt(ji,jj) 1520 2771 hbl(ji,jj) = hbl(ji,jj) + MAX( rn2(ji,jj,jk) , 0._wp ) * e3w(ji,jj,jk,Kmm) … … 1523 2774 ! 1524 2775 DO_2D( 1, 1, 1, 1 ) 1525 iiki = imld_rst(ji,jj) 1526 hbl (ji,jj) = gdepw(ji,jj,iiki ,Kmm) * ssmask(ji,jj) ! Turbocline depth 2776 iiki = MAX(4,imld_rst(ji,jj)) 2777 hbl (ji,jj) = gdepw(ji,jj,iiki,Kmm ) ! Turbocline depth 2778 dh (ji,jj) = e3t(ji,jj,iiki-1,Kmm ) ! Turbocline depth 1527 2779 END_2D 1528 hbl = MAX(hbl,epsln) 1529 hbli(:,:) = hbl(:,:) 1530 DEALLOCATE( imld_rst ) 2780 1531 2781 WRITE(numout,*) ' ===>>>> : hbl computed from stratification' 2782 2783 IF( ln_osm_mle ) THEN 2784 hmle(:,:) = hbl(:,:) ! Initialise MLE depth. 2785 WRITE(numout,*) ' ===>>>> : hmle set = to hbl' 2786 END IF 2787 2788 ww(:,:,:) = 0._wp 2789 WRITE(numout,*) ' ===>>>> : wn not in restart file, set to zero initially' 1532 2790 END SUBROUTINE osm_rst 1533 2791 … … 1560 2818 ENDIF 1561 2819 1562 ! add non-local temperature and salinity flux1563 2820 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 1564 2821 pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) & … … 1570 2827 END_3D 1571 2828 1572 1573 ! save the non-local tracer flux trends for diagnostic 2829 ! save the non-local tracer flux trends for diagnostics 1574 2830 IF( l_trdtra ) THEN 1575 2831 ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:) 1576 2832 ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) - ztrds(:,:,:) 1577 !!bug gm jpttdzdf ==> jpttosm 1578 CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_ zdf, ztrdt )1579 CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_ zdf, ztrds )2833 2834 CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_osm, ztrdt ) 2835 CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_osm, ztrds ) 1580 2836 DEALLOCATE( ztrdt ) ; DEALLOCATE( ztrds ) 1581 2837 ENDIF … … 1643 2899 1644 2900 !!====================================================================== 2901 1645 2902 END MODULE zdfosm
Note: See TracChangeset
for help on using the changeset viewer.
