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, ycdragm, ycdragh)

    ! 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(:, 2:) ! (knon, 2:klev) coefficient, vitesse

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

    real, intent(out):: ycdragm(:), ycdragh(:) ! (knon)

    ! 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, ycdragm, ycdragh) 

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