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 !! 3.3 ! 2010-10 (C. Bricaud) add in the reference !!---------------------------------------------------------------------- #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 zdf_gls_init ! 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 under stable stratif. (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) :: rn_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) :: rcm_sf = 0.73_wp ! Shear free turbulence parameters REAL(wp) :: ra_sf = -2.0_wp ! Must be negative -2 < ra_sf < -1 REAL(wp) :: rl_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) !!-------------------------------------------------------------------- ! 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) = rsbc_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(rsbc_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) = rc03 * 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. + re2 * ( 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 rsc_psi*zwall/rsc_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 rsc_psi/rsc_psi0 IF (ln_sigpsi) THEN zsigpsi = MIN(1., zesh2/eps(ji,jj,jk)) ! 0. <= zsigpsi <= 1. zwall_psi(ji,jj,jk) = rsc_psi / (zsigpsi * rsc_psi + (1.-zsigpsi) * rsc_psi0 / MAX(zwall(ji,jj,jk),1.)) ELSE zwall_psi(ji,jj,jk) = 1. ENDIF ! ! building the matrix zcof = rfact_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) = rsc_psi/rsc_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( rsbc_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( rsbc_tke1 * ustars2(:,:) * ( (zhsro(:,:)+fsdepw(:,:,2))/zhsro(:,:) )**ra_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( rc02r * 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( rc02r * 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( rsbc_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(:,:) = rsbc_tke3 * ustars2(:,:)**1.5 * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5*ra_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( rc02r * 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( rc02r * 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( rc02r * 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( rc02r * 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)) / (rc0 * 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) = rc02 * en(ji,jj,jk) * mxln(ji,jj,jk)**rnn 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)) ! rpsi3 = dir*rpsi3m+(1-dir)*rpsi3p ! ! shear prod. - stratif. destruction prod = rpsi1 * zratio * shear(ji,jj,jk) ! ! stratif. destruction buoy = rpsi3 * zratio * (- avt(ji,jj,jk) * rn2(ji,jj,jk)) ! ! shear prod. - stratif. destruction diss = rpsi2 * 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 = rfact_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(:,:) = rl_sf * zhsro(:,:) psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * 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))**(rmm*ra_sf+rnn) ) & & / zhsro(:,:)**(rmm*ra_sf) psi (:,:,2) = rsbc_psi1 * ustars2(:,:)**rmm * 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) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * 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) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * 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(:,:) = rl_sf * zhsro(:,:) psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * 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))**(rmm*ra_sf+rnn-1.) / zhsro(:,:)**(rmm*ra_sf) zflxs(:,:) = rsbc_psi3 * ( zwall_psi(:,:,1)*avm(:,:,1) + zwall_psi(:,:,2)*avm(:,:,2) ) & & * en(:,:,1)**rmm * zdep psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) ! ELSE ! No wave induced mixing ! zdep(:,:) = vkarmn * zhsro(:,:) psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * 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(:,:) = rsbc_psi2 * ( avm(:,:,1) + avm(:,:,2) ) * en(:,:,1)**rmm * zdep**(rnn-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) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn 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) = rc0**rpp * en(ji,jj,ibot )**rmm * zdep(ji,jj)**rnn 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) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn ! 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 = rsbc_psi2 * ( avm(ji,jj,ibot) + avm(ji,jj,ibotm1) ) * & & (0.5*(en(ji,jj,ibot)+en(ji,jj,ibotm1)))**rmm * zdep(ji,jj)**(rnn-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) = rc03 * 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) = rc04 * 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) = rc0**(3.+rpp/rnn) * en(ji,jj,jk)**(1.5+rmm/rnn) * psi(ji,jj,jk)**(-1./rnn) 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) = rc03 * 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( rn_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, rgh0 ) gh = MAX( gh, rghmin ) ! Stability functions from Kantha and Clayson (if C2=C3=0 => Galperin) sh = ra2*( 1.-6.*ra1/rb1 ) / ( 1.-3.*ra2*gh*(6.*ra1+rb2*( 1.-rc3 ) ) ) sm = ( rb1**(-1./3.) + ( 18*ra1*ra1 + 9.*ra1*ra2*(1.-rc2) )*sh*gh ) / (1.-9.*ra1*ra2*gh) ! ! Store stability function in avmu and avmv avmu(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk) avmv(ji,jj,jk) = rc_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, rgh0 ) gh = MAX( gh, rghmin ) gh = gh*rf6 ! 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*rf6 gm = MIN ( (rd0 - rd1*gh + rd3*gh**2.)/(rd2-rd4*gh) , gm ) ! Stability functions from Canuto rcff = rd0 - rd1*gh +rd2*gm + rd3*gh**2. - rd4*gh*gm + rd5*gm**2. sm = (rs0 - rs1*gh + rs2*gm) / rcff sh = (rs4 - rs5*gh + rs6*gm) / rcff ! ! Store stability function in avmu and avmv avmu(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk) avmv(ji,jj,jk) = rc_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) = rcm_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) = rc0 / 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 !! !!---------------------------------------------------------------------- INTEGER :: jk ! dummy loop indices REAL(wp):: zcr ! local scalar !! NAMELIST/namzdf_gls/rn_emin, rn_epsmin, ln_length_lim, & & rn_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) rn_clim_galp = ', rn_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' rpp = 0. rmm = 1. rnn = 1. rsc_tke = 1.96 rsc_psi = 1.96 rpsi1 = 0.9 rpsi3p = 1. rpsi2 = 0.5 ! SELECT CASE ( nn_stab_func ) ! CASE( 0, 1 ) ! G88 or KC stability functions rpsi3m = 2.53 CASE( 2 ) ! Canuto A stability functions rpsi3m = 2.38 CASE( 3 ) ! Canuto B stability functions rpsi3m = 2.38 ! caution : constant not identified END SELECT ! CASE( 1 ) ! k-eps ! IF(lwp) WRITE(numout,*) 'The choosen closure is k-eps' rpp = 3. rmm = 1.5 rnn = -1. rsc_tke = 1. rsc_psi = 1.3 ! Schmidt number for psi rpsi1 = 1.44 rpsi3p = 1. rpsi2 = 1.92 ! SELECT CASE ( nn_stab_func ) ! CASE( 0, 1 ) ! G88 or KC stability functions rpsi3m = -0.52 CASE( 2 ) ! Canuto A stability functions rpsi3m = -0.629 CASE( 3 ) ! Canuto B stability functions rpsi3m = -0.566 END SELECT ! CASE( 2 ) ! k-omega ! IF(lwp) WRITE(numout,*) 'The choosen closure is k-omega' rpp = -1. rmm = 0.5 rnn = -1. rsc_tke = 2. rsc_psi = 2. rpsi1 = 0.555 rpsi3p = 1. rpsi2 = 0.833 ! SELECT CASE ( nn_stab_func ) ! CASE( 0, 1 ) ! G88 or KC stability functions rpsi3m = -0.58 CASE( 2 ) ! Canuto A stability functions rpsi3m = -0.64 CASE( 3 ) ! Canuto B stability functions rpsi3m = -0.64 ! caution : constant not identified END SELECT ! CASE( 3 ) ! gen ! IF(lwp) WRITE(numout,*) 'The choosen closure is generic' rpp = 2. rmm = 1. rnn = -0.67 rsc_tke = 0.8 rsc_psi = 1.07 rpsi1 = 1. rpsi3p = 1. rpsi2 = 1.22 ! SELECT CASE ( nn_stab_func ) ! CASE( 0, 1 ) ! G88 or KC stability functions rpsi3m = 0.1 CASE( 2 ) ! Canuto A stability functions rpsi3m = 0.05 CASE( 3 ) ! Canuto B stability functions rpsi3m = 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' rc2 = 0. rc3 = 0. rc_diff = 1. rc0 = 0.5544 rcm_sf = 0.9884 rghmin = -0.28 rgh0 = 0.0233 rghcri = 0.02 ! CASE ( 1 ) ! Kantha-Clayson stability functions ! IF(lwp) WRITE(numout,*) 'Stability functions from Kantha-Clayson' rc2 = 0.7 rc3 = 0.2 rc_diff = 1. rc0 = 0.5544 rcm_sf = 0.9884 rghmin = -0.28 rgh0 = 0.0233 rghcri = 0.02 ! CASE ( 2 ) ! Canuto A stability functions ! IF(lwp) WRITE(numout,*) 'Stability functions from Canuto A' rs0 = 1.5*rl1*rl5**2. rs1 = -rl4*(rl6+rl7) + 2.*rl4*rl5*(rl1-(1./3.)*rl2-rl3)+1.5*rl1*rl5*rl8 rs2 = -(3./8.)*rl1*(rl6**2.-rl7**2.) rs4 = 2.*rl5 rs5 = 2.*rl4 rs6 = (2./3.)*rl5*(3.*rl3**2.-rl2**2.)-0.5*rl5*rl1*(3.*rl3-rl2)+0.75*rl1*(rl6-rl7) rd0 = 3*rl5**2. rd1 = rl5*(7.*rl4+3.*rl8) rd2 = rl5**2.*(3.*rl3**2.-rl2**2.)-0.75*(rl6**2.-rl7**2.) rd3 = rl4*(4.*rl4+3.*rl8) rd4 = rl4*(rl2*rl6-3.*rl3*rl7-rl5*(rl2**2.-rl3**2.))+rl5*rl8*(3.*rl3**2.-rl2**2.) rd5 = 0.25*(rl2**2.-3.*rl3**2.)*(rl6**2.-rl7**2.) rc0 = 0.5268 rf6 = 8. / (rc0**6.) rc_diff = SQRT(2.)/(rc0**3.) rcm_sf = 0.7310 rghmin = -0.28 rgh0 = 0.0329 rghcri = 0.03 ! CASE ( 3 ) ! Canuto B stability functions ! IF(lwp) WRITE(numout,*) 'Stability functions from Canuto B' rs0 = 1.5*rm1*rm5**2. rs1 = -rm4*(rm6+rm7) + 2.*rm4*rm5*(rm1-(1./3.)*rm2-rm3)+1.5*rm1*rm5*rm8 rs2 = -(3./8.)*rm1*(rm6**2.-rm7**2.) rs4 = 2.*rm5 rs5 = 2.*rm4 rs6 = (2./3.)*rm5*(3.*rm3**2.-rm2**2.)-0.5*rm5*rm1*(3.*rm3-rm2)+0.75*rm1*(rm6-rm7) rd0 = 3*rm5**2. rd1 = rm5*(7.*rm4+3.*rm8) rd2 = rm5**2.*(3.*rm3**2.-rm2**2.)-0.75*(rm6**2.-rm7**2.) rd3 = rm4*(4.*rm4+3.*rm8) rd4 = rm4*(rm2*rm6-3.*rm3*rm7-rm5*(rm2**2.-rm3**2.))+rm5*rm8*(3.*rm3**2.-rm2**2.) rd5 = 0.25*(rm2**2.-3.*rm3**2.)*(rm6**2.-rm7**2.) rc0 = 0.5268 !! rc0 = 0.5540 (Warner ...) to verify ! rf6 = 8. / (rc0**6.) rc_diff = SQRT(2.)/(rc0**3.) rcm_sf = 0.7470 rghmin = -0.28 rgh0 = 0.0444 rghcri = 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*rsc_tke) * rcm_sf /vkarmn rsc_psi0 = vkarmn**2/(rpsi2*rcm_sf**2) * & & ( rnn**2 - 4./3.*zcr*rnn*rmm - 1./3.*zcr*rnn & & + 2./9.*rmm*zcr**2 + 4./9.*zcr**2*rmm**2) ELSE rsc_psi0 = rsc_psi ENDIF ! Shear free turbulence parameters: ! ra_sf = -4.*rnn*SQRT(rsc_tke) / ( (1.+4.*rmm)*SQRT(rsc_tke) & & - SQRT(rsc_tke + 24.*rsc_psi0*rpsi2 ) ) rl_sf = rc0 * SQRT(rc0/rcm_sf) * SQRT( ( (1. + 4.*rmm + 8.*rmm**2)*rsc_tke & & + 12. * rsc_psi0*rpsi2 - (1. + 4.*rmm) & & *SQRT(rsc_tke*(rsc_tke & & + 24.*rsc_psi0*rpsi2)) ) & & /(12.*rnn**2.) & & ) ! Control print ! IF(lwp) THEN WRITE(numout,*) WRITE(numout,*) 'Limit values' WRITE(numout,*) '~~~~~~~~~~~~' WRITE(numout,*) 'Parameter m = ',rmm WRITE(numout,*) 'Parameter n = ',rnn WRITE(numout,*) 'Parameter p = ',rpp WRITE(numout,*) 'rpsi1 = ',rpsi1 WRITE(numout,*) 'rpsi2 = ',rpsi2 WRITE(numout,*) 'rpsi3m = ',rpsi3m WRITE(numout,*) 'rpsi3p = ',rpsi3p WRITE(numout,*) 'rsc_tke = ',rsc_tke WRITE(numout,*) 'rsc_psi = ',rsc_psi WRITE(numout,*) 'rsc_psi0 = ',rsc_psi0 WRITE(numout,*) 'rc0 = ',rc0 WRITE(numout,*) WRITE(numout,*) 'Shear free turbulence parameters:' WRITE(numout,*) 'rcm_sf = ',rcm_sf WRITE(numout,*) 'ra_sf = ',ra_sf WRITE(numout,*) 'rl_sf = ',rl_sf WRITE(numout,*) ENDIF ! Constants initialization rc02r = 1. / rc0**2. rc02 = rc0**2._wp rc03 = rc0**3._wp rc04 = rc0**4._wp rc03_sqrt2_galp = rc03 / SQRT(2._wp) / rn_clim_galp rsbc_mb = 0.5 * (15.8*rn_crban)**(2./3.) ! Surf. bound. cond. from Mellor and Blumberg rsbc_std = 3.75 ! Surf. bound. cond. standard (prod=diss) rsbc_tke1 = (-rsc_tke*rn_crban/(rcm_sf*ra_sf*rl_sf))**(2./3.) ! k_eps = 53. Dirichlet + Wave breaking rsbc_tke2 = 0.5 / rau0 rsbc_tke3 = rdt * rn_crban ! Neumann + Wave breaking rsbc_zs = rn_charn/grav ! Charnock formula rsbc_psi1 = rc0**rpp * rsbc_tke1**rmm * rl_sf**rnn ! Dirichlet + Wave breaking rsbc_psi2 = -0.5 * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking rsbc_psi3 = -0.5 * rdt * rc0**rpp * rl_sf**rnn / rsc_psi * (rnn + rmm*ra_sf) ! Neumann + Wave breaking rfact_tke = -0.5 / rsc_tke * rdt ! Cst used for the Diffusion term of tke rfact_psi = -0.5 / rsc_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 SUBROUTINE gls_rst( kt,cdrw ) ! Empty routine INTEGER , INTENT(in) :: kt ! ocean time-step CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag WRITE(*,*) 'gls_Rst: You should not have seen this print! error?', kt, cdrw END SUBROUTINE gls_rst #endif !!====================================================================== END MODULE zdfgls