Newer
Older
MODULE drag3
USE param
USE inputdata
type, public:: data_in
real, public:: ws, dt, st, dq, cflh, z0in
end type
type, public:: data_outdef
real, public:: zl, ri, re, lnzuzt, zu, ztout, rith, cm, ch, ct, ckt
end type
type, public:: data_par
integer, public :: it=10
end type
type, public:: data_lutyp
integer, public :: lu_indx=1
end type
! SUBROUTINE surf_fluxMAS(mas1_w, mas1_dt, mas1_st, mas1_dq, out, par1, lu1)
SUBROUTINE surf_fluxMAS(mas1_w, mas1_dt, mas1_st, mas1_dq, mas1_cflh, mas1_z0in, &
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
masout_zl, masout_ri, masout_re, masout_lnzuzt, masout_zu, masout_ztout,&
masout_rith, masout_cm, masout_ch, masout_ct, masout_ckt,&
par1, lu1,numst)
real, dimension (numst) :: mas1_w
real, dimension (numst) :: mas1_dt
real, dimension (numst) :: mas1_st
real, dimension (numst) :: mas1_dq
real, dimension (numst) :: mas1_cflh
real, dimension (numst) :: mas1_z0in
real, dimension (numst) :: masout_zl
real, dimension (numst) :: masout_ri
real, dimension (numst) :: masout_re
real, dimension (numst) :: masout_lnzuzt
real, dimension (numst) :: masout_zu
real, dimension (numst) :: masout_ztout
real, dimension (numst) :: masout_rith
real, dimension (numst) :: masout_cm
real, dimension (numst) :: masout_ch
real, dimension (numst) :: masout_ct
real, dimension (numst) :: masout_ckt
type (data_in) in
type (data_outdef) out
type (data_par) par1
type (data_lutyp) lu1
integer i,numst
do i=1,numst
in%ws=mas1_w(i)
in%dt=mas1_dt(i)
in%st=mas1_st(i)
in%dq=mas1_dq(i)
in%cflh=mas1_cflh(i)
in%z0in=mas1_z0in(i)
!if ((in%ws == in%ws).and.(in%dt == in%dt).and.(in%st == in%st).and.(in%dq == in%dq)) then
CALL surf_flux(in, out, par1, lu1)
!end if
masout_zl(i)=out%zl
masout_ri(i)=out%ri
masout_re(i)=out%re
masout_lnzuzt(i)=out%lnzuzt
masout_zu(i)=out%zu
masout_ztout(i)=out%ztout
masout_rith(i)=out%rith
masout_cm(i)=out%cm
masout_ch(i)=out%ch
masout_ct(i)=out%ct
masout_ckt(i)=out%ckt
enddo
END SUBROUTINE surf_fluxMAS
SUBROUTINE surf_flux(in, out, par, lu)
type (data_in) , intent(in) :: in
type (data_outdef) out
type (data_par) par
type (data_lutyp) lu
real ws, dt, dq, cflh, z0in
integer it
integer lu_indx
real zl, ri, re, lnzuzt, zu, ztin, ri_th, cm, ch, ct, ckt
real z0, d3, d0max, U10m, a1, y1, cimin, f, c1, u2, h0, u3, x7, d0, d00, zt, h00, ft0, an4, an5
real al, Tsurf, Rib, q4, t4, u, r1, f0, f4, am, o, dd, x1, y0, x0, y10, a2ch, x10, p1, p0, h, d1, f1
real d, c4, c1min, c0, c, b1, an0
integer i, j
ws=in%ws
dt=in%dt
Tsurf=in%st
dq=in%dq
cflh=in%cflh
z0in=in%z0in
it=par%it
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
u=ws
t4=dt
q4=dq
h=cflh
z0=z0in
AN5=(Pr_t_inf_inv / Pr_t_0_inv)**4
D1=(2.0E0*alpha_h-AN5*alpha_m-SQRT((AN5*alpha_m)**2+4.0E0*AN5*alpha_h*(alpha_h-alpha_m)))/(2.0E0*alpha_h**2)
Y10=(1.0E0-alpha_m*D1)**.25E0
X10=(1.0E0-alpha_h*D1)**.5E0
P1=2.0E0*ATAN(Y10)+ALOG((Y10-1.0E0)/(Y10+1.0E0))
P0=ALOG((X10-1.0E0)/(X10+1.0E0))
d3=0.0e0
d0max=2.0e0
if(z0 < 0.0e0) then
!> @brief Test - definition z0 of sea surface
!> @details if lu_indx=2.or.lu_indx=3 call z0sea module (Ramil_Dasha)
d0max = 8.0e0
U10m=u
a1=0.0e0
y1=25.0e0
c1min=alog(h1/1.0e0)/kappa
do i=1,it
f=a2-2.0e0*alog(U10m)
do j=1,it
c1=(f+2.0e0*alog(y1))/kappa
! looks like the check should use U10 instead of U
! but note that a1 should be set = 0 in cycle beforehand
if(u <= 8.0e0) a1=alog(1.0e0+a3*((y1/U10m)**3))/kappa
c1=c1-a1
c1=max(c1,c1min)
y1=c1
end do
z0=h1*exp(-c1*kappa)
z0=max(z0,0.000015e0)
U10m=u*alog(h1/z0)/(alog(h/z0))
end do
h0=h/z0
u3=U10m/c1
else
! ......parameters from viscosity sublayer......
h0=h/z0
u3=u*kappa/alog(h0)
end if
x7=u3*z0/an
if(x7 <= x8) then
d0=an1*alog(al1*x7)+an2
else
d0=al2*(x7**0.45e0)
end if
! ......humidity stratification and ri-number......
al=g/Tsurf
d0=min(d0,d0max)
Rib=al*h*(t4+0.61e0*Tsurf*q4)/u**2
d00=d0
zt=z0/exp(d00)
h00=h/zt
ft0=alog(h00)
! ......definition of r-prim......
an4=d1/h0
an5=d1/h00
if (abs(d0) < 1.0e-10) an5=an4
an5=sqrt(1.0e0-g0*an5)
an4=(1.0e0-alpha_m*an4)**0.25e0
f0=alog((x10-1.0e0)*(an5+1.0e0)/((x10+1.0e0)*(an5-1.0e0)))/Pr_t_0_inv
f4=2.0e0*(atan(y10)-atan(an4))+alog((y10-1.0e0)*(an4+1.0e0)/((y10+1.0e0)*(an4-1.0e0)))
r1=d1*f0/(f4*f4)
! ......definition of dz,ta,fu,fo,fiu,fio......
if (Rib > 0.0e0) then
! ......stable stratification......
!write (*,*) 'stable'
Rib=min(Rib,r0)
f=alog(h0)
f1=d0/f
a1=b4*Rib
a2ch=(f1+1.0e0)/Pr_t_0_inv-2.0e0*a1
d3=f*(sqrt(a2ch**2+4.0e0*a1*(1.0e0-a1))-a2ch)/(2.0e0*b4*(1.0e0-a1))
f1=b4*d3
f4=f+f1
f0=(f+d0)/Pr_t_0_inv+f1
o=1.0e0/Pr_t_0_inv+f1
am=1.0e0+f1
else if (Rib < r1) then
! ......strong instability.....
!write (*,*) 'instability'
d3=d1
do i=1,it+1
d=d3/h0
dd=d3/h00
if (abs(d0) < 1.0e-10) dd=d
a1=(d1/d3)**(1.0e0/3.0e0)
x0=sqrt(1.0e0-g0*dd)
y0=(1.0e0-alpha_m*d)**0.25e0
c=alog((x0+1.0e0)/(x0-1.0e0))
b1=-2.0e0*atan(y0)+alog((y0+1.0e0)/(y0-1.0e0))
f=3.0e0*(1.0e0-a1)
f4=f/y10+p1+b1
f0=(f/x10+p0+c)/Pr_t_0_inv
if (i == it+1) exit
d3=Rib*f4**2/f0
end do
am=a1/y10
o=a1/(Pr_t_0_inv*x10)
else if (Rib > -0.001e0) then
! ......nearly neutral......
write (*,*) 'neutral'
f4=alog(h0)
f0=ft0/Pr_t_0_inv
if (abs(d0) < 1.0e-10) f0=f4/Pr_t_0_inv
am=1.0e0
o=1.0e0/Pr_t_0_inv
else
! ......week and semistrong instability......
!write (*,*) 'semistrong'
f1=alog(h0)
if (abs(d0) < 1.0e-10) ft0=f1
d3=Rib*Pr_t_0_inv*f1**2/ft0
do i=1,it+1
d=d3/h0
dd=d3/h00
if (abs(d0) < 1.0e-10) dd=d
y1=(1.0e0-alpha_m*d3)**0.25e0
x1=sqrt(1.0e0-g0*d3)
y0=(1.0e0-alpha_m*d)**0.25e0
x0=sqrt(1.0e0-g0*dd)
y0=max(y0,1.000001e0)
x0=max(x0,1.000001e0)
f4=alog((y1-1.0e0)*(y0+1.0e0)/((y1+1.0e0)*(y0-1.0e0)))+2.0e0*(atan(y1)-atan(y0))
f0=alog((x1-1.0e0)*(x0+1.0e0)/((x1+1.0e0)*(x0-1.0e0)))/Pr_t_0_inv
if (i == it + 1) exit
d3=Rib*f4**2/f0
end do
am=(1.0e0-alpha_m*d3)**(-0.25e0)
o=1.0e0/(Pr_t_0_inv*sqrt(1.0e0-g0*d3))
end if
! ......computation of cu,co,k(h),alft
c4=kappa/f4
c0=kappa/f0
an4=kappa*c4*u*h/am
an0=am/o
! ......exit......
out%zl=d3
out%ri=Rib
out%re=x7
out%lnzuzt=d00
out%zu=z0
out%ztout=zt
out%rith=r1
out%cm=c4
out%ch=c0
out%ct=an4
out%ckt=an0
return
END SUBROUTINE surf_flux
END MODULE drag3