/[MITgcm]/MITgcm/pkg/bulk_force/bulkf_formula_lanl.F
ViewVC logotype

Contents of /MITgcm/pkg/bulk_force/bulkf_formula_lanl.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph

Revision 1.3 - (show annotations) (download)
Thu Oct 9 04:19:19 2003 UTC (20 years, 8 months ago) by edhill
Branch: MAIN
CVS Tags: checkpoint51k_post, checkpoint51o_pre, checkpoint51l_post, checkpoint52, checkpoint51t_post, checkpoint51n_post, checkpoint51s_post, checkpoint51n_pre, checkpoint52b_pre, checkpoint51l_pre, checkpoint51q_post, checkpoint51r_post, checkpoint51i_post, checkpoint52a_pre, checkpoint51o_post, checkpoint52a_post, ecco_c52_e35, checkpoint51m_post, checkpoint51p_post, checkpoint51u_post
Branch point for: branch-nonh, tg2-branch, checkpoint51n_branch
Changes since 1.2: +3 -1 lines 
 o first check-in for the "branch-genmake2" merge
 o verification suite as run on shelley (gcc 3.2.2):

Wed Oct  8 23:42:29 EDT 2003
                T           S           U           V
G D M    c        m  s        m  s        m  s        m  s
E p a R  g  m  m  e  .  m  m  e  .  m  m  e  .  m  m  e  .
N n k u  2  i  a  a  d  i  a  a  d  i  a  a  d  i  a  a  d
2 d e n  d  n  x  n  .  n  x  n  .  n  x  n  .  n  x  n  .

OPTFILE=NONE

