source: LMDZ5/trunk/libf/phylmd/cdrag.F90 @ 2232

Last change on this file since 2232 was 2232, checked in by fhourdin, 9 years ago

1) Fusion des procédures clcdrag.F90 et coefcdrag.F90

dans une procédure unique cdrag.F90.
Les deux anciennes sont obsolètes mais maintenues dans le
code quelques mois pour d'éventuels tests.

2) modification du nom du module arth.F90 en arth_m.F90 pour

avoir le même nom de fichier et de module (règle adoptée
pour LMDZ et utilisée par makelmdz).

1) merging clcdrag.F90 and coefcdrag.F90 in one procedure cdrag.F90
2) renaming arth.F90 in arth_m.F90

  1. Wang/FH
File size: 9.6 KB
Line 
1!
2!$Id: cdrag.F90 ???? 2015-??-?? ??:??:??Z ? $
3!
4 SUBROUTINE cdrag( knon,  nsrf,   &
5     speed, t1,    q1,    zgeop1, &
6     psol,  tsurf, qsurf, rugos,  &
7     pcfm,  pcfh,  zri,   pref )
8
9  USE dimphy
10  USE indice_sol_mod
11  IMPLICIT NONE
12! ================================================================= c
13!
14! Objet : calcul des cdrags pour le moment (pcfm) et
15!         les flux de chaleur sensible et latente (pcfh).   
16!
17! Modified histroy:
18!   27-Jan-2014: richardson number inconsistant between
19!   coefcdrag.F90 and clcdrag.F90, Fuxing WANG wrote this subroutine
20!   by merging coefcdrag and clcdrag.
21!
22! References:
23!   Louis, J. F., 1979: A parametric model of vertical eddy fluxes in the
24!     atmosphere. Boundary-Layer Meteorology. 01/1979; 17(2):187-202.
25!   Louis, J. F., Tiedtke, M. and Geleyn, J. F., 1982: `A short history of the
26!     operational PBL parametrization at ECMWF'. Workshop on boundary layer
27!     parametrization, November 1981, ECMWF, Reading, England.
28!     Page: 19. Equations in Table 1.
29!   Anton Beljaars. May 1992. The parametrization of the planetary boundary layer. 
30!     European Centre for Medium-Range Weather Forecasts.
31!     Equations: 110-113. Page 40.
32!   Miller,M.J., A.C.M.Beljaars, T.N.Palmer. 1992. The sensitivity of the ECMWF
33!     model to the parameterization of evaporation from the tropical oceans. J.
34!     Climate, 5:418-434.
35!
36! ================================================================= c
37!
38! knon----input-I- nombre de points pour un type de surface
39! nsrf----input-I- indice pour le type de surface; voir indicesol.h
40! speed---input-R- module du vent au 1er niveau du modele
41! t1------input-R- temperature de l'air au 1er niveau du modele
42! q1------input-R- humidite de l'air au 1er niveau du modele
43! zgeop---input-R- geopotentiel au 1er niveau du modele
44! psol----input-R- pression au sol
45! tsurf---input-R- temperature de l'air a la surface
46! qsurf---input-R- humidite de l'air a la surface
47! rugos---input-R- rugosite
48!! u1, v1 are removed, speed is used. Fuxing WANG, 04/03/2015,
49!! u1------input-R- vent zonal au 1er niveau du modele
50!! v1------input-R- vent meridien au 1er niveau du modele
51!
52! pcfm---output-R- cdrag pour le moment
53! pcfh---output-R- cdrag pour les flux de chaleur latente et sensible
54! zri----output-R- Richardson number
55! pref---output-R- pression au niveau zgeop/RG
56!
57! Parameters:
58! ckap-----Karman constant
59! cb,cc,cd-parameters in Louis et al., 1982
60! ================================================================= c
61!
62!
63! Parametres d'entree
64!*****************************************************************
65  INTEGER, INTENT(IN)                      :: knon, nsrf
66  REAL, DIMENSION(klon), INTENT(IN)        :: speed ! module du vent au 1er niveau du modele
67  REAL, DIMENSION(klon), INTENT(IN)        :: zgeop1! geopotentiel au 1er niveau du modele
68  REAL, DIMENSION(klon), INTENT(IN)        :: psol  ! pression au sol
69  REAL, DIMENSION(klon), INTENT(IN)        :: t1    ! temperature at 1st level
70  REAL, DIMENSION(klon), INTENT(IN)        :: q1    ! humidity at 1st level
71  REAL, DIMENSION(klon), INTENT(IN)        :: tsurf ! Surface temperature (K)
72  REAL, DIMENSION(klon), INTENT(IN)        :: qsurf ! Surface humidity (Kg/Kg)
73  REAL, DIMENSION(klon), INTENT(IN)        :: rugos ! Rugosity at surface (m)
74!  paprs, pplay u1, v1: to be deleted
75!  they were in the old clcdrag. Fuxing WANG, 04/03/2015
76!  REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs
77!  REAL, DIMENSION(klon,klev), INTENT(IN)   :: pplay
78!  REAL, DIMENSION(klon), INTENT(IN)        :: u1, v1
79
80! Parametres de sortie
81!******************************************************************
82  REAL, DIMENSION(klon), INTENT(OUT)       :: pcfm  ! Drag coefficient for heat flux
83  REAL, DIMENSION(klon), INTENT(OUT)       :: pcfh  ! Drag coefficient for momentum
84  REAL, DIMENSION(klon), INTENT(OUT)       :: zri   ! Richardson number
85  REAL, DIMENSION(klon), INTENT(OUT)       :: pref  ! Pression au niveau zgeop/RG
86
87! Parametres local
88  INTEGER                                  :: ng_q1    ! Number of grids that q1 < 0.0
89  INTEGER                                  :: ng_qsurf ! Number of grids that qsurf < 0.0
90!  zgeop1, psol: to be deleted, they are inputs now. Fuxing WANG, 04/03/2015
91!  REAL, DIMENSION(klon)                    :: zgeop1! geopotentiel au 1er niveau du modele
92!  REAL, DIMENSION(klon)                    :: psol  ! pression au sol
93!
94! ================================================================= c
95!
96  INCLUDE "YOMCST.h"
97  INCLUDE "YOETHF.h"
98!  INCLUDE "indicesol.h"
99  INCLUDE "clesphys.h"
100  INCLUDE "iniprint.h"
101!
102! Quelques constantes et options:
103!!$PB      REAL, PARAMETER :: ckap=0.35, cb=5.0, cc=5.0, cd=5.0, cepdu2=(0.1)**2
104  REAL, PARAMETER :: CKAP=0.40, CB=5.0, CC=5.0, CD=5.0, CEPDU2 = (0.1)**2
105!
106! Variables locales :
107  INTEGER               :: i
108  REAL                  :: zdu2, ztsolv
109  REAL                  :: ztvd, zscf
110  REAL                  :: zucf, zcr
111  REAL                  :: friv, frih
112  REAL, DIMENSION(klon) :: zcfm1, zcfm2 ! Drag coefficient for momentum
113  REAL, DIMENSION(klon) :: zcfh1, zcfh2 ! Drag coefficient for heat flux
114  LOGICAL, PARAMETER    :: zxli=.FALSE. ! calcul des cdrags selon Laurent Li
115  REAL, DIMENSION(klon) :: zcdn         ! Drag coefficient in neutral conditions
116!
117! Fonctions thermodynamiques et fonctions d'instabilite
118  REAL                  :: fsta, fins, x
119  fsta(x) = 1.0 / (1.0+10.0*x*(1+8.0*x))
120  fins(x) = SQRT(1.0-18.0*x)
121
122! ================================================================= c
123! Fuxing WANG, 04/03/2015, delete the calculation of zgeop1
124! (le geopotentiel du premier couche de modele).
125! zgeop1 is an input ivariable in this subroutine.
126!  DO i = 1, knon
127!     zgeop1(i) = RD * t1(i) / (0.5*(paprs(i,1)+pplay(i,1))) &
128!          * (paprs(i,1)-pplay(i,1))
129!  END DO
130! ================================================================= c
131!
132! Fuxing WANG, 04/03/2015
133! To check if there are negative q1, qsurf values.
134  ng_q1 = 0     ! Initialization
135  ng_qsurf = 0  ! Initialization
136  DO i = 1, knon
137     IF (q1(i).LT.0.0)     ng_q1 = ng_q1 + 1
138     IF (qsurf(i).LT.0.0)  ng_qsurf = ng_qsurf + 1
139  ENDDO
140  IF (ng_q1.GT.0) THEN
141      WRITE(lunout,*)" *** Warning: Negative q1(humidity at 1st level) values in cdrag.F90 !"
142      WRITE(lunout,*)" The total number of the grids is: ", ng_q1
143      WRITE(lunout,*)" The negative q1 is set to zero "
144!      abort_message="voir ci-dessus"
145!      CALL abort_gcm(modname,abort_message,1)
146  ENDIF
147  IF (ng_qsurf.GT.0) THEN
148      WRITE(lunout,*)" *** Warning: Negative qsurf(humidity at surface) values in cdrag.F90 !"
149      WRITE(lunout,*)" The total number of the grids is: ", ng_qsurf
150      WRITE(lunout,*)" The negative qsurf is set to zero "
151!      abort_message="voir ci-dessus"
152!      CALL abort_gcm(modname,abort_message,1)
153  ENDIF
154
155! Calculer le frottement au sol (Cdrag)
156  DO i = 1, knon
157!------------------------------------------------------------
158! u1, v1 are replaced by speed. Fuxing WANG, 04/03/2015,
159!     zdu2 = MAX(CEPDU2,u1(i)**2+v1(i)**2)
160!------------------------------------------------------------
161     zdu2 = MAX(CEPDU2, speed(i)**2)
162!     psol(i) = paprs(i,1)
163     pref(i) = EXP(LOG(psol(i)) - zgeop1(i)/(RD*t1(i)* &
164                 (1.+ RETV * max(q1(i),0.0))))  ! negative q1 set to zero
165!------------ the old calculations in clcdrag----------------
166!     ztsolv = tsurf(i) * (1.0+RETV*qsurf(i))
167!     ztvd = (t1(i)+zgeop1(i)/RCPD/(1.+RVTMP2*q1(i))) &
168!          *(1.+RETV*q1(i))
169!------------------------------------------------------------
170! Fuxing WANG, 04/03/2015, in this revised version,
171! the negative qsurf and q1 are set to zero (as in coefcdrag)
172     ztsolv = tsurf(i) * (1.0+RETV*max(qsurf(i),0.0)) ! negative qsurf set to zero
173     ztvd = (t1(i)+zgeop1(i)/RCPD/(1.+RVTMP2*q1(i))) &
174          *(1.+RETV*max(q1(i),0.0)) ! negative q1 set to zero
175     zri(i) = zgeop1(i)*(ztvd-ztsolv)/(zdu2*ztvd)
176     zcdn(i) = (CKAP/LOG(1.+zgeop1(i)/(RG*rugos(i))))**2
177
178     IF (zri(i) .GT. 0.) THEN      ! situation stable
179        zri(i) = MIN(20.,zri(i))
180        IF (.NOT.zxli) THEN
181           zscf = SQRT(1.+CD*ABS(zri(i)))
182           friv = AMAX1(1. / (1.+2.*CB*zri(i)/ZSCF), f_ri_cd_min)
183           zcfm1(i) = zcdn(i) * friv
184           frih = AMAX1(1./ (1.+3.*CB*zri(i)*ZSCF), f_ri_cd_min )
185!!$ PB     zcfh1(i) = zcdn(i) * frih
186!!$ PB     zcfh1(i) = f_cdrag_stable * zcdn(i) * frih
187           zcfh1(i) = f_cdrag_ter * zcdn(i) * frih
188           IF(nsrf.EQ.is_oce) zcfh1(i) = f_cdrag_oce * zcdn(i) * frih
189!!$ PB
190           pcfm(i) = zcfm1(i)
191           pcfh(i) = zcfh1(i)
192        ELSE
193           pcfm(i) = zcdn(i)* fsta(zri(i))
194           pcfh(i) = zcdn(i)* fsta(zri(i))
195        ENDIF
196     ELSE                          ! situation instable
197        IF (.NOT.zxli) THEN
198           zucf = 1./(1.+3.0*CB*CC*zcdn(i)*SQRT(ABS(zri(i)) &
199                *(1.0+zgeop1(i)/(RG*rugos(i)))))
200           zcfm2(i) = zcdn(i)*amax1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min)
201!!$ PB     zcfh2(i) = zcdn(i)*amax1((1.-3.0*cb*zri(i)*zucf),f_ri_cd_min)
202           zcfh2(i) = f_cdrag_ter*zcdn(i)*amax1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min)
203           pcfm(i) = zcfm2(i)
204           pcfh(i) = zcfh2(i)
205        ELSE
206           pcfm(i) = zcdn(i)* fins(zri(i))
207           pcfh(i) = zcdn(i)* fins(zri(i))
208        ENDIF
209! cdrah sur l'ocean cf. Miller et al. (1992)
210        zcr = (0.0016/(zcdn(i)*SQRT(zdu2)))*ABS(ztvd-ztsolv)**(1./3.)
211        IF(nsrf.EQ.is_oce) pcfh(i) =f_cdrag_oce* zcdn(i)*(1.0+zcr**1.25)**(1./1.25)
212     ENDIF
213  END DO
214
215! ================================================================= c
216     
217  ! IM cf JLD : on seuille cdrag_m et cdrag_h
218  IF (nsrf == is_oce) THEN
219     DO i=1,knon
220        pcfm(i)=MIN(pcfm(i),cdmmax)
221        pcfh(i)=MIN(pcfh(i),cdhmax)
222     END DO
223  END IF
224
225END SUBROUTINE cdrag
Note: See TracBrowser for help on using the repository browser.