-
Notifications
You must be signed in to change notification settings - Fork 237
/
cheapaml_lanl_flux.F
217 lines (195 loc) · 7.1 KB
/
cheapaml_lanl_flux.F
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
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
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
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
#include "CHEAPAML_OPTIONS.h"
#undef ALLOW_THSICE
CBOP
C !ROUTINE: CHEAPAML_LANL_FLUX
C !INTERFACE:
SUBROUTINE CHEAPAML_LANL_FLUX(
I i,j,bi,bj,
O fsha, flha, evp, xolw, ssqt, q100 )
C !DESCRIPTION:
C ==================================================================
C SUBROUTINE cheapaml_LANL_flux
C ==================================================================
C o compute surface fluxes using LANL algorithms
C !USES:
IMPLICIT NONE
C === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "DYNVARS.h"
#include "GRID.h"
c#include "FFIELDS.h"
c#ifdef ALLOW_THSICE
c#include "THSICE_VARS.h"
c#endif
#include "CHEAPAML.h"
C !INPUT PARAMETERS:
INTEGER i,j,bi,bj
C !OUTPUT PARAMETERS:
_RL fsha, flha, evp, xolw, ssqt, q100
CEOP
C Output:
C ustress, vstress - wind stress
C fsha - sensible heat flux
C flha - latent heat flux
C xolw - oceanic upwelled long wave radiation
C ssqt - sat. specific humidity at atm layer top
C Input:
C uwind, vwind - mean wind speed (m/s)
C Tair - mean air temperature (K) at height ht (m)
C theta(k=1) - sea surface temperature (C)
C Qair - Specific humidity kg/kg
C Solar - short wave net solar flux at surface (W/m^2)
C Tr - relaxation profile for temperature on boundaries (C)
C qr - relaxation profile for specific humidity (kg/kg)
C i,j,bi,bj - indices of data
C !LOCAL VARIABLES:
C iceornot :: variables to include seaice effect
INTEGER iceornot
_RL deltaTm
_RL uss,usm,uw,vw
_RL cheapaml_BulkCdn
_RL to
_RL t
_RL t0,QaR
_RL ssq, q
_RL deltap, delq, pt, psx100, z100ol
_RL rdn,ren,rhn,zice,zref
_RL rd,re,rh,tta,ttas,toa,ttt
_RL ustar,tstar,qstar,ht,hu,hq
_RL aln,cdalton,czol,psim_fac
_RL huol,stable,xsq,x,psimh,psixh
_RL clha, csha
INTEGER niter_bulk,iter
C useful values
C hardwire atmospheric relative humidity at 80%
QaR=0.8 _d 0
C factor to compute rainfall from specific humidity
C inverse of time step
deltaTm=1. _d 0/deltaT
C reference values to compute turbulent flux
ht=zt
hq=zq
hu=zu
zref = zt
zice=.0005 _d 0
aln = log(ht/zref)
C for iterating on turbulence
niter_bulk = 5
cdalton = 0.0346000 _d 0
czol = zref*xkar*gravity
psim_fac=5. _d 0
C determine wind stress
IF(.NOT.useStressOption)THEN
if (maskC(i,j,1,bi,bj).ne.0. _d 0) then
#ifdef ALLOW_THSICE
if (ICEMASK(i,j,bi,bj).gt.0. _d 0) then
if (snowheight(i,j,bi,bj).gt.3. _d -1) then
iceornot=2
else
iceornot=1
endif
else
iceornot=0
endif
#else
iceornot=0
#endif
uw=uwind(i,j,bi,bj)
vw=vwind(i,j,bi,bj)
uss=sqrt(uw**2+vw**2)
usm=max(uss,1. _d 0)
cheapaml_BulkCdn = cdrag_1/usm + cdrag_2 + cdrag_3*usm
ustress(i,j,bi,bj)= rhoa*cheapaml_BulkCdn*uss*uw
vstress(i,j,bi,bj)= rhoa*cheapaml_BulkCdn*uss*vw
else
usm=0. _d 0
ustress(i,j,bi,bj) = 0. _d 0
vstress(i,j,bi,bj) = 0. _d 0
endif
C wind stress computed
ENDIF
C diabatic and freshwater flux forcing
to=theta(i,j,1,bi,bj)
t=Tair(i,j,bi,bj)
toa=to+Celsius2K
tta=t+Celsius2K
ttas=tta+gamma_blk*zref
ttt=tta-( cheaphgrid(i,j,bi,bj)- zref)*gamma_blk
pt=p0*(1-gamma_blk*cheaphgrid(i,j,bi,bj)/ttas)
& **(gravity/gamma_blk/gasR)
C specific humidities
ssq= ssq0*exp( lath*(ssq1-ssq2/toa) ) / p0
ssqt = ssq0*exp( lath*(ssq1-ssq2/ttt) ) / pt
C saturation no more at the top:
ssqt = 0.7 _d 0*ssq
if (useFreshWaterFlux) then
q=qair(i,j,bi,bj)
else
q=QaR * ssq
endif
C adjust temperature from reference height to formula height
deltap = t - to + gamma_blk*(zref-ht)
delq = q - ssq
ttas = tta+gamma_blk*(zref-ht)
t0 = ttas*(1. _d 0 + humid_fac*q)
C initialize estimate exchange coefficients
rdn=xkar/(log(zref/zice))
rhn=rdn
ren=rdn
C calculate turbulent scales
ustar=rdn*usm
tstar=rhn*deltap
qstar=ren*delq
C iteration with psi-functions to find transfer coefficients
do iter=1,niter_bulk
huol = czol/ustar**2 *(tstar/t0 +
& qstar/(1. _d 0/humid_fac+q))
huol = sign( min(abs(huol),10. _d 0), huol)
stable = 5. _d -1 + sign(5. _d -1 , huol)
xsq = max(sqrt(abs(1. _d 0 - 16. _d 0*huol)),1. _d 0)
x = sqrt(xsq)
psimh = -5. _d 0*huol*stable + (1. _d 0-stable)*
& (2. _d 0*log(5. _d -1*(1. _d 0+x)) +
& 2. _d 0*log(5. _d -1*(1. _d 0+xsq)) -
& 2. _d 0*atan(x) + pi*.5 _d 0)
psixh = -5. _d 0*huol*stable + (1. _d 0-stable)*
& (2. _d 0*log(5. _d -1*(1. _d 0+xsq)))
C Update the transfer coefficients
rd = rdn/(1. _d 0 + rdn*(aln-psimh)/xkar)
rh = rhn/(1. _d 0 + rhn*(aln-psixh)/xkar)
re = rh
C Update ustar, tstar, qstar using updated, shifted coefficients.
ustar = rd*usm
qstar = re*delq
tstar = rh*deltap
enddo
usm=max(uss,0.5 _d 0)
csha = rhoa*cpair*usm*rh*rd
clha = rhoa*lath*usm*re*rd
fsha = csha*deltap
flha = clha*delq
evp = -flha/lath
C the sensible and latent heat fluxes, fsha and flha,
C are computed so that positive values are downward.
C the convention for cheapaml is upward fluxes are positive,
C so they must be multiplied by -1
fsha=-fsha
flha=-flha
C oceanic upwelled long wave
xolw=stefan*(toa)**4
C compute specific humidity at 100m
huol = czol/ustar**2 *(tstar/t0 +
& qstar/(1. _d 0/humid_fac+q))
huol = sign( min(abs(huol),10. _d 0), huol)
stable = 5. _d -1 + sign(5. _d -1 , huol)
z100ol = 100. _d 0 *xkar*gravity/ustar**2 *(tstar/t0
& + qstar/(1. _d 0/humid_fac+q))
xsq = max(sqrt(abs(1. _d 0 - 16. _d 0*z100ol)),1. _d 0)
x = sqrt(xsq)
psx100 = -5. _d 0*z100ol*stable + (1. _d 0-stable)*
& (2. _d 0*log(5. _d -1*(1. _d 0+xsq)))
q100=ssq+qstar*(dlog(100. _d 0/zice)-psx100)
RETURN
END