Sources/phylmd/coefkz.f

module coefkz_m

  IMPLICIT none

contains

  SUBROUTINE coefkz(nsrf, knon, paprs, pplay, ksta, ksta_ter, ts, rugos, u, v, &
       t, q, qsurf, coefm, coefh)

    ! Authors: F. Hourdin, M. Forichon, Z. X. Li (LMD/CNRS)
    ! date: 1993/09/22
    ! Objet : calculer le coefficient de frottement du sol ("Cdrag") et les
    ! coefficients d'échange turbulent dans l'atmosphère.

    USE indicesol, ONLY: is_oce
    USE dimphy, ONLY: klev, klon
    USE conf_gcm_m, ONLY: prt_level
    USE suphec_m, ONLY: rcpd, rd, retv, rg, rkappa, rlstt, rlvtt, rtt
    USE yoethf_m, ONLY: r2es, r5ies, r5les, rvtmp2
    USE fcttre, ONLY: dqsatl, dqsats, foede, foeew, qsatl, qsats, thermcep
    USE conf_phys_m, ONLY: iflag_pbl
    use clcdrag_m, only: clcdrag

    ! Arguments:

    integer, intent(in):: nsrf ! indicateur de la nature du sol
    INTEGER, intent(in):: knon ! nombre de points a traiter

    REAL, intent(in):: paprs(klon, klev+1)
    ! pression a chaque intercouche (en Pa)

    real, intent(in):: pplay(klon, klev)
    ! pression au milieu de chaque couche (en Pa)

    REAL, intent(in):: ksta, ksta_ter
    REAL, intent(in):: ts(klon) ! temperature du sol (en Kelvin)
    REAL, intent(in):: rugos(klon) ! longeur de rugosite (en m)
    REAL, intent(in):: u(klon, klev), v(klon, klev) ! wind
    REAL, intent(in):: t(klon, klev) ! temperature (K)
    real, intent(in):: q(klon, klev) ! vapeur d'eau (kg/kg)
    real, intent(in):: qsurf(klon) 
    REAL, intent(out):: coefm(:, :) ! (knon, klev) coefficient, vitesse

    real, intent(out):: coefh(:, :) ! (knon, klev) 
    ! coefficient, chaleur et humidité

    ! Local:

    INTEGER itop(knon) ! 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.
    REAL, PARAMETER:: cc = 5.
    REAL, PARAMETER:: cd = 5.
    REAL, PARAMETER:: clam = 160.
    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)
    REAL, PARAMETER:: mixlen = 35. ! constante contrôlant longueur de mélange
    INTEGER, PARAMETER:: isommet = klev ! 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

    INTEGER i, k, kk
    REAL zgeop(klon, klev)
    REAL zmgeom(klon)
    REAL ri(klon)
    REAL l2(klon)

    REAL u1(klon), v1(klon), t1(klon), q1(klon), z1(klon)

    REAL zdphi, zdu2, ztvd, ztvu, cdn
    REAL scf
    REAL zt, zq, zdelta, zcvm5, zcor, zqs, zfr, zdqs
    REAL z2geomf, zalh2, alm2, zscfh, scfm
    REAL, PARAMETER:: t_coup = 273.15
    REAL gamt(2:klev) ! contre-gradient pour la chaleur sensible: Kelvin/metre

    !--------------------------------------------------------------------

    ! 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

    IF ( nsrf .NE. is_oce ) THEN
       kstable = ksta_ter
    ELSE
       kstable = ksta
    ENDIF

    ! Calculer les géopotentiels 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, .false., u1, v1, t1, q1, z1, ts, qsurf, &
         rugos, coefm(:, 1), coefh(:, 1)) 

    ! Calculer les coefficients turbulents dans l'atmosphere

    itop = isommet

    loop_vertical: DO k = 2, isommet
       loop_horiz: 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):

          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))
             ri(i) =zmgeom(i)*(ztvd-ztvu)/(zdu2*0.5*(ztvd+ztvu))
             ri(i) = ri(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
             ri(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)))
             ri(i) = ri(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:

          cdn = SQRT(zdu2) / zmgeom(i) * RG

          IF (opt_ec) THEN
             z2geomf = zgeop(i, k-1)+zgeop(i, k)
             alm2 = (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 (ri(i) < 0.) THEN
                ! situation instable
                scf = ((zgeop(i, k)/zgeop(i, k-1))**(1./3.)-1.)**3 &
                     / (zmgeom(i)/RG)**3 / (zgeop(i, k-1)/RG)
                scf = SQRT(-ri(i)*scf)
                scfm = 1.0 / (1.0+3.0*cb*cc*alm2*scf)
                zscfh = 1.0 / (1.0+3.0*cb*cc*zalh2*scf)
                coefm(i, k) = cdn * alm2 * (1. - 2. * cb * ri(i) * scfm)
                coefh(i, k) = cdn*zalh2*(1.-3.0*cb*ri(i)*zscfh)
             ELSE
                ! situation stable
                scf = SQRT(1.+cd*ri(i))
                coefm(i, k) = cdn * alm2 / (1. + 2. * cb * ri(i) / scf)
                coefh(i, k) = cdn*zalh2/(1.+3.0*cb*ri(i)*scf)
             ENDIF
          ELSE
             l2(i) = (mixlen*MAX(0.0, (paprs(i, k)-paprs(i, itop(i)+1)) &
                  /(paprs(i, 2)-paprs(i, itop(i)+1)) ))**2
             coefm(i, k) = sqrt(max(cdn**2 * (ric - ri(i)) / ric, kstable))
             coefm(i, k)= l2(i) * coefm(i, k)
             coefh(i, k) = coefm(i, k) / prandtl ! h et m different
          ENDIF
       ENDDO loop_horiz
    ENDDO loop_vertical

    ! Au-delà du sommet, pas de diffusion turbulente :
    forall (i = 1: knon)
       coefh(i, itop(i) + 1:) = 0.
       coefm(i, itop(i) + 1:) = 0.
    END forall

  END SUBROUTINE coefkz

end module coefkz_m
1	guez	47	module coefkz_m
2	guez	3
3	guez	40	IMPLICIT none
4	guez	3
5	guez	47	contains
6	guez	3
7	guez	47	SUBROUTINE coefkz(nsrf, knon, paprs, pplay, ksta, ksta_ter, ts, rugos, u, v, &
8	guez	62	t, q, qsurf, coefm, coefh)
9	guez	40
10	guez	62	! Authors: F. Hourdin, M. Forichon, Z. X. Li (LMD/CNRS)
11	guez	47	! date: 1993/09/22
12			! Objet : calculer le coefficient de frottement du sol ("Cdrag") et les
13			! coefficients d'échange turbulent dans l'atmosphère.
14	guez	40
15	guez	62	USE indicesol, ONLY: is_oce
16	guez	71	USE dimphy, ONLY: klev, klon
17	guez	62	USE conf_gcm_m, ONLY: prt_level
18			USE suphec_m, ONLY: rcpd, rd, retv, rg, rkappa, rlstt, rlvtt, rtt
19			USE yoethf_m, ONLY: r2es, r5ies, r5les, rvtmp2
20			USE fcttre, ONLY: dqsatl, dqsats, foede, foeew, qsatl, qsats, thermcep
21			USE conf_phys_m, ONLY: iflag_pbl
22			use clcdrag_m, only: clcdrag
23	guez	40
24	guez	47	! Arguments:
25	guez	40
26	guez	47	integer, intent(in):: nsrf ! indicateur de la nature du sol
27			INTEGER, intent(in):: knon ! nombre de points a traiter
28	guez	40
29	guez	47	REAL, intent(in):: paprs(klon, klev+1)
30			! pression a chaque intercouche (en Pa)
31	guez	40
32	guez	47	real, intent(in):: pplay(klon, klev)
33			! pression au milieu de chaque couche (en Pa)
34	guez	40
35	guez	47	REAL, intent(in):: ksta, ksta_ter
36			REAL, intent(in):: ts(klon) ! temperature du sol (en Kelvin)
37			REAL, intent(in):: rugos(klon) ! longeur de rugosite (en m)
38			REAL, intent(in):: u(klon, klev), v(klon, klev) ! wind
39			REAL, intent(in):: t(klon, klev) ! temperature (K)
40			real, intent(in):: q(klon, klev) ! vapeur d'eau (kg/kg)
41			real, intent(in):: qsurf(klon)
42	guez	62	REAL, intent(out):: coefm(:, :) ! (knon, klev) coefficient, vitesse
43	guez	40
44	guez	62	real, intent(out):: coefh(:, :) ! (knon, klev)
45			! coefficient, chaleur et humidité
46
47	guez	47	! Local:
48	guez	40
49	guez	62	INTEGER itop(knon) ! numero de couche du sommet de la couche limite
50	guez	40
51	guez	47	! Quelques constantes et options:
52	guez	40
53	guez	47	REAL, PARAMETER:: cepdu2 =0.1**2
54			REAL, PARAMETER:: CKAP = 0.4
55			REAL, PARAMETER:: cb = 5.
56			REAL, PARAMETER:: cc = 5.
57			REAL, PARAMETER:: cd = 5.
58			REAL, PARAMETER:: clam = 160.
59	guez	62	REAL, PARAMETER:: ratqs = 0.05 ! largeur de distribution de vapeur d'eau
60
61			LOGICAL, PARAMETER:: richum = .TRUE.
62			! utilise le nombre de Richardson humide
63
64	guez	47	REAL, PARAMETER:: ric = 0.4 ! nombre de Richardson critique
65			REAL, PARAMETER:: prandtl = 0.4
66	guez	40
67	guez	47	REAL kstable ! diffusion minimale (situation stable)
68	guez	62	REAL, PARAMETER:: mixlen = 35. ! constante contrôlant longueur de mélange
69			INTEGER, PARAMETER:: isommet = klev ! sommet de la couche limite
70	guez	40
71	guez	47	LOGICAL, PARAMETER:: tvirtu = .TRUE.
72			! calculer Ri d'une maniere plus performante
73	guez	40
74	guez	47	LOGICAL, PARAMETER:: opt_ec = .FALSE.
75			! formule du Centre Europeen dans l'atmosphere
76	guez	40
77	guez	47	INTEGER i, k, kk
78			REAL zgeop(klon, klev)
79			REAL zmgeom(klon)
80	guez	62	REAL ri(klon)
81			REAL l2(klon)
82	guez	40
83	guez	47	REAL u1(klon), v1(klon), t1(klon), q1(klon), z1(klon)
84	guez	40
85	guez	62	REAL zdphi, zdu2, ztvd, ztvu, cdn
86			REAL scf
87	guez	47	REAL zt, zq, zdelta, zcvm5, zcor, zqs, zfr, zdqs
88	guez	62	REAL z2geomf, zalh2, alm2, zscfh, scfm
89	guez	47	REAL, PARAMETER:: t_coup = 273.15
90			REAL gamt(2:klev) ! contre-gradient pour la chaleur sensible: Kelvin/metre
91	guez	40
92	guez	47	!--------------------------------------------------------------------
93	guez	40
94	guez	47	! Prescrire la valeur de contre-gradient
95			if (iflag_pbl.eq.1) then
96			DO k = 3, klev
97			gamt(k) = -1.0E-03
98			ENDDO
99			gamt(2) = -2.5E-03
100			else
101			DO k = 2, klev
102			gamt(k) = 0.0
103			ENDDO
104			ENDIF
105	guez	40
106	guez	47	IF ( nsrf .NE. is_oce ) THEN
107			kstable = ksta_ter
108			ELSE
109			kstable = ksta
110			ENDIF
111	guez	40
112	guez	47	! Calculer les géopotentiels de chaque couche
113			DO i = 1, knon
114			zgeop(i, 1) = RD * t(i, 1) / (0.5 * (paprs(i, 1) + pplay(i, 1))) &
115			* (paprs(i, 1) - pplay(i, 1))
116			ENDDO
117			DO k = 2, klev
118			DO i = 1, knon
119			zgeop(i, k) = zgeop(i, k-1) &
120			+ RD * 0.5*(t(i, k-1)+t(i, k)) / paprs(i, k) &
121			* (pplay(i, k-1)-pplay(i, k))
122			ENDDO
123			ENDDO
124	guez	40
125	guez	47	! Calculer le frottement au sol (Cdrag)
126	guez	40
127	guez	47	DO i = 1, knon
128			u1(i) = u(i, 1)
129			v1(i) = v(i, 1)
130			t1(i) = t(i, 1)
131			q1(i) = q(i, 1)
132			z1(i) = zgeop(i, 1)
133			ENDDO
134	guez	40
135	guez	47	CALL clcdrag(klon, knon, nsrf, .false., u1, v1, t1, q1, z1, ts, qsurf, &
136	guez	62	rugos, coefm(:, 1), coefh(:, 1))
137	guez	40
138	guez	47	! Calculer les coefficients turbulents dans l'atmosphere
139	guez	40
140	guez	62	itop = isommet
141	guez	40
142	guez	47	loop_vertical: DO k = 2, isommet
143			loop_horiz: DO i = 1, knon
144			zdu2 = MAX(cepdu2, (u(i, k)-u(i, k-1))**2 &
145			+(v(i, k)-v(i, k-1))**2)
146			zmgeom(i) = zgeop(i, k)-zgeop(i, k-1)
147			zdphi =zmgeom(i) / 2.0
148			zt = (t(i, k)+t(i, k-1)) * 0.5
149			zq = (q(i, k)+q(i, k-1)) * 0.5
150	guez	40
151	guez	47	! calculer Qs et dQs/dT:
152	guez	40
153	guez	47	IF (thermcep) THEN
154			zdelta = MAX(0., SIGN(1., RTT-zt))
155			zcvm5 = R5LESRLVTT/RCPD/(1.0+RVTMP2zq)*(1.-zdelta) &
156			+ R5IESRLSTT/RCPD/(1.0+RVTMP2zq)*zdelta
157			zqs = R2ES * FOEEW(zt, zdelta) / pplay(i, k)
158			zqs = MIN(0.5, zqs)
159			zcor = 1./(1.-RETV*zqs)
160			zqs = zqs*zcor
161			zdqs = FOEDE(zt, zdelta, zcvm5, zqs, zcor)
162			ELSE
163	guez	62	IF (zt < t_coup) THEN
164	guez	47	zqs = qsats(zt) / pplay(i, k)
165			zdqs = dqsats(zt, zqs)
166			ELSE
167			zqs = qsatl(zt) / pplay(i, k)
168			zdqs = dqsatl(zt, zqs)
169			ENDIF
170			ENDIF
171	guez	40
172	guez	47	! calculer la fraction nuageuse (processus humide):
173	guez	40
174	guez	47	zfr = (zq+ratqszq-zqs) / (2.0ratqs*zq)
175			zfr = MAX(0.0, MIN(1.0, zfr))
176			IF (.NOT.richum) zfr = 0.0
177	guez	40
178	guez	47	! calculer le nombre de Richardson:
179	guez	40
180	guez	47	IF (tvirtu) THEN
181			ztvd =( t(i, k) &
182			+ zdphi/RCPD/(1.+RVTMP2*zq) &
183			( (1.-zfr) + zfr(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) &
184			)(1.+RETVq(i, k))
185			ztvu =( t(i, k-1) &
186			- zdphi/RCPD/(1.+RVTMP2*zq) &
187			( (1.-zfr) + zfr(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) &
188			)(1.+RETVq(i, k-1))
189	guez	62	ri(i) =zmgeom(i)(ztvd-ztvu)/(zdu20.5*(ztvd+ztvu))
190			ri(i) = ri(i) &
191	guez	47	+ zmgeom(i)zmgeom(i)/RGgamt(k) &
192			(paprs(i, k)/101325.0)*RKAPPA &
193			/(zdu20.5(ztvd+ztvu))
194			ELSE
195			! calcul de Ridchardson compatible LMD5
196	guez	62	ri(i) =(RCPD*(t(i, k)-t(i, k-1)) &
197	guez	47	-RD0.5(t(i, k)+t(i, k-1))/paprs(i, k) &
198			*(pplay(i, k)-pplay(i, k-1)) &
199			)zmgeom(i)/(zdu20.5RCPD(t(i, k-1)+t(i, k)))
200	guez	62	ri(i) = ri(i) + &
201	guez	47	zmgeom(i)zmgeom(i)gamt(k)/RG &
202			(paprs(i, k)/101325.0)*RKAPPA &
203			/(zdu20.5(t(i, k-1)+t(i, k)))
204			ENDIF
205	guez	40
206	guez	47	! finalement, les coefficients d'echange sont obtenus:
207	guez	40
208	guez	62	cdn = SQRT(zdu2) / zmgeom(i) * RG
209	guez	40
210	guez	47	IF (opt_ec) THEN
211			z2geomf = zgeop(i, k-1)+zgeop(i, k)
212	guez	62	alm2 = (0.5ckap/RGz2geomf &
213	guez	47	/(1.+0.5ckap/rg/clamz2geomf))**2
214			zalh2 = (0.5ckap/rgz2geomf &
215			/(1.+0.5ckap/RG/(clamSQRT(1.5cd))z2geomf))**2
216	guez	62	IF (ri(i) < 0.) THEN
217	guez	47	! situation instable
218	guez	62	scf = ((zgeop(i, k)/zgeop(i, k-1))(1./3.)-1.)3 &
219	guez	47	/ (zmgeom(i)/RG)**3 / (zgeop(i, k-1)/RG)
220	guez	62	scf = SQRT(-ri(i)*scf)
221			scfm = 1.0 / (1.0+3.0cbccalm2scf)
222			zscfh = 1.0 / (1.0+3.0cbcczalh2scf)
223			coefm(i, k) = cdn * alm2 * (1. - 2. * cb * ri(i) * scfm)
224			coefh(i, k) = cdnzalh2(1.-3.0cbri(i)*zscfh)
225	guez	47	ELSE
226			! situation stable
227	guez	62	scf = SQRT(1.+cd*ri(i))
228			coefm(i, k) = cdn * alm2 / (1. + 2. * cb * ri(i) / scf)
229			coefh(i, k) = cdnzalh2/(1.+3.0cbri(i)scf)
230	guez	47	ENDIF
231			ELSE
232	guez	62	l2(i) = (mixlen*MAX(0.0, (paprs(i, k)-paprs(i, itop(i)+1)) &
233	guez	47	/(paprs(i, 2)-paprs(i, itop(i)+1)) ))**2
234	guez	62	coefm(i, k) = sqrt(max(cdn*2 (ric - ri(i)) / ric, kstable))
235			coefm(i, k)= l2(i) * coefm(i, k)
236			coefh(i, k) = coefm(i, k) / prandtl ! h et m different
237	guez	47	ENDIF
238			ENDDO loop_horiz
239			ENDDO loop_vertical
240	guez	40
241	guez	62	! Au-delà du sommet, pas de diffusion turbulente :
242			forall (i = 1: knon)
243			coefh(i, itop(i) + 1:) = 0.
244			coefm(i, itop(i) + 1:) = 0.
245			END forall
246	guez	40
247	guez	47	END SUBROUTINE coefkz
248	guez	40
249	guez	47	end module coefkz_m