/[MITgcm]/MITgcm/pkg/exf/exf_getffields.F
ViewVC logotype

Contents of /MITgcm/pkg/exf/exf_getffields.F

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


Revision 1.2 - (show annotations) (download)
Mon Jul 30 20:41:20 2001 UTC (22 years, 10 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint43a-release1mods, checkpoint40pre7, checkpoint40pre6, checkpoint40pre9, checkpoint40pre8, release1_b1, checkpoint43, release1-branch_tutorials, checkpoint40pre4, release1-branch-end, checkpoint40pre5, ecco-branch-mod1, release1_beta1, checkpoint42, checkpoint40, checkpoint41, release1-branch_branchpoint
Branch point for: release1-branch, release1, ecco-branch, release1_coupled
Changes since 1.1: +5 -14 lines
Updated to c40 code.

1 c $Header: /u/gcmpack/models/MITgcmUV/pkg/exf/exf_getffields.F,v 1.1 2001/05/14 22:08:40 heimbach Exp $
2
3 #include "EXF_CPPOPTIONS.h"
4
5 subroutine exf_GetFFields(
6 I mycurrenttime,
7 I mycurrentiter,
8 I mythid
9 & )
10
11 c ==================================================================
12 c SUBROUTINE exf_GetFFields
13 c ==================================================================
14 c
15 c o Get the surface fluxes either from file or as derived from bulk
16 c formulae that use the atmospheric state.
17 c
18 c The calculation of the bulk surface fluxes has been adapted from
19 c the NCOM model which uses the formulae given in Large and Pond
20 c (1981 & 1982 )
21 c
22 c
23 c Header taken from NCOM version: ncom1.4.1
24 c -----------------------------------------
25 c
26 c Following procedures and coefficients in Large and Pond
27 c (1981 ; 1982)
28 c
29 c Output: Bulk estimates of the turbulent surface fluxes.
30 c -------
31 c
32 c hs - sensible heat flux (W/m^2), into ocean
33 c hl - latent heat flux (W/m^2), into ocean
34 c
35 c Input:
36 c ------
37 c
38 c us - mean wind speed (m/s) at height hu (m)
39 c th - mean air temperature (K) at height ht (m)
40 c qh - mean air humidity (kg/kg) at height hq (m)
41 c sst - sea surface temperature (K)
42 c tk0 - Kelvin temperature at 0 Celsius (K)
43 c
44 c Assume 1) a neutral 10m drag coefficient =
45 c
46 c cdn = .0027/u10 + .000142 + .0000764 u10
47 c
48 c 2) a neutral 10m stanton number =
49 c
50 c ctn = .0327 sqrt(cdn), unstable
51 c ctn = .0180 sqrt(cdn), stable
52 c
53 c 3) a neutral 10m dalton number =
54 c
55 c cen = .0346 sqrt(cdn)
56 c
57 c 4) the saturation humidity of air at
58 c
59 c t(k) = exf_BulkqSat(t) (kg/m^3)
60 c
61 c Note: 1) here, tstar = <wt>/u*, and qstar = <wq>/u*.
62 c 2) wind speeds should all be above a minimum speed,
63 c say 0.5 m/s.
64 c 3) with optional iteration loop, niter=3, should suffice.
65 c 4) this version is for analyses inputs with hu = 10m and
66 c ht = hq.
67 c 5) sst enters in Celsius.
68 c
69 c ------------------------------------
70 c
71 c By setting CPP options in the header file *EXF_CPPOPTIONS.h* it
72 c is possible to combine data sets in four different ways:
73 c
74 c The following options are available:
75 c
76 c ALLOW_ATM_TEMP (UAT)
77 c ALLOW_ATM_WIND (UAW)
78 c
79 c
80 c UAT | UAW | action
81 c ----------------------------------------------------
82 c undefined | undefined | Use surface fluxes.
83 c undefined | defined | Assume cdn(u) given to
84 c | | infer the wind stress.
85 c defined | undefined | Compute wind field from
86 c | | given stress assuming a
87 c | | linear relation.
88 c defined | defined | Use the bulk formulae.
89 c ----------------------------------------------------
90 c
91 c Implementations of the bulk formulae exist for the follwing
92 c versions of the MITgcm:
93 c
94 c MITgcm : Patrick Heimbach
95 c MITgcmUV: Christian Eckert
96 c
97 c started: Christian Eckert eckert@mit.edu 27-Aug-1999
98 c
99 c changed: Christian Eckert eckert@mit.edu 14-Jan-2000
100 c
101 c - restructured the original version in order to have a
102 c better interface to the MITgcmUV.
103 c
104 c Christian Eckert eckert@mit.edu 12-Feb-2000
105 c
106 c - Changed Routine names (package prefix: exf_)
107 c
108 c Patrick Heimbach, heimbach@mit.edu 04-May-2000
109 c
110 c - changed the handling of precip and sflux with respect
111 c to CPP options ALLOW_BULKFORMULAE and ALLOW_ATM_TEMP
112 c
113 c - included some CPP flags ALLOW_BULKFORMULAE to make
114 c sure ALLOW_ATM_TEMP, ALLOW_ATM_WIND are used only in
115 c conjunction with defined ALLOW_BULKFORMULAE
116 c
117 c - statement functions discarded
118 c
119 c Ralf.Giering@FastOpt.de 25-Mai-2000
120 c
121 c - total rewrite using new subroutines
122 c
123 c ==================================================================
124 c SUBROUTINE exf_GetFFields
125 c ==================================================================
126
127 implicit none
128
129 c == global variables ==
130
131 #include "EEPARAMS.h"
132 #include "SIZE.h"
133 #include "PARAMS.h"
134 #include "DYNVARS.h"
135 #include "GRID.h"
136
137 #include "exf_fields.h"
138 #include "exf_constants.h"
139
140 c == routine arguments ==
141
142 integer mythid
143 integer mycurrentiter
144 _RL mycurrenttime
145
146 c == local variables ==
147
148 integer bi,bj
149 integer i,j,k
150
151 #ifdef ALLOW_BULKFORMULAE
152
153 #ifdef ALLOW_ATM_TEMP
154 integer iter
155 _RL aln
156 _RL delq
157 _RL deltap
158 _RL hqol
159 _RL htol
160 _RL huol
161 _RL psimh
162 _RL psixh
163 _RL qstar
164 _RL rd
165 _RL re
166 _RL rdn
167 _RL rh
168 _RL ssq
169 _RL stable
170 _RL tau
171 _RL tstar
172 _RL t0
173 _RL ustar
174 _RL uzn
175 _RL xsq
176 _RL x
177 _RL evap(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
178 #endif /* ALLOW_ATM_TEMP */
179
180 _RL us
181 _RL cw
182 _RL sw
183 _RL sh
184 _RL hs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
185 _RL hl(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
186 _RL hfl
187
188 #endif /* ALLOW_BULKFORMULAE */
189
190 c == external functions ==
191
192 integer ilnblnk
193 external ilnblnk
194
195 #ifdef ALLOW_BULKFORMULAE
196 _RL exf_BulkqSat
197 external exf_BulkqSat
198 _RL exf_BulkCdn
199 external exf_BulkCdn
200 _RL exf_BulkRhn
201 external exf_BulkRhn
202 #endif /* ALLOW_BULKFORMULAE */
203
204 c == end of interface ==
205
206 c determine forcing field records
207
208 #ifdef ALLOW_BULKFORMULAE
209
210 c Determine where we are in time and set counters, flags and
211 c the linear interpolation factors accordingly.
212 #ifdef ALLOW_ATM_TEMP
213 c Atmospheric temperature.
214 call exf_set_atemp( atemp
215 & , mycurrenttime, mycurrentiter, mythid )
216
217 c Atmospheric humidity.
218 call exf_set_aqh( aqh
219 & , mycurrenttime, mycurrentiter, mythid )
220
221 c Net long wave radiative flux.
222 call exf_set_lwflux( lwflux
223 & , mycurrenttime, mycurrentiter, mythid )
224
225 c Net short wave radiative flux.
226 call exf_set_swflux( swflux
227 & , mycurrenttime, mycurrentiter, mythid )
228
229 c Precipitation.
230 call exf_set_precip( precip
231 & , mycurrenttime, mycurrentiter, mythid )
232
233 aln = log(ht/zref)
234
235 #else
236
237 c Atmospheric heat flux.
238 call exf_set_hflux( hflux
239 & , mycurrenttime, mycurrentiter, mythid )
240
241 c Salt flux.
242 call exf_set_sflux( sflux
243 & , mycurrenttime, mycurrentiter, mythid )
244
245 #ifdef ALLOW_KPP
246 c Net short wave radiative flux.
247 call exf_set_swflux( swflux
248 & , mycurrenttime, mycurrentiter, mythid )
249
250 #endif /* ALLOW_KPP */
251
252 #endif /* ALLOW_ATM_TEMP */
253
254 #ifdef ALLOW_ATM_WIND
255 c Zonal wind.
256 call exf_set_uwind( uwind
257 & , mycurrenttime, mycurrentiter, mythid )
258
259 c Meridional wind.
260 call exf_set_vwind( vwind
261 & , mycurrenttime, mycurrentiter, mythid )
262
263 #else
264 c Zonal wind stress.
265 call exf_set_ustress( ustress
266 & , mycurrenttime, mycurrentiter, mythid )
267
268 c Meridional wind stress.
269 call exf_set_vstress( vstress
270 & , mycurrenttime, mycurrentiter, mythid )
271
272 #endif /* ALLOW_ATM_WIND */
273
274 #else /* ALLOW_BULKFORMULAE undefined */
275 c Atmospheric heat flux.
276 call exf_set_hflux( hflux
277 & , mycurrenttime, mycurrentiter, mythid )
278
279 c Salt flux.
280 call exf_set_sflux( sflux
281 & , mycurrenttime, mycurrentiter, mythid )
282
283 c Zonal wind stress.
284 call exf_set_ustress( ustress
285 & , mycurrenttime, mycurrentiter, mythid )
286
287 c Meridional wind stress.
288 call exf_set_vstress( vstress
289 & , mycurrenttime, mycurrentiter, mythid )
290
291 #ifdef ALLOW_KPP
292 c Net short wave radiative flux.
293 call exf_set_swflux( swflux
294 & , mycurrenttime, mycurrentiter, mythid )
295 #endif /* ALLOW_KPP */
296
297 #endif /* ALLOW_BULKFORMULAE */
298
299 c Loop over tiles.
300 do bj = mybylo(mythid),mybyhi(mythid)
301 do bi = mybxlo(mythid),mybxhi(mythid)
302 k = 1
303
304 do j = 1-oly,sny+oly
305 do i = 1-olx,snx+olx
306
307 #ifdef ALLOW_BULKFORMULAE
308
309 c Compute the turbulent surface fluxes.
310 c (Bulk formulae estimates)
311
312 #ifdef ALLOW_ATM_WIND
313 c Wind speed and direction.
314 us = sqrt(uwind(i,j,bi,bj)*uwind(i,j,bi,bj) +
315 & vwind(i,j,bi,bj)*vwind(i,j,bi,bj))
316 if ( us .ne. 0. _d 0 ) then
317 cw = uwind(i,j,bi,bj)/us
318 sw = vwind(i,j,bi,bj)/us
319 else
320 cw = 0. _d 0
321 sw = 0. _d 0
322 endif
323 sh = max(us,umin)
324 #else
325 #ifdef ALLOW_ATM_TEMP
326
327 c The variables us, sh and rdn have to be computed from given
328 c wind stresses inverting relationship for neutral drag coeff.
329 c cdn.
330 c The inversion is based on linear and quadratic form of
331 c cdn(umps); ustar can be directly computed from stress;
332
333 ustar = sqrt(ustress(i,j,bi,bj)*ustress(i,j,bi,bj) +
334 & vstress(i,j,bi,bj)*vstress(i,j,bi,bj))/
335 & atmrho
336 cw = ustress(i,j,bi,bj)/ustar
337 sw = ustress(i,j,bi,bj)/ustar
338
339 if ( ustar .eq. 0. _d 0 ) then
340 us = 0. _d 0
341 else if ( ustar .lt. ustofu11 ) then
342 tmp1 = -cquadrag_2/cquadrag_1/2
343 tmp2 = sqrt(tmp1*tmp1 + ustar*ustar/cquadrag_1)
344 us = sqrt(tmp1 + tmp2)
345 else
346 tmp3 = clindrag_2/clindrag_1/3
347 tmp4 = ustar*ustar/clindrag_1/2 - tmp3**3
348 tmp5 = sqrt(ustar*ustar/clindrag_1*
349 & (ustar*ustar/clindrag_1/4 - tmp3**3))
350 us = (tmp4 + tmp5)**(1/3) +
351 & tmp3**2 * (tmp4 + tmp5)**(-1/3) - tmp3
352 endif
353
354 if ( us .ne. 0 ) then
355 rdn = ustar/us
356 else
357 rdn = 0. _d 0
358 end if
359
360 sh = max(us,umin)
361 #endif /* ALLOW_ATM_TEMP */
362 #endif /* ALLOW_ATM_WIND */
363
364 #ifdef ALLOW_ATM_TEMP
365 c Initial guess: z/l=0.0; hu=ht=hq=z
366 c Iterations: converge on z/l and hence the fluxes.
367 c t0 : virtual temperature (K)
368 c ssq : sea surface humidity (kg/kg)
369 c deltap : potential temperature diff (K)
370
371 if ( atemp(i,j,bi,bj) .ne. 0. _d 0 ) then
372 t0 = atemp(i,j,bi,bj)*
373 & (exf_one + humid_fac*aqh(i,j,bi,bj))
374 ssq = saltsat*
375 & exf_BulkqSat(theta(i,j,1,bi,bj) + cen2kel)/
376 & atmrho
377 deltap = atemp(i,j,bi,bj) + gamma_blk*ht -
378 & theta(i,j,1,bi,bj) - cen2kel
379 delq = aqh(i,j,bi,bj) - ssq
380 stable = exf_half + sign(exf_half, deltap)
381 rdn = sqrt(exf_BulkCdn(sh))
382 ustar = rdn*sh
383 tstar = exf_BulkRhn(stable)*deltap
384 qstar = cdalton*delq
385
386 do iter = 1,niter_bulk
387
388 huol = czol*(tstar/t0 +
389 & qstar/(exf_one/humid_fac+aqh(i,j,bi,bj)))/
390 & ustar**2
391 huol = max(huol,zolmin)
392 stable = exf_half + sign(exf_half, huol)
393 htol = huol*ht/hu
394 hqol = huol*hq/hu
395
396 c Evaluate all stability functions assuming hq = ht.
397 xsq = max(sqrt(abs(exf_one - 16.*huol)),exf_one)
398 x = sqrt(xsq)
399 psimh = -psim_fac*huol*stable +
400 & (exf_one - stable)*
401 & log((exf_one + x*(exf_two + x))*
402 & (exf_one + xsq)/8.) - exf_two*atan(x) +
403 & pi*exf_half
404 xsq = max(sqrt(abs(exf_one - 16.*htol)),exf_one)
405 psixh = -psim_fac*htol*stable + (exf_one - stable)*
406 & exf_two*log((exf_one + xsq)/exf_two)
407
408 c Shift wind speed using old coefficient
409 c rd = rdn/(exf_one + rdn/karman*(log(hu/zref) - psimh) )
410 rd = rdn/(exf_one - rdn/karman*psimh )
411 uzn = max(sh*rd/rdn, umin)
412
413 c Update the transfer coefficients at 10 meters
414 c and neutral stability.
415 rdn = sqrt(exf_BulkCdn(uzn))
416
417 c Shift all coefficients to the measurement height
418 c and stability.
419 c rd = rdn/(exf_one + rdn/karman*(log(hu/zref) - psimh))
420 rd = rdn/(exf_one - rdn/karman*psimh)
421 rh = exf_BulkRhn(stable)/(exf_one +
422 & exf_BulkRhn(stable)/
423 & karman*(aln - psixh))
424 re = cdalton/(exf_one + cdalton/karman*(aln - psixh))
425
426 c Update ustar, tstar, qstar using updated, shifted
427 c coefficients.
428 ustar = rd*sh
429 qstar = re*delq
430 tstar = rh*deltap
431
432 enddo
433
434 tau = atmrho*ustar**2
435 tau = tau*us/sh
436 hs(i,j,bi,bj) = atmcp*tau*tstar/ustar
437 hl(i,j,bi,bj) = flamb*tau*qstar/ustar
438
439 evap(i,j,bi,bj) = tau*qstar/ustar
440
441 ustress(i,j,bi,bj) = tau*cw
442 vstress(i,j,bi,bj) = tau*sw
443 ce hflux(i,j,bi,bj) = atmcp*tau*tstar/ustar +
444 ce & flamb*tau*qstar/ustar
445 else
446 ustress(i,j,bi,bj) = 0. _d 0
447 vstress(i,j,bi,bj) = 0. _d 0
448 hflux (i,j,bi,bj) = 0. _d 0
449 hs(i,j,bi,bj) = 0. _d 0
450 hl(i,j,bi,bj) = 0. _d 0
451 endif
452
453 #else
454 #ifdef ALLOW_ATM_WIND
455 ustress(i,j,bi,bj) = atmrho*exf_BulkCdn(sh)*us*
456 & uwind(i,j,bi,bj)
457 vstress(i,j,bi,bj) = atmrho*exf_BulkCdn(sh)*us*
458 & vwind(i,j,bi,bj)
459 #endif /* ALLOW_ATM_WIND */
460 #endif /* ALLOW_ATM_TEMP */
461 enddo
462 enddo
463 enddo
464 enddo
465
466 c Add all contributions.
467 k = 1
468 do bj = mybylo(mythid),mybyhi(mythid)
469 do bi = mybxlo(mythid),mybxhi(mythid)
470 do j = 1,sny
471 do i = 1,snx
472 c Net surface heat flux.
473 #ifdef ALLOW_ATM_TEMP
474 hfl = 0. _d 0
475 hfl = hfl - hs(i,j,bi,bj)
476 hfl = hfl - hl(i,j,bi,bj)
477 hfl = hfl + lwflux(i,j,bi,bj)
478 #ifndef ALLOW_KPP
479 hfl = hfl + swflux(i,j,bi,bj)
480 #endif /* ALLOW_KPP undef */
481 c Heat flux:
482 hflux(i,j,bi,bj) = hfl*maskc(i,j,k,bi,bj)
483 c Salt flux from Precipitation and Evaporation.
484 sflux(i,j,bi,bj) = precip(i,j,bi,bj) - evap(i,j,bi,bj)
485 #endif /* ALLOW_ATM_TEMP */
486
487 #else
488 hflux(i,j,bi,bj) = hflux(i,j,bi,bj)*maskc(i,j,k,bi,bj)
489 sflux(i,j,bi,bj) = sflux(i,j,bi,bj)*maskc(i,j,k,bi,bj)
490 #endif /* ALLOW_BULKFORMULAE */
491 enddo
492 enddo
493 enddo
494 enddo
495
496 c Update the tile edges.
497 _EXCH_XY_R8(hflux, mythid)
498 _EXCH_XY_R8(sflux, mythid)
499 _EXCH_XY_R8(ustress, mythid)
500 _EXCH_XY_R8(vstress, mythid)
501
502 #ifdef ALLOW_KPP
503 _EXCH_XY_R8(swflux, mythid)
504 #endif /* ALLOW_KPP */
505
506 end

  ViewVC Help
Powered by ViewVC 1.1.22