Sources/phylmd/coefkz.f

module coefkz_m

  IMPLICIT none

contains

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