source: LMDZ6/trunk/libf/phylmd/cdrag.F90 @ 3513

Last change on this file since 3513 was 2868, checked in by fhourdin, 8 years ago

Option pour les fonctions de stabilité pour les coefficients
d'échanges en surface (Etienne Vignon)

  • Property svn:keywords set to Id
File size: 12.1 KB
Line 
1!
2!$Id: cdrag.F90 2868 2017-05-02 11:16:06Z fhourdin $
3!
4 SUBROUTINE cdrag(knon,  nsrf,   &
5     speed, t1,    q1,    zgeop1, &
6     psol,  tsurf, qsurf, z0m, z0h,  &
7     cdm,  cdh,  zri,   pref)
8
9  USE dimphy
10  USE indice_sol_mod
11  USE print_control_mod, ONLY: lunout, prt_level
12  USE ioipsl_getin_p_mod, ONLY : getin_p
13
14  IMPLICIT NONE
15! ================================================================= c
16!
17! Objet : calcul des cdrags pour le moment (pcfm) et
18!         les flux de chaleur sensible et latente (pcfh) d'apr??s
19!         Louis 1982, Louis 1979, King et al 2001
20!         ou Zilitinkevich et al 2002  pour les cas stables, Louis 1979
21!         et 1982 pour les cas instables
22!
23! Modified history:
24!  writting on the 20/05/2016
25!  modified on the 13/12/2016 to be adapted to LMDZ6
26!
27! References:
28!   Louis, J. F., 1979: A parametric model of vertical eddy fluxes in the
29!     atmosphere. Boundary-Layer Meteorology. 01/1979; 17(2):187-202.
30!   Louis, J. F., Tiedtke, M. and Geleyn, J. F., 1982: `A short history of the
31!     operational PBL parametrization at ECMWF'. Workshop on boundary layer
32!     parametrization, November 1981, ECMWF, Reading, England.
33!     Page: 19. Equations in Table 1.
34!   Mascart P, Noilhan J, Giordani H 1995.A MODIFIED PARAMETERIZATION OF FLUX-PROFILE RELATIONSHIPS
35!    IN THE SURFACE LAYER USING DIFFERENT ROUGHNESS LENGTH VALUES FOR HEAT AND MOMENTUM
36!    Boundary-Layer Meteorology 72: 331-344
37!   Anton Beljaars. May 1992. The parametrization of the planetary boundary layer. 
38!     European Centre for Medium-Range Weather Forecasts.
39!     Equations: 110-113. Page 40.
40!   Miller,M.J., A.C.M.Beljaars, T.N.Palmer. 1992. The sensitivity of the ECMWF
41!     model to the parameterization of evaporation from the tropical oceans. J.
42!     Climate, 5:418-434.
43!   King J.C, Connolley, W.M ad Derbyshire S.H. 2001, Sensitivity of Modelled Antarctic climate
44!   to surface and boundary-layer flux parametrizations
45!   QJRMS, 127, pp 779-794
46!
47! ================================================================= c
48! ================================================================= c
49! On choisit le couple de fonctions de correction avec deux flags:
50! Un pour les cas instables, un autre pour les cas stables
51!
52! iflag_corr_insta:
53!                   1: Louis 1979 avec les modifications de Mascart 1995 (z0/= z0h)
54!                   2: Louis 1982
55!                   3: Laurent Li
56!
57! iflag_corr_sta:
58!                   1: Louis 1979 avec les modifications de Mascart 1995 (z0/= z0h)
59!                   2: Louis 1982
60!                   3: Laurent Li
61!                   4: King  2001 (SHARP)
62!                   5: MO 1st order theory (allow collapse of turbulence)
63!           
64!
65!*****************************************************************
66! Parametres d'entree
67!*****************************************************************
68 
69  INTEGER, INTENT(IN)                      :: knon, nsrf ! nombre de points de grille sur l'horizontal + type de surface
70  REAL, DIMENSION(klon), INTENT(IN)        :: speed ! module du vent au 1er niveau du modele
71  REAL, DIMENSION(klon), INTENT(IN)        :: zgeop1! geopotentiel au 1er niveau du modele
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  REAL, DIMENSION(klon), INTENT(IN)        :: t1  ! Temperature au premier niveau (K)
76  REAL, DIMENSION(klon), INTENT(IN)        :: q1  ! humidite specifique au premier niveau (kg/kg)
77  REAL, DIMENSION(klon), INTENT(IN)        :: psol ! pression au sol
78
79
80
81! Parametres de sortie
82!******************************************************************
83  REAL, DIMENSION(klon), INTENT(OUT)       :: cdm  ! Drag coefficient for heat flux
84  REAL, DIMENSION(klon), INTENT(OUT)       :: cdh  ! 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! Variables Locales
89!******************************************************************
90 
91
92  INCLUDE "YOMCST.h"
93  INCLUDE "YOETHF.h"
94  INCLUDE "clesphys.h"
95
96
97  REAL, PARAMETER       :: CKAP=0.40, CKAPT=0.42
98  REAL                   CEPDU2
99  REAL                   ALPHA
100  REAL                   CB,CC,CD,C2,C3
101  REAL                   MU, CM, CH, B, CMstar, CHstar
102  REAL                   PM, PH, BPRIME
103  REAL                   C
104  INTEGER                ng_q1                      ! Number of grids that q1 < 0.0
105  INTEGER                ng_qsurf                   ! Number of grids that qsurf < 0.0
106  INTEGER                i
107  REAL                   zdu2, ztsolv
108  REAL                   ztvd, zscf
109  REAL                   zucf, zcr
110  REAL                   friv, frih
111  REAL, DIMENSION(klon) :: FM, FH                   ! stability functions
112  REAL, DIMENSION(klon) :: cdmn, cdhn               ! Drag coefficient in neutral conditions
113  REAL zzzcd
114
115  LOGICAL, SAVE :: firstcall = .TRUE.
116  !$OMP THREADPRIVATE(firstcall)
117  INTEGER, SAVE :: iflag_corr_sta
118  !$OMP THREADPRIVATE(iflag_corr_sta)
119  INTEGER, SAVE :: iflag_corr_insta
120  !$OMP THREADPRIVATE(iflag_corr_insta)
121
122!===================================================================c
123! Valeurs numeriques des constantes
124!===================================================================c
125
126
127! Minimum du carre du vent
128
129 CEPDU2 = (0.1)**2
130
131! Louis 1982
132
133  CB=5.0
134  CC=5.0
135  CD=5.0
136
137
138! King 2001
139
140  C2=0.25
141  C3=0.0625
142
143 
144! Louis 1979
145
146  BPRIME=4.7
147  B=9.4
148 
149
150!MO
151
152  ALPHA=5.0
153 
154
155! ================================================================= c
156! Tests avant de commencer
157! Fuxing WANG, 04/03/2015
158! To check if there are negative q1, qsurf values.
159!====================================================================c
160  ng_q1 = 0     ! Initialization
161  ng_qsurf = 0  ! Initialization
162  DO i = 1, knon
163     IF (q1(i).LT.0.0)     ng_q1 = ng_q1 + 1
164     IF (qsurf(i).LT.0.0)  ng_qsurf = ng_qsurf + 1
165  ENDDO
166  IF (ng_q1.GT.0 .and. prt_level > 5) THEN
167      WRITE(lunout,*)" *** Warning: Negative q1(humidity at 1st level) values in cdrag.F90 !"
168      WRITE(lunout,*)" The total number of the grids is: ", ng_q1
169      WRITE(lunout,*)" The negative q1 is set to zero "
170!      abort_message="voir ci-dessus"
171!      CALL abort_physic(modname,abort_message,1)
172  ENDIF
173  IF (ng_qsurf.GT.0 .and. prt_level > 5) THEN
174      WRITE(lunout,*)" *** Warning: Negative qsurf(humidity at surface) values in cdrag.F90 !"
175      WRITE(lunout,*)" The total number of the grids is: ", ng_qsurf
176      WRITE(lunout,*)" The negative qsurf is set to zero "
177!      abort_message="voir ci-dessus"
178!      CALL abort_physic(modname,abort_message,1)
179  ENDIF
180
181
182
183!=============================================================================c
184! Calcul du cdrag
185!=============================================================================c
186
187! On choisit les fonctions de stabilite utilisees au premier appel
188!**************************************************************************
189  IF (firstcall) THEN
190   iflag_corr_sta=2
191   iflag_corr_insta=2
192 
193   CALL getin_p('iflag_corr_sta',iflag_corr_sta)
194   CALL getin_p('iflag_corr_insta',iflag_corr_insta)
195
196   firstcall = .FALSE.
197 ENDIF
198
199!xxxxxxxxxxxxxxxxxxxxxxx
200  DO i = 1, knon        ! Boucle sur l'horizontal
201!xxxxxxxxxxxxxxxxxxxxxxx
202
203
204! calculs preliminaires:
205!***********************
206     
207
208     zdu2 = MAX(CEPDU2, speed(i)**2)
209     pref(i) = EXP(LOG(psol(i)) - zgeop1(i)/(RD*t1(i)* &
210                 (1.+ RETV * max(q1(i),0.0))))           ! negative q1 set to zero
211     ztsolv = tsurf(i) * (1.0+RETV*max(qsurf(i),0.0))    ! negative qsurf set to zero
212     ztvd = (t1(i)+zgeop1(i)/RCPD/(1.+RVTMP2*q1(i))) &
213          *(1.+RETV*max(q1(i),0.0))                      ! negative q1 set to zero
214     zri(i) = zgeop1(i)*(ztvd-ztsolv)/(zdu2*ztvd)
215
216
217! Coefficients CD neutres : k^2/ln(z/z0) et k^2/(ln(z/z0)*ln(z/z0h)):
218!********************************************************************
219
220     zzzcd=CKAP/LOG(1.+zgeop1(i)/(RG*z0m(i)))
221     cdmn(i) = zzzcd*zzzcd
222     cdhn(i) = zzzcd*(CKAP/LOG(1.+zgeop1(i)/(RG*z0h(i))))
223
224
225! Calcul des fonctions de stabilit?? FMs, FHs, FMi, FHi :
226!*******************************************************
227
228!''''''''''''''
229! Cas instables
230!''''''''''''''
231
232 IF (zri(i) .LT. 0.) THEN   
233
234
235        SELECT CASE (iflag_corr_insta)
236   
237        CASE (1) ! Louis 1979 + Mascart 1995
238
239           MU=LOG(MAX(z0m(i)/z0h(i),0.01))
240           CMstar=6.8741+2.6933*MU-0.3601*(MU**2)+0.0154*(MU**3)
241           PM=0.5233-0.0815*MU+0.0135*(MU**2)-0.001*(MU**3)
242           CHstar=3.2165+4.3431*MU+0.536*(MU**2)-0.0781*(MU**3)
243           PH=0.5802-0.1571*MU+0.0327*(MU**2)-0.0026*(MU**3)
244           CH=CHstar*B*CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) &
245            & * CKAPT/LOG(z0h(i)+zgeop1(i)/(RG*z0h(i)))       &
246            & * ((zgeop1(i)/(RG*z0h(i)))**PH)
247           CM=CMstar*B*CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) &
248            & *CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i)))         &
249            & * ((zgeop1(i)/(RG*z0m(i)))**PM)
250         
251
252
253
254           FM(i)=1.-B*zri(i)/(1.+CM*SQRT(ABS(zri(i))))
255           FH(i)=1.-B*zri(i)/(1.+CH*SQRT(ABS(zri(i))))
256
257        CASE (2) ! Louis 1982
258
259           zucf = 1./(1.+3.0*CB*CC*cdmn(i)*SQRT(ABS(zri(i)) &
260                *(1.0+zgeop1(i)/(RG*z0m(i)))))
261           FM(i) = AMAX1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min)
262           FH(i) = AMAX1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min)
263
264
265        CASE (3) ! Laurent Li
266
267           
268           FM(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min)
269           FH(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min)
270
271
272
273         CASE default ! Louis 1982
274
275           zucf = 1./(1.+3.0*CB*CC*cdmn(i)*SQRT(ABS(zri(i)) &
276                *(1.0+zgeop1(i)/(RG*z0m(i)))))
277           FM(i) = AMAX1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min)
278           FH(i) = AMAX1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min)
279
280
281         END SELECT
282
283
284
285! Calcul des drags
286
287
288       cdm(i)=cdmn(i)*FM(i)
289       cdh(i)=f_cdrag_ter*cdhn(i)*FH(i)
290
291
292! Traitement particulier des cas oceaniques
293! on applique Miller et al 1992 en l'absence de gustiness
294
295  IF (nsrf == is_oce) THEN
296!        cdh(i)=f_cdrag_oce*cdhn(i)*FH(i)
297
298        IF(iflag_gusts==0) THEN
299           zcr = (0.0016/(cdmn(i)*SQRT(zdu2)))*ABS(ztvd-ztsolv)**(1./3.)
300           cdh(i) =f_cdrag_oce* cdhn(i)*(1.0+zcr**1.25)**(1./1.25)
301        ENDIF
302
303
304        cdm(i)=MIN(cdm(i),cdmmax)
305        cdh(i)=MIN(cdh(i),cdhmax)
306
307  END IF
308
309
310
311 ELSE                         
312
313!'''''''''''''''
314! Cas stables :
315!'''''''''''''''
316        zri(i) = MIN(20.,zri(i))
317
318       SELECT CASE (iflag_corr_sta)
319   
320        CASE (1) ! Louis 1979 + Mascart 1995
321   
322           FM(i)=MAX(1./((1+BPRIME*zri(i))**2),f_ri_cd_min)
323           FH(i)=FM(i)
324
325
326        CASE (2) ! Louis 1982
327           
328           zscf = SQRT(1.+CD*ABS(zri(i)))
329           FM(i)= AMAX1(1. / (1.+2.*CB*zri(i)/zscf), f_ri_cd_min)
330           FH(i)= AMAX1(1./ (1.+3.*CB*zri(i)*zscf), f_ri_cd_min )
331
332
333        CASE (3) ! Laurent Li
334   
335           FM(i)=MAX(1.0 / (1.0+10.0*zri(i)*(1+8.0*zri(i))),f_ri_cd_min)
336           FH(i)=FM(i)
337
338
339        CASE (4)  ! King 2001
340         
341           if (zri(i) .LT. C2/2.) then
342           FM(i)=MAX((1.-zri(i)/C2)**2,f_ri_cd_min)
343           FH(i)=  FM(i)
344
345
346           else
347           FM(i)=MAX(C3*((C2/zri(i))**2),f_ri_cd_min)
348           FH(i)= FM(i)
349           endif
350
351
352        CASE (5) ! MO
353   
354          if (zri(i) .LT. 1./alpha) then
355
356           FM(i)=MAX((1.-alpha*zri(i))**2,f_ri_cd_min)
357           FH(i)=FM(i)
358
359           else
360
361
362           FM(i)=MAX(1E-7,f_ri_cd_min)
363           FH(i)=FM(i)
364
365          endif
366
367
368
369
370
371        CASE default ! Louis 1982
372
373           zscf = SQRT(1.+CD*ABS(zri(i)))
374           FM(i)= AMAX1(1. / (1.+2.*CB*zri(i)/zscf), f_ri_cd_min)
375           FH(i)= AMAX1(1./ (1.+3.*CB*zri(i)*zscf), f_ri_cd_min )
376
377
378
379   END SELECT
380
381! Calcul des drags
382
383
384       cdm(i)=cdmn(i)*FM(i)
385       cdh(i)=f_cdrag_ter*cdhn(i)*FH(i)
386
387       IF(nsrf.EQ.is_oce) THEN
388
389        cdh(i)=f_cdrag_oce*cdhn(i)*FH(i)
390        cdm(i)=MIN(cdm(i),cdmmax)
391        cdh(i)=MIN(cdh(i),cdhmax)
392
393      ENDIF
394
395
396 ENDIF
397
398
399
400
401!xxxxxxxxxxx
402  END DO   !  Fin de la boucle sur l'horizontal
403!xxxxxxxxxxx
404! ================================================================= c
405     
406 
407
408END SUBROUTINE cdrag
409
410
411!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Note: See TracBrowser for help on using the repository browser.