!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE yamada4(ni, nsrf, ngrid, dt, g, rconst, plev, temp, zlev, zlay, u, v, teta, & cd, tke, eps, km, kn, kq, ustar, iflag_pbl, drgpro) USE dimphy, only : klev,klon USE phys_local_var_mod, only: wprime USE yamada_ini_mod, only : new_yamada4,yamada4_num,hboville USE yamada_ini_mod, only : prt_level, lunout,pbl_lmixmin_alpha,b1,kap,viscom,viscoh USE yamada_ini_mod, only : ric, yun,ydeux,lmixmin,iflag_vdif_q2 IMPLICIT NONE ! ************************************************************************************************ ! ! yamada4: subroutine qui calcule le transfert turbulent avec une fermeture d'ordre 2 ou 2.5 ! ! Reference: Simulation of nocturnal drainage flows by a q2l Turbulence Closure Model ! Yamada T. ! J. Atmos. Sci, 40, 91-106, 1983 ! !************************************************************************************************ ! Input : !'====== ! ni: indice horizontal sur la grille de base, non restreinte ! nsrf: type de surface ! ngrid: nombre de mailles concern??es sur l'horizontal ! dt : pas de temps ! g : g ! rconst: constante de l'air sec ! zlev : altitude a chaque niveau (interface inferieure de la couche ! de meme indice) ! zlay : altitude au centre de chaque couche ! u,v : vitesse au centre de chaque couche ! (en entree : la valeur au debut du pas de temps) ! teta : temperature potentielle virtuelle au centre de chaque couche ! (en entree : la valeur au debut du pas de temps) ! cd : cdrag pour la quantit?? de mouvement ! (en entree : la valeur au debut du pas de temps) ! ustar: vitesse de friction calcul??e par une formule diagnostique ! iflag_pbl: flag pour choisir des options du sch??ma de turbulence ! ! iflag_pbl doit valoir entre 6 et 9 ! l=6, on prend systematiquement une longueur d'equilibre ! iflag_pbl=6 : MY 2.0 ! iflag_pbl=7 : MY 2.0.Fournier ! iflag_pbl=8/9 : MY 2.5 ! iflag_pbl=8 with special obsolete treatments for convergence ! with Cmpi5 NPv3.1 simulations ! iflag_pbl=10/11 : New scheme M2 and N2 explicit and dissiptation exact ! iflag_pbl=12 = 11 with vertical diffusion off q2 ! ! 2013/04/01 (FH hourdin@lmd.jussieu.fr) ! Correction for very stable PBLs (iflag_pbl=10 and 11) ! iflag_pbl=8 converges numerically with NPv3.1 ! iflag_pbl=11 -> the model starts with NP from start files created by ce0l ! -> the model can run with longer time-steps. ! 2016/11/30 (EV etienne.vignon@lmd.ipsl.fr) ! On met tke (=q2/2) en entr??e plut??t que q2 ! On corrige l'update de la tke ! 2020/10/01 (EV) ! On ajoute la dissipation de la TKE en diagnostique de sortie ! ! Inpout/Output : !============== ! tke : tke au bas de chaque couche ! (en entree : la valeur au debut du pas de temps) ! (en sortie : la valeur a la fin du pas de temps) ! Outputs: !========== ! eps: tke dissipation rate ! km : diffusivite turbulente de quantite de mouvement (au bas de chaque ! couche) ! (en sortie : la valeur a la fin du pas de temps) ! kn : diffusivite turbulente des scalaires (au bas de chaque couche) ! (en sortie : la valeur a la fin du pas de temps) ! !....................................................................... !======================================================================= ! Declarations: !======================================================================= ! Inputs/Outputs !---------------- REAL dt, g, rconst REAL plev(klon, klev+1), temp(klon, klev) REAL ustar(klon) REAL kmin, qmin, pblhmin(klon), coriol(klon) REAL zlev(klon, klev+1) REAL zlay(klon, klev) REAL u(klon, klev) REAL v(klon, klev) REAL teta(klon, klev) REAL cd(klon) REAL tke(klon, klev+1) REAL eps(klon,klev+1) REAL unsdz(klon, klev) REAL unsdzdec(klon, klev+1) REAL kn(klon, klev+1) REAL km(klon, klev+1) INTEGER iflag_pbl, ngrid, nsrf INTEGER ni(klon) !FC REAL drgpro(klon,klev) REAL winds(klon,klev) ! Local !------- REAL q2(klon, klev+1) REAL kmpre(klon, klev+1), tmp2, qpre REAL mpre(klon, klev+1) REAL kq(klon, klev+1) REAL ff(klon, klev+1), delta(klon, klev+1) REAL aa(klon, klev+1), aa0, aa1 INTEGER nlay, nlev INTEGER ig, jg, k REAL ri, zrif, zalpha, zsm, zsn REAL rif(klon, klev+1), sm(klon, klev+1), alpha(klon, klev) REAL m2(klon, klev+1), dz(klon, klev+1), zq, n2(klon, klev+1) REAL dtetadz(klon, klev+1) REAL m2cstat, mcstat, kmcstat REAL l(klon, klev+1) REAL zz(klon, klev+1) INTEGER iter REAL dissip(klon,klev), tkeprov,tkeexp, shear(klon,klev), buoy(klon,klev) REAL :: disseff REAL,SAVE :: rifc !$OMP THREADPRIVATE(rifc) REAL,SAVE :: seuilsm, seuilalpha !$OMP THREADPRIVATE(seuilsm, seuilalpha) REAL frif, falpha, fsm REAL rino(klon, klev+1), smyam(klon, klev), styam(klon, klev), & lyam(klon, klev), knyam(klon, klev), w2yam(klon, klev), t2yam(klon, klev) CHARACTER (len = 20) :: modname = 'yamada4' CHARACTER (len = 80) :: abort_message ! Fonctions utilis??es !-------------------- frif(ri) = 0.6588*(ri+0.1776-sqrt(ri*ri-0.3221*ri+0.03156)) falpha(ri) = 1.318*(0.2231-ri)/(0.2341-ri) fsm(ri) = 1.96*(0.1912-ri)*(0.2341-ri)/((1.-ri)*(0.2231-ri)) IF (new_yamada4) THEN ! Corrections et reglages issus du travail de these d'Etienne Vignon. rifc=frif(ric) seuilsm=fsm(frif(ric)) seuilalpha=falpha(frif(ric)) ELSE rifc=0.191 seuilalpha=1.12 seuilsm=0.085 ENDIF !=============================================================================== ! Flags, tests et ??valuations de constantes !=============================================================================== ! On utilise ou non la routine de Holstalg Boville pour les cas tres stables IF (.NOT. (iflag_pbl>=6 .AND. iflag_pbl<=12)) THEN abort_message='probleme de coherence dans appel a MY' CALL abort_physic(modname,abort_message,1) END IF nlay = klev nlev = klev + 1 !======================================================================== ! Calcul des increments verticaux !========================================================================= ! Attention: zlev n'est pas declare a nlev DO ig = 1, ngrid zlev(ig, nlev) = zlay(ig, nlay) + (zlay(ig,nlay)-zlev(ig,nlev-1)) END DO DO k = 1, nlay DO ig = 1, ngrid unsdz(ig, k) = 1.E+0/(zlev(ig,k+1)-zlev(ig,k)) END DO END DO DO ig = 1, ngrid unsdzdec(ig, 1) = 1.E+0/(zlay(ig,1)-zlev(ig,1)) END DO DO k = 2, nlay DO ig = 1, ngrid unsdzdec(ig, k) = 1.E+0/(zlay(ig,k)-zlay(ig,k-1)) END DO END DO DO ig = 1, ngrid unsdzdec(ig, nlay+1) = 1.E+0/(zlev(ig,nlay+1)-zlay(ig,nlay)) END DO !========================================================================= ! Richardson number and stability functions !========================================================================= ! initialize arrays: m2(1:ngrid, :) = 0.0 sm(1:ngrid, :) = 0.0 rif(1:ngrid, :) = 0.0 !------------------------------------------------------------ DO k = 2, klev DO ig = 1, ngrid dz(ig, k) = zlay(ig, k) - zlay(ig, k-1) m2(ig, k) = ((u(ig,k)-u(ig,k-1))**2+(v(ig,k)-v(ig, & k-1))**2)/(dz(ig,k)*dz(ig,k)) dtetadz(ig, k) = (teta(ig,k)-teta(ig,k-1))/dz(ig, k) n2(ig, k) = g*2.*dtetadz(ig, k)/(teta(ig,k-1)+teta(ig,k)) ri = n2(ig, k)/max(m2(ig,k), 1.E-10) IF (ri<ric) THEN rif(ig, k) = frif(ri) ELSE rif(ig, k) = rifc END IF if (new_yamada4) then alpha(ig, k) = max(falpha(rif(ig,k)),seuilalpha) sm(ig, k) = max(fsm(rif(ig,k)),seuilsm) else IF (rif(ig,k)<0.16) THEN alpha(ig, k) = falpha(rif(ig,k)) sm(ig, k) = fsm(rif(ig,k)) ELSE alpha(ig, k) = seuilalpha sm(ig, k) = seuilsm END IF end if zz(ig, k) = b1*m2(ig, k)*(1.-rif(ig,k))*sm(ig, k) END DO END DO !======================================================================= ! DIFFERENT TYPES DE SCHEMA de YAMADA !======================================================================= ! On commence par calculer q2 a partir de la tke IF (new_yamada4) THEN DO k=1,klev+1 q2(1:ngrid,k)=tke(1:ngrid,k)*ydeux ENDDO ELSE DO k=1,klev+1 q2(1:ngrid,k)=tke(1:ngrid,k) ENDDO ENDIF ! ==================================================================== ! Computing the mixing length ! ==================================================================== CALL mixinglength(ni,nsrf,ngrid,iflag_pbl,pbl_lmixmin_alpha,lmixmin,zlay,zlev,u,v,q2,n2, l) !-------------- ! Yamada 2.0 !-------------- IF (iflag_pbl==6) THEN DO k = 2, klev q2(1:ngrid, k) = l(1:ngrid, k)**2*zz(1:ngrid, k) END DO !------------------ ! Yamada 2.Fournier !------------------ ELSE IF (iflag_pbl==7) THEN ! Calcul de l, km, au pas precedent !.................................... DO k = 2, klev DO ig = 1, ngrid delta(ig, k) = q2(ig, k)/(l(ig,k)**2*sm(ig,k)) kmpre(ig, k) = l(ig, k)*sqrt(q2(ig,k))*sm(ig, k) mpre(ig, k) = sqrt(m2(ig,k)) END DO END DO DO k = 2, klev - 1 DO ig = 1, ngrid m2cstat = max(alpha(ig,k)*n2(ig,k)+delta(ig,k)/b1, 1.E-12) mcstat = sqrt(m2cstat) ! Puis on ecrit la valeur de q qui annule l'equation de m supposee en q3 !......................................................................... IF (k==2) THEN kmcstat = 1.E+0/mcstat*(unsdz(ig,k)*kmpre(ig,k+1)*mpre(ig,k+1)+ & unsdz(ig,k-1)*cd(ig)*(sqrt(u(ig,3)**2+v(ig,3)**2)-mcstat/unsdzdec & (ig,k)-mpre(ig,k+1)/unsdzdec(ig,k+1))**2)/(unsdz(ig,k)+unsdz(ig,k & -1)) ELSE kmcstat = 1.E+0/mcstat*(unsdz(ig,k)*kmpre(ig,k+1)*mpre(ig,k+1)+ & unsdz(ig,k-1)*kmpre(ig,k-1)*mpre(ig,k-1))/ & (unsdz(ig,k)+unsdz(ig,k-1)) END IF tmp2 = kmcstat/(sm(ig,k)/q2(ig,k))/l(ig, k) q2(ig, k) = max(tmp2, 1.E-12)**(2./3.) END DO END DO ! ------------------------ ! Yamada 2.5 a la Didi !------------------------- ELSE IF (iflag_pbl==8 .OR. iflag_pbl==9) THEN ! Calcul de l, km, au pas precedent !.................................... DO k = 2, klev DO ig = 1, ngrid delta(ig, k) = q2(ig, k)/(l(ig,k)**2*sm(ig,k)) IF (delta(ig,k)<1.E-20) THEN delta(ig, k) = 1.E-20 END IF km(ig, k) = l(ig, k)*sqrt(q2(ig,k))*sm(ig, k) aa0 = (m2(ig,k)-alpha(ig,k)*n2(ig,k)-delta(ig,k)/b1) aa1 = (m2(ig,k)*(1.-rif(ig,k))-delta(ig,k)/b1) aa(ig, k) = aa1*dt/(delta(ig,k)*l(ig,k)) qpre = sqrt(q2(ig,k)) IF (aa(ig,k)>0.) THEN q2(ig, k) = (qpre+aa(ig,k)*qpre*qpre)**2 ELSE q2(ig, k) = (qpre/(1.-aa(ig,k)*qpre))**2 END IF ! else ! iflag_pbl=9 ! if (aa(ig,k)*qpre.gt.0.9) then ! q2(ig,k)=(qpre*10.)**2 ! else ! q2(ig,k)=(qpre/(1.-aa(ig,k)*qpre))**2 ! endif ! endif q2(ig, k) = min(max(q2(ig,k),1.E-10), 1.E4) END DO END DO ELSE IF (iflag_pbl>=10) THEN shear(:,:)=0. buoy(:,:)=0. dissip(:,:)=0. km(:,:)=0. IF (yamada4_num>=1) THEN DO k = 2, klev - 1 DO ig=1,ngrid q2(ig, k) = min(max(q2(ig,k),1.E-10), 1.E4) km(ig, k) = l(ig, k)*sqrt(q2(ig,k))*sm(ig, k) shear(ig,k)=km(ig, k)*m2(ig, k) buoy(ig,k)=km(ig, k)*m2(ig, k)*(-1.*rif(ig,k)) ! dissip(ig,k)=min(max(((sqrt(q2(ig,k)))**3)/(b1*l(ig,k)),1.E-12),1.E4) dissip(ig,k)=((sqrt(q2(ig,k)))**3)/(b1*l(ig,k)) ENDDO ENDDO IF (yamada4_num==1) THEN ! Schema du MAR tel quel DO k = 2, klev - 1 DO ig=1,ngrid tkeprov=q2(ig,k)/ydeux tkeprov= tkeprov* & & (tkeprov+dt*(shear(ig,k)+max(0.,buoy(ig,k))))/ & & (tkeprov+dt*((-1.)*min(0.,buoy(ig,k))+dissip(ig,k)+drgpro(ig,k)*tkeprov)) q2(ig,k)=tkeprov*ydeux ENDDO ENDDO ELSE IF (yamada4_num==2) THEN ! version modifiee avec integration exacte pour la dissipation DO k = 2, klev - 1 DO ig=1,ngrid tkeprov=q2(ig,k)/ydeux disseff=dissip(ig,k)-min(0.,buoy(ig,k)) tkeprov = tkeprov/(1.+dt*disseff/(2.*tkeprov))**2 tkeprov= tkeprov+dt*(shear(ig,k)+max(0.,buoy(ig,k))) q2(ig,k)=tkeprov*ydeux ! En cas stable, on traite la flotabilite comme la ! dissipation, en supposant que buoy/q2^3 est constant. ! Puis on prend la solution exacte ENDDO ENDDO ELSE IF (yamada4_num==3) THEN ! version modifiee avec integration exacte pour la dissipation DO k = 2, klev - 1 DO ig=1,ngrid tkeprov=q2(ig,k)/ydeux disseff=dissip(ig,k)-min(0.,buoy(ig,k)) tkeprov=tkeprov*exp(-dt*disseff/tkeprov) tkeprov= tkeprov+dt*(shear(ig,k)+max(0.,buoy(ig,k))) q2(ig,k)=tkeprov*ydeux ! En cas stable, on traite la flotabilite comme la ! dissipation, en supposant que buoy/q2^3 est constant. ! Puis on prend la solution exacte ENDDO ENDDO ELSE IF (yamada4_num==4) THEN ! version modifiee avec integration exacte pour la dissipation DO k = 2, klev - 1 DO ig=1,ngrid tkeprov=q2(ig,k)/ydeux tkeprov= tkeprov+dt*(shear(ig,k)+max(0.,buoy(ig,k))) tkeprov= tkeprov* & & tkeprov/ & & (tkeprov+dt*((-1.)*min(0.,buoy(ig,k))+dissip(ig,k))) q2(ig,k)=tkeprov*ydeux ! En cas stable, on traite la flotabilite comme la ! dissipation, en supposant que buoy/q2^3 est constant. ! Puis on prend la solution exacte ENDDO ENDDO !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! Attention, yamada4_num=5 est inexacte car néglige les termes de flottabilité !! en conditions instables !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ELSE IF (yamada4_num==5) THEN ! version modifiee avec integration exacte pour la dissipation DO k = 2, klev - 1 DO ig=1,ngrid tkeprov=q2(ig,k)/ydeux !FC on ajoute la dissipation due aux arbres disseff=dissip(ig,k)-min(0.,buoy(ig,k)) + drgpro(ig,k)*tkeprov tkeexp=exp(-dt*disseff/tkeprov) ! on prend en compte la tke cree par les arbres winds(ig,k)=sqrt(u(ig,k)**2+v(ig,k)**2) tkeprov= (shear(ig,k)+ & & drgpro(ig,k)*(winds(ig,k))**3)*tkeprov/disseff*(1.-tkeexp)+tkeprov*tkeexp q2(ig,k)=tkeprov*ydeux ! En cas stable, on traite la flotabilite comme la ! dissipation, en supposant que buoy/q2^3 est constant. ! Puis on prend la solution exacte ENDDO ENDDO ELSE IF (yamada4_num==6) THEN ! version modifiee avec integration exacte pour la dissipation DO k = 2, klev - 1 DO ig=1,ngrid ! En cas stable, on traite la flotabilite comme la ! dissipation, en supposant que dissipeff/TKE est constant. ! Puis on prend la solution exacte ! With drag and dissipation from high vegetation (EV & FC, 05/10/2020) winds(ig,k)=sqrt(u(ig,k)**2+v(ig,k)**2) tkeprov=q2(ig,k)/ydeux tkeprov=tkeprov+max(buoy(ig,k)+shear(ig,k)+drgpro(ig,k)*(winds(ig,k))**3,0.)*dt disseff=dissip(ig,k)-min(0.,buoy(ig,k)+shear(ig,k)+drgpro(ig,k)*(winds(ig,k))**3) + drgpro(ig,k)*tkeprov tkeexp=exp(-dt*disseff/tkeprov) tkeprov= tkeprov*tkeexp q2(ig,k)=tkeprov*ydeux ENDDO ENDDO ENDIF DO k = 2, klev - 1 DO ig=1,ngrid q2(ig, k) = min(max(q2(ig,k),1.E-10), 1.E4) ENDDO ENDDO ELSE DO k = 2, klev - 1 km(1:ngrid, k) = l(1:ngrid, k)*sqrt(q2(1:ngrid,k))*sm(1:ngrid, k) q2(1:ngrid, k) = q2(1:ngrid, k) + ydeux*dt*km(1:ngrid, k)*m2(1:ngrid, k)*(1.-rif(1:ngrid,k)) ! q2(1:ngrid, k) = q2(1:ngrid, k) + dt*km(1:ngrid, k)*m2(1:ngrid, k)*(1.-rif(1:ngrid,k)) q2(1:ngrid, k) = min(max(q2(1:ngrid,k),1.E-10), 1.E4) q2(1:ngrid, k) = 1./(1./sqrt(q2(1:ngrid,k))+dt/(yun*l(1:ngrid,k)*b1)) ! q2(1:ngrid, k) = 1./(1./sqrt(q2(1:ngrid,k))+dt/(2*l(1:ngrid,k)*b1)) q2(1:ngrid, k) = q2(1:ngrid, k)*q2(1:ngrid, k) END DO ENDIF ELSE abort_message='Cas nom prevu dans yamada4' CALL abort_physic(modname,abort_message,1) END IF ! Fin du cas 8 ! ==================================================================== ! Calcul des coefficients de melange ! ==================================================================== DO k = 2, klev DO ig = 1, ngrid zq = sqrt(q2(ig,k)) km(ig, k) = l(ig, k)*zq*sm(ig, k) ! For momentum kn(ig, k) = km(ig, k)*alpha(ig, k) ! For scalars kq(ig, k) = l(ig, k)*zq*0.2 ! For TKE END DO END DO !==================================================================== ! Transport diffusif vertical de la TKE par la TKE !==================================================================== ! initialize near-surface and top-layer mixing coefficients !........................................................... kq(1:ngrid, 1) = kq(1:ngrid, 2) ! constant (ie no gradient) near the surface kq(1:ngrid, klev+1) = 0 ! zero at the top ! Transport diffusif vertical de la TKE. !....................................... IF (iflag_vdif_q2==1) THEN q2(1:ngrid, 1) = q2(1:ngrid, 2) CALL vdif_q2(dt, g, rconst, ngrid, plev, temp, kq, q2) END IF !==================================================================== ! Traitement particulier pour les cas tres stables, introduction d'une ! longueur de m??lange minimale !==================================================================== ! ! Reference: Local versus Nonlocal boundary-layer diffusion in a global climate model ! Holtslag A.A.M. and Boville B.A. ! J. Clim., 6, 1825-1842, 1993 IF (hboville) THEN IF (prt_level>1) THEN WRITE(lunout,*) 'YAMADA4 0' END IF DO ig = 1, ngrid coriol(ig) = 1.E-4 pblhmin(ig) = 0.07*ustar(ig)/max(abs(coriol(ig)), 2.546E-5) END DO IF (1==1) THEN IF (iflag_pbl==8 .OR. iflag_pbl==10) THEN DO k = 2, klev DO ig = 1, ngrid IF (teta(ig,2)>teta(ig,1)) THEN qmin = ustar(ig)*(max(1.-zlev(ig,k)/pblhmin(ig),0.))**2 kmin = kap*zlev(ig, k)*qmin ELSE kmin = -1. ! kmin n'est utilise que pour les SL stables. END IF IF (kn(ig,k)<kmin .OR. km(ig,k)<kmin) THEN kn(ig, k) = kmin km(ig, k) = kmin kq(ig, k) = kmin ! la longueur de melange est suposee etre l= kap z ! K=l q Sm d'ou q2=(K/l Sm)**2 q2(ig, k) = (qmin/sm(ig,k))**2 END IF END DO END DO ELSE DO k = 2, klev DO ig = 1, ngrid IF (teta(ig,2)>teta(ig,1)) THEN qmin = ustar(ig)*(max(1.-zlev(ig,k)/pblhmin(ig),0.))**2 kmin = kap*zlev(ig, k)*qmin ELSE kmin = -1. ! kmin n'est utilise que pour les SL stables. END IF IF (kn(ig,k)<kmin .OR. km(ig,k)<kmin) THEN kn(ig, k) = kmin km(ig, k) = kmin kq(ig, k) = kmin ! la longueur de melange est suposee etre l= kap z ! K=l q Sm d'ou q2=(K/l Sm)**2 sm(ig, k) = 1. alpha(ig, k) = 1. q2(ig, k) = min((qmin/sm(ig,k))**2, 10.) zq = sqrt(q2(ig,k)) km(ig, k) = l(ig, k)*zq*sm(ig, k) kn(ig, k) = km(ig, k)*alpha(ig, k) kq(ig, k) = l(ig, k)*zq*0.2 END IF END DO END DO END IF END IF END IF ! hboville ! Ajout d'une viscosite moleculaire km(1:ngrid,2:klev)=km(1:ngrid,2:klev)+viscom kn(1:ngrid,2:klev)=kn(1:ngrid,2:klev)+viscoh kq(1:ngrid,2:klev)=kq(1:ngrid,2:klev)+viscoh IF (prt_level>1) THEN WRITE(lunout,*)'YAMADA4 1' END IF !(prt_level>1) THEN !====================================================== ! Estimations de w'2 et T'2 d'apres Abdela et McFarlane !====================================================== ! ! Reference: A New Second-Order Turbulence Closure Scheme for the Planetary Boundary Layer ! Abdella K and McFarlane N ! J. Atmos. Sci., 54, 1850-1867, 1997 ! Diagnostique pour stokage !.......................... IF (1==0) THEN rino = rif smyam(1:ngrid, 1) = 0. styam(1:ngrid, 1) = 0. lyam(1:ngrid, 1) = 0. knyam(1:ngrid, 1) = 0. w2yam(1:ngrid, 1) = 0. t2yam(1:ngrid, 1) = 0. smyam(1:ngrid, 2:klev) = sm(1:ngrid, 2:klev) styam(1:ngrid, 2:klev) = sm(1:ngrid, 2:klev)*alpha(1:ngrid, 2:klev) lyam(1:ngrid, 2:klev) = l(1:ngrid, 2:klev) knyam(1:ngrid, 2:klev) = kn(1:ngrid, 2:klev) ! Calcul de w'2 et T'2 !....................... w2yam(1:ngrid, 2:klev) = q2(1:ngrid, 2:klev)*0.24 + & lyam(1:ngrid, 2:klev)*5.17*kn(1:ngrid, 2:klev)*n2(1:ngrid, 2:klev)/ & sqrt(q2(1:ngrid,2:klev)) t2yam(1:ngrid, 2:klev) = 9.1*kn(1:ngrid, 2:klev)* & dtetadz(1:ngrid, 2:klev)**2/sqrt(q2(1:ngrid,2:klev))* & lyam(1:ngrid, 2:klev) END IF !============================================================================ ! Mise a jour de la tke !============================================================================ IF (new_yamada4) THEN DO k=1,klev+1 tke(1:ngrid,k)=q2(1:ngrid,k)/ydeux ENDDO ELSE DO k=1,klev+1 tke(1:ngrid,k)=q2(1:ngrid,k) ENDDO ENDIF !============================================================================ ! Diagnostique de la dissipation et vitesse verticale !============================================================================ ! Diagnostics eps(:,:)=0. wprime(1:ngrid,:,nsrf)=0. DO k=2,klev DO ig=1,ngrid eps(ig,k)=dissip(ig,k) jg=ni(ig) wprime(jg,k,nsrf)=sqrt(MAX(1./3*q2(ig,k),0.)) ENDDO ENDDO !============================================================================= RETURN END SUBROUTINE yamada4 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE vdif_q2(timestep, gravity, rconst, ngrid, plev, temp, kmy, q2) USE dimphy, only : klev,klon IMPLICIT NONE ! vdif_q2: subroutine qui calcule la diffusion de la TKE par la TKE ! avec un schema implicite en temps avec ! inversion d'un syst??me tridiagonal ! ! Reference: Description of the interface with the surface and ! the computation of the turbulet diffusion in LMDZ ! Technical note on LMDZ ! Dufresne, J-L, Ghattas, J. and Grandpeix, J-Y ! !============================================================================ ! Declarations !============================================================================ REAL plev(klon, klev+1) REAL temp(klon, klev) REAL timestep REAL gravity, rconst REAL kstar(klon, klev+1), zz REAL kmy(klon, klev+1) REAL q2(klon, klev+1) REAL deltap(klon, klev+1) REAL denom(klon, klev+1), alpha(klon, klev+1), beta(klon, klev+1) INTEGER ngrid INTEGER i, k !========================================================================= ! Calcul !========================================================================= DO k = 1, klev DO i = 1, ngrid zz = (plev(i,k)+plev(i,k+1))*gravity/(rconst*temp(i,k)) kstar(i, k) = 0.125*(kmy(i,k+1)+kmy(i,k))*zz*zz/ & (plev(i,k)-plev(i,k+1))*timestep END DO END DO DO k = 2, klev DO i = 1, ngrid deltap(i, k) = 0.5*(plev(i,k-1)-plev(i,k+1)) END DO END DO DO i = 1, ngrid deltap(i, 1) = 0.5*(plev(i,1)-plev(i,2)) deltap(i, klev+1) = 0.5*(plev(i,klev)-plev(i,klev+1)) denom(i, klev+1) = deltap(i, klev+1) + kstar(i, klev) alpha(i, klev+1) = deltap(i, klev+1)*q2(i, klev+1)/denom(i, klev+1) beta(i, klev+1) = kstar(i, klev)/denom(i, klev+1) END DO DO k = klev, 2, -1 DO i = 1, ngrid denom(i, k) = deltap(i, k) + (1.-beta(i,k+1))*kstar(i, k) + & kstar(i, k-1) alpha(i, k) = (q2(i,k)*deltap(i,k)+kstar(i,k)*alpha(i,k+1))/denom(i, k) beta(i, k) = kstar(i, k-1)/denom(i, k) END DO END DO ! Si on recalcule q2(1) !....................... IF (1==0) THEN DO i = 1, ngrid denom(i, 1) = deltap(i, 1) + (1-beta(i,2))*kstar(i, 1) q2(i, 1) = (q2(i,1)*deltap(i,1)+kstar(i,1)*alpha(i,2))/denom(i, 1) END DO END IF DO k = 2, klev + 1 DO i = 1, ngrid q2(i, k) = alpha(i, k) + beta(i, k)*q2(i, k-1) END DO END DO RETURN END SUBROUTINE vdif_q2 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE vdif_q2e(timestep, gravity, rconst, ngrid, plev, temp, kmy, q2) USE dimphy, only : klev,klon IMPLICIT NONE ! vdif_q2e: subroutine qui calcule la diffusion de TKE par la TKE ! avec un schema explicite en temps !==================================================== ! Declarations !==================================================== REAL plev(klon, klev+1) REAL temp(klon, klev) REAL timestep REAL gravity, rconst REAL kstar(klon, klev+1), zz REAL kmy(klon, klev+1) REAL q2(klon, klev+1) REAL deltap(klon, klev+1) REAL denom(klon, klev+1), alpha(klon, klev+1), beta(klon, klev+1) INTEGER ngrid INTEGER i, k !================================================== ! Calcul !================================================== DO k = 1, klev DO i = 1, ngrid zz = (plev(i,k)+plev(i,k+1))*gravity/(rconst*temp(i,k)) kstar(i, k) = 0.125*(kmy(i,k+1)+kmy(i,k))*zz*zz/ & (plev(i,k)-plev(i,k+1))*timestep END DO END DO DO k = 2, klev DO i = 1, ngrid deltap(i, k) = 0.5*(plev(i,k-1)-plev(i,k+1)) END DO END DO DO i = 1, ngrid deltap(i, 1) = 0.5*(plev(i,1)-plev(i,2)) deltap(i, klev+1) = 0.5*(plev(i,klev)-plev(i,klev+1)) END DO DO k = klev, 2, -1 DO i = 1, ngrid q2(i, k) = q2(i, k) + (kstar(i,k)*(q2(i,k+1)-q2(i, & k))-kstar(i,k-1)*(q2(i,k)-q2(i,k-1)))/deltap(i, k) END DO END DO DO i = 1, ngrid q2(i, 1) = q2(i, 1) + (kstar(i,1)*(q2(i,2)-q2(i,1)))/deltap(i, 1) q2(i, klev+1) = q2(i, klev+1) + (-kstar(i,klev)*(q2(i,klev+1)-q2(i, & klev)))/deltap(i, klev+1) END DO RETURN END SUBROUTINE vdif_q2e !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE mixinglength(ni, nsrf, ngrid,iflag_pbl,pbl_lmixmin_alpha,lmixmin,zlay,zlev,u,v,q2,n2, lmix) USE dimphy, only : klev,klon USE yamada_ini_mod, only : l0 USE phys_state_var_mod, only: zstd, zsig, zmea USE phys_local_var_mod, only: l_mixmin, l_mix USE yamada_ini_mod, only : kap, kapb ! zstd: ecart type de la'altitud e sous-maille ! zmea: altitude moyenne sous maille ! zsig: pente moyenne de le maille USE geometry_mod, only: cell_area ! aire_cell: aire de la maille IMPLICIT NONE !************************************************************************* ! Subrourine qui calcule la longueur de m??lange dans le sch??ma de turbulence ! avec la formule de Blackadar ! Calcul d'un minimum en fonction de l'orographie sous-maille: ! L'id??e est la suivante: plus il y a de relief, plus il y a du m??lange ! induit par les circulations meso et submeso ??chelles. ! ! References: * The vertical distribution of wind and turbulent exchange in a neutral atmosphere ! Blackadar A.K. ! J. Geophys. Res., 64, No 8, 1962 ! ! * An evaluation of neutral and convective planetary boundary-layer parametrisations relative ! to large eddy simulations ! Ayotte K et al ! Boundary Layer Meteorology, 79, 131-175, 1996 ! ! ! * Local Similarity in the Stable Boundary Layer and Mixing length Approaches: consistency of concepts ! Van de Wiel B.J.H et al ! Boundary-Lay Meteorol, 128, 103-166, 2008 ! ! ! Histoire: !---------- ! * premi??re r??daction, Etienne et Frederic, 09/06/2016 ! ! *********************************************************************** !================================================================== ! Declarations !================================================================== ! Inputs !------- INTEGER ni(klon) ! indice sur la grille original (non restreinte) INTEGER nsrf ! Type de surface INTEGER ngrid ! Nombre de points concern??s sur l'horizontal INTEGER iflag_pbl ! Choix du sch??ma de turbulence REAL pbl_lmixmin_alpha ! on active ou non le calcul de la longueur de melange minimum en fonction du relief REAL lmixmin ! Minimum absolu de la longueur de m??lange REAL zlay(klon, klev) ! altitude du centre de la couche REAL zlev(klon, klev+1) ! atitude de l'interface inf??rieure de la couche REAL u(klon, klev) ! vitesse du vent zonal REAL v(klon, klev) ! vitesse du vent meridional REAL q2(klon, klev+1) ! energie cin??tique turbulente REAL n2(klon, klev+1) ! frequence de Brunt-Vaisala !In/out !------- ! Outputs !--------- REAL lmix(klon, klev+1) ! Longueur de melange ! Local !------- INTEGER ig,jg, k REAL h_oro(klon) REAL hlim(klon) REAL zq REAL sq(klon), sqz(klon) REAL fl, zzz, zl0, zq2, zn2 REAL famorti, zzzz, zh_oro, zhlim REAL l1(klon, klev+1), l2(klon,klev+1) REAL winds(klon, klev) REAL xcell REAL zstdslope(klon) REAL lmax REAL l2strat, l2neutre, extent REAL l2limit(klon) !=============================================================== ! Fonctions utiles !=============================================================== ! Calcul de l suivant la formule de Blackadar 1962 adapt??e par Ayotte 1996 !.......................................................................... fl(zzz, zl0, zq2, zn2) = max(min(l0(ig)*kap*zlev(ig, & k)/(kap*zlev(ig,k)+l0(ig)),0.5*sqrt(q2(ig,k))/sqrt( & max(n2(ig,k),1.E-10))), 1.E-5) ! Fonction d'amortissement de la turbulence au dessus de la montagne ! On augmente l'amortissement en diminuant la valeur de hlim (extent) dans le code !..................................................................... famorti(zzzz, zh_oro, zhlim)=(-1.)*ATAN((zzzz-zh_oro)/(zhlim-zh_oro))*2./3.1416+1. IF (ngrid==0) RETURN !===================================================================== ! CALCUL de la LONGUEUR de m??lange suivant BLACKADAR: l1 !===================================================================== l1(1:ngrid,:)=0. IF (iflag_pbl==8 .OR. iflag_pbl==10) THEN ! Iterative computation of l0 ! This version is kept for iflag_pbl only for convergence ! with NPv3.1 Cmip5 simulations !................................................................... DO ig = 1, ngrid sq(ig) = 1.E-10 sqz(ig) = 1.E-10 END DO DO k = 2, klev - 1 DO ig = 1, ngrid zq = sqrt(q2(ig,k)) sqz(ig) = sqz(ig) + zq*zlev(ig, k)*(zlay(ig,k)-zlay(ig,k-1)) sq(ig) = sq(ig) + zq*(zlay(ig,k)-zlay(ig,k-1)) END DO END DO DO ig = 1, ngrid l0(ig) = 0.2*sqz(ig)/sq(ig) END DO DO k = 2, klev DO ig = 1, ngrid l1(ig, k) = fl(zlev(ig,k), l0(ig), q2(ig,k), n2(ig,k)) END DO END DO ELSE ! In all other case, the assymptotic mixing length l0 is imposed (150m) !...................................................................... l0(1:ngrid) = 150. DO k = 2, klev DO ig = 1, ngrid l1(ig, k) = fl(zlev(ig,k), l0(ig), q2(ig,k), n2(ig,k)) END DO END DO END IF !=========================================================================================== ! CALCUL d'une longueur de melange minimum en fonctions de la topographie sous maille: l2 ! si pbl_lmixmin_alpha=TRUE et si on se trouve sur de la terre ( pas actif sur les ! glacier, la glace de mer et les oc??ans) !=========================================================================================== l2(1:ngrid,:)=0.0 l_mixmin(1:ngrid,:,nsrf)=0. l_mix(1:ngrid,:,nsrf)=0. hlim(1:ngrid)=0. IF (nsrf .EQ. 1) THEN ! coefficients !-------------- extent=2. ! On ??tend l'impact du relief jusqu'?? extent*h, extent >1. lmax=150. ! Longueur de m??lange max dans l'absolu ! calculs !--------- DO ig=1,ngrid ! On calcule la hauteur du relief !................................. ! On ne peut pas prendre zstd seulement pour caracteriser le relief sous maille ! car sur un terrain pentu mais sans relief, zstd est non nul (comme en Antarctique, C. Genthon) ! On corrige donc zstd avec l'ecart type de l'altitude dans la maille sans relief ! (en gros, une maille de taille xcell avec une pente constante zstdslope) jg=ni(ig) ! IF (zsig(jg) .EQ. 0.) THEN ! zstdslope(ig)=0. ! ELSE ! xcell=sqrt(cell_area(jg)) ! zstdslope(ig)=max((xcell*zsig(jg)-zmea(jg))**3 /(3.*zsig(jg)),0.) ! zstdslope(ig)=sqrt(zstdslope(ig)) ! END IF ! h_oro(ig)=max(zstd(jg)-zstdslope(ig),0.) ! Hauteur du relief h_oro(ig)=zstd(jg) hlim(ig)=extent*h_oro(ig) ENDDO l2limit(1:ngrid)=0. DO k=2,klev DO ig=1,ngrid winds(ig,k)=sqrt(u(ig,k)**2+v(ig,k)**2) IF (zlev(ig,k) .LE. h_oro(ig)) THEN ! sous l'orographie l2strat= kapb*pbl_lmixmin_alpha*winds(ig,k)/sqrt(max(n2(ig,k),1.E-10)) ! si stratifi??, amplitude d'oscillation * kappab (voir Van de Wiel et al 2008) l2neutre=kap*zlev(ig,k)*h_oro(ig)/(kap*zlev(ig,k)+h_oro(ig)) ! Dans le cas neutre, formule de blackadar. tend asymptotiquement vers h l2neutre=MIN(l2neutre,lmax) ! On majore par lmax l2limit(ig)=MIN(l2neutre,l2strat) ! Calcule de l2 (minimum de la longueur en cas neutre et celle en situation stratifi??e) l2(ig,k)=l2limit(ig) ELSE IF (zlev(ig,k) .LE. hlim(ig)) THEN ! Si on est au dessus des montagnes, mais affect?? encore par elles ! Au dessus des montagnes, on prend la l2limit au sommet des montagnes ! (la derni??re calcul??e dans la boucle k, vu que k est un indice croissant avec z) ! et on multiplie l2limit par une fonction qui d??croit entre h et hlim l2(ig,k)=l2limit(ig)*famorti(zlev(ig,k),h_oro(ig), hlim(ig)) ELSE ! Au dessus de extent*h, on prend l2=l0 l2(ig,k)=0. END IF ENDDO ENDDO ENDIF ! pbl_lmixmin_alpha !================================================================================== ! On prend le max entre la longueur de melange de blackadar et celle calcul??e ! en fonction de la topographie !=================================================================================== DO k=1,klev+1 DO ig=1,ngrid lmix(ig,k)=MAX(MAX(l1(ig,k), l2(ig,k)),lmixmin) ENDDO ENDDO ! Diagnostics DO k=1,klev+1 DO ig=1,ngrid jg=ni(ig) l_mix(jg,k,nsrf)=lmix(ig,k) l_mixmin(jg,k,nsrf)=MAX(l2(ig,k),lmixmin) ENDDO ENDDO END SUBROUTINE mixinglength