- Timestamp:
- 2020-05-14T21:46:00+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/OCE/ZDF/zdfgls.F90 (modified) (25 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
- Property svn:externals
-
old new 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 9 # SETTE 10 ^/utils/CI/sette@HEAD sette
-
- Property svn:externals
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/ZDF/zdfgls.F90
r12178 r12928 104 104 105 105 !! * Substitutions 106 # include " vectopt_loop_substitute.h90"106 # include "do_loop_substitute.h90" 107 107 !!---------------------------------------------------------------------- 108 108 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 124 124 125 125 126 SUBROUTINE zdf_gls( kt, p_sh2, p_avm, p_avt )126 SUBROUTINE zdf_gls( kt, Kbb, Kmm, p_sh2, p_avm, p_avt ) 127 127 !!---------------------------------------------------------------------- 128 128 !! *** ROUTINE zdf_gls *** … … 134 134 !! 135 135 INTEGER , INTENT(in ) :: kt ! ocean time step 136 INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices 136 137 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: p_sh2 ! shear production term 137 138 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: p_avm, p_avt ! momentum and tracer Kz (w-points) … … 166 167 167 168 ! Compute surface, top and bottom friction at T-points 168 DO jj = 2, jpjm1 169 DO ji = fs_2, fs_jpim1 ! vector opt. 170 ! 171 ! surface friction 172 ustar2_surf(ji,jj) = r1_rau0 * taum(ji,jj) * tmask(ji,jj,1) 173 ! 169 DO_2D_00_00 170 ! 171 ! surface friction 172 ustar2_surf(ji,jj) = r1_rho0 * taum(ji,jj) * tmask(ji,jj,1) 173 ! 174 174 !!gm Rq we may add here r_ke0(_top/_bot) ? ==>> think about that... 175 ! bottom friction (explicit before friction) 176 zmsku = ( 2._wp - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 177 zmskv = ( 2._wp - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) ! (CAUTION: CdU<0) 178 ustar2_bot(ji,jj) = - rCdU_bot(ji,jj) * SQRT( ( zmsku*( ub(ji,jj,mbkt(ji,jj))+ub(ji-1,jj,mbkt(ji,jj)) ) )**2 & 179 & + ( zmskv*( vb(ji,jj,mbkt(ji,jj))+vb(ji,jj-1,mbkt(ji,jj)) ) )**2 ) 180 END DO 181 END DO 175 ! bottom friction (explicit before friction) 176 zmsku = ( 2._wp - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 177 zmskv = ( 2._wp - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) ! (CAUTION: CdU<0) 178 ustar2_bot(ji,jj) = - rCdU_bot(ji,jj) * SQRT( ( zmsku*( uu(ji,jj,mbkt(ji,jj),Kbb)+uu(ji-1,jj,mbkt(ji,jj),Kbb) ) )**2 & 179 & + ( zmskv*( vv(ji,jj,mbkt(ji,jj),Kbb)+vv(ji,jj-1,mbkt(ji,jj),Kbb) ) )**2 ) 180 END_2D 182 181 IF( ln_isfcav ) THEN !top friction 183 DO jj = 2, jpjm1 184 DO ji = fs_2, fs_jpim1 ! vector opt. 185 zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) ) 186 zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) ) ! (CAUTION: CdU<0) 187 ustar2_top(ji,jj) = - rCdU_top(ji,jj) * SQRT( ( zmsku*( ub(ji,jj,mikt(ji,jj))+ub(ji-1,jj,mikt(ji,jj)) ) )**2 & 188 & + ( zmskv*( vb(ji,jj,mikt(ji,jj))+vb(ji,jj-1,mikt(ji,jj)) ) )**2 ) 189 END DO 190 END DO 182 DO_2D_00_00 183 zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) ) 184 zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) ) ! (CAUTION: CdU<0) 185 ustar2_top(ji,jj) = - rCdU_top(ji,jj) * SQRT( ( zmsku*( uu(ji,jj,mikt(ji,jj),Kbb)+uu(ji-1,jj,mikt(ji,jj),Kbb) ) )**2 & 186 & + ( zmskv*( vv(ji,jj,mikt(ji,jj),Kbb)+vv(ji,jj-1,mikt(ji,jj),Kbb) ) )**2 ) 187 END_2D 191 188 ENDIF 192 189 … … 206 203 END SELECT 207 204 ! 208 DO jk = 2, jpkm1 !== Compute dissipation rate ==! 209 DO jj = 1, jpjm1 210 DO ji = 1, jpim1 211 eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / hmxl_n(ji,jj,jk) 212 END DO 213 END DO 214 END DO 205 DO_3D_10_10( 2, jpkm1 ) 206 eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / hmxl_n(ji,jj,jk) 207 END_3D 215 208 216 209 ! Save tke at before time step … … 219 212 220 213 IF( nn_clos == 0 ) THEN ! Mellor-Yamada 221 DO jk = 2, jpkm1 222 DO jj = 2, jpjm1 223 DO ji = fs_2, fs_jpim1 ! vector opt. 224 zup = hmxl_n(ji,jj,jk) * gdepw_n(ji,jj,mbkt(ji,jj)+1) 225 zdown = vkarmn * gdepw_n(ji,jj,jk) * ( -gdepw_n(ji,jj,jk) + gdepw_n(ji,jj,mbkt(ji,jj)+1) ) 226 zcoef = ( zup / MAX( zdown, rsmall ) ) 227 zwall (ji,jj,jk) = ( 1._wp + re2 * zcoef*zcoef ) * tmask(ji,jj,jk) 228 END DO 229 END DO 230 END DO 214 DO_3D_00_00( 2, jpkm1 ) 215 zup = hmxl_n(ji,jj,jk) * gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) 216 zdown = vkarmn * gdepw(ji,jj,jk,Kmm) * ( -gdepw(ji,jj,jk,Kmm) + gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) ) 217 zcoef = ( zup / MAX( zdown, rsmall ) ) 218 zwall (ji,jj,jk) = ( 1._wp + re2 * zcoef*zcoef ) * tmask(ji,jj,jk) 219 END_3D 231 220 ENDIF 232 221 … … 244 233 ! Warning : after this step, en : right hand side of the matrix 245 234 246 DO jk = 2, jpkm1 247 DO jj = 2, jpjm1 248 DO ji = 2, jpim1 249 ! 250 buoy = - p_avt(ji,jj,jk) * rn2(ji,jj,jk) ! stratif. destruction 251 ! 252 diss = eps(ji,jj,jk) ! dissipation 253 ! 254 zdir = 0.5_wp + SIGN( 0.5_wp, p_sh2(ji,jj,jk) + buoy ) ! zdir =1(=0) if shear(ji,jj,jk)+buoy >0(<0) 255 ! 256 zesh2 = zdir*(p_sh2(ji,jj,jk)+buoy)+(1._wp-zdir)*p_sh2(ji,jj,jk) ! production term 257 zdiss = zdir*(diss/en(ji,jj,jk)) +(1._wp-zdir)*(diss-buoy)/en(ji,jj,jk) ! dissipation term 235 DO_3D_00_00( 2, jpkm1 ) 236 ! 237 buoy = - p_avt(ji,jj,jk) * rn2(ji,jj,jk) ! stratif. destruction 238 ! 239 diss = eps(ji,jj,jk) ! dissipation 240 ! 241 zdir = 0.5_wp + SIGN( 0.5_wp, p_sh2(ji,jj,jk) + buoy ) ! zdir =1(=0) if shear(ji,jj,jk)+buoy >0(<0) 242 ! 243 zesh2 = zdir*(p_sh2(ji,jj,jk)+buoy)+(1._wp-zdir)*p_sh2(ji,jj,jk) ! production term 244 zdiss = zdir*(diss/en(ji,jj,jk)) +(1._wp-zdir)*(diss-buoy)/en(ji,jj,jk) ! dissipation term 258 245 !!gm better coding, identical results 259 246 ! zesh2 = p_sh2(ji,jj,jk) + zdir*buoy ! production term 260 247 ! zdiss = ( diss - (1._wp-zdir)*buoy ) / en(ji,jj,jk) ! dissipation term 261 248 !!gm 262 ! 263 ! Compute a wall function from 1. to rsc_psi*zwall/rsc_psi0 264 ! Note that as long that Dirichlet boundary conditions are NOT set at the first and last levels (GOTM style) 265 ! there is no need to set a boundary condition for zwall_psi at the top and bottom boundaries. 266 ! Otherwise, this should be rsc_psi/rsc_psi0 267 IF( ln_sigpsi ) THEN 268 zsigpsi = MIN( 1._wp, zesh2 / eps(ji,jj,jk) ) ! 0. <= zsigpsi <= 1. 269 zwall_psi(ji,jj,jk) = rsc_psi / & 270 & ( zsigpsi * rsc_psi + (1._wp-zsigpsi) * rsc_psi0 / MAX( zwall(ji,jj,jk), 1._wp ) ) 271 ELSE 272 zwall_psi(ji,jj,jk) = 1._wp 273 ENDIF 274 ! 275 ! building the matrix 276 zcof = rfact_tke * tmask(ji,jj,jk) 277 ! ! lower diagonal, in fact not used for jk = 2 (see surface conditions) 278 zd_lw(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk ) + p_avm(ji,jj,jk-1) ) / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk) ) 279 ! ! upper diagonal, in fact not used for jk = ibotm1 (see bottom conditions) 280 zd_up(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk ) ) / ( e3t_n(ji,jj,jk ) * e3w_n(ji,jj,jk) ) 281 ! ! diagonal 282 zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) + rdt * zdiss * wmask(ji,jj,jk) 283 ! ! right hand side in en 284 en(ji,jj,jk) = en(ji,jj,jk) + rdt * zesh2 * wmask(ji,jj,jk) 285 END DO 286 END DO 287 END DO 249 ! 250 ! Compute a wall function from 1. to rsc_psi*zwall/rsc_psi0 251 ! Note that as long that Dirichlet boundary conditions are NOT set at the first and last levels (GOTM style) 252 ! there is no need to set a boundary condition for zwall_psi at the top and bottom boundaries. 253 ! Otherwise, this should be rsc_psi/rsc_psi0 254 IF( ln_sigpsi ) THEN 255 zsigpsi = MIN( 1._wp, zesh2 / eps(ji,jj,jk) ) ! 0. <= zsigpsi <= 1. 256 zwall_psi(ji,jj,jk) = rsc_psi / & 257 & ( zsigpsi * rsc_psi + (1._wp-zsigpsi) * rsc_psi0 / MAX( zwall(ji,jj,jk), 1._wp ) ) 258 ELSE 259 zwall_psi(ji,jj,jk) = 1._wp 260 ENDIF 261 ! 262 ! building the matrix 263 zcof = rfact_tke * tmask(ji,jj,jk) 264 ! ! lower diagonal, in fact not used for jk = 2 (see surface conditions) 265 zd_lw(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk ) + p_avm(ji,jj,jk-1) ) / ( e3t(ji,jj,jk-1,Kmm) * e3w(ji,jj,jk,Kmm) ) 266 ! ! upper diagonal, in fact not used for jk = ibotm1 (see bottom conditions) 267 zd_up(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk ) ) / ( e3t(ji,jj,jk ,Kmm) * e3w(ji,jj,jk,Kmm) ) 268 ! ! diagonal 269 zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) + rn_Dt * zdiss * wmask(ji,jj,jk) 270 ! ! right hand side in en 271 en(ji,jj,jk) = en(ji,jj,jk) + rn_Dt * zesh2 * wmask(ji,jj,jk) 272 END_3D 288 273 ! 289 274 zdiag(:,:,jpk) = 1._wp … … 306 291 ! 307 292 ! One level below 308 en (:,:,2) = MAX( rc02r * ustar2_surf(:,:) * ( 1._wp + rsbc_tke1 * ((zhsro(:,:)+gdepw _n(:,:,2)) &293 en (:,:,2) = MAX( rc02r * ustar2_surf(:,:) * ( 1._wp + rsbc_tke1 * ((zhsro(:,:)+gdepw(:,:,2,Kmm)) & 309 294 & / zhsro(:,:) )**(1.5_wp*ra_sf) )**(2._wp/3._wp) , rn_emin ) 310 295 zd_lw(:,:,2) = 0._wp … … 325 310 zdiag(:,:,2) = zdiag(:,:,2) + zd_lw(:,:,2) ! Remove zd_lw from zdiag 326 311 zd_lw(:,:,2) = 0._wp 327 zkar (:,:) = (rl_sf + (vkarmn-rl_sf)*(1.-EXP(-rtrans*gdept _n(:,:,1)/zhsro(:,:)) ))312 zkar (:,:) = (rl_sf + (vkarmn-rl_sf)*(1.-EXP(-rtrans*gdept(:,:,1,Kmm)/zhsro(:,:)) )) 328 313 zflxs(:,:) = rsbc_tke2 * ustar2_surf(:,:)**1.5_wp * zkar(:,:) & 329 & * ( ( zhsro(:,:)+gdept _n(:,:,1) ) / zhsro(:,:) )**(1.5_wp*ra_sf)330 !!gm why not : * ( 1._wp + gdept _n(:,:,1) / zhsro(:,:) )**(1.5_wp*ra_sf)331 en(:,:,2) = en(:,:,2) + zflxs(:,:) / e3w _n(:,:,2)314 & * ( ( zhsro(:,:)+gdept(:,:,1,Kmm) ) / zhsro(:,:) )**(1.5_wp*ra_sf) 315 !!gm why not : * ( 1._wp + gdept(:,:,1,Kmm) / zhsro(:,:) )**(1.5_wp*ra_sf) 316 en(:,:,2) = en(:,:,2) + zflxs(:,:) / e3w(:,:,2,Kmm) 332 317 ! 333 318 ! … … 342 327 ! ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = rn_lmin 343 328 ! ! Balance between the production and the dissipation terms 344 DO jj = 2, jpjm1 345 DO ji = fs_2, fs_jpim1 ! vector opt. 329 DO_2D_00_00 346 330 !!gm This means that bottom and ocean w-level above have a specified "en" value. Sure ???? 347 331 !! With thick deep ocean level thickness, this may be quite large, no ??? 348 332 !! in particular in ocean cavities where top stratification can be large... 349 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 350 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 333 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 334 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 335 ! 336 z_en = MAX( rc02r * ustar2_bot(ji,jj), rn_emin ) 337 ! 338 ! Dirichlet condition applied at: 339 ! Bottom level (ibot) & Just above it (ibotm1) 340 zd_lw(ji,jj,ibot) = 0._wp ; zd_lw(ji,jj,ibotm1) = 0._wp 341 zd_up(ji,jj,ibot) = 0._wp ; zd_up(ji,jj,ibotm1) = 0._wp 342 zdiag(ji,jj,ibot) = 1._wp ; zdiag(ji,jj,ibotm1) = 1._wp 343 en (ji,jj,ibot) = z_en ; en (ji,jj,ibotm1) = z_en 344 END_2D 345 ! 346 IF( ln_isfcav) THEN ! top boundary (ocean cavity) 347 DO_2D_00_00 348 itop = mikt(ji,jj) ! k top w-point 349 itopp1 = mikt(ji,jj) + 1 ! k+1 1st w-point below the top one 350 ! ! mask at the ocean surface points 351 z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) ) 351 352 ! 352 z_en = MAX( rc02r * ustar2_bot(ji,jj), rn_emin ) 353 ! 353 !!gm TO BE VERIFIED !!! 354 354 ! Dirichlet condition applied at: 355 ! Bottom level (ibot) & Just above it (ibotm1) 356 zd_lw(ji,jj,ibot) = 0._wp ; zd_lw(ji,jj,ibotm1) = 0._wp 357 zd_up(ji,jj,ibot) = 0._wp ; zd_up(ji,jj,ibotm1) = 0._wp 358 zdiag(ji,jj,ibot) = 1._wp ; zdiag(ji,jj,ibotm1) = 1._wp 359 en (ji,jj,ibot) = z_en ; en (ji,jj,ibotm1) = z_en 360 END DO 361 END DO 362 ! 363 IF( ln_isfcav) THEN ! top boundary (ocean cavity) 364 DO jj = 2, jpjm1 365 DO ji = fs_2, fs_jpim1 ! vector opt. 366 itop = mikt(ji,jj) ! k top w-point 367 itopp1 = mikt(ji,jj) + 1 ! k+1 1st w-point below the top one 368 ! ! mask at the ocean surface points 369 z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) ) 370 ! 371 !!gm TO BE VERIFIED !!! 372 ! Dirichlet condition applied at: 373 ! top level (itop) & Just below it (itopp1) 374 zd_lw(ji,jj,itop) = 0._wp ; zd_lw(ji,jj,itopp1) = 0._wp 375 zd_up(ji,jj,itop) = 0._wp ; zd_up(ji,jj,itopp1) = 0._wp 376 zdiag(ji,jj,itop) = 1._wp ; zdiag(ji,jj,itopp1) = 1._wp 377 en (ji,jj,itop) = z_en ; en (ji,jj,itopp1) = z_en 378 END DO 379 END DO 355 ! top level (itop) & Just below it (itopp1) 356 zd_lw(ji,jj,itop) = 0._wp ; zd_lw(ji,jj,itopp1) = 0._wp 357 zd_up(ji,jj,itop) = 0._wp ; zd_up(ji,jj,itopp1) = 0._wp 358 zdiag(ji,jj,itop) = 1._wp ; zdiag(ji,jj,itopp1) = 1._wp 359 en (ji,jj,itop) = z_en ; en (ji,jj,itopp1) = z_en 360 END_2D 380 361 ENDIF 381 362 ! 382 363 CASE ( 1 ) ! Neumman boundary condition 383 364 ! 384 DO jj = 2, jpjm1 385 DO ji = fs_2, fs_jpim1 ! vector opt. 386 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 387 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 388 ! 389 z_en = MAX( rc02r * ustar2_bot(ji,jj), rn_emin ) 365 DO_2D_00_00 366 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 367 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 368 ! 369 z_en = MAX( rc02r * ustar2_bot(ji,jj), rn_emin ) 370 ! 371 ! Bottom level Dirichlet condition: 372 ! Bottom level (ibot) & Just above it (ibotm1) 373 ! Dirichlet ! Neumann 374 zd_lw(ji,jj,ibot) = 0._wp ! ! Remove zd_up from zdiag 375 zdiag(ji,jj,ibot) = 1._wp ; zdiag(ji,jj,ibotm1) = zdiag(ji,jj,ibotm1) + zd_up(ji,jj,ibotm1) 376 zd_up(ji,jj,ibot) = 0._wp ; zd_up(ji,jj,ibotm1) = 0._wp 377 END_2D 378 IF( ln_isfcav) THEN ! top boundary (ocean cavity) 379 DO_2D_00_00 380 itop = mikt(ji,jj) ! k top w-point 381 itopp1 = mikt(ji,jj) + 1 ! k+1 1st w-point below the top one 382 ! ! mask at the ocean surface points 383 z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) ) 390 384 ! 391 385 ! Bottom level Dirichlet condition: 392 386 ! Bottom level (ibot) & Just above it (ibotm1) 393 387 ! Dirichlet ! Neumann 394 zd_lw(ji,jj,ibot) = 0._wp ! ! Remove zd_up from zdiag 395 zdiag(ji,jj,ibot) = 1._wp ; zdiag(ji,jj,ibotm1) = zdiag(ji,jj,ibotm1) + zd_up(ji,jj,ibotm1) 396 zd_up(ji,jj,ibot) = 0._wp ; zd_up(ji,jj,ibotm1) = 0._wp 397 END DO 398 END DO 399 IF( ln_isfcav) THEN ! top boundary (ocean cavity) 400 DO jj = 2, jpjm1 401 DO ji = fs_2, fs_jpim1 ! vector opt. 402 itop = mikt(ji,jj) ! k top w-point 403 itopp1 = mikt(ji,jj) + 1 ! k+1 1st w-point below the top one 404 ! ! mask at the ocean surface points 405 z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) ) 406 ! 407 ! Bottom level Dirichlet condition: 408 ! Bottom level (ibot) & Just above it (ibotm1) 409 ! Dirichlet ! Neumann 410 zd_lw(ji,jj,itop) = 0._wp ! ! Remove zd_up from zdiag 411 zdiag(ji,jj,itop) = 1._wp ; zdiag(ji,jj,itopp1) = zdiag(ji,jj,itopp1) + zd_up(ji,jj,itopp1) 412 zd_up(ji,jj,itop) = 0._wp ; zd_up(ji,jj,itopp1) = 0._wp 413 END DO 414 END DO 388 zd_lw(ji,jj,itop) = 0._wp ! ! Remove zd_up from zdiag 389 zdiag(ji,jj,itop) = 1._wp ; zdiag(ji,jj,itopp1) = zdiag(ji,jj,itopp1) + zd_up(ji,jj,itopp1) 390 zd_up(ji,jj,itop) = 0._wp ; zd_up(ji,jj,itopp1) = 0._wp 391 END_2D 415 392 ENDIF 416 393 ! … … 420 397 ! ---------------------------------------------------------- 421 398 ! 422 DO jk = 2, jpkm1 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 423 DO jj = 2, jpjm1 424 DO ji = fs_2, fs_jpim1 ! vector opt. 425 zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 426 END DO 427 END DO 428 END DO 429 DO jk = 2, jpk ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 430 DO jj = 2, jpjm1 431 DO ji = fs_2, fs_jpim1 ! vector opt. 432 zd_lw(ji,jj,jk) = en(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) * zd_lw(ji,jj,jk-1) 433 END DO 434 END DO 435 END DO 436 DO jk = jpk-1, 2, -1 ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 437 DO jj = 2, jpjm1 438 DO ji = fs_2, fs_jpim1 ! vector opt. 439 en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 440 END DO 441 END DO 442 END DO 399 DO_3D_00_00( 2, jpkm1 ) 400 zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 401 END_3D 402 DO_3D_00_00( 2, jpk ) 403 zd_lw(ji,jj,jk) = en(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) * zd_lw(ji,jj,jk-1) 404 END_3D 405 DO_3DS_00_00( jpk-1, 2, -1 ) 406 en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 407 END_3D 443 408 ! ! set the minimum value of tke 444 409 en(:,:,:) = MAX( en(:,:,:), rn_emin ) … … 453 418 ! 454 419 CASE( 0 ) ! k-kl (Mellor-Yamada) 455 DO jk = 2, jpkm1 456 DO jj = 2, jpjm1 457 DO ji = fs_2, fs_jpim1 ! vector opt. 458 psi(ji,jj,jk) = eb(ji,jj,jk) * hmxl_b(ji,jj,jk) 459 END DO 460 END DO 461 END DO 420 DO_3D_00_00( 2, jpkm1 ) 421 psi(ji,jj,jk) = eb(ji,jj,jk) * hmxl_b(ji,jj,jk) 422 END_3D 462 423 ! 463 424 CASE( 1 ) ! k-eps 464 DO jk = 2, jpkm1 465 DO jj = 2, jpjm1 466 DO ji = fs_2, fs_jpim1 ! vector opt. 467 psi(ji,jj,jk) = eps(ji,jj,jk) 468 END DO 469 END DO 470 END DO 425 DO_3D_00_00( 2, jpkm1 ) 426 psi(ji,jj,jk) = eps(ji,jj,jk) 427 END_3D 471 428 ! 472 429 CASE( 2 ) ! k-w 473 DO jk = 2, jpkm1 474 DO jj = 2, jpjm1 475 DO ji = fs_2, fs_jpim1 ! vector opt. 476 psi(ji,jj,jk) = SQRT( eb(ji,jj,jk) ) / ( rc0 * hmxl_b(ji,jj,jk) ) 477 END DO 478 END DO 479 END DO 430 DO_3D_00_00( 2, jpkm1 ) 431 psi(ji,jj,jk) = SQRT( eb(ji,jj,jk) ) / ( rc0 * hmxl_b(ji,jj,jk) ) 432 END_3D 480 433 ! 481 434 CASE( 3 ) ! generic 482 DO jk = 2, jpkm1 483 DO jj = 2, jpjm1 484 DO ji = fs_2, fs_jpim1 ! vector opt. 485 psi(ji,jj,jk) = rc02 * eb(ji,jj,jk) * hmxl_b(ji,jj,jk)**rnn 486 END DO 487 END DO 488 END DO 435 DO_3D_00_00( 2, jpkm1 ) 436 psi(ji,jj,jk) = rc02 * eb(ji,jj,jk) * hmxl_b(ji,jj,jk)**rnn 437 END_3D 489 438 ! 490 439 END SELECT … … 497 446 ! Warning : after this step, en : right hand side of the matrix 498 447 499 DO jk = 2, jpkm1 500 DO jj = 2, jpjm1 501 DO ji = fs_2, fs_jpim1 ! vector opt. 502 ! 503 ! psi / k 504 zratio = psi(ji,jj,jk) / eb(ji,jj,jk) 505 ! 506 ! psi3+ : stable : B=-KhN²<0 => N²>0 if rn2>0 zdir = 1 (stable) otherwise zdir = 0 (unstable) 507 zdir = 0.5_wp + SIGN( 0.5_wp, rn2(ji,jj,jk) ) 508 ! 509 rpsi3 = zdir * rpsi3m + ( 1._wp - zdir ) * rpsi3p 510 ! 511 ! shear prod. - stratif. destruction 512 prod = rpsi1 * zratio * p_sh2(ji,jj,jk) 513 ! 514 ! stratif. destruction 515 buoy = rpsi3 * zratio * (- p_avt(ji,jj,jk) * rn2(ji,jj,jk) ) 516 ! 517 ! shear prod. - stratif. destruction 518 diss = rpsi2 * zratio * zwall(ji,jj,jk) * eps(ji,jj,jk) 519 ! 520 zdir = 0.5_wp + SIGN( 0.5_wp, prod + buoy ) ! zdir =1(=0) if shear(ji,jj,jk)+buoy >0(<0) 521 ! 522 zesh2 = zdir * ( prod + buoy ) + (1._wp - zdir ) * prod ! production term 523 zdiss = zdir * ( diss / psi(ji,jj,jk) ) + (1._wp - zdir ) * (diss-buoy) / psi(ji,jj,jk) ! dissipation term 524 ! 525 ! building the matrix 526 zcof = rfact_psi * zwall_psi(ji,jj,jk) * tmask(ji,jj,jk) 527 ! ! lower diagonal 528 zd_lw(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk ) + p_avm(ji,jj,jk-1) ) / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk) ) 529 ! ! upper diagonal 530 zd_up(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk ) ) / ( e3t_n(ji,jj,jk ) * e3w_n(ji,jj,jk) ) 531 ! ! diagonal 532 zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) + rdt * zdiss * wmask(ji,jj,jk) 533 ! ! right hand side in psi 534 psi(ji,jj,jk) = psi(ji,jj,jk) + rdt * zesh2 * wmask(ji,jj,jk) 535 END DO 536 END DO 537 END DO 448 DO_3D_00_00( 2, jpkm1 ) 449 ! 450 ! psi / k 451 zratio = psi(ji,jj,jk) / eb(ji,jj,jk) 452 ! 453 ! psi3+ : stable : B=-KhN²<0 => N²>0 if rn2>0 zdir = 1 (stable) otherwise zdir = 0 (unstable) 454 zdir = 0.5_wp + SIGN( 0.5_wp, rn2(ji,jj,jk) ) 455 ! 456 rpsi3 = zdir * rpsi3m + ( 1._wp - zdir ) * rpsi3p 457 ! 458 ! shear prod. - stratif. destruction 459 prod = rpsi1 * zratio * p_sh2(ji,jj,jk) 460 ! 461 ! stratif. destruction 462 buoy = rpsi3 * zratio * (- p_avt(ji,jj,jk) * rn2(ji,jj,jk) ) 463 ! 464 ! shear prod. - stratif. destruction 465 diss = rpsi2 * zratio * zwall(ji,jj,jk) * eps(ji,jj,jk) 466 ! 467 zdir = 0.5_wp + SIGN( 0.5_wp, prod + buoy ) ! zdir =1(=0) if shear(ji,jj,jk)+buoy >0(<0) 468 ! 469 zesh2 = zdir * ( prod + buoy ) + (1._wp - zdir ) * prod ! production term 470 zdiss = zdir * ( diss / psi(ji,jj,jk) ) + (1._wp - zdir ) * (diss-buoy) / psi(ji,jj,jk) ! dissipation term 471 ! 472 ! building the matrix 473 zcof = rfact_psi * zwall_psi(ji,jj,jk) * tmask(ji,jj,jk) 474 ! ! lower diagonal 475 zd_lw(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk ) + p_avm(ji,jj,jk-1) ) / ( e3t(ji,jj,jk-1,Kmm) * e3w(ji,jj,jk,Kmm) ) 476 ! ! upper diagonal 477 zd_up(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk ) ) / ( e3t(ji,jj,jk ,Kmm) * e3w(ji,jj,jk,Kmm) ) 478 ! ! diagonal 479 zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) + rn_Dt * zdiss * wmask(ji,jj,jk) 480 ! ! right hand side in psi 481 psi(ji,jj,jk) = psi(ji,jj,jk) + rn_Dt * zesh2 * wmask(ji,jj,jk) 482 END_3D 538 483 ! 539 484 zdiag(:,:,jpk) = 1._wp … … 554 499 ! 555 500 ! One level below 556 zkar (:,:) = (rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdepw _n(:,:,2)/zhsro(:,:) )))557 zdep (:,:) = (zhsro(:,:) + gdepw _n(:,:,2)) * zkar(:,:)501 zkar (:,:) = (rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdepw(:,:,2,Kmm)/zhsro(:,:) ))) 502 zdep (:,:) = (zhsro(:,:) + gdepw(:,:,2,Kmm)) * zkar(:,:) 558 503 psi (:,:,2) = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 559 504 zd_lw(:,:,2) = 0._wp … … 575 520 ! 576 521 ! Set psi vertical flux at the surface: 577 zkar (:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdept _n(:,:,1)/zhsro(:,:) )) ! Lengh scale slope578 zdep (:,:) = ((zhsro(:,:) + gdept _n(:,:,1)) / zhsro(:,:))**(rmm*ra_sf)522 zkar (:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdept(:,:,1,Kmm)/zhsro(:,:) )) ! Lengh scale slope 523 zdep (:,:) = ((zhsro(:,:) + gdept(:,:,1,Kmm)) / zhsro(:,:))**(rmm*ra_sf) 579 524 zflxs(:,:) = (rnn + rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:))*(1._wp + rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 580 525 zdep (:,:) = rsbc_psi1 * (zwall_psi(:,:,1)*p_avm(:,:,1)+zwall_psi(:,:,2)*p_avm(:,:,2)) * & 581 & ustar2_surf(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + gdept _n(:,:,1))**(rnn-1.)526 & ustar2_surf(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + gdept(:,:,1,Kmm))**(rnn-1.) 582 527 zflxs(:,:) = zdep(:,:) * zflxs(:,:) 583 psi (:,:,2) = psi(:,:,2) + zflxs(:,:) / e3w _n(:,:,2)528 psi (:,:,2) = psi(:,:,2) + zflxs(:,:) / e3w(:,:,2,Kmm) 584 529 ! 585 530 END SELECT … … 596 541 ! ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = vkarmn * r_z0_bot 597 542 ! ! Balance between the production and the dissipation terms 598 DO jj = 2, jpjm1 599 DO ji = fs_2, fs_jpim1 ! vector opt. 600 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 601 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 602 zdep(ji,jj) = vkarmn * r_z0_bot 603 psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn 604 zd_lw(ji,jj,ibot) = 0._wp 605 zd_up(ji,jj,ibot) = 0._wp 606 zdiag(ji,jj,ibot) = 1._wp 607 ! 608 ! Just above last level, Dirichlet condition again (GOTM like) 609 zdep(ji,jj) = vkarmn * ( r_z0_bot + e3t_n(ji,jj,ibotm1) ) 610 psi (ji,jj,ibotm1) = rc0**rpp * en(ji,jj,ibot )**rmm * zdep(ji,jj)**rnn 611 zd_lw(ji,jj,ibotm1) = 0._wp 612 zd_up(ji,jj,ibotm1) = 0._wp 613 zdiag(ji,jj,ibotm1) = 1._wp 614 END DO 615 END DO 543 DO_2D_00_00 544 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 545 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 546 zdep(ji,jj) = vkarmn * r_z0_bot 547 psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn 548 zd_lw(ji,jj,ibot) = 0._wp 549 zd_up(ji,jj,ibot) = 0._wp 550 zdiag(ji,jj,ibot) = 1._wp 551 ! 552 ! Just above last level, Dirichlet condition again (GOTM like) 553 zdep(ji,jj) = vkarmn * ( r_z0_bot + e3t(ji,jj,ibotm1,Kmm) ) 554 psi (ji,jj,ibotm1) = rc0**rpp * en(ji,jj,ibot )**rmm * zdep(ji,jj)**rnn 555 zd_lw(ji,jj,ibotm1) = 0._wp 556 zd_up(ji,jj,ibotm1) = 0._wp 557 zdiag(ji,jj,ibotm1) = 1._wp 558 END_2D 616 559 ! 617 560 CASE ( 1 ) ! Neumman boundary condition 618 561 ! 619 DO jj = 2, jpjm1 620 DO ji = fs_2, fs_jpim1 ! vector opt. 621 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 622 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 623 ! 624 ! Bottom level Dirichlet condition: 625 zdep(ji,jj) = vkarmn * r_z0_bot 626 psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn 627 ! 628 zd_lw(ji,jj,ibot) = 0._wp 629 zd_up(ji,jj,ibot) = 0._wp 630 zdiag(ji,jj,ibot) = 1._wp 631 ! 632 ! Just above last level: Neumann condition with flux injection 633 zdiag(ji,jj,ibotm1) = zdiag(ji,jj,ibotm1) + zd_up(ji,jj,ibotm1) ! Remove zd_up from zdiag 634 zd_up(ji,jj,ibotm1) = 0. 635 ! 636 ! Set psi vertical flux at the bottom: 637 zdep(ji,jj) = r_z0_bot + 0.5_wp*e3t_n(ji,jj,ibotm1) 638 zflxb = rsbc_psi2 * ( p_avm(ji,jj,ibot) + p_avm(ji,jj,ibotm1) ) & 639 & * (0.5_wp*(en(ji,jj,ibot)+en(ji,jj,ibotm1)))**rmm * zdep(ji,jj)**(rnn-1._wp) 640 psi(ji,jj,ibotm1) = psi(ji,jj,ibotm1) + zflxb / e3w_n(ji,jj,ibotm1) 641 END DO 642 END DO 562 DO_2D_00_00 563 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 564 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 565 ! 566 ! Bottom level Dirichlet condition: 567 zdep(ji,jj) = vkarmn * r_z0_bot 568 psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn 569 ! 570 zd_lw(ji,jj,ibot) = 0._wp 571 zd_up(ji,jj,ibot) = 0._wp 572 zdiag(ji,jj,ibot) = 1._wp 573 ! 574 ! Just above last level: Neumann condition with flux injection 575 zdiag(ji,jj,ibotm1) = zdiag(ji,jj,ibotm1) + zd_up(ji,jj,ibotm1) ! Remove zd_up from zdiag 576 zd_up(ji,jj,ibotm1) = 0. 577 ! 578 ! Set psi vertical flux at the bottom: 579 zdep(ji,jj) = r_z0_bot + 0.5_wp*e3t(ji,jj,ibotm1,Kmm) 580 zflxb = rsbc_psi2 * ( p_avm(ji,jj,ibot) + p_avm(ji,jj,ibotm1) ) & 581 & * (0.5_wp*(en(ji,jj,ibot)+en(ji,jj,ibotm1)))**rmm * zdep(ji,jj)**(rnn-1._wp) 582 psi(ji,jj,ibotm1) = psi(ji,jj,ibotm1) + zflxb / e3w(ji,jj,ibotm1,Kmm) 583 END_2D 643 584 ! 644 585 END SELECT … … 647 588 ! ---------------- 648 589 ! 649 DO jk = 2, jpkm1 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 650 DO jj = 2, jpjm1 651 DO ji = fs_2, fs_jpim1 ! vector opt. 652 zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 653 END DO 654 END DO 655 END DO 656 DO jk = 2, jpk ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 657 DO jj = 2, jpjm1 658 DO ji = fs_2, fs_jpim1 ! vector opt. 659 zd_lw(ji,jj,jk) = psi(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) * zd_lw(ji,jj,jk-1) 660 END DO 661 END DO 662 END DO 663 DO jk = jpk-1, 2, -1 ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 664 DO jj = 2, jpjm1 665 DO ji = fs_2, fs_jpim1 ! vector opt. 666 psi(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * psi(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 667 END DO 668 END DO 669 END DO 590 DO_3D_00_00( 2, jpkm1 ) 591 zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 592 END_3D 593 DO_3D_00_00( 2, jpk ) 594 zd_lw(ji,jj,jk) = psi(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) * zd_lw(ji,jj,jk-1) 595 END_3D 596 DO_3DS_00_00( jpk-1, 2, -1 ) 597 psi(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * psi(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 598 END_3D 670 599 671 600 ! Set dissipation … … 675 604 ! 676 605 CASE( 0 ) ! k-kl (Mellor-Yamada) 677 DO jk = 1, jpkm1 678 DO jj = 2, jpjm1 679 DO ji = fs_2, fs_jpim1 ! vector opt. 680 eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / MAX( psi(ji,jj,jk), rn_epsmin) 681 END DO 682 END DO 683 END DO 606 DO_3D_00_00( 1, jpkm1 ) 607 eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / MAX( psi(ji,jj,jk), rn_epsmin) 608 END_3D 684 609 ! 685 610 CASE( 1 ) ! k-eps 686 DO jk = 1, jpkm1 687 DO jj = 2, jpjm1 688 DO ji = fs_2, fs_jpim1 ! vector opt. 689 eps(ji,jj,jk) = psi(ji,jj,jk) 690 END DO 691 END DO 692 END DO 611 DO_3D_00_00( 1, jpkm1 ) 612 eps(ji,jj,jk) = psi(ji,jj,jk) 613 END_3D 693 614 ! 694 615 CASE( 2 ) ! k-w 695 DO jk = 1, jpkm1 696 DO jj = 2, jpjm1 697 DO ji = fs_2, fs_jpim1 ! vector opt. 698 eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk) 699 END DO 700 END DO 701 END DO 616 DO_3D_00_00( 1, jpkm1 ) 617 eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk) 618 END_3D 702 619 ! 703 620 CASE( 3 ) ! generic … … 705 622 zex1 = ( 1.5_wp + rmm/rnn ) 706 623 zex2 = -1._wp / rnn 707 DO jk = 1, jpkm1 708 DO jj = 2, jpjm1 709 DO ji = fs_2, fs_jpim1 ! vector opt. 710 eps(ji,jj,jk) = zcoef * en(ji,jj,jk)**zex1 * psi(ji,jj,jk)**zex2 711 END DO 712 END DO 713 END DO 624 DO_3D_00_00( 1, jpkm1 ) 625 eps(ji,jj,jk) = zcoef * en(ji,jj,jk)**zex1 * psi(ji,jj,jk)**zex2 626 END_3D 714 627 ! 715 628 END SELECT … … 717 630 ! Limit dissipation rate under stable stratification 718 631 ! -------------------------------------------------- 719 DO jk = 1, jpkm1 ! Note that this set boundary conditions on hmxl_n at the same time 720 DO jj = 2, jpjm1 721 DO ji = fs_2, fs_jpim1 ! vector opt. 722 ! limitation 723 eps (ji,jj,jk) = MAX( eps(ji,jj,jk), rn_epsmin ) 724 hmxl_n(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk) 725 ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated) 726 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 727 IF( ln_length_lim ) hmxl_n(ji,jj,jk) = MIN( rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), hmxl_n(ji,jj,jk) ) 728 END DO 729 END DO 730 END DO 632 DO_3D_00_00( 1, jpkm1 ) 633 ! limitation 634 eps (ji,jj,jk) = MAX( eps(ji,jj,jk), rn_epsmin ) 635 hmxl_n(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk) 636 ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated) 637 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 638 IF( ln_length_lim ) hmxl_n(ji,jj,jk) = MIN( rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), hmxl_n(ji,jj,jk) ) 639 END_3D 731 640 732 641 ! … … 737 646 ! 738 647 CASE ( 0 , 1 ) ! Galperin or Kantha-Clayson stability functions 739 DO jk = 2, jpkm1 740 DO jj = 2, jpjm1 741 DO ji = fs_2, fs_jpim1 ! vector opt. 742 ! zcof = l²/q² 743 zcof = hmxl_b(ji,jj,jk) * hmxl_b(ji,jj,jk) / ( 2._wp*eb(ji,jj,jk) ) 744 ! Gh = -N²l²/q² 745 gh = - rn2(ji,jj,jk) * zcof 746 gh = MIN( gh, rgh0 ) 747 gh = MAX( gh, rghmin ) 748 ! Stability functions from Kantha and Clayson (if C2=C3=0 => Galperin) 749 sh = ra2*( 1._wp-6._wp*ra1/rb1 ) / ( 1.-3.*ra2*gh*(6.*ra1+rb2*( 1._wp-rc3 ) ) ) 750 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) 751 ! 752 ! Store stability function in zstt and zstm 753 zstt(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk) 754 zstm(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk) 755 END DO 756 END DO 757 END DO 648 DO_3D_00_00( 2, jpkm1 ) 649 ! zcof = l²/q² 650 zcof = hmxl_b(ji,jj,jk) * hmxl_b(ji,jj,jk) / ( 2._wp*eb(ji,jj,jk) ) 651 ! Gh = -N²l²/q² 652 gh = - rn2(ji,jj,jk) * zcof 653 gh = MIN( gh, rgh0 ) 654 gh = MAX( gh, rghmin ) 655 ! Stability functions from Kantha and Clayson (if C2=C3=0 => Galperin) 656 sh = ra2*( 1._wp-6._wp*ra1/rb1 ) / ( 1.-3.*ra2*gh*(6.*ra1+rb2*( 1._wp-rc3 ) ) ) 657 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) 658 ! 659 ! Store stability function in zstt and zstm 660 zstt(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk) 661 zstm(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk) 662 END_3D 758 663 ! 759 664 CASE ( 2, 3 ) ! Canuto stability functions 760 DO jk = 2, jpkm1 761 DO jj = 2, jpjm1 762 DO ji = fs_2, fs_jpim1 ! vector opt. 763 ! zcof = l²/q² 764 zcof = hmxl_b(ji,jj,jk)*hmxl_b(ji,jj,jk) / ( 2._wp * eb(ji,jj,jk) ) 765 ! Gh = -N²l²/q² 766 gh = - rn2(ji,jj,jk) * zcof 767 gh = MIN( gh, rgh0 ) 768 gh = MAX( gh, rghmin ) 769 gh = gh * rf6 770 ! Gm = M²l²/q² Shear number 771 shr = p_sh2(ji,jj,jk) / MAX( p_avm(ji,jj,jk), rsmall ) 772 gm = MAX( shr * zcof , 1.e-10 ) 773 gm = gm * rf6 774 gm = MIN ( (rd0 - rd1*gh + rd3*gh*gh) / (rd2-rd4*gh) , gm ) 775 ! Stability functions from Canuto 776 rcff = rd0 - rd1*gh +rd2*gm + rd3*gh*gh - rd4*gh*gm + rd5*gm*gm 777 sm = (rs0 - rs1*gh + rs2*gm) / rcff 778 sh = (rs4 - rs5*gh + rs6*gm) / rcff 779 ! 780 ! Store stability function in zstt and zstm 781 zstt(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk) 782 zstm(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk) 783 END DO 784 END DO 785 END DO 665 DO_3D_00_00( 2, jpkm1 ) 666 ! zcof = l²/q² 667 zcof = hmxl_b(ji,jj,jk)*hmxl_b(ji,jj,jk) / ( 2._wp * eb(ji,jj,jk) ) 668 ! Gh = -N²l²/q² 669 gh = - rn2(ji,jj,jk) * zcof 670 gh = MIN( gh, rgh0 ) 671 gh = MAX( gh, rghmin ) 672 gh = gh * rf6 673 ! Gm = M²l²/q² Shear number 674 shr = p_sh2(ji,jj,jk) / MAX( p_avm(ji,jj,jk), rsmall ) 675 gm = MAX( shr * zcof , 1.e-10 ) 676 gm = gm * rf6 677 gm = MIN ( (rd0 - rd1*gh + rd3*gh*gh) / (rd2-rd4*gh) , gm ) 678 ! Stability functions from Canuto 679 rcff = rd0 - rd1*gh +rd2*gm + rd3*gh*gh - rd4*gh*gm + rd5*gm*gm 680 sm = (rs0 - rs1*gh + rs2*gm) / rcff 681 sh = (rs4 - rs5*gh + rs6*gm) / rcff 682 ! 683 ! Store stability function in zstt and zstm 684 zstt(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk) 685 zstm(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk) 686 END_3D 786 687 ! 787 688 END SELECT … … 794 695 ! default value, in case jpk > mbkt(ji,jj)+1. Not needed but avoid a bug when looking for undefined values (-fpe0) 795 696 zstm(:,:,jpk) = 0. 796 DO jj = 2, jpjm1 ! update bottom with good values 797 DO ji = fs_2, fs_jpim1 ! vector opt. 798 zstm(ji,jj,mbkt(ji,jj)+1) = zstm(ji,jj,mbkt(ji,jj)) 799 END DO 800 END DO 697 DO_2D_00_00 698 zstm(ji,jj,mbkt(ji,jj)+1) = zstm(ji,jj,mbkt(ji,jj)) 699 END_2D 801 700 802 701 zstt(:,:, 1) = wmask(:,:, 1) ! default value not needed but avoid a bug when looking for undefined values (-fpe0) … … 811 710 ! later overwritten by surface/bottom boundaries conditions, so we don't really care of p_avm(:,:1) and p_avm(:,:jpk) 812 711 ! for zd_lw and zd_up but they have to be defined to avoid a bug when looking for undefined values (-fpe0) 813 DO jk = 1, jpk 814 DO jj = 2, jpjm1 815 DO ji = fs_2, fs_jpim1 ! vector opt. 816 zsqen = SQRT( 2._wp * en(ji,jj,jk) ) * hmxl_n(ji,jj,jk) 817 zavt = zsqen * zstt(ji,jj,jk) 818 zavm = zsqen * zstm(ji,jj,jk) 819 p_avt(ji,jj,jk) = MAX( zavt, avtb(jk) ) * wmask(ji,jj,jk) ! apply mask for zdfmxl routine 820 p_avm(ji,jj,jk) = MAX( zavm, avmb(jk) ) ! Note that avm is not masked at the surface and the bottom 821 END DO 822 END DO 823 END DO 712 DO_3D_00_00( 1, jpk ) 713 zsqen = SQRT( 2._wp * en(ji,jj,jk) ) * hmxl_n(ji,jj,jk) 714 zavt = zsqen * zstt(ji,jj,jk) 715 zavm = zsqen * zstm(ji,jj,jk) 716 p_avt(ji,jj,jk) = MAX( zavt, avtb(jk) ) * wmask(ji,jj,jk) ! apply mask for zdfmxl routine 717 p_avm(ji,jj,jk) = MAX( zavm, avmb(jk) ) ! Note that avm is not masked at the surface and the bottom 718 END_3D 824 719 p_avt(:,:,1) = 0._wp 825 720 ! 826 IF( ln_ctl) THEN721 IF(sn_cfctl%l_prtctl) THEN 827 722 CALL prt_ctl( tab3d_1=en , clinfo1=' gls - e: ', tab3d_2=p_avt, clinfo2=' t: ', kdim=jpk) 828 723 CALL prt_ctl( tab3d_1=p_avm, clinfo1=' gls - m: ', kdim=jpk ) … … 857 752 !!---------------------------------------------------------- 858 753 ! 859 REWIND( numnam_ref ) ! Namelist namzdf_gls in reference namelist : Vertical eddy diffivity and viscosity using gls turbulent closure scheme860 754 READ ( numnam_ref, namzdf_gls, IOSTAT = ios, ERR = 901) 861 755 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_gls in reference namelist' ) 862 756 863 REWIND( numnam_cfg ) ! Namelist namzdf_gls in configuration namelist : Vertical eddy diffivity and viscosity using gls turbulent closure scheme864 757 READ ( numnam_cfg, namzdf_gls, IOSTAT = ios, ERR = 902 ) 865 758 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namzdf_gls in configuration namelist' ) … … 1114 1007 rc04 = rc03 * rc0 1115 1008 rsbc_tke1 = -3._wp/2._wp*rn_crban*ra_sf*rl_sf ! Dirichlet + Wave breaking 1116 rsbc_tke2 = r dt * rn_crban / rl_sf ! Neumann + Wave breaking1009 rsbc_tke2 = rn_Dt * rn_crban / rl_sf ! Neumann + Wave breaking 1117 1010 zcr = MAX(rsmall, rsbc_tke1**(1./(-ra_sf*3._wp/2._wp))-1._wp ) 1118 1011 rtrans = 0.2_wp / zcr ! Ad. inverse transition length between log and wave layer 1119 1012 rsbc_zs1 = rn_charn/grav ! Charnock formula for surface roughness 1120 1013 rsbc_zs2 = rn_frac_hs / 0.85_wp / grav * 665._wp ! Rascle formula for surface roughness 1121 rsbc_psi1 = -0.5_wp * r dt * rc0**(rpp-2._wp*rmm) / rsc_psi1122 rsbc_psi2 = -0.5_wp * r dt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking1123 ! 1124 rfact_tke = -0.5_wp / rsc_tke * r dt ! Cst used for the Diffusion term of tke1125 rfact_psi = -0.5_wp / rsc_psi * r dt ! Cst used for the Diffusion term of tke1014 rsbc_psi1 = -0.5_wp * rn_Dt * rc0**(rpp-2._wp*rmm) / rsc_psi 1015 rsbc_psi2 = -0.5_wp * rn_Dt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking 1016 ! 1017 rfact_tke = -0.5_wp / rsc_tke * rn_Dt ! Cst used for the Diffusion term of tke 1018 rfact_psi = -0.5_wp / rsc_psi * rn_Dt ! Cst used for the Diffusion term of tke 1126 1019 ! 1127 1020 ! !* Wall proximity function
Note: See TracChangeset
for help on using the changeset viewer.
