source: LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/cdrag_mod.F90 @ 4688

Last change on this file since 4688 was 3886, checked in by Laurent Fairhead, 4 years ago

By request, re-introduction of two 'bugs' so that running without the iflag_new_t2mq2m
gives exactly the same results as CMIP6 previous simulations
  • Property svn:executable set to *
File size: 24.9 KB
Line 
1!
2MODULE cdrag_mod
3!
4! This module contains some procedures for calculation of the cdrag
5! coefficients for turbulent diffusion at surface
6!
7  IMPLICIT NONE
8
9CONTAINS
10!
11!****************************************************************************************
12!
13!r original routine svn3623
14!
15 SUBROUTINE cdrag(knon,  nsrf,   &
16     speed, t1,    q1,    zgeop1, &
17     psol,  tsurf, qsurf, z0m, z0h,  &
18     cdm,  cdh,  zri,   pref)
19
20  USE dimphy
21  USE indice_sol_mod
22  USE print_control_mod, ONLY: lunout, prt_level
23  USE ioipsl_getin_p_mod, ONLY : getin_p
24
25  IMPLICIT NONE
26! ================================================================= c
27!
28! Objet : calcul des cdrags pour le moment (pcfm) et
29!         les flux de chaleur sensible et latente (pcfh) d'apr??s
30!         Louis 1982, Louis 1979, King et al 2001
31!         ou Zilitinkevich et al 2002  pour les cas stables, Louis 1979
32!         et 1982 pour les cas instables
33!
34! Modified history:
35!  writting on the 20/05/2016
36!  modified on the 13/12/2016 to be adapted to LMDZ6
37!
38! References:
39!   Louis, J. F., 1979: A parametric model of vertical eddy fluxes in the
40!     atmosphere. Boundary-Layer Meteorology. 01/1979; 17(2):187-202.
41!   Louis, J. F., Tiedtke, M. and Geleyn, J. F., 1982: `A short history of the
42!     operational PBL parametrization at ECMWF'. Workshop on boundary layer
43!     parametrization, November 1981, ECMWF, Reading, England.
44!     Page: 19. Equations in Table 1.
45!   Mascart P, Noilhan J, Giordani H 1995.A MODIFIED PARAMETERIZATION OF FLUX-PROFILE RELATIONSHIPS
46!    IN THE SURFACE LAYER USING DIFFERENT ROUGHNESS LENGTH VALUES FOR HEAT AND MOMENTUM
47!    Boundary-Layer Meteorology 72: 331-344
48!   Anton Beljaars. May 1992. The parametrization of the planetary boundary layer. 
49!     European Centre for Medium-Range Weather Forecasts.
50!     Equations: 110-113. Page 40.
51!   Miller,M.J., A.C.M.Beljaars, T.N.Palmer. 1992. The sensitivity of the ECMWF
52!     model to the parameterization of evaporation from the tropical oceans. J.
53!     Climate, 5:418-434.
54!   King J.C, Connolley, W.M ad Derbyshire S.H. 2001, Sensitivity of Modelled Antarctic climate
55!   to surface and boundary-layer flux parametrizations
56!   QJRMS, 127, pp 779-794
57!
58! ================================================================= c
59! ================================================================= c
60! On choisit le couple de fonctions de correction avec deux flags:
61! Un pour les cas instables, un autre pour les cas stables
62!
63! iflag_corr_insta:
64!                   1: Louis 1979 avec les modifications de Mascart 1995 (z0/= z0h)
65!                   2: Louis 1982
66!                   3: Laurent Li
67!
68! iflag_corr_sta:
69!                   1: Louis 1979 avec les modifications de Mascart 1995 (z0/= z0h)
70!                   2: Louis 1982
71!                   3: Laurent Li
72!                   4: King  2001 (SHARP)
73!                   5: MO 1st order theory (allow collapse of turbulence)
74!           
75!
76!*****************************************************************
77! Parametres d'entree
78!*****************************************************************
79 
80  INTEGER, INTENT(IN)                      :: knon, nsrf ! nombre de points de grille sur l'horizontal + type de surface
81  REAL, DIMENSION(klon), INTENT(IN)        :: speed ! module du vent au 1er niveau du modele
82  REAL, DIMENSION(klon), INTENT(IN)        :: zgeop1! geopotentiel au 1er niveau du modele
83  REAL, DIMENSION(klon), INTENT(IN)        :: tsurf ! Surface temperature (K)
84  REAL, DIMENSION(klon), INTENT(IN)        :: qsurf ! Surface humidity (Kg/Kg)
85  REAL, DIMENSION(klon), INTENT(IN)        :: z0m, z0h ! Rugosity at surface (m)
86  REAL, DIMENSION(klon), INTENT(IN)        :: t1  ! Temperature au premier niveau (K)
87  REAL, DIMENSION(klon), INTENT(IN)        :: q1  ! humidite specifique au premier niveau (kg/kg)
88  REAL, DIMENSION(klon), INTENT(IN)        :: psol ! pression au sol
89
90
91
92! Parametres de sortie
93!******************************************************************
94  REAL, DIMENSION(klon), INTENT(OUT)       :: cdm   ! Drag coefficient for momentum
95  REAL, DIMENSION(klon), INTENT(OUT)       :: cdh   ! Drag coefficient for heat flux
96  REAL, DIMENSION(klon), INTENT(OUT)       :: zri   ! Richardson number
97  REAL, DIMENSION(klon), INTENT(OUT)       :: pref  ! Pression au niveau zgeop/RG
98
99! Variables Locales
100!******************************************************************
101 
102
103  INCLUDE "YOMCST.h"
104  INCLUDE "YOETHF.h"
105  INCLUDE "clesphys.h"
106
107
108  REAL, PARAMETER       :: CKAP=0.40, CKAPT=0.42
109  REAL                   CEPDU2
110  REAL                   ALPHA
111  REAL                   CB,CC,CD,C2,C3
112  REAL                   MU, CM, CH, B, CMstar, CHstar
113  REAL                   PM, PH, BPRIME
114  REAL                   C
115  INTEGER                ng_q1                      ! Number of grids that q1 < 0.0
116  INTEGER                ng_qsurf                   ! Number of grids that qsurf < 0.0
117  INTEGER                i
118  REAL                   zdu2, ztsolv
119  REAL                   ztvd, zscf
120  REAL                   zucf, zcr
121  REAL                   friv, frih
122  REAL, DIMENSION(klon) :: FM, FH                   ! stability functions
123  REAL, DIMENSION(klon) :: cdmn, cdhn               ! Drag coefficient in neutral conditions
124  REAL zzzcd
125
126  LOGICAL, SAVE :: firstcall = .TRUE.
127  !$OMP THREADPRIVATE(firstcall)
128  INTEGER, SAVE :: iflag_corr_sta
129  !$OMP THREADPRIVATE(iflag_corr_sta)
130  INTEGER, SAVE :: iflag_corr_insta
131  !$OMP THREADPRIVATE(iflag_corr_insta)
132
133  integer, save :: iflag_new_t2mq2m
134  !$OMP THREADPRIVATE(iflag_new_t2mq2m)
135
136!===================================================================c
137! Valeurs numeriques des constantes
138!===================================================================c
139
140
141! Minimum du carre du vent
142
143 CEPDU2 = (0.1)**2
144
145! Louis 1982
146
147  CB=5.0
148  CC=5.0
149  CD=5.0
150
151
152! King 2001
153
154  C2=0.25
155  C3=0.0625
156
157 
158! Louis 1979
159
160  BPRIME=4.7
161  B=9.4
162 
163
164!MO
165
166  ALPHA=5.0
167 
168
169! ================================================================= c
170! Tests avant de commencer
171! Fuxing WANG, 04/03/2015
172! To check if there are negative q1, qsurf values.
173!====================================================================c
174  ng_q1 = 0     ! Initialization
175  ng_qsurf = 0  ! Initialization
176  DO i = 1, knon
177     IF (q1(i).LT.0.0)     ng_q1 = ng_q1 + 1
178     IF (qsurf(i).LT.0.0)  ng_qsurf = ng_qsurf + 1
179  ENDDO
180  IF (ng_q1.GT.0 .and. prt_level > 5) THEN
181      WRITE(lunout,*)" *** Warning: Negative q1(humidity at 1st level) values in cdrag.F90 !"
182      WRITE(lunout,*)" The total number of the grids is: ", ng_q1
183      WRITE(lunout,*)" The negative q1 is set to zero "
184!      abort_message="voir ci-dessus"
185!      CALL abort_physic(modname,abort_message,1)
186  ENDIF
187  IF (ng_qsurf.GT.0 .and. prt_level > 5) THEN
188      WRITE(lunout,*)" *** Warning: Negative qsurf(humidity at surface) values in cdrag.F90 !"
189      WRITE(lunout,*)" The total number of the grids is: ", ng_qsurf
190      WRITE(lunout,*)" The negative qsurf is set to zero "
191!      abort_message="voir ci-dessus"
192!      CALL abort_physic(modname,abort_message,1)
193  ENDIF
194
195
196
197!=============================================================================c
198! Calcul du cdrag
199!=============================================================================c
200
201! On choisit les fonctions de stabilite utilisees au premier appel
202!**************************************************************************
203  IF (firstcall) THEN
204   iflag_corr_sta=2
205   iflag_corr_insta=2
206   iflag_new_t2mq2m=0
207
208   CALL getin_p('iflag_corr_sta',iflag_corr_sta)
209   CALL getin_p('iflag_corr_insta',iflag_corr_insta)
210   CALL getin_p('iflag_new_t2mq2m',iflag_new_t2mq2m)
211   firstcall = .FALSE.
212 ENDIF
213
214!xxxxxxxxxxxxxxxxxxxxxxx
215  DO i = 1, knon        ! Boucle sur l'horizontal
216!xxxxxxxxxxxxxxxxxxxxxxx
217
218
219! calculs preliminaires:
220!***********************
221     
222
223     zdu2 = MAX(CEPDU2, speed(i)**2)
224     pref(i) = EXP(LOG(psol(i)) - zgeop1(i)/(RD*t1(i)* &
225                 (1.+ RETV * max(q1(i),0.0))))           ! negative q1 set to zero
226     ztsolv = tsurf(i) * (1.0+RETV*max(qsurf(i),0.0))    ! negative qsurf set to zero
227     IF (iflag_new_t2mq2m==0) THEN
228        ztvd = (t1(i)+zgeop1(i)/RCPD/(1.+RVTMP2*q1(i))) &
229             *(1.+RETV*max(q1(i),0.0))                   ! negative q1 set to zero
230     ELSE
231        ztvd = (t1(i)+zgeop1(i)/RCPD/(1.+RVTMP2*max(q1(i),0.0))) &
232             *(1.+RETV*max(q1(i),0.0))                   ! negative q1 set to zero
233     ENDIF
234     
235     zri(i) = zgeop1(i)*(ztvd-ztsolv)/(zdu2*ztvd)
236
237
238! Coefficients CD neutres : k^2/ln(z/z0) et k^2/(ln(z/z0)*ln(z/z0h)):
239!********************************************************************
240
241     zzzcd=CKAP/LOG(1.+zgeop1(i)/(RG*z0m(i)))
242     cdmn(i) = zzzcd*zzzcd
243     cdhn(i) = zzzcd*(CKAP/LOG(1.+zgeop1(i)/(RG*z0h(i))))
244
245
246! Calcul des fonctions de stabilit?? FMs, FHs, FMi, FHi :
247!*******************************************************
248
249!''''''''''''''
250! Cas instables
251!''''''''''''''
252
253 IF (zri(i) .LT. 0.) THEN   
254
255
256        SELECT CASE (iflag_corr_insta)
257   
258        CASE (1) ! Louis 1979 + Mascart 1995
259
260           MU=LOG(MAX(z0m(i)/z0h(i),0.01))
261           CMstar=6.8741+2.6933*MU-0.3601*(MU**2)+0.0154*(MU**3)
262           PM=0.5233-0.0815*MU+0.0135*(MU**2)-0.001*(MU**3)
263           CHstar=3.2165+4.3431*MU+0.536*(MU**2)-0.0781*(MU**3)
264           PH=0.5802-0.1571*MU+0.0327*(MU**2)-0.0026*(MU**3)
265           CH=CHstar*B*CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) &
266            & * CKAPT/LOG(z0h(i)+zgeop1(i)/(RG*z0h(i)))       &
267            & * ((zgeop1(i)/(RG*z0h(i)))**PH)
268           CM=CMstar*B*CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) &
269            & *CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i)))         &
270            & * ((zgeop1(i)/(RG*z0m(i)))**PM)
271         
272
273
274
275           FM(i)=1.-B*zri(i)/(1.+CM*SQRT(ABS(zri(i))))
276           FH(i)=1.-B*zri(i)/(1.+CH*SQRT(ABS(zri(i))))
277
278        CASE (2) ! Louis 1982
279
280           zucf = 1./(1.+3.0*CB*CC*cdmn(i)*SQRT(ABS(zri(i)) &
281                *(1.0+zgeop1(i)/(RG*z0m(i)))))
282           FM(i) = AMAX1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min)
283           FH(i) = AMAX1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min)
284
285
286        CASE (3) ! Laurent Li
287
288           
289           FM(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min)
290           FH(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min)
291
292
293
294         CASE default ! Louis 1982
295
296           zucf = 1./(1.+3.0*CB*CC*cdmn(i)*SQRT(ABS(zri(i)) &
297                *(1.0+zgeop1(i)/(RG*z0m(i)))))
298           FM(i) = AMAX1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min)
299           FH(i) = AMAX1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min)
300
301
302         END SELECT
303
304
305
306! Calcul des drags
307
308
309       cdm(i)=cdmn(i)*FM(i)
310       cdh(i)=f_cdrag_ter*cdhn(i)*FH(i)
311
312
313! Traitement particulier des cas oceaniques
314! on applique Miller et al 1992 en l'absence de gustiness
315
316  IF (nsrf == is_oce) THEN
317!        cdh(i)=f_cdrag_oce*cdhn(i)*FH(i)
318
319        IF(iflag_gusts==0) THEN
320           zcr = (0.0016/(cdmn(i)*SQRT(zdu2)))*ABS(ztvd-ztsolv)**(1./3.)
321           cdh(i) =f_cdrag_oce* cdhn(i)*(1.0+zcr**1.25)**(1./1.25)
322        ENDIF
323
324
325        cdm(i)=MIN(cdm(i),cdmmax)
326        cdh(i)=MIN(cdh(i),cdhmax)
327
328  END IF
329
330
331
332 ELSE                         
333
334!'''''''''''''''
335! Cas stables :
336!'''''''''''''''
337        zri(i) = MIN(20.,zri(i))
338
339       SELECT CASE (iflag_corr_sta)
340   
341        CASE (1) ! Louis 1979 + Mascart 1995
342   
343           FM(i)=MAX(1./((1+BPRIME*zri(i))**2),f_ri_cd_min)
344           FH(i)=FM(i)
345
346
347        CASE (2) ! Louis 1982
348           
349           zscf = SQRT(1.+CD*ABS(zri(i)))
350           FM(i)= AMAX1(1. / (1.+2.*CB*zri(i)/zscf), f_ri_cd_min)
351           FH(i)= AMAX1(1./ (1.+3.*CB*zri(i)*zscf), f_ri_cd_min )
352
353
354        CASE (3) ! Laurent Li
355   
356           FM(i)=MAX(1.0 / (1.0+10.0*zri(i)*(1+8.0*zri(i))),f_ri_cd_min)
357           FH(i)=FM(i)
358
359
360        CASE (4)  ! King 2001
361         
362           if (zri(i) .LT. C2/2.) then
363           FM(i)=MAX((1.-zri(i)/C2)**2,f_ri_cd_min)
364           FH(i)=  FM(i)
365
366
367           else
368           FM(i)=MAX(C3*((C2/zri(i))**2),f_ri_cd_min)
369           FH(i)= FM(i)
370           endif
371
372
373        CASE (5) ! MO
374   
375          if (zri(i) .LT. 1./alpha) then
376
377           FM(i)=MAX((1.-alpha*zri(i))**2,f_ri_cd_min)
378           FH(i)=FM(i)
379
380           else
381
382
383           FM(i)=MAX(1E-7,f_ri_cd_min)
384           FH(i)=FM(i)
385
386          endif
387
388
389
390
391
392        CASE default ! Louis 1982
393
394           zscf = SQRT(1.+CD*ABS(zri(i)))
395           FM(i)= AMAX1(1. / (1.+2.*CB*zri(i)/zscf), f_ri_cd_min)
396           FH(i)= AMAX1(1./ (1.+3.*CB*zri(i)*zscf), f_ri_cd_min )
397
398
399
400   END SELECT
401
402! Calcul des drags
403
404
405       cdm(i)=cdmn(i)*FM(i)
406       cdh(i)=f_cdrag_ter*cdhn(i)*FH(i)
407
408       IF(nsrf.EQ.is_oce) THEN
409
410        cdh(i)=f_cdrag_oce*cdhn(i)*FH(i)
411        cdm(i)=MIN(cdm(i),cdmmax)
412        cdh(i)=MIN(cdh(i),cdhmax)
413
414      ENDIF
415
416
417 ENDIF
418
419!xxxxxxxxxxx
420  END DO   !  Fin de la boucle sur l'horizontal
421!xxxxxxxxxxx
422! ================================================================= c
423     
424END SUBROUTINE cdrag
425
426!
427 SUBROUTINE cdragn_ri1(knon,  nsrf,   &
428     speed, t1,    q1,    zgeop1, &
429     psol,  tsurf, qsurf, z0m, z0h,  &
430     ri1, iri1, &
431     cdm,  cdh,  zri,   pref)
432
433  USE dimphy
434  USE indice_sol_mod
435  USE print_control_mod, ONLY: lunout, prt_level
436  USE ioipsl_getin_p_mod, ONLY : getin_p
437
438  IMPLICIT NONE
439! ================================================================= c
440!
441! Objet : calcul des cdrags pour le moment (pcfm) et
442!         les flux de chaleur sensible et latente (pcfh) d'apr??s
443!         Louis 1982, Louis 1979, King et al 2001
444!         ou Zilitinkevich et al 2002  pour les cas stables, Louis 1979
445!         et 1982 pour les cas instables
446!
447! Modified history:
448!  writting on the 20/05/2016
449!  modified on the 13/12/2016 to be adapted to LMDZ6
450!
451! References:
452!   Louis, J. F., 1979: A parametric model of vertical eddy fluxes in the
453!     atmosphere. Boundary-Layer Meteorology. 01/1979; 17(2):187-202.
454!   Louis, J. F., Tiedtke, M. and Geleyn, J. F., 1982: `A short history of the
455!     operational PBL parametrization at ECMWF'. Workshop on boundary layer
456!     parametrization, November 1981, ECMWF, Reading, England.
457!     Page: 19. Equations in Table 1.
458!   Mascart P, Noilhan J, Giordani H 1995.A MODIFIED PARAMETERIZATION OF FLUX-PROFILE RELATIONSHIPS
459!    IN THE SURFACE LAYER USING DIFFERENT ROUGHNESS LENGTH VALUES FOR HEAT AND MOMENTUM
460!    Boundary-Layer Meteorology 72: 331-344
461!   Anton Beljaars. May 1992. The parametrization of the planetary boundary layer. 
462!     European Centre for Medium-Range Weather Forecasts.
463!     Equations: 110-113. Page 40.
464!   Miller,M.J., A.C.M.Beljaars, T.N.Palmer. 1992. The sensitivity of the ECMWF
465!     model to the parameterization of evaporation from the tropical oceans. J.
466!     Climate, 5:418-434.
467!   King J.C, Connolley, W.M ad Derbyshire S.H. 2001, Sensitivity of Modelled Antarctic climate
468!   to surface and boundary-layer flux parametrizations
469!   QJRMS, 127, pp 779-794
470!
471! ================================================================= c
472! ================================================================= c
473! On choisit le couple de fonctions de correction avec deux flags:
474! Un pour les cas instables, un autre pour les cas stables
475!
476! iflag_corr_insta:
477!                   1: Louis 1979 avec les modifications de Mascart 1995 (z0/= z0h)
478!                   2: Louis 1982
479!                   3: Laurent Li
480!
481! iflag_corr_sta:
482!                   1: Louis 1979 avec les modifications de Mascart 1995 (z0/= z0h)
483!                   2: Louis 1982
484!                   3: Laurent Li
485!                   4: King  2001 (SHARP)
486!                   5: MO 1st order theory (allow collapse of turbulence)
487!           
488!
489!*****************************************************************
490! Parametres d'entree
491!*****************************************************************
492 
493  INTEGER, INTENT(IN)                      :: knon, nsrf ! nombre de points de grille sur l'horizontal + type de surface
494  REAL, DIMENSION(klon), INTENT(IN)        :: speed ! module du vent au 1er niveau du modele
495  REAL, DIMENSION(klon), INTENT(IN)        :: zgeop1! geopotentiel au 1er niveau du modele
496  REAL, DIMENSION(klon), INTENT(IN)        :: tsurf ! Surface temperature (K)
497  REAL, DIMENSION(klon), INTENT(IN)        :: qsurf ! Surface humidity (Kg/Kg)
498  REAL, DIMENSION(klon), INTENT(IN)        :: z0m, z0h ! Rugosity at surface (m)
499  REAL, DIMENSION(klon), INTENT(IN)        :: ri1   ! Richardson 1ere couche
500  INTEGER, INTENT(IN)                      :: iri1   !
501  REAL, DIMENSION(klon), INTENT(IN)        :: t1  ! Temperature au premier niveau (K)
502  REAL, DIMENSION(klon), INTENT(IN)        :: q1  ! humidite specifique au premier niveau (kg/kg)
503  REAL, DIMENSION(klon), INTENT(IN)        :: psol ! pression au sol
504
505
506
507! Parametres de sortie
508!******************************************************************
509  REAL, DIMENSION(klon), INTENT(OUT)       :: cdm  ! Drag coefficient for heat flux
510  REAL, DIMENSION(klon), INTENT(OUT)       :: cdh  ! Drag coefficient for momentum
511  REAL, DIMENSION(klon), INTENT(OUT)       :: zri   ! Richardson number
512  REAL, DIMENSION(klon), INTENT(OUT)       :: pref  ! Pression au niveau zgeop/RG
513! Variables Locales
514!******************************************************************
515 
516
517  INCLUDE "YOMCST.h"
518  INCLUDE "YOETHF.h"
519  INCLUDE "clesphys.h"
520
521
522  REAL, PARAMETER       :: CKAP=0.40, CKAPT=0.42
523  REAL                   CEPDU2
524  REAL                   ALPHA
525  REAL                   CB,CC,CD,C2,C3
526  REAL                   MU, CM, CH, B, CMstar, CHstar
527  REAL                   PM, PH, BPRIME
528  REAL                   C
529  INTEGER                ng_q1                      ! Number of grids that q1 < 0.0
530  INTEGER                ng_qsurf                   ! Number of grids that qsurf < 0.0
531  INTEGER                i
532  REAL                   zdu2, ztsolv
533  REAL                   ztvd, zscf
534  REAL                   zucf, zcr
535  REAL                   friv, frih
536  REAL, DIMENSION(klon) :: FM, FH                   ! stability functions
537  REAL, DIMENSION(klon) :: cdmn, cdhn               ! Drag coefficient in neutral conditions
538  REAL zzzcd
539
540  LOGICAL, SAVE :: firstcall = .TRUE.
541  !$OMP THREADPRIVATE(firstcall)
542  INTEGER, SAVE :: iflag_corr_sta
543  !$OMP THREADPRIVATE(iflag_corr_sta)
544  INTEGER, SAVE :: iflag_corr_insta
545  !$OMP THREADPRIVATE(iflag_corr_insta)
546
547!===================================================================c
548! Valeurs numeriques des constantes
549!===================================================================c
550
551
552! Minimum du carre du vent
553
554 CEPDU2 = (0.1)**2
555! Louis 1982
556
557  CB=5.0
558  CC=5.0
559  CD=5.0
560
561
562! King 2001
563
564  C2=0.25
565  C3=0.0625
566
567 
568! Louis 1979
569
570  BPRIME=4.7
571  B=9.4
572 
573
574!MO
575
576  ALPHA=5.0
577 
578
579! ================================================================= c
580! Tests avant de commencer
581! Fuxing WANG, 04/03/2015
582! To check if there are negative q1, qsurf values.
583!====================================================================c
584  ng_q1 = 0     ! Initialization
585  ng_qsurf = 0  ! Initialization
586  DO i = 1, knon
587     IF (q1(i).LT.0.0)     ng_q1 = ng_q1 + 1
588     IF (qsurf(i).LT.0.0)  ng_qsurf = ng_qsurf + 1
589  ENDDO
590  IF (ng_q1.GT.0 .and. prt_level > 5) THEN
591      WRITE(lunout,*)" *** Warning: Negative q1(humidity at 1st level) values in cdrag.F90 !"
592      WRITE(lunout,*)" The total number of the grids is: ", ng_q1
593      WRITE(lunout,*)" The negative q1 is set to zero "
594!      abort_message="voir ci-dessus"
595!      CALL abort_physic(modname,abort_message,1)
596  ENDIF
597  IF (ng_qsurf.GT.0 .and. prt_level > 5) THEN
598      WRITE(lunout,*)" *** Warning: Negative qsurf(humidity at surface) values in cdrag.F90 !"
599      WRITE(lunout,*)" The total number of the grids is: ", ng_qsurf
600      WRITE(lunout,*)" The negative qsurf is set to zero "
601!      abort_message="voir ci-dessus"
602!      CALL abort_physic(modname,abort_message,1)
603  ENDIF
604
605
606
607!=============================================================================c
608! Calcul du cdrag
609!=============================================================================c
610
611! On choisit les fonctions de stabilite utilisees au premier appel
612!**************************************************************************
613  IF (firstcall) THEN
614   iflag_corr_sta=2
615   iflag_corr_insta=2
616 
617   CALL getin_p('iflag_corr_sta',iflag_corr_sta)
618   CALL getin_p('iflag_corr_insta',iflag_corr_insta)
619
620   firstcall = .FALSE.
621 ENDIF
622
623!xxxxxxxxxxxxxxxxxxxxxxx
624  DO i = 1, knon        ! Boucle sur l'horizontal
625!xxxxxxxxxxxxxxxxxxxxxxx
626
627
628! calculs preliminaires:
629!***********************
630     
631
632     zdu2 = MAX(CEPDU2, speed(i)**2)
633     pref(i) = EXP(LOG(psol(i)) - zgeop1(i)/(RD*t1(i)* &
634                 (1.+ RETV * max(q1(i),0.0))))           ! negative q1 set to zero
635     ztsolv = tsurf(i) * (1.0+RETV*max(qsurf(i),0.0))    ! negative qsurf set to zero
636     ztvd = (t1(i)+zgeop1(i)/RCPD/(1.+RVTMP2*max(q1(i),0.0))) &
637          *(1.+RETV*max(q1(i),0.0))                      ! negative q1 set to zero
638     zri(i) = zgeop1(i)*(ztvd-ztsolv)/(zdu2*ztvd)
639
640     IF (iri1.EQ.1) THEN
641       zri(i) = ri1(i)
642     ENDIF
643
644! Coefficients CD neutres : k^2/ln(z/z0) et k^2/(ln(z/z0)*ln(z/z0h)):
645!********************************************************************
646
647     zzzcd=CKAP/LOG(1.+zgeop1(i)/(RG*z0m(i)))
648     cdmn(i) = zzzcd*zzzcd
649     cdhn(i) = zzzcd*(CKAP/LOG(1.+zgeop1(i)/(RG*z0h(i))))
650
651
652! Calcul des fonctions de stabilit?? FMs, FHs, FMi, FHi :
653!*******************************************************
654
655!''''''''''''''
656! Cas instables
657!''''''''''''''
658
659 IF (zri(i) .LT. 0.) THEN   
660
661
662        SELECT CASE (iflag_corr_insta)
663   
664        CASE (1) ! Louis 1979 + Mascart 1995
665
666           MU=LOG(MAX(z0m(i)/z0h(i),0.01))
667           CMstar=6.8741+2.6933*MU-0.3601*(MU**2)+0.0154*(MU**3)
668           PM=0.5233-0.0815*MU+0.0135*(MU**2)-0.001*(MU**3)
669           CHstar=3.2165+4.3431*MU+0.536*(MU**2)-0.0781*(MU**3)
670           PH=0.5802-0.1571*MU+0.0327*(MU**2)-0.0026*(MU**3)
671           CH=CHstar*B*CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) &
672            & * CKAPT/LOG(z0h(i)+zgeop1(i)/(RG*z0h(i)))       &
673            & * ((zgeop1(i)/(RG*z0h(i)))**PH)
674           CM=CMstar*B*CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) &
675            & *CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i)))         &
676            & * ((zgeop1(i)/(RG*z0m(i)))**PM)
677         
678
679
680
681           FM(i)=1.-B*zri(i)/(1.+CM*SQRT(ABS(zri(i))))
682           FH(i)=1.-B*zri(i)/(1.+CH*SQRT(ABS(zri(i))))
683
684        CASE (2) ! Louis 1982
685
686           zucf = 1./(1.+3.0*CB*CC*cdmn(i)*SQRT(ABS(zri(i)) &
687                *(1.0+zgeop1(i)/(RG*z0m(i)))))
688           FM(i) = AMAX1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min)
689           FH(i) = AMAX1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min)
690
691
692        CASE (3) ! Laurent Li
693
694           
695           FM(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min)
696           FH(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min)
697
698
699
700         CASE default ! Louis 1982
701
702           zucf = 1./(1.+3.0*CB*CC*cdmn(i)*SQRT(ABS(zri(i)) &
703                *(1.0+zgeop1(i)/(RG*z0m(i)))))
704           FM(i) = AMAX1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min)
705           FH(i) = AMAX1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min)
706
707
708         END SELECT
709
710
711
712! Calcul des drags
713
714
715       cdm(i)=cdmn(i)*FM(i)
716       cdh(i)=f_cdrag_ter*cdhn(i)*FH(i)
717
718
719! Traitement particulier des cas oceaniques
720! on applique Miller et al 1992 en l'absence de gustiness
721
722  IF (nsrf == is_oce) THEN
723!        cdh(i)=f_cdrag_oce*cdhn(i)*FH(i)
724
725        IF(iflag_gusts==0) THEN
726           zcr = (0.0016/(cdmn(i)*SQRT(zdu2)))*ABS(ztvd-ztsolv)**(1./3.)
727           cdh(i) =f_cdrag_oce* cdhn(i)*(1.0+zcr**1.25)**(1./1.25)
728        ENDIF
729
730
731        cdm(i)=MIN(cdm(i),cdmmax)
732        cdh(i)=MIN(cdh(i),cdhmax)
733
734  END IF
735
736
737
738 ELSE                         
739
740!'''''''''''''''
741! Cas stables :
742!'''''''''''''''
743        zri(i) = MIN(20.,zri(i))
744
745       SELECT CASE (iflag_corr_sta)
746   
747        CASE (1) ! Louis 1979 + Mascart 1995
748   
749           FM(i)=MAX(1./((1+BPRIME*zri(i))**2),f_ri_cd_min)
750           FH(i)=FM(i)
751
752
753        CASE (2) ! Louis 1982
754           
755           zscf = SQRT(1.+CD*ABS(zri(i)))
756           FM(i)= AMAX1(1. / (1.+2.*CB*zri(i)/zscf), f_ri_cd_min)
757           FH(i)= AMAX1(1./ (1.+3.*CB*zri(i)*zscf), f_ri_cd_min )
758
759
760        CASE (3) ! Laurent Li
761   
762           FM(i)=MAX(1.0 / (1.0+10.0*zri(i)*(1+8.0*zri(i))),f_ri_cd_min)
763           FH(i)=FM(i)
764
765
766        CASE (4)  ! King 2001
767         
768           if (zri(i) .LT. C2/2.) then
769           FM(i)=MAX((1.-zri(i)/C2)**2,f_ri_cd_min)
770           FH(i)=  FM(i)
771
772
773           else
774           FM(i)=MAX(C3*((C2/zri(i))**2),f_ri_cd_min)
775           FH(i)= FM(i)
776           endif
777
778
779        CASE (5) ! MO
780   
781          if (zri(i) .LT. 1./alpha) then
782
783           FM(i)=MAX((1.-alpha*zri(i))**2,f_ri_cd_min)
784           FH(i)=FM(i)
785
786           else
787
788
789           FM(i)=MAX(1E-7,f_ri_cd_min)
790           FH(i)=FM(i)
791
792          endif
793
794
795
796
797
798        CASE default ! Louis 1982
799
800           zscf = SQRT(1.+CD*ABS(zri(i)))
801           FM(i)= AMAX1(1. / (1.+2.*CB*zri(i)/zscf), f_ri_cd_min)
802           FH(i)= AMAX1(1./ (1.+3.*CB*zri(i)*zscf), f_ri_cd_min )
803
804
805
806   END SELECT
807
808! Calcul des drags
809
810
811       cdm(i)=cdmn(i)*FM(i)
812       cdh(i)=f_cdrag_ter*cdhn(i)*FH(i)
813
814       IF(nsrf.EQ.is_oce) THEN
815
816        cdh(i)=f_cdrag_oce*cdhn(i)*FH(i)
817        cdm(i)=MIN(cdm(i),cdmmax)
818        cdh(i)=MIN(cdh(i),cdhmax)
819
820      ENDIF
821
822
823 ENDIF
824
825!xxxxxxxxxxx
826  END DO   !  Fin de la boucle sur l'horizontal
827!xxxxxxxxxxx
828! ================================================================= c
829     
830END SUBROUTINE cdragn_ri1
831
832END MODULE cdrag_mod
Note: See TracBrowser for help on using the repository browser.