/[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.1 - (show annotations) (download)
Mon May 14 22:08:40 2001 UTC (23 years, 1 month ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint40pre3, checkpoint40pre1, checkpoint40pre2, checkpoint39
Added external forcing package.
Not presently supported by mitgcm, i.e. disabled by default.

1 c $Header: /u/gcmpack/development/heimbach/ecco_env/pkg/exf/exf_getffields.F,v 1.6 2001/05/08 13:26:33 ralf 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 _RL maskc(1-olx:snx+olx,1-oly:sny+oly)
191
192 c == external functions ==
193
194 integer ilnblnk
195 external ilnblnk
196
197 #ifdef ALLOW_BULKFORMULAE
198 _RL exf_BulkqSat
199 external exf_BulkqSat
200 _RL exf_BulkCdn
201 external exf_BulkCdn
202 _RL exf_BulkRhn
203 external exf_BulkRhn
204 #endif /* ALLOW_BULKFORMULAE */
205
206 c == end of interface ==
207
208 c determine forcing field records
209
210 #ifdef ALLOW_BULKFORMULAE
211
212 c Determine where we are in time and set counters, flags and
213 c the linear interpolation factors accordingly.
214 #ifdef ALLOW_ATM_TEMP
215 c Atmospheric temperature.
216 call exf_set_atemp( atemp
217 & , mycurrenttime, mycurrentiter, mythid )
218
219 c Atmospheric humidity.
220 call exf_set_aqh( aqh
221 & , mycurrenttime, mycurrentiter, mythid )
222
223 c Net long wave radiative flux.
224 call exf_set_lwflux( lwflux
225 & , mycurrenttime, mycurrentiter, mythid )
226
227 c Net short wave radiative flux.
228 call exf_set_swflux( swflux
229 & , mycurrenttime, mycurrentiter, mythid )
230
231 c Precipitation.
232 call exf_set_precip( precip
233 & , mycurrenttime, mycurrentiter, mythid )
234
235 aln = log(ht/zref)
236
237 #else
238
239 c Atmospheric heat flux.
240 call exf_set_hflux( hflux
241 & , mycurrenttime, mycurrentiter, mythid )
242
243 c Salt flux.
244 call exf_set_sflux( sflux
245 & , mycurrenttime, mycurrentiter, mythid )
246
247 #ifdef ALLOW_KPP
248 c Net short wave radiative flux.
249 call exf_set_swflux( swflux
250 & , mycurrenttime, mycurrentiter, mythid )
251
252 #endif /* ALLOW_KPP */
253
254 #endif /* ALLOW_ATM_TEMP */
255
256 #ifdef ALLOW_ATM_WIND
257 c Zonal wind.
258 call exf_set_uwind( uwind
259 & , mycurrenttime, mycurrentiter, mythid )
260
261 c Meridional wind.
262 call exf_set_vwind( vwind
263 & , mycurrenttime, mycurrentiter, mythid )
264
265 #else
266 c Zonal wind stress.
267 call exf_set_ustress( ustress
268 & , mycurrenttime, mycurrentiter, mythid )
269
270 c Meridional wind stress.
271 call exf_set_vstress( vstress
272 & , mycurrenttime, mycurrentiter, mythid )
273
274 #endif /* ALLOW_ATM_WIND */
275
276 #else /* ALLOW_BULKFORMULAE undefined */
277 c Atmospheric heat flux.
278 call exf_set_hflux( hflux
279 & , mycurrenttime, mycurrentiter, mythid )
280
281 c Salt flux.
282 call exf_set_sflux( sflux
283 & , mycurrenttime, mycurrentiter, mythid )
284
285 c Zonal wind stress.
286 call exf_set_ustress( ustress
287 & , mycurrenttime, mycurrentiter, mythid )
288
289 c Meridional wind stress.
290 call exf_set_vstress( vstress
291 & , mycurrenttime, mycurrentiter, mythid )
292
293 #ifdef ALLOW_KPP
294 c Net short wave radiative flux.
295 call exf_set_swflux( swflux
296 & , mycurrenttime, mycurrentiter, mythid )
297 #endif /* ALLOW_KPP */
298
299 #endif /* ALLOW_BULKFORMULAE */
300
301 c Loop over tiles.
302 do bj = mybylo(mythid),mybyhi(mythid)
303 do bi = mybxlo(mythid),mybxhi(mythid)
304 k = 1
305
306 c Calculate mask for tracer cells (0 => land, 1 => water).
307 do j = 1-oly,sny+oly
308 do i = 1-olx,snx+olx
309 maskc(i,j) = exf_one
310 if (_hFacC(i,j,k,bi,bj) .eq. 0.) maskc(i,j) = 0. _d 0
311 enddo
312 enddo
313
314 do j = 1-oly,sny+oly
315 do i = 1-olx,snx+olx
316
317 #ifdef ALLOW_BULKFORMULAE
318
319 c Compute the turbulent surface fluxes.
320 c (Bulk formulae estimates)
321
322 #ifdef ALLOW_ATM_WIND
323 c Wind speed and direction.
324 us = sqrt(uwind(i,j,bi,bj)*uwind(i,j,bi,bj) +
325 & vwind(i,j,bi,bj)*vwind(i,j,bi,bj))
326 if ( us .ne. 0. _d 0 ) then
327 cw = uwind(i,j,bi,bj)/us
328 sw = vwind(i,j,bi,bj)/us
329 else
330 cw = 0. _d 0
331 sw = 0. _d 0
332 endif
333 sh = max(us,umin)
334 #else
335 #ifdef ALLOW_ATM_TEMP
336
337 c The variables us, sh and rdn have to be computed from given
338 c wind stresses inverting relationship for neutral drag coeff.
339 c cdn.
340 c The inversion is based on linear and quadratic form of
341 c cdn(umps); ustar can be directly computed from stress;
342
343 ustar = sqrt(ustress(i,j,bi,bj)*ustress(i,j,bi,bj) +
344 & vstress(i,j,bi,bj)*vstress(i,j,bi,bj))/
345 & atmrho
346 cw = ustress(i,j,bi,bj)/ustar
347 sw = ustress(i,j,bi,bj)/ustar
348
349 if ( ustar .eq. 0. _d 0 ) then
350 us = 0. _d 0
351 else if ( ustar .lt. ustofu11 ) then
352 tmp1 = -cquadrag_2/cquadrag_1/2
353 tmp2 = sqrt(tmp1*tmp1 + ustar*ustar/cquadrag_1)
354 us = sqrt(tmp1 + tmp2)
355 else
356 tmp3 = clindrag_2/clindrag_1/3
357 tmp4 = ustar*ustar/clindrag_1/2 - tmp3**3
358 tmp5 = sqrt(ustar*ustar/clindrag_1*
359 & (ustar*ustar/clindrag_1/4 - tmp3**3))
360 us = (tmp4 + tmp5)**(1/3) +
361 & tmp3**2 * (tmp4 + tmp5)**(-1/3) - tmp3
362 endif
363
364 if ( us .ne. 0 ) then
365 rdn = ustar/us
366 else
367 rdn = 0. _d 0
368 end if
369
370 sh = max(us,umin)
371 #endif /* ALLOW_ATM_TEMP */
372 #endif /* ALLOW_ATM_WIND */
373
374 #ifdef ALLOW_ATM_TEMP
375 c Initial guess: z/l=0.0; hu=ht=hq=z
376 c Iterations: converge on z/l and hence the fluxes.
377 c t0 : virtual temperature (K)
378 c ssq : sea surface humidity (kg/kg)
379 c deltap : potential temperature diff (K)
380
381 if ( atemp(i,j,bi,bj) .ne. 0. _d 0 ) then
382 t0 = atemp(i,j,bi,bj)*
383 & (exf_one + humid_fac*aqh(i,j,bi,bj))
384 ssq = saltsat*
385 & exf_BulkqSat(theta(i,j,1,bi,bj) + cen2kel)/
386 & atmrho
387 deltap = atemp(i,j,bi,bj) + gamma_blk*ht -
388 & theta(i,j,1,bi,bj) - cen2kel
389 delq = aqh(i,j,bi,bj) - ssq
390 stable = exf_half + sign(exf_half, deltap)
391 rdn = sqrt(exf_BulkCdn(sh))
392 ustar = rdn*sh
393 tstar = exf_BulkRhn(stable)*deltap
394 qstar = cdalton*delq
395
396 do iter = 1,niter_bulk
397
398 huol = czol*(tstar/t0 +
399 & qstar/(exf_one/humid_fac+aqh(i,j,bi,bj)))/
400 & ustar**2
401 huol = max(huol,zolmin)
402 stable = exf_half + sign(exf_half, huol)
403 htol = huol*ht/hu
404 hqol = huol*hq/hu
405
406 c Evaluate all stability functions assuming hq = ht.
407 xsq = max(sqrt(abs(exf_one - 16.*huol)),exf_one)
408 x = sqrt(xsq)
409 psimh = -psim_fac*huol*stable +
410 & (exf_one - stable)*
411 & log((exf_one + x*(exf_two + x))*
412 & (exf_one + xsq)/8.) - exf_two*atan(x) +
413 & pi*exf_half
414 xsq = max(sqrt(abs(exf_one - 16.*htol)),exf_one)
415 psixh = -psim_fac*htol*stable + (exf_one - stable)*
416 & exf_two*log((exf_one + xsq)/exf_two)
417
418 c Shift wind speed using old coefficient
419 c rd = rdn/(exf_one + rdn/karman*(log(hu/zref) - psimh) )
420 rd = rdn/(exf_one - rdn/karman*psimh )
421 uzn = max(sh*rd/rdn, umin)
422
423 c Update the transfer coefficients at 10 meters
424 c and neutral stability.
425 rdn = sqrt(exf_BulkCdn(uzn))
426
427 c Shift all coefficients to the measurement height
428 c and stability.
429 c rd = rdn/(exf_one + rdn/karman*(log(hu/zref) - psimh))
430 rd = rdn/(exf_one - rdn/karman*psimh)
431 rh = exf_BulkRhn(stable)/(exf_one +
432 & exf_BulkRhn(stable)/
433 & karman*(aln - psixh))
434 re = cdalton/(exf_one + cdalton/karman*(aln - psixh))
435
436 c Update ustar, tstar, qstar using updated, shifted
437 c coefficients.
438 ustar = rd*sh
439 qstar = re*delq
440 tstar = rh*deltap
441
442 enddo
443
444 tau = atmrho*ustar**2
445 tau = tau*us/sh
446 hs(i,j,bi,bj) = atmcp*tau*tstar/ustar
447 hl(i,j,bi,bj) = flamb*tau*qstar/ustar
448
449 evap(i,j,bi,bj) = tau*qstar/ustar
450
451 ustress(i,j,bi,bj) = tau*cw
452 vstress(i,j,bi,bj) = tau*sw
453 ce hflux(i,j,bi,bj) = atmcp*tau*tstar/ustar +
454 ce & flamb*tau*qstar/ustar
455 else
456 ustress(i,j,bi,bj) = 0. _d 0
457 vstress(i,j,bi,bj) = 0. _d 0
458 hflux (i,j,bi,bj) = 0. _d 0
459 hs(i,j,bi,bj) = 0. _d 0
460 hl(i,j,bi,bj) = 0. _d 0
461 endif
462
463 #else
464 #ifdef ALLOW_ATM_WIND
465 ustress(i,j,bi,bj) = atmrho*exf_BulkCdn(sh)*us*
466 & uwind(i,j,bi,bj)
467 vstress(i,j,bi,bj) = atmrho*exf_BulkCdn(sh)*us*
468 & vwind(i,j,bi,bj)
469 #endif /* ALLOW_ATM_WIND */
470 #endif /* ALLOW_ATM_TEMP */
471 enddo
472 enddo
473 enddo
474 enddo
475
476 c Add all contributions.
477 do bj = mybylo(mythid),mybyhi(mythid)
478 do bi = mybxlo(mythid),mybxhi(mythid)
479 do j = 1,sny
480 do i = 1,snx
481 c Net surface heat flux.
482 #ifdef ALLOW_ATM_TEMP
483 hfl = 0. _d 0
484 hfl = hfl - hs(i,j,bi,bj)
485 hfl = hfl - hl(i,j,bi,bj)
486 hfl = hfl + lwflux(i,j,bi,bj)
487 #ifndef ALLOW_KPP
488 hfl = hfl + swflux(i,j,bi,bj)
489 #endif /* ALLOW_KPP undef */
490 c Heat flux:
491 hflux(i,j,bi,bj) = hfl*maskc(i,j)
492 c Salt flux from Precipitation and Evaporation.
493 sflux(i,j,bi,bj) = precip(i,j,bi,bj) - evap(i,j,bi,bj)
494 #endif /* ALLOW_ATM_TEMP */
495
496 #else
497 hflux(i,j,bi,bj) = hflux(i,j,bi,bj)*maskc(i,j)
498 sflux(i,j,bi,bj) = sflux(i,j,bi,bj)*maskc(i,j)
499 #endif /* ALLOW_BULKFORMULAE */
500 enddo
501 enddo
502 enddo
503 enddo
504
505 c Update the tile edges.
506 _EXCH_XY_R8(hflux, mythid)
507 _EXCH_XY_R8(sflux, mythid)
508 _EXCH_XY_R8(ustress, mythid)
509 _EXCH_XY_R8(vstress, mythid)
510
511 #ifdef ALLOW_KPP
512 _EXCH_XY_R8(swflux, mythid)
513 #endif /* ALLOW_KPP */
514
515 end

  ViewVC Help
Powered by ViewVC 1.1.22