source: LMDZ5/branches/LMDZ5-DOFOCO/libf/phylmd/coefcdrag.F90 @ 4400

Last change on this file since 4400 was 1940, checked in by jghattas, 11 years ago

Bug corrections done by Fuxing Wang.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 5.3 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! Quelques constantes :
55      REAL, parameter :: RKAR=0.40, CB=5.0, CC=5.0, CD=5.0, cepdu2=(0.1)**2
56!
57! Variables locales :
58      INTEGER :: i
59      REAL, dimension(klon) :: zdu2, zdphi, ztsolv, ztvd
60      REAL, dimension(klon) :: zscf, friv, frih, zucf, zcr
61      REAL, dimension(klon) :: zcfm1, zcfh1
62      REAL, dimension(klon) :: zcfm2, zcfh2
63      REAL, dimension(klon) :: trm0, trm1
64!-------------------------------------------------------------------------
65      REAL :: fsta, fins, x
66      fsta(x) = 1.0 / (1.0+10.0*x*(1+8.0*x))
67      fins(x) = SQRT(1.0-18.0*x)
68!-------------------------------------------------------------------------
69!
70      DO i = 1, knon
71!
72       zdphi(i) = zgeop(i)
73       zdu2(i) = max(cepdu2,speed(i)**2)
74       pref(i) = exp(log(psol(i)) - zdphi(i)/(RD*t(i)* &
75                 (1.+ RETV * max(q(i),0.0))))
76       ztsolv(i) = ts(i)
77       trm0(i) = 1. + RETV * max(qsurf(i),0.0)
78       trm1(i) = 1. + RETV * max(q(i),0.0)
79       ztsolv(i) = ztsolv(i) * trm0(i)
80       ztvd(i) = (t(i)+zdphi(i)/RCPD/(1.+RVTMP2*q(i))) &
81          *(1.+RETV*q(i))
82       zri1(i) = zdphi(i)*(ztvd(i)-ztsolv(i))/(zdu2(i)*ztvd(i))
83!
84! on teste zri1 par rapport au Richardson de la 1ere couche ri1
85!
86!IM +++
87       IF(1.EQ.0) THEN
88       IF (okri) THEN
89         IF (ri1(i).GE.0.0.AND.zri1(i).LT.0.0) THEN
90           zri1(i) = ri1(i)
91         ELSE IF(ri1(i).LT.0.0.AND.zri1(i).GE.0.0) THEN
92           zri1(i) = ri1(i)
93         ENDIF
94       ENDIF
95       ENDIF
96!IM ---
97!
98       cdran(i) = (RKAR/log(1.+zdphi(i)/(RG*rugos(i))))**2
99
100       IF (zri1(i) .ge. 0.) THEN
101!
102! situation stable : pour eviter les inconsistances dans les cas
103! tres stables on limite zri1 a 20. cf Hess et al. (1995)
104!
105         zri1(i) = min(20.,zri1(i))
106!
107         IF (.NOT.zxli) THEN
108           zscf(i) = SQRT(1.+CD*ABS(zri1(i)))
109           friv(i) = max(1. / (1.+2.*CB*zri1(i)/ zscf(i)), 0.1)
110           zcfm1(i) = cdran(i) * friv(i)
111           frih(i) = max(1./ (1.+3.*CB*zri1(i)*zscf(i)), 0.1 )
112           zcfh1(i) = cdran(i) * frih(i)
113           cdram(i) = zcfm1(i)
114           cdrah(i) = zcfh1(i)
115         ELSE
116           cdram(i) = cdran(i)* fsta(zri1(i))
117           cdrah(i) = cdran(i)* fsta(zri1(i))
118         ENDIF
119!
120       ELSE
121!
122! situation instable
123!
124         IF (.NOT.zxli) THEN
125           zucf(i) = 1./(1.+3.0*CB*CC*cdran(i)*SQRT(ABS(zri1(i)) &
126                 *(1.0+zdphi(i)/(RG*rugos(i)))))
127           zcfm2(i) = cdran(i)*max((1.-2.0*CB*zri1(i)*zucf(i)),0.1)
128           zcfh2(i) = cdran(i)*max((1.-3.0*CB*zri1(i)*zucf(i)),0.1)
129           cdram(i) = zcfm2(i)
130           cdrah(i) = zcfh2(i)
131         ELSE
132           cdram(i) = cdran(i)* fins(zri1(i))
133           cdrah(i) = cdran(i)* fins(zri1(i))
134         ENDIF
135!
136! cdrah sur l'ocean cf. Miller et al. (1992)
137!
138         zcr(i) = (0.0016/(cdran(i)*SQRT(zdu2(i))))*ABS(ztvd(i)-ztsolv(i)) &
139               **(1./3.)
140         IF (nsrf.EQ.is_oce) cdrah(i) = cdran(i)*(1.0+zcr(i)**1.25) &
141                  **(1./1.25)
142       ENDIF
143!
144      END DO
145      RETURN
146      END SUBROUTINE coefcdrag
Note: See TracBrowser for help on using the repository browser.