source: LMDZ6/branches/contrails/libf/phylmd/cdrag_mod.f90 @ 5431

Last change on this file since 5431 was 5285, checked in by abarral, 2 months ago

As discussed internally, remove generic ONLY: ... for new _mod_h modules

  • Property svn:executable set to *
  • Property svn:keywords set to Id
File size: 25.7 KB
Line 
1!
2!$Id: cdrag_mod.f90 5285 2024-10-28 13:33:29Z aborella $
3!
4!
5MODULE cdrag_mod
6!
7! This module contains some procedures for calculation of the cdrag
8! coefficients for turbulent diffusion at surface
9!
10  IMPLICIT NONE
11
12CONTAINS
13!
14!****************************************************************************************
15!
16!r original routine svn3623
17!
18 SUBROUTINE cdrag(knon,  nsrf,   &
19     speed, t1,    q1,    zgeop1, &
20     psol, pblh,  tsurf, qsurf, z0m, z0h,  &
21     ri_in, iri_in, &
22     cdm,  cdh,  zri,   pref, prain, tsol , pat1)
23
24  USE dimphy
25  USE coare_cp_mod, ONLY: coare_cp
26  USE coare30_flux_cnrm_mod, ONLY: coare30_flux_cnrm
27  USE indice_sol_mod
28  USE print_control_mod, ONLY: lunout, prt_level
29  USE ioipsl_getin_p_mod, ONLY : getin_p
30  USE lmdz_atke_turbulence_ini, ONLY : smmin, ric, cinf, cepsilon, pr_slope, pr_asym, pr_neut
31  USE yomcst_mod_h
32  USE clesphys_mod_h
33  USE yoethf_mod_h
34  IMPLICIT NONE
35! ================================================================= c
36!
37! Objet : calcul des cdrags pour le moment (pcfm) et
38!         les flux de chaleur sensible et latente (pcfh) d'apr??s
39!         Louis 1982, Louis 1979, King et al 2001
40!         ou Zilitinkevich et al 2002  pour les cas stables, Louis 1979
41!         et 1982 pour les cas instables
42!
43! Modified history:
44!  writting on the 20/05/2016
45!  modified on the 13/12/2016 to be adapted to LMDZ6
46!
47! References:
48!   Louis, J. F., 1979: A parametric model of vertical eddy fluxes in the
49!     atmosphere. Boundary-Layer Meteorology. 01/1979; 17(2):187-202.
50!   Louis, J. F., Tiedtke, M. and Geleyn, J. F., 1982: `A short history of the
51!     operational PBL parametrization at ECMWF'. Workshop on boundary layer
52!     parametrization, November 1981, ECMWF, Reading, England.
53!     Page: 19. Equations in Table 1.
54!   Mascart P, Noilhan J, Giordani H 1995.A MODIFIED PARAMETERIZATION OF FLUX-PROFILE RELATIONSHIPS
55!    IN THE SURFACE LAYER USING DIFFERENT ROUGHNESS LENGTH VALUES FOR HEAT AND MOMENTUM
56!    Boundary-Layer Meteorology 72: 331-344
57!   Anton Beljaars. May 1992. The parametrization of the planetary boundary layer. 
58!     European Centre for Medium-Range Weather Forecasts.
59!     Equations: 110-113. Page 40.
60!   Miller,M.J., A.C.M.Beljaars, T.N.Palmer. 1992. The sensitivity of the ECMWF
61!     model to the parameterization of evaporation from the tropical oceans. J.
62!     Climate, 5:418-434.
63!   King J.C, Connolley, W.M ad Derbyshire S.H. 2001, Sensitivity of Modelled Antarctic climate
64!   to surface and boundary-layer flux parametrizations
65!   QJRMS, 127, pp 779-794
66!
67! ================================================================= c
68! ================================================================= c
69! On choisit le couple de fonctions de correction avec deux flags:
70! Un pour les cas instables, un autre pour les cas stables
71!
72! iflag_corr_insta:
73!                   1: Louis 1979 avec les modifications de Mascart 1995 (z0/= z0h)
74!                   2: Louis 1982
75!                   3: Laurent Li
76!
77! iflag_corr_sta:
78!                   1: Louis 1979 avec les modifications de Mascart 1995 (z0/= z0h)
79!                   2: Louis 1982
80!                   3: Laurent Li
81!                   4: King  2001 (SHARP)
82!                   5: MO 1st order theory (allow collapse of turbulence)
83!           
84!
85!*****************************************************************
86! Parametres d'entree
87!*****************************************************************
88 
89  INTEGER, INTENT(IN)                      :: knon, nsrf ! nombre de points de grille sur l'horizontal + type de surface
90  REAL, DIMENSION(klon), INTENT(IN)        :: speed    ! module du vent au 1er niveau du modele
91  REAL, DIMENSION(klon), INTENT(IN)        :: zgeop1   ! geopotentiel au 1er niveau du modele
92  REAL, DIMENSION(klon), INTENT(IN)        :: tsurf    ! Surface temperature (K)
93  REAL, DIMENSION(klon), INTENT(IN)        :: qsurf    ! Surface humidity (Kg/Kg)
94  REAL, DIMENSION(klon), INTENT(INOUT)     :: z0m, z0h ! Rugosity at surface (m)
95  REAL, DIMENSION(klon), INTENT(IN)        :: ri_in ! Input Richardson 1st layer for first guess calculations of screen var.
96  INTEGER, INTENT(IN)                      :: iri_in! iflag to activate cdrag calculation using ri1
97  REAL, DIMENSION(klon), INTENT(IN)        :: t1       ! Temperature au premier niveau (K)
98  REAL, DIMENSION(klon), INTENT(IN)        :: q1       ! humidite specifique au premier niveau (kg/kg)
99  REAL, DIMENSION(klon), INTENT(IN)        :: psol     ! pression au sol
100
101!------------------ Rajout pour COARE (OT2018) --------------------
102  REAL, DIMENSION(klon), INTENT(INOUT)     :: pblh  !hauteur de CL
103  REAL, DIMENSION(klon), INTENT(IN)        :: prain !rapport au precipitation
104  REAL, DIMENSION(klon), INTENT(IN)        :: tsol  !SST imposé sur la surface oceanique
105  REAL, DIMENSION(klon), INTENT(IN)        :: pat1  !pression premier lev
106
107
108
109! Parametres de sortie
110!******************************************************************
111  REAL, DIMENSION(klon), INTENT(OUT)       :: cdm   ! Drag coefficient for momentum
112  REAL, DIMENSION(klon), INTENT(OUT)       :: cdh   ! Drag coefficient for heat flux
113  REAL, DIMENSION(klon), INTENT(OUT)       :: zri   ! Richardson number
114  REAL, DIMENSION(klon), INTENT(INOUT)       :: pref  ! Pression au niveau zgeop/RG
115
116! Variables Locales
117!******************************************************************
118
119  REAL, PARAMETER       :: CKAP=0.40, CKAPT=0.42
120  REAL                   CEPDU2
121  REAL                   ALPHA
122  REAL                   CB,CC,CD,C2,C3
123  REAL                   MU, CM, CH, B, CMstar, CHstar
124  REAL                   PM, PH, BPRIME
125  INTEGER                ng_q1                      ! Number of grids that q1 < 0.0
126  INTEGER                ng_qsurf                   ! Number of grids that qsurf < 0.0
127  INTEGER                i, k
128  REAL                   zdu2, ztsolv
129  REAL                   ztvd, zscf
130  REAL                   zucf, zcr
131  REAL, DIMENSION(klon) :: FM, FH                   ! stability functions
132  REAL, DIMENSION(klon) :: cdmn, cdhn               ! Drag coefficient in neutral conditions
133  REAL zzzcd
134  REAL, DIMENSION(klon) :: sm, prandtl              ! Stability function from atke turbulence scheme
135  REAL                  ri0, ri1, cn                ! to have dimensionless term in sm and prandtl
136
137!------------------ Rajout (OT2018) --------------------
138!------------------ Rajout pour les appelles BULK (OT) --------------------
139  REAL, DIMENSION(klon,2) :: rugos_itm
140  REAL, DIMENSION(klon,2) :: rugos_ith
141  REAL, PARAMETER         :: tol_it_rugos=1.e-4
142  REAL, DIMENSION(3)      :: coeffs
143  LOGICAL                 :: mixte
144  REAL                    :: z_0m
145  REAL                    :: z_0h
146  REAL, DIMENSION(klon)   :: cdmm
147  REAL, DIMENSION(klon)   :: cdhh
148
149!------------------RAJOUT POUR ECUME -------------------
150
151  REAL, DIMENSION(klon) :: PQSAT
152  REAL, DIMENSION(klon) :: PSFTH
153  REAL, DIMENSION(klon) :: PFSTQ
154  REAL, DIMENSION(klon) :: PUSTAR
155  REAL, DIMENSION(klon) :: PCD      ! Drag coefficient for momentum
156  REAL, DIMENSION(klon) :: PCDN     ! Drag coefficient for momentum
157  REAL, DIMENSION(klon) :: PCH      ! Drag coefficient for momentum
158  REAL, DIMENSION(klon) :: PCE      ! Drag coefficient for momentum
159  REAL, DIMENSION(klon) :: PRI
160  REAL, DIMENSION(klon) :: PRESA
161  REAL, DIMENSION(klon) :: PSSS
162
163  LOGICAL               :: OPRECIP
164  LOGICAL               :: OPWEBB
165  LOGICAL               :: OPERTFLUX
166  LOGICAL               :: LPRECIP
167  LOGICAL               :: LPWG
168
169
170
171  LOGICAL, SAVE :: firstcall = .TRUE.
172  !$OMP THREADPRIVATE(firstcall)
173  INTEGER, SAVE :: iflag_corr_sta
174  !$OMP THREADPRIVATE(iflag_corr_sta)
175  INTEGER, SAVE :: iflag_corr_insta
176  !$OMP THREADPRIVATE(iflag_corr_insta)
177  LOGICAL, SAVE :: ok_cdrag_iter
178  !$OMP THREADPRIVATE(ok_cdrag_iter)
179
180!===================================================================c
181! Valeurs numeriques des constantes
182!===================================================================c
183
184
185! Minimum du carre du vent
186
187 CEPDU2 = (0.1)**2
188
189! Louis 1982
190
191  CB=5.0
192  CC=5.0
193  CD=5.0
194
195
196! King 2001
197
198  C2=0.25
199  C3=0.0625
200
201 
202! Louis 1979
203
204  BPRIME=4.7
205  B=9.4
206 
207
208!MO
209
210  ALPHA=5.0
211
212! Consistent with atke scheme
213  cn=(1./sqrt(cepsilon))**(2./3.)
214  ri0=2./rpi*(cinf - cn)*ric/cn
215  ri1=-2./rpi * (pr_asym - pr_neut)/pr_slope
216
217
218! ================================================================= c
219! Tests avant de commencer
220! Fuxing WANG, 04/03/2015
221! To check if there are negative q1, qsurf values.
222!====================================================================c
223  ng_q1 = 0     ! Initialization
224  ng_qsurf = 0  ! Initialization
225  DO i = 1, knon
226     IF (q1(i).LT.0.0)     ng_q1 = ng_q1 + 1
227     IF (qsurf(i).LT.0.0)  ng_qsurf = ng_qsurf + 1
228  ENDDO
229  IF (ng_q1.GT.0 .and. prt_level > 5) THEN
230      WRITE(lunout,*)" *** Warning: Negative q1(humidity at 1st level) values in cdrag.F90 !"
231      WRITE(lunout,*)" The total number of the grids is: ", ng_q1
232      WRITE(lunout,*)" The negative q1 is set to zero "
233!      abort_message="voir ci-dessus"
234!      CALL abort_physic(modname,abort_message,1)
235  ENDIF
236  IF (ng_qsurf.GT.0 .and. prt_level > 5) THEN
237      WRITE(lunout,*)" *** Warning: Negative qsurf(humidity at surface) values in cdrag.F90 !"
238      WRITE(lunout,*)" The total number of the grids is: ", ng_qsurf
239      WRITE(lunout,*)" The negative qsurf is set to zero "
240!      abort_message="voir ci-dessus"
241!      CALL abort_physic(modname,abort_message,1)
242  ENDIF
243
244
245
246!=============================================================================c
247! Calcul du cdrag
248!=============================================================================c
249
250! On choisit les fonctions de stabilite utilisees au premier appel
251!**************************************************************************
252 IF (firstcall) THEN
253   iflag_corr_sta=2
254   iflag_corr_insta=2
255   ok_cdrag_iter = .FALSE.
256 
257   CALL getin_p('iflag_corr_sta',iflag_corr_sta)
258   CALL getin_p('iflag_corr_insta',iflag_corr_insta)
259   CALL getin_p('ok_cdrag_iter',ok_cdrag_iter)
260
261   firstcall = .FALSE.
262 ENDIF
263
264!------------------ Rajout (OT2018) --------------------
265!---------    Rajout pour itération sur rugosité     ----------------
266  rugos_itm(:,2) = z0m
267  rugos_itm(:,1) = 3.*tol_it_rugos*z0m
268
269  rugos_ith(:,2) = z0h !cp nouveau rugos_it
270  rugos_ith(:,1) = 3.*tol_it_rugos*z0h
271!--------------------------------------------------------------------
272
273!xxxxxxxxxxxxxxxxxxxxxxx
274  DO i = 1, knon        ! Boucle sur l'horizontal
275!xxxxxxxxxxxxxxxxxxxxxxx
276
277
278! calculs preliminaires:
279!***********************
280     
281     zdu2 = MAX(CEPDU2, speed(i)**2)
282     pref(i) = EXP(LOG(psol(i)) - zgeop1(i)/(RD*t1(i)* &
283          (1.+ RETV * max(q1(i),0.0))))           ! negative q1 set to zero
284     ztsolv = tsurf(i) * (1.0+RETV*max(qsurf(i),0.0))    ! negative qsurf set to zero
285     ztvd = (t1(i)+zgeop1(i)/RCPD/(1.+RVTMP2*q1(i))) &
286          *(1.+RETV*max(q1(i),0.0))                      ! negative q1 set to zero
287     
288!------------------ Rajout (OT2018) --------------------
289!--------------   ON DUPLIQUE POUR METTRE ITERATION SUR OCEAN    ------------------------     
290     IF (iri_in.EQ.1) THEN
291       zri(i) = ri_in(i)
292     ENDIF
293
294     IF (nsrf == is_oce) THEN
295       
296!------------------  Pour Core 2 choix Core Pur ou Core Mixte  --------------------
297       IF ((choix_bulk > 1 .AND. choix_bulk < 4) .AND. (nsrf .eq. is_oce)) THEN
298         IF(choix_bulk .eq. 2) THEN
299           mixte = .false.
300         ELSE
301            mixte = .true.
302         ENDIF
303         call clc_core_cp ( sqrt(zdu2),t1(i)-tsurf(i),q1(i)-qsurf(i),t1(i),q1(i),&
304             zgeop1(i)/RG, zgeop1(i)/RG,  zgeop1(i)/RG,&
305             psol(i),nit_bulk,mixte,&
306             coeffs,z_0m,z_0h)
307         cdmm(i) = coeffs(1)
308         cdhh(i) = coeffs(2)
309         cdm(i)=cdmm(i)
310         cdh(i)=cdhh(i)
311         write(*,*) "clc_core cd ch",cdmm(i),cdhh(i)
312
313!------------------                 Pour ECUME                 --------------------
314       ELSE IF ((choix_bulk .eq. 4) .and. (nsrf .eq. is_oce)) THEN
315         OPRECIP = .false.
316         OPWEBB  = .false.
317         OPERTFLUX = .false.
318         IF (nsrf .eq. is_oce) THEN
319           PSSS = 0.0
320         ENDIF
321         call ini_csts
322         call ecumev6_flux( z_0m,t1(i),tsurf(i),&
323             q1(i),qsurf(i),sqrt(zdu2),zgeop1(i)/RG,PSSS,zgeop1(i)/RG,&
324             psol(i),pat1(i), OPRECIP, OPWEBB,&
325             PSFTH,PFSTQ,PUSTAR,PCD,PCDN,&
326             PCH,PCE,PRI,PRESA,prain,&
327             z_0h,OPERTFLUX,coeffs)
328
329         cdmm(i) = coeffs(1)
330         cdhh(i) = coeffs(2)
331         cdm(i)=cdmm(i)
332         cdh(i)=cdhh(i)
333   
334         write(*,*) "ecume cd ch",cdmm(i),cdhh(i)
335
336!------------------                 Pour COARE CNRM                 --------------------
337       ELSE IF ((choix_bulk .eq. 5) .and. (nsrf .eq. is_oce)) THEN
338         LPRECIP = .false.
339         LPWG    = .false.
340         call ini_csts
341         block
342           real, dimension(1) :: z0m_1d, z_0h_1d, sqrt_zdu2_1d, zgeop1_rg_1d  ! convert scalar to 1D for call
343           z0m_1d = z0m
344           z_0h_1d = z0h
345           sqrt_zdu2_1d = sqrt(zdu2)
346           zgeop1_rg_1d=zgeop1(i)/RG
347           call coare30_flux_cnrm(z0m_1d,t1(i),tsurf(i), q1(i),  &
348               sqrt_zdu2_1d,zgeop1_rg_1d,zgeop1_rg_1d,psol(i),qsurf(i),PQSAT, &
349               PSFTH,PFSTQ,PUSTAR,PCD,PCDN,PCH,PCE,PRI, &
350               PRESA,prain,pat1(i),z_0h_1d, LPRECIP, LPWG, coeffs)
351
352         end block
353         cdmm(i) = coeffs(1)
354         cdhh(i) = coeffs(2)
355         cdm(i)=cdmm(i)
356         cdh(i)=cdhh(i)
357         write(*,*) "coare CNRM cd ch",cdmm(i),cdhh(i)
358
359!------------------                 Pour COARE Maison                 --------------------
360       ELSE IF ((choix_bulk .eq. 1) .and. (nsrf .eq. is_oce)) THEN
361         IF ( pblh(i) .eq. 0. ) THEN
362           pblh(i) = 1500.
363         ENDIF
364         write(*,*) "debug size",size(coeffs)
365         call coare_cp(sqrt(zdu2),t1(i)-tsurf(i),q1(i)-qsurf(i),&
366               t1(i),q1(i),&
367               zgeop1(i)/RG,zgeop1(i)/RG,zgeop1(i)/RG,&
368               psol(i), pblh(i),&
369               nit_bulk,&
370               coeffs,z_0m,z_0h)
371         cdmm(i) = coeffs(1)
372         cdhh(i) = coeffs(3)
373         cdm(i)=cdmm(i)
374         cdh(i)=cdhh(i)
375         write(*,*) "coare cd ch",cdmm(i),cdhh(i)
376       ELSE
377!------------------                 Pour La param LMDZ (ocean)              --------------------
378         zri(i) = zgeop1(i)*(ztvd-ztsolv)/(zdu2*ztvd)
379         IF (iri_in.EQ.1) THEN
380           zri(i) = ri_in(i)
381         ENDIF
382       ENDIF
383     
384
385!----------------------- Rajout des itérations --------------
386       DO  k=1,nit_bulk
387
388! Coefficients CD neutres : k^2/ln(z/z0) et k^2/(ln(z/z0)*ln(z/z0h)):
389!********************************************************************
390         zzzcd=CKAP/LOG(1.+zgeop1(i)/(RG*rugos_itm(i,2)))
391         cdmn(i) = zzzcd*zzzcd
392         cdhn(i) = zzzcd*(CKAP/LOG(1.+zgeop1(i)/(RG*rugos_ith(i,2))))
393
394! Calcul des fonctions de stabilite FMs, FHs, FMi, FHi :
395!*******************************************************
396         IF (zri(i) .LT. 0.) THEN   
397           SELECT CASE (iflag_corr_insta)
398             CASE (1) ! Louis 1979 + Mascart 1995
399                  MU=LOG(MAX(z0m(i)/z0h(i),0.01))
400                  CMstar=6.8741+2.6933*MU-0.3601*(MU**2)+0.0154*(MU**3)
401                  PM=0.5233-0.0815*MU+0.0135*(MU**2)-0.001*(MU**3)
402                  CHstar=3.2165+4.3431*MU+0.536*(MU**2)-0.0781*(MU**3)
403                  PH=0.5802-0.1571*MU+0.0327*(MU**2)-0.0026*(MU**3)
404                  CH=CHstar*B*CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) &
405                     & * CKAPT/LOG(z0h(i)+zgeop1(i)/(RG*z0h(i)))       &
406                     & * ((zgeop1(i)/(RG*z0h(i)))**PH)
407                  CM=CMstar*B*CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) &
408                     & *CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i)))         &
409                     & * ((zgeop1(i)/(RG*z0m(i)))**PM)
410                  FM(i)=1.-B*zri(i)/(1.+CM*SQRT(ABS(zri(i))))
411                  FH(i)=1.-B*zri(i)/(1.+CH*SQRT(ABS(zri(i))))
412             CASE (2) ! Louis 1982
413                  zucf = 1./(1.+3.0*CB*CC*cdmn(i)*SQRT(ABS(zri(i)) &
414                        *(1.0+zgeop1(i)/(RG*z0m(i)))))
415                  FM(i) = AMAX1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min)
416                  FH(i) = AMAX1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min)
417             CASE (3) ! Laurent Li
418                  FM(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min)
419                  FH(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min)
420             CASE (6) ! Consistent with turbulence scheme (in stationary case) derived in atke (2023)
421                  sm(i) = 2./rpi * (cinf - cn) * atan(-zri(i)/ri0) + cn
422                  prandtl(i) = -2./rpi * (pr_asym - pr_neut) * atan(zri(i)/ri1) + pr_neut
423                  FM(i) = MAX((sm(i)**(3./2.) * sqrt(cepsilon) * (1 - zri(i) / prandtl(i))**(1./2.)),f_ri_cd_min)
424                  FH(i) = MAX((FM(i) / prandtl(i)), f_ri_cd_min)
425             CASE default ! Louis 1982
426                  zucf = 1./(1.+3.0*CB*CC*cdmn(i)*SQRT(ABS(zri(i)) &
427                         *(1.0+zgeop1(i)/(RG*z0m(i)))))
428                  FM(i) = AMAX1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min)
429                  FH(i) = AMAX1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min)
430           END SELECT
431! Calcul des drags
432           cdmm(i)=cdmn(i)*FM(i)
433           cdhh(i)=f_cdrag_ter*cdhn(i)*FH(i)
434! Traitement particulier des cas oceaniques
435! on applique Miller et al 1992 en l'absence de gustiness
436           IF (nsrf == is_oce) THEN
437!            cdh(i)=f_cdrag_oce*cdhn(i)*FH(i)
438             IF (iflag_gusts==0) THEN
439               zcr = (0.0016/(cdmn(i)*SQRT(zdu2)))*ABS(ztvd-ztsolv)**(1./3.)
440               cdhh(i) =f_cdrag_oce* cdhn(i)*(1.0+zcr**1.25)**(1./1.25)
441             ENDIF
442             cdmm(i)=MIN(cdmm(i),cdmmax)
443             cdhh(i)=MIN(cdhh(i),cdhmax)
444!             write(*,*) "LMDZ cd ch",cdmm(i),cdhh(i)
445           END IF
446         ELSE                         
447
448!'''''''''''''''
449! Cas stables :
450!'''''''''''''''
451           zri(i) = MIN(20.,zri(i))
452           SELECT CASE (iflag_corr_sta)
453             CASE (1) ! Louis 1979 + Mascart 1995
454                  FM(i)=MAX(1./((1+BPRIME*zri(i))**2),f_ri_cd_min)
455                  FH(i)=FM(i)
456             CASE (2) ! Louis 1982
457                  zscf = SQRT(1.+CD*ABS(zri(i)))
458                  FM(i)= AMAX1(1. / (1.+2.*CB*zri(i)/zscf), f_ri_cd_min)
459                  FH(i)= AMAX1(1./ (1.+3.*CB*zri(i)*zscf), f_ri_cd_min )
460             CASE (3) ! Laurent Li
461                  FM(i)=MAX(1.0 / (1.0+10.0*zri(i)*(1+8.0*zri(i))),f_ri_cd_min)
462                  FH(i)=FM(i)
463             CASE (4)  ! King 2001
464                  IF (zri(i) .LT. C2/2.) THEN
465                    FM(i)=MAX((1.-zri(i)/C2)**2,f_ri_cd_min)
466                    FH(i)=  FM(i)
467                  ELSE
468                    FM(i)=MAX(C3*((C2/zri(i))**2),f_ri_cd_min)
469                    FH(i)= FM(i)
470                  ENDIF
471             CASE (5) ! MO
472                  if (zri(i) .LT. 1./alpha) then
473                    FM(i)=MAX((1.-alpha*zri(i))**2,f_ri_cd_min)
474                    FH(i)=FM(i)
475                  else
476                    FM(i)=MAX(1E-7,f_ri_cd_min)
477                    FH(i)=FM(i)
478                  endif
479             CASE (6) ! Consistent with turbulence scheme (in stationary case) derived in atke (2023)
480                  sm(i) = MAX(smmin,cn*(1.-zri(i)/ric))
481                           ! prandlt expression from venayagamoorthy and stretch 2010, Li et al 2019
482                  prandtl(i) = pr_neut*exp(-pr_slope/pr_neut*zri(i)+zri(i)/pr_neut) &
483                                + zri(i) * pr_slope
484                  FM(i) = MAX((sm(i)**(3./2.) * sqrt(cepsilon) * (1 - zri(i) / prandtl(i))**(1./2.)),f_ri_cd_min)
485                  FH(i) = MAX((FM(i) / prandtl(i)), f_ri_cd_min)
486             CASE default ! Louis 1982
487                  zscf = SQRT(1.+CD*ABS(zri(i)))
488                  FM(i)= AMAX1(1. / (1.+2.*CB*zri(i)/zscf), f_ri_cd_min)
489                  FH(i)= AMAX1(1./ (1.+3.*CB*zri(i)*zscf), f_ri_cd_min )
490           END SELECT
491
492! Calcul des drags
493           
494           cdmm(i)=cdmn(i)*FM(i)
495           cdhh(i)=f_cdrag_ter*cdhn(i)*FH(i)
496           
497           IF (choix_bulk == 0) THEN
498             cdm(i)=cdmn(i)*FM(i)
499             cdh(i)=f_cdrag_ter*cdhn(i)*FH(i)
500           ENDIF
501
502           IF (nsrf.EQ.is_oce) THEN
503             cdhh(i)=f_cdrag_oce*cdhn(i)*FH(i)
504             cdmm(i)=MIN(cdmm(i),cdmmax)
505             cdhh(i)=MIN(cdhh(i),cdhmax)
506           ENDIF
507           IF (ok_cdrag_iter) THEN
508             rugos_itm(i,1) = rugos_itm(i,2)
509             rugos_ith(i,1) = rugos_ith(i,2)
510             rugos_itm(i,2) =  0.018*cdmm(i) * (speed(i))/RG  &
511                              +  0.11*14e-6 / SQRT(cdmm(i) * zdu2)
512
513!---------- Version SEPARATION DES Z0 ----------------------
514             IF (iflag_z0_oce==0) THEN
515               rugos_ith(i,2) = rugos_itm(i,2)
516             ELSE IF (iflag_z0_oce==1) THEN
517               rugos_ith(i,2) = 0.40*14e-6 / SQRT(cdmm(i) * zdu2)
518             ENDIF
519           ENDIF
520         ENDIF
521         IF (ok_cdrag_iter) THEN
522           rugos_itm(i,2) = MAX(1.5e-05,rugos_itm(i,2))
523           rugos_ith(i,2) = MAX(1.5e-05,rugos_ith(i,2))
524         ENDIF
525       ENDDO
526       IF (nsrf.EQ.is_oce) THEN
527         cdm(i)=MIN(cdmm(i),cdmmax)
528         cdh(i)=MIN(cdhh(i),cdhmax)
529       ENDIF
530       z0m = rugos_itm(:,2)
531       z0h = rugos_ith(:,2)
532     ELSE ! (nsrf == is_oce)
533       zri(i) = zgeop1(i)*(ztvd-ztsolv)/(zdu2*ztvd)
534       IF (iri_in.EQ.1) THEN
535         zri(i) = ri_in(i)
536       ENDIF
537
538! Coefficients CD neutres : k^2/ln(z/z0) et k^2/(ln(z/z0)*ln(z/z0h)):
539!********************************************************************
540       zzzcd=CKAP/LOG(1.+zgeop1(i)/(RG*z0m(i)))
541       cdmn(i) = zzzcd*zzzcd
542       cdhn(i) = zzzcd*(CKAP/LOG(1.+zgeop1(i)/(RG*z0h(i))))
543
544
545! Calcul des fonctions de stabilit?? FMs, FHs, FMi, FHi :
546!*******************************************************
547!''''''''''''''
548! Cas instables
549!''''''''''''''
550       IF (zri(i) .LT. 0.) THEN   
551         SELECT CASE (iflag_corr_insta)
552           CASE (1) ! Louis 1979 + Mascart 1995
553                MU=LOG(MAX(z0m(i)/z0h(i),0.01))
554                CMstar=6.8741+2.6933*MU-0.3601*(MU**2)+0.0154*(MU**3)
555                PM=0.5233-0.0815*MU+0.0135*(MU**2)-0.001*(MU**3)
556                CHstar=3.2165+4.3431*MU+0.536*(MU**2)-0.0781*(MU**3)
557                PH=0.5802-0.1571*MU+0.0327*(MU**2)-0.0026*(MU**3)
558                CH=CHstar*B*CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) &
559                   & * CKAPT/LOG(z0h(i)+zgeop1(i)/(RG*z0h(i)))       &
560                   & * ((zgeop1(i)/(RG*z0h(i)))**PH)
561                CM=CMstar*B*CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) &
562                   & *CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i)))         &
563                   & * ((zgeop1(i)/(RG*z0m(i)))**PM)
564                FM(i)=1.-B*zri(i)/(1.+CM*SQRT(ABS(zri(i))))
565                FH(i)=1.-B*zri(i)/(1.+CH*SQRT(ABS(zri(i))))
566           CASE (2) ! Louis 1982
567                zucf = 1./(1.+3.0*CB*CC*cdmn(i)*SQRT(ABS(zri(i)) &
568                       *(1.0+zgeop1(i)/(RG*z0m(i)))))
569                FM(i) = AMAX1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min)
570                FH(i) = AMAX1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min)
571           CASE (3) ! Laurent Li
572                FM(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min)
573                FH(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min)
574           CASE (6) ! Consistent with turbulence scheme (in stationary case) derived in atke (2023)
575                 sm(i) = 2./rpi * (cinf - cn) * atan(-zri(i)/ri0) + cn
576                 prandtl(i) = -2./rpi * (pr_asym - pr_neut) * atan(zri(i)/ri1) + pr_neut
577                 FM(i) = MAX((sm(i)**(3./2.) * sqrt(cepsilon) * (1 - zri(i) / prandtl(i))**(1./2.)),f_ri_cd_min)
578                 FH(i) = MAX((FM(i) / prandtl(i)), f_ri_cd_min)
579           CASE default ! Louis 1982
580                zucf = 1./(1.+3.0*CB*CC*cdmn(i)*SQRT(ABS(zri(i)) &
581                      *(1.0+zgeop1(i)/(RG*z0m(i)))))
582                FM(i) = AMAX1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min)
583                FH(i) = AMAX1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min)
584         END SELECT
585! Calcul des drags
586         cdm(i)=cdmn(i)*FM(i)
587         cdh(i)=f_cdrag_ter*cdhn(i)*FH(i)
588       ELSE                         
589!'''''''''''''''
590! Cas stables :
591!'''''''''''''''
592         zri(i) = MIN(20.,zri(i))
593         SELECT CASE (iflag_corr_sta)
594           CASE (1) ! Louis 1979 + Mascart 1995
595                FM(i)=MAX(1./((1+BPRIME*zri(i))**2),f_ri_cd_min)
596                FH(i)=FM(i)
597           CASE (2) ! Louis 1982
598                zscf = SQRT(1.+CD*ABS(zri(i)))
599                FM(i)= AMAX1(1. / (1.+2.*CB*zri(i)/zscf), f_ri_cd_min)
600                FH(i)= AMAX1(1./ (1.+3.*CB*zri(i)*zscf), f_ri_cd_min )
601           CASE (3) ! Laurent Li
602                FM(i)=MAX(1.0 / (1.0+10.0*zri(i)*(1+8.0*zri(i))),f_ri_cd_min)
603                FH(i)=FM(i)
604           CASE (4)  ! King 2001
605                if (zri(i) .LT. C2/2.) then
606                  FM(i)=MAX((1.-zri(i)/C2)**2,f_ri_cd_min)
607                  FH(i)=  FM(i)
608                else
609                  FM(i)=MAX(C3*((C2/zri(i))**2),f_ri_cd_min)
610                  FH(i)= FM(i)
611                endif
612           CASE (5) ! MO
613                if (zri(i) .LT. 1./alpha) then
614                  FM(i)=MAX((1.-alpha*zri(i))**2,f_ri_cd_min)
615                  FH(i)=FM(i)
616                else
617                  FM(i)=MAX(1E-7,f_ri_cd_min)
618                  FH(i)=FM(i)
619                endif
620           CASE (6) ! Consistent with turbulence scheme (in stationary case) derived in atke (2023)
621                sm(i) = MAX(0.,cn*(1.-zri(i)/ric))
622                prandtl(i) = pr_neut + zri(i) * pr_slope
623                FM(i) = MAX((sm(i)**(3./2.) * sqrt(cepsilon) * (1 - zri(i) / prandtl(i))**(1./2.)),f_ri_cd_min)
624                FH(i) = MAX((FM(i) / prandtl(i)), f_ri_cd_min)
625           CASE default ! Louis 1982
626                zscf = SQRT(1.+CD*ABS(zri(i)))
627                FM(i)= AMAX1(1. / (1.+2.*CB*zri(i)/zscf), f_ri_cd_min)
628                FH(i)= AMAX1(1./ (1.+3.*CB*zri(i)*zscf), f_ri_cd_min )
629         END SELECT
630! Calcul des drags
631         cdm(i)=cdmn(i)*FM(i)
632         cdh(i)=f_cdrag_ter*cdhn(i)*FH(i)
633       ENDIF
634     ENDIF ! fin du if (nsrf == is_oce)
635  END DO   !  Fin de la boucle sur l'horizontal
636
637END SUBROUTINE cdrag
638
639END MODULE cdrag_mod
Note: See TracBrowser for help on using the repository browser.