LMDZ
cdrag.F90
Go to the documentation of this file.
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 
231 END SUBROUTINE cdrag
iflag_gusts
!$Id ok_orolf LOGICAL ok_limitvrai LOGICAL ok_all_xml INTEGER iflag_ener_conserv REAL solaire RCFC12 RCFC12_act CFC12_ppt!IM ajout CFMIP2 CMIP5 LOGICAL ok_4xCO2atm RCFC12_per CFC12_ppt_per!OM correction du bilan d eau global!OM Correction sur precip KE REAL cvl_corr!OM Fonte calotte dans bilan eau LOGICAL ok_lic_melt!IM simulateur ISCCP INTEGER overlap!IM seuils cdrh REAL cdhmax!IM param stabilite s terres et en dehors REAL f_ri_cd_min!IM MAFo pmagic evap0!Frottement au f_cdrag_oce REAL f_z0qh_oce REAL z0h_seaice INTEGER iflag_gusts
Definition: clesphys.h:46
false
!$Id itapm1 ENDIF!IM on interpole les champs sur les niveaux STD de pression!IM a chaque pas de temps de la physique c!positionnement de l argument logique a false c!pour ne pas recalculer deux fois la meme chose!c!a cet effet un appel a plevel_new a ete deplace c!a la fin de la serie d appels c!la boucle DO nlevSTD a ete internalisee c!dans d ou la creation de cette routine c c!CALL false
Definition: calcul_STDlev.h:26
cdrag
subroutine cdrag(knon, nsrf, speed, t1, q1, zgeop1, psol, tsurf, qsurf, z0m, z0h, pcfm, pcfh, zri, pref)
Definition: cdrag.F90:8
print_control_mod
Definition: print_control_mod.F90:2
dimphy
Definition: dimphy.F90:1
indice_sol_mod::is_oce
integer, parameter is_oce
Definition: indice_sol_mod.F90:3
cdmmax
!$Id ok_orolf LOGICAL ok_limitvrai LOGICAL ok_all_xml INTEGER iflag_ener_conserv REAL solaire RCFC12 RCFC12_act CFC12_ppt!IM ajout CFMIP2 CMIP5 LOGICAL ok_4xCO2atm RCFC12_per CFC12_ppt_per!OM correction du bilan d eau global!OM Correction sur precip KE REAL cvl_corr!OM Fonte calotte dans bilan eau LOGICAL ok_lic_melt!IM simulateur ISCCP INTEGER overlap!IM seuils cdrh REAL cdmmax
Definition: clesphys.h:23
indice_sol_mod
Definition: indice_sol_mod.F90:1
lunout
!$Header!gestion des impressions de sorties et de débogage la sortie standard prt_level COMMON comprint lunout
Definition: iniprint.h:7
rg
real rg
Definition: comcstphy.h:1