- Timestamp:
- 2010-11-16T14:03:39+01:00 (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90
r2395 r2397 5 5 !! turbulent closure parameterization 6 6 !!====================================================================== 7 !! History : 3.0 ! 2009-09(G. Reffray) Original code8 !! 3.3 ! 2010-10 (C. Bricaud) add in the reference7 !! History : 3.0 ! 2009-09 (G. Reffray) Original code 8 !! 3.3 ! 2010-10 (C. Bricaud) Add in the reference 9 9 !!---------------------------------------------------------------------- 10 10 #if defined key_zdfgls || defined key_esopa … … 34 34 35 35 PUBLIC zdf_gls ! routine called in step module 36 PUBLIC zdf_gls_init ! routine called in stepmodule36 PUBLIC zdf_gls_init ! routine called in opa module 37 37 PUBLIC gls_rst ! routine called in step module 38 38 … … 46 46 ! !!! ** Namelist namzdf_gls ** 47 47 LOGICAL :: ln_crban = .FALSE. ! =T use Craig and Banner scheme 48 LOGICAL :: ln_length_lim = .FALSE. ! use limit on the dissipation rate under stable stratif. (Galperin et al., 1988) 49 LOGICAL :: ln_sigpsi = .FALSE. ! Activate Burchard (2003) modification for k-eps closure AND wave breaking mixing 50 REAL(wp) :: rn_epsmin = 1.e-12_wp ! minimum value of dissipation (m2/s3) 51 REAL(wp) :: rn_emin = 1.e-6_wp ! minimum value of TKE (m2/s2) 48 LOGICAL :: ln_length_lim = .FALSE. ! use limit on the dissipation rate under stable stratification (Galperin et al. 1988) 49 LOGICAL :: ln_sigpsi = .FALSE. ! Activate Burchard (2003) modification for k-eps closure & wave breaking mixing 52 50 INTEGER :: nn_tkebc_surf = 0 ! TKE surface boundary condition (=0/1) 53 51 INTEGER :: nn_tkebc_bot = 0 ! TKE bottom boundary condition (=0/1) … … 57 55 INTEGER :: nn_clos = 0 ! closure 0/1/2/3 MY82/k-eps/k-w/gen 58 56 REAL(wp) :: rn_clim_galp = 0.53_wp ! Holt 2008 value for k-eps: 0.267 59 REAL(wp) :: hsro = 0.003_wp ! Minimum surface roughness 57 REAL(wp) :: rn_epsmin = 1.e-12_wp ! minimum value of dissipation (m2/s3) 58 REAL(wp) :: rn_emin = 1.e-6_wp ! minimum value of TKE (m2/s2) 60 59 REAL(wp) :: rn_charn = 2.e+5_wp ! Charnock constant for surface breaking waves mixing : 1400. (standard) or 2.e5 (Stacey value) 61 60 REAL(wp) :: rn_crban = 100._wp ! Craig and Banner constant for surface breaking waves mixing 62 61 63 REAL(wp) :: rcm_sf = 0.73_wp ! Shear free turbulence parameters 64 REAL(wp) :: ra_sf = -2.0_wp ! Must be negative -2 < ra_sf < -1 65 REAL(wp) :: rl_sf = 0.2_wp ! 0 <rl_sf<vkarmn 66 REAL(wp) :: rghmin = -0.28 67 REAL(wp) :: rgh0 = 0.0329 68 REAL(wp) :: rghcri = 0.03 62 REAL(wp) :: hsro = 0.003_wp ! Minimum surface roughness 63 REAL(wp) :: rcm_sf = 0.73_wp ! Shear free turbulence parameters 64 REAL(wp) :: ra_sf = -2.0_wp ! Must be negative -2 < ra_sf < -1 65 REAL(wp) :: rl_sf = 0.2_wp ! 0 <rl_sf<vkarmn 66 REAL(wp) :: rghmin = -0.28_wp 67 REAL(wp) :: rgh0 = 0.0329_wp 68 REAL(wp) :: rghcri = 0.03_wp 69 69 REAL(wp) :: ra1 = 0.92_wp 70 70 REAL(wp) :: ra2 = 0.74_wp … … 88 88 REAL(wp) :: rm7 = 0.0_wp 89 89 REAL(wp) :: rm8 = 0.318_wp 90 REAL(wp) :: rc02 91 REAL(wp) :: rc02r 92 REAL(wp) :: rc03 93 REAL(wp) :: rc04 94 REAL(wp) :: rc03_sqrt2_galp 95 REAL(wp) :: rsbc_mb 96 REAL(wp) :: rsbc_std 97 REAL(wp) :: rsbc_tke1 98 REAL(wp) :: rsbc_tke2 99 REAL(wp) :: rsbc_tke3 100 REAL(wp) :: rsbc_psi1 101 REAL(wp) :: rsbc_psi2 102 REAL(wp) :: rsbc_psi3 103 REAL(wp) :: rsbc_zs 104 REAL(wp) :: rfact_tke, rfact_psi 105 REAL(wp) :: rc0, rc2, rc3, rf6, rcff, rc_diff 106 REAL(wp) :: rs0, rs1, rs2, rs4, rs5, rs6 107 REAL(wp) :: rd0, rd1, rd2, rd3, rd4, rd5 108 REAL(wp) :: rsc_tke, rsc_psi, rpsi1, rpsi2, rpsi3, rsc_psi0 109 REAL(wp) :: rpsi3m, rpsi3p, rpp, rmm, rnn 90 91 REAL(wp) :: rc02, rc02r, rc03, rc04 ! coefficients deduced from above parameters 92 REAL(wp) :: rc03_sqrt2_galp ! - - - - 93 REAL(wp) :: rsbc_tke1, rsbc_tke2, rsbc_tke3, rfact_tke ! - - - - 94 REAL(wp) :: rsbc_psi1, rsbc_psi2, rsbc_psi3, rfact_psi ! - - - - 95 REAL(wp) :: rsbc_mb , rsbc_std , rsbc_zs ! - - - - 96 REAL(wp) :: rc0, rc2, rc3, rf6, rcff, rc_diff ! - - - - 97 REAL(wp) :: rs0, rs1, rs2, rs4, rs5, rs6 ! - - - - 98 REAL(wp) :: rd0, rd1, rd2, rd3, rd4, rd5 ! - - - - 99 REAL(wp) :: rsc_tke, rsc_psi, rpsi1, rpsi2, rpsi3, rsc_psi0 ! - - - - 100 REAL(wp) :: rpsi3m, rpsi3p, rpp, rmm, rnn ! - - - - 110 101 111 102 !! * Substitutions … … 124 115 !! 125 116 !! ** Purpose : Compute the vertical eddy viscosity and diffusivity 126 !! coefficients using a 2.5turbulent closure scheme.117 !! coefficients using the GLS turbulent closure scheme. 127 118 !!---------------------------------------------------------------------- 128 119 USE oce, z_elem_a => ua ! use ua as workspace … … 133 124 INTEGER, INTENT(in) :: kt ! ocean time step 134 125 INTEGER :: ji, jj, jk, ibot, ibotm1, dir ! dummy loop arguments 135 ! 136 REAL(wp) :: zesh2, zsigpsi ! temporary scalars 137 REAL(wp) :: ztx2, zty2, zup, zdown, zcof ! - 138 REAL(wp) :: zratio, zrn2, zflxb, sh ! - - 139 REAL(wp) :: prod, buoy, diss, zdiss, sm ! - - 140 REAL(wp) :: gh, gm, shr, dif, zsqen, zav ! - - 141 REAL(wp), DIMENSION(jpi,jpj) :: zdep ! 142 REAL(wp), DIMENSION(jpi,jpj) :: zflxs ! Turbulence fluxed induced by internal waves 143 REAL(wp), DIMENSION(jpi,jpj) :: zhsro ! Surface roughness (surface waves) 144 REAL(wp), DIMENSION(jpi,jpj,jpk) :: eb ! tke at time before 145 REAL(wp), DIMENSION(jpi,jpj,jpk) :: mxlb ! mixing length at time before 146 REAL(wp), DIMENSION(jpi,jpj,jpk) :: shear ! vertical shear 147 REAL(wp), DIMENSION(jpi,jpj,jpk) :: eps ! dissipation rate 148 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwall_psi ! Wall function for psi schmidt number in the wb case 149 ! (if ln_sigpsi.AND.ln_crban) 126 REAL(wp) :: zesh2, zsigpsi, zcoef, zex1, zex2 ! local scalars 127 REAL(wp) :: ztx2, zty2, zup, zdown, zcof ! - - 128 REAL(wp) :: zratio, zrn2, zflxb, sh ! - - 129 REAL(wp) :: prod, buoy, diss, zdiss, sm ! - - 130 REAL(wp) :: gh, gm, shr, dif, zsqen, zav ! - - 131 REAL(wp), DIMENSION(jpi,jpj) :: zdep ! 132 REAL(wp), DIMENSION(jpi,jpj) :: zflxs ! Turbulence fluxed induced by internal waves 133 REAL(wp), DIMENSION(jpi,jpj) :: zhsro ! Surface roughness (surface waves) 134 REAL(wp), DIMENSION(jpi,jpj,jpk) :: eb ! tke at time before 135 REAL(wp), DIMENSION(jpi,jpj,jpk) :: mxlb ! mixing length at time before 136 REAL(wp), DIMENSION(jpi,jpj,jpk) :: shear ! vertical shear 137 REAL(wp), DIMENSION(jpi,jpj,jpk) :: eps ! dissipation rate 138 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwall_psi ! Wall function use in the wb case (ln_sigpsi.AND.ln_crban=T) 150 139 !!-------------------------------------------------------------------- 151 140 152 141 ! Preliminary computing 153 142 154 ustars2 = 0. ; ustarb2 = 0. ; psi = 0. ; zwall_psi = 0.143 ustars2 = 0._wp ; ustarb2 = 0._wp ; psi = 0._wp ; zwall_psi = 0._wp 155 144 156 145 ! Compute wind stress at T-points … … 162 151 ! wind stress 163 152 ! squared surface velocity scale 164 ztx2 = 2. * (utau(ji-1,jj )*umask(ji-1,jj,1) + utau(ji,jj)*umask(ji,jj,1)) / & 165 & MAX(1., umask(ji-1,jj,1) + umask(ji,jj,1)) 166 zty2 = 2. * (vtau(ji ,jj-1)*vmask(ji,jj-1,1) + vtau(ji,jj)*vmask(ji,jj,1)) / & 167 & MAX(1., vmask(ji,jj-1,1) + vmask(ji,jj,1)) 168 ustars2(ji,jj) = rsbc_tke2 * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) 153 ustars2(ji,jj) = rau0r * taum(ji,jj) * tmask(ji,jj,1) 169 154 ! 170 155 ! bottom friction 171 ztx2 = 2. * (wbotu(ji-1,jj )*umask(ji-1,jj,1) + wbotu(ji,jj)*umask(ji,jj,1)) /&172 & MAX(1., umask(ji-1,jj,1) + umask(ji,jj,1))173 zty2 = 2. * (wbotv(ji ,jj-1)*vmask(ji,jj-1,1) + wbotv(ji,jj)*vmask(ji,jj,1)) /&174 & MAX(1., vmask(ji,jj-1,1) + vmask(ji,jj,1))175 ustarb2(ji,jj) = 0.5 *SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)176 END DO177 END DO156 ztx2 = ( wbotu(ji-1,jj)*umask(ji-1,jj,1) + wbotu(ji,jj)*umask(ji,jj,1) ) & 157 & / MAX( 1., umask(ji-1,jj,1) + umask(ji,jj,1) ) 158 zty2 = ( wbotv(ji,jj-1)*vmask(ji,jj-1,1) + wbotv(ji,jj)*vmask(ji,jj,1) ) & 159 & / MAX( 1., vmask(ji,jj-1,1) + vmask(ji,jj,1) ) 160 ustarb2(ji,jj) = SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) 161 END DO 162 END DO 178 163 179 164 ! In case of breaking surface waves mixing, 180 165 ! Compute surface roughness length according to Charnock formula: 181 IF (ln_crban) THEN 182 zhsro(:,:) = MAX(rsbc_zs * ustars2(:,:), hsro) 183 ELSE 184 zhsro(:,:) = hsro 166 IF( ln_crban ) THEN ; zhsro(:,:) = MAX(rsbc_zs * ustars2(:,:), hsro) 167 ELSE ; zhsro(:,:) = hsro 185 168 ENDIF 186 169 … … 198 181 & * fse3vw_b(ji,jj,jk) ) 199 182 eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT(en(ji,jj,jk)) / mxln(ji,jj,jk) 200 END DO201 END DO202 END DO183 END DO 184 END DO 185 END DO 203 186 ! 204 187 ! Lateral boundary conditions (avmu,avmv) (sign unchanged) 205 CALL lbc_lnk( avmu, 'U', 1. ) ; 188 CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. ) 206 189 207 190 ! Save tke at before time step … … 209 192 mxlb(:,:,:) = mxln(:,:,:) 210 193 211 IF ( nn_clos .EQ. 0 ) THEN! Mellor-Yamada194 IF( nn_clos == 0 ) THEN ! Mellor-Yamada 212 195 DO jk = 2, jpkm1 213 196 DO jj = 2, jpjm1 214 197 DO ji = fs_2, fs_jpim1 ! vector opt. 215 198 zup = mxln(ji,jj,jk) * fsdepw(ji,jj,mbathy(ji,jj)) 216 zdown = vkarmn * fsdepw(ji,jj,jk) * (-fsdepw(ji,jj,jk)+fsdepw(ji,jj,mbathy(ji,jj))) 217 zwall (ji,jj,jk) = ( 1. + re2 * ( zup / MAX( zdown, rsmall ) )**2. ) * tmask(ji,jj,jk) 218 ENDDO 219 ENDDO 220 ENDDO 199 zdown = vkarmn * fsdepw(ji,jj,jk) * ( -fsdepw(ji,jj,jk) + fsdepw(ji,jj,mbathy(ji,jj)) ) 200 zcoef = ( zup / MAX( zdown, rsmall ) ) 201 zwall (ji,jj,jk) = ( 1._wp + re2 * zcoef*zcoef ) * tmask(ji,jj,jk) 202 END DO 203 END DO 204 END DO 221 205 ENDIF 222 206 … … 248 232 diss = eps(ji,jj,jk) 249 233 ! 250 dir = 0.5 + sign(0.5,shear(ji,jj,jk)+buoy) ! dir =1(=0) if shear(ji,jj,jk)+buoy >0(<0)251 ! 252 zesh2 = dir*(shear(ji,jj,jk)+buoy)+(1 -dir)*shear(ji,jj,jk) ! production term253 zdiss = dir*(diss/en(ji,jj,jk)) +(1 -dir)*(diss-buoy)/en(ji,jj,jk) ! dissipation term234 dir = 0.5_wp + SIGN( 0.5_wp, shear(ji,jj,jk) + buoy ) ! dir =1(=0) if shear(ji,jj,jk)+buoy >0(<0) 235 ! 236 zesh2 = dir*(shear(ji,jj,jk)+buoy)+(1._wp-dir)*shear(ji,jj,jk) ! production term 237 zdiss = dir*(diss/en(ji,jj,jk)) +(1._wp-dir)*(diss-buoy)/en(ji,jj,jk) ! dissipation term 254 238 ! 255 239 ! Compute a wall function from 1. to rsc_psi*zwall/rsc_psi0 … … 257 241 ! there is no need to set a boundary condition for zwall_psi at the top and bottom boundaries. 258 242 ! Otherwise, this should be rsc_psi/rsc_psi0 259 IF (ln_sigpsi) THEN260 zsigpsi = MIN(1., zesh2/eps(ji,jj,jk))! 0. <= zsigpsi <= 1.261 zwall_psi(ji,jj,jk) = rsc_psi / (zsigpsi * rsc_psi + (1.-zsigpsi) * rsc_psi0 / MAX(zwall(ji,jj,jk),1.))243 IF( ln_sigpsi ) THEN 244 zsigpsi = MIN( 1._wp, zesh2 / eps(ji,jj,jk) ) ! 0. <= zsigpsi <= 1. 245 zwall_psi(ji,jj,jk) = rsc_psi / ( zsigpsi * rsc_psi + (1._wp-zsigpsi) * rsc_psi0 / MAX( zwall(ji,jj,jk), 1._wp ) ) 262 246 ELSE 263 zwall_psi(ji,jj,jk) = 1.247 zwall_psi(ji,jj,jk) = 1._wp 264 248 ENDIF 265 249 ! … … 276 260 ! 277 261 ! diagonal 278 z_elem_b(ji,jj,jk) = 1. - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) &279 & + rdt * zdiss * tmask(ji,jj,jk)262 z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) & 263 & + rdt * zdiss * tmask(ji,jj,jk) 280 264 ! 281 265 ! right hand side in en … … 285 269 END DO 286 270 ! 287 z_elem_b(:,:,jpk) = 1. 271 z_elem_b(:,:,jpk) = 1._wp 288 272 ! 289 273 ! Set surface condition on zwall_psi (1 at the bottom) 290 IF (ln_sigpsi) THEN 291 DO jj = 2, jpjm1 292 DO ji = fs_2, fs_jpim1 ! vector opt. 293 zwall_psi(ji,jj,1) = rsc_psi/rsc_psi0 294 END DO 295 END DO 274 IF( ln_sigpsi ) THEN 275 zcoef = rsc_psi / rsc_psi0 276 DO jj = 2, jpjm1 277 DO ji = fs_2, fs_jpim1 ! vector opt. 278 zwall_psi(ji,jj,1) = zcoef 279 END DO 280 END DO 296 281 ENDIF 297 282 … … 302 287 ! 303 288 CASE ( 0 ) ! Dirichlet case 304 ! 305 IF (ln_crban) THEN ! Wave induced mixing case 306 ! ! en(1) = q2(1) = 0.5 * (15.8 * Ccb)^(2/3) * u*^2 307 ! ! balance between the production and the dissipation terms including the wave effect 308 ! 309 en(:,:,1) = MAX( rsbc_tke1 * ustars2(:,:), rn_emin ) 310 z_elem_a(:,:,1) = en(:,:,1) 311 z_elem_c(:,:,1) = 0. 312 z_elem_b(:,:,1) = 1. 313 ! 314 ! one level below 315 en(:,:,2) = MAX( rsbc_tke1 * ustars2(:,:) * ( (zhsro(:,:)+fsdepw(:,:,2))/zhsro(:,:) )**ra_sf, rn_emin ) 316 z_elem_a(:,:,2) = 0. 317 z_elem_c(:,:,2) = 0. 318 z_elem_b(:,:,2) = 1. 319 ! 320 ELSE ! No wave induced mixing case 321 ! ! en(1) = u*^2/C0^2 & l(1) = K*zs 322 ! ! balance between the production and the dissipation terms 323 ! 324 en(:,:,1) = MAX( rc02r * ustars2(:,:), rn_emin ) 325 z_elem_a(:,:,1) = en(:,:,1) 326 z_elem_c(:,:,1) = 0. 327 z_elem_b(:,:,1) = 1. 328 ! 329 ! one level below 330 en(:,:,2) = MAX( rc02r * ustars2(:,:), rn_emin ) 331 z_elem_a(:,:,2) = 0. 332 z_elem_c(:,:,2) = 0. 333 z_elem_b(:,:,2) = 1. 334 ! 335 ENDIF 336 ! 289 ! 290 IF (ln_crban) THEN ! Wave induced mixing case 291 ! ! en(1) = q2(1) = 0.5 * (15.8 * Ccb)^(2/3) * u*^2 292 ! ! balance between the production and the dissipation terms including the wave effect 293 en(:,:,1) = MAX( rsbc_tke1 * ustars2(:,:), rn_emin ) 294 z_elem_a(:,:,1) = en(:,:,1) 295 z_elem_c(:,:,1) = 0._wp 296 z_elem_b(:,:,1) = 1._wp 297 ! 298 ! one level below 299 en(:,:,2) = MAX( rsbc_tke1 * ustars2(:,:) * ( (zhsro(:,:)+fsdepw(:,:,2))/zhsro(:,:) )**ra_sf, rn_emin ) 300 z_elem_a(:,:,2) = 0._wp 301 z_elem_c(:,:,2) = 0._wp 302 z_elem_b(:,:,2) = 1._wp 303 ! 304 ELSE ! No wave induced mixing case 305 ! ! en(1) = u*^2/C0^2 & l(1) = K*zs 306 ! ! balance between the production and the dissipation terms 307 en(:,:,1) = MAX( rc02r * ustars2(:,:), rn_emin ) 308 z_elem_a(:,:,1) = en(:,:,1) 309 z_elem_c(:,:,1) = 0._wp 310 z_elem_b(:,:,1) = 1._wp 311 ! 312 ! one level below 313 en(:,:,2) = MAX( rc02r * ustars2(:,:), rn_emin ) 314 z_elem_a(:,:,2) = 0._wp 315 z_elem_c(:,:,2) = 0._wp 316 z_elem_b(:,:,2) = 1._wp 317 ! 318 ENDIF 319 ! 337 320 CASE ( 1 ) ! Neumann boundary condition on d(e)/dz 338 !339 IF (ln_crban) THEN ! Shear free case: d(e)/dz= Fw340 !341 ! Dirichlet conditions at k=1 (Cosmetic)342 en(:,:,1) = MAX( rsbc_tke1 * ustars2(:,:), rn_emin )343 z_elem_a(:,:,1) = en(:,:,1)344 z_elem_c(:,:,1) = 0.345 z_elem_b(:,:,1) = 1.346 ! at k=2, set de/dz=Fw347 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b348 z_elem_a(:,:,2) = 0.349 zflxs(:,:) = rsbc_tke3 * ustars2(:,:)**1.5 * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5*ra_sf)350 en(:,:,2) = en(:,:,2) + zflxs(:,:)/fse3w(:,:,2)351 !352 ELSE ! No wave induced mixing case: d(e)/dz=0.353 !354 ! Dirichlet conditions at k=1 (Cosmetic)355 en(:,:,1) = MAX( rc02r * ustars2(:,:), rn_emin )356 z_elem_a(:,:,1) = en(:,:,1)357 z_elem_c(:,:,1) = 0.358 z_elem_b(:,:,1) = 1.359 ! at k=2 set de/dz=0.:360 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2)! Remove z_elem_a from z_elem_b361 z_elem_a(:,:,2) = 0.362 !363 ENDIF364 !321 ! 322 IF (ln_crban) THEN ! Shear free case: d(e)/dz= Fw 323 ! 324 ! Dirichlet conditions at k=1 (Cosmetic) 325 en(:,:,1) = MAX( rsbc_tke1 * ustars2(:,:), rn_emin ) 326 z_elem_a(:,:,1) = en(:,:,1) 327 z_elem_c(:,:,1) = 0._wp 328 z_elem_b(:,:,1) = 1._wp 329 ! at k=2, set de/dz=Fw 330 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 331 z_elem_a(:,:,2) = 0._wp 332 zflxs(:,:) = rsbc_tke3 * ustars2(:,:)**1.5_wp * ( (zhsro(:,:)+fsdept(:,:,1) ) / zhsro(:,:) )**(1.5*ra_sf) 333 en(:,:,2) = en(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 334 ! 335 ELSE ! No wave induced mixing case: d(e)/dz=0. 336 ! 337 ! Dirichlet conditions at k=1 (Cosmetic) 338 en(:,:,1) = MAX( rc02r * ustars2(:,:), rn_emin ) 339 z_elem_a(:,:,1) = en(:,:,1) 340 z_elem_c(:,:,1) = 0._wp 341 z_elem_b(:,:,1) = 1._wp 342 ! at k=2 set de/dz=0.: 343 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 344 z_elem_a(:,:,2) = 0._wp 345 ! 346 ENDIF 347 ! 365 348 END SELECT 366 349 … … 371 354 ! 372 355 CASE ( 0 ) ! Dirichlet 373 ! ! en(ibot) = u*^2 / Co2 and mxln(ibot) = rn_lmin 374 ! ! Balance between the production and the dissipation terms 375 ! 356 ! ! en(ibot) = u*^2 / Co2 and mxln(ibot) = rn_lmin 357 ! ! Balance between the production and the dissipation terms 376 358 !CDIR NOVERRCHK 377 DO jj = 2, jpjm1359 DO jj = 2, jpjm1 378 360 !CDIR NOVERRCHK 379 DO ji = fs_2, fs_jpim1 ! vector opt.380 ibot = mbathy(ji,jj)381 ibotm1 = ibot-1382 !383 ! Bottom level Dirichlet condition:384 z_elem_a(ji,jj,ibot ) = 0.385 z_elem_c(ji,jj,ibot ) = 0.386 z_elem_b(ji,jj,ibot ) = 1.387 en(ji,jj,ibot ) = MAX( rc02r * ustarb2(ji,jj), rn_emin )388 !389 ! Just above last level, Dirichlet condition again390 z_elem_a(ji,jj,ibotm1) = 0.391 z_elem_c(ji,jj,ibotm1) = 0.392 z_elem_b(ji,jj,ibotm1) = 1.393 en(ji,jj,ibotm1) = MAX( rc02r * ustarb2(ji,jj), rn_emin )394 END DO395 END DO396 !361 DO ji = fs_2, fs_jpim1 ! vector opt. 362 ibot = mbathy(ji,jj) 363 ibotm1 = ibot-1 364 ! 365 ! Bottom level Dirichlet condition: 366 z_elem_a(ji,jj,ibot ) = 0._wp 367 z_elem_c(ji,jj,ibot ) = 0._wp 368 z_elem_b(ji,jj,ibot ) = 1._wp 369 en(ji,jj,ibot ) = MAX( rc02r * ustarb2(ji,jj), rn_emin ) 370 ! 371 ! Just above last level, Dirichlet condition again 372 z_elem_a(ji,jj,ibotm1) = 0._wp 373 z_elem_c(ji,jj,ibotm1) = 0._wp 374 z_elem_b(ji,jj,ibotm1) = 1._wp 375 en(ji,jj,ibotm1) = MAX( rc02r * ustarb2(ji,jj), rn_emin ) 376 END DO 377 END DO 378 ! 397 379 CASE ( 1 ) ! Neumman boundary condition 398 !380 ! 399 381 !CDIR NOVERRCHK 400 DO jj = 2, jpjm1382 DO jj = 2, jpjm1 401 383 !CDIR NOVERRCHK 402 DO ji = fs_2, fs_jpim1 ! vector opt.403 ibot = mbathy(ji,jj)404 ibotm1 = ibot-1405 !406 ! Bottom level Dirichlet condition:407 z_elem_a(ji,jj,ibot) = 0.408 z_elem_c(ji,jj,ibot) = 0.409 z_elem_b(ji,jj,ibot) = 1.410 en(ji,jj,ibot) = MAX( rc02r * ustarb2(ji,jj), rn_emin )411 !412 ! Just above last level: Neumann condition413 z_elem_b(ji,jj,ibotm1) = z_elem_b(ji,jj,ibotm1) + z_elem_c(ji,jj,ibotm1)! Remove z_elem_c from z_elem_b414 z_elem_c(ji,jj,ibotm1) = 0.415 END DO416 END DO417 !384 DO ji = fs_2, fs_jpim1 ! vector opt. 385 ibot = mbathy(ji,jj) 386 ibotm1 = ibot-1 387 ! 388 ! Bottom level Dirichlet condition: 389 z_elem_a(ji,jj,ibot) = 0._wp 390 z_elem_c(ji,jj,ibot) = 0._wp 391 z_elem_b(ji,jj,ibot) = 1._wp 392 en(ji,jj,ibot) = MAX( rc02r * ustarb2(ji,jj), rn_emin ) 393 ! 394 ! Just above last level: Neumann condition 395 z_elem_b(ji,jj,ibotm1) = z_elem_b(ji,jj,ibotm1) + z_elem_c(ji,jj,ibotm1) ! Remove z_elem_c from z_elem_b 396 z_elem_c(ji,jj,ibotm1) = 0._wp 397 END DO 398 END DO 399 ! 418 400 END SELECT 419 401 … … 442 424 END DO 443 425 END DO 444 ! 445 ! set the minimum value of tke 426 ! ! set the minimum value of tke 446 427 en(:,:,:) = MAX( en(:,:,:), rn_emin ) 447 428 … … 455 436 ! 456 437 CASE( 0 ) ! k-kl (Mellor-Yamada) 457 ! 458 DO jk = 2, jpkm1 459 DO jj = 2, jpjm1 460 DO ji = fs_2, fs_jpim1 ! vector opt. 461 psi(ji,jj,jk) = en(ji,jj,jk) * mxln(ji,jj,jk) 462 ENDDO 463 ENDDO 464 ENDDO 465 ! 438 DO jk = 2, jpkm1 439 DO jj = 2, jpjm1 440 DO ji = fs_2, fs_jpim1 ! vector opt. 441 psi(ji,jj,jk) = en(ji,jj,jk) * mxln(ji,jj,jk) 442 END DO 443 END DO 444 END DO 445 ! 466 446 CASE( 1 ) ! k-eps 467 ! 468 DO jk = 2, jpkm1 469 DO jj = 2, jpjm1 470 DO ji = fs_2, fs_jpim1 ! vector opt. 471 psi(ji,jj,jk) = eps(ji,jj,jk) 472 ENDDO 473 ENDDO 474 ENDDO 475 ! 447 DO jk = 2, jpkm1 448 DO jj = 2, jpjm1 449 DO ji = fs_2, fs_jpim1 ! vector opt. 450 psi(ji,jj,jk) = eps(ji,jj,jk) 451 END DO 452 END DO 453 END DO 454 ! 476 455 CASE( 2 ) ! k-w 477 ! 478 DO jk = 2, jpkm1 479 DO jj = 2, jpjm1 480 DO ji = fs_2, fs_jpim1 ! vector opt. 481 psi(ji,jj,jk) = SQRT(en(ji,jj,jk)) / (rc0 * mxln(ji,jj,jk)) 482 ENDDO 483 ENDDO 484 ENDDO 485 ! 486 CASE( 3 ) ! gen 487 ! 488 DO jk = 2, jpkm1 489 DO jj = 2, jpjm1 490 DO ji = fs_2, fs_jpim1 ! vector opt. 491 psi(ji,jj,jk) = rc02 * en(ji,jj,jk) * mxln(ji,jj,jk)**rnn 492 ENDDO 493 ENDDO 494 ENDDO 495 ! 456 DO jk = 2, jpkm1 457 DO jj = 2, jpjm1 458 DO ji = fs_2, fs_jpim1 ! vector opt. 459 psi(ji,jj,jk) = SQRT( en(ji,jj,jk) ) / ( rc0 * mxln(ji,jj,jk) ) 460 END DO 461 END DO 462 END DO 463 ! 464 CASE( 3 ) ! generic 465 DO jk = 2, jpkm1 466 DO jj = 2, jpjm1 467 DO ji = fs_2, fs_jpim1 ! vector opt. 468 psi(ji,jj,jk) = rc02 * en(ji,jj,jk) * mxln(ji,jj,jk)**rnn 469 END DO 470 END DO 471 END DO 472 ! 496 473 END SELECT 497 474 ! … … 511 488 ! 512 489 ! psi3+ : stable : B=-KhN²<0 => N²>0 if rn2>0 dir = 1 (stable) otherwise dir = 0 (unstable) 513 dir = 0.5 + sign(0.5,rn2(ji,jj,jk))514 ! 515 rpsi3 = dir *rpsi3m+(1-dir)*rpsi3p490 dir = 0.5_wp + SIGN( 0.5_wp, rn2(ji,jj,jk) ) 491 ! 492 rpsi3 = dir * rpsi3m + ( 1._wp - dir ) * rpsi3p 516 493 ! 517 494 ! shear prod. - stratif. destruction … … 519 496 ! 520 497 ! stratif. destruction 521 buoy = rpsi3 * zratio * (- avt(ji,jj,jk) * rn2(ji,jj,jk) )498 buoy = rpsi3 * zratio * (- avt(ji,jj,jk) * rn2(ji,jj,jk) ) 522 499 ! 523 500 ! shear prod. - stratif. destruction 524 501 diss = rpsi2 * zratio * zwall(ji,jj,jk) * eps(ji,jj,jk) 525 502 ! 526 dir = 0.5 + sign(0.5,prod+buoy) ! dir =1(=0) if shear(ji,jj,jk)+buoy >0(<0)527 ! 528 zesh2 = dir *(prod+buoy) +(1-dir)*prod! production term529 zdiss = dir *(diss/psi(ji,jj,jk)) +(1-dir)*(diss-buoy)/psi(ji,jj,jk) ! dissipation term503 dir = 0.5_wp + SIGN( 0.5_wp, prod + buoy ) ! dir =1(=0) if shear(ji,jj,jk)+buoy >0(<0) 504 ! 505 zesh2 = dir * ( prod + buoy ) + (1._wp - dir ) * prod ! production term 506 zdiss = dir * ( diss / psi(ji,jj,jk) ) + (1._wp - dir ) * (diss-buoy) / psi(ji,jj,jk) ! dissipation term 530 507 ! 531 508 ! building the matrix … … 538 515 & / ( fse3t(ji,jj,jk ) * fse3w(ji,jj,jk) ) 539 516 ! diagonal 540 z_elem_b(ji,jj,jk) = 1. - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) &541 & + rdt * zdiss * tmask(ji,jj,jk)517 z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) & 518 & + rdt * zdiss * tmask(ji,jj,jk) 542 519 ! 543 520 ! right hand side in psi … … 547 524 END DO 548 525 ! 549 z_elem_b(:,:,jpk) = 1. 526 z_elem_b(:,:,jpk) = 1._wp 550 527 551 528 ! Surface boundary condition on psi … … 555 532 ! 556 533 CASE ( 0 ) ! Dirichlet boundary conditions 557 !558 IF (ln_crban) THEN! Wave induced mixing case559 ! ! en(1) = q2(1) = 0.5 * (15.8 * Ccb)^(2/3) * u*^2560 ! ! balance between the production and the dissipation terms including the wave effect561 !562 zdep(:,:) = rl_sf * zhsro(:,:)563 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)564 z_elem_a(:,:,1) = psi(:,:,1)565 z_elem_c(:,:,1) = 0.566 z_elem_b(:,:,1) = 1.567 !568 ! one level below569 zdep(:,:) = ( (zhsro(:,:) + fsdepw(:,:,2))**(rmm*ra_sf+rnn) ) &570 & / zhsro(:,:)**(rmm*ra_sf)571 psi (:,:,2) = rsbc_psi1 * ustars2(:,:)**rmm * zdep(:,:) * tmask(:,:,1)572 z_elem_a(:,:,2) = 0.573 z_elem_c(:,:,2) = 0.574 z_elem_b(:,:,2) = 1.575 !576 ELSE ! No wave induced mixing case577 ! ! en(1) = u*^2/C0^2 & l(1) = K*zs578 ! ! balance between the production and the dissipation terms579 !580 zdep(:,:) = vkarmn * zhsro(:,:)581 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)582 z_elem_a(:,:,1) = psi(:,:,1)583 z_elem_c(:,:,1) = 0.584 z_elem_b(:,:,1) = 1.585 !586 ! one level below587 zdep(:,:) = vkarmn * ( zhsro(:,:) + fsdepw(:,:,2) )588 psi (:,:,2) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)589 z_elem_a(:,:,2) = 0.590 z_elem_c(:,:,2) = 0.591 z_elem_b(:,:,2) = 1.592 !593 ENDIF594 !534 ! 535 IF( ln_crban ) THEN ! Wave induced mixing case 536 ! ! en(1) = q2(1) = 0.5 * (15.8 * Ccb)^(2/3) * u*^2 537 ! ! balance between the production and the dissipation terms including the wave effect 538 zdep(:,:) = rl_sf * zhsro(:,:) 539 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 540 z_elem_a(:,:,1) = psi(:,:,1) 541 z_elem_c(:,:,1) = 0._wp 542 z_elem_b(:,:,1) = 1._wp 543 ! 544 ! one level below 545 zex1 = (rmm*ra_sf+rnn) 546 zex2 = (rmm*ra_sf) 547 zdep(:,:) = ( (zhsro(:,:) + fsdepw(:,:,2))**zex1 ) / zhsro(:,:)**zex2 548 psi (:,:,2) = rsbc_psi1 * ustars2(:,:)**rmm * zdep(:,:) * tmask(:,:,1) 549 z_elem_a(:,:,2) = 0._wp 550 z_elem_c(:,:,2) = 0._wp 551 z_elem_b(:,:,2) = 1._wp 552 ! 553 ELSE ! No wave induced mixing case 554 ! ! en(1) = u*^2/C0^2 & l(1) = K*zs 555 ! ! balance between the production and the dissipation terms 556 ! 557 zdep(:,:) = vkarmn * zhsro(:,:) 558 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 559 z_elem_a(:,:,1) = psi(:,:,1) 560 z_elem_c(:,:,1) = 0._wp 561 z_elem_b(:,:,1) = 1._wp 562 ! 563 ! one level below 564 zdep(:,:) = vkarmn * ( zhsro(:,:) + fsdepw(:,:,2) ) 565 psi (:,:,2) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 566 z_elem_a(:,:,2) = 0._wp 567 z_elem_c(:,:,2) = 0._wp 568 z_elem_b(:,:,2) = 1. 569 ! 570 ENDIF 571 ! 595 572 CASE ( 1 ) ! Neumann boundary condition on d(psi)/dz 596 !597 IF (ln_crban) THEN ! Wave induced mixing case598 !599 zdep(:,:) = rl_sf * zhsro(:,:)600 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)601 z_elem_a(:,:,1) = psi(:,:,1)602 z_elem_c(:,:,1) = 0.603 z_elem_b(:,:,1) = 1.604 !605 ! Neumann condition at k=2606 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b607 z_elem_a(:,:,2) = 0.608 !609 ! Set psi vertical flux at the surface:610 zdep(:,:) = (zhsro(:,:) + fsdept(:,:,1))**(rmm*ra_sf+rnn-1.) / zhsro(:,:)**(rmm*ra_sf)611 zflxs(:,:) = rsbc_psi3 * ( zwall_psi(:,:,1)*avm(:,:,1) + zwall_psi(:,:,2)*avm(:,:,2) ) &612 &* en(:,:,1)**rmm * zdep613 psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2)614 !573 ! 574 IF( ln_crban ) THEN ! Wave induced mixing case 575 ! 576 zdep(:,:) = rl_sf * zhsro(:,:) 577 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 578 z_elem_a(:,:,1) = psi(:,:,1) 579 z_elem_c(:,:,1) = 0._wp 580 z_elem_b(:,:,1) = 1._wp 581 ! 582 ! Neumann condition at k=2 583 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 584 z_elem_a(:,:,2) = 0._wp 585 ! 586 ! Set psi vertical flux at the surface: 587 zdep(:,:) = (zhsro(:,:) + fsdept(:,:,1))**(rmm*ra_sf+rnn-1._wp) / zhsro(:,:)**(rmm*ra_sf) 588 zflxs(:,:) = rsbc_psi3 * ( zwall_psi(:,:,1)*avm(:,:,1) + zwall_psi(:,:,2)*avm(:,:,2) ) & 589 & * en(:,:,1)**rmm * zdep 590 psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 591 ! 615 592 ELSE ! No wave induced mixing 616 !617 zdep(:,:) = vkarmn * zhsro(:,:)618 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)619 z_elem_a(:,:,1) = psi(:,:,1)620 z_elem_c(:,:,1) = 0.621 z_elem_b(:,:,1) = 1.622 !623 ! Neumann condition at k=2624 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b625 z_elem_a(ji,jj,2) = 0.626 !627 ! Set psi vertical flux at the surface:628 zdep(:,:)= zhsro(:,:) + fsdept(:,:,1)629 zflxs(:,:) = rsbc_psi2 * ( avm(:,:,1) + avm(:,:,2) ) * en(:,:,1)**rmm * zdep**(rnn-1.)630 psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2)631 !632 ENDIF633 !593 ! 594 zdep(:,:) = vkarmn * zhsro(:,:) 595 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 596 z_elem_a(:,:,1) = psi(:,:,1) 597 z_elem_c(:,:,1) = 0._wp 598 z_elem_b(:,:,1) = 1._wp 599 ! 600 ! Neumann condition at k=2 601 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 602 z_elem_a(ji,jj,2) = 0._wp 603 ! 604 ! Set psi vertical flux at the surface: 605 zdep(:,:) = zhsro(:,:) + fsdept(:,:,1) 606 zflxs(:,:) = rsbc_psi2 * ( avm(:,:,1) + avm(:,:,2) ) * en(:,:,1)**rmm * zdep**(rnn-1._wp) 607 psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 608 ! 609 ENDIF 610 ! 634 611 END SELECT 635 612 … … 639 616 SELECT CASE ( nn_psibc_bot ) 640 617 ! 641 !642 618 CASE ( 0 ) ! Dirichlet 643 ! ! en(ibot) = u*^2 / Co2 and mxln(ibot) = vkarmn * rn_hbro644 ! ! Balance between the production and the dissipation terms619 ! ! en(ibot) = u*^2 / Co2 and mxln(ibot) = vkarmn * rn_hbro 620 ! ! Balance between the production and the dissipation terms 645 621 !CDIR NOVERRCHK 646 DO jj = 2, jpjm1622 DO jj = 2, jpjm1 647 623 !CDIR NOVERRCHK 648 DO ji = fs_2, fs_jpim1 ! vector opt. 649 ibot = mbathy(ji,jj) 650 ibotm1 = ibot-1 651 zdep(ji,jj) = vkarmn * rn_hbro 652 psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn 653 z_elem_a(ji,jj,ibot) = 0. 654 z_elem_c(ji,jj,ibot) = 0. 655 z_elem_b(ji,jj,ibot) = 1. 656 ! 657 ! Just above last level, Dirichlet condition again (GOTM like) 658 zdep(ji,jj) = vkarmn * (rn_hbro + fse3t(ji,jj,ibotm1)) 659 psi (ji,jj,ibotm1) = rc0**rpp * en(ji,jj,ibot )**rmm * zdep(ji,jj)**rnn 660 z_elem_a(ji,jj,ibotm1) = 0. 661 z_elem_c(ji,jj,ibotm1) = 0. 662 z_elem_b(ji,jj,ibotm1) = 1. 663 END DO 664 END DO 665 ! 666 ! 624 DO ji = fs_2, fs_jpim1 ! vector opt. 625 ibot = mbathy(ji,jj) 626 ibotm1 = ibot-1 627 zdep(ji,jj) = vkarmn * rn_hbro 628 psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn 629 z_elem_a(ji,jj,ibot) = 0._wp 630 z_elem_c(ji,jj,ibot) = 0._wp 631 z_elem_b(ji,jj,ibot) = 1._wp 632 ! 633 ! Just above last level, Dirichlet condition again (GOTM like) 634 zdep(ji,jj) = vkarmn * (rn_hbro + fse3t(ji,jj,ibotm1)) 635 psi (ji,jj,ibotm1) = rc0**rpp * en(ji,jj,ibot )**rmm * zdep(ji,jj)**rnn 636 z_elem_a(ji,jj,ibotm1) = 0._wp 637 z_elem_c(ji,jj,ibotm1) = 0._wp 638 z_elem_b(ji,jj,ibotm1) = 1._wp 639 END DO 640 END DO 641 ! 667 642 CASE ( 1 ) ! Neumman boundary condition 668 !643 ! 669 644 !CDIR NOVERRCHK 670 DO jj = 2, jpjm1645 DO jj = 2, jpjm1 671 646 !CDIR NOVERRCHK 672 DO ji = fs_2, fs_jpim1 ! vector opt.673 ibot = mbathy(ji,jj)674 ibotm1 = ibot-1675 !676 ! Bottom level Dirichlet condition:677 zdep(ji,jj) = vkarmn * rn_hbro678 psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn679 !680 z_elem_a(ji,jj,ibot) = 0.681 z_elem_c(ji,jj,ibot) = 0.682 z_elem_b(ji,jj,ibot) = 1.683 !684 ! Just above last level: Neumann condition with flux injection685 z_elem_b(ji,jj,ibotm1) = z_elem_b(ji,jj,ibotm1) + z_elem_c(ji,jj,ibotm1) ! Remove z_elem_c from z_elem_b686 z_elem_c(ji,jj,ibotm1) = 0.687 !688 ! Set psi vertical flux at the bottom:689 zdep(ji,jj) = rn_hbro + 0.5*fse3t(ji,jj,ibotm1)690 zflxb = rsbc_psi2 * ( avm(ji,jj,ibot) + avm(ji,jj,ibotm1) ) *&691 & (0.5*(en(ji,jj,ibot)+en(ji,jj,ibotm1)))**rmm * zdep(ji,jj)**(rnn-1.)692 psi(ji,jj,ibotm1) = psi(ji,jj,ibotm1) + zflxb / fse3w(ji,jj,ibotm1)693 END DO694 END DO695 !647 DO ji = fs_2, fs_jpim1 ! vector opt. 648 ibot = mbathy(ji,jj) 649 ibotm1 = ibot-1 650 ! 651 ! Bottom level Dirichlet condition: 652 zdep(ji,jj) = vkarmn * rn_hbro 653 psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn 654 ! 655 z_elem_a(ji,jj,ibot) = 0._wp 656 z_elem_c(ji,jj,ibot) = 0._wp 657 z_elem_b(ji,jj,ibot) = 1._wp 658 ! 659 ! Just above last level: Neumann condition with flux injection 660 z_elem_b(ji,jj,ibotm1) = z_elem_b(ji,jj,ibotm1) + z_elem_c(ji,jj,ibotm1) ! Remove z_elem_c from z_elem_b 661 z_elem_c(ji,jj,ibotm1) = 0. 662 ! 663 ! Set psi vertical flux at the bottom: 664 zdep(ji,jj) = rn_hbro + 0.5_wp*fse3t(ji,jj,ibotm1) 665 zflxb = rsbc_psi2 * ( avm(ji,jj,ibot) + avm(ji,jj,ibotm1) ) & 666 & * (0.5_wp*(en(ji,jj,ibot)+en(ji,jj,ibotm1)))**rmm * zdep(ji,jj)**(rnn-1._wp) 667 psi(ji,jj,ibotm1) = psi(ji,jj,ibotm1) + zflxb / fse3w(ji,jj,ibotm1) 668 END DO 669 END DO 670 ! 696 671 END SELECT 697 672 … … 727 702 ! 728 703 CASE( 0 ) ! k-kl (Mellor-Yamada) 729 ! 730 DO jk = 1, jpkm1 731 DO jj = 2, jpjm1 732 DO ji = fs_2, fs_jpim1 ! vector opt. 733 eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * en(ji,jj,jk) * SQRT(en(ji,jj,jk)) / psi(ji,jj,jk) 734 ENDDO 735 ENDDO 736 ENDDO 737 ! 704 DO jk = 1, jpkm1 705 DO jj = 2, jpjm1 706 DO ji = fs_2, fs_jpim1 ! vector opt. 707 eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / psi(ji,jj,jk) 708 END DO 709 END DO 710 END DO 711 ! 738 712 CASE( 1 ) ! k-eps 739 ! 740 DO jk = 1, jpkm1 741 DO jj = 2, jpjm1 742 DO ji = fs_2, fs_jpim1 ! vector opt. 743 eps(ji,jj,jk) = psi(ji,jj,jk) 744 ENDDO 745 ENDDO 746 ENDDO 747 ! 713 DO jk = 1, jpkm1 714 DO jj = 2, jpjm1 715 DO ji = fs_2, fs_jpim1 ! vector opt. 716 eps(ji,jj,jk) = psi(ji,jj,jk) 717 END DO 718 END DO 719 END DO 720 ! 748 721 CASE( 2 ) ! k-w 749 ! 750 DO jk = 1, jpkm1 751 DO jj = 2, jpjm1 752 DO ji = fs_2, fs_jpim1 ! vector opt. 753 eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk) 754 ENDDO 755 ENDDO 756 ENDDO 757 ! 758 CASE( 3 ) ! gen 759 ! 760 DO jk = 1, jpkm1 761 DO jj = 2, jpjm1 762 DO ji = fs_2, fs_jpim1 ! vector opt. 763 eps(ji,jj,jk) = rc0**(3.+rpp/rnn) * en(ji,jj,jk)**(1.5+rmm/rnn) * psi(ji,jj,jk)**(-1./rnn) 764 ENDDO 765 ENDDO 766 ENDDO 767 ! 722 DO jk = 1, jpkm1 723 DO jj = 2, jpjm1 724 DO ji = fs_2, fs_jpim1 ! vector opt. 725 eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk) 726 END DO 727 END DO 728 END DO 729 ! 730 CASE( 3 ) ! generic 731 zcoef = rc0**( 3._wp + rpp/rnn ) 732 zex1 = ( 1.5_wp + rmm/rnn ) 733 zex2 = -1._wp / rnn 734 DO jk = 1, jpkm1 735 DO jj = 2, jpjm1 736 DO ji = fs_2, fs_jpim1 ! vector opt. 737 eps(ji,jj,jk) = zcoef * en(ji,jj,jk)**zex1 * psi(ji,jj,jk)**zex2 738 END DO 739 END DO 740 END DO 741 ! 768 742 END SELECT 769 743 … … 775 749 ! limitation 776 750 eps(ji,jj,jk) = MAX( eps(ji,jj,jk), rn_epsmin ) 777 mxln(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk)) / eps(ji,jj,jk)751 mxln(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk) 778 752 ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated) 779 753 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 780 mxln(ji,jj,jk) = MIN( rn_clim_galp*SQRT( 2.*en(ji,jj,jk)/zrn2 ), mxln(ji,jj,jk))754 mxln(ji,jj,jk) = MIN( rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), mxln(ji,jj,jk) ) 781 755 END DO 782 756 END DO … … 790 764 ! 791 765 CASE ( 0 , 1 ) ! Galperin or Kantha-Clayson stability functions 792 ! 793 DO jk = 2, jpkm1 794 DO jj = 2, jpjm1 795 DO ji = fs_2, fs_jpim1 ! vector opt. 796 ! zcof = l²/q² 797 zcof = mxlb(ji,jj,jk)**2. / ( 2.*eb(ji,jj,jk) ) 798 ! Gh = -N²l²/q² 799 gh = - rn2(ji,jj,jk) * zcof 800 gh = MIN( gh, rgh0 ) 801 gh = MAX( gh, rghmin ) 802 ! Stability functions from Kantha and Clayson (if C2=C3=0 => Galperin) 803 sh = ra2*( 1.-6.*ra1/rb1 ) / ( 1.-3.*ra2*gh*(6.*ra1+rb2*( 1.-rc3 ) ) ) 804 sm = ( rb1**(-1./3.) + ( 18*ra1*ra1 + 9.*ra1*ra2*(1.-rc2) )*sh*gh ) / (1.-9.*ra1*ra2*gh) 805 ! 806 ! Store stability function in avmu and avmv 807 avmu(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk) 808 avmv(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk) 809 END DO 810 END DO 811 END DO 812 ! 766 DO jk = 2, jpkm1 767 DO jj = 2, jpjm1 768 DO ji = fs_2, fs_jpim1 ! vector opt. 769 ! zcof = l²/q² 770 zcof = mxlb(ji,jj,jk) * mxlb(ji,jj,jk) / ( 2._wp*eb(ji,jj,jk) ) 771 ! Gh = -N²l²/q² 772 gh = - rn2(ji,jj,jk) * zcof 773 gh = MIN( gh, rgh0 ) 774 gh = MAX( gh, rghmin ) 775 ! Stability functions from Kantha and Clayson (if C2=C3=0 => Galperin) 776 sh = ra2*( 1._wp-6._wp*ra1/rb1 ) / ( 1.-3.*ra2*gh*(6.*ra1+rb2*( 1._wp-rc3 ) ) ) 777 sm = ( rb1**(-1._wp/3._wp) + ( 18._wp*ra1*ra1 + 9._wp*ra1*ra2*(1._wp-rc2) )*sh*gh ) / (1._wp-9._wp*ra1*ra2*gh) 778 ! 779 ! Store stability function in avmu and avmv 780 avmu(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk) 781 avmv(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk) 782 END DO 783 END DO 784 END DO 785 ! 813 786 CASE ( 2, 3 ) ! Canuto stability functions 814 ! 815 DO jk = 2, jpkm1 816 DO jj = 2, jpjm1 817 DO ji = fs_2, fs_jpim1 ! vector opt. 818 ! zcof = l²/q² 819 zcof = mxlb(ji,jj,jk)**2. / ( 2.*eb(ji,jj,jk) ) 820 ! Gh = -N²l²/q² 821 gh = - rn2(ji,jj,jk) * zcof 822 gh = MIN( gh, rgh0 ) 823 gh = MAX( gh, rghmin ) 824 gh = gh*rf6 825 ! Gm = M²l²/q² Shear number 826 shr = shear(ji,jj,jk)/MAX(avm(ji,jj,jk), rsmall) 827 gm = shr * zcof 828 gm = MAX(gm, 1.e-10) 829 gm = gm*rf6 830 gm = MIN ( (rd0 - rd1*gh + rd3*gh**2.)/(rd2-rd4*gh) , gm ) 831 ! Stability functions from Canuto 832 rcff = rd0 - rd1*gh +rd2*gm + rd3*gh**2. - rd4*gh*gm + rd5*gm**2. 833 sm = (rs0 - rs1*gh + rs2*gm) / rcff 834 sh = (rs4 - rs5*gh + rs6*gm) / rcff 835 ! 836 ! Store stability function in avmu and avmv 837 avmu(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk) 838 avmv(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk) 839 END DO 840 END DO 841 END DO 842 ! 787 DO jk = 2, jpkm1 788 DO jj = 2, jpjm1 789 DO ji = fs_2, fs_jpim1 ! vector opt. 790 ! zcof = l²/q² 791 zcof = mxlb(ji,jj,jk)*mxlb(ji,jj,jk) / ( 2._wp * eb(ji,jj,jk) ) 792 ! Gh = -N²l²/q² 793 gh = - rn2(ji,jj,jk) * zcof 794 gh = MIN( gh, rgh0 ) 795 gh = MAX( gh, rghmin ) 796 gh = gh * rf6 797 ! Gm = M²l²/q² Shear number 798 shr = shear(ji,jj,jk) / MAX( avm(ji,jj,jk), rsmall ) 799 gm = MAX( shr * zcof , 1.e-10 ) 800 gm = gm * rf6 801 gm = MIN ( (rd0 - rd1*gh + rd3*gh*gh) / (rd2-rd4*gh) , gm ) 802 ! Stability functions from Canuto 803 rcff = rd0 - rd1*gh +rd2*gm + rd3*gh*gh - rd4*gh*gm + rd5*gm*gm 804 sm = (rs0 - rs1*gh + rs2*gm) / rcff 805 sh = (rs4 - rs5*gh + rs6*gm) / rcff 806 ! 807 ! Store stability function in avmu and avmv 808 avmu(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk) 809 avmv(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk) 810 END DO 811 END DO 812 END DO 813 ! 843 814 END SELECT 844 815 845 816 ! Boundary conditions on stability functions for momentum (Neumann): 846 817 ! Lines below are useless if GOTM style Dirichlet conditions are used 818 zcoef = rcm_sf / SQRT( 2._wp ) 847 819 DO jj = 2, jpjm1 848 820 DO ji = fs_2, fs_jpim1 ! vector opt. 849 avmv(ji,jj,1) = rcm_sf / SQRT(2.) 850 END DO 851 END DO 821 avmv(ji,jj,1) = zcoef 822 END DO 823 END DO 824 zcoef = rc0 / SQRT( 2._wp ) 852 825 DO jj = 2, jpjm1 853 826 DO ji = fs_2, fs_jpim1 ! vector opt. 854 ibot=mbathy(ji,jj) 855 ibotm1=ibot-1 856 avmv(ji,jj,ibot) = rc0 / SQRT(2.) 827 ibot = mbathy(ji,jj) 828 avmv(ji,jj,ibot) = zcoef 857 829 END DO 858 830 END DO … … 863 835 DO jj = 2, jpjm1 864 836 DO ji = fs_2, fs_jpim1 ! vector opt. 865 zsqen = SQRT(2.*en(ji,jj,jk)) * mxln(ji,jj,jk) 866 zav = zsqen * avmu(ji,jj,jk) 867 avt(ji,jj,jk) = MAX(zav,avtb(jk))*tmask(ji,jj,jk) ! apply mask for zdfmxl routine 868 zav = zsqen * avmv(ji,jj,jk) 869 avm(ji,jj,jk) = MAX(zav,avmb(jk)) ! Note that avm is not masked at the surface and the bottom 870 END DO 871 END DO 872 END DO 873 837 zsqen = SQRT( 2._wp * en(ji,jj,jk) ) * mxln(ji,jj,jk) 838 zav = zsqen * avmu(ji,jj,jk) 839 avt(ji,jj,jk) = MAX( zav, avtb(jk) )*tmask(ji,jj,jk) ! apply mask for zdfmxl routine 840 zav = zsqen * avmv(ji,jj,jk) 841 avm(ji,jj,jk) = MAX( zav, avmb(jk) ) ! Note that avm is not masked at the surface and the bottom 842 END DO 843 END DO 844 END DO 874 845 ! 875 846 ! Lateral boundary conditions (sign unchanged) 876 ! 877 avt(:,:,1) = 0. 847 avt(:,:,1) = 0._wp 878 848 CALL lbc_lnk( avm, 'W', 1. ) ; CALL lbc_lnk( avt, 'W', 1. ) 879 849 … … 881 851 DO jj = 2, jpjm1 882 852 DO ji = fs_2, fs_jpim1 ! vector opt. 883 avmu(ji,jj,jk) = ( avm (ji,jj,jk)*tmask(ji,jj,jk) + avm (ji+1,jj ,jk)*tmask(ji+1,jj ,jk) ) * umask(ji,jj,jk) & 884 & / MAX( 1., tmask(ji,jj,jk) + tmask(ji+1,jj ,jk) ) 885 avmv(ji,jj,jk) = ( avm (ji,jj,jk)*tmask(ji,jj,jk) + avm (ji ,jj+1,jk)*tmask(ji ,jj+1,jk) ) * vmask(ji,jj,jk) & 886 & / MAX( 1., tmask(ji,jj,jk) + tmask(ji ,jj+1,jk) ) 887 END DO 888 END DO 889 END DO 890 891 avmu(:,:,1) = 0. 892 avmv(:,:,1) = 0. 893 894 CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. ) ! Lateral boundary conditions 853 avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj ,jk) ) * umask(ji,jj,jk) 854 avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji ,jj+1,jk) ) * vmask(ji,jj,jk) 855 END DO 856 END DO 857 END DO 858 avmu(:,:,1) = 0._wp ; avmv(:,:,1) = 0._wp ! set surface to zero 859 CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. ) ! Lateral boundary conditions 895 860 896 861 IF(ln_ctl) THEN … … 918 883 !! 919 884 !!---------------------------------------------------------------------- 885 USE dynzdf_exp 886 USE trazdf_exp 887 ! 920 888 INTEGER :: jk ! dummy loop indices 921 889 REAL(wp):: zcr ! local scalar … … 929 897 !!---------------------------------------------------------- 930 898 931 ! Read Namelist namzdf_gls 932 ! ------------------------ 933 REWIND ( numnam ) 899 REWIND ( numnam ) !* Read Namelist namzdf_gls 934 900 READ ( numnam, namzdf_gls ) 935 901 936 ! Parameter control and print 937 ! --------------------------- 938 ! Control print 939 IF(lwp) THEN 902 IF(lwp) THEN !* Control print 940 903 WRITE(numout,*) 941 904 WRITE(numout,*) 'zdf_gls_init : gls turbulent closure scheme' 942 905 WRITE(numout,*) '~~~~~~~~~~~~' 943 WRITE(numout,*) ' 944 WRITE(numout,*) ' minimum value of en rn_emin= ', rn_emin945 WRITE(numout,*) ' minimum value of eps rn_epsmin= ', rn_epsmin946 WRITE(numout,*) ' Surface roughness (m) hsro= ', hsro947 WRITE(numout,*) ' Bottom roughness (m) rn_hbro= ', rn_hbro948 WRITE(numout,*) ' Limit dissipation rate under stable stratification ln_length_lim = ',ln_length_lim949 WRITE(numout,*) ' Galperin length scale limitation coef (Standard: 0.53, Holt: 0.26) rn_clim_galp= ', rn_clim_galp950 WRITE(numout,*) ' TKE Surface boundary conditionnn_tkebc_surf = ', nn_tkebc_surf951 WRITE(numout,*) ' TKE Bottom boundary conditionnn_tkebc_bot = ', nn_tkebc_bot952 WRITE(numout,*) ' PSI Surface boundary conditionnn_psibc_surf = ', nn_psibc_surf953 WRITE(numout,*) ' PSI Bottom boundary conditionnn_psibc_bot = ', nn_psibc_bot954 WRITE(numout,*) ' Craig and Banner scheme ln_crban= ', ln_crban955 WRITE(numout,*) ' Modify psi Schmidt number (wb case) ln_sigpsi= ', ln_sigpsi956 WRITE(numout,*) ' Craig and Banner coef. rn_crban= ', rn_crban957 WRITE(numout,*) ' Charnock coef. rn_charn= ', rn_charn958 WRITE(numout,*) ' Stability functions nn_stab_func = ',nn_stab_func959 WRITE(numout,*) ' Closure nn_clos = ',nn_clos906 WRITE(numout,*) ' Namelist namzdf_gls : set gls mixing parameters' 907 WRITE(numout,*) ' minimum value of en rn_emin = ', rn_emin 908 WRITE(numout,*) ' minimum value of eps rn_epsmin = ', rn_epsmin 909 WRITE(numout,*) ' Surface roughness (m) hsro = ', hsro 910 WRITE(numout,*) ' Bottom roughness (m) rn_hbro = ', rn_hbro 911 WRITE(numout,*) ' Limit dissipation rate under stable stratif. ln_length_lim = ', ln_length_lim 912 WRITE(numout,*) ' Galperin limit (Standard: 0.53, Holt: 0.26) rn_clim_galp = ', rn_clim_galp 913 WRITE(numout,*) ' TKE Surface boundary condition nn_tkebc_surf = ', nn_tkebc_surf 914 WRITE(numout,*) ' TKE Bottom boundary condition nn_tkebc_bot = ', nn_tkebc_bot 915 WRITE(numout,*) ' PSI Surface boundary condition nn_psibc_surf = ', nn_psibc_surf 916 WRITE(numout,*) ' PSI Bottom boundary condition nn_psibc_bot = ', nn_psibc_bot 917 WRITE(numout,*) ' Craig and Banner scheme ln_crban = ', ln_crban 918 WRITE(numout,*) ' Modify psi Schmidt number (wb case) ln_sigpsi = ', ln_sigpsi 919 WRITE(numout,*) ' Craig and Banner coefficient rn_crban = ', rn_crban 920 WRITE(numout,*) ' Charnock coefficient rn_charn = ', rn_charn 921 WRITE(numout,*) ' Stability functions nn_stab_func = ', nn_stab_func 922 WRITE(numout,*) ' Type of closure nn_clos = ', nn_clos 960 923 WRITE(numout,*) 961 924 ENDIF 962 925 963 ! Check of some namelist values926 ! !* Check of some namelist values 964 927 IF( nn_tkebc_surf < 0 .OR. nn_tkebc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_tkebc_surf is 0 or 1' ) 965 928 IF( nn_psibc_surf < 0 .OR. nn_psibc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_psibc_surf is 0 or 1' ) 966 IF( nn_tkebc_bot < 0 .OR. nn_tkebc_bot> 1 ) CALL ctl_stop( 'bad flag: nn_tkebc_bot is 0 or 1' )967 IF( nn_psibc_bot < 0 .OR. nn_psibc_bot> 1 ) CALL ctl_stop( 'bad flag: nn_psibc_bot is 0 or 1' )968 IF( nn_stab_func < 0 .OR. nn_stab_func > 3) CALL ctl_stop( 'bad flag: nn_stab_func is 0, 1, 2 and 3' )969 IF( nn_clos < 0 .OR. nn_clos > 3) CALL ctl_stop( 'bad flag: nn_clos is 0, 1, 2 or 3' )929 IF( nn_tkebc_bot < 0 .OR. nn_tkebc_bot > 1 ) CALL ctl_stop( 'bad flag: nn_tkebc_bot is 0 or 1' ) 930 IF( nn_psibc_bot < 0 .OR. nn_psibc_bot > 1 ) CALL ctl_stop( 'bad flag: nn_psibc_bot is 0 or 1' ) 931 IF( nn_stab_func < 0 .OR. nn_stab_func > 3 ) CALL ctl_stop( 'bad flag: nn_stab_func is 0, 1, 2 and 3' ) 932 IF( nn_clos < 0 .OR. nn_clos > 3 ) CALL ctl_stop( 'bad flag: nn_clos is 0, 1, 2 or 3' ) 970 933 971 934 ! Initialisation of the parameters for the choosen closure … … 975 938 ! 976 939 CASE( 0 ) ! k-kl (Mellor-Yamada) 977 !978 IF(lwp) WRITE(numout,*) 'The choosen closure is k-kl closed to the classical Mellor-Yamada'979 rpp = 0.980 rmm = 1.981 rnn = 1.982 rsc_tke = 1.96983 rsc_psi = 1.96984 rpsi1 = 0.9985 rpsi3p = 1.986 rpsi2 = 0.5987 !940 ! 941 IF(lwp) WRITE(numout,*) 'The choosen closure is k-kl closed to the classical Mellor-Yamada' 942 rpp = 0._wp 943 rmm = 1._wp 944 rnn = 1._wp 945 rsc_tke = 1.96_wp 946 rsc_psi = 1.96_wp 947 rpsi1 = 0.9_wp 948 rpsi3p = 1._wp 949 rpsi2 = 0.5_wp 950 ! 988 951 SELECT CASE ( nn_stab_func ) 989 ! 990 CASE( 0, 1 ) ! G88 or KC stability functions 991 rpsi3m = 2.53 992 CASE( 2 ) ! Canuto A stability functions 993 rpsi3m = 2.38 994 CASE( 3 ) ! Canuto B stability functions 995 rpsi3m = 2.38 ! caution : constant not identified 952 CASE( 0, 1 ) ; rpsi3m = 2.53_wp ! G88 or KC stability functions 953 CASE( 2 ) ; rpsi3m = 2.38_wp ! Canuto A stability functions 954 CASE( 3 ) ; rpsi3m = 2.38 ! Canuto B stability functions (caution : constant not identified) 996 955 END SELECT 997 !956 ! 998 957 CASE( 1 ) ! k-eps 999 ! 1000 IF(lwp) WRITE(numout,*) 'The choosen closure is k-eps' 1001 rpp = 3. 1002 rmm = 1.5 1003 rnn = -1. 1004 rsc_tke = 1. 1005 rsc_psi = 1.3 ! Schmidt number for psi 1006 rpsi1 = 1.44 1007 rpsi3p = 1. 1008 rpsi2 = 1.92 1009 ! 1010 SELECT CASE ( nn_stab_func ) 1011 ! 1012 CASE( 0, 1 ) ! G88 or KC stability functions 1013 rpsi3m = -0.52 1014 CASE( 2 ) ! Canuto A stability functions 1015 rpsi3m = -0.629 1016 CASE( 3 ) ! Canuto B stability functions 1017 rpsi3m = -0.566 1018 END SELECT 1019 ! 958 ! 959 IF(lwp) WRITE(numout,*) 'The choosen closure is k-eps' 960 rpp = 3._wp 961 rmm = 1.5_wp 962 rnn = -1._wp 963 rsc_tke = 1._wp 964 rsc_psi = 1.3_wp ! Schmidt number for psi 965 rpsi1 = 1.44_wp 966 rpsi3p = 1._wp 967 rpsi2 = 1.92_wp 968 ! 969 SELECT CASE ( nn_stab_func ) 970 CASE( 0, 1 ) ; rpsi3m = -0.52_wp ! G88 or KC stability functions 971 CASE( 2 ) ; rpsi3m = -0.629_wp ! Canuto A stability functions 972 CASE( 3 ) ; rpsi3m = -0.566 ! Canuto B stability functions 973 END SELECT 974 ! 1020 975 CASE( 2 ) ! k-omega 1021 ! 1022 IF(lwp) WRITE(numout,*) 'The choosen closure is k-omega' 1023 rpp = -1. 1024 rmm = 0.5 1025 rnn = -1. 1026 rsc_tke = 2. 1027 rsc_psi = 2. 1028 rpsi1 = 0.555 1029 rpsi3p = 1. 1030 rpsi2 = 0.833 1031 ! 1032 SELECT CASE ( nn_stab_func ) 1033 ! 1034 CASE( 0, 1 ) ! G88 or KC stability functions 1035 rpsi3m = -0.58 1036 CASE( 2 ) ! Canuto A stability functions 1037 rpsi3m = -0.64 1038 CASE( 3 ) ! Canuto B stability functions 1039 rpsi3m = -0.64 ! caution : constant not identified 1040 END SELECT 1041 ! 1042 CASE( 3 ) ! gen 1043 ! 1044 IF(lwp) WRITE(numout,*) 'The choosen closure is generic' 1045 rpp = 2. 1046 rmm = 1. 1047 rnn = -0.67 1048 rsc_tke = 0.8 1049 rsc_psi = 1.07 1050 rpsi1 = 1. 1051 rpsi3p = 1. 1052 rpsi2 = 1.22 1053 ! 1054 SELECT CASE ( nn_stab_func ) 1055 ! 1056 CASE( 0, 1 ) ! G88 or KC stability functions 1057 rpsi3m = 0.1 1058 CASE( 2 ) ! Canuto A stability functions 1059 rpsi3m = 0.05 1060 CASE( 3 ) ! Canuto B stability functions 1061 rpsi3m = 0.05 ! caution : constant not identified 1062 END SELECT 1063 ! 976 ! 977 IF(lwp) WRITE(numout,*) 'The choosen closure is k-omega' 978 rpp = -1._wp 979 rmm = 0.5_wp 980 rnn = -1._wp 981 rsc_tke = 2._wp 982 rsc_psi = 2._wp 983 rpsi1 = 0.555_wp 984 rpsi3p = 1._wp 985 rpsi2 = 0.833_wp 986 ! 987 SELECT CASE ( nn_stab_func ) 988 CASE( 0, 1 ) ; rpsi3m = -0.58_wp ! G88 or KC stability functions 989 CASE( 2 ) ; rpsi3m = -0.64_wp ! Canuto A stability functions 990 CASE( 3 ) ; rpsi3m = -0.64_wp ! Canuto B stability functions caution : constant not identified) 991 END SELECT 992 ! 993 CASE( 3 ) ! generic 994 ! 995 IF(lwp) WRITE(numout,*) 'The choosen closure is generic' 996 rpp = 2._wp 997 rmm = 1._wp 998 rnn = -0.67_wp 999 rsc_tke = 0.8_wp 1000 rsc_psi = 1.07_wp 1001 rpsi1 = 1._wp 1002 rpsi3p = 1._wp 1003 rpsi2 = 1.22_wp 1004 ! 1005 SELECT CASE ( nn_stab_func ) 1006 CASE( 0, 1 ) ; rpsi3m = 0.1_wp ! G88 or KC stability functions 1007 CASE( 2 ) ; rpsi3m = 0.05_wp ! Canuto A stability functions 1008 CASE( 3 ) ; rpsi3m = 0.05_wp ! Canuto B stability functions caution : constant not identified) 1009 END SELECT 1010 ! 1064 1011 END SELECT 1065 1012 … … 1070 1017 ! 1071 1018 CASE ( 0 ) ! Galperin stability functions 1072 !1073 IF(lwp) WRITE(numout,*) 'Stability functions from Galperin'1074 rc2 = 0.1075 rc3 = 0.1076 rc_diff = 1.1077 rc0 = 0.55441078 rcm_sf = 0.98841079 rghmin = -0.281080 rgh0 = 0.02331081 rghcri = 0.021082 !1019 ! 1020 IF(lwp) WRITE(numout,*) 'Stability functions from Galperin' 1021 rc2 = 0._wp 1022 rc3 = 0._wp 1023 rc_diff = 1._wp 1024 rc0 = 0.5544_wp 1025 rcm_sf = 0.9884_wp 1026 rghmin = -0.28_wp 1027 rgh0 = 0.0233_wp 1028 rghcri = 0.02_wp 1029 ! 1083 1030 CASE ( 1 ) ! Kantha-Clayson stability functions 1084 !1085 IF(lwp) WRITE(numout,*) 'Stability functions from Kantha-Clayson'1086 rc2 = 0.71087 rc3 = 0.21088 rc_diff = 1.1089 rc0 = 0.55441090 rcm_sf = 0.98841091 rghmin = -0.281092 rgh0 = 0.02331093 rghcri = 0.021094 !1031 ! 1032 IF(lwp) WRITE(numout,*) 'Stability functions from Kantha-Clayson' 1033 rc2 = 0.7_wp 1034 rc3 = 0.2_wp 1035 rc_diff = 1._wp 1036 rc0 = 0.5544_wp 1037 rcm_sf = 0.9884_wp 1038 rghmin = -0.28_wp 1039 rgh0 = 0.0233_wp 1040 rghcri = 0.02_wp 1041 ! 1095 1042 CASE ( 2 ) ! Canuto A stability functions 1096 ! 1097 IF(lwp) WRITE(numout,*) 'Stability functions from Canuto A' 1098 rs0 = 1.5*rl1*rl5**2. 1099 rs1 = -rl4*(rl6+rl7) + 2.*rl4*rl5*(rl1-(1./3.)*rl2-rl3)+1.5*rl1*rl5*rl8 1100 rs2 = -(3./8.)*rl1*(rl6**2.-rl7**2.) 1101 rs4 = 2.*rl5 1102 rs5 = 2.*rl4 1103 rs6 = (2./3.)*rl5*(3.*rl3**2.-rl2**2.)-0.5*rl5*rl1*(3.*rl3-rl2)+0.75*rl1*(rl6-rl7) 1104 rd0 = 3*rl5**2. 1105 rd1 = rl5*(7.*rl4+3.*rl8) 1106 rd2 = rl5**2.*(3.*rl3**2.-rl2**2.)-0.75*(rl6**2.-rl7**2.) 1107 rd3 = rl4*(4.*rl4+3.*rl8) 1108 rd4 = rl4*(rl2*rl6-3.*rl3*rl7-rl5*(rl2**2.-rl3**2.))+rl5*rl8*(3.*rl3**2.-rl2**2.) 1109 rd5 = 0.25*(rl2**2.-3.*rl3**2.)*(rl6**2.-rl7**2.) 1110 rc0 = 0.5268 1111 rf6 = 8. / (rc0**6.) 1112 rc_diff = SQRT(2.)/(rc0**3.) 1113 rcm_sf = 0.7310 1114 rghmin = -0.28 1115 rgh0 = 0.0329 1116 rghcri = 0.03 1117 ! 1043 ! 1044 IF(lwp) WRITE(numout,*) 'Stability functions from Canuto A' 1045 rs0 = 1.5_wp * rl1 * rl5*rl5 1046 rs1 = -rl4*(rl6+rl7) + 2._wp*rl4*rl5*(rl1-(1._wp/3._wp)*rl2-rl3) + 1.5_wp*rl1*rl5*rl8 1047 rs2 = -(3._wp/8._wp) * rl1*(rl6*rl6-rl7*rl7) 1048 rs4 = 2._wp * rl5 1049 rs5 = 2._wp * rl4 1050 rs6 = (2._wp/3._wp) * rl5 * ( 3._wp*rl3*rl3 - rl2*rl2 ) - 0.5_wp * rl5*rl1 * (3._wp*rl3-rl2) & 1051 & + 0.75_wp * rl1 * ( rl6 - rl7 ) 1052 rd0 = 3._wp * rl5*rl5 1053 rd1 = rl5 * ( 7._wp*rl4 + 3._wp*rl8 ) 1054 rd2 = rl5*rl5 * ( 3._wp*rl3*rl3 - rl2*rl2 ) - 0.75_wp*(rl6*rl6 - rl7*rl7 ) 1055 rd3 = rl4 * ( 4._wp*rl4 + 3._wp*rl8) 1056 rd4 = rl4 * ( rl2 * rl6 - 3._wp*rl3*rl7 - rl5*(rl2*rl2 - rl3*rl3 ) ) + rl5*rl8 * ( 3._wp*rl3*rl3 - rl2*rl2 ) 1057 rd5 = 0.25_wp * ( rl2*rl2 - 3._wp *rl3*rl3 ) * ( rl6*rl6 - rl7*rl7 ) 1058 rc0 = 0.5268_wp 1059 rf6 = 8._wp / (rc0**6._wp) 1060 rc_diff = SQRT(2._wp) / (rc0**3._wp) 1061 rcm_sf = 0.7310_wp 1062 rghmin = -0.28_wp 1063 rgh0 = 0.0329_wp 1064 rghcri = 0.03_wp 1065 ! 1118 1066 CASE ( 3 ) ! Canuto B stability functions 1119 !1120 IF(lwp) WRITE(numout,*) 'Stability functions from Canuto B'1121 rs0 = 1.5*rm1*rm5**2.1122 rs1 = -rm4*(rm6+rm7) + 2.*rm4*rm5*(rm1-(1./3.)*rm2-rm3)+1.5*rm1*rm5*rm81123 rs2 = -(3./8.)*rm1*(rm6**2.-rm7**2.)1124 rs4 = 2.*rm51125 rs5 = 2.*rm41126 rs6 = (2./3.)*rm5*(3.*rm3**2.-rm2**2.)-0.5*rm5*rm1*(3.*rm3-rm2)+0.75*rm1*(rm6-rm7)1127 rd0 = 3*rm5**2.1128 rd1 = rm5*(7.*rm4+3.*rm8)1129 rd2 = rm5**2.*(3.*rm3**2.-rm2**2.)-0.75*(rm6**2.-rm7**2.)1130 rd3 = rm4*(4.*rm4+3.*rm8)1131 rd4 = rm4*(rm2*rm6-3.*rm3*rm7-rm5*(rm2**2.-rm3**2.))+rm5*rm8*(3.*rm3**2.-rm2**2.)1132 rd5 = 0.25*(rm2**2.-3.*rm3**2.)*(rm6**2.-rm7**2.)1133 rc0 = 0.5268 !! rc0 = 0.5540(Warner ...) to verify !1134 rf6 = 8. / (rc0**6.)1135 rc_diff = SQRT(2.)/(rc0**3.)1136 rcm_sf = 0.74701137 rghmin = -0.281138 rgh0 = 0.04441139 rghcri = 0.04141140 !1067 ! 1068 IF(lwp) WRITE(numout,*) 'Stability functions from Canuto B' 1069 rs0 = 1.5_wp * rm1 * rm5*rm5 1070 rs1 = -rm4 * (rm6+rm7) + 2._wp * rm4*rm5*(rm1-(1._wp/3._wp)*rm2-rm3) + 1.5_wp * rm1*rm5*rm8 1071 rs2 = -(3._wp/8._wp) * rm1 * (rm6*rm6-rm7*rm7 ) 1072 rs4 = 2._wp * rm5 1073 rs5 = 2._wp * rm4 1074 rs6 = (2._wp/3._wp) * rm5 * (3._wp*rm3*rm3-rm2*rm2) - 0.5_wp * rm5*rm1*(3._wp*rm3-rm2) + 0.75_wp * rm1*(rm6-rm7) 1075 rd0 = 3._wp * rm5*rm5 1076 rd1 = rm5 * (7._wp*rm4 + 3._wp*rm8) 1077 rd2 = rm5*rm5 * (3._wp*rm3*rm3 - rm2*rm2) - 0.75_wp * (rm6*rm6 - rm7*rm7) 1078 rd3 = rm4 * ( 4._wp*rm4 + 3._wp*rm8 ) 1079 rd4 = rm4 * ( rm2*rm6 -3._wp*rm3*rm7 - rm5*(rm2*rm2 - rm3*rm3) ) + rm5 * rm8 * ( 3._wp*rm3*rm3 - rm2*rm2 ) 1080 rd5 = 0.25_wp * ( rm2*rm2 - 3._wp*rm3*rm3 ) * ( rm6*rm6 - rm7*rm7 ) 1081 rc0 = 0.5268_wp !! rc0 = 0.5540_wp (Warner ...) to verify ! 1082 rf6 = 8._wp / ( rc0**6._wp ) 1083 rc_diff = SQRT(2._wp)/(rc0**3.) 1084 rcm_sf = 0.7470_wp 1085 rghmin = -0.28_wp 1086 rgh0 = 0.0444_wp 1087 rghcri = 0.0414_wp 1088 ! 1141 1089 END SELECT 1142 1090 1143 ! Set Schmidt number for psi diffusion 1144 ! In the wave breaking case 1091 ! Set Schmidt number for psi diffusion in the wave breaking case 1145 1092 ! See equation 13 of Carniel et al, Ocean modelling, 30, 225-239, 2009 1146 1093 ! or equation (17) of Burchard, JPO, 31, 3133-3145, 2001 1147 IF ((ln_sigpsi).AND.(ln_crban)) THEN1148 zcr = SQRT( 1.5*rsc_tke) * rcm_sf /vkarmn1149 rsc_psi0 = vkarmn* *2/(rpsi2*rcm_sf**2) *&1150 & ( rnn**2 - 4./3.*zcr*rnn*rmm - 1./3.*zcr*rnn&1151 & + 2./9.*rmm*zcr**2 + 4./9.*zcr**2*rmm**2)1094 IF( ln_sigpsi .AND. ln_crban ) THEN 1095 zcr = SQRT( 1.5_wp*rsc_tke ) * rcm_sf / vkarmn 1096 rsc_psi0 = vkarmn*vkarmn / ( rpsi2 * rcm_sf*rcm_sf ) & 1097 & * ( rnn*rnn - 4._wp/3._wp * zcr*rnn*rmm - 1._wp/3._wp * zcr*rnn & 1098 & + 2._wp/9._wp * rmm * zcr*zcr + 4._wp/9._wp * zcr*zcr * rmm*rmm ) 1152 1099 ELSE 1153 1100 rsc_psi0 = rsc_psi … … 1156 1103 ! Shear free turbulence parameters: 1157 1104 ! 1158 ra_sf = -4.*rnn*SQRT(rsc_tke) / ( (1.+4.*rmm)*SQRT(rsc_tke) & 1159 & - SQRT(rsc_tke + 24.*rsc_psi0*rpsi2 ) ) 1160 rl_sf = rc0 * SQRT(rc0/rcm_sf) * SQRT( ( (1. + 4.*rmm + 8.*rmm**2)*rsc_tke & 1161 & + 12. * rsc_psi0*rpsi2 - (1. + 4.*rmm) & 1162 & *SQRT(rsc_tke*(rsc_tke & 1163 & + 24.*rsc_psi0*rpsi2)) ) & 1164 & /(12.*rnn**2.) & 1165 & ) 1166 1167 ! Control print 1168 ! 1169 IF(lwp) THEN 1105 ra_sf = -4._wp * rnn * SQRT( rsc_tke ) / ( (1._wp+4._wp*rmm) * SQRT( rsc_tke ) & 1106 & - SQRT(rsc_tke + 24._wp*rsc_psi0*rpsi2 ) ) 1107 rl_sf = rc0 * SQRT( rc0 / rcm_sf ) & 1108 & * SQRT( ( (1._wp + 4._wp*rmm + 8._wp*rmm*rmm) * rsc_tke & 1109 & + 12._wp * rsc_psi0 * rpsi2 & 1110 & - (1._wp + 4._wp*rmm) * SQRT( rsc_tke*(rsc_tke+ 24._wp*rsc_psi0*rpsi2) ) ) & 1111 & / ( 12._wp*rnn*rnn ) ) 1112 1113 ! 1114 IF(lwp) THEN ! Control print 1170 1115 WRITE(numout,*) 1171 1116 WRITE(numout,*) 'Limit values' … … 1191 1136 1192 1137 ! Constants initialization 1193 rc02r = 1. / rc0**2. 1194 rc02 = rc0**2._wp 1195 rc03 = rc0**3._wp 1196 rc04 = rc0**4._wp 1138 rc02 = rc0 * rc0 ; rc02r = 1. / rc02 1139 rc03 = rc02 * rc0 1140 rc04 = rc03 * rc0 1197 1141 rc03_sqrt2_galp = rc03 / SQRT(2._wp) / rn_clim_galp 1198 rsbc_mb = 0.5 * (15.8*rn_crban)**(2./3.)! Surf. bound. cond. from Mellor and Blumberg1199 rsbc_std = 3.75 ! Surf. bound. cond. standard (prod=diss)1200 rsbc_tke1 = (-rsc_tke*rn_crban/(rcm_sf*ra_sf*rl_sf))**(2. /3.)! k_eps = 53. Dirichlet + Wave breaking1201 rsbc_tke2 = 0.5 / rau01142 rsbc_mb = 0.5_wp * (15.8_wp*rn_crban)**(2._wp/3._wp) ! Surf. bound. cond. from Mellor and Blumberg 1143 rsbc_std = 3.75_wp ! Surf. bound. cond. standard (prod=diss) 1144 rsbc_tke1 = (-rsc_tke*rn_crban/(rcm_sf*ra_sf*rl_sf))**(2._wp/3._wp) ! k_eps = 53. Dirichlet + Wave breaking 1145 rsbc_tke2 = 0.5_wp / rau0 1202 1146 rsbc_tke3 = rdt * rn_crban ! Neumann + Wave breaking 1203 rsbc_zs = rn_charn /grav! Charnock formula1147 rsbc_zs = rn_charn / grav ! Charnock formula 1204 1148 rsbc_psi1 = rc0**rpp * rsbc_tke1**rmm * rl_sf**rnn ! Dirichlet + Wave breaking 1205 rsbc_psi2 = -0.5 * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking1206 rsbc_psi3 = -0.5 * rdt * rc0**rpp * rl_sf**rnn / rsc_psi * (rnn + rmm*ra_sf) ! Neumann + Wave breaking1207 rfact_tke = -0.5 / rsc_tke * rdt ! Cst used for the Diffusion term of tke1208 rfact_psi = -0.5 / rsc_psi * rdt ! Cst used for the Diffusion term of tke1209 1210 ! Wall proximity function1149 rsbc_psi2 = -0.5_wp * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking 1150 rsbc_psi3 = -0.5_wp * rdt * rc0**rpp * rl_sf**rnn / rsc_psi * (rnn + rmm*ra_sf) ! Neumann + Wave breaking 1151 rfact_tke = -0.5_wp / rsc_tke * rdt ! Cst used for the Diffusion term of tke 1152 rfact_psi = -0.5_wp / rsc_psi * rdt ! Cst used for the Diffusion term of tke 1153 1154 ! !* Wall proximity function 1211 1155 zwall (:,:,:) = 1._wp * tmask(:,:,:) 1212 1156 1213 ! !* set vertical eddy coef. to the background value1157 ! !* set vertical eddy coef. to the background value 1214 1158 DO jk = 1, jpk 1215 1159 avt (:,:,jk) = avtb(jk) * tmask(:,:,jk) … … 1218 1162 avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk) 1219 1163 END DO 1220 ! !* read or initialize all required files1164 ! !* read or initialize all required files 1221 1165 CALL gls_rst( nit000, 'READ' ) 1222 1166 ! … … 1236 1180 INTEGER , INTENT(in) :: kt ! ocean time-step 1237 1181 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 1238 ! !1182 ! 1239 1183 INTEGER :: jit, jk ! dummy loop indices 1240 1184 INTEGER :: id1, id2, id3, id4, id5, id6, id7, id8 … … 1327 1271 WRITE(*,*) 'zdf_gls: You should not have seen this print! error?', kt 1328 1272 END SUBROUTINE zdf_gls 1329 SUBROUTINE zdf_gls_init ! Empty routine 1330 END SUBROUTINE zdf_gls_init 1331 SUBROUTINE gls_rst( kt,cdrw ) ! Empty routine 1273 SUBROUTINE gls_rst( kt, cdrw ) ! Empty routine 1332 1274 INTEGER , INTENT(in) :: kt ! ocean time-step 1333 1275 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 1334 WRITE(*,*) 'gls_ Rst: You should not have seen this print! error?', kt, cdrw1276 WRITE(*,*) 'gls_rst: You should not have seen this print! error?', kt, cdrw 1335 1277 END SUBROUTINE gls_rst 1336 1278 #endif … … 1338 1280 !!====================================================================== 1339 1281 END MODULE zdfgls 1282
Note: See TracChangeset
for help on using the changeset viewer.