Y Y Y Y 13 16 16 16  0 16 16 16 16 16 16 16 16 13 12  0  0 pass  adjustment.128x64x1
Y Y Y Y 16 16 16 16  0 16 16 16 16 16 16  0  0 16 16  0  0 pass  adjustment.cs-32x32x1
Y Y Y Y 16 16 16 16  0 16 16 16 16 16 16 22  0 16 16 22  0 pass  adjust_nlfs.cs-32x32x1
Y Y Y Y -- 13 13 16 16 13 13 13 13 16 16 16 16 16 16 16 16 N/O   advect_cs
Y Y Y Y -- 22 16 16 16 16 16 16 13 16 16 16 16 16 16 16 16 N/O   advect_xy
Y Y Y Y -- 13 16 13 16 16 16 16 16 16 16 22 16 16 16 16 16 N/O   advect_xz
Y Y Y Y 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 pass  aim.5l_cs
Y Y Y Y 14 16 16 16 16 16 16 16 16 13 16 16 16 16 16 13 16 pass  aim.5l_Equatorial_Channel
Y Y Y Y 16 16 16 16 16 16 16 16 16 16 16 13 16 16 13 13 16 pass  aim.5l_LatLon
Y Y Y Y 13 16 16 16 16 16 16 16 16 16 13 12 13 13 16 13 16 pass  exp0
Y Y Y Y 14 16 16 16 16 16 16 16 22 16 16 16 13 16 16 22 16 pass  exp1
Y Y Y Y 13 13 16 13 16 16 16 16 16 13 13 16 16 13 13 13 13 pass  exp2
Y Y Y Y 16 16 16 16 16 16 16 16 22 16 16 16 16 16 16 16 16 pass  exp4
Y Y Y Y 16 16 16 16 16 16 16 16 16 16 16 22 16 16 16 22 16 pass  exp5
Y Y Y Y 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 pass  front_relax
Y Y Y Y 14 16 16 13 13 16 16 13 13 16 13 13 16 12 13 13 16 pass  global_ocean.90x40x15
Y Y Y Y 10 16 16 13 13 16 13 16 16 13 13 13 13 16 16 13 16 FAIL  global_ocean.cs32x15
Y Y Y Y  6 11 12 13 13 12 13 16 13  9  9  9  9 10  9  9 11 FAIL  global_ocean_pressure
Y Y Y Y 14 16 16 13 16 16 16 13 13 13 13 13 16 12 16 13 16 pass  global_with_exf
Y Y Y Y 14 16 16 16 16 16 16 16 16 11 13 22 13 16 16  9 16 pass  hs94.128x64x5
Y Y Y Y 13 16 16 16 16 16 16 16 16 11 16 16 16 13 16 22 13 pass  hs94.1x64x5
Y Y Y Y 14 16 16 16 16 16 16 16 16 13 16 13 13 16 16 22 13 pass  hs94.cs-32x32x5
Y Y Y Y 10 10 16 13 13 16 16 16 22 16 13 13 13 13 13 22 13 FAIL  ideal_2D_oce
Y Y Y Y  8 16 16 16 16 16 16 16 16 13 13  8 16 16 16 16 16 FAIL  internal_wave
Y Y Y Y 14 16 16 16 16 16 16 16 16 13 13 22 13 13 13 22 16 pass  inverted_barometer
Y Y Y Y 12 16 16 16 16 16 16 16 16 16 13 12 13 13 13 13 13 FAIL  lab_sea
Y Y Y Y 11 16 16 16 16 16 16 16 13 13 13 12 13 16 13 12 13 FAIL  natl_box
Y Y Y Y 16 16 16 16 16 16 16 16 22 16 16 16 16 16 16 16 16 pass  plume_on_slope
Y Y Y Y 13 16 16 16 16 13 16 16 16 16 16 16 16 13 16 16 16 pass  solid-body.cs-32x32x1
1 C $Header: /u/u3/gcmpack/MITgcm/pkg/bulk_force/bulkf_formula_lanl.F,v 1.2.4.1 2003/10/02 18:18:33 adcroft Exp $
2 C $Name: $
3 c swd -- bulkf formula used in bulkf and ice pkgs
4
5 c taken from exf package
6 #include "BULK_FORCE_OPTIONS.h"
7 subroutine bulkf_formula_lanl(
8 I uw, vw, us, Ta, Qa, nc, tsf_in,
9 I flwupa, flha, fsha, df0dT,
10 I ust, vst, evp, ssq, iceornot, windread
11 & )
12
13 IMPLICIT NONE
14 c
15 #include "EEPARAMS.h"
16 #include "SIZE.h"
17 #include "PARAMS.h"
18 #include "DYNVARS.h"
19 #include "GRID.h"
20 #include "FFIELDS.h"
21 #ifdef ALLOW_BULK_FORCE
22 #include "BULKF_ICE_CONSTANTS.h"
23 #include "BULKF.h"
24 #endif
25
26 c
27 c calculate bulk formula fluxes over open ocean or seaice
28 c
29 c input
30 _RL us ! wind speed
31 _RL uw ! zonal wind speed (at grid center)
32 _RL vw ! meridional wind speed (at grid center)
33 _RL Ta ! air temperature at ht
34 _RL Qa ! specific humidity at ht
35 _RL tsf_in ! surface temperature (either ice or open water)
36 _RL nc ! fraction cloud cover
37 integer iceornot ! 0=open water, 1=ice cover
38 logical windread !
39 c output
40 _RL flwupa ! upward long wave radiation
41 _RL flha ! latent heat flux
42 _RL fsha ! sensible heat flux
43 _RL df0dT ! derivative of heat flux with respect to Tsf
44 _RL ust ! zonal wind stress (at grid center)
45 _RL vst ! meridional wind stress (at grid center)
46 _RL evp ! evaporation rate (over open water)
47 _RL ssq ! surface specific humidity (kg/kg)
48 c
49 c local variables
50 _RL tsf ! surface temperature in K
51 _RL ht ! reference height (m)
52 _RL hq ! reference height for humidity (m)
53 _RL hu ! reference height for wind speed (m)
54 _RL zref ! reference height
55 _RL czol
56 _RL usm ! wind speed limited
57 _RL t0 ! virtual temperature (K)
58 _RL deltap ! potential temperature diff (K)
59 _RL delq ! specific humidity difference
60 _RL stable
61 _RL rdn ,ren, rhn
62 _RL ustar
63 _RL tstar
64 _RL qstar
65 _RL huol
66 _RL htol
67 _RL hqol
68 _RL xsq
69 _RL x
70 _RL re
71 _RL rh
72 _RL tau
73 _RL psimh
74 _RL psixh
75 _RL rd
76 _RL uzn
77 _RL shn
78 _RL aln
79 _RL cdalton
80 _RL dflhdT ! derivative of latent heat with respect to T
81 _RL dfshdT ! derivative of sensible heat with respect to T
82 _RL dflwupdT ! derivative of long wave with respect to T
83 _RL mixratio
84 _RL ea
85 _RL psim_fac
86 _RL umin
87 _RL lath
88 _RL csha
89 _RL clha
90 _RL zice
91 integer niter_bulk, iter
92
93 #ifdef ALLOW_BULK_FORCE
94
95 c --- external functions ---
96 _RL exf_BulkCdn
97 external exf_BulkCdn
98 _RL exf_BulkqSat
99 external exf_BulkqSat
100 _RL exf_BulkRhn
101 external exf_BulkRhn
102
103 cQQQQQQQQQQQQQQQQQQQQq
104 c -- compute turbulent surface fluxes
105 ht = 2.d0
106 hq = 2.d0
107 hu = 10.d0
108 zref = 10.d0
109 zice = 0.0005
110 aln = log(ht/zref)
111 niter_bulk = 5
112 cdalton = 0.0346000
113 czol = zref*xkar*gravity
114 psim_fac=5.d0
115 umin=1.d0
116 c
117 lath=Lvap
118 if (iceornot.gt.0) lath=Lvap+Lfresh
119 Tsf=Tsf_in+Tf0kel
120 c Wind speed
121 if (us.eq.0) then
122 us = sqrt(uw*uw + vw*vw)
123 endif
124 usm = max(us,umin)
125 cQQQ try to limit drag
126 cQQ usm = min(usm,5.d0)
127 c
128 t0 = Ta*(1.d0 + humid_fac*Qa)
129 cQQ ssq = 0.622*6.11*exp(22.47*(1.d0-Tf0kel/tsf))/p0
130 ssq =3.797915*exp(lath*(7.93252d-6-2.166847d-3/Tsf))/p0
131 cBB debugging
132 cBB print*,'ice, ssq', ssq
133 c
134 deltap = ta - tsf + gamma_blk*ht
135 delq = Qa - ssq
136 c
137 c initialize estimate exchange coefficients
138 rdn=xkar/(log(zref/zice))
139 rhn=rdn
140 ren=rdn
141 c calculate turbulent scales
142 ustar=rdn*usm
143 tstar=rhn*deltap
144 qstar=ren*delq
145 c
146 c interation with psi-functions to find transfer coefficients
147 do iter=1,niter_bulk
148 huol = czol/ustar**2 *(tstar/t0 +
149 & qstar/(1.d0/humid_fac+Qa))
150 huol = sign( min(abs(huol),10.d0), huol)
151 stable = 5. _d -1 + sign(5. _d -1 , huol)
152 xsq = max(sqrt(abs(1.d0 - 16.d0*huol)),1.d0)
153 x = sqrt(xsq)
154 psimh = -5.d0*huol*stable + (1.d0-stable)*
155 & (2.d0*log(5.e-1*(1.d0+x)) +
156 & 2.d0*log(5.e-1*(1.d0+xsq)) -
157 & 2.d0*atan(x) + pi*.5d0)
158 psixh = -5.d0*huol*stable + (1.d0-stable)*
159 & (2.d0*log(5.e-1*(1.d0+xsq)))
160
161 c Update the transfer coefficients
162
163 rd = rdn/(1.d0 + rdn*(aln-psimh)/xkar)
164 rh = rhn/(1.d0 + rhn*(aln-psixh)/xkar)
165 re = rh
166 c Update ustar, tstar, qstar using updated, shifted coefficients.
167 ustar = rd*usm
168 qstar = re*delq
169 tstar = rh*deltap
170 enddo
171 c
172 tau = rhoa*ustar**2
173 tau = tau*us/usm
174 csha = rhoa*cpair*us*rh*rd
175 clha = rhoa*lath*us*re*rd
176 c
177 fsha = csha*deltap
178 flha = clha*delq
179 evp = flha/lath/rhosw
180 c
181 c up long wave radiation
182 cQQ mixratio=Qa/(1-Qa)
183 cQQ ea=p0*mixratio/(0.62197+mixratio)
184 cQQ flwupa=-0.985*stefan*tsf**4
185 cQQ & *(0.39-0.05*sqrt(ea))
186 cQQ & *(1-0.6*nc**2)
187 if (iceornot.eq.0) then
188 flwupa=-ocean_emissivity*stefan*tsf**4
189 else
190 if (iceornot.eq.2) then
191 flwupa=-snow_emissivity*stefan*tsf**4
192 else
193 flwupa=-ice_emissivity*stefan*tsf**4
194 endif
195 endif
196 #ifdef ALLOW_THERM_SEAICE
197 cswdcou -- remove ---
198 cQQ if (iceornot.gt.0) then
199 cswdcou --- end remove ---
200 c derivatives with respect to Tsf
201 cQQ dflhdT = -clha*Tf0kel*ssq*22.47/(tsf**2)
202 c dflhdT = -clha*Lath*ssq/(Rvap*tsf**2)
203 dflhdT = -clha*ssq*Lath*2.166847E-3/(Tsf**2)
204 dfshdT = -csha
205 if (iceornot.eq.0) then
206 dflwupdT=-4.*ocean_emissivity*stefan*tsf**3
207 else
208 if (iceornot.eq.2) then
209 dflwupdT=-4.*snow_emissivity*stefan*tsf**3
210 else
211 dflwupdT=-4.*ice_emissivity*stefan*tsf**3
212 endif
213 endif
214 cQQ dflwupdT=-4.*0.985*stefan*tsf**3
215 cQQ & *(0.39-0.05*sqrt(ea))
216 cQQ & *(1-0.6*nc**2)
217 c total derivative with respect to surface temperature
218 df0dT=dflwupdT+dfshdT+dflhdT
219 cBB
220 cBB print*,'derivatives:',df0dT,dflwupdT,dfshdT,dflhdT
221 cswdcou -- remove ---
222 cQQ endif
223 cswdcou
224 #endif /*ALLOW_THERM_SEAICE*/
225 c
226 c wind stress at center points
227 if (.NOT.windread) then
228 ust = rhoa*exf_BulkCdn(usm)*us*uw
229 vst = rhoa*exf_BulkCdn(usm)*us*vw
230 endif
231 #endif /*ALLOW_BULK_FORCE*/
232
233 return
234 end
  ViewVC Help
Powered by ViewVC 1.1.22