! ! $Header$ ! 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) ! USE dimphy IMPLICIT NONE CONTAINS ! !**************************************************************************************** ! SUBROUTINE coef_diff_turb(dtime, nsrf, knon, ni, & ypaprs, ypplay, yu, yv, yq, yt, yts, yrugos, yqsurf, & ycoefm, ycoefh ,yq2) ! ! Calculate coefficients(ycoefm, ycoefh) for turbulent diffusion in the ! atmosphere and coefficients for turbulent diffusion at surface - cdrag ! (ycoefm(:,1), ycoefh(:,1)). ! ! ! 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, yrugos, yqsurf REAL, DIMENSION(klon,klev+1) :: yq2 ! Output arguments !**************************************************************************************** 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) :: y_cd_h, y_cd_m REAL, DIMENSION(klon) :: yustar ! Include !**************************************************************************************** INCLUDE "clesphys.h" INCLUDE "indicesol.h" INCLUDE "iniprint.h" INCLUDE "compbl.h" INCLUDE "YOETHF.h" INCLUDE "YOMCST.h" !**************************************************************************************** ! Start calculation ! - Initilalize output variables !**************************************************************************************** ycoefm(:,:) = 0.0 ycoefh(:,:) = 0.0 !**************************************************************************************** ! Methode 1 : ! !**************************************************************************************** CALL coefkz(nsrf, knon, ypaprs, ypplay, & ksta, ksta_ter, & yts, yrugos, yu, yv, yt, yq, & yqsurf, & ycoefm, ycoefh) !**************************************************************************************** ! Methode 2 : ! !**************************************************************************************** IF (iflag_pbl.EQ.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 !**************************************************************************************** ! IM cf JLD : on seuille ycoefm et ycoefh ! !**************************************************************************************** IF (nsrf.EQ.is_oce) THEN DO j=1,knon ycoefm(j,1)=MIN(ycoefm(j,1),cdmmax) ycoefh(j,1)=MIN(ycoefh(j,1),cdhmax) 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,ycoefm, & 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 ! Methode 3 : !**************************************************************************************** IF (iflag_pbl.GE.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 : ! y_cd_h(1:knon) = ycoefm(1:knon,1) ! y_cd_m(1:knon) = ycoefh(1:knon,1) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! y_cd_m(1:knon) = ycoefm(1:knon,1) y_cd_h(1:knon) = ycoefh(1:knon,1) CALL ustarhb(knon,yu,yv,y_cd_m, yustar) IF (prt_level > 9) THEN WRITE(lunout,*) 'USTAR = ',yustar ENDIF ! iflag_pbl peut etre utilise comme longuer de melange IF (iflag_pbl.GE.11) THEN CALL vdif_kcay(knon,dtime,RG,RD,ypaprs,yt, & yzlev,yzlay,yu,yv,yteta, & y_cd_m,yq2,q2diag,ykmm,ykmn,yustar, & iflag_pbl) ELSE CALL yamada4(knon,dtime,RG,RD,ypaprs,yt, & yzlev,yzlay,yu,yv,yteta, & y_cd_m,yq2,ykmm,ykmn,ykmq,yustar, & iflag_pbl) ENDIF ycoefm(1:knon,1)=y_cd_m(1:knon) ycoefh(1:knon,1)=y_cd_h(1:knon) ycoefm(1:knon,2:klev)=ykmm(1:knon,2:klev) ycoefh(1:knon,2:klev)=ykmn(1:knon,2:klev) ENDIF !(iflag_pbl.ge.3) END SUBROUTINE coef_diff_turb ! !**************************************************************************************** ! SUBROUTINE coefkz(nsrf, knon, paprs, pplay, & ksta, ksta_ter, & ts, rugos, & u,v,t,q, & qsurf, & pcfm, pcfh) !====================================================================== ! 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) ! rugos----input-R- longeur de rugosite (en m) ! 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 "FCTTRE.h" INCLUDE "iniprint.h" INCLUDE "indicesol.h" INCLUDE "compbl.h" INCLUDE "YOMCST.h" ! ! Arguments: ! INTEGER, INTENT(IN) :: knon, nsrf 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) :: rugos REAL, DIMENSION(klon,klev), INTENT(OUT) :: pcfm, pcfh ! Local variables ! numero de couche du sommet de la couche limite INTEGER, DIMENSION(klon) :: itop ! ! Quelques constantes et options: ! REAL cepdu2, ckap, cb, cc, cd, clam PARAMETER (cepdu2 =(0.1)**2) PARAMETER (CKAP=0.4) PARAMETER (cb=5.0) PARAMETER (cc=5.0) PARAMETER (cd=5.0) PARAMETER (clam=160.0) REAL ratqs ! largeur de distribution de vapeur d'eau PARAMETER (ratqs=0.05) LOGICAL richum ! utilise le nombre de Richardson humide PARAMETER (richum=.TRUE.) REAL ric ! nombre de Richardson critique PARAMETER(ric=0.4) REAL prandtl PARAMETER (prandtl=0.4) REAL kstable ! diffusion minimale (situation stable) ! GKtest ! PARAMETER (kstable=1.0e-10) REAL ksta, ksta_ter !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 mixlen ! constante controlant longueur de melange PARAMETER (mixlen=35.0) INTEGER isommet ! le sommet de la couche limite LOGICAL tvirtu ! calculer Ri d'une maniere plus performante PARAMETER (tvirtu=.TRUE.) LOGICAL opt_ec ! formule du Centre Europeen dans l'atmosphere PARAMETER (opt_ec=.FALSE.) ! ! Variables locales: INTEGER i, k !IM 120704 REAL zgeop(klon,klev) REAL zmgeom(klon) REAL zri(klon) REAL zl2(klon) REAL u1(klon), v1(klon), t1(klon), q1(klon), z1(klon) REAL pcfm1(klon), pcfh1(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 t_coup PARAMETER (t_coup=273.15) !IM LOGICAL check PARAMETER (check=.FALSE.) ! ! contre-gradient pour la chaleur sensible: Kelvin/metre REAL gamt(2:klev) REAL qsurf(klon) LOGICAL, SAVE :: appel1er !$OMP THREADPRIVATE(appel1er) ! ! Fonctions thermodynamiques et fonctions d'instabilite REAL fsta, fins, x LOGICAL zxli ! utiliser un jeu de fonctions simples PARAMETER (zxli=.FALSE.) fsta(x) = 1.0 / (1.0+10.0*x*(1+8.0*x)) fins(x) = SQRT(1.0-18.0*x) DATA appel1er /.TRUE./ 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.EQ.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 .NE. 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 le frottement au sol (Cdrag) ! DO i = 1, knon u1(i) = u(i,1) v1(i) = v(i,1) t1(i) = t(i,1) q1(i) = q(i,1) z1(i) = zgeop(i,1) ENDDO CALL clcdrag(klon, knon, nsrf, zxli, & u1, v1, t1, q1, z1, & ts, qsurf, rugos, & pcfm1, pcfh1) !IM ts, qsurf, rugos, DO i = 1, knon pcfm(i,1)=pcfm1(i) pcfh(i,1)=pcfh1(i) 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 .LT. 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): ! zfr = (zq+ratqs*zq-zqs) / (2.0*ratqs*zq) 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).LT.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 .LE. 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) !====================================================================== ! 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 prandtl PARAMETER (prandtl=0.4) REAL kstable PARAMETER (kstable=0.002) ! PARAMETER (kstable=0.001) REAL mixlen ! constante controlant longueur de melange PARAMETER (mixlen=35.0) REAL seuil ! au-dela l'inversion est consideree trop faible PARAMETER (seuil=-0.02) ! 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 "indicesol.h" 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).GT.0.8*paprs(i,1) .AND. & zdthdp.LT.zdthmin(i) ) THEN zdthmin(i) = zdthdp invb(i) = k ENDIF ENDDO ENDDO ! ! Introduire une diffusion: ! IF ( nsrf.EQ.is_oce ) THEN DO k = 2, klev DO i = 1, knon !IM cf FH/GK IF ( (nsrf.NE.is_oce) .OR. ! si ce n'est pas sur l'ocean !IM cf FH/GK . (invb(i).EQ.klev) .OR. ! s'il n'y a pas d'inversion !IM cf JLD/ GKtest TERkz2 ! IF ( (nsrf.EQ.is_ter) .OR. ! si on est sur la terre ! fin GKtest ! s'il n'y a pas d'inversion ou si l'inversion est trop faible ! IF ( (nsrf.EQ.is_oce) .AND. & IF ( (invb(i).EQ.klev) .OR. (zdthmin(i).GT.seuil) ) 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