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 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
    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
    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, zcvm5, zcor, zqs, zfr, zdqs
    logical zdelta
    REAL z2geomf, zalh2, alm2, zscfh, scfm
    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:

          zdelta = RTT >=zt
          zcvm5 = merge(R5IES * RLSTT, R5LES * RLVTT, zdelta) / RCPD &
               / (1. + RVTMP2*zq)
          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)

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