source: LMDZ5/branches/LMDZ5V1.0-dev/libf/phylmd/clcdrag.F90 @ 3900

Last change on this file since 3900 was 1279, checked in by Laurent Fairhead, 15 years ago

Merged LMDZ4-dev branch changes r1241:1278 into the trunk
Running trunk and LMDZ4-dev in LMDZOR configuration on local
machine (sequential) and SX8 (4-proc) yields identical results
(restart and restartphy are identical binarily)
Log history from r1241 to r1278 is available by switching to
source:LMDZ4/branches/LMDZ4-dev-20091210

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 4.8 KB
Line 
1!
2!$Id: clcdrag.F90 1279 2009-12-10 09:02:56Z evignon $
3!
4SUBROUTINE clcdrag(knon, nsrf, paprs, pplay,&
5     u1, v1, t1, q1, &
6     tsurf, qsurf, rugos, &
7     pcfm, pcfh)
8
9  USE dimphy
10  IMPLICIT NONE
11! ================================================================= c
12!
13! Objet : calcul des cdrags pour le moment (pcfm) et
14!         les flux de chaleur sensible et latente (pcfh).   
15!
16! ================================================================= c
17!
18! knon----input-I- nombre de points pour un type de surface
19! nsrf----input-I- indice pour le type de surface; voir indicesol.h
20! u1-------input-R- vent zonal au 1er niveau du modele
21! v1-------input-R- vent meridien au 1er niveau du modele
22! t1-------input-R- temperature de l'air au 1er niveau du modele
23! q1-------input-R- humidite de l'air au 1er niveau du modele
24! tsurf------input-R- temperature de l'air a la surface
25! qsurf---input-R- humidite de l'air a la surface
26! rugos---input-R- rugosite
27!
28! pcfm---output-R- cdrag pour le moment
29! pcfh---output-R- cdrag pour les flux de chaleur latente et sensible
30!
31  INTEGER, INTENT(IN)                      :: knon, nsrf
32  REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs
33  REAL, DIMENSION(klon,klev), INTENT(IN)   :: pplay
34  REAL, DIMENSION(klon), INTENT(IN)        :: u1, v1, t1, q1
35  REAL, DIMENSION(klon), INTENT(IN)        :: tsurf, qsurf
36  REAL, DIMENSION(klon), INTENT(IN)        :: rugos
37  REAL, DIMENSION(klon), INTENT(OUT)       :: pcfm, pcfh
38!
39! ================================================================= c
40!
41  INCLUDE "YOMCST.h"
42  INCLUDE "YOETHF.h"
43  INCLUDE "indicesol.h"
44  INCLUDE "clesphys.h"
45!
46! Quelques constantes et options:
47!!$PB      REAL, PARAMETER :: ckap=0.35, cb=5.0, cc=5.0, cd=5.0, cepdu2=(0.1)**2
48  REAL, PARAMETER :: ckap=0.40, cb=5.0, cc=5.0, cd=5.0, cepdu2=(0.1)**2
49!
50! Variables locales :
51  INTEGER               :: i
52  REAL                  :: zdu2, ztsolv
53  REAL                  :: ztvd, zscf
54  REAL                  :: zucf, zcr
55  REAL                  :: friv, frih
56  REAL, DIMENSION(klon) :: zcfm1, zcfm2
57  REAL, DIMENSION(klon) :: zcfh1, zcfh2
58  REAL, DIMENSION(klon) :: zcdn
59  REAL, DIMENSION(klon) :: zri
60  REAL, DIMENSION(klon) :: zgeop1       ! geopotentiel au 1er niveau du modele
61  LOGICAL, PARAMETER    :: zxli=.FALSE. ! calcul des cdrags selon Laurent Li
62!
63! Fonctions thermodynamiques et fonctions d'instabilite
64  REAL                  :: fsta, fins, x
65  fsta(x) = 1.0 / (1.0+10.0*x*(1+8.0*x))
66  fins(x) = SQRT(1.0-18.0*x)
67
68! ================================================================= c
69!
70! Calculer le geopotentiel du premier couche de modele
71!
72  DO i = 1, knon
73     zgeop1(i) = RD * t1(i) / (0.5*(paprs(i,1)+pplay(i,1))) &
74          * (paprs(i,1)-pplay(i,1))
75  END DO
76! ================================================================= c
77!
78! Calculer le frottement au sol (Cdrag)
79!
80  DO i = 1, knon
81     zdu2 = MAX(cepdu2,u1(i)**2+v1(i)**2)
82     ztsolv = tsurf(i) * (1.0+RETV*qsurf(i))
83     ztvd = (t1(i)+zgeop1(i)/RCPD/(1.+RVTMP2*q1(i))) &
84          *(1.+RETV*q1(i))
85     zri(i) = zgeop1(i)*(ztvd-ztsolv)/(zdu2*ztvd)
86     zcdn(i) = (ckap/LOG(1.+zgeop1(i)/(RG*rugos(i))))**2
87
88!!$        IF (zri(i) .ge. 0.) THEN      ! situation stable
89     IF (zri(i) .GT. 0.) THEN      ! situation stable
90        zri(i) = MIN(20.,zri(i))
91        IF (.NOT.zxli) THEN
92           zscf = SQRT(1.+cd*ABS(zri(i)))
93           FRIV = AMAX1(1. / (1.+2.*CB*zri(i)/ZSCF), 0.1)
94           zcfm1(i) = zcdn(i) * FRIV
95           FRIH = AMAX1(1./ (1.+3.*CB*zri(i)*ZSCF), 0.1 )
96!!$  PB          zcfh1(i) = zcdn(i) * FRIH
97!!$ PB           zcfh1(i) = f_cdrag_stable * zcdn(i) * FRIH
98           zcfh1(i) = f_cdrag_ter * zcdn(i) * FRIH
99           IF(nsrf.EQ.is_oce) zcfh1(i) = f_cdrag_oce * zcdn(i) * FRIH
100!!$ PB
101           pcfm(i) = zcfm1(i)
102           pcfh(i) = zcfh1(i)
103        ELSE
104           pcfm(i) = zcdn(i)* fsta(zri(i))
105           pcfh(i) = zcdn(i)* fsta(zri(i))
106        ENDIF
107     ELSE                          ! situation instable
108        IF (.NOT.zxli) THEN
109           zucf = 1./(1.+3.0*cb*cc*zcdn(i)*SQRT(ABS(zri(i)) &
110                *(1.0+zgeop1(i)/(RG*rugos(i)))))
111           zcfm2(i) = zcdn(i)*amax1((1.-2.0*cb*zri(i)*zucf),0.1)
112!!$PB            zcfh2(i) = zcdn(i)*amax1((1.-3.0*cb*zri(i)*zucf),0.1)
113           zcfh2(i) = f_cdrag_ter*zcdn(i)*amax1((1.-3.0*cb*zri(i)*zucf),0.1)
114           pcfm(i) = zcfm2(i)
115           pcfh(i) = zcfh2(i)
116        ELSE
117           pcfm(i) = zcdn(i)* fins(zri(i))
118           pcfh(i) = zcdn(i)* fins(zri(i))
119        ENDIF
120        zcr = (0.0016/(zcdn(i)*SQRT(zdu2)))*ABS(ztvd-ztsolv)**(1./3.)
121        IF(nsrf.EQ.is_oce) pcfh(i) =f_cdrag_oce* zcdn(i)*(1.0+zcr**1.25)**(1./1.25)
122     ENDIF
123  END DO
124
125! ================================================================= c
126     
127  ! IM cf JLD : on seuille cdrag_m et cdrag_h
128  IF (nsrf == is_oce) THEN
129     DO i=1,knon
130        pcfm(i)=MIN(pcfm(i),cdmmax)
131        pcfh(i)=MIN(pcfh(i),cdhmax)
132     END DO
133  END IF
134
135END SUBROUTINE clcdrag
Note: See TracBrowser for help on using the repository browser.