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

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

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

Revision 1.2 - (show annotations) (download)
Thu Oct 9 04:19:19 2003 UTC (20 years, 8 months ago) by edhill
Branch: MAIN
CVS Tags: checkpoint51l_post, checkpoint51k_post, checkpoint51o_post, checkpoint51q_post, checkpoint52a_pre, checkpoint51r_post, checkpoint52, checkpoint52b_pre, checkpoint51o_pre, checkpoint51t_post, checkpoint52a_post, checkpoint51i_post, checkpoint51n_pre, checkpoint51p_post, checkpoint51n_post, checkpoint51l_pre, checkpoint51u_post, checkpoint51m_post, ecco_c52_e35, checkpoint51s_post
Branch point for: checkpoint51n_branch, branch-nonh, tg2-branch
Changes since 1.1: +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.F,v 1.1.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 #include "BULK_FORCE_OPTIONS.h"
6 subroutine bulkf_formula(
7 I uw, vw, us, Ta, Qa, nc, tsf_in,
8 I flwup, flha, fsha, df0dT,
9 I ust, vst, evp, ssq, iceornot, windread
10 & )
11
12 IMPLICIT NONE
13 c
14 #include "EEPARAMS.h"
15 #include "SIZE.h"
16 #include "PARAMS.h"
17 #include "DYNVARS.h"
18 #include "GRID.h"
19 #include "FFIELDS.h"
20 #ifdef ALLOW_BULKFORMULA
21 #include "BULKF_ICE_CONSTANTS.h"
22 #include "BULKF.h"
23 #endif
24
25 c
26 c calculate bulk formula fluxes over open ocean or seaice
27 c
28 c input
29 _RL us ! wind speed
30 _RL uw ! zonal wind speed (at grid center)
31 _RL vw ! meridional wind speed (at grid center)
32 _RL Ta ! air temperature at ht
33 _RL Qa ! specific humidity at ht
34 _RL tsf_in ! surface temperature (either ice or open water)
35 _RL nc ! fraction cloud cover
36 integer iceornot ! 0=open water, 1=ice cover
37 logical windread !
38 c output
39 _RL flwup ! upward long wave radiation
40 _RL flha ! latent heat flux
41 _RL fsha ! sensible heat flux
42 _RL df0dT ! derivative of heat flux with respect to Tsf
43 _RL ust ! zonal wind stress (at grid center)
44 _RL vst ! meridional wind stress (at grid center)
45 _RL evp ! evaporation rate (over open water)
46 _RL ssq ! surface specific humidity (kg/kg)
47 c
48 c local variables
49 _RL tsf ! surface temperature in K
50 _RL ht ! reference height (m)
51 _RL usm ! wind speed limited
52 _RL t0 ! virtual temperature (K)
53 _RL deltap ! potential temperature diff (K)
54 _RL delq ! specific humidity difference
55 _RL stable
56 _RL lhcoef
57 _RL shcoef
58 _RL dflhdT ! derivative of latent heat with respect to T
59 _RL dfshdT ! derivative of sensible heat with respect to T
60 _RL dflwupdT ! derivative of long wave with respect to T
61 _RL mixratio
62 _RL ea
63
64 #ifdef ALLOW_BULKFORMULA
65
66 c --- external functions ---
67 _RL exf_BulkCdn
68 external exf_BulkCdn
69
70 Qa=max(Qa,1e-7)
71 Tsf=Tsf_in+Tf0kel
72 ht=10
73 c Wind speed
74 if (us.eq.0) then
75 us = sqrt(uw*uw + vw*vw)
76 endif
77 usm = max(us,1.d-5)
78 cQQQ try to limit drag
79 usm = min(usm,5.d0)
80 c
81 t0 = Ta*(1.d0 + humid_fac*Qa)
82 ssq = qcoef*exp(22.47*(1.d0-Tf0kel/tsf))
83 cBB debugging
84 cBB print*,'ice, ssq', qcoef, ssq
85 c
86 deltap = ta - tsf + gamma_blk*ht
87 delq = Qa - ssq
88 stable = 5.d-1 + sign(5.d-1, deltap)
89 cQQ -- check all of this
90 c latent heat coeff. following Gill 2.4.6 and 2.6.2 formula
91 lhcoef =Lvap*rhoa*usm*1.5e-3 !QQ name coeff.
92 flha = lhcoef*delq
93 evp= lhcoef*delq/Lvap/rhosw
94 c sensible heat coeff. following Gill 2.4.5 formula
95 if (stable.gt.0) then
96 shcoef = 0.83e-3*rhoa*usm*cpair
97 else
98 shcoef = 1.10e-3*rhoa*usm*cpair
99 endif
100 fsha = shcoef*deltap
101 c up long wave radiation
102 mixratio=Qa/(1-Qa)
103 ea=p0*mixratio/(0.62197+mixratio)
104 flwup=-0.985*stefan*tsf**4
105 & *(0.39-0.05*sqrt(ea))
106 & *(1-0.6*nc**2)
107 #ifdef ALLOW_TSEAICE
108 if (iceornot.eq.1) then
109 c derivatives with respect to Tsf
110 dflhdT = -lhcoef*Tf0kel*ssq*22.47/(tsf**2)
111 dfshdT = -shcoef
112 dflwupdT=-4.*0.985*stefan*tsf**3
113 & *(0.39-0.05*sqrt(ea))
114 c total derivative with respect to surface temperature
115 df0dT=dflwupdT+dfshdT+dflhdT
116 cBB
117 cBB print*,'derivatives:',df0dT,dflwupdT,dfshdT,dflhdT
118 endif
119 #endif /*ALLOW_TSEAICE*/
120 c
121 c wind stress at center points
122 if (.NOT.windread) then
123 ust = rhoa*exf_BulkCdn(usm)*us*uw
124 vst = rhoa*exf_BulkCdn(usm)*us*vw
125 endif
126 #endif /*ALLOW_BULKFORMULA*/
127
128 return
129 end
  ViewVC Help
Powered by ViewVC 1.1.22