source: LMDZ5/branches/AI-cosp/libf/phylmd/cdrag.F90 @ 5360

Last change on this file since 5360 was 2311, checked in by Ehouarn Millour, 9 years ago

Further modifications to enforce physics/dynamics separation:

  • moved iniprint.h and misc_mod back to dyn3d_common, as these should only be used by dynamics.
  • created print_control_mod in the physics to store flags prt_level, lunout, debug to be local to physics (should be used rather than iniprint.h)
  • created abort_physic.F90 , which does the same job as abort_gcm() did, but should be used instead when in physics.
  • reactivated inifis (turned it into a module, inifis_mod.F90) to initialize physical constants and print_control_mod flags.

EM
File size: 9.9 KB
Line 
1!
2!$Id: cdrag.F90 ???? 2015-??-?? ??:??:??Z ? $
3!
4 SUBROUTINE cdrag( knon,  nsrf,   &
5     speed, t1,    q1,    zgeop1, &
6     psol,  tsurf, qsurf, z0m, z0h,  &
7     pcfm,  pcfh,  zri,   pref )
8
9  USE dimphy
10  USE indice_sol_mod
11  USE print_control_mod, ONLY: lunout
12  IMPLICIT NONE
13! ================================================================= c
14!
15! Objet : calcul des cdrags pour le moment (pcfm) et
16!         les flux de chaleur sensible et latente (pcfh).   
17!
18! Modified histroy:
19!   27-Jan-2014: richardson number inconsistant between
20!   coefcdrag.F90 and clcdrag.F90, Fuxing WANG wrote this subroutine
21!   by merging coefcdrag and clcdrag.
22!
23! References:
24!   Louis, J. F., 1979: A parametric model of vertical eddy fluxes in the
25!     atmosphere. Boundary-Layer Meteorology. 01/1979; 17(2):187-202.
26!   Louis, J. F., Tiedtke, M. and Geleyn, J. F., 1982: `A short history of the
27!     operational PBL parametrization at ECMWF'. Workshop on boundary layer
28!     parametrization, November 1981, ECMWF, Reading, England.
29!     Page: 19. Equations in Table 1.
30!   Anton Beljaars. May 1992. The parametrization of the planetary boundary layer. 
31!     European Centre for Medium-Range Weather Forecasts.
32!     Equations: 110-113. Page 40.
33!   Miller,M.J., A.C.M.Beljaars, T.N.Palmer. 1992. The sensitivity of the ECMWF
34!     model to the parameterization of evaporation from the tropical oceans. J.
35!     Climate, 5:418-434.
36!
37! ================================================================= c
38!
39! knon----input-I- nombre de points pour un type de surface
40! nsrf----input-I- indice pour le type de surface; voir indicesol.h
41! speed---input-R- module du vent au 1er niveau du modele
42! t1------input-R- temperature de l'air au 1er niveau du modele
43! q1------input-R- humidite de l'air au 1er niveau du modele
44! zgeop---input-R- geopotentiel au 1er niveau du modele
45! psol----input-R- pression au sol
46! tsurf---input-R- temperature de l'air a la surface
47! qsurf---input-R- humidite de l'air a la surface
48! z0m, z0h---input-R- rugosite
49!! u1, v1 are removed, speed is used. Fuxing WANG, 04/03/2015,
50!! u1------input-R- vent zonal au 1er niveau du modele
51!! v1------input-R- vent meridien au 1er niveau du modele
52!
53! pcfm---output-R- cdrag pour le moment
54! pcfh---output-R- cdrag pour les flux de chaleur latente et sensible
55! zri----output-R- Richardson number
56! pref---output-R- pression au niveau zgeop/RG
57!
58! Parameters:
59! ckap-----Karman constant
60! cb,cc,cd-parameters in Louis et al., 1982
61! ================================================================= c
62!
63!
64! Parametres d'entree
65!*****************************************************************
66  INTEGER, INTENT(IN)                      :: knon, nsrf
67  REAL, DIMENSION(klon), INTENT(IN)        :: speed ! module du vent au 1er niveau du modele
68  REAL, DIMENSION(klon), INTENT(IN)        :: zgeop1! geopotentiel au 1er niveau du modele
69  REAL, DIMENSION(klon), INTENT(IN)        :: psol  ! pression au sol
70  REAL, DIMENSION(klon), INTENT(IN)        :: t1    ! temperature at 1st level
71  REAL, DIMENSION(klon), INTENT(IN)        :: q1    ! humidity at 1st level
72  REAL, DIMENSION(klon), INTENT(IN)        :: tsurf ! Surface temperature (K)
73  REAL, DIMENSION(klon), INTENT(IN)        :: qsurf ! Surface humidity (Kg/Kg)
74  REAL, DIMENSION(klon), INTENT(IN)        :: z0m, z0h ! Rugosity at surface (m)
75!  paprs, pplay u1, v1: to be deleted
76!  they were in the old clcdrag. Fuxing WANG, 04/03/2015
77!  REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs
78!  REAL, DIMENSION(klon,klev), INTENT(IN)   :: pplay
79!  REAL, DIMENSION(klon), INTENT(IN)        :: u1, v1
80
81! Parametres de sortie
82!******************************************************************
83  REAL, DIMENSION(klon), INTENT(OUT)       :: pcfm  ! Drag coefficient for heat flux
84  REAL, DIMENSION(klon), INTENT(OUT)       :: pcfh  ! Drag coefficient for momentum
85  REAL, DIMENSION(klon), INTENT(OUT)       :: zri   ! Richardson number
86  REAL, DIMENSION(klon), INTENT(OUT)       :: pref  ! Pression au niveau zgeop/RG
87
88! Parametres local
89  INTEGER                                  :: ng_q1    ! Number of grids that q1 < 0.0
90  INTEGER                                  :: ng_qsurf ! Number of grids that qsurf < 0.0
91!  zgeop1, psol: to be deleted, they are inputs now. Fuxing WANG, 04/03/2015
92!  REAL, DIMENSION(klon)                    :: zgeop1! geopotentiel au 1er niveau du modele
93!  REAL, DIMENSION(klon)                    :: psol  ! pression au sol
94!
95! ================================================================= c
96!
97  INCLUDE "YOMCST.h"
98  INCLUDE "YOETHF.h"
99!  INCLUDE "indicesol.h"
100  INCLUDE "clesphys.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_m, zcdn_h         ! 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_physic(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_physic(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
177
178! Coefficients CD neutres pour m et h
179     zcdn_m(i) = (CKAP/LOG(1.+zgeop1(i)/(RG*z0m(i))))**2
180     zcdn_h(i) = (CKAP/LOG(1.+zgeop1(i)/(RG*z0h(i))))**2
181
182     IF (zri(i) .GT. 0.) THEN      ! situation stable
183        zri(i) = MIN(20.,zri(i))
184        IF (.NOT.zxli) THEN
185           zscf = SQRT(1.+CD*ABS(zri(i)))
186           friv = AMAX1(1. / (1.+2.*CB*zri(i)/ZSCF), f_ri_cd_min)
187           zcfm1(i) = zcdn_m(i) * friv
188           frih = AMAX1(1./ (1.+3.*CB*zri(i)*ZSCF), f_ri_cd_min )
189!!$ PB     zcfh1(i) = zcdn(i) * frih
190!!$ PB     zcfh1(i) = f_cdrag_stable * zcdn(i) * frih
191           zcfh1(i) = f_cdrag_ter * zcdn_h(i) * frih
192           IF(nsrf.EQ.is_oce) zcfh1(i) = f_cdrag_oce * zcdn_h(i) * frih
193!!$ PB
194           pcfm(i) = zcfm1(i)
195           pcfh(i) = zcfh1(i)
196        ELSE
197           pcfm(i) = zcdn_m(i)* fsta(zri(i))
198           pcfh(i) = zcdn_h(i)* fsta(zri(i))
199        ENDIF
200     ELSE                          ! situation instable
201        IF (.NOT.zxli) THEN
202           zucf = 1./(1.+3.0*CB*CC*zcdn_m(i)*SQRT(ABS(zri(i)) &
203                *(1.0+zgeop1(i)/(RG*z0m(i)))))
204           zcfm2(i) = zcdn_m(i)*amax1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min)
205!!$ PB     zcfh2(i) = zcdn_h(i)*amax1((1.-3.0*cb*zri(i)*zucf),f_ri_cd_min)
206           zcfh2(i) = f_cdrag_ter*zcdn_h(i)*amax1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min)
207           pcfm(i) = zcfm2(i)
208           pcfh(i) = zcfh2(i)
209        ELSE
210           pcfm(i) = zcdn_m(i)* fins(zri(i))
211           pcfh(i) = zcdn_h(i)* fins(zri(i))
212        ENDIF
213        IF(iflag_gusts==0) THEN
214! cdrah sur l'ocean cf. Miller et al. (1992) - only active when gustiness parameterization is not active
215           zcr = (0.0016/(zcdn_m(i)*SQRT(zdu2)))*ABS(ztvd-ztsolv)**(1./3.)
216           IF(nsrf.EQ.is_oce) pcfh(i) =f_cdrag_oce* zcdn_h(i)*(1.0+zcr**1.25)**(1./1.25)
217        ENDIF
218     ENDIF
219  END DO
220
221! ================================================================= c
222     
223  ! IM cf JLD : on seuille cdrag_m et cdrag_h
224  IF (nsrf == is_oce) THEN
225     DO i=1,knon
226        pcfm(i)=MIN(pcfm(i),cdmmax)
227        pcfh(i)=MIN(pcfh(i),cdhmax)
228     END DO
229  END IF
230
231END SUBROUTINE cdrag
Note: See TracBrowser for help on using the repository browser.