source: dynamico_lmdz/simple_physics/phyparam/param/vdif_cd.F @ 4183

Last change on this file since 4183 was 4176, checked in by dubos, 6 years ago

simple_physics : copy code from FH

File size: 2.6 KB
Line 
1      SUBROUTINE vdif_cd( ngrid,pz0,pg,pz,pu,pv,pts,ph,pcdv,pcdh)
2      IMPLICIT NONE
3c=======================================================================
4c
5c   Subject: computation of the surface drag coefficient using the
6c   -------  approch developed by Loui for ECMWF.
7c
8c   Author: Frederic Hourdin  15 /10 /93
9c   -------
10c
11c   Arguments:
12c   ----------
13c
14c   inputs:
15c   ------
16c     ngrid            size of the horizontal grid
17c     pg               gravity (m s -2)
18c     pz(ngrid)        height of the first atmospheric layer
19c     pu(ngrid)        u component of the wind in that layer
20c     pv(ngrid)        v component of the wind in that layer
21c     pts(ngrid)       surfacte temperature
22c     ph(ngrid)        potential temperature T*(p/ps)^kappa
23c
24c   outputs:
25c   --------
26c     pcdv(ngrid)      Cd for the wind
27c     pcdh(ngrid)      Cd for potential temperature
28c
29c=======================================================================
30c
31c-----------------------------------------------------------------------
32c   Declarations:
33c   -------------
34
35c   Arguments:
36c   ----------
37
38      INTEGER ngrid,nlay
39      REAL pz0(ngrid)
40      REAL pg,pz(ngrid)
41      REAL pu(ngrid),pv(ngrid)
42      REAL pts(ngrid),ph(ngrid)
43      REAL pcdv(ngrid),pcdh(ngrid)
44
45c   Local:
46c   ------
47
48      INTEGER ig
49
50      REAL zu2,z1,zri,zcd0,zz
51
52      REAL karman,b,c,d,c2b,c3bc,c3b,z0,umin2
53      LOGICAL firstcal
54      DATA karman,b,c,d,umin2/.4,5.,5.,5.,1.e-12/
55      DATA firstcal/.true./
56      SAVE b,c,d,karman,c2b,c3bc,c3b,firstcal,umin2
57
58c-----------------------------------------------------------------------
59c   couche de surface:
60c   ------------------
61
62c     DO ig=1,ngrid
63c        zu2=pu(ig)*pu(ig)+pv(ig)*pv(ig)+umin2
64c        pcdv(ig)=pz0(ig)*(1.+sqrt(zu2))
65c        pcdh(ig)=pcdv(ig)
66c     ENDDO
67c     RETURN
68
69      IF (firstcal) THEN
70         c2b=2.*b
71         c3bc=3.*b*c
72         c3b=3.*b
73         firstcal=.false.
74      ENDIF
75
76c!!!! WARNING, verifier la formule originale de Louis!
77      DO ig=1,ngrid
78         zu2=pu(ig)*pu(ig)+pv(ig)*pv(ig)+umin2
79         zri=pg*pz(ig)*(ph(ig)-pts(ig))/(ph(ig)*zu2)
80         z1=1.+pz(ig)/pz0(ig)
81         zcd0=karman/log(z1)
82         zcd0=zcd0*zcd0*sqrt(zu2)
83         IF(zri.LT.0.) THEN
84            z1=b*zri/(1.+c3bc*zcd0*sqrt(-z1*zri))
85            pcdv(ig)=zcd0*(1.-2.*z1)
86            pcdh(ig)=zcd0*(1.-3.*z1)
87         ELSE
88            zz=sqrt(1.+d*zri)
89            pcdv(ig)=zcd0/(1.+c2b*zri/zz)
90            pcdh(ig)=zcd0/(1.+c3b*zri*zz)
91         ENDIF
92      ENDDO
93
94c-----------------------------------------------------------------------
95
96      RETURN
97      END
Note: See TracBrowser for help on using the repository browser.