- Timestamp:
- 2020-10-01T13:33:30+02:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfgls.F90
r13295 r13553 19 19 USE dom_oce ! ocean space and time domain 20 20 USE domvvl ! ocean space and time domain : variable volume layer 21 USE zdfdrg , ONLY : ln_drg_OFF ! top/bottom free-slip flag 21 22 USE zdfdrg , ONLY : r_z0_top , r_z0_bot ! top/bottom roughness 22 23 USE zdfdrg , ONLY : rCdU_top , rCdU_bot ! top/bottom friction … … 53 54 INTEGER :: nn_bc_bot ! bottom boundary condition (=0/1) 54 55 INTEGER :: nn_z0_met ! Method for surface roughness computation 56 INTEGER :: nn_z0_ice ! Roughness accounting for sea ice 55 57 INTEGER :: nn_stab_func ! stability functions G88, KC or Canuto (=0/1/2) 56 58 INTEGER :: nn_clos ! closure 0/1/2/3 MY82/k-eps/k-w/gen … … 61 63 REAL(wp) :: rn_crban ! Craig and Banner constant for surface breaking waves mixing 62 64 REAL(wp) :: rn_hsro ! Minimum surface roughness 65 REAL(wp) :: rn_hsri ! Ice ocean roughness 63 66 REAL(wp) :: rn_frac_hs ! Fraction of wave height as surface roughness (if nn_z0_met > 1) 64 67 … … 152 155 REAL(wp), DIMENSION(jpi,jpj) :: zflxs ! Turbulence fluxed induced by internal waves 153 156 REAL(wp), DIMENSION(jpi,jpj) :: zhsro ! Surface roughness (surface waves) 157 REAL(wp), DIMENSION(jpi,jpj) :: zice_fra ! Tapering of wave breaking under sea ice 154 158 REAL(wp), DIMENSION(jpi,jpj,jpk) :: eb ! tke at time before 155 159 REAL(wp), DIMENSION(jpi,jpj,jpk) :: hmxl_b ! mixing length at time before … … 167 171 ustar2_bot (:,:) = 0._wp 168 172 173 SELECT CASE ( nn_z0_ice ) 174 CASE( 0 ) ; zice_fra(:,:) = 0._wp 175 CASE( 1 ) ; zice_fra(:,:) = TANH( fr_i(:,:) * 10._wp ) 176 CASE( 2 ) ; zice_fra(:,:) = fr_i(:,:) 177 CASE( 3 ) ; zice_fra(:,:) = MIN( 4._wp * fr_i(:,:) , 1._wp ) 178 END SELECT 179 169 180 ! Compute surface, top and bottom friction at T-points 170 DO_2D( 0, 0, 0, 0 ) 171 ! 172 ! surface friction 173 ustar2_surf(ji,jj) = r1_rho0 * taum(ji,jj) * tmask(ji,jj,1) 174 ! 175 !!gm Rq we may add here r_ke0(_top/_bot) ? ==>> think about that... 176 ! bottom friction (explicit before friction) 177 zmsku = ( 2._wp - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 178 zmskv = ( 2._wp - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) ! (CAUTION: CdU<0) 179 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 & 180 & + ( zmskv*( vv(ji,jj,mbkt(ji,jj),Kbb)+vv(ji,jj-1,mbkt(ji,jj),Kbb) ) )**2 ) 181 DO_2D( 0, 0, 0, 0 ) !== surface ocean friction 182 ustar2_surf(ji,jj) = r1_rho0 * taum(ji,jj) * tmask(ji,jj,1) ! surface friction 181 183 END_2D 182 IF( ln_isfcav ) THEN !top friction 183 DO_2D( 0, 0, 0, 0 ) 184 zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) ) 185 zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) ) ! (CAUTION: CdU<0) 186 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 & 187 & + ( zmskv*( vv(ji,jj,mikt(ji,jj),Kbb)+vv(ji,jj-1,mikt(ji,jj),Kbb) ) )**2 ) 184 ! 185 !!gm Rq we may add here r_ke0(_top/_bot) ? ==>> think about that... 186 ! 187 IF( .NOT.ln_drg_OFF ) THEN !== top/bottom friction (explicit before friction) 188 DO_2D( 0, 0, 0, 0 ) ! bottom friction (explicit before friction) 189 zmsku = ( 2._wp - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 190 zmskv = ( 2._wp - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) ! (CAUTION: CdU<0) 191 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 & 192 & + ( zmskv*( vv(ji,jj,mbkt(ji,jj),Kbb)+vv(ji,jj-1,mbkt(ji,jj),Kbb) ) )**2 ) 188 193 END_2D 194 IF( ln_isfcav ) THEN 195 DO_2D( 0, 0, 0, 0 ) ! top friction 196 zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) ) 197 zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) ) ! (CAUTION: CdU<0) 198 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 & 199 & + ( zmskv*( vv(ji,jj,mikt(ji,jj),Kbb)+vv(ji,jj-1,mikt(ji,jj),Kbb) ) )**2 ) 200 END_2D 201 ENDIF 189 202 ENDIF 190 203 … … 204 217 END SELECT 205 218 ! 206 DO_3D( 1, 0, 1, 0, 2, jpkm1 ) 219 ! adapt roughness where there is sea ice 220 zhsro(:,:) = ( (1._wp-zice_fra(:,:)) * zhsro(:,:) + zice_fra(:,:) * rn_hsri )*tmask(:,:,1) + (1._wp - tmask(:,:,1))*rn_hsro 221 ! 222 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) !== Compute dissipation rate ==! 207 223 eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / hmxl_n(ji,jj,jk) 208 224 END_3D … … 288 304 CASE ( 0 ) ! Dirichlet boundary condition (set e at k=1 & 2) 289 305 ! First level 290 en (:,:,1) = MAX( rn_emin , rc02r * ustar2_surf(:,:) * (1._wp + rsbc_tke1)**r2_3 )306 en (:,:,1) = MAX( rn_emin , rc02r * ustar2_surf(:,:) * (1._wp + (1._wp-zice_fra(:,:))*rsbc_tke1)**r2_3 ) 291 307 zd_lw(:,:,1) = en(:,:,1) 292 308 zd_up(:,:,1) = 0._wp … … 294 310 ! 295 311 ! One level below 296 en (:,:,2) = MAX( rc02r * ustar2_surf(:,:) * ( 1._wp + rsbc_tke1 * ((zhsro(:,:)+gdepw(:,:,2,Kmm))&297 & / zhsro(:,:) )**(1.5_wp*ra_sf) )**(2._wp/3._wp) 312 en (:,:,2) = MAX( rc02r * ustar2_surf(:,:) * ( 1._wp + (1._wp-zice_fra(:,:))*rsbc_tke1 * ((zhsro(:,:)+gdepw(:,:,2,Kmm)) & 313 & / zhsro(:,:) )**(1.5_wp*ra_sf) )**(2._wp/3._wp) , rn_emin ) 298 314 zd_lw(:,:,2) = 0._wp 299 315 zd_up(:,:,2) = 0._wp … … 304 320 ! 305 321 ! Dirichlet conditions at k=1 306 en (:,:,1) = MAX( rc02r * ustar2_surf(:,:) * (1._wp + rsbc_tke1)**r2_3 , rn_emin )322 en (:,:,1) = MAX( rc02r * ustar2_surf(:,:) * (1._wp + (1._wp-zice_fra(:,:))*rsbc_tke1)**r2_3 , rn_emin ) 307 323 zd_lw(:,:,1) = en(:,:,1) 308 324 zd_up(:,:,1) = 0._wp … … 311 327 ! at k=2, set de/dz=Fw 312 328 !cbr 313 zdiag(:,:,2) = zdiag(:,:,2) + zd_lw(:,:,2) ! Remove zd_lw from zdiag 314 zd_lw(:,:,2) = 0._wp 329 DO_2D( 0, 0, 0, 0 ) ! zdiag zd_lw not defined/used on the halo 330 zdiag(ji,jj,2) = zdiag(ji,jj,2) + zd_lw(ji,jj,2) ! Remove zd_lw from zdiag 331 zd_lw(ji,jj,2) = 0._wp 332 END_2D 315 333 zkar (:,:) = (rl_sf + (vkarmn-rl_sf)*(1.-EXP(-rtrans*gdept(:,:,1,Kmm)/zhsro(:,:)) )) 316 zflxs(:,:) = rsbc_tke2 * ustar2_surf(:,:)**1.5_wp * zkar(:,:) &334 zflxs(:,:) = rsbc_tke2 * (1._wp-zice_fra(:,:)) * ustar2_surf(:,:)**1.5_wp * zkar(:,:) & 317 335 & * ( ( zhsro(:,:)+gdept(:,:,1,Kmm) ) / zhsro(:,:) )**(1.5_wp*ra_sf) 318 336 !!gm why not : * ( 1._wp + gdept(:,:,1,Kmm) / zhsro(:,:) )**(1.5_wp*ra_sf) … … 400 418 ! ---------------------------------------------------------- 401 419 ! 402 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 420 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 403 421 zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 404 422 END_3D 405 DO_3D( 0, 0, 0, 0, 2, jpk )423 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 406 424 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) 407 425 END_3D 408 DO_3DS( 0, 0, 0, 0, jpk -1, 2, -1 )426 DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 ) ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 409 427 en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 410 428 END_3D … … 521 539 ! 522 540 ! Neumann condition at k=2 523 zdiag(:,:,2) = zdiag(:,:,2) + zd_lw(:,:,2) ! Remove zd_lw from zdiag 524 zd_lw(:,:,2) = 0._wp 541 DO_2D( 0, 0, 0, 0 ) ! zdiag zd_lw not defined/used on the halo 542 zdiag(ji,jj,2) = zdiag(ji,jj,2) + zd_lw(ji,jj,2) ! Remove zd_lw from zdiag 543 zd_lw(ji,jj,2) = 0._wp 544 END_2D 525 545 ! 526 546 ! Set psi vertical flux at the surface: 527 547 zkar (:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdept(:,:,1,Kmm)/zhsro(:,:) )) ! Lengh scale slope 528 548 zdep (:,:) = ((zhsro(:,:) + gdept(:,:,1,Kmm)) / zhsro(:,:))**(rmm*ra_sf) 529 zflxs(:,:) = (rnn + rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:))*(1._wp + rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 549 zflxs(:,:) = (rnn + (1._wp-zice_fra(:,:))*rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:)) & 550 & *(1._wp + (1._wp-zice_fra(:,:))*rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 530 551 zdep (:,:) = rsbc_psi1 * (zwall_psi(:,:,1)*p_avm(:,:,1)+zwall_psi(:,:,2)*p_avm(:,:,2)) * & 531 552 & ustar2_surf(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + gdept(:,:,1,Kmm))**(rnn-1.) … … 593 614 ! ---------------- 594 615 ! 595 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 616 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 596 617 zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 597 618 END_3D 598 DO_3D( 0, 0, 0, 0, 2, jpk )619 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 599 620 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) 600 621 END_3D 601 DO_3DS( 0, 0, 0, 0, jpk -1, 2, -1 )622 DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 ) ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 602 623 psi(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * psi(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 603 624 END_3D … … 635 656 ! Limit dissipation rate under stable stratification 636 657 ! -------------------------------------------------- 637 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 658 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! Note that this set boundary conditions on hmxl_n at the same time 638 659 ! limitation 639 660 eps (ji,jj,jk) = MAX( eps(ji,jj,jk), rn_epsmin ) … … 700 721 ! default value, in case jpk > mbkt(ji,jj)+1. Not needed but avoid a bug when looking for undefined values (-fpe0) 701 722 zstm(:,:,jpk) = 0. 702 DO_2D( 0, 0, 0, 0 ) 723 DO_2D( 0, 0, 0, 0 ) ! update bottom with good values 703 724 zstm(ji,jj,mbkt(ji,jj)+1) = zstm(ji,jj,mbkt(ji,jj)) 704 725 END_2D … … 750 771 REAL(wp):: zcr ! local scalar 751 772 !! 752 NAMELIST/namzdf_gls/rn_emin, rn_epsmin, ln_length_lim, &753 & rn_clim_galp, ln_sigpsi, rn_hsro, 754 & rn_crban, rn_charn, rn_frac_hs, &755 & nn_bc_surf, nn_bc_bot, nn_z0_met, 773 NAMELIST/namzdf_gls/rn_emin, rn_epsmin, ln_length_lim, & 774 & rn_clim_galp, ln_sigpsi, rn_hsro, rn_hsri, & 775 & rn_crban, rn_charn, rn_frac_hs, & 776 & nn_bc_surf, nn_bc_bot, nn_z0_met, nn_z0_ice, & 756 777 & nn_stab_func, nn_clos 757 778 !!---------------------------------------------------------- … … 779 800 WRITE(numout,*) ' Charnock coefficient rn_charn = ', rn_charn 780 801 WRITE(numout,*) ' Surface roughness formula nn_z0_met = ', nn_z0_met 802 WRITE(numout,*) ' surface wave breaking under ice nn_z0_ice = ', nn_z0_ice 803 SELECT CASE( nn_z0_ice ) 804 CASE( 0 ) ; WRITE(numout,*) ' ==>>> no impact of ice cover on surface wave breaking' 805 CASE( 1 ) ; WRITE(numout,*) ' ==>>> roughness uses rn_hsri and is weigthed by 1-TANH( fr_i(:,:) * 10 )' 806 CASE( 2 ) ; WRITE(numout,*) ' ==>>> roughness uses rn_hsri and is weighted by 1-fr_i(:,:)' 807 CASE( 3 ) ; WRITE(numout,*) ' ==>>> roughness uses rn_hsri and is weighted by 1-MIN( 1, 4 * fr_i(:,:) )' 808 CASE DEFAULT 809 CALL ctl_stop( 'zdf_gls_init: wrong value for nn_z0_ice, should be 0,1,2, or 3') 810 END SELECT 781 811 WRITE(numout,*) ' Wave height frac. (used if nn_z0_met=2) rn_frac_hs = ', rn_frac_hs 782 812 WRITE(numout,*) ' Stability functions nn_stab_func = ', nn_stab_func 783 813 WRITE(numout,*) ' Type of closure nn_clos = ', nn_clos 784 814 WRITE(numout,*) ' Surface roughness (m) rn_hsro = ', rn_hsro 815 WRITE(numout,*) ' Ice-ocean roughness (used if nn_z0_ice/=0) rn_hsri = ', rn_hsri 785 816 WRITE(numout,*) 786 817 WRITE(numout,*) ' Namelist namdrg_top/_bot: used values:'
Note: See TracChangeset
for help on using the changeset viewer.
