source: LMDZ5/branches/testing/libf/phylmd/coefcdrag.F90 @ 3046

Last change on this file since 3046 was 2408, checked in by Laurent Fairhead, 9 years ago

Merged trunk changes r2298:2396 into testing branch

  • 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: 6.0 KB
RevLine 
[524]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)
[1795]8
9      USE indice_sol_mod
10
[524]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
[1795]23! nsrf----input-I- indice pour le type de surface; voir indice_sol_mod.F90
[524]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!
[793]52      include "YOMCST.h"
53      include "YOETHF.h"
[2056]54      INCLUDE "clesphys.h"
[524]55! Quelques constantes :
[1061]56      REAL, parameter :: RKAR=0.40, CB=5.0, CC=5.0, CD=5.0, cepdu2=(0.1)**2
[524]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
[2160]65
[2258]66  CHARACTER (LEN=80) :: abort_message
67  CHARACTER (LEN=20) :: modname = 'coefcdra'
68
69
70!
71
72
[524]73!-------------------------------------------------------------------------
74      REAL :: fsta, fins, x
75      fsta(x) = 1.0 / (1.0+10.0*x*(1+8.0*x))
76      fins(x) = SQRT(1.0-18.0*x)
77!-------------------------------------------------------------------------
[2258]78
79  abort_message='obsolete, remplace par cdrag, use at you own risk'
[2408]80  CALL abort_physic(modname,abort_message,1)
[2258]81
[524]82!
83      DO i = 1, knon
84!
85       zdphi(i) = zgeop(i)
[1061]86       zdu2(i) = max(cepdu2,speed(i)**2)
[524]87       pref(i) = exp(log(psol(i)) - zdphi(i)/(RD*t(i)* &
88                 (1.+ RETV * max(q(i),0.0))))
89       ztsolv(i) = ts(i)
[2056]90!       ztvd(i) = t(i) * (psol(i)/pref(i))**RKAPPA
[2160]91!       ztvd(i) = (t(i)+zdphi(i)/RCPD/(1.+RVTMP2*q(i))) &
92!          *(1.+RETV*q(i))
93       ztvd(i) = (t(i)+zdphi(i)/RCPD/(1.+RVTMP2*q(i)))
[524]94       trm0(i) = 1. + RETV * max(qsurf(i),0.0)
95       trm1(i) = 1. + RETV * max(q(i),0.0)
96       ztsolv(i) = ztsolv(i) * trm0(i)
[2160]97       ztvd(i) = ztvd(i) * trm1(i)
[524]98       zri1(i) = zdphi(i)*(ztvd(i)-ztsolv(i))/(zdu2(i)*ztvd(i))
99!
100! on teste zri1 par rapport au Richardson de la 1ere couche ri1
101!
102!IM +++
103       IF(1.EQ.0) THEN
104       IF (okri) THEN
105         IF (ri1(i).GE.0.0.AND.zri1(i).LT.0.0) THEN
106           zri1(i) = ri1(i)
107         ELSE IF(ri1(i).LT.0.0.AND.zri1(i).GE.0.0) THEN
108           zri1(i) = ri1(i)
109         ENDIF
110       ENDIF
111       ENDIF
112!IM ---
113!
114       cdran(i) = (RKAR/log(1.+zdphi(i)/(RG*rugos(i))))**2
115
116       IF (zri1(i) .ge. 0.) THEN
117!
118! situation stable : pour eviter les inconsistances dans les cas
119! tres stables on limite zri1 a 20. cf Hess et al. (1995)
120!
121         zri1(i) = min(20.,zri1(i))
122!
123         IF (.NOT.zxli) THEN
124           zscf(i) = SQRT(1.+CD*ABS(zri1(i)))
[2160]125           friv(i) = max(1. / (1.+2.*CB*zri1(i)/ zscf(i)), f_ri_cd_min)
[524]126           zcfm1(i) = cdran(i) * friv(i)
[2160]127           frih(i) = max(1./ (1.+3.*CB*zri1(i)*zscf(i)), f_ri_cd_min )
[2056]128!           zcfh1(i) = cdran(i) * frih(i)
129           zcfh1(i) = f_cdrag_ter*cdran(i) * frih(i)
130           IF(nsrf.EQ.is_oce) zcfh1(i)=f_cdrag_oce*cdran(i)*frih(i)
[524]131           cdram(i) = zcfm1(i)
132           cdrah(i) = zcfh1(i)
133         ELSE
134           cdram(i) = cdran(i)* fsta(zri1(i))
135           cdrah(i) = cdran(i)* fsta(zri1(i))
136         ENDIF
137!
138       ELSE
139!
140! situation instable
141!
142         IF (.NOT.zxli) THEN
143           zucf(i) = 1./(1.+3.0*CB*CC*cdran(i)*SQRT(ABS(zri1(i)) &
144                 *(1.0+zdphi(i)/(RG*rugos(i)))))
[2160]145           zcfm2(i) = cdran(i)*max((1.-2.0*CB*zri1(i)*zucf(i)),f_ri_cd_min)
146!           zcfh2(i) = cdran(i)*max((1.-3.0*CB*zri1(i)*zucf(i)),f_ri_cd_min)
147           zcfh2(i) = f_cdrag_ter*cdran(i)*max((1.-3.0*CB*zri1(i)*zucf(i)),f_ri_cd_min)
[524]148           cdram(i) = zcfm2(i)
149           cdrah(i) = zcfh2(i)
150         ELSE
151           cdram(i) = cdran(i)* fins(zri1(i))
152           cdrah(i) = cdran(i)* fins(zri1(i))
153         ENDIF
154!
155! cdrah sur l'ocean cf. Miller et al. (1992)
156!
157         zcr(i) = (0.0016/(cdran(i)*SQRT(zdu2(i))))*ABS(ztvd(i)-ztsolv(i)) &
158               **(1./3.)
[2056]159!         IF (nsrf.EQ.is_oce) cdrah(i) = cdran(i)*(1.0+zcr(i)**1.25) &
160!                  **(1./1.25)
161         IF (nsrf.EQ.is_oce) cdrah(i)=f_cdrag_oce*cdran(i)*(1.0+zcr(i)**1.25) &
[524]162                  **(1./1.25)
163       ENDIF
164!
165      END DO
166      RETURN
167      END SUBROUTINE coefcdrag
Note: See TracBrowser for help on using the repository browser.