Sources/phylmd/coefkz.f

module coefkz_m

  IMPLICIT none

contains

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

    ! Authors: F. Hourdin, M. Forichon, Z. X. Li (LMD/CNRS)
    ! Date: September 22nd, 1993

    ! Objet : calculer les coefficients d'échange turbulent dans
    ! l'atmosphère.

    USE conf_phys_m, ONLY: iflag_pbl
    USE dimphy, ONLY: klev
    USE fcttre, ONLY: foede, foeew
    USE indicesol, ONLY: is_oce
    USE suphec_m, ONLY: rcpd, rd, retv, rg, rkappa, rlstt, rlvtt, rtt
    USE yoethf_m, ONLY: r2es, r5ies, r5les, rvtmp2

    integer, intent(in):: nsrf ! indicateur de la nature du sol

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

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

    REAL, intent(in):: ksta, ksta_ter
    REAL, intent(in):: ts(:) ! (knon) temperature du sol (en Kelvin)
    REAL, intent(in):: u(:, :), v(:, :) ! (knon, klev) wind
    REAL, intent(in):: t(:, :) ! (knon, klev) temperature (K)
    real, intent(in):: q(:, :) ! (knon, klev) vapeur d'eau (kg/kg)
    REAL, intent(in):: zgeop(:, :) ! (knon, klev)
    REAL, intent(out):: coefm(:, 2:) ! (knon, 2:klev) coefficient, vitesse

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

    ! Local:

    INTEGER knon ! nombre de points a traiter

    INTEGER itop(size(ts)) ! (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 zmgeom(size(ts))
    REAL ri(size(ts))
    REAL l2(size(ts))
    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

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

    knon = size(ts)

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