source: LMDZ6/branches/cirrus/libf/phylmd/cdrag_mod.F90 @ 5434

Last change on this file since 5434 was 5202, checked in by Laurent Fairhead, 3 months ago

Updating cirrus branch to trunk revision 5171

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