- Timestamp:
- 2017-06-06T15:55:44+02:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90
r8093 r8143 8 8 !! 3.3 ! 2010-10 (C. Bricaud) Add in the reference 9 9 !! 4.0 ! 2017-04 (G. Madec) remove CPP keys & avm at t-point only 10 !! - ! 2017-05 (G. Madec) add top friction as boundary condition 10 11 !!---------------------------------------------------------------------- 11 12 … … 19 20 USE domvvl ! ocean space and time domain : variable volume layer 20 21 USE zdf_oce ! ocean vertical physics 21 !!gm old22 USE zdfbfr , ONLY : rn_tfrz0, rn_bfrz0 ! top/bottom roughness23 !!gm new24 22 USE zdfdrg , ONLY : r_z0_top , r_z0_bot ! top/bottom roughness 25 23 USE zdfdrg , ONLY : rCdU_top , rCdU_bot ! top/bottom friction 26 !!gm27 24 USE sbc_oce ! surface boundary condition: ocean 28 25 USE phycst ! physical constants … … 32 29 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 33 30 USE lib_mpp ! MPP manager 34 USE wrk_nemo ! work arrays35 31 USE prtctl ! Print control 36 32 USE in_out_manager ! I/O manager … … 149 145 REAL(wp) :: prod, buoy, diss, zdiss, sm ! - - 150 146 REAL(wp) :: gh, gm, shr, dif, zsqen, zavt, zavm ! - - 147 REAL(wp) :: zmsku, zmskv ! - - 151 148 REAL(wp), DIMENSION(jpi,jpj) :: zdep 152 149 REAL(wp), DIMENSION(jpi,jpj) :: zkar … … 158 155 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwall_psi ! Wall function use in the wb case (ln_sigpsi) 159 156 REAL(wp), DIMENSION(jpi,jpj,jpk) :: psi ! psi at time now 160 REAL(wp), DIMENSION(jpi,jpj,jpk) :: z_elem_a ! element of the first matrix diagonal 161 REAL(wp), DIMENSION(jpi,jpj,jpk) :: z_elem_b ! element of the second matrix diagonal 162 REAL(wp), DIMENSION(jpi,jpj,jpk) :: z_elem_c ! element of the third matrix diagonal 157 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zd_lw, zd_up, zdiag ! lower, upper and diagonal of the matrix 163 158 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zstt, zstm ! stability function on tracer and momentum 164 159 !!-------------------------------------------------------------------- … … 171 166 ustar2_top (:,:) = 0._wp ; zwall_psi(:,:,:) = 0._wp 172 167 ustar2_bot (:,:) = 0._wp 173 174 175 168 176 169 ! Compute surface, top and bottom friction at T-points … … 181 174 ustar2_surf(ji,jj) = r1_rau0 * taum(ji,jj) * tmask(ji,jj,1) 182 175 ! 183 !!gm old184 ! bottom friction (explicit before friction)185 ! Note that we chose here not to bound the friction as in dynbfr)186 ztx2 = ( bfrua(ji,jj) * ub(ji,jj,mbku(ji,jj)) + bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj)) ) &187 & * ( 1._wp - 0.5_wp * umask(ji,jj,1) * umask(ji-1,jj,1) )188 zty2 = ( bfrva(ji,jj) * vb(ji,jj,mbkv(ji,jj)) + bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1)) ) &189 & * ( 1._wp - 0.5_wp * vmask(ji,jj,1) * vmask(ji,jj-1,1) )190 ustar2_bot(ji,jj) = SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)191 END DO192 END DO193 !!gm new194 176 !!gm Rq we may add here r_ke0(_top/_bot) ? ==>> think about that... 195 !! bottom friction (explicit before friction)196 !zmsku = ( 2._wp - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) )197 !zmskv = ( 2._wp - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) ! (CAUTION: CdU<0)198 ! ustar2_bot(ji,jj) = - rCdU_bot(ji,jj) * SQRT( ( zmsku*( ub(ji,jj,mbkt(ji,jj))+ub(ji-1,jj,mbkt(ji,jj)) ) )**2 199 !& + ( zmskv*( vb(ji,jj,mbkt(ji,jj))+vb(ji,jj-1,mbkt(ji,jj)) ) )**2 )200 !END DO201 !END DO202 !IF( ln_isfcav ) THEN !top friction203 !DO jj = 2, jpjm1204 !DO ji = fs_2, fs_jpim1 ! vector opt.205 !zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) )206 !zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) ) ! (CAUTION: CdU<0)207 ! ustar2_top(ji,jj) = - rCdU_top(ji,jj) * SQRT( ( zmsku*( ub(ji,jj,mikt(ji,jj))+ub(ji-1,jj,mikt(ji,jj)) ) )**2 208 !& + ( zmskv*( vb(ji,jj,mikt(ji,jj))+vb(ji,jj-1,mikt(ji,jj)) ) )**2 )209 !END DO210 !END DO211 !ENDIF212 !!gm 177 ! bottom friction (explicit before friction) 178 zmsku = ( 2._wp - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 179 zmskv = ( 2._wp - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) ! (CAUTION: CdU<0) 180 ustar2_bot(ji,jj) = - rCdU_bot(ji,jj) * SQRT( ( zmsku*( ub(ji,jj,mbkt(ji,jj))+ub(ji-1,jj,mbkt(ji,jj)) ) )**2 & 181 & + ( zmskv*( vb(ji,jj,mbkt(ji,jj))+vb(ji,jj-1,mbkt(ji,jj)) ) )**2 ) 182 END DO 183 END DO 184 IF( ln_isfcav ) THEN !top friction 185 DO jj = 2, jpjm1 186 DO ji = fs_2, fs_jpim1 ! vector opt. 187 zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) ) 188 zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) ) ! (CAUTION: CdU<0) 189 ustar2_top(ji,jj) = - rCdU_top(ji,jj) * SQRT( ( zmsku*( ub(ji,jj,mikt(ji,jj))+ub(ji-1,jj,mikt(ji,jj)) ) )**2 & 190 & + ( zmskv*( vb(ji,jj,mikt(ji,jj))+vb(ji,jj-1,mikt(ji,jj)) ) )**2 ) 191 END DO 192 END DO 193 ENDIF 194 213 195 SELECT CASE ( nn_z0_met ) !== Set surface roughness length ==! 214 196 CASE ( 0 ) ! Constant roughness … … 261 243 ! The surface boundary condition are set after 262 244 ! The bottom boundary condition are also set after. In standard e(bottom)=0. 263 ! z _elem_b : diagonal z_elem_c : upper diagonal z_elem_a: lower diagonal245 ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal 264 246 ! Warning : after this step, en : right hand side of the matrix 265 247 … … 296 278 zcof = rfact_tke * tmask(ji,jj,jk) 297 279 ! ! lower diagonal 298 z _elem_a(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) )280 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) ) 299 281 ! ! upper diagonal 300 z _elem_c(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) )282 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) ) 301 283 ! ! diagonal 302 z _elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) + rdt * zdiss * wmask(ji,jj,jk)284 zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) + rdt * zdiss * wmask(ji,jj,jk) 303 285 ! ! right hand side in en 304 286 en(ji,jj,jk) = en(ji,jj,jk) + rdt * zesh2 * wmask(ji,jj,jk) … … 307 289 END DO 308 290 ! 309 z _elem_b(:,:,jpk) = 1._wp291 zdiag(:,:,jpk) = 1._wp 310 292 ! 311 293 ! Set surface condition on zwall_psi (1 at the bottom) … … 318 300 SELECT CASE ( nn_bc_surf ) 319 301 ! 320 CASE ( 0 ) ! Dirichlet case302 CASE ( 0 ) ! Dirichlet boundary condition (set e at k=1 & 2) 321 303 ! First level 322 en(:,:,1) = rc02r * ustar2_surf(:,:) * (1._wp + rsbc_tke1)**r2_3 323 en(:,:,1) = MAX(en(:,:,1), rn_emin) 324 z_elem_a(:,:,1) = en(:,:,1) 325 z_elem_c(:,:,1) = 0._wp 326 z_elem_b(:,:,1) = 1._wp 304 en (:,:,1) = MAX( rn_emin , rc02r * ustar2_surf(:,:) * (1._wp + rsbc_tke1)**r2_3 ) 305 zd_lw(:,:,1) = en(:,:,1) 306 zd_up(:,:,1) = 0._wp 307 zdiag(:,:,1) = 1._wp 327 308 ! 328 309 ! One level below 329 en(:,:,2) = rc02r * ustar2_surf(:,:) * (1._wp + rsbc_tke1 * ((zhsro(:,:)+gdepw_n(:,:,2)) & 330 & / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp) 331 en(:,:,2) = MAX(en(:,:,2), rn_emin ) 332 z_elem_a(:,:,2) = 0._wp 333 z_elem_c(:,:,2) = 0._wp 334 z_elem_b(:,:,2) = 1._wp 335 ! 336 ! 337 CASE ( 1 ) ! Neumann boundary condition on d(e)/dz 310 en (:,:,2) = MAX( rc02r * ustar2_surf(:,:) * ( 1._wp + rsbc_tke1 * ((zhsro(:,:)+gdepw_n(:,:,2)) & 311 & / zhsro(:,:) )**(1.5_wp*ra_sf) )**(2._wp/3._wp) , rn_emin ) 312 zd_lw(:,:,2) = 0._wp 313 zd_up(:,:,2) = 0._wp 314 zdiag(:,:,2) = 1._wp 315 ! 316 ! 317 CASE ( 1 ) ! Neumann boundary condition (set d(e)/dz) 338 318 ! 339 319 ! Dirichlet conditions at k=1 340 en(:,:,1) = rc02r * ustar2_surf(:,:) * (1._wp + rsbc_tke1)**r2_3 341 en(:,:,1) = MAX(en(:,:,1), rn_emin) 342 z_elem_a(:,:,1) = en(:,:,1) 343 z_elem_c(:,:,1) = 0._wp 344 z_elem_b(:,:,1) = 1._wp 320 en (:,:,1) = MAX( rc02r * ustar2_surf(:,:) * (1._wp + rsbc_tke1)**r2_3 , rn_emin ) 321 zd_lw(:,:,1) = en(:,:,1) 322 zd_up(:,:,1) = 0._wp 323 zdiag(:,:,1) = 1._wp 345 324 ! 346 325 ! at k=2, set de/dz=Fw 347 326 !cbr 348 z _elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b349 z _elem_a(:,:,2) = 0._wp350 zkar (:,:)= (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans*gdept_n(:,:,1)/zhsro(:,:)) ))351 zflxs(:,:) 352 & * ((zhsro(:,:)+gdept_n(:,:,1)) / zhsro(:,:))**(1.5_wp*ra_sf)353 354 en(:,:,2) = en(:,:,2) + zflxs(:,:) /e3w_n(:,:,2)327 zdiag(:,:,2) = zdiag(:,:,2) + zd_lw(:,:,2) ! Remove zd_lw from zdiag 328 zd_lw(:,:,2) = 0._wp 329 zkar (:,:) = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans*gdept_n(:,:,1)/zhsro(:,:)) )) 330 zflxs(:,:) = rsbc_tke2 * ustar2_surf(:,:)**1.5_wp * zkar(:,:) & 331 & * ( ( zhsro(:,:)+gdept_n(:,:,1) ) / zhsro(:,:) )**(1.5_wp*ra_sf) 332 !!gm why not : * ( 1._wp + gdept_n(:,:,1) / zhsro(:,:) )**(1.5_wp*ra_sf) 333 en(:,:,2) = en(:,:,2) + zflxs(:,:) / e3w_n(:,:,2) 355 334 ! 356 335 ! … … 377 356 ! Dirichlet condition applied at: 378 357 ! Bottom level (ibot) & Just above it (ibotm1) 379 z _elem_a(ji,jj,ibot) = 0._wp ; z_elem_a(ji,jj,ibotm1) = 0._wp380 z _elem_c(ji,jj,ibot) = 0._wp ; z_elem_c(ji,jj,ibotm1) = 0._wp381 z _elem_b(ji,jj,ibot) = 1._wp ; z_elem_b(ji,jj,ibotm1) = 1._wp382 en (ji,jj,ibot) = z_en ; en(ji,jj,ibotm1) = z_en383 END DO 384 END DO 385 !!gm new 386 !IF( ln_isfcav) THEN ! top boundary (ocean cavity)387 !DO jj = 2, jpjm1388 !DO ji = fs_2, fs_jpim1 ! vector opt.389 !itop = mikt(ji,jj) ! k top w-point390 !itopp1 = mikt(ji,jj) + 1 ! k+1 1st w-point below the top one391 !! ! mask at the ocean surface points392 !z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) )393 !!394 ! ! Dirichlet condition applied at: 395 ! ! top level (itop) & Just below it (itopp1)396 ! z_elem_a(ji,jj,itop) = 0._wp ; z_elem_a(ji,jj,ipotm1) = 0._wp 397 ! z_elem_c(ji,jj,itop) = 0._wp ; z_elem_c(ji,jj,itopp1) = 0._wp398 ! z_elem_b(ji,jj,itop) = 1._wp ; z_elem_b(ji,jj,itopp1) = 1._wp399 ! en (ji,jj,itop) = zen ; en (ji,jj,itopp1) = z_en 400 ! END DO 401 !END DO402 ! ENDIF 403 !!gm 358 zd_lw(ji,jj,ibot) = 0._wp ; zd_lw(ji,jj,ibotm1) = 0._wp 359 zd_up(ji,jj,ibot) = 0._wp ; zd_up(ji,jj,ibotm1) = 0._wp 360 zdiag(ji,jj,ibot) = 1._wp ; zdiag(ji,jj,ibotm1) = 1._wp 361 en (ji,jj,ibot) = z_en ; en (ji,jj,ibotm1) = z_en 362 END DO 363 END DO 364 ! 365 IF( ln_isfcav) THEN ! top boundary (ocean cavity) 366 DO jj = 2, jpjm1 367 DO ji = fs_2, fs_jpim1 ! vector opt. 368 itop = mikt(ji,jj) ! k top w-point 369 itopp1 = mikt(ji,jj) + 1 ! k+1 1st w-point below the top one 370 ! ! mask at the ocean surface points 371 z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) ) 372 ! 373 !!gm TO BE VERIFIED !!! 374 ! Dirichlet condition applied at: 375 ! top level (itop) & Just below it (itopp1) 376 zd_lw(ji,jj,itop) = 0._wp ; zd_lw(ji,jj,itopp1) = 0._wp 377 zd_up(ji,jj,itop) = 0._wp ; zd_up(ji,jj,itopp1) = 0._wp 378 zdiag(ji,jj,itop) = 1._wp ; zdiag(ji,jj,itopp1) = 1._wp 379 en (ji,jj,itop) = z_en ; en (ji,jj,itopp1) = z_en 380 END DO 381 END DO 382 ENDIF 404 383 ! 405 384 CASE ( 1 ) ! Neumman boundary condition … … 415 394 ! Bottom level (ibot) & Just above it (ibotm1) 416 395 ! Dirichlet ! Neumann 417 z_elem_a(ji,jj,ibot) = 0._wp ! ! Remove z_elem_c from z_elem_b 418 z_elem_b(ji,jj,ibot) = 1._wp ; z_elem_b(ji,jj,ibotm1) = z_elem_b(ji,jj,ibotm1) + z_elem_c(ji,jj,ibotm1) 419 z_elem_c(ji,jj,ibot) = 0._wp ; z_elem_c(ji,jj,ibotm1) = 0._wp 420 END DO 421 END DO 422 !!gm new 423 ! IF( ln_isfcav) THEN ! top boundary (ocean cavity) 424 ! DO jj = 2, jpjm1 425 ! DO ji = fs_2, fs_jpim1 ! vector opt. 426 ! itop = mikt(ji,jj) ! k top w-point 427 ! itopp1 = mikt(ji,jj) + 1 ! k+1 1st w-point below the top one 428 ! ! ! mask at the ocean surface points 429 ! z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) ) 430 ! ! 431 ! ! Bottom level Dirichlet condition: 432 ! ! Bottom level (ibot) & Just above it (ibotm1) 433 ! ! Dirichlet ! Neumann 434 ! z_elem_a(ji,jj,itop) = 0._wp ! ! Remove z_elem_c from z_elem_b 435 ! z_elem_b(ji,jj,itop) = 1._wp ; z_elem_b(ji,jj,itopp1) = z_elem_b(ji,jj,itopp1) + z_elem_c(ji,jj,itopp1) 436 ! z_elem_c(ji,jj,itop) = 0._wp ; z_elem_c(ji,jj,itopp1) = 0._wp 437 ! END DO 438 ! END DO 439 ! ENDIF 440 !!gm 396 zd_lw(ji,jj,ibot) = 0._wp ! ! Remove zd_up from zdiag 397 zdiag(ji,jj,ibot) = 1._wp ; zdiag(ji,jj,ibotm1) = zdiag(ji,jj,ibotm1) + zd_up(ji,jj,ibotm1) 398 zd_up(ji,jj,ibot) = 0._wp ; zd_up(ji,jj,ibotm1) = 0._wp 399 END DO 400 END DO 401 IF( ln_isfcav) THEN ! top boundary (ocean cavity) 402 DO jj = 2, jpjm1 403 DO ji = fs_2, fs_jpim1 ! vector opt. 404 itop = mikt(ji,jj) ! k top w-point 405 itopp1 = mikt(ji,jj) + 1 ! k+1 1st w-point below the top one 406 ! ! mask at the ocean surface points 407 z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) ) 408 ! 409 ! Bottom level Dirichlet condition: 410 ! Bottom level (ibot) & Just above it (ibotm1) 411 ! Dirichlet ! Neumann 412 zd_lw(ji,jj,itop) = 0._wp ! ! Remove zd_up from zdiag 413 zdiag(ji,jj,itop) = 1._wp ; zdiag(ji,jj,itopp1) = zdiag(ji,jj,itopp1) + zd_up(ji,jj,itopp1) 414 zd_up(ji,jj,itop) = 0._wp ; zd_up(ji,jj,itopp1) = 0._wp 415 END DO 416 END DO 417 ENDIF 441 418 ! 442 419 END SELECT … … 448 425 DO jj = 2, jpjm1 449 426 DO ji = fs_2, fs_jpim1 ! vector opt. 450 z _elem_b(ji,jj,jk) = z_elem_b(ji,jj,jk) - z_elem_a(ji,jj,jk) * z_elem_c(ji,jj,jk-1) / z_elem_b(ji,jj,jk-1)427 zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 451 428 END DO 452 429 END DO … … 455 432 DO jj = 2, jpjm1 456 433 DO ji = fs_2, fs_jpim1 ! vector opt. 457 z _elem_a(ji,jj,jk) = en(ji,jj,jk) - z_elem_a(ji,jj,jk) / z_elem_b(ji,jj,jk-1) * z_elem_a(ji,jj,jk-1)434 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) 458 435 END DO 459 436 END DO … … 462 439 DO jj = 2, jpjm1 463 440 DO ji = fs_2, fs_jpim1 ! vector opt. 464 en(ji,jj,jk) = ( z _elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) * en(ji,jj,jk+1) ) / z_elem_b(ji,jj,jk)441 en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 465 442 END DO 466 443 END DO … … 519 496 ! Resolution of a tridiagonal linear system by a "methode de chasse" 520 497 ! computation from level 2 to jpkm1 (e(1) already computed and e(jpk)=0 ). 521 ! z _elem_b : diagonal z_elem_c : upper diagonal z_elem_a: lower diagonal498 ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal 522 499 ! Warning : after this step, en : right hand side of the matrix 523 500 … … 543 520 diss = rpsi2 * zratio * zwall(ji,jj,jk) * eps(ji,jj,jk) 544 521 ! 545 zdir = 0.5_wp + SIGN( 0.5_wp, prod + buoy ) ! zdir =1(=0) if shear(ji,jj,jk)+buoy >0(<0)522 zdir = 0.5_wp + SIGN( 0.5_wp, prod + buoy ) ! zdir =1(=0) if shear(ji,jj,jk)+buoy >0(<0) 546 523 ! 547 524 zesh2 = zdir * ( prod + buoy ) + (1._wp - zdir ) * prod ! production term … … 551 528 zcof = rfact_psi * zwall_psi(ji,jj,jk) * tmask(ji,jj,jk) 552 529 ! ! lower diagonal 553 z _elem_a(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) )530 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) ) 554 531 ! ! upper diagonal 555 z _elem_c(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) )532 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) ) 556 533 ! ! diagonal 557 z _elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) + rdt * zdiss * wmask(ji,jj,jk)534 zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) + rdt * zdiss * wmask(ji,jj,jk) 558 535 ! ! right hand side in psi 559 536 psi(ji,jj,jk) = psi(ji,jj,jk) + rdt * zesh2 * wmask(ji,jj,jk) … … 562 539 END DO 563 540 ! 564 z _elem_b(:,:,jpk) = 1._wp541 zdiag(:,:,jpk) = 1._wp 565 542 566 543 ! Surface boundary condition on psi … … 574 551 zdep (:,:) = zhsro(:,:) * rl_sf ! Cosmetic 575 552 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 576 z _elem_a(:,:,1) = psi(:,:,1)577 z _elem_c(:,:,1) = 0._wp578 z _elem_b(:,:,1) = 1._wp553 zd_lw(:,:,1) = psi(:,:,1) 554 zd_up(:,:,1) = 0._wp 555 zdiag(:,:,1) = 1._wp 579 556 ! 580 557 ! One level below … … 582 559 zdep (:,:) = (zhsro(:,:) + gdepw_n(:,:,2)) * zkar(:,:) 583 560 psi (:,:,2) = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 584 z _elem_a(:,:,2) = 0._wp585 z _elem_c(:,:,2) = 0._wp586 z _elem_b(:,:,2) = 1._wp561 zd_lw(:,:,2) = 0._wp 562 zd_up(:,:,2) = 0._wp 563 zdiag(:,:,2) = 1._wp 587 564 ! 588 565 CASE ( 1 ) ! Neumann boundary condition on d(psi)/dz … … 591 568 zdep (:,:) = zhsro(:,:) * rl_sf 592 569 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 593 z _elem_a(:,:,1) = psi(:,:,1)594 z _elem_c(:,:,1) = 0._wp595 z _elem_b(:,:,1) = 1._wp570 zd_lw(:,:,1) = psi(:,:,1) 571 zd_up(:,:,1) = 0._wp 572 zdiag(:,:,1) = 1._wp 596 573 ! 597 574 ! Neumann condition at k=2 598 z _elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b599 z _elem_a(:,:,2) = 0._wp575 zdiag(:,:,2) = zdiag(:,:,2) + zd_lw(:,:,2) ! Remove zd_lw from zdiag 576 zd_lw(:,:,2) = 0._wp 600 577 ! 601 578 ! Set psi vertical flux at the surface: … … 613 590 ! -------------------------------- 614 591 ! 615 SELECT CASE ( nn_bc_bot ) 592 !!gm should be done for ISF (top boundary cond.) 593 !!gm so, totally new staff needed ===>>> think about that ! 594 ! 595 SELECT CASE ( nn_bc_bot ) ! bottom boundary 616 596 ! 617 597 CASE ( 0 ) ! Dirichlet 618 ! ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = vkarmn * r n_bfrz0598 ! ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = vkarmn * r_z0_bot 619 599 ! ! Balance between the production and the dissipation terms 620 600 DO jj = 2, jpjm1 … … 622 602 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 623 603 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 624 zdep(ji,jj) = vkarmn * r n_bfrz0604 zdep(ji,jj) = vkarmn * r_z0_bot 625 605 psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn 626 z _elem_a(ji,jj,ibot) = 0._wp627 z _elem_c(ji,jj,ibot) = 0._wp628 z _elem_b(ji,jj,ibot) = 1._wp606 zd_lw(ji,jj,ibot) = 0._wp 607 zd_up(ji,jj,ibot) = 0._wp 608 zdiag(ji,jj,ibot) = 1._wp 629 609 ! 630 610 ! Just above last level, Dirichlet condition again (GOTM like) 631 zdep(ji,jj) = vkarmn * ( r n_bfrz0+ e3t_n(ji,jj,ibotm1) )611 zdep(ji,jj) = vkarmn * ( r_z0_bot + e3t_n(ji,jj,ibotm1) ) 632 612 psi (ji,jj,ibotm1) = rc0**rpp * en(ji,jj,ibot )**rmm * zdep(ji,jj)**rnn 633 z _elem_a(ji,jj,ibotm1) = 0._wp634 z _elem_c(ji,jj,ibotm1) = 0._wp635 z _elem_b(ji,jj,ibotm1) = 1._wp613 zd_lw(ji,jj,ibotm1) = 0._wp 614 zd_up(ji,jj,ibotm1) = 0._wp 615 zdiag(ji,jj,ibotm1) = 1._wp 636 616 END DO 637 617 END DO … … 645 625 ! 646 626 ! Bottom level Dirichlet condition: 647 zdep(ji,jj) = vkarmn * r n_bfrz0627 zdep(ji,jj) = vkarmn * r_z0_bot 648 628 psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn 649 629 ! 650 z _elem_a(ji,jj,ibot) = 0._wp651 z _elem_c(ji,jj,ibot) = 0._wp652 z _elem_b(ji,jj,ibot) = 1._wp630 zd_lw(ji,jj,ibot) = 0._wp 631 zd_up(ji,jj,ibot) = 0._wp 632 zdiag(ji,jj,ibot) = 1._wp 653 633 ! 654 634 ! Just above last level: Neumann condition with flux injection 655 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_b656 z _elem_c(ji,jj,ibotm1) = 0.635 zdiag(ji,jj,ibotm1) = zdiag(ji,jj,ibotm1) + zd_up(ji,jj,ibotm1) ! Remove zd_up from zdiag 636 zd_up(ji,jj,ibotm1) = 0. 657 637 ! 658 638 ! Set psi vertical flux at the bottom: 659 zdep(ji,jj) = r n_bfrz0+ 0.5_wp*e3t_n(ji,jj,ibotm1)639 zdep(ji,jj) = r_z0_bot + 0.5_wp*e3t_n(ji,jj,ibotm1) 660 640 zflxb = rsbc_psi2 * ( p_avm(ji,jj,ibot) + p_avm(ji,jj,ibotm1) ) & 661 641 & * (0.5_wp*(en(ji,jj,ibot)+en(ji,jj,ibotm1)))**rmm * zdep(ji,jj)**(rnn-1._wp) … … 672 652 DO jj = 2, jpjm1 673 653 DO ji = fs_2, fs_jpim1 ! vector opt. 674 z _elem_b(ji,jj,jk) = z_elem_b(ji,jj,jk) - z_elem_a(ji,jj,jk) * z_elem_c(ji,jj,jk-1) / z_elem_b(ji,jj,jk-1)654 zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 675 655 END DO 676 656 END DO … … 679 659 DO jj = 2, jpjm1 680 660 DO ji = fs_2, fs_jpim1 ! vector opt. 681 z _elem_a(ji,jj,jk) = psi(ji,jj,jk) - z_elem_a(ji,jj,jk) / z_elem_b(ji,jj,jk-1) * z_elem_a(ji,jj,jk-1)661 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) 682 662 END DO 683 663 END DO … … 686 666 DO jj = 2, jpjm1 687 667 DO ji = fs_2, fs_jpim1 ! vector opt. 688 psi(ji,jj,jk) = ( z _elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) * psi(ji,jj,jk+1) ) / z_elem_b(ji,jj,jk)668 psi(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * psi(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 689 669 END DO 690 670 END DO … … 814 794 zstm(:,:,1) = zstm(:,:,2) 815 795 816 !!gm should be done for ISF (top boundary cond.)817 796 DO jj = 2, jpjm1 818 797 DO ji = fs_2, fs_jpim1 ! vector opt. … … 820 799 END DO 821 800 END DO 801 !!gm should be done for ISF (top boundary cond.) 802 !!gm so, totally new staff needed!!gm 822 803 823 804 ! Compute diffusivities/viscosities … … 904 885 WRITE(numout,*) ' Type of closure nn_clos = ', nn_clos 905 886 WRITE(numout,*) ' Surface roughness (m) rn_hsro = ', rn_hsro 906 !!gm old 907 WRITE(numout,*) ' top roughness (m) (nambfr namelist) rn_tfrz0 = ', rn_tfrz0 908 WRITE(numout,*) ' Bottom roughness (m) (nambfr namelist) rn_bfrz0 = ', rn_bfrz0 909 !!gm new 910 ! WRITE(numout,*) ' Namelist namdrg_top/_bot used values:' 911 ! WRITE(numout,*) ' top ocean cavity roughness (m) rn_z0(_top) = ', r_z0_top 912 ! WRITE(numout,*) ' Bottom seafloor roughness (m) rn_z0(_bot) = ', r_z0_bot 913 !!gm 887 WRITE(numout,*) 888 WRITE(numout,*) ' Namelist namdrg_top/_bot: used values:' 889 WRITE(numout,*) ' top ocean cavity roughness (m) rn_z0(_top) = ', r_z0_top 890 WRITE(numout,*) ' Bottom seafloor roughness (m) rn_z0(_bot) = ', r_z0_bot 914 891 WRITE(numout,*) 915 892 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
