source: LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmd/coare_cp.F90 @ 5072

Last change on this file since 5072 was 4661, checked in by Laurent Fairhead, 15 months ago

Inclusion and merge of Olivier's modifications
File size: 6.4 KB
Line 
1real function psit_30(zet)
2
3  real, intent(in) :: zet
4
5  if(zet<0) then
6     x=(1.-(15*zet))**.5
7     psik=2*log((1+x)/2)
8     x=(1.-(34.15*zet))**.3333
9     psic=1.5*log((1.+x+x*x)/3.)-sqrt(3.)*atan((1.+2.*x)/sqrt(3.))+4.*atan(1.)/sqrt(3.)
10     f=zet*zet/(1+zet*zet)
11     psit_30=(1-f)*psik+f*psic   
12
13  else
14     c=min(50.,.35*zet)
15     psit_30=-((1.+2./3.*zet)**1.5+.6667*(zet-14.28)/exp(c)+8.525)
16  endif
17end FUNCTION psit_30
18
19real function psiuo(zet)
20
21  real, intent(in) :: zet
22
23  if (zet<0) then
24
25     x=(1.-15.*zet)**.25
26     psik=2.*log((1.+x)/2.)+log((1.+x*x)/2.)-2.*atan(x)+2.*atan(1.)
27     x=(1.-10.15*zet)**.3333
28     psic=1.5*log((1.+x+x*x)/3.)-sqrt(3.)*atan((1.+2.*x)/sqrt(3.))+4.*atan(1.)/sqrt(3.)
29     f=zet*zet/(1+zet*zet)
30     psiuo=(1-f)*psik+f*psic                                               
31  else
32     c=min(50.,.35*zet)
33     psiuo=-((1+1.0*zet)**1.0+.667*(zet-14.28)/exp(c)+8.525)
34  endif
35END FUNCTION psiuo
36
37
38real function gen_q1(ts,p,dq)
39
40  real, intent(in) :: ts,p,dq
41  real es,qs
42
43  es = 6.112*exp(17.502*(ts-273.15)/(ts-32.18))*.98*(1.0007+3.46e-8*p)
44  qs = es*62.197/(p-37.8*es)
45
46  q1 = dq + qs 
47
48end function gen_q1
49
50
51
52
53subroutine coare_cp(du,dt,dq,&
54     t,q,&
55     zu,zt,zq,&
56     p,zi,&
57     nits,&
58     coeffs,rugosm,rugosh)
59  !     zo,tau,hsb,hlb,&
60  !     var)
61  !version with shortened iteration  modified Rt and Rq
62  !uses wave information wave period in s and wave ht in m
63  !no wave, standard coare 2.6 charnock:  jwave=0
64  !Oost et al.  zo=50/2/pi L (u*/c)**4.5 if jwave=1
65  !taylor and yelland  zo=1200 h*(L/h)**4.5 jwave=2
66  !
67  USE MODD_CSTS,       ONLY : XKARMAN, XG, XSTEFAN, XRD, XRV, XPI, &
68                            XLVTT, XCL, XCPD, XCPV, XRHOLW, XTT, &
69                            XP00
70
71  IMPLICIT NONE
72
73  real, intent(in) :: du,dt,dq,t,q
74  real, intent(in) :: zu,zt,zq,p,zi
75  integer, intent(in) :: nits
76  ! real, dimension (nits), intent(out) :: zo,tau,hsb,hlb
77  ! real, dimension(nits,3), intent(out) :: var
78  real, dimension(3), intent(out) :: coeffs
79  real, intent(out) :: rugosm
80  real, intent(out) :: rugosh
81
82  real, parameter ::  beta=1.2, von=.4, fdg = 1. ,&
83       tdk = 273.16, pi = 3.141593, grav = 9.82,&
84       rgas = 287.1
85
86  integer, dimension(3) :: shape_input
87
88  real bf, cc, rhoa, visa,&
89       u10, ut, uts, ut0, ug,&
90       cd10, ch10, ct, ct10,&
91       charn, ribu,&
92       l, l10, &
93       ribcu,rr,&
94       zet,zetu, zom10, zoh10, zot, zoq, zot10,&
95       cd, ch, le, cpa, cpv,&
96       usr,qsr,tsr,&
97       zom, t_c
98
99!, ZTWAVE, ZHWAVE, ZCWAVE, ZLWAVE
100
101
102
103
104  real old_usr, old_tsr, old_qsr,tmp
105
106  real, external :: psit_30, psiuo, grv
107  integer i,j,k
108!---------------- Rajout pour prendre en compte différent Z0 --------------------------------!
109!  INTEGER :: NGRVWAVES        ! Pour le choix du z0
110!  NGRVWAVES = 2
111!---------------------------------------------------------------------------------------------
112!-------------------- Attention Modif réalisée pas SURRR -------------------------------------
113!---------------------------------------------------------------------------------------------
114
115  Ribcu=-zu/zi/.004/Beta**3
116  cpa=1004.67
117
118  t_c = t - 273.15
119
120
121  Le=(2.501-.00237*(t_c-dt))*1e6
122
123
124  visa=1.326e-5*(1+6.542e-3*(t_c)+8.301e-6*t_c**2-4.84e-9*t_c**3)
125
126
127  cpv=cpa*(1+0.84*Q)
128  rhoa=P/(Rgas*t*(1+0.61*Q))
129
130
131
132
133
134
135  ug=.2
136
137  ut=sqrt(du*du+ug*ug)
138  ut=MAX(ut , 0.1 * MIN(10.,zu) )
139
140  u10=ut*log(10/1e-4)/log(zu/1e-4)
141  usr=.035*u10                              ! turbulent friction velocity (m/s), including gustiness
142  zom10=0.011*usr*usr/grav+0.11*visa/usr     ! roughness length for u (smith 88)
143  Cd10=(von/log(10/zom10))**2
144!  zoh10=0.40*visa/usr+1.4e-5
145!  Ch10=((von**2)/((log(10/zom10))*(log(10/zoh10)))) !ammener à devenir ca
146  Ch10=0.00115
147  Ct10=Ch10/sqrt(Cd10)
148  zot10=10/exp(von/Ct10)                    ! roughness length for t
149  Cd=(von/log(zu/zom10))**2
150  Ct=von/log(zt/zot10)
151  CC=von*Ct/Cd
152
153  ut0 = ut
154
155
156  ut = ut0
157  Ribu=grav*zu/t*(dt+.61*t*dq)/ut**2
158
159  if (Ribu .LT. 0) then
160     zetu=CC*Ribu/(1+Ribu/Ribcu)
161  else
162     zetu=CC*Ribu*(1+27/9*Ribu/CC)
163  endif
164
165  L10=zu/zetu
166
167  ! if (zetu .GT. 50) then
168  !    nits=1
169  ! endif
170
171  usr=ut*von/(log(zu/zom10)-psiuo(zetu))
172  tsr=dt*von*fdg/(log(zt/zot10)-psit_30(zt/L10))
173  qsr=dq*von*fdg/(log(zq/zot10)-psit_30(zq/L10))
174  !  tkt=.001
175
176  ! charnock constant - lin par morceau - constant
177  if ( ut <= 10. ) then
178     charn=0.011
179  else
180     if (ut .GT. 18) then
181        charn=0.018
182     else
183        charn=0.011+(ut-10)/(18-10)*(0.018-0.011)
184     endif
185  endif
186
187!  ZHWAVE = 0.018*ut*ut*(1.+0.015*ut)
188!  ZTWAVE = 0.729*ut
189!  ZCWAVE = XG*ZTWAVE/(2.*pi)
190!  ZLWAVE = ZTWAVE*ZCWAVE
191
192
193  !***************  bulk loop ************
194  do i=1, nits
195
196     zet=von*grav*zu/t*(tsr*(1+0.61*Q)+.61*t*qsr)/(usr*usr)/(1+0.61*Q)
197
198!     IF (NGRVWAVES==0) THEN
199!       zom = charn*usr*usr/XG + 0.11*visa/usr !Smith 1988
200!     ELSE IF (NGRVWAVES==1) THEN
201!       zom = (50./(2.*pi))*ZLWAVE*(usr/ZCWAVE)**4.5 &
202!               + 0.11*visa/usr                       !Oost et al. 2002 
203!     ELSE IF (NGRVWAVES==2) THEN
204!       zom = 1200.*ZHWAVE*(ZHWAVE/ZLWAVE)**4.5 &
205!               + 0.11*visa/usr                       !Taulor and Yelland 2001 
206!     ENDIF
207
208     zom=charn*usr*usr/grav+0.11*visa/usr
209
210     rr=zom*usr/visa
211     L=zu/zet
212     zoq=min(1.15e-4,5.5e-5/rr**.6)  ! a modifier
213     zot=zoq                         ! a modifier
214!     zot=0.40*visa/usr+1.4e-5
215
216     old_usr = usr
217     old_tsr = tsr
218     old_qsr = qsr
219
220     usr=ut*von/(log(zu/zom)-psiuo(zu/L))
221     tsr=dt*von*fdg/(log(zt/zot)-psit_30(zt/L))
222     qsr=dq*von*fdg/(log(zq/zoq)-psit_30(zq/L))
223
224     Bf=-grav/t*usr*(tsr+.61*t*qsr)
225
226     if (Bf .GT. 0) then
227        ug=Beta*(Bf*zi)**.333
228     else
229        ug=.2
230     endif
231
232     ut=sqrt(du*du+ug*ug)
233     ut=MAX(ut , 0.1 * MIN(10.,zu) )
234
235
236  enddo !bulk iter loop
237
238  ! coeffs(m1,m2,m3,m4,m5,1)=rhoa*usr*usr*du/ut !stress
239  ! coeffs(m1,m2,m3,m4,m5,2)=rhoa*cpa*usr*tsr
240  ! coeffs(m1,m2,m3,m4,m5,3)=rhoa*Le*usr*qsr
241
242  tmp = (von/(log(zu/zom)-psiuo(zu/L)) ) * ut / du
243
244  coeffs = [tmp**2,&
245       von*fdg/(log(zt/zot)-psit_30(zt/L)) * tmp,&
246       von*fdg/(log(zq/zoq)-psit_30(zq/L)) * tmp]
247
248  rugosm = zom
249  rugosh = zot
250
251
252end subroutine coare_cp
Note: See TracBrowser for help on using the repository browser.