da_mo_correction.inc
References to this file elsewhere.
1 function da_mo_correction (ho, po, to, qo, &
2 hm, pm, tm, qm, um ,vm, &
3 roughness) RESULT (correc)
4
5 !----------------------------------------------------------------------------
6 ! Purpose: Compute the correction factor to convert surface wind into 40m
7 ! wind using similarity laws.
8 !
9 ! Reference:
10 ! ---------
11 ! A description of the fifth generation Penn State/NCAR Mesoscale Model
12 ! Grell et al. 1994, page 29-30 and 80-83.
13 !----------------------------------------------------------------------------
14
15 implicit none
16
17 real, intent (in) :: ho, po, to, qo
18 real, intent (in) :: hm, pm, tm, qm, um, vm
19 real, intent (in) :: roughness
20
21 real :: correc, winfac
22
23 real :: thm , tho, tvm, tvo, thvm, thvo, rcp
24 real :: za, Vc2, Va2, V2
25 ! FIX? real :: hdif, rib, rll, hll, zint
26 real :: hdif, rib, hll, zint
27
28 ! Height difference (in m) above wich correction is applied
29
30 real, parameter :: hmax = 10.0, hs_max = 40.0
31
32 ! Default Roughness length im m (for land, set to 0 if on sea)
33
34 real, parameter :: zint0 = 0.2
35
36 rcp = gas_constant/cp
37
38 ! 0.0 initialize correction factor to 1
39 ! =====================================
40
41 correc = 1.0
42 winfac = 1.0
43
44 ! 1. height difference and roughness length
45 ! ==========================================
46
47 ! 1.1 Height difference
48 ! -----------------
49
50 hdif = hm - ho
51
52 ! 1.2 No correction if height difference is less than hmax.
53 ! -----------------------------------------------------
54
55 if (hdif <= hmax) then
56 return
57 end if
58
59 ! too many
60 ! if (trace_use) call da_trace_entry("da_mo_correction")
61
62 ! 1.3 Compute the roughness length based upon season and land use
63 ! -----------------------------------------------------------
64
65 zint = roughness
66
67 if (zint < 0.0001) zint = 0.0001
68
69 ! 2. potential temperature at model surface and observation location
70 ! ===================================================================
71
72 ! 2.1 potential temperature on model surface
73 ! --------------------------------------
74
75 thm = tm * (1000.0 / (pm/100.0)) ** rcp
76
77 ! 2.2 potential temperature at observation location
78 ! ---------------------------------------------
79
80 tho = to * (1000.0 / (po/100.0)) ** rcp
81
82 ! 3. virtual temperature at model surface and observation location
83 ! ===================================================================
84
85 ! 3.1 Virtual temperature on model surface
86 ! -------------------------------------
87
88 tvm = tm * (1.0 + 0.608 * qm)
89
90 ! 3.2 Virtual temperature at observation location
91 ! -------------------------------------------
92
93 tvo = to * (1.0 + 0.608 * qo)
94
95 ! 4. potential virtual temperature at model surface and observation location
96 ! ===========================================================================
97
98 ! 4.1 potential virtual temperature on model surface
99 ! ----------------------------------------------
100
101 thvm = tvm * (1000.0 / (pm/100.0)) ** rcp
102
103 ! 4.2 potential virtual temperature at observation location
104 ! -----------------------------------------------------
105
106 thvo = tvo * (1000.0 / (po/100.0)) ** rcp
107
108
109 ! 5. bulk richardson number and moni-obukov length
110 ! =================================================
111
112 ! 5.1 Pre-calculations
113 ! ----------------
114
115 za = hm - ho
116
117 ! Because the following formula is derived based on
118 ! the surface layer height is hs_max=40m. Under
119 ! free convection, we assume that atmospheric state
120 ! above the surface layer is fully mixed, and the
121 ! wind at the lowest sigma half level is same as the
122 ! wind at top of the surface layer.
123
124 if (za > hs_max) za = hs_max
125
126 Va2 = um*um + vm*vm
127 Vc2 = 4.0 * MAX ((tho - thm), 0.0)
128
129 V2 = Va2 + Vc2
130
131 ! 5.2 Bulk richardson number
132 ! ----------------------
133
134 rib = (gravity * za / thm) * (thvm - thvo) / V2
135
136 ! 5.3 Monin-obukov length
137 ! -------------------
138
139 hll = rib * LOG (za/zint)
140
141 ! FIX? is this right?
142 ! rll = za / hll
143
144 ! 5.4 Ratio PBL height/Monin Obukov length: za/rll
145 ! ------------------------------------
146
147 hll = ABS (hll)
148
149
150 ! 6. CORRECTION FACTOR BASED UPON REGIME
151 ! =======================================
152
153 ! 6.1 Stable conditions (REGIME 1)
154 ! ---------------------------
155
156 ! correc = 1.0 ! rib > 0.2
157
158 ! 6.2 Mechanically driven turbulence (REGIME 2)
159 ! ------------------------------------------
160
161 ! correc = 1.0 ! 0.0 =< rib <= 0.2
162
163 correc = 1.0
164
165 if (rib < 0.0) then
166
167 ! 6.3 Unstable Forced convection (REGIME 3)
168 ! -------------------------------------
169
170 ! correc = 1.0 ! hll <= 1.5
171
172
173 ! 6.4 Free convection (REGIME 4)
174 ! --------------------------
175
176 if (hll > 1.5) then
177 if (zint < 0.2) then
178 correc = 1.000 + 0.320 * zint ** 0.200
179 else if (zint < 0.0) then
180 correc = 1.169 + 0.315 * zint
181 end if
182
183 ! 6.4.1 The factor depended on Za (MM5, param.F)
184
185 winfac = 0.5*(log(za/0.05)/log(40.0/0.05) &
186 + log(za)/log(40.0))
187
188 correc = correc * winfac
189 end if
190 end if
191
192 ! if (trace_use) call da_trace_exit("da_mo_correction")
193
194 end function da_mo_correction
195
196