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

Annotation 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 - (hide annotations) (download)
Thu Oct 9 04:19:19 2003 UTC (20 years, 7 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 edhill 1.3 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 cheisey 1.1 c swd -- bulkf formula used in bulkf and ice pkgs
4    
5     c taken from exf package
6 edhill 1.3 #include "BULK_FORCE_OPTIONS.h"
7 cheisey 1.1 subroutine bulkf_formula_lanl(
8     I uw, vw, us, Ta, Qa, nc, tsf_in,
9 cheisey 1.2 I flwupa, flha, fsha, df0dT,
10 cheisey 1.1 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 cheisey 1.2 _RL flwupa ! upward long wave radiation
41 cheisey 1.1 _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 cheisey 1.2 cQQ flwupa=-0.985*stefan*tsf**4
185 cheisey 1.1 cQQ & *(0.39-0.05*sqrt(ea))
186     cQQ & *(1-0.6*nc**2)
187     if (iceornot.eq.0) then
188 cheisey 1.2 flwupa=-ocean_emissivity*stefan*tsf**4
189 cheisey 1.1 else
190     if (iceornot.eq.2) then
191 cheisey 1.2 flwupa=-snow_emissivity*stefan*tsf**4
192 cheisey 1.1 else
193 cheisey 1.2 flwupa=-ice_emissivity*stefan*tsf**4
194 cheisey 1.1 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