source: LMDZ6/trunk/libf/phylmd/cdrag_mod.F90 @ 4440

Last change on this file since 4440 was 3817, checked in by musat, 4 years ago

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