/[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.4 - (show annotations) (download)
Tue Nov 12 20:34:41 2002 UTC (21 years, 7 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint46n_post, checkpoint47c_post, checkpoint47d_pre, checkpoint47a_post, checkpoint47d_post, checkpoint47b_post, checkpoint47
Changes since 1.3: +14 -11 lines
Merging from release1_p8:
o exf:
  updated external forcing package
  - bug fixes carried over from ecco-branch
    (missing OBCS_OPTIONS.h in two routines)
  - enable easy to use "no forcing".
  - added exf I/O for atmospheric loading
  - added exf I/O for runoff data
  - transfered scaling between exf <-> MITgcm to exf namelist
  - removing old exfa stuff

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

  ViewVC Help
Powered by ViewVC 1.1.22