phylmd/Interface_surf/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:: ratqs = 0.05 ! largeur de distribution de vapeur d'eau
    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
    INTEGER i, k
    REAL zmgeom(size(ts))
    REAL ri(size(ts))
    REAL l2(size(ts))
    REAL zdphi, zdu2, ztvd, ztvu, cdn
    REAL zt, zq, zcvm5, zcor, zqs, zfr, zdqs
    logical zdelta
    REAL gamt(2:klev) ! contre-gradient pour la chaleur sensible: Kelvin/metre

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

    knon = size(ts)

    ! Prescrire la valeur de contre-gradient
    if (iflag_pbl == 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

    kstable = merge(ksta, ksta_ter, nsrf == is_oce)

    ! Calculer les coefficients turbulents dans l'atmosphere

    itop = isommet

    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:
          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))

          ! calculer le nombre de Richardson:
          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))

          ! finalement, les coefficients d'echange sont obtenus:

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

          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
       ENDDO
    ENDDO

    ! 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	279	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	guez	62	REAL, PARAMETER:: ratqs = 0.05 ! largeur de distribution de vapeur d'eau
52	guez	47	REAL, PARAMETER:: ric = 0.4 ! nombre de Richardson critique
53			REAL, PARAMETER:: prandtl = 0.4
54	guez	40
55	guez	47	REAL kstable ! diffusion minimale (situation stable)
56	guez	62	REAL, PARAMETER:: mixlen = 35. ! constante contrôlant longueur de mélange
57			INTEGER, PARAMETER:: isommet = klev ! sommet de la couche limite
58	guez	105	INTEGER i, k
59	guez	248	REAL zmgeom(size(ts))
60			REAL ri(size(ts))
61			REAL l2(size(ts))
62	guez	62	REAL zdphi, zdu2, ztvd, ztvu, cdn
63	guez	103	REAL zt, zq, zcvm5, zcor, zqs, zfr, zdqs
64			logical zdelta
65	guez	47	REAL gamt(2:klev) ! contre-gradient pour la chaleur sensible: Kelvin/metre
66	guez	40
67	guez	47	!--------------------------------------------------------------------
68	guez	40
69	guez	248	knon = size(ts)
70	guez	208
71	guez	47	! Prescrire la valeur de contre-gradient
72	guez	279	if (iflag_pbl == 1) then
73	guez	47	DO k = 3, klev
74	guez	248	gamt(k) = - 1E-3
75	guez	47	ENDDO
76	guez	248	gamt(2) = - 2.5E-3
77	guez	47	else
78			DO k = 2, klev
79			gamt(k) = 0.0
80			ENDDO
81			ENDIF
82	guez	40
83	guez	288	kstable = merge(ksta, ksta_ter, nsrf == is_oce)
84	guez	40
85	guez	47	! Calculer les coefficients turbulents dans l'atmosphere
86	guez	40
87	guez	62	itop = isommet
88	guez	40
89	guez	288	DO k = 2, isommet
90			DO i = 1, knon
91	guez	248	zdu2 = MAX(cepdu2, (u(i, k) - u(i, k - 1))**2 &
92			+ (v(i, k) - v(i, k - 1))**2)
93			zmgeom(i) = zgeop(i, k) - zgeop(i, k - 1)
94			zdphi = zmgeom(i) / 2.0
95			zt = (t(i, k) + t(i, k - 1)) * 0.5
96			zq = (q(i, k) + q(i, k - 1)) * 0.5
97	guez	40
98	guez	47	! calculer Qs et dQs/dT:
99	guez	207	zdelta = RTT >=zt
100			zcvm5 = merge(R5IES * RLSTT, R5LES * RLVTT, zdelta) / RCPD &
101	guez	248	/ (1. + RVTMP2 * zq)
102	guez	207	zqs = R2ES * FOEEW(zt, zdelta) / pplay(i, k)
103			zqs = MIN(0.5, zqs)
104	guez	248	zcor = 1./(1. - RETV * zqs)
105			zqs = zqs * zcor
106	guez	207	zdqs = FOEDE(zt, zdelta, zcvm5, zqs, zcor)
107	guez	40
108	guez	47	! calculer la fraction nuageuse (processus humide):
109	guez	248	zfr = (zq + ratqs * zq - zqs) / (2.0 * ratqs * zq)
110	guez	47	zfr = MAX(0.0, MIN(1.0, zfr))
111	guez	40
112	guez	279	! calculer le nombre de Richardson:
113	guez	288	ztvd = (t(i, k) &
114			+ zdphi/RCPD/(1. + RVTMP2 * zq) &
115			* ((1. - zfr) + zfr * (1. + RLVTT * zqs/RD/zt)/(1. + zdqs)) &
116			) * (1. + RETV * q(i, k))
117			ztvu = (t(i, k - 1) &
118			- zdphi/RCPD/(1. + RVTMP2 * zq) &
119			* ((1. - zfr) + zfr * (1. + RLVTT * zqs/RD/zt)/(1. + zdqs)) &
120			) * (1. + RETV * q(i, k - 1))
121			ri(i) = zmgeom(i) * (ztvd - ztvu)/(zdu2 * 0.5 * (ztvd + ztvu))
122			ri(i) = ri(i) &
123			+ zmgeom(i) * zmgeom(i)/RG * gamt(k) &
124			* (paprs(i, k)/101325.0)**RKAPPA &
125			/(zdu2 * 0.5 * (ztvd + ztvu))
126	guez	40
127	guez	47	! finalement, les coefficients d'echange sont obtenus:
128	guez	40
129	guez	62	cdn = SQRT(zdu2) / zmgeom(i) * RG
130	guez	40
131	guez	288	l2(i) = (mixlen * MAX(0.0, (paprs(i, k) - paprs(i, itop(i) + 1)) &
132			/(paprs(i, 2) - paprs(i, itop(i) + 1))))**2
133			coefm(i, k) = sqrt(max(cdn*2 (ric - ri(i)) / ric, kstable))
134			coefm(i, k) = l2(i) * coefm(i, k)
135			coefh(i, k) = coefm(i, k) / prandtl ! h et m different
136			ENDDO
137			ENDDO
138	guez	40
139	guez	62	! Au-delà du sommet, pas de diffusion turbulente :
140			forall (i = 1: knon)
141			coefh(i, itop(i) + 1:) = 0.
142			coefm(i, itop(i) + 1:) = 0.
143			END forall
144	guez	40
145	guez	47	END SUBROUTINE coefkz
146	guez	40
147	guez	47	end module coefkz_m