MODULE zdfgls !!====================================================================== !! *** MODULE zdfgls *** !! Ocean physics: vertical mixing coefficient computed from the gls !! turbulent closure parameterization !!====================================================================== !! History : 3.0 ! 2009-09 (G. Reffray) : Original code !!---------------------------------------------------------------------- #if defined key_zdfgls || defined key_esopa !!---------------------------------------------------------------------- !! 'key_zdfgls' Generic Length Scale vertical physics !!---------------------------------------------------------------------- !!---------------------------------------------------------------------- !! zdf_gls : update momentum and tracer Kz from a gls scheme !! zdf_gls_init : initialization, namelist read, and parameters control !! gls_rst : read/write gls restart in ocean restart file !!---------------------------------------------------------------------- USE oce ! ocean dynamics and active tracers USE dom_oce ! ocean space and time domain USE domvvl ! ocean space and time domain : variable volume layer USE zdf_oce ! ocean vertical physics USE sbc_oce ! surface boundary condition: ocean USE phycst ! physical constants USE zdfmxl ! mixed layer USE restart ! only for lrst_oce USE lbclnk ! ocean lateral boundary conditions (or mpp link) USE prtctl ! Print control USE in_out_manager ! I/O manager USE iom ! I/O manager library USE zdfbfr, ONLY : rn_hbro, wbotu, wbotv ! bottom roughness and bottom stresses IMPLICIT NONE PRIVATE PUBLIC zdf_gls ! routine called in step module PUBLIC gls_rst ! routine called in step module LOGICAL , PUBLIC, PARAMETER :: lk_zdfgls = .TRUE. !: TKE vertical mixing flag REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: en !: now turbulent kinetic energy REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: mxln !: now mixing length REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: zwall !: wall function REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: ustars2 !: Squared surface velocity scale at T-points REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: ustarb2 !: Squared bottom velocity scale at T-points ! !!! ** Namelist namzdf_gls ** LOGICAL :: ln_crban = .FALSE. ! =T use Craig and Banner scheme LOGICAL :: ln_length_lim = .FALSE. ! use limit on the dissipation rate understable stratification (Galperin et al., 1988) LOGICAL :: ln_sigpsi = .FALSE. ! Activate Burchard (2003) modification for k-eps closure AND wave breaking mixing REAL(wp) :: rn_epsmin = 1.e-12_wp ! minimum value of dissipation (m2/s3) REAL(wp) :: rn_emin = 1.e-6_wp ! minimum value of TKE (m2/s2) INTEGER :: nn_tkebc_surf = 0 ! TKE surface boundary condition (=0/1) INTEGER :: nn_tkebc_bot = 0 ! TKE bottom boundary condition (=0/1) INTEGER :: nn_psibc_surf = 0 ! PSI surface boundary condition (=0/1) INTEGER :: nn_psibc_bot = 0 ! PSI bottom boundary condition (=0/1) INTEGER :: nn_stab_func = 0 ! stability functions G88, KC or Canuto (=0/1/2) INTEGER :: nn_clos = 0 ! closure 0/1/2/3 MY82/k-eps/k-w/gen REAL(wp) :: clim_galp = 0.53_wp ! Holt 2008 value for k-eps: 0.267 REAL(wp) :: hsro = 0.003_wp ! Minimum surface roughness REAL(wp) :: rn_charn = 2.e+5_wp ! Charnock constant for surface breaking waves mixing : 1400. (standard) or 2.e5 (Stacey value) REAL(wp) :: rn_crban = 100._wp ! Craig and Banner constant for surface breaking waves mixing REAL(wp) :: rn_cm_sf = 0.73_wp ! Shear free turbulence parameters REAL(wp) :: rn_a_sf = -2.0_wp ! Must be negative -2 < rn_a_sf < -1 REAL(wp) :: rn_l_sf = 0.2_wp ! 0 ua ! use ua as workspace USE oce, z_elem_b => va ! use va as workspace USE oce, z_elem_c => ta ! use ta as workspace USE oce, psi => sa ! use sa as workspace ! INTEGER, INTENT(in) :: kt ! ocean time step INTEGER :: ji, jj, jk, ibot, ibotm1, dir ! dummy loop arguments ! REAL(wp) :: zesh2, zsigpsi ! temporary scalars REAL(wp) :: ztx2, zty2, zup, zdown, zcof ! - REAL(wp) :: zratio, zrn2, zflxb, sh ! - - REAL(wp) :: prod, buoy, diss, zdiss, sm ! - - REAL(wp) :: gh, gm, shr, dif, zsqen, zav ! - - REAL(wp), DIMENSION(jpi,jpj) :: zdep ! REAL(wp), DIMENSION(jpi,jpj) :: zflxs ! Turbulence fluxed induced by internal waves REAL(wp), DIMENSION(jpi,jpj) :: zhsro ! Surface roughness (surface waves) REAL(wp), DIMENSION(jpi,jpj,jpk) :: eb ! tke at time before REAL(wp), DIMENSION(jpi,jpj,jpk) :: mxlb ! mixing length at time before REAL(wp), DIMENSION(jpi,jpj,jpk) :: shear ! vertical shear REAL(wp), DIMENSION(jpi,jpj,jpk) :: eps ! dissipation rate REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwall_psi ! Wall function for psi schmidt number in the wb case ! (if ln_sigpsi.AND.ln_crban) !!-------------------------------------------------------------------- IF( kt == nit000 ) CALL zdf_gls_init ! Initialization (first time-step only) !!-------------------------------------------------------------------- ! Preliminary computing ustars2 = 0. ; ustarb2 = 0. ; psi = 0. ; zwall_psi = 0. ! Compute wind stress at T-points !CDIR NOVERRCHK DO jj = 2, jpjm1 !CDIR NOVERRCHK DO ji = fs_2, fs_jpim1 ! vector opt. ! ! wind stress ! squared surface velocity scale ztx2 = 2. * (utau(ji-1,jj )*umask(ji-1,jj,1) + utau(ji,jj)*umask(ji,jj,1)) / & & MAX(1., umask(ji-1,jj,1) + umask(ji,jj,1)) zty2 = 2. * (vtau(ji ,jj-1)*vmask(ji,jj-1,1) + vtau(ji,jj)*vmask(ji,jj,1)) / & & MAX(1., vmask(ji,jj-1,1) + vmask(ji,jj,1)) ustars2(ji,jj) = rn_sbc_tke2 * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) ! ! bottom friction ztx2 = 2. * (wbotu(ji-1,jj )*umask(ji-1,jj,1) + wbotu(ji,jj)*umask(ji,jj,1)) / & & MAX(1., umask(ji-1,jj,1) + umask(ji,jj,1)) zty2 = 2. * (wbotv(ji ,jj-1)*vmask(ji,jj-1,1) + wbotv(ji,jj)*vmask(ji,jj,1)) / & & MAX(1., vmask(ji,jj-1,1) + vmask(ji,jj,1)) ustarb2(ji,jj) = 0.5 * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) ENDDO ENDDO ! In case of breaking surface waves mixing, ! Compute surface roughness length according to Charnock formula: IF (ln_crban) THEN zhsro(:,:) = MAX(rn_sbc_zs * ustars2(:,:), hsro) ELSE zhsro(:,:) = hsro ENDIF ! Compute shear and dissipation rate DO jk = 2, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) ) & & * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) ) & & / ( fse3uw_n(ji,jj,jk) & & * fse3uw_b(ji,jj,jk) ) avmv(ji,jj,jk) = avmv(ji,jj,jk) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) ) & & * ( vb(ji,jj,jk-1) - vb(ji,jj,jk) ) & & / ( fse3vw_n(ji,jj,jk) & & * fse3vw_b(ji,jj,jk) ) eps(ji,jj,jk) = rn_c03 * en(ji,jj,jk) * SQRT(en(ji,jj,jk)) / mxln(ji,jj,jk) ENDDO ENDDO ENDDO ! ! Lateral boundary conditions (avmu,avmv) (sign unchanged) CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. ) ! Save tke at before time step eb (:,:,:) = en (:,:,:) mxlb(:,:,:) = mxln(:,:,:) IF ( nn_clos .EQ. 0 ) THEN ! Mellor-Yamada DO jk = 2, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zup = mxln(ji,jj,jk) * fsdepw(ji,jj,mbathy(ji,jj)) zdown = vkarmn * fsdepw(ji,jj,jk) * (-fsdepw(ji,jj,jk)+fsdepw(ji,jj,mbathy(ji,jj))) zwall (ji,jj,jk) = ( 1. + rn_e2 * ( zup / MAX( zdown, rsmall ) )**2. ) * tmask(ji,jj,jk) ENDDO ENDDO ENDDO ENDIF !!---------------------------------!! !! Equation to prognostic k !! !!---------------------------------!! ! ! Now Turbulent kinetic energy (output in en) ! ------------------------------- ! Resolution of a tridiagonal linear system by a "methode de chasse" ! computation from level 2 to jpkm1 (e(1) computed after and e(jpk)=0 ). ! The surface boundary condition are set after ! The bottom boundary condition are also set after. In standard e(bottom)=0. ! z_elem_b : diagonal z_elem_c : upper diagonal z_elem_a : lower diagonal ! Warning : after this step, en : right hand side of the matrix DO jk = 2, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ! ! shear prod. at w-point weightened by mask shear(ji,jj,jk) = ( avmu(ji-1,jj,jk) + avmu(ji,jj,jk) ) / MAX( 1.e0 , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) & & + ( avmv(ji,jj-1,jk) + avmv(ji,jj,jk) ) / MAX( 1.e0 , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) ) ! ! stratif. destruction buoy = - avt(ji,jj,jk) * rn2(ji,jj,jk) ! ! shear prod. - stratif. destruction diss = eps(ji,jj,jk) ! dir = 0.5 + sign(0.5,shear(ji,jj,jk)+buoy) ! dir =1(=0) if shear(ji,jj,jk)+buoy >0(<0) ! zesh2 = dir*(shear(ji,jj,jk)+buoy)+(1-dir)*shear(ji,jj,jk) ! production term zdiss = dir*(diss/en(ji,jj,jk)) +(1-dir)*(diss-buoy)/en(ji,jj,jk) ! dissipation term ! ! Compute a wall function from 1. to rn_sc_psi*zwall/rn_sc_psi0 ! Note that as long that Dirichlet boundary conditions are NOT set at the first and last levels (GOTM style) ! there is no need to set a boundary condition for zwall_psi at the top and bottom boundaries. ! Otherwise, this should be rn_sc_psi/rn_sc_psi0 IF (ln_sigpsi) THEN zsigpsi = MIN(1., zesh2/eps(ji,jj,jk)) ! 0. <= zsigpsi <= 1. zwall_psi(ji,jj,jk) = rn_sc_psi / (zsigpsi * rn_sc_psi + (1.-zsigpsi) * rn_sc_psi0 / MAX(zwall(ji,jj,jk),1.)) ELSE zwall_psi(ji,jj,jk) = 1. ENDIF ! ! building the matrix zcof = zfact_tke * tmask(ji,jj,jk) ! ! lower diagonal z_elem_a(ji,jj,jk) = zcof * ( avm (ji,jj,jk ) + avm (ji,jj,jk-1) ) & & / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk ) ) ! ! upper diagonal z_elem_c(ji,jj,jk) = zcof * ( avm (ji,jj,jk+1) + avm (ji,jj,jk ) ) & & / ( fse3t(ji,jj,jk ) * fse3w(ji,jj,jk) ) ! ! diagonal z_elem_b(ji,jj,jk) = 1. - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) & & + rdt * zdiss * tmask(ji,jj,jk) ! ! right hand side in en en(ji,jj,jk) = en(ji,jj,jk) + rdt * zesh2 * tmask(ji,jj,jk) END DO END DO END DO ! z_elem_b(:,:,jpk) = 1. ! ! Set surface condition on zwall_psi (1 at the bottom) IF (ln_sigpsi) THEN DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zwall_psi(ji,jj,1) = rn_sc_psi/rn_sc_psi0 END DO END DO ENDIF ! Surface boundary condition on tke ! --------------------------------- ! SELECT CASE ( nn_tkebc_surf ) ! CASE ( 0 ) ! Dirichlet case ! IF (ln_crban) THEN ! Wave induced mixing case ! ! en(1) = q2(1) = 0.5 * (15.8 * Ccb)^(2/3) * u*^2 ! ! balance between the production and the dissipation terms including the wave effect ! en(:,:,1) = MAX( rn_sbc_tke1 * ustars2(:,:), rn_emin ) z_elem_a(:,:,1) = en(:,:,1) z_elem_c(:,:,1) = 0. z_elem_b(:,:,1) = 1. ! ! one level below en(:,:,2) = MAX( rn_sbc_tke1 * ustars2(:,:) * ( (zhsro(:,:)+fsdepw(:,:,2))/zhsro(:,:) )**rn_a_sf, rn_emin ) z_elem_a(:,:,2) = 0. z_elem_c(:,:,2) = 0. z_elem_b(:,:,2) = 1. ! ELSE ! No wave induced mixing case ! ! en(1) = u*^2/C0^2 & l(1) = K*zs ! ! balance between the production and the dissipation terms ! en(:,:,1) = MAX( rn_c02r * ustars2(:,:), rn_emin ) z_elem_a(:,:,1) = en(:,:,1) z_elem_c(:,:,1) = 0. z_elem_b(:,:,1) = 1. ! ! one level below en(:,:,2) = MAX( rn_c02r * ustars2(:,:), rn_emin ) z_elem_a(:,:,2) = 0. z_elem_c(:,:,2) = 0. z_elem_b(:,:,2) = 1. ! ENDIF ! CASE ( 1 ) ! Neumann boundary condition on d(e)/dz ! IF (ln_crban) THEN ! Shear free case: d(e)/dz= Fw ! ! Dirichlet conditions at k=1 (Cosmetic) en(:,:,1) = MAX( rn_sbc_tke1 * ustars2(:,:), rn_emin ) z_elem_a(:,:,1) = en(:,:,1) z_elem_c(:,:,1) = 0. z_elem_b(:,:,1) = 1. ! at k=2, set de/dz=Fw z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b z_elem_a(:,:,2) = 0. zflxs(:,:) = rn_sbc_tke3 * ustars2(:,:)**1.5 * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5*rn_a_sf) en(:,:,2) = en(:,:,2) + zflxs(:,:)/fse3w(:,:,2) ! ELSE ! No wave induced mixing case: d(e)/dz=0. ! ! Dirichlet conditions at k=1 (Cosmetic) en(:,:,1) = MAX( rn_c02r * ustars2(:,:), rn_emin ) z_elem_a(:,:,1) = en(:,:,1) z_elem_c(:,:,1) = 0. z_elem_b(:,:,1) = 1. ! at k=2 set de/dz=0.: z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b z_elem_a(:,:,2) = 0. ! ENDIF ! END SELECT ! Bottom boundary condition on tke ! -------------------------------- ! SELECT CASE ( nn_tkebc_bot ) ! CASE ( 0 ) ! Dirichlet ! ! en(ibot) = u*^2 / Co2 and mxln(ibot) = rn_lmin ! ! Balance between the production and the dissipation terms ! !CDIR NOVERRCHK DO jj = 2, jpjm1 !CDIR NOVERRCHK DO ji = fs_2, fs_jpim1 ! vector opt. ibot = mbathy(ji,jj) ibotm1 = ibot-1 ! ! Bottom level Dirichlet condition: z_elem_a(ji,jj,ibot ) = 0. z_elem_c(ji,jj,ibot ) = 0. z_elem_b(ji,jj,ibot ) = 1. en(ji,jj,ibot ) = MAX( rn_c02r * ustarb2(ji,jj), rn_emin ) ! ! Just above last level, Dirichlet condition again z_elem_a(ji,jj,ibotm1) = 0. z_elem_c(ji,jj,ibotm1) = 0. z_elem_b(ji,jj,ibotm1) = 1. en(ji,jj,ibotm1) = MAX( rn_c02r * ustarb2(ji,jj), rn_emin ) END DO END DO ! CASE ( 1 ) ! Neumman boundary condition ! !CDIR NOVERRCHK DO jj = 2, jpjm1 !CDIR NOVERRCHK DO ji = fs_2, fs_jpim1 ! vector opt. ibot = mbathy(ji,jj) ibotm1 = ibot-1 ! ! Bottom level Dirichlet condition: z_elem_a(ji,jj,ibot) = 0. z_elem_c(ji,jj,ibot) = 0. z_elem_b(ji,jj,ibot) = 1. en(ji,jj,ibot) = MAX( rn_c02r * ustarb2(ji,jj), rn_emin ) ! ! Just above last level: Neumann condition 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 z_elem_c(ji,jj,ibotm1) = 0. END DO END DO ! END SELECT ! Matrix inversion (en prescribed at surface and the bottom) ! ---------------------------------------------------------- ! DO jk = 2, jpkm1 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. 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) END DO END DO END DO DO jk = 2, jpk ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. 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) END DO END DO END DO DO jk = jpk-1, 2, -1 ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. 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) END DO END DO END DO ! ! set the minimum value of tke en(:,:,:) = MAX( en(:,:,:), rn_emin ) !!----------------------------------------!! !! Solve prognostic equation for psi !! !!----------------------------------------!! ! Set psi to previous time step value ! SELECT CASE ( nn_clos ) ! CASE( 0 ) ! k-kl (Mellor-Yamada) ! DO jk = 2, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. psi(ji,jj,jk) = en(ji,jj,jk) * mxln(ji,jj,jk) ENDDO ENDDO ENDDO ! CASE( 1 ) ! k-eps ! DO jk = 2, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. psi(ji,jj,jk) = eps(ji,jj,jk) ENDDO ENDDO ENDDO ! CASE( 2 ) ! k-w ! DO jk = 2, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. psi(ji,jj,jk) = SQRT(en(ji,jj,jk)) / (rn_c0 * mxln(ji,jj,jk)) ENDDO ENDDO ENDDO ! CASE( 3 ) ! gen ! DO jk = 2, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. psi(ji,jj,jk) = rn_c02 * en(ji,jj,jk) * mxln(ji,jj,jk)**znn ENDDO ENDDO ENDDO ! END SELECT ! ! Now gls (output in psi) ! ------------------------------- ! Resolution of a tridiagonal linear system by a "methode de chasse" ! computation from level 2 to jpkm1 (e(1) already computed and e(jpk)=0 ). ! z_elem_b : diagonal z_elem_c : upper diagonal z_elem_a : lower diagonal ! Warning : after this step, en : right hand side of the matrix DO jk = 2, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ! ! psi / k zratio = psi(ji,jj,jk) / eb(ji,jj,jk) ! ! psi3+ : stable : B=-KhN²<0 => N²>0 if rn2>0 dir = 1 (stable) otherwise dir = 0 (unstable) dir = 0.5 + sign(0.5,rn2(ji,jj,jk)) ! rn_psi3 = dir*rn_psi3m+(1-dir)*rn_psi3p ! ! shear prod. - stratif. destruction prod = rn_psi1 * zratio * shear(ji,jj,jk) ! ! stratif. destruction buoy = rn_psi3 * zratio * (- avt(ji,jj,jk) * rn2(ji,jj,jk)) ! ! shear prod. - stratif. destruction diss = rn_psi2 * zratio * zwall(ji,jj,jk) * eps(ji,jj,jk) ! dir = 0.5 + sign(0.5,prod+buoy) ! dir =1(=0) if shear(ji,jj,jk)+buoy >0(<0) ! zesh2 = dir*(prod+buoy) +(1-dir)*prod ! production term zdiss = dir*(diss/psi(ji,jj,jk)) +(1-dir)*(diss-buoy)/psi(ji,jj,jk) ! dissipation term ! ! building the matrix zcof = zfact_psi * zwall_psi(ji,jj,jk) * tmask(ji,jj,jk) ! lower diagonal z_elem_a(ji,jj,jk) = zcof * ( avm (ji,jj,jk ) + avm (ji,jj,jk-1) ) & & / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk ) ) ! upper diagonal z_elem_c(ji,jj,jk) = zcof * ( avm (ji,jj,jk+1) + avm (ji,jj,jk ) ) & & / ( fse3t(ji,jj,jk ) * fse3w(ji,jj,jk) ) ! diagonal z_elem_b(ji,jj,jk) = 1. - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) & & + rdt * zdiss * tmask(ji,jj,jk) ! ! right hand side in psi psi(ji,jj,jk) = psi(ji,jj,jk) + rdt * zesh2 * tmask(ji,jj,jk) END DO END DO END DO ! z_elem_b(:,:,jpk) = 1. ! Surface boundary condition on psi ! --------------------------------- ! SELECT CASE ( nn_psibc_surf ) ! CASE ( 0 ) ! Dirichlet boundary conditions ! IF (ln_crban) THEN ! Wave induced mixing case ! ! en(1) = q2(1) = 0.5 * (15.8 * Ccb)^(2/3) * u*^2 ! ! balance between the production and the dissipation terms including the wave effect ! zdep(:,:) = rn_l_sf * zhsro(:,:) psi (:,:,1) = rn_c0**zpp * en(:,:,1)**zmm * zdep(:,:)**znn * tmask(:,:,1) z_elem_a(:,:,1) = psi(:,:,1) z_elem_c(:,:,1) = 0. z_elem_b(:,:,1) = 1. ! ! one level below zdep(:,:) = ( (zhsro(:,:) + fsdepw(:,:,2))**(zmm*rn_a_sf+znn) ) & & / zhsro(:,:)**(zmm*rn_a_sf) psi (:,:,2) = rn_sbc_psi1 * ustars2(:,:)**zmm * zdep(:,:) * tmask(:,:,1) z_elem_a(:,:,2) = 0. z_elem_c(:,:,2) = 0. z_elem_b(:,:,2) = 1. ! ELSE ! No wave induced mixing case ! ! en(1) = u*^2/C0^2 & l(1) = K*zs ! ! balance between the production and the dissipation terms ! zdep(:,:) = vkarmn * zhsro(:,:) psi (:,:,1) = rn_c0**zpp * en(:,:,1)**zmm * zdep(:,:)**znn * tmask(:,:,1) z_elem_a(:,:,1) = psi(:,:,1) z_elem_c(:,:,1) = 0. z_elem_b(:,:,1) = 1. ! ! one level below zdep(:,:) = vkarmn * ( zhsro(:,:) + fsdepw(:,:,2) ) psi (:,:,2) = rn_c0**zpp * en(:,:,1)**zmm * zdep(:,:)**znn * tmask(:,:,1) z_elem_a(:,:,2) = 0. z_elem_c(:,:,2) = 0. z_elem_b(:,:,2) = 1. ! ENDIF ! CASE ( 1 ) ! Neumann boundary condition on d(psi)/dz ! IF (ln_crban) THEN ! Wave induced mixing case ! zdep(:,:) = rn_l_sf * zhsro(:,:) psi (:,:,1) = rn_c0**zpp * en(:,:,1)**zmm * zdep(:,:)**znn * tmask(:,:,1) z_elem_a(:,:,1) = psi(:,:,1) z_elem_c(:,:,1) = 0. z_elem_b(:,:,1) = 1. ! ! Neumann condition at k=2 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b z_elem_a(:,:,2) = 0. ! ! Set psi vertical flux at the surface: zdep(:,:) = (zhsro(:,:) + fsdept(:,:,1))**(zmm*rn_a_sf+znn-1.) / zhsro(:,:)**(zmm*rn_a_sf) zflxs(:,:) = rn_sbc_psi3 * ( zwall_psi(:,:,1)*avm(:,:,1) + zwall_psi(:,:,2)*avm(:,:,2) ) & & * en(:,:,1)**zmm * zdep psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) ! ELSE ! No wave induced mixing ! zdep(:,:) = vkarmn * zhsro(:,:) psi (:,:,1) = rn_c0**zpp * en(:,:,1)**zmm * zdep(:,:)**znn * tmask(:,:,1) z_elem_a(:,:,1) = psi(:,:,1) z_elem_c(:,:,1) = 0. z_elem_b(:,:,1) = 1. ! ! Neumann condition at k=2 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b z_elem_a(ji,jj,2) = 0. ! ! Set psi vertical flux at the surface: zdep(:,:) = zhsro(:,:) + fsdept(:,:,1) zflxs(:,:) = rn_sbc_psi2 * ( avm(:,:,1) + avm(:,:,2) ) * en(:,:,1)**zmm * zdep**(znn-1.) psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) ! ENDIF ! END SELECT ! Bottom boundary condition on psi ! -------------------------------- ! SELECT CASE ( nn_psibc_bot ) ! ! CASE ( 0 ) ! Dirichlet ! ! en(ibot) = u*^2 / Co2 and mxln(ibot) = vkarmn * rn_hbro ! ! Balance between the production and the dissipation terms !CDIR NOVERRCHK DO jj = 2, jpjm1 !CDIR NOVERRCHK DO ji = fs_2, fs_jpim1 ! vector opt. ibot = mbathy(ji,jj) ibotm1 = ibot-1 zdep(ji,jj) = vkarmn * rn_hbro psi (ji,jj,ibot) = rn_c0**zpp * en(ji,jj,ibot)**zmm * zdep(ji,jj)**znn z_elem_a(ji,jj,ibot) = 0. z_elem_c(ji,jj,ibot) = 0. z_elem_b(ji,jj,ibot) = 1. ! ! Just above last level, Dirichlet condition again (GOTM like) zdep(ji,jj) = vkarmn * (rn_hbro + fse3t(ji,jj,ibotm1)) psi (ji,jj,ibotm1) = rn_c0**zpp * en(ji,jj,ibot )**zmm * zdep(ji,jj)**znn z_elem_a(ji,jj,ibotm1) = 0. z_elem_c(ji,jj,ibotm1) = 0. z_elem_b(ji,jj,ibotm1) = 1. END DO END DO ! ! CASE ( 1 ) ! Neumman boundary condition ! !CDIR NOVERRCHK DO jj = 2, jpjm1 !CDIR NOVERRCHK DO ji = fs_2, fs_jpim1 ! vector opt. ibot = mbathy(ji,jj) ibotm1 = ibot-1 ! ! Bottom level Dirichlet condition: zdep(ji,jj) = vkarmn * rn_hbro psi (ji,jj,ibot) = rn_c0**zpp * en(ji,jj,ibot)**zmm * zdep(ji,jj)**znn ! z_elem_a(ji,jj,ibot) = 0. z_elem_c(ji,jj,ibot) = 0. z_elem_b(ji,jj,ibot) = 1. ! ! Just above last level: Neumann condition with flux injection 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 z_elem_c(ji,jj,ibotm1) = 0. ! ! Set psi vertical flux at the bottom: zdep(ji,jj) = rn_hbro + 0.5*fse3t(ji,jj,ibotm1) zflxb = rn_sbc_psi2 * ( avm(ji,jj,ibot) + avm(ji,jj,ibotm1) ) * & & (0.5*(en(ji,jj,ibot)+en(ji,jj,ibotm1)))**zmm * zdep(ji,jj)**(znn-1.) psi(ji,jj,ibotm1) = psi(ji,jj,ibotm1) + zflxb / fse3w(ji,jj,ibotm1) END DO END DO ! END SELECT ! Matrix inversion ! ---------------- ! DO jk = 2, jpkm1 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. 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) END DO END DO END DO DO jk = 2, jpk ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. 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) END DO END DO END DO DO jk = jpk-1, 2, -1 ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. 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) END DO END DO END DO ! Set dissipation !---------------- SELECT CASE ( nn_clos ) ! CASE( 0 ) ! k-kl (Mellor-Yamada) ! DO jk = 1, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. eps(ji,jj,jk) = rn_c03 * en(ji,jj,jk) * en(ji,jj,jk) * SQRT(en(ji,jj,jk)) / psi(ji,jj,jk) ENDDO ENDDO ENDDO ! CASE( 1 ) ! k-eps ! DO jk = 1, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. eps(ji,jj,jk) = psi(ji,jj,jk) ENDDO ENDDO ENDDO ! CASE( 2 ) ! k-w ! DO jk = 1, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. eps(ji,jj,jk) = rn_c04 * en(ji,jj,jk) * psi(ji,jj,jk) ENDDO ENDDO ENDDO ! CASE( 3 ) ! gen ! DO jk = 1, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. eps(ji,jj,jk) = rn_c0**(3.+zpp/znn) * en(ji,jj,jk)**(1.5+zmm/znn) * psi(ji,jj,jk)**(-1./znn) ENDDO ENDDO ENDDO ! END SELECT ! Limit dissipation rate under stable stratification ! -------------------------------------------------- DO jk = 1, jpkm1 ! Note that this set boundary conditions on mxln at the same time DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ! limitation eps(ji,jj,jk) = MAX( eps(ji,jj,jk), rn_epsmin ) mxln(ji,jj,jk) = rn_c03 * en(ji,jj,jk) * SQRT(en(ji,jj,jk)) / eps(ji,jj,jk) ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated) zrn2 = MAX( rn2(ji,jj,jk), rsmall ) mxln(ji,jj,jk) = MIN( clim_galp*SQRT( 2.*en(ji,jj,jk)/zrn2 ), mxln(ji,jj,jk) ) END DO END DO END DO ! ! Stability function and vertical viscosity and diffusivity ! --------------------------------------------------------- ! SELECT CASE ( nn_stab_func ) ! CASE ( 0 , 1 ) ! Galperin or Kantha-Clayson stability functions ! DO jk = 2, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ! zcof = l²/q² zcof = mxlb(ji,jj,jk)**2. / ( 2.*eb(ji,jj,jk) ) ! Gh = -N²l²/q² gh = - rn2(ji,jj,jk) * zcof ! gh = MIN(gh, gh-(gh-rn_ghcri)**2./(gh + rn_gh0 - 2.0*rn_ghcri)) gh = MIN( gh, rn_gh0 ) gh = MAX( gh, rn_ghmin ) ! Stability functions from Kantha and Clayson (if C2=C3=0 => Galperin) sh = rn_a2*( 1.-6.*rn_a1/rn_b1 ) / ( 1.-3.*rn_a2*gh*(6.*rn_a1+rn_b2*( 1.-rn_c3 ) ) ) sm = ( rn_b1**(-1./3.) + ( 18*rn_a1*rn_a1 + 9.*rn_a1*rn_a2*(1.-rn_c2) )*sh*gh ) / (1.-9.*rn_a1*rn_a2*gh) ! ! Store stability function in avmu and avmv avmu(ji,jj,jk) = rn_c_diff * sh * tmask(ji,jj,jk) avmv(ji,jj,jk) = rn_c_diff * sm * tmask(ji,jj,jk) END DO END DO END DO ! CASE ( 2, 3 ) ! Canuto stability functions ! DO jk = 2, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ! zcof = l²/q² zcof = mxlb(ji,jj,jk)**2. / ( 2.*eb(ji,jj,jk) ) ! Gh = -N²l²/q² gh = - rn2(ji,jj,jk) * zcof ! gh = MIN(gh, gh-(gh-rn_ghcri)**2./(gh + rn_gh0 - 2.0*rn_ghcri)) gh = MIN( gh, rn_gh0 ) gh = MAX( gh, rn_ghmin ) gh = gh*rn_f6 ! Gm = M²l²/q² Shear number shr = shear(ji,jj,jk)/MAX(avm(ji,jj,jk), rsmall) gm = shr * zcof gm = MAX(gm, 1.e-10) gm = gm*rn_f6 gm = MIN ( (rn_d0 - rn_d1*gh + rn_d3*gh**2.)/(rn_d2-rn_d4*gh) , gm ) ! Stability functions from Canuto rn_cff = rn_d0 - rn_d1*gh +rn_d2*gm + rn_d3*gh**2. - rn_d4*gh*gm + rn_d5*gm**2. sm = (rn_s0 - rn_s1*gh + rn_s2*gm) / rn_cff sh = (rn_s4 - rn_s5*gh + rn_s6*gm) / rn_cff ! ! Store stability function in avmu and avmv avmu(ji,jj,jk) = rn_c_diff * sh * tmask(ji,jj,jk) avmv(ji,jj,jk) = rn_c_diff * sm * tmask(ji,jj,jk) END DO END DO END DO ! END SELECT ! Boundary conditions on stability functions for momentum (Neumann): ! Lines below are useless if GOTM style Dirichlet conditions are used DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. avmv(ji,jj,1) = rn_cm_sf / SQRT(2.) END DO END DO DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ibot=mbathy(ji,jj) ibotm1=ibot-1 avmv(ji,jj,ibot) = rn_c0 / SQRT(2.) END DO END DO ! Compute diffusivities/viscosities ! The computation below could be restrained to jk=2 to jpkm1 if GOTM style Dirichlet conditions are used DO jk = 1, jpk DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zsqen = SQRT(2.*en(ji,jj,jk)) * mxln(ji,jj,jk) zav = zsqen * avmu(ji,jj,jk) avt(ji,jj,jk) = MAX(zav,avtb(jk))*tmask(ji,jj,jk) ! apply mask for zdfmxl routine zav = zsqen * avmv(ji,jj,jk) avm(ji,jj,jk) = MAX(zav,avmb(jk)) ! Note that avm is not masked at the surface and the bottom END DO END DO END DO ! ! Lateral boundary conditions (sign unchanged) ! avt(:,:,1) = 0. CALL lbc_lnk( avm, 'W', 1. ) ; CALL lbc_lnk( avt, 'W', 1. ) DO jk = 2, jpkm1 !* vertical eddy viscosity at u- and v-points DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. 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) & & / MAX( 1., tmask(ji,jj,jk) + tmask(ji+1,jj ,jk) ) 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) & & / MAX( 1., tmask(ji,jj,jk) + tmask(ji ,jj+1,jk) ) END DO END DO END DO avmu(:,:,1) = 0. avmv(:,:,1) = 0. CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. ) ! Lateral boundary conditions IF(ln_ctl) THEN CALL prt_ctl( tab3d_1=en , clinfo1=' gls - e: ', tab3d_2=avt, clinfo2=' t: ', ovlap=1, kdim=jpk) CALL prt_ctl( tab3d_1=avmu, clinfo1=' gls - u: ', mask1=umask, & & tab3d_2=avmv, clinfo2= ' v: ', mask2=vmask, ovlap=1, kdim=jpk ) ENDIF ! END SUBROUTINE zdf_gls SUBROUTINE zdf_gls_init !!---------------------------------------------------------------------- !! *** ROUTINE zdf_gls_init *** !! !! ** Purpose : Initialization of the vertical eddy diffivity and !! viscosity when using a gls turbulent closure scheme !! !! ** Method : Read the namzdf_gls namelist and check the parameters !! called at the first timestep (nit000) !! !! ** input : Namlist namzdf_gls !! !! ** Action : Increase by 1 the nstop flag is setting problem encounter !! !!---------------------------------------------------------------------- USE dynzdf_exp USE trazdf_exp ! # if defined key_vectopt_memory INTEGER :: ji, jj, jk ! dummy loop indices # else INTEGER :: jk ! dummy loop indices # endif REAL(wp):: zcr !! NAMELIST/namzdf_gls/rn_emin, rn_epsmin, ln_length_lim, & & clim_galp, ln_crban, ln_sigpsi, & & rn_crban, rn_charn, & & nn_tkebc_surf, nn_tkebc_bot, & & nn_psibc_surf, nn_psibc_bot, & & nn_stab_func, nn_clos !!---------------------------------------------------------- ! Read Namelist namzdf_gls ! ------------------------ REWIND ( numnam ) READ ( numnam, namzdf_gls ) ! Parameter control and print ! --------------------------- ! Control print IF(lwp) THEN WRITE(numout,*) WRITE(numout,*) 'zdf_gls_init : gls turbulent closure scheme' WRITE(numout,*) '~~~~~~~~~~~~' WRITE(numout,*) ' Namelist namzdf_gls : set gls mixing parameters' WRITE(numout,*) ' minimum value of en rn_emin = ', rn_emin WRITE(numout,*) ' minimum value of eps rn_epsmin = ', rn_epsmin WRITE(numout,*) ' Surface roughness (m) hsro = ', hsro WRITE(numout,*) ' Bottom roughness (m) rn_hbro = ', rn_hbro WRITE(numout,*) ' Limit dissipation rate under stable stratification ln_length_lim = ',ln_length_lim WRITE(numout,*) ' Galperin length scale limitation coef (Standard: 0.53, Holt: 0.26) clim_galp = ', clim_galp WRITE(numout,*) ' TKE Surface boundary condition nn_tkebc_surf = ', nn_tkebc_surf WRITE(numout,*) ' TKE Bottom boundary condition nn_tkebc_bot = ', nn_tkebc_bot WRITE(numout,*) ' PSI Surface boundary condition nn_psibc_surf = ', nn_psibc_surf WRITE(numout,*) ' PSI Bottom boundary condition nn_psibc_bot = ', nn_psibc_bot WRITE(numout,*) ' Craig and Banner scheme ln_crban = ', ln_crban WRITE(numout,*) ' Modify psi Schmidt number (wb case) ln_sigpsi = ', ln_sigpsi WRITE(numout,*) ' Craig and Banner coef. rn_crban = ', rn_crban WRITE(numout,*) ' Charnock coef. rn_charn = ', rn_charn WRITE(numout,*) ' Stability functions nn_stab_func = ',nn_stab_func WRITE(numout,*) ' Closure nn_clos = ',nn_clos WRITE(numout,*) ENDIF ! Check of some namelist values IF( nn_tkebc_surf < 0 .OR. nn_tkebc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_tkebc_surf is 0 or 1' ) IF( nn_psibc_surf < 0 .OR. nn_psibc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_psibc_surf is 0 or 1' ) IF( nn_tkebc_bot < 0 .OR. nn_tkebc_bot > 1 ) CALL ctl_stop( 'bad flag: nn_tkebc_bot is 0 or 1' ) IF( nn_psibc_bot < 0 .OR. nn_psibc_bot > 1 ) CALL ctl_stop( 'bad flag: nn_psibc_bot is 0 or 1' ) IF( nn_stab_func < 0 .OR. nn_stab_func > 3) CALL ctl_stop( 'bad flag: nn_stab_func is 0, 1, 2 and 3' ) IF( nn_clos < 0 .OR. nn_clos > 3) CALL ctl_stop( 'bad flag: nn_clos is 0, 1, 2 or 3' ) ! Initialisation of the parameters for the choosen closure ! -------------------------------------------------------- ! SELECT CASE ( nn_clos ) ! CASE( 0 ) ! k-kl (Mellor-Yamada) ! IF(lwp) WRITE(numout,*) 'The choosen closure is k-kl closed to the classical Mellor-Yamada' zpp = 0. zmm = 1. znn = 1. rn_sc_tke = 1.96 rn_sc_psi = 1.96 rn_psi1 = 0.9 rn_psi3p = 1. rn_psi2 = 0.5 ! SELECT CASE ( nn_stab_func ) ! CASE( 0, 1 ) ! G88 or KC stability functions rn_psi3m = 2.53 CASE( 2 ) ! Canuto A stability functions rn_psi3m = 2.38 CASE( 3 ) ! Canuto B stability functions rn_psi3m = 2.38 ! caution : constant not identified END SELECT ! CASE( 1 ) ! k-eps ! IF(lwp) WRITE(numout,*) 'The choosen closure is k-eps' zpp = 3. zmm = 1.5 znn = -1. rn_sc_tke = 1. rn_sc_psi = 1.3 ! Schmidt number for psi rn_psi1 = 1.44 rn_psi3p = 1. rn_psi2 = 1.92 ! SELECT CASE ( nn_stab_func ) ! CASE( 0, 1 ) ! G88 or KC stability functions rn_psi3m = -0.52 CASE( 2 ) ! Canuto A stability functions rn_psi3m = -0.629 CASE( 3 ) ! Canuto B stability functions rn_psi3m = -0.566 END SELECT ! CASE( 2 ) ! k-omega ! IF(lwp) WRITE(numout,*) 'The choosen closure is k-omega' zpp = -1. zmm = 0.5 znn = -1. rn_sc_tke = 2. rn_sc_psi = 2. rn_psi1 = 0.555 rn_psi3p = 1. rn_psi2 = 0.833 ! SELECT CASE ( nn_stab_func ) ! CASE( 0, 1 ) ! G88 or KC stability functions rn_psi3m = -0.58 CASE( 2 ) ! Canuto A stability functions rn_psi3m = -0.64 CASE( 3 ) ! Canuto B stability functions rn_psi3m = -0.64 ! caution : constant not identified END SELECT ! CASE( 3 ) ! gen ! IF(lwp) WRITE(numout,*) 'The choosen closure is generic' zpp = 2. zmm = 1. znn = -0.67 rn_sc_tke = 0.8 rn_sc_psi = 1.07 rn_psi1 = 1. rn_psi3p = 1. rn_psi2 = 1.22 ! SELECT CASE ( nn_stab_func ) ! CASE( 0, 1 ) ! G88 or KC stability functions rn_psi3m = 0.1 CASE( 2 ) ! Canuto A stability functions rn_psi3m = 0.05 CASE( 3 ) ! Canuto B stability functions rn_psi3m = 0.05 ! caution : constant not identified END SELECT ! END SELECT ! Initialisation of the parameters of the stability functions ! ----------------------------------------------------------- ! SELECT CASE ( nn_stab_func ) ! CASE ( 0 ) ! Galperin stability functions ! IF(lwp) WRITE(numout,*) 'Stability functions from Galperin' rn_c2 = 0. rn_c3 = 0. rn_c_diff = 1. rn_c0 = 0.5544 rn_cm_sf = 0.9884 rn_ghmin = -0.28 rn_gh0 = 0.0233 rn_ghcri = 0.02 ! CASE ( 1 ) ! Kantha-Clayson stability functions ! IF(lwp) WRITE(numout,*) 'Stability functions from Kantha-Clayson' rn_c2 = 0.7 rn_c3 = 0.2 rn_c_diff = 1. rn_c0 = 0.5544 rn_cm_sf = 0.9884 rn_ghmin = -0.28 rn_gh0 = 0.0233 rn_ghcri = 0.02 ! CASE ( 2 ) ! Canuto A stability functions ! IF(lwp) WRITE(numout,*) 'Stability functions from Canuto A' rn_s0 = 1.5*rn_l1*rn_l5**2. rn_s1 = -rn_l4*(rn_l6+rn_l7) + 2.*rn_l4*rn_l5*(rn_l1-(1./3.)*rn_l2-rn_l3)+1.5*rn_l1*rn_l5*rn_l8 rn_s2 = -(3./8.)*rn_l1*(rn_l6**2.-rn_l7**2.) rn_s4 = 2.*rn_l5 rn_s5 = 2.*rn_l4 rn_s6 = (2./3.)*rn_l5*(3.*rn_l3**2.-rn_l2**2.)-0.5*rn_l5*rn_l1*(3.*rn_l3-rn_l2)+0.75*rn_l1*(rn_l6-rn_l7) rn_d0 = 3*rn_l5**2. rn_d1 = rn_l5*(7.*rn_l4+3.*rn_l8) rn_d2 = rn_l5**2.*(3.*rn_l3**2.-rn_l2**2.)-0.75*(rn_l6**2.-rn_l7**2.) rn_d3 = rn_l4*(4.*rn_l4+3.*rn_l8) rn_d4 = rn_l4*(rn_l2*rn_l6-3.*rn_l3*rn_l7-rn_l5*(rn_l2**2.-rn_l3**2.))+rn_l5*rn_l8*(3.*rn_l3**2.-rn_l2**2.) rn_d5 = 0.25*(rn_l2**2.-3.*rn_l3**2.)*(rn_l6**2.-rn_l7**2.) rn_c0 = 0.5268 rn_f6 = 8. / (rn_c0**6.) rn_c_diff = SQRT(2.)/(rn_c0**3.) rn_cm_sf = 0.7310 rn_ghmin = -0.28 rn_gh0 = 0.0329 rn_ghcri = 0.03 ! CASE ( 3 ) ! Canuto B stability functions ! IF(lwp) WRITE(numout,*) 'Stability functions from Canuto B' rn_s0 = 1.5*rn_m1*rn_m5**2. rn_s1 = -rn_m4*(rn_m6+rn_m7) + 2.*rn_m4*rn_m5*(rn_m1-(1./3.)*rn_m2-rn_m3)+1.5*rn_m1*rn_m5*rn_m8 rn_s2 = -(3./8.)*rn_m1*(rn_m6**2.-rn_m7**2.) rn_s4 = 2.*rn_m5 rn_s5 = 2.*rn_m4 rn_s6 = (2./3.)*rn_m5*(3.*rn_m3**2.-rn_m2**2.)-0.5*rn_m5*rn_m1*(3.*rn_m3-rn_m2)+0.75*rn_m1*(rn_m6-rn_m7) rn_d0 = 3*rn_m5**2. rn_d1 = rn_m5*(7.*rn_m4+3.*rn_m8) rn_d2 = rn_m5**2.*(3.*rn_m3**2.-rn_m2**2.)-0.75*(rn_m6**2.-rn_m7**2.) rn_d3 = rn_m4*(4.*rn_m4+3.*rn_m8) rn_d4 = rn_m4*(rn_m2*rn_m6-3.*rn_m3*rn_m7-rn_m5*(rn_m2**2.-rn_m3**2.))+rn_m5*rn_m8*(3.*rn_m3**2.-rn_m2**2.) rn_d5 = 0.25*(rn_m2**2.-3.*rn_m3**2.)*(rn_m6**2.-rn_m7**2.) rn_c0 = 0.5268 !! rn_c0 = 0.5540 (Warner ...) to verify ! rn_f6 = 8. / (rn_c0**6.) rn_c_diff = SQRT(2.)/(rn_c0**3.) rn_cm_sf = 0.7470 rn_ghmin = -0.28 rn_gh0 = 0.0444 rn_ghcri = 0.0414 ! END SELECT ! Set Schmidt number for psi diffusion ! In the wave breaking case ! See equation 13 of Carniel et al, Ocean modelling, 30, 225-239, 2009 ! or equation (17) of Burchard, JPO, 31, 3133-3145, 2001 IF ((ln_sigpsi).AND.(ln_crban)) THEN zcr = SQRT(1.5*rn_sc_tke) * rn_cm_sf /vkarmn rn_sc_psi0 = vkarmn**2/(rn_psi2*rn_cm_sf**2) * & & ( znn**2 - 4./3.*zcr*znn*zmm - 1./3.*zcr*znn & & + 2./9.*zmm*zcr**2 + 4./9.*zcr**2*zmm**2) ELSE rn_sc_psi0 = rn_sc_psi ENDIF ! Shear free turbulence parameters: ! rn_a_sf = -4.*znn*SQRT(rn_sc_tke) / ( (1.+4.*zmm)*SQRT(rn_sc_tke) & & - SQRT(rn_sc_tke + 24.*rn_sc_psi0*rn_psi2 ) ) rn_l_sf = rn_c0 * SQRT(rn_c0/rn_cm_sf) * SQRT( ( (1. + 4.*zmm + 8.*zmm**2)*rn_sc_tke & & + 12. * rn_sc_psi0*rn_psi2 - (1. + 4.*zmm) & & *SQRT(rn_sc_tke*(rn_sc_tke & & + 24.*rn_sc_psi0*rn_psi2)) ) & & /(12.*znn**2.) & & ) ! Control print ! IF(lwp) THEN WRITE(numout,*) WRITE(numout,*) 'Limit values' WRITE(numout,*) '~~~~~~~~~~~~' WRITE(numout,*) 'Parameter m = ',zmm WRITE(numout,*) 'Parameter n = ',znn WRITE(numout,*) 'Parameter p = ',zpp WRITE(numout,*) 'rn_psi1 = ',rn_psi1 WRITE(numout,*) 'rn_psi2 = ',rn_psi2 WRITE(numout,*) 'rn_psi3m = ',rn_psi3m WRITE(numout,*) 'rn_psi3p = ',rn_psi3p WRITE(numout,*) 'rn_sc_tke = ',rn_sc_tke WRITE(numout,*) 'rn_sc_psi = ',rn_sc_psi WRITE(numout,*) 'rn_sc_psi0 = ',rn_sc_psi0 WRITE(numout,*) 'rn_c0 = ',rn_c0 WRITE(numout,*) WRITE(numout,*) 'Shear free turbulence parameters:' WRITE(numout,*) 'rn_cm_sf = ',rn_cm_sf WRITE(numout,*) 'rn_a_sf = ',rn_a_sf WRITE(numout,*) 'rn_l_sf = ',rn_l_sf WRITE(numout,*) ENDIF ! Constants initialization rn_c02r = 1. / rn_c0**2. rn_c02 = rn_c0**2._wp rn_c03 = rn_c0**3._wp rn_c04 = rn_c0**4._wp rn_c03_sqrt2_galp = rn_c03 / SQRT(2._wp) / clim_galp rn_sbc_mb = 0.5 * (15.8*rn_crban)**(2./3.) ! Surf. bound. cond. from Mellor and Blumberg rn_sbc_std = 3.75 ! Surf. bound. cond. standard (prod=diss) rn_sbc_tke1 = (-rn_sc_tke*rn_crban/(rn_cm_sf*rn_a_sf*rn_l_sf))**(2./3.) ! k_eps = 53. Dirichlet + Wave breaking rn_sbc_tke2 = 0.5 / rau0 rn_sbc_tke3 = rdt * rn_crban ! Neumann + Wave breaking rn_sbc_zs = rn_charn/grav ! Charnock formula rn_sbc_psi1 = rn_c0**zpp * rn_sbc_tke1**zmm * rn_l_sf**znn ! Dirichlet + Wave breaking rn_sbc_psi2 = -0.5 * rdt * rn_c0**zpp * znn * vkarmn**znn / rn_sc_psi ! Neumann + NO Wave breaking rn_sbc_psi3 = -0.5 * rdt * rn_c0**zpp * rn_l_sf**znn / rn_sc_psi * (znn + zmm*rn_a_sf) ! Neumann + Wave breaking zfact_tke = -0.5 / rn_sc_tke * rdt ! Cst used for the Diffusion term of tke zfact_psi = -0.5 / rn_sc_psi * rdt ! Cst used for the Diffusion term of tke ! Wall proximity function zwall (:,:,:) = 1._wp * tmask(:,:,:) ! !* set vertical eddy coef. to the background value DO jk = 1, jpk avt (:,:,jk) = avtb(jk) * tmask(:,:,jk) avm (:,:,jk) = avmb(jk) * tmask(:,:,jk) avmu(:,:,jk) = avmb(jk) * umask(:,:,jk) avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk) END DO ! !* read or initialize all required files CALL gls_rst( nit000, 'READ' ) ! END SUBROUTINE zdf_gls_init SUBROUTINE gls_rst( kt, cdrw ) !!--------------------------------------------------------------------- !! *** ROUTINE ts_rst *** !! !! ** Purpose : Read or write TKE file (en) in restart file !! !! ** Method : use of IOM library !! if the restart does not contain TKE, en is either !! set to rn_emin or recomputed (nn_igls/=0) !!---------------------------------------------------------------------- INTEGER , INTENT(in) :: kt ! ocean time-step CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag ! INTEGER :: jit, jk ! dummy loop indices INTEGER :: id1, id2, id3, id4, id5, id6, id7, id8 INTEGER :: ji, jj, ikbu, ikbv, ikbum1, ikbvm1 REAL(wp):: cbx, cby !!---------------------------------------------------------------------- ! IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise ! ! --------------- IF( ln_rstart ) THEN !* Read the restart file id1 = iom_varid( numror, 'en' , ldstop = .FALSE. ) id2 = iom_varid( numror, 'avt' , ldstop = .FALSE. ) id3 = iom_varid( numror, 'avm' , ldstop = .FALSE. ) id4 = iom_varid( numror, 'avmu' , ldstop = .FALSE. ) id5 = iom_varid( numror, 'avmv' , ldstop = .FALSE. ) id6 = iom_varid( numror, 'mxln' , ldstop = .FALSE. ) id7 = iom_varid( numror, 'wbotu', ldstop = .FALSE. ) id8 = iom_varid( numror, 'wbotv', ldstop = .FALSE. ) ! IF( MIN( id1, id2, id3, id4, id5, id6, id7, id8 ) > 0 ) THEN ! all required arrays exist CALL iom_get( numror, jpdom_autoglo, 'en' , en ) CALL iom_get( numror, jpdom_autoglo, 'avt' , avt ) CALL iom_get( numror, jpdom_autoglo, 'avm' , avm ) CALL iom_get( numror, jpdom_autoglo, 'avmu' , avmu ) CALL iom_get( numror, jpdom_autoglo, 'avmv' , avmv ) CALL iom_get( numror, jpdom_autoglo, 'mxln' , mxln ) CALL iom_get( numror, jpdom_autoglo, 'wbotu' , wbotu ) CALL iom_get( numror, jpdom_autoglo, 'wbotv' , wbotv ) ELSE IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without gls scheme, en and mxln computed by iterative loop' IF(lwp) WRITE(numout,*) ' ===>>>> : The bottom stresses are estimated' en (:,:,:) = rn_emin mxln(:,:,:) = 0.001 ! Initialize bottom stresses DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ikbu = MIN( mbathy(ji+1,jj ), mbathy(ji,jj) ) ikbum1 = MAX( ikbu-1, 1 ) ikbv = MIN( mbathy(ji,jj+1), mbathy(ji,jj) ) ikbvm1 = MAX( ikbv-1, 1 ) cbx = avmu(ji,jj,ikbu) / fse3uw(ji,jj,ikbu) cby = avmv(ji,jj,ikbv) / fse3vw(ji,jj,ikbv) wbotu(ji,jj) = -cbx * un(ji,jj,ikbum1)*umask(ji,jj,1) wbotv(ji,jj) = -cby * vn(ji,jj,ikbvm1)*vmask(ji,jj,1) END DO END DO DO jit = nit000 + 1, nit000 + 10 ; CALL zdf_gls( jit ) ; END DO ENDIF ELSE !* Start from rest IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of en and mxln by background values' IF(lwp) WRITE(numout,*) ' ===>>>> : The bottom stresses are estimated' en (:,:,:) = rn_emin mxln(:,:,:) = 0.001 ! Initialize bottom stresses DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ikbu = MIN( mbathy(ji+1,jj ), mbathy(ji,jj) ) ikbum1 = MAX( ikbu-1, 1 ) ikbv = MIN( mbathy(ji,jj+1), mbathy(ji,jj) ) ikbvm1 = MAX( ikbv-1, 1 ) cbx = avmu(ji,jj,ikbu) / fse3uw(ji,jj,ikbu) cby = avmv(ji,jj,ikbv) / fse3vw(ji,jj,ikbv) wbotu(ji,jj) = -cbx * un(ji,jj,ikbum1)*umask(ji,jj,1) wbotv(ji,jj) = -cby * vn(ji,jj,ikbvm1)*vmask(ji,jj,1) END DO END DO ENDIF ! ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file ! ! ------------------- IF(lwp) WRITE(numout,*) '---- gls-rst ----' CALL iom_rstput( kt, nitrst, numrow, 'en' , en ) CALL iom_rstput( kt, nitrst, numrow, 'avt' , avt ) CALL iom_rstput( kt, nitrst, numrow, 'avm' , avm ) CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu ) CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv ) CALL iom_rstput( kt, nitrst, numrow, 'mxln' , mxln ) ! ENDIF ! END SUBROUTINE gls_rst #else !!---------------------------------------------------------------------- !! Dummy module : NO TKE scheme !!---------------------------------------------------------------------- LOGICAL, PUBLIC, PARAMETER :: lk_zdfgls = .FALSE. !: TKE flag CONTAINS SUBROUTINE zdf_gls( kt ) ! Empty routine WRITE(*,*) 'zdf_gls: You should not have seen this print! error?', kt END SUBROUTINE zdf_gls #endif !!====================================================================== END MODULE zdfgls