source: LMDZ6/trunk/libf/phylmd/coare_cp_mod.f90 @ 5444

Last change on this file since 5444 was 5268, checked in by abarral, 2 months ago

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