MODULE coef_diff_turb_mod ! This module contains some procedures for calculation of the coefficients of the ! turbulent diffusion in the atmosphere and coefficients for turbulent diffusion ! at surface(cdrag) IMPLICIT NONE CONTAINS !**************************************************************************************** SUBROUTINE coef_diff_turb(dtime, nsrf, knon, ni, & ypaprs, ypplay, yu, yv, yq, yt, yts, yqsurf, ycdragm, & ycoefm, ycoefh ,yq2, yeps, ydrgpro) USE dimphy USE indice_sol_mod USE lmdz_print_control, ONLY: prt_level, lunout ! Calculate coefficients(ycoefm, ycoefh) for turbulent diffusion in the ! atmosphere ! NB! No values are calculated between surface and the first model layer. ! ycoefm(:,1) and ycoefh(:,1) are not valid !!! ! Input arguments !**************************************************************************************** REAL, INTENT(IN) :: dtime INTEGER, INTENT(IN) :: nsrf, knon INTEGER, DIMENSION(klon), INTENT(IN) :: ni REAL, DIMENSION(klon,klev+1), INTENT(IN) :: ypaprs REAL, DIMENSION(klon,klev), INTENT(IN) :: ypplay REAL, DIMENSION(klon,klev), INTENT(IN) :: yu, yv REAL, DIMENSION(klon,klev), INTENT(IN) :: yq, yt REAL, DIMENSION(klon), INTENT(IN) :: yts, yqsurf REAL, DIMENSION(klon), INTENT(IN) :: ycdragm !FC REAL, DIMENSION(klon,klev), INTENT(IN) :: ydrgpro ! InOutput arguments !**************************************************************************************** REAL, DIMENSION(klon,klev+1), INTENT(INOUT):: yq2 ! Output arguments !**************************************************************************************** REAL, DIMENSION(klon,klev+1), INTENT(OUT) :: yeps REAL, DIMENSION(klon,klev), INTENT(OUT) :: ycoefh REAL, DIMENSION(klon,klev), INTENT(OUT) :: ycoefm ! Other local variables !**************************************************************************************** INTEGER :: k, i, j REAL, DIMENSION(klon,klev) :: ycoefm0, ycoefh0, yzlay, yteta REAL, DIMENSION(klon,klev+1) :: yzlev, q2diag, ykmm, ykmn, ykmq REAL, DIMENSION(klon) :: yustar ! Include !**************************************************************************************** INCLUDE "clesphys.h" INCLUDE "compbl.h" INCLUDE "YOETHF.h" INCLUDE "YOMCST.h" ykmm = 0 !ym missing init ykmn = 0 !ym missing init ykmq = 0 !ym missing init yeps(:,:) = 0. !**************************************************************************************** ! Calcul de coefficients de diffusion turbulent de l'atmosphere : ! ycoefm(:,2:klev), ycoefh(:,2:klev) !**************************************************************************************** CALL coefkz(nsrf, knon, ypaprs, ypplay, & ksta, ksta_ter, & yts, yu, yv, yt, yq, & yqsurf, & ycoefm, ycoefh) !**************************************************************************************** ! Eventuelle recalcule des coeffeicients de diffusion turbulent de l'atmosphere : ! ycoefm(:,2:klev), ycoefh(:,2:klev) !**************************************************************************************** IF (iflag_pbl==1) THEN CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, & ycoefm0, ycoefh0) DO k = 2, klev DO i = 1, knon ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k)) ycoefh(i,k) = MAX(ycoefh(i,k),ycoefh0(i,k)) ENDDO ENDDO ENDIF !**************************************************************************************** ! Calcul d'une diffusion minimale pour les conditions tres stables !**************************************************************************************** IF (ok_kzmin) THEN CALL coefkzmin(knon,ypaprs,ypplay,yu,yv,yt,yq,ycdragm, & ycoefm0,ycoefh0) DO k = 2, klev DO i = 1, knon ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k)) ycoefh(i,k) = MAX(ycoefh(i,k),ycoefh0(i,k)) ENDDO ENDDO ENDIF !**************************************************************************************** ! MELLOR ET YAMADA adapte a Mars Richard Fournier et Frederic Hourdin !**************************************************************************************** IF (iflag_pbl>=3) THEN yzlay(1:knon,1)= & RD*yt(1:knon,1)/(0.5*(ypaprs(1:knon,1)+ypplay(1:knon,1))) & *(ypaprs(1:knon,1)-ypplay(1:knon,1))/RG DO k=2,klev DO i = 1, knon yzlay(i,k)= & yzlay(i,k-1)+RD*0.5*(yt(i,k-1)+yt(i,k)) & /ypaprs(i,k)*(ypplay(i,k-1)-ypplay(i,k))/RG END DO END DO DO k=1,klev DO i = 1, knon yteta(i,k)= & yt(i,k)*(ypaprs(i,1)/ypplay(i,k))**RKAPPA & *(1.+0.61*yq(i,k)) END DO END DO yzlev(1:knon,1)=0. yzlev(1:knon,klev+1)=2.*yzlay(1:knon,klev)-yzlay(1:knon,klev-1) DO k=2,klev DO i = 1, knon yzlev(i,k)=0.5*(yzlay(i,k)+yzlay(i,k-1)) END DO END DO !!$!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!$! Pour memoire, le papier Hourdin et al. 2002 a ete obtenur avec un !!$! bug sur les coefficients de surface : !!$! ycdragh(1:knon) = ycoefm(1:knon,1) !!$! ycdragm(1:knon) = ycoefh(1:knon,1) !!$!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Normalement, on peut passer dans les codes avec knon=0 ! Mais ca fait planter le replay. ! En attendant une réécriture, on a joute des if (Fredho) IF ( klon>1 .OR. (klon==1 .AND. knon==1) ) THEN CALL ustarhb(knon,klev,knon,yu,yv,ycdragm, yustar) endif IF (prt_level > 9) THEN WRITE(lunout,*) 'USTAR = ',(yustar(i),i=1,knon) ENDIF ! iflag_pbl peut etre utilise comme longuer de melange IF (iflag_pbl>=31) THEN IF ( klon>1 .OR. (klon==1 .AND. knon==1) ) THEN CALL vdif_kcay(knon,klev,knon,dtime,RG,RD,ypaprs,yt, & yzlev,yzlay,yu,yv,yteta, & ycdragm,yq2,q2diag,ykmm,ykmn,yustar, & iflag_pbl) endif ELSE IF (iflag_pbl<20) THEN CALL yamada4(ni,nsrf,knon,dtime,RG,RD,ypaprs,yt, & yzlev,yzlay,yu,yv,yteta, & ycdragm,yq2,yeps,ykmm,ykmn,ykmq,yustar, & iflag_pbl,ydrgpro) !FC ENDIF ycoefm(1:knon,2:klev)=ykmm(1:knon,2:klev) ycoefh(1:knon,2:klev)=ykmn(1:knon,2:klev) ELSE ! No TKE for Standard Physics yq2=0. ENDIF !(iflag_pbl.ge.3) END SUBROUTINE coef_diff_turb !**************************************************************************************** SUBROUTINE coefkz(nsrf, knon, paprs, pplay, & ksta, ksta_ter, & ts, & u,v,t,q, & qsurf, & pcfm, pcfh) USE dimphy USE indice_sol_mod USE lmdz_print_control, ONLY: prt_level, lunout !====================================================================== ! Auteur(s) F. Hourdin, M. Forichon, Z.X. Li (LMD/CNRS) date: 19930922 ! (une version strictement identique a l'ancien modele) ! Objet: calculer le coefficient du frottement du sol (Cdrag) et les ! coefficients d'echange turbulent dans l'atmosphere. ! Arguments: ! nsrf-----input-I- indicateur de la nature du sol ! knon-----input-I- nombre de points a traiter ! paprs----input-R- pregssion a chaque intercouche (en Pa) ! pplay----input-R- pression au milieu de chaque couche (en Pa) ! ts-------input-R- temperature du sol (en Kelvin) ! u--------input-R- vitesse u ! v--------input-R- vitesse v ! t--------input-R- temperature (K) ! q--------input-R- vapeur d'eau (kg/kg) ! pcfm-----output-R- coefficients a calculer (vitesse) ! pcfh-----output-R- coefficients a calculer (chaleur et humidite) !====================================================================== INCLUDE "YOETHF.h" INCLUDE "YOMCST.h" INCLUDE "FCTTRE.h" INCLUDE "compbl.h" ! Arguments: INTEGER, INTENT(IN) :: knon, nsrf REAL, INTENT(IN) :: ksta, ksta_ter REAL, DIMENSION(klon), INTENT(IN) :: ts REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay REAL, DIMENSION(klon,klev), INTENT(IN) :: u, v, t, q REAL, DIMENSION(klon), INTENT(IN) :: qsurf REAL, DIMENSION(klon,klev), INTENT(OUT) :: pcfm, pcfh ! Local variables: INTEGER, DIMENSION(klon) :: itop ! numero de couche du sommet de la couche limite ! Quelques constantes et options: REAL, PARAMETER :: cepdu2=0.1**2 REAL, PARAMETER :: CKAP=0.4 REAL, PARAMETER :: cb=5.0 REAL, PARAMETER :: cc=5.0 REAL, PARAMETER :: cd=5.0 REAL, PARAMETER :: clam=160.0 REAL, PARAMETER :: ratqs=0.05 ! largeur de distribution de vapeur d'eau LOGICAL, PARAMETER :: richum=.TRUE. ! utilise le nombre de Richardson humide REAL, PARAMETER :: ric=0.4 ! nombre de Richardson critique REAL, PARAMETER :: prandtl=0.4 REAL kstable ! diffusion minimale (situation stable) ! GKtest ! PARAMETER (kstable=1.0e-10) !IM: 261103 REAL kstable_ter, kstable_sinon !IM: 211003 cf GK PARAMETER (kstable_ter = 1.0e-6) !IM: 261103 PARAMETER (kstable_ter = 1.0e-8) !IM: 261103 PARAMETER (kstable_ter = 1.0e-10) !IM: 261103 PARAMETER (kstable_sinon = 1.0e-10) ! fin GKtest REAL, PARAMETER :: mixlen=35.0 ! constante controlant longueur de melange INTEGER isommet ! le sommet de la couche limite LOGICAL, PARAMETER :: tvirtu=.TRUE. ! calculer Ri d'une maniere plus performante LOGICAL, PARAMETER :: opt_ec=.FALSE.! formule du Centre Europeen dans l'atmosphere ! Variables locales: INTEGER i, k !IM 120704 REAL zgeop(klon,klev) REAL zmgeom(klon) REAL zri(klon) REAL zl2(klon) REAL zdphi, zdu2, ztvd, ztvu, zcdn REAL zscf REAL zt, zq, zdelta, zcvm5, zcor, zqs, zfr, zdqs REAL z2geomf, zalh2, zalm2, zscfh, zscfm REAL, PARAMETER :: t_coup=273.15 LOGICAL, PARAMETER :: check=.FALSE. ! contre-gradient pour la chaleur sensible: Kelvin/metre REAL gamt(2:klev) LOGICAL, SAVE :: appel1er=.TRUE. !$OMP THREADPRIVATE(appel1er) ! Fonctions thermodynamiques et fonctions d'instabilite REAL fsta, fins, x fsta(x) = 1.0 / (1.0+10.0*x*(1+8.0*x)) fins(x) = SQRT(1.0-18.0*x) isommet=klev IF (appel1er) THEN IF (prt_level > 9) THEN WRITE(lunout,*)'coefkz, opt_ec:', opt_ec WRITE(lunout,*)'coefkz, richum:', richum IF (richum) WRITE(lunout,*)'coefkz, ratqs:', ratqs WRITE(lunout,*)'coefkz, isommet:', isommet WRITE(lunout,*)'coefkz, tvirtu:', tvirtu appel1er = .FALSE. ENDIF ENDIF ! Initialiser les sorties DO k = 1, klev DO i = 1, knon pcfm(i,k) = 0.0 pcfh(i,k) = 0.0 ENDDO ENDDO DO i = 1, knon itop(i) = 0 ENDDO ! Prescrire la valeur de contre-gradient IF (iflag_pbl==1) THEN DO k = 3, klev gamt(k) = -1.0E-03 ENDDO gamt(2) = -2.5E-03 ELSE DO k = 2, klev gamt(k) = 0.0 ENDDO ENDIF !IM cf JLD/ GKtest IF ( nsrf /= is_oce ) THEN !IM 261103 kstable = kstable_ter kstable = ksta_ter ELSE !IM 261103 kstable = kstable_sinon kstable = ksta ENDIF !IM cf JLD/ GKtest fin ! Calculer les geopotentiels de chaque couche DO i = 1, knon zgeop(i,1) = RD * t(i,1) / (0.5*(paprs(i,1)+pplay(i,1))) & * (paprs(i,1)-pplay(i,1)) ENDDO DO k = 2, klev DO i = 1, knon zgeop(i,k) = zgeop(i,k-1) & + RD * 0.5*(t(i,k-1)+t(i,k)) / paprs(i,k) & * (pplay(i,k-1)-pplay(i,k)) ENDDO ENDDO ! Calculer les coefficients turbulents dans l'atmosphere DO i = 1, knon itop(i) = isommet ENDDO DO k = 2, isommet DO i = 1, knon zdu2=MAX(cepdu2,(u(i,k)-u(i,k-1))**2 & +(v(i,k)-v(i,k-1))**2) zmgeom(i)=zgeop(i,k)-zgeop(i,k-1) zdphi =zmgeom(i) / 2.0 zt = (t(i,k)+t(i,k-1)) * 0.5 zq = (q(i,k)+q(i,k-1)) * 0.5 ! Calculer Qs et dQs/dT: IF (thermcep) THEN zdelta = MAX(0.,SIGN(1.,RTT-zt)) zcvm5 = R5LES*RLVTT/RCPD/(1.0+RVTMP2*zq)*(1.-zdelta) & + R5IES*RLSTT/RCPD/(1.0+RVTMP2*zq)*zdelta zqs = R2ES * FOEEW(zt,zdelta) / pplay(i,k) zqs = MIN(0.5,zqs) zcor = 1./(1.-RETV*zqs) zqs = zqs*zcor zdqs = FOEDE(zt,zdelta,zcvm5,zqs,zcor) ELSE IF (zt < t_coup) THEN zqs = qsats(zt) / pplay(i,k) zdqs = dqsats(zt,zqs) ELSE zqs = qsatl(zt) / pplay(i,k) zdqs = dqsatl(zt,zqs) ENDIF ENDIF ! calculer la fraction nuageuse (processus humide): IF (zq /= 0.) THEN zfr = (zq+ratqs*zq-zqs) / (2.0*ratqs*zq) else zfr = 0. end if zfr = MAX(0.0,MIN(1.0,zfr)) IF (.NOT.richum) zfr = 0.0 ! calculer le nombre de Richardson: IF (tvirtu) THEN ztvd =( t(i,k) & + zdphi/RCPD/(1.+RVTMP2*zq) & *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) & )*(1.+RETV*q(i,k)) ztvu =( t(i,k-1) & - zdphi/RCPD/(1.+RVTMP2*zq) & *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) & )*(1.+RETV*q(i,k-1)) zri(i) =zmgeom(i)*(ztvd-ztvu)/(zdu2*0.5*(ztvd+ztvu)) zri(i) = zri(i) & + zmgeom(i)*zmgeom(i)/RG*gamt(k) & *(paprs(i,k)/101325.0)**RKAPPA & /(zdu2*0.5*(ztvd+ztvu)) ELSE ! calcul de Ridchardson compatible LMD5 zri(i) =(RCPD*(t(i,k)-t(i,k-1)) & -RD*0.5*(t(i,k)+t(i,k-1))/paprs(i,k) & *(pplay(i,k)-pplay(i,k-1)) & )*zmgeom(i)/(zdu2*0.5*RCPD*(t(i,k-1)+t(i,k))) zri(i) = zri(i) + & zmgeom(i)*zmgeom(i)*gamt(k)/RG & *(paprs(i,k)/101325.0)**RKAPPA & /(zdu2*0.5*(t(i,k-1)+t(i,k))) ENDIF ! finalement, les coefficients d'echange sont obtenus: zcdn=SQRT(zdu2) / zmgeom(i) * RG IF (opt_ec) THEN z2geomf=zgeop(i,k-1)+zgeop(i,k) zalm2=(0.5*ckap/RG*z2geomf & /(1.+0.5*ckap/rg/clam*z2geomf))**2 zalh2=(0.5*ckap/rg*z2geomf & /(1.+0.5*ckap/RG/(clam*SQRT(1.5*cd))*z2geomf))**2 IF (zri(i)<0.0) THEN ! situation instable zscf = ((zgeop(i,k)/zgeop(i,k-1))**(1./3.)-1.)**3 & / (zmgeom(i)/RG)**3 / (zgeop(i,k-1)/RG) zscf = SQRT(-zri(i)*zscf) zscfm = 1.0 / (1.0+3.0*cb*cc*zalm2*zscf) zscfh = 1.0 / (1.0+3.0*cb*cc*zalh2*zscf) pcfm(i,k)=zcdn*zalm2*(1.-2.0*cb*zri(i)*zscfm) pcfh(i,k)=zcdn*zalh2*(1.-3.0*cb*zri(i)*zscfh) ELSE ! situation stable zscf=SQRT(1.+cd*zri(i)) pcfm(i,k)=zcdn*zalm2/(1.+2.0*cb*zri(i)/zscf) pcfh(i,k)=zcdn*zalh2/(1.+3.0*cb*zri(i)*zscf) ENDIF ELSE zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,itop(i)+1)) & /(paprs(i,2)-paprs(i,itop(i)+1)) ))**2 pcfm(i,k)=SQRT(MAX(zcdn*zcdn*(ric-zri(i))/ric, kstable)) pcfm(i,k)= zl2(i)* pcfm(i,k) pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different ENDIF ENDDO ENDDO ! Au-dela du sommet, pas de diffusion turbulente: DO i = 1, knon IF (itop(i)+1 <= klev) THEN DO k = itop(i)+1, klev pcfh(i,k) = 0.0 pcfm(i,k) = 0.0 ENDDO ENDIF ENDDO END SUBROUTINE coefkz !**************************************************************************************** SUBROUTINE coefkz2(nsrf, knon, paprs, pplay,t, & pcfm, pcfh) USE dimphy USE indice_sol_mod !====================================================================== ! J'introduit un peu de diffusion sauf dans les endroits ! ou une forte inversion est presente ! On peut dire qu'il represente la convection peu profonde ! Arguments: ! nsrf-----input-I- indicateur de la nature du sol ! knon-----input-I- nombre de points a traiter ! paprs----input-R- pression a chaque intercouche (en Pa) ! pplay----input-R- pression au milieu de chaque couche (en Pa) ! t--------input-R- temperature (K) ! pcfm-----output-R- coefficients a calculer (vitesse) ! pcfh-----output-R- coefficients a calculer (chaleur et humidite) !====================================================================== ! Arguments: INTEGER, INTENT(IN) :: knon, nsrf REAL, DIMENSION(klon, klev+1), INTENT(IN) :: paprs REAL, DIMENSION(klon, klev), INTENT(IN) :: pplay REAL, DIMENSION(klon, klev), INTENT(IN) :: t(klon,klev) REAL, DIMENSION(klon, klev), INTENT(OUT) :: pcfm, pcfh ! Quelques constantes et options: REAL, PARAMETER :: prandtl=0.4 REAL, PARAMETER :: kstable=0.002 ! REAL, PARAMETER :: kstable=0.001 REAL, PARAMETER :: mixlen=35.0 ! constante controlant longueur de melange REAL, PARAMETER :: seuil=-0.02 ! au-dela l'inversion est consideree trop faible ! PARAMETER (seuil=-0.04) ! PARAMETER (seuil=-0.06) ! PARAMETER (seuil=-0.09) ! Variables locales: INTEGER i, k, invb(knon) REAL zl2(knon) REAL zdthmin(knon), zdthdp INCLUDE "YOMCST.h" ! Initialiser les sorties DO k = 1, klev DO i = 1, knon pcfm(i,k) = 0.0 pcfh(i,k) = 0.0 ENDDO ENDDO ! Chercher la zone d'inversion forte DO i = 1, knon invb(i) = klev zdthmin(i)=0.0 ENDDO DO k = 2, klev/2-1 DO i = 1, knon zdthdp = (t(i,k)-t(i,k+1))/(pplay(i,k)-pplay(i,k+1)) & - RD * 0.5*(t(i,k)+t(i,k+1))/RCPD/paprs(i,k+1) zdthdp = zdthdp * 100.0 IF (pplay(i,k)>0.8*paprs(i,1) .AND. & zdthdpseuil) ) THEN zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,klev+1)) & /(paprs(i,2)-paprs(i,klev+1)) ))**2 pcfm(i,k)= zl2(i)* kstable pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different ENDIF ENDDO ENDDO ENDIF END SUBROUTINE coefkz2 !**************************************************************************************** END MODULE coef_diff_turb_mod