Sources/phylmd/coefkz.f

module coefkz_m

  IMPLICIT none

contains

  SUBROUTINE coefkz(nsrf, 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: September 22nd, 1993
    ! Objet : calculer le coefficient de frottement du sol ("Cdrag") et les
    ! coefficients d'échange turbulent dans l'atmosphère.

    use clcdrag_m, only: clcdrag
    USE conf_phys_m, ONLY: iflag_pbl
    USE dimphy, ONLY: klev, klon
    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(:, :) ! (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(:) ! (knon) temperature du sol (en Kelvin)
    REAL, intent(in):: rugos(:) ! (klon) longeur de rugosite (en m)
    REAL, intent(in):: u(:, :), 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(:) ! (knon) 
    REAL, intent(out):: coefm(:, :) ! (knon, klev) coefficient, vitesse

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

    ! Local:

    INTEGER knon ! nombre de points a traiter
    INTEGER itop(size(coefm, 1)) ! (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

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

    knon = size(coefm, 1)

    ! 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(nsrf, 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	208	SUBROUTINE coefkz(nsrf, paprs, pplay, ksta, ksta_ter, ts, rugos, u, v, t, &
8			q, qsurf, 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	47	! 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	208	use clcdrag_m, only: clcdrag
16			USE conf_phys_m, ONLY: iflag_pbl
17			USE dimphy, ONLY: klev, klon
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	208	REAL, intent(in):: paprs(:, :) ! (klon, klev+1)
26	guez	47	! pression a chaque intercouche (en Pa)
27	guez	40
28	guez	208	real, intent(in):: pplay(:, :) ! (klon, 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	208	REAL, intent(in):: rugos(:) ! (klon) longeur de rugosite (en m)
34			REAL, intent(in):: u(:, :), v(:, :) ! (klon, klev) wind
35			REAL, intent(in):: t(:, :) ! (klon, klev) temperature (K)
36			real, intent(in):: q(:, :) ! (klon, klev) vapeur d'eau (kg/kg)
37	guez	221	real, intent(in):: qsurf(:) ! (knon)
38	guez	62	REAL, intent(out):: coefm(:, :) ! (knon, klev) coefficient, vitesse
39	guez	40
40	guez	62	real, intent(out):: coefh(:, :) ! (knon, klev)
41			! coefficient, chaleur et humidité
42
43	guez	47	! Local:
44	guez	40
45	guez	208	INTEGER knon ! nombre de points a traiter
46	guez	229	INTEGER itop(size(coefm, 1)) ! (knon) numero de couche du sommet
47			! de la couche limite
48	guez	40
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	47	REAL zgeop(klon, klev)
77			REAL zmgeom(klon)
78	guez	62	REAL ri(klon)
79			REAL l2(klon)
80	guez	40
81	guez	47	REAL u1(klon), v1(klon), t1(klon), q1(klon), z1(klon)
82	guez	40
83	guez	62	REAL zdphi, zdu2, ztvd, ztvu, cdn
84			REAL scf
85	guez	103	REAL zt, zq, zcvm5, zcor, zqs, zfr, zdqs
86			logical zdelta
87	guez	62	REAL z2geomf, zalh2, alm2, zscfh, scfm
88	guez	47	REAL gamt(2:klev) ! contre-gradient pour la chaleur sensible: Kelvin/metre
89	guez	40
90	guez	47	!--------------------------------------------------------------------
91	guez	40
92	guez	208	knon = size(coefm, 1)
93
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	221	CALL clcdrag(nsrf, u1, v1, t1, q1, z1, ts, qsurf, rugos, coefm(:, 1), &
136			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	207	zdelta = RTT >=zt
154			zcvm5 = merge(R5IES * RLSTT, R5LES * RLVTT, zdelta) / RCPD &
155			/ (1. + RVTMP2*zq)
156			zqs = R2ES * FOEEW(zt, zdelta) / pplay(i, k)
157			zqs = MIN(0.5, zqs)
158			zcor = 1./(1.-RETV*zqs)
159			zqs = zqs*zcor
160			zdqs = FOEDE(zt, zdelta, zcvm5, zqs, zcor)
161	guez	40
162	guez	47	! calculer la fraction nuageuse (processus humide):
163	guez	40
164	guez	47	zfr = (zq+ratqszq-zqs) / (2.0ratqs*zq)
165			zfr = MAX(0.0, MIN(1.0, zfr))
166			IF (.NOT.richum) zfr = 0.0
167	guez	40
168	guez	47	! calculer le nombre de Richardson:
169	guez	40
170	guez	47	IF (tvirtu) THEN
171			ztvd =( t(i, k) &
172			+ zdphi/RCPD/(1.+RVTMP2*zq) &
173			( (1.-zfr) + zfr(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) &
174			)(1.+RETVq(i, k))
175			ztvu =( t(i, k-1) &
176			- zdphi/RCPD/(1.+RVTMP2*zq) &
177			( (1.-zfr) + zfr(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) &
178			)(1.+RETVq(i, k-1))
179	guez	62	ri(i) =zmgeom(i)(ztvd-ztvu)/(zdu20.5*(ztvd+ztvu))
180			ri(i) = ri(i) &
181	guez	47	+ zmgeom(i)zmgeom(i)/RGgamt(k) &
182			(paprs(i, k)/101325.0)*RKAPPA &
183			/(zdu20.5(ztvd+ztvu))
184			ELSE
185			! calcul de Ridchardson compatible LMD5
186	guez	62	ri(i) =(RCPD*(t(i, k)-t(i, k-1)) &
187	guez	47	-RD0.5(t(i, k)+t(i, k-1))/paprs(i, k) &
188			*(pplay(i, k)-pplay(i, k-1)) &
189			)zmgeom(i)/(zdu20.5RCPD(t(i, k-1)+t(i, k)))
190	guez	62	ri(i) = ri(i) + &
191	guez	47	zmgeom(i)zmgeom(i)gamt(k)/RG &
192			(paprs(i, k)/101325.0)*RKAPPA &
193			/(zdu20.5(t(i, k-1)+t(i, k)))
194			ENDIF
195	guez	40
196	guez	47	! finalement, les coefficients d'echange sont obtenus:
197	guez	40
198	guez	62	cdn = SQRT(zdu2) / zmgeom(i) * RG
199	guez	40
200	guez	47	IF (opt_ec) THEN
201			z2geomf = zgeop(i, k-1)+zgeop(i, k)
202	guez	62	alm2 = (0.5ckap/RGz2geomf &
203	guez	47	/(1.+0.5ckap/rg/clamz2geomf))**2
204			zalh2 = (0.5ckap/rgz2geomf &
205			/(1.+0.5ckap/RG/(clamSQRT(1.5cd))z2geomf))**2
206	guez	62	IF (ri(i) < 0.) THEN
207	guez	47	! situation instable
208	guez	62	scf = ((zgeop(i, k)/zgeop(i, k-1))(1./3.)-1.)3 &
209	guez	47	/ (zmgeom(i)/RG)**3 / (zgeop(i, k-1)/RG)
210	guez	62	scf = SQRT(-ri(i)*scf)
211			scfm = 1.0 / (1.0+3.0cbccalm2scf)
212			zscfh = 1.0 / (1.0+3.0cbcczalh2scf)
213			coefm(i, k) = cdn * alm2 * (1. - 2. * cb * ri(i) * scfm)
214			coefh(i, k) = cdnzalh2(1.-3.0cbri(i)*zscfh)
215	guez	47	ELSE
216			! situation stable
217	guez	62	scf = SQRT(1.+cd*ri(i))
218			coefm(i, k) = cdn * alm2 / (1. + 2. * cb * ri(i) / scf)
219			coefh(i, k) = cdnzalh2/(1.+3.0cbri(i)scf)
220	guez	47	ENDIF
221			ELSE
222	guez	62	l2(i) = (mixlen*MAX(0.0, (paprs(i, k)-paprs(i, itop(i)+1)) &
223	guez	47	/(paprs(i, 2)-paprs(i, itop(i)+1)) ))**2
224	guez	62	coefm(i, k) = sqrt(max(cdn*2 (ric - ri(i)) / ric, kstable))
225			coefm(i, k)= l2(i) * coefm(i, k)
226			coefh(i, k) = coefm(i, k) / prandtl ! h et m different
227	guez	47	ENDIF
228			ENDDO loop_horiz
229			ENDDO loop_vertical
230	guez	40
231	guez	62	! Au-delà du sommet, pas de diffusion turbulente :
232			forall (i = 1: knon)
233			coefh(i, itop(i) + 1:) = 0.
234			coefm(i, itop(i) + 1:) = 0.
235			END forall
236	guez	40
237	guez	47	END SUBROUTINE coefkz
238	guez	40
239	guez	47	end module coefkz_m