source: LMDZ5/trunk/libf/phylmd/coefcdrag.F90 @ 2201

Last change on this file since 2201 was 2158, checked in by Laurent Fairhead, 10 years ago

Modification pour éviter de faire des calculs avec un q négatif qui peut entraîner
des calculs de t2m irréalistes (raison du plantage avec ORCHIDEE en mode debug a priori)


Modifications to prevent calculations with a negative humidity that would lead to
unrealistic calculations for temperature

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 5.8 KB
Line 
1!
2!
3!
4      SUBROUTINE coefcdrag (klon, knon, nsrf, zxli, &
5                            speed, t, q, zgeop, psol, &
6                            ts, qsurf, rugos, okri, ri1, &
7                            cdram, cdrah, cdran, zri1, pref)
8
9      USE indice_sol_mod
10
11      IMPLICIT none
12!-------------------------------------------------------------------------
13! Objet : calcul des cdrags pour le moment (cdram) et les flux de chaleur
14!         sensible et latente (cdrah), du cdrag neutre (cdran),
15!         du nombre de Richardson entre la surface et le niveau de reference
16!         (zri1) et de la pression au niveau de reference (pref).   
17!
18! I. Musat, 01.07.2002
19!-------------------------------------------------------------------------
20!
21! klon----input-I- dimension de la grille physique (= nb_pts_latitude X nb_pts_longitude)
22! knon----input-I- nombre de points pour un type de surface
23! nsrf----input-I- indice pour le type de surface; voir indice_sol_mod.F90
24! zxli----input-L- TRUE si calcul des cdrags selon Laurent Li
25! speed---input-R- module du vent au 1er niveau du modele
26! t-------input-R- temperature de l'air au 1er niveau du modele
27! q-------input-R- humidite de l'air au 1er niveau du modele
28! zgeop---input-R- geopotentiel au 1er niveau du modele
29! psol----input-R- pression au sol
30! ts------input-R- temperature de l'air a la surface
31! qsurf---input-R- humidite de l'air a la surface
32! rugos---input-R- rugosite
33! okri----input-L- TRUE si on veut tester le nb. Richardson entre la sfce
34!                  et zref par rapport au Ri entre la sfce et la 1ere couche
35! ri1-----input-R- nb. Richardson entre la surface et la 1ere couche
36!
37! cdram--output-R- cdrag pour le moment
38! cdrah--output-R- cdrag pour les flux de chaleur latente et sensible
39! cdran--output-R- cdrag neutre
40! zri1---output-R- nb. Richardson entre la surface et la couche zgeop/RG
41! pref---output-R- pression au niveau zgeop/RG
42!
43      INTEGER, intent(in) :: klon, knon, nsrf
44      LOGICAL, intent(in) :: zxli
45      REAL, dimension(klon), intent(in) :: speed, t, q, zgeop, psol
46      REAL, dimension(klon), intent(in) :: ts, qsurf, rugos, ri1
47      LOGICAL, intent(in) :: okri   
48!
49      REAL, dimension(klon), intent(out) :: cdram, cdrah, cdran, zri1, pref
50!-------------------------------------------------------------------------
51!
52      include "YOMCST.h"
53      include "YOETHF.h"
54      INCLUDE "clesphys.h"
55! Quelques constantes :
56      REAL, parameter :: RKAR=0.40, CB=5.0, CC=5.0, CD=5.0, cepdu2=(0.1)**2
57!
58! Variables locales :
59      INTEGER :: i
60      REAL, dimension(klon) :: zdu2, zdphi, ztsolv, ztvd
61      REAL, dimension(klon) :: zscf, friv, frih, zucf, zcr
62      REAL, dimension(klon) :: zcfm1, zcfh1
63      REAL, dimension(klon) :: zcfm2, zcfh2
64      REAL, dimension(klon) :: trm0, trm1
65
66!-------------------------------------------------------------------------
67      REAL :: fsta, fins, x
68      fsta(x) = 1.0 / (1.0+10.0*x*(1+8.0*x))
69      fins(x) = SQRT(1.0-18.0*x)
70!-------------------------------------------------------------------------
71!
72      DO i = 1, knon
73!
74       zdphi(i) = zgeop(i)
75       zdu2(i) = max(cepdu2,speed(i)**2)
76       pref(i) = exp(log(psol(i)) - zdphi(i)/(RD*t(i)* &
77                 (1.+ RETV * max(q(i),0.0))))
78       ztsolv(i) = ts(i)
79!       ztvd(i) = t(i) * (psol(i)/pref(i))**RKAPPA
80!       ztvd(i) = (t(i)+zdphi(i)/RCPD/(1.+RVTMP2*q(i))) &
81!          *(1.+RETV*q(i))
82       ztvd(i) = (t(i)+zdphi(i)/RCPD/(1.+RVTMP2*q(i)))
83       trm0(i) = 1. + RETV * max(qsurf(i),0.0)
84       trm1(i) = 1. + RETV * max(q(i),0.0)
85       ztsolv(i) = ztsolv(i) * trm0(i)
86       ztvd(i) = ztvd(i) * trm1(i)
87       zri1(i) = zdphi(i)*(ztvd(i)-ztsolv(i))/(zdu2(i)*ztvd(i))
88!
89! on teste zri1 par rapport au Richardson de la 1ere couche ri1
90!
91!IM +++
92       IF(1.EQ.0) THEN
93       IF (okri) THEN
94         IF (ri1(i).GE.0.0.AND.zri1(i).LT.0.0) THEN
95           zri1(i) = ri1(i)
96         ELSE IF(ri1(i).LT.0.0.AND.zri1(i).GE.0.0) THEN
97           zri1(i) = ri1(i)
98         ENDIF
99       ENDIF
100       ENDIF
101!IM ---
102!
103       cdran(i) = (RKAR/log(1.+zdphi(i)/(RG*rugos(i))))**2
104
105       IF (zri1(i) .ge. 0.) THEN
106!
107! situation stable : pour eviter les inconsistances dans les cas
108! tres stables on limite zri1 a 20. cf Hess et al. (1995)
109!
110         zri1(i) = min(20.,zri1(i))
111!
112         IF (.NOT.zxli) THEN
113           zscf(i) = SQRT(1.+CD*ABS(zri1(i)))
114           friv(i) = max(1. / (1.+2.*CB*zri1(i)/ zscf(i)), f_ri_cd_min)
115           zcfm1(i) = cdran(i) * friv(i)
116           frih(i) = max(1./ (1.+3.*CB*zri1(i)*zscf(i)), f_ri_cd_min )
117!           zcfh1(i) = cdran(i) * frih(i)
118           zcfh1(i) = f_cdrag_ter*cdran(i) * frih(i)
119           IF(nsrf.EQ.is_oce) zcfh1(i)=f_cdrag_oce*cdran(i)*frih(i)
120           cdram(i) = zcfm1(i)
121           cdrah(i) = zcfh1(i)
122         ELSE
123           cdram(i) = cdran(i)* fsta(zri1(i))
124           cdrah(i) = cdran(i)* fsta(zri1(i))
125         ENDIF
126!
127       ELSE
128!
129! situation instable
130!
131         IF (.NOT.zxli) THEN
132           zucf(i) = 1./(1.+3.0*CB*CC*cdran(i)*SQRT(ABS(zri1(i)) &
133                 *(1.0+zdphi(i)/(RG*rugos(i)))))
134           zcfm2(i) = cdran(i)*max((1.-2.0*CB*zri1(i)*zucf(i)),f_ri_cd_min)
135!           zcfh2(i) = cdran(i)*max((1.-3.0*CB*zri1(i)*zucf(i)),f_ri_cd_min)
136           zcfh2(i) = f_cdrag_ter*cdran(i)*max((1.-3.0*CB*zri1(i)*zucf(i)),f_ri_cd_min)
137           cdram(i) = zcfm2(i)
138           cdrah(i) = zcfh2(i)
139         ELSE
140           cdram(i) = cdran(i)* fins(zri1(i))
141           cdrah(i) = cdran(i)* fins(zri1(i))
142         ENDIF
143!
144! cdrah sur l'ocean cf. Miller et al. (1992)
145!
146         zcr(i) = (0.0016/(cdran(i)*SQRT(zdu2(i))))*ABS(ztvd(i)-ztsolv(i)) &
147               **(1./3.)
148!         IF (nsrf.EQ.is_oce) cdrah(i) = cdran(i)*(1.0+zcr(i)**1.25) &
149!                  **(1./1.25)
150         IF (nsrf.EQ.is_oce) cdrah(i)=f_cdrag_oce*cdran(i)*(1.0+zcr(i)**1.25) &
151                  **(1./1.25)
152       ENDIF
153!
154      END DO
155      RETURN
156      END SUBROUTINE coefcdrag
Note: See TracBrowser for help on using the repository browser.