source: LMDZ5/trunk/libf/phylmd/clcdrag.F90 @ 2201

Last change on this file since 2201 was 2126, checked in by fhourdin, 10 years ago

Introduction du cas Dice couplé atmosphere/surface
+ nouveau paramètre de contrôle f_ri_cd_min, seuil minimum
sur la fonction f(Ri) en facteur du coefficient de traîné neutre.
Par défaut : f_ri_cd_min=0.1 (comme avant)

Introduction of the coupled atmosphere/surface dice case.
+ new parameter f_ri_cd_min, minimum threshold on the f(Ri) fonction
that multiply the neutral drag coefficient.

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 4.9 KB
RevLine 
[524]1!
[1279]2!$Id: clcdrag.F90 2126 2014-10-15 00:03:57Z crio $
3!
[1067]4SUBROUTINE clcdrag(knon, nsrf, paprs, pplay,&
5     u1, v1, t1, q1, &
6     tsurf, qsurf, rugos, &
7     pcfm, pcfh)
8
9  USE dimphy
[1785]10  USE indice_sol_mod
11
[1067]12  IMPLICIT NONE
[524]13! ================================================================= c
14!
15! Objet : calcul des cdrags pour le moment (pcfm) et
16!         les flux de chaleur sensible et latente (pcfh).   
17!
18! ================================================================= c
19!
20! knon----input-I- nombre de points pour un type de surface
[1785]21! nsrf----input-I- indice pour le type de surface; voir indice_sol_mod.F90
[1067]22! u1-------input-R- vent zonal au 1er niveau du modele
23! v1-------input-R- vent meridien au 1er niveau du modele
24! t1-------input-R- temperature de l'air au 1er niveau du modele
25! q1-------input-R- humidite de l'air au 1er niveau du modele
26! tsurf------input-R- temperature de l'air a la surface
[524]27! qsurf---input-R- humidite de l'air a la surface
28! rugos---input-R- rugosite
29!
30! pcfm---output-R- cdrag pour le moment
31! pcfh---output-R- cdrag pour les flux de chaleur latente et sensible
32!
[1071]33  INTEGER, INTENT(IN)                      :: knon, nsrf
34  REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs
35  REAL, DIMENSION(klon,klev), INTENT(IN)   :: pplay
36  REAL, DIMENSION(klon), INTENT(IN)        :: u1, v1, t1, q1
37  REAL, DIMENSION(klon), INTENT(IN)        :: tsurf, qsurf
38  REAL, DIMENSION(klon), INTENT(IN)        :: rugos
39  REAL, DIMENSION(klon), INTENT(OUT)       :: pcfm, pcfh
40!
[524]41! ================================================================= c
42!
[1071]43  INCLUDE "YOMCST.h"
44  INCLUDE "YOETHF.h"
45  INCLUDE "clesphys.h"
[524]46!
47! Quelques constantes et options:
48!!$PB      REAL, PARAMETER :: ckap=0.35, cb=5.0, cc=5.0, cd=5.0, cepdu2=(0.1)**2
[1067]49  REAL, PARAMETER :: ckap=0.40, cb=5.0, cc=5.0, cd=5.0, cepdu2=(0.1)**2
[524]50!
51! Variables locales :
[1067]52  INTEGER               :: i
53  REAL                  :: zdu2, ztsolv
54  REAL                  :: ztvd, zscf
55  REAL                  :: zucf, zcr
56  REAL                  :: friv, frih
57  REAL, DIMENSION(klon) :: zcfm1, zcfm2
58  REAL, DIMENSION(klon) :: zcfh1, zcfh2
59  REAL, DIMENSION(klon) :: zcdn
60  REAL, DIMENSION(klon) :: zri
61  REAL, DIMENSION(klon) :: zgeop1       ! geopotentiel au 1er niveau du modele
62  LOGICAL, PARAMETER    :: zxli=.FALSE. ! calcul des cdrags selon Laurent Li
[524]63!
64! Fonctions thermodynamiques et fonctions d'instabilite
[1067]65  REAL                  :: fsta, fins, x
66  fsta(x) = 1.0 / (1.0+10.0*x*(1+8.0*x))
67  fins(x) = SQRT(1.0-18.0*x)
68
[524]69! ================================================================= c
70!
[1067]71! Calculer le geopotentiel du premier couche de modele
72!
73  DO i = 1, knon
74     zgeop1(i) = RD * t1(i) / (0.5*(paprs(i,1)+pplay(i,1))) &
75          * (paprs(i,1)-pplay(i,1))
76  END DO
77! ================================================================= c
78!
[524]79! Calculer le frottement au sol (Cdrag)
80!
[1067]81  DO i = 1, knon
82     zdu2 = MAX(cepdu2,u1(i)**2+v1(i)**2)
83     ztsolv = tsurf(i) * (1.0+RETV*qsurf(i))
84     ztvd = (t1(i)+zgeop1(i)/RCPD/(1.+RVTMP2*q1(i))) &
85          *(1.+RETV*q1(i))
86     zri(i) = zgeop1(i)*(ztvd-ztsolv)/(zdu2*ztvd)
87     zcdn(i) = (ckap/LOG(1.+zgeop1(i)/(RG*rugos(i))))**2
88
[524]89!!$        IF (zri(i) .ge. 0.) THEN      ! situation stable
[1067]90     IF (zri(i) .GT. 0.) THEN      ! situation stable
91        zri(i) = MIN(20.,zri(i))
92        IF (.NOT.zxli) THEN
93           zscf = SQRT(1.+cd*ABS(zri(i)))
[2126]94           FRIV = AMAX1(1. / (1.+2.*CB*zri(i)/ZSCF), f_ri_cd_min)
[1067]95           zcfm1(i) = zcdn(i) * FRIV
[2126]96           FRIH = AMAX1(1./ (1.+3.*CB*zri(i)*ZSCF), f_ri_cd_min )
[524]97!!$  PB          zcfh1(i) = zcdn(i) * FRIH
[1279]98!!$ PB           zcfh1(i) = f_cdrag_stable * zcdn(i) * FRIH
99           zcfh1(i) = f_cdrag_ter * zcdn(i) * FRIH
100           IF(nsrf.EQ.is_oce) zcfh1(i) = f_cdrag_oce * zcdn(i) * FRIH
101!!$ PB
[1067]102           pcfm(i) = zcfm1(i)
103           pcfh(i) = zcfh1(i)
104        ELSE
105           pcfm(i) = zcdn(i)* fsta(zri(i))
106           pcfh(i) = zcdn(i)* fsta(zri(i))
107        ENDIF
108     ELSE                          ! situation instable
109        IF (.NOT.zxli) THEN
110           zucf = 1./(1.+3.0*cb*cc*zcdn(i)*SQRT(ABS(zri(i)) &
111                *(1.0+zgeop1(i)/(RG*rugos(i)))))
[2126]112           zcfm2(i) = zcdn(i)*amax1((1.-2.0*cb*zri(i)*zucf),f_ri_cd_min)
113!!$PB            zcfh2(i) = zcdn(i)*amax1((1.-3.0*cb*zri(i)*zucf),f_ri_cd_min)
114           zcfh2(i) = f_cdrag_ter*zcdn(i)*amax1((1.-3.0*cb*zri(i)*zucf),f_ri_cd_min)
[1067]115           pcfm(i) = zcfm2(i)
116           pcfh(i) = zcfh2(i)
117        ELSE
118           pcfm(i) = zcdn(i)* fins(zri(i))
119           pcfh(i) = zcdn(i)* fins(zri(i))
[524]120        ENDIF
[1067]121        zcr = (0.0016/(zcdn(i)*SQRT(zdu2)))*ABS(ztvd-ztsolv)**(1./3.)
[1279]122        IF(nsrf.EQ.is_oce) pcfh(i) =f_cdrag_oce* zcdn(i)*(1.0+zcr**1.25)**(1./1.25)
[1067]123     ENDIF
124  END DO
125
126! ================================================================= c
127     
128  ! IM cf JLD : on seuille cdrag_m et cdrag_h
129  IF (nsrf == is_oce) THEN
130     DO i=1,knon
131        pcfm(i)=MIN(pcfm(i),cdmmax)
132        pcfh(i)=MIN(pcfh(i),cdhmax)
133     END DO
134  END IF
135
136END SUBROUTINE clcdrag
Note: See TracBrowser for help on using the repository browser.