pkg/exf/exf_getffields.F

c $Header: /u/gcmpack/development/heimbach/ecco_env/pkg/exf/exf_getffields.F,v 1.13 2002/01/11 17:34:46 heimbach Exp $

#include "EXF_CPPOPTIONS.h"

      subroutine exf_GetFFields( mycurrenttime, mycurrentiter, mythid )

c     ==================================================================
c     SUBROUTINE exf_GetFFields
c     ==================================================================
c
c     o Get the surface fluxes either from file or as derived from bulk
c       formulae that use the atmospheric state.
c
c       The calculation of the bulk surface fluxes has been adapted from
c       the NCOM model which uses the formulae given in Large and Pond
c       (1981 & 1982 )
c
c
c       Header taken from NCOM version: ncom1.4.1
c       -----------------------------------------
c
c       Following procedures and coefficients in Large and Pond
c       (1981 ; 1982)
c
c       Output: Bulk estimates of the turbulent surface fluxes.
c       -------
c
c       hs  - sensible heat flux  (W/m^2), into ocean
c       hl  - latent   heat flux  (W/m^2), into ocean
c
c       Input:
c       ------
c
c       us  - mean wind speed (m/s)     at height hu (m)
c       th  - mean air temperature (K)  at height ht (m)
c       qh  - mean air humidity (kg/kg) at height hq (m)
c       sst - sea surface temperature (K)
c       tk0 - Kelvin temperature at 0 Celsius (K)
c
c       Assume 1) a neutral 10m drag coefficient =
c
c                 cdn = .0027/u10 + .000142 + .0000764 u10
c
c              2) a neutral 10m stanton number =
c
c                 ctn = .0327 sqrt(cdn), unstable
c                 ctn = .0180 sqrt(cdn), stable
c
c              3) a neutral 10m dalton number =
c
c                 cen = .0346 sqrt(cdn)
c
c              4) the saturation humidity of air at
c
c                 t(k) = exf_BulkqSat(t)  (kg/m^3)
c
c       Note:  1) here, tstar = <wt>/u*, and qstar = <wq>/u*.
c              2) wind speeds should all be above a minimum speed,
c                 say 0.5 m/s.
c              3) with optional iteration loop, niter=3, should suffice.
c              4) this version is for analyses inputs with hu = 10m and
c                 ht = hq.
c              5) sst enters in Celsius.
c
c       ------------------------------------
c
c       By setting CPP options in the header file *EXF_CPPOPTIONS.h* it
c       is possible to combine data sets in four different ways:
c
c       The following options are available:
c
c         ALLOW_ATM_TEMP (UAT)
c         ALLOW_ATM_WIND (UAW)
c
c
c             UAT    |    UAW      |         action
c         ----------------------------------------------------
c         undefined  |  undefined  |  Use surface fluxes.
c         undefined  |    defined  |  Assume cdn(u) given to
c                    |             |  infer the wind stress.
c           defined  |  undefined  |  Compute wind field from
c                    |             |  given stress assuming a
c                    |             |  linear relation.
c           defined  |    defined  |  Use the bulk formulae.
c         ----------------------------------------------------
c
c       Implementations of the bulk formulae exist for the follwing
c       versions of the MITgcm:
c
c       MITgcm  : Patrick Heimbach
c       MITgcmUV: Christian Eckert
c
c       started: Christian Eckert eckert@mit.edu 27-Aug-1999
c
c       changed: Christian Eckert eckert@mit.edu 14-Jan-2000
c
c              - restructured the original version in order to have a
c                better interface to the MITgcmUV.
c
c              Christian Eckert eckert@mit.edu  12-Feb-2000
c
c              - Changed Routine names (package prefix: exf_)
c
c              Patrick Heimbach, heimbach@mit.edu  04-May-2000
c
c              - changed the handling of precip and sflux with respect
c                to CPP options ALLOW_BULKFORMULAE and ALLOW_ATM_TEMP
c
c              - included some CPP flags ALLOW_BULKFORMULAE to make
c                sure ALLOW_ATM_TEMP, ALLOW_ATM_WIND are used only in
c                conjunction with defined ALLOW_BULKFORMULAE
c
c              - statement functions discarded
c
c              Ralf.Giering@FastOpt.de 25-Mai-2000
c
c              - total rewrite using new subroutines
c
c     ==================================================================
c     SUBROUTINE exf_GetFFields
c     ==================================================================

      implicit none

c     == global variables ==

#include "EEPARAMS.h"
#include "SIZE.h"
#include "PARAMS.h"
#include "DYNVARS.h"
#include "GRID.h"

#include "exf_fields.h"
#include "exf_constants.h"

#ifdef ALLOW_AUTODIFF_TAMC
#include "tamc.h"
#endif

c     == routine arguments ==

      integer mythid
      integer mycurrentiter
      _RL     mycurrenttime

c     == local variables ==

      integer bi,bj
      integer i,j,k

#ifdef ALLOW_BULKFORMULAE

#ifdef ALLOW_ATM_TEMP
      integer iter
      _RL     aln
      _RL     delq
      _RL     deltap
      _RL     hqol
      _RL     htol
      _RL     huol
      _RL     psimh
      _RL     psixh
      _RL     qstar
      _RL     rd
      _RL     re
      _RL     rdn
      _RL     rh
      _RL     ssttmp
      _RL     ssq
      _RL     stable
      _RL     tstar
      _RL     t0
      _RL     ustar
      _RL     uzn
      _RL     shn
      _RL     xsq
      _RL     x
      _RL     tau
      _RL     evap(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
#ifdef ALLOW_AUTODIFF_TAMC
      integer ikey_1
      integer ikey_2
#endif
#endif /* ALLOW_ATM_TEMP */

      _RL     ustmp
      _RL     us
      _RL     cw
      _RL     sw
      _RL     sh
      _RL     hs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
      _RL     hl(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
      _RL     hfl

#endif /* ALLOW_BULKFORMULAE */

c     == external functions ==

      integer  ilnblnk
      external ilnblnk

#ifdef ALLOW_BULKFORMULAE
      _RL       exf_BulkqSat
      external  exf_BulkqSat
      _RL       exf_BulkCdn
      external  exf_BulkCdn
      _RL       exf_BulkRhn
      external  exf_BulkRhn
#endif /* ALLOW_BULKFORMULAE */

c     == end of interface ==

c     determine forcing field records

#ifdef ALLOW_BULKFORMULAE

#if (defined (ALLOW_ATM_TEMP) || defined (ALLOW_ATM_WIND))
      aln = log(ht/zref)
#endif

c     Determine where we are in time and set counters, flags and
c     the linear interpolation factors accordingly.
#ifdef ALLOW_ATM_TEMP
c     Atmospheric temperature.
      call exf_set_atemp( mycurrenttime, mycurrentiter, mythid )

c     Atmospheric humidity.
      call exf_set_aqh( mycurrenttime, mycurrentiter, mythid )

c     Net long wave radiative flux.
      call exf_set_lwflux( mycurrenttime, mycurrentiter, mythid )

c     Net short wave radiative flux.
      call exf_set_swflux( mycurrenttime, mycurrentiter, mythid )

c     Precipitation.
      call exf_set_precip( mycurrenttime, mycurrentiter, mythid )

#ifdef ALLOW_ATEMP_CONTROL
      call ctrl_getatemp ( mycurrenttime, mycurrentiter, mythid )
#endif

#ifdef ALLOW_AQH_CONTROL
      call ctrl_getaqh ( mycurrenttime, mycurrentiter, mythid )
#endif

#else

c     Atmospheric heat flux.
      call exf_set_hflux( mycurrenttime, mycurrentiter, mythid )

c     Salt flux.
      call exf_set_sflux( mycurrenttime, mycurrentiter, mythid )

#ifdef ALLOW_KPP
c     Net short wave radiative flux.
      call exf_set_swflux( mycurrenttime, mycurrentiter, mythid )

#endif /* ALLOW_KPP */

#endif /* ALLOW_ATM_TEMP */

#ifdef ALLOW_ATM_WIND
c     Zonal wind.
      call exf_set_uwind( mycurrenttime, mycurrentiter, mythid )

c     Meridional wind.
      call exf_set_vwind( mycurrenttime, mycurrentiter, mythid )

#ifdef ALLOW_UWIND_CONTROL
      call ctrl_getuwind ( mycurrenttime, mycurrentiter, mythid )
#endif

#ifdef ALLOW_VWIND_CONTROL
      call ctrl_getvwind ( mycurrenttime, mycurrentiter, mythid )
#endif

#else
c     Zonal wind stress.
      call exf_set_ustress( mycurrenttime, mycurrentiter, mythid )

c     Meridional wind stress.
      call exf_set_vstress( mycurrenttime, mycurrentiter, mythid )

#endif /* ALLOW_ATM_WIND */

#else /* ALLOW_BULKFORMULAE undefined */
c     Atmospheric heat flux.
      call exf_set_hflux( mycurrenttime, mycurrentiter, mythid )

c     Salt flux.
      call exf_set_sflux( mycurrenttime, mycurrentiter, mythid )

c     Zonal wind stress.
      call exf_set_ustress( mycurrenttime, mycurrentiter, mythid )

c     Meridional wind stress.
      call exf_set_vstress( mycurrenttime, mycurrentiter, mythid )

#ifdef ALLOW_KPP
c     Net short wave radiative flux.
      call exf_set_swflux( mycurrenttime, mycurrentiter, mythid )
#endif

#endif /* ALLOW_BULKFORMULAE */

c     Loop over tiles.
#ifdef ALLOW_AUTODIFF_TAMC
C--   HPF directive to help TAMC
CHPF$ INDEPENDENT
#endif /* ALLOW_AUTODIFF_TAMC */
      do bj = mybylo(mythid),mybyhi(mythid)
#ifdef ALLOW_AUTODIFF_TAMC
C--    HPF directive to help TAMC
CHPF$  INDEPENDENT
#endif
        do bi = mybxlo(mythid),mybxhi(mythid)

          k = 1

          do j = 1-oly,sny+oly
            do i = 1-olx,snx+olx

#ifdef ALLOW_BULKFORMULAE

#ifdef ALLOW_AUTODIFF_TAMC
               act1 = bi - myBxLo(myThid)
               max1 = myBxHi(myThid) - myBxLo(myThid) + 1
               act2 = bj - myByLo(myThid)
               max2 = myByHi(myThid) - myByLo(myThid) + 1
               act3 = myThid - 1
               max3 = nTx*nTy
               act4 = ikey_dynamics - 1

               ikey_1 = i
     &              + sNx*(j-1)
     &              + sNx*sNy*act1
     &              + sNx*sNy*max1*act2
     &              + sNx*sNy*max1*max2*act3
     &              + sNx*sNy*max1*max2*max3*act4
#endif

c             Compute the turbulent surface fluxes.
c                  (Bulk formulae estimates)

#ifdef ALLOW_ATM_WIND
c             Wind speed and direction.
              ustmp = uwind(i,j,bi,bj)*uwind(i,j,bi,bj) +
     &                vwind(i,j,bi,bj)*vwind(i,j,bi,bj)
              if ( ustmp .ne. 0. _d 0 ) then
                us = sqrt(ustmp)
                cw = uwind(i,j,bi,bj)/us
                sw = vwind(i,j,bi,bj)/us
              else
                us = 0. _d 0
                cw = 0. _d 0
                sw = 0. _d 0
              endif
              sh = max(us,umin)
#else
#ifdef ALLOW_ATM_TEMP

c             The variables us, sh and rdn have to be computed from
c             given wind stresses inverting relationship for neutral
c             drag coeff. cdn.
c             The inversion is based on linear and quadratic form of
c             cdn(umps); ustar can be directly computed from stress;

              ustmp = ustress(i,j,bi,bj)*ustress(i,j,bi,bj) + 
     &                vstress(i,j,bi,bj)*vstress(i,j,bi,bj)
              if ( ustmp .ne. 0. _d 0 ) then
                ustar = sqrt(ustmp/atmrho)
                cw = ustress(i,j,bi,bj)/sqrt(ustmp)
                sw = vstress(i,j,bi,bj)/sqrt(ustmp)
              else
                 ustar = 0. _d 0
                 cw    = 0. _d 0
                 sw    = 0. _d 0
              endif

              if ( ustar .eq. 0. _d 0 ) then
                us = 0. _d 0
              else if ( ustar .lt. ustofu11 ) then
                tmp1 = -cquadrag_2/cquadrag_1/2
                tmp2 = sqrt(tmp1*tmp1 + ustar*ustar/cquadrag_1)
                us   = sqrt(tmp1 + tmp2)
              else
                tmp3 = clindrag_2/clindrag_1/3
                tmp4 = ustar*ustar/clindrag_1/2 - tmp3**3
                tmp5 = sqrt(ustar*ustar/clindrag_1*
     &                      (ustar*ustar/clindrag_1/4 - tmp3**3))
                us   = (tmp4 + tmp5)**(1/3) +
     &                 tmp3**2 * (tmp4 + tmp5)**(-1/3) - tmp3
              endif

              if ( us .ne. 0 ) then
                 rdn = ustar/us
              else
                 rdn = 0. _d 0
              end if

              sh    = max(us,umin)
#endif /* ALLOW_ATM_TEMP */
#endif /* ALLOW_ATM_WIND */

#ifdef ALLOW_ATM_TEMP

c             Initial guess: z/l=0.0; hu=ht=hq=z
c             Iterations:    converge on z/l and hence the fluxes.
c             t0     : virtual temperature (K)
c             ssq    : sea surface humidity (kg/kg)
c             deltap : potential temperature diff (K)

              if ( atemp(i,j,bi,bj) .ne. 0. _d 0 ) then
                t0     = atemp(i,j,bi,bj)*
     &                   (exf_one + humid_fac*aqh(i,j,bi,bj))
                ssttmp = theta(i,j,k,bi,bj)
                ssq    = saltsat*
     &                   exf_BulkqSat(ssttmp + cen2kel)/
     &                   atmrho
                deltap = atemp(i,j,bi,bj)   + gamma_blk*ht -
     &                   ssttmp - cen2kel
                delq   = aqh(i,j,bi,bj) - ssq
                stable = exf_half + sign(exf_half, deltap)
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE sh     = comlev1_exf_1, key = ikey_1
#endif
                rdn    = sqrt(exf_BulkCdn(sh))
                ustar  = rdn*sh
                tstar  = exf_BulkRhn(stable)*deltap 
                qstar  = cdalton*delq 

                do iter = 1,niter_bulk

#ifdef ALLOW_AUTODIFF_TAMC
                   ikey_2 = iter
     &                  + niter_bulk*(i-1)
     &                  + sNx*niter_bulk*(j-1)
     &                  + sNx*niter_bulk*sNy*act1
     &                  + sNx*niter_bulk*sNy*max1*act2
     &                  + sNx*niter_bulk*sNy*max1*max2*act3
     &                  + sNx*niter_bulk*sNy*max1*max2*max3*act4
#endif

#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE rdn    = comlev1_exf_2, key = ikey_2
CADJ STORE ustar  = comlev1_exf_2, key = ikey_2
CADJ STORE qstar  = comlev1_exf_2, key = ikey_2
CADJ STORE tstar  = comlev1_exf_2, key = ikey_2
CADJ STORE sh     = comlev1_exf_2, key = ikey_2
CADJ STORE us     = comlev1_exf_2, key = ikey_2
#endif

                  huol   = czol*(tstar/t0 +
     &                     qstar/(exf_one/humid_fac+aqh(i,j,bi,bj)))/
     &                           ustar**2
                  huol   = max(huol,zolmin)
                  stable = exf_half + sign(exf_half, huol)
                  htol   = huol*ht/hu
                  hqol   = huol*hq/hu

c                 Evaluate all stability functions assuming hq = ht.
                  xsq    = max(sqrt(abs(exf_one - 16.*huol)),exf_one)
                   x     = sqrt(xsq)
                  psimh  = -psim_fac*huol*stable +
     &                     (exf_one - stable)*
     &                     log((exf_one + x*(exf_two + x))*
     &                     (exf_one + xsq)/8.) - exf_two*atan(x) +
     &                     pi*exf_half
                  xsq    = max(sqrt(abs(exf_one - 16.*htol)),exf_one)
                  psixh  = -psim_fac*htol*stable + (exf_one - stable)*
     &                     exf_two*log((exf_one + xsq)/exf_two)

c                 Shift wind speed using old coefficient
ccc                  rd   = rdn/(exf_one + rdn/karman*
ccc     &                 (log(hu/zref) - psimh) )
                  rd   = rdn/(exf_one - rdn/karman*psimh )
                  shn  = sh*rd/rdn
                  uzn  = max(shn, umin)

c                 Update the transfer coefficients at 10 meters
c                 and neutral stability.

                  rdn = sqrt(exf_BulkCdn(uzn))

c                 Shift all coefficients to the measurement height
c                 and stability.
c                 rd = rdn/(exf_one + rdn/karman*(log(hu/zref) - psimh))
                  rd = rdn/(exf_one - rdn/karman*psimh)
                  rh = exf_BulkRhn(stable)/(exf_one +
     &                                      exf_BulkRhn(stable)/
     &                                     karman*(aln - psixh))
                  re = cdalton/(exf_one + cdalton/karman*(aln - psixh))

c                 Update ustar, tstar, qstar using updated, shifted
c                 coefficients.
                  ustar = rd*sh  
                  qstar = re*delq 
                  tstar = rh*deltap
                  tau   = atmrho*ustar**2
                  tau   = tau*us/sh

                enddo

#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE ustar  = comlev1_exf_1, key = ikey_1
CADJ STORE qstar  = comlev1_exf_1, key = ikey_1
CADJ STORE tstar  = comlev1_exf_1, key = ikey_1
CADJ STORE tau    = comlev1_exf_1, key = ikey_1
CADJ STORE cw     = comlev1_exf_1, key = ikey_1
CADJ STORE sw     = comlev1_exf_1, key = ikey_1
#endif

                hs(i,j,bi,bj)      = atmcp*tau*tstar/ustar
                hl(i,j,bi,bj)      = flamb*tau*qstar/ustar

                evap(i,j,bi,bj)    = tau*qstar/ustar

                ustress(i,j,bi,bj) = tau*cw
                vstress(i,j,bi,bj) = tau*sw
ce              hflux(i,j,bi,bj)   = atmcp*tau*tstar/ustar +
ce     &                             flamb*tau*qstar/ustar
              else
                ustress(i,j,bi,bj) = 0. _d 0
                vstress(i,j,bi,bj) = 0. _d 0
                hflux  (i,j,bi,bj) = 0. _d 0
                hs(i,j,bi,bj)      = 0. _d 0
                hl(i,j,bi,bj)      = 0. _d 0
              endif

#else
#ifdef ALLOW_ATM_WIND
              ustress(i,j,bi,bj) = atmrho*exf_BulkCdn(sh)*us*
     &                             uwind(i,j,bi,bj)
              vstress(i,j,bi,bj) = atmrho*exf_BulkCdn(sh)*us*
     &                             vwind(i,j,bi,bj)
#endif /* ALLOW_ATM_WIND */
#endif /* ALLOW_ATM_TEMP */
            enddo
          enddo
        enddo
      enddo

c     Add all contributions.
      do bj = mybylo(mythid),mybyhi(mythid)
        do bi = mybxlo(mythid),mybxhi(mythid)
          do j = 1,sny
            do i = 1,snx
c             Net surface heat flux.
#ifdef ALLOW_ATM_TEMP
              hfl = 0. _d 0
              hfl = hfl - hs(i,j,bi,bj)
              hfl = hfl - hl(i,j,bi,bj)
              hfl = hfl + lwflux(i,j,bi,bj)
#ifndef ALLOW_KPP
              hfl = hfl + swflux(i,j,bi,bj)
#endif /* ALLOW_KPP undef */
c             Heat flux:
              hflux(i,j,bi,bj) = hfl*maskc(i,j,1,bi,bj)
c             Salt flux from Precipitation and Evaporation.
              sflux(i,j,bi,bj) = precip(i,j,bi,bj) - evap(i,j,bi,bj)
#endif /* ALLOW_ATM_TEMP */

#else
              hflux(i,j,bi,bj) = hflux(i,j,bi,bj)*maskc(i,j,1,bi,bj)
              sflux(i,j,bi,bj) = sflux(i,j,bi,bj)*maskc(i,j,1,bi,bj)
#endif /* ALLOW_BULKFORMULAE */
            enddo
          enddo
        enddo
      enddo

c     Update the tile edges.
      _EXCH_XY_R8(hflux,   mythid)
      _EXCH_XY_R8(sflux,   mythid)
      _EXCH_XY_R8(ustress, mythid)
      _EXCH_XY_R8(vstress, mythid)

#ifdef ALLOW_KPP
      _EXCH_XY_R8(swflux, mythid)
#endif /* ALLOW_KPP */

      end
1	heimbach	1.3	c $Header: /u/gcmpack/development/heimbach/ecco_env/pkg/exf/exf_getffields.F,v 1.13 2002/01/11 17:34:46 heimbach Exp $
2	heimbach	1.1
3			#include "EXF_CPPOPTIONS.h"
4
5	heimbach	1.3	subroutine exf_GetFFields( mycurrenttime, mycurrentiter, mythid )
6	heimbach	1.1
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
97			c - restructured the original version in order to have a
98			c better interface to the MITgcmUV.
99			c
100			c Christian Eckert eckert@mit.edu 12-Feb-2000
101			c
102			c - Changed Routine names (package prefix: exf_)
103			c
104			c Patrick Heimbach, heimbach@mit.edu 04-May-2000
105			c
106			c - changed the handling of precip and sflux with respect
107			c to CPP options ALLOW_BULKFORMULAE and ALLOW_ATM_TEMP
108			c
109			c - included some CPP flags ALLOW_BULKFORMULAE to make
110			c sure ALLOW_ATM_TEMP, ALLOW_ATM_WIND are used only in
111			c conjunction with defined ALLOW_BULKFORMULAE
112			c
113			c - statement functions discarded
114			c
115			c Ralf.Giering@FastOpt.de 25-Mai-2000
116			c
117			c - total rewrite using new subroutines
118			c
119			c ==================================================================
120			c SUBROUTINE exf_GetFFields
121			c ==================================================================
122
123			implicit none
124
125			c == global variables ==
126
127			#include "EEPARAMS.h"
128			#include "SIZE.h"
129			#include "PARAMS.h"
130			#include "DYNVARS.h"
131			#include "GRID.h"
132
133			#include "exf_fields.h"
134			#include "exf_constants.h"
135
136	heimbach	1.3	#ifdef ALLOW_AUTODIFF_TAMC
137			#include "tamc.h"
138			#endif
139
140	heimbach	1.1	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	heimbach	1.3	_RL ssttmp
169	heimbach	1.1	_RL ssq
170			_RL stable
171			_RL tstar
172			_RL t0
173			_RL ustar
174			_RL uzn
175	heimbach	1.3	_RL shn
176	heimbach	1.1	_RL xsq
177			_RL x
178	heimbach	1.3	_RL tau
179	heimbach	1.1	_RL evap(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
180	heimbach	1.3	#ifdef ALLOW_AUTODIFF_TAMC
181			integer ikey_1
182			integer ikey_2
183			#endif
184	heimbach	1.1	#endif /* ALLOW_ATM_TEMP */
185
186	heimbach	1.3	_RL ustmp
187	heimbach	1.1	_RL us
188			_RL cw
189			_RL sw
190			_RL sh
191			_RL hs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
192			_RL hl(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
193			_RL hfl
194
195			#endif /* ALLOW_BULKFORMULAE */
196
197			c == external functions ==
198
199			integer ilnblnk
200			external ilnblnk
201
202			#ifdef ALLOW_BULKFORMULAE
203			_RL exf_BulkqSat
204			external exf_BulkqSat
205			_RL exf_BulkCdn
206			external exf_BulkCdn
207			_RL exf_BulkRhn
208			external exf_BulkRhn
209			#endif /* ALLOW_BULKFORMULAE */
210
211			c == end of interface ==
212
213			c determine forcing field records
214
215			#ifdef ALLOW_BULKFORMULAE
216
217	heimbach	1.3	#if (defined (ALLOW_ATM_TEMP) \|\| defined (ALLOW_ATM_WIND))
218			aln = log(ht/zref)
219			#endif
220
221	heimbach	1.1	c Determine where we are in time and set counters, flags and
222			c the linear interpolation factors accordingly.
223			#ifdef ALLOW_ATM_TEMP
224			c Atmospheric temperature.
225	heimbach	1.3	call exf_set_atemp( mycurrenttime, mycurrentiter, mythid )
226	heimbach	1.1
227			c Atmospheric humidity.
228	heimbach	1.3	call exf_set_aqh( mycurrenttime, mycurrentiter, mythid )
229	heimbach	1.1
230			c Net long wave radiative flux.
231	heimbach	1.3	call exf_set_lwflux( mycurrenttime, mycurrentiter, mythid )
232	heimbach	1.1
233			c Net short wave radiative flux.
234	heimbach	1.3	call exf_set_swflux( mycurrenttime, mycurrentiter, mythid )
235	heimbach	1.1
236			c Precipitation.
237	heimbach	1.3	call exf_set_precip( mycurrenttime, mycurrentiter, mythid )
238	heimbach	1.1
239	heimbach	1.3	#ifdef ALLOW_ATEMP_CONTROL
240			call ctrl_getatemp ( mycurrenttime, mycurrentiter, mythid )
241			#endif
242
243			#ifdef ALLOW_AQH_CONTROL
244			call ctrl_getaqh ( mycurrenttime, mycurrentiter, mythid )
245			#endif
246	heimbach	1.1
247			#else
248
249			c Atmospheric heat flux.
250	heimbach	1.3	call exf_set_hflux( mycurrenttime, mycurrentiter, mythid )
251	heimbach	1.1
252			c Salt flux.
253	heimbach	1.3	call exf_set_sflux( mycurrenttime, mycurrentiter, mythid )
254	heimbach	1.1
255			#ifdef ALLOW_KPP
256			c Net short wave radiative flux.
257	heimbach	1.3	call exf_set_swflux( mycurrenttime, mycurrentiter, mythid )
258	heimbach	1.1
259			#endif /* ALLOW_KPP */
260
261			#endif /* ALLOW_ATM_TEMP */
262
263			#ifdef ALLOW_ATM_WIND
264			c Zonal wind.
265	heimbach	1.3	call exf_set_uwind( mycurrenttime, mycurrentiter, mythid )
266	heimbach	1.1
267			c Meridional wind.
268	heimbach	1.3	call exf_set_vwind( mycurrenttime, mycurrentiter, mythid )
269
270			#ifdef ALLOW_UWIND_CONTROL
271			call ctrl_getuwind ( mycurrenttime, mycurrentiter, mythid )
272			#endif
273
274			#ifdef ALLOW_VWIND_CONTROL
275			call ctrl_getvwind ( mycurrenttime, mycurrentiter, mythid )
276			#endif
277	heimbach	1.1
278			#else
279			c Zonal wind stress.
280	heimbach	1.3	call exf_set_ustress( mycurrenttime, mycurrentiter, mythid )
281	heimbach	1.1
282			c Meridional wind stress.
283	heimbach	1.3	call exf_set_vstress( mycurrenttime, mycurrentiter, mythid )
284	heimbach	1.1
285			#endif /* ALLOW_ATM_WIND */
286
287			#else /* ALLOW_BULKFORMULAE undefined */
288			c Atmospheric heat flux.
289	heimbach	1.3	call exf_set_hflux( mycurrenttime, mycurrentiter, mythid )
290	heimbach	1.1
291			c Salt flux.
292	heimbach	1.3	call exf_set_sflux( mycurrenttime, mycurrentiter, mythid )
293	heimbach	1.1
294			c Zonal wind stress.
295	heimbach	1.3	call exf_set_ustress( mycurrenttime, mycurrentiter, mythid )
296	heimbach	1.1
297			c Meridional wind stress.
298	heimbach	1.3	call exf_set_vstress( mycurrenttime, mycurrentiter, mythid )
299	heimbach	1.1
300			#ifdef ALLOW_KPP
301			c Net short wave radiative flux.
302	heimbach	1.3	call exf_set_swflux( mycurrenttime, mycurrentiter, mythid )
303			#endif
304	heimbach	1.1
305			#endif /* ALLOW_BULKFORMULAE */
306
307			c Loop over tiles.
308	heimbach	1.3	#ifdef ALLOW_AUTODIFF_TAMC
309			C-- HPF directive to help TAMC
310			CHPF$ INDEPENDENT
311			#endif /* ALLOW_AUTODIFF_TAMC */
312	heimbach	1.1	do bj = mybylo(mythid),mybyhi(mythid)
313	heimbach	1.3	#ifdef ALLOW_AUTODIFF_TAMC
314			C-- HPF directive to help TAMC
315			CHPF$ INDEPENDENT
316			#endif
317	heimbach	1.1	do bi = mybxlo(mythid),mybxhi(mythid)
318	heimbach	1.3
319	heimbach	1.1	k = 1
320
321			do j = 1-oly,sny+oly
322			do i = 1-olx,snx+olx
323
324			#ifdef ALLOW_BULKFORMULAE
325
326	heimbach	1.3	#ifdef ALLOW_AUTODIFF_TAMC
327			act1 = bi - myBxLo(myThid)
328			max1 = myBxHi(myThid) - myBxLo(myThid) + 1
329			act2 = bj - myByLo(myThid)
330			max2 = myByHi(myThid) - myByLo(myThid) + 1
331			act3 = myThid - 1
332			max3 = nTx*nTy
333			act4 = ikey_dynamics - 1
334
335			ikey_1 = i
336			& + sNx*(j-1)
337			& + sNxsNyact1
338			& + sNxsNymax1*act2
339			& + sNxsNymax1max2act3
340			& + sNxsNymax1max2max3*act4
341			#endif
342
343	heimbach	1.1	c Compute the turbulent surface fluxes.
344			c (Bulk formulae estimates)
345
346			#ifdef ALLOW_ATM_WIND
347			c Wind speed and direction.
348	heimbach	1.3	ustmp = uwind(i,j,bi,bj)*uwind(i,j,bi,bj) +
349			& vwind(i,j,bi,bj)*vwind(i,j,bi,bj)
350			if ( ustmp .ne. 0. _d 0 ) then
351			us = sqrt(ustmp)
352	heimbach	1.1	cw = uwind(i,j,bi,bj)/us
353			sw = vwind(i,j,bi,bj)/us
354			else
355	heimbach	1.3	us = 0. _d 0
356	heimbach	1.1	cw = 0. _d 0
357			sw = 0. _d 0
358			endif
359			sh = max(us,umin)
360			#else
361			#ifdef ALLOW_ATM_TEMP
362
363	heimbach	1.3	c The variables us, sh and rdn have to be computed from
364			c given wind stresses inverting relationship for neutral
365			c drag coeff. cdn.
366	heimbach	1.1	c The inversion is based on linear and quadratic form of
367			c cdn(umps); ustar can be directly computed from stress;
368
369	heimbach	1.3	ustmp = ustress(i,j,bi,bj)*ustress(i,j,bi,bj) +
370			& vstress(i,j,bi,bj)*vstress(i,j,bi,bj)
371			if ( ustmp .ne. 0. _d 0 ) then
372			ustar = sqrt(ustmp/atmrho)
373			cw = ustress(i,j,bi,bj)/sqrt(ustmp)
374			sw = vstress(i,j,bi,bj)/sqrt(ustmp)
375			else
376			ustar = 0. _d 0
377			cw = 0. _d 0
378			sw = 0. _d 0
379			endif
380	heimbach	1.1
381			if ( ustar .eq. 0. _d 0 ) then
382			us = 0. _d 0
383			else if ( ustar .lt. ustofu11 ) then
384			tmp1 = -cquadrag_2/cquadrag_1/2
385			tmp2 = sqrt(tmp1tmp1 + ustarustar/cquadrag_1)
386			us = sqrt(tmp1 + tmp2)
387			else
388			tmp3 = clindrag_2/clindrag_1/3
389			tmp4 = ustarustar/clindrag_1/2 - tmp3*3
390			tmp5 = sqrt(ustarustar/clindrag_1
391			& (ustarustar/clindrag_1/4 - tmp3*3))
392			us = (tmp4 + tmp5)**(1/3) +
393			& tmp3*2 (tmp4 + tmp5)**(-1/3) - tmp3
394			endif
395
396			if ( us .ne. 0 ) then
397			rdn = ustar/us
398			else
399			rdn = 0. _d 0
400			end if
401
402			sh = max(us,umin)
403			#endif /* ALLOW_ATM_TEMP */
404			#endif /* ALLOW_ATM_WIND */
405
406			#ifdef ALLOW_ATM_TEMP
407	heimbach	1.3
408	heimbach	1.1	c Initial guess: z/l=0.0; hu=ht=hq=z
409			c Iterations: converge on z/l and hence the fluxes.
410			c t0 : virtual temperature (K)
411			c ssq : sea surface humidity (kg/kg)
412			c deltap : potential temperature diff (K)
413
414			if ( atemp(i,j,bi,bj) .ne. 0. _d 0 ) then
415			t0 = atemp(i,j,bi,bj)*
416			& (exf_one + humid_fac*aqh(i,j,bi,bj))
417	heimbach	1.3	ssttmp = theta(i,j,k,bi,bj)
418	heimbach	1.1	ssq = saltsat*
419	heimbach	1.3	& exf_BulkqSat(ssttmp + cen2kel)/
420	heimbach	1.1	& atmrho
421			deltap = atemp(i,j,bi,bj) + gamma_blk*ht -
422	heimbach	1.3	& ssttmp - cen2kel
423	heimbach	1.1	delq = aqh(i,j,bi,bj) - ssq
424			stable = exf_half + sign(exf_half, deltap)
425	heimbach	1.3	#ifdef ALLOW_AUTODIFF_TAMC
426			CADJ STORE sh = comlev1_exf_1, key = ikey_1
427			#endif
428	heimbach	1.1	rdn = sqrt(exf_BulkCdn(sh))
429			ustar = rdn*sh
430			tstar = exf_BulkRhn(stable)*deltap
431			qstar = cdalton*delq
432
433			do iter = 1,niter_bulk
434
435	heimbach	1.3	#ifdef ALLOW_AUTODIFF_TAMC
436			ikey_2 = iter
437			& + niter_bulk*(i-1)
438			& + sNxniter_bulk(j-1)
439			& + sNxniter_bulksNy*act1
440			& + sNxniter_bulksNymax1act2
441			& + sNxniter_bulksNymax1max2*act3
442			& + sNxniter_bulksNymax1max2max3act4
443			#endif
444
445			#ifdef ALLOW_AUTODIFF_TAMC
446			CADJ STORE rdn = comlev1_exf_2, key = ikey_2
447			CADJ STORE ustar = comlev1_exf_2, key = ikey_2
448			CADJ STORE qstar = comlev1_exf_2, key = ikey_2
449			CADJ STORE tstar = comlev1_exf_2, key = ikey_2
450			CADJ STORE sh = comlev1_exf_2, key = ikey_2
451			CADJ STORE us = comlev1_exf_2, key = ikey_2
452			#endif
453
454	heimbach	1.1	huol = czol*(tstar/t0 +
455			& qstar/(exf_one/humid_fac+aqh(i,j,bi,bj)))/
456			& ustar**2
457			huol = max(huol,zolmin)
458			stable = exf_half + sign(exf_half, huol)
459			htol = huol*ht/hu
460			hqol = huol*hq/hu
461
462			c Evaluate all stability functions assuming hq = ht.
463			xsq = max(sqrt(abs(exf_one - 16.*huol)),exf_one)
464			x = sqrt(xsq)
465			psimh = -psim_fachuolstable +
466			& (exf_one - stable)*
467			& log((exf_one + x(exf_two + x))
468			& (exf_one + xsq)/8.) - exf_two*atan(x) +
469			& pi*exf_half
470			xsq = max(sqrt(abs(exf_one - 16.*htol)),exf_one)
471			psixh = -psim_fachtolstable + (exf_one - stable)*
472			& exf_two*log((exf_one + xsq)/exf_two)
473
474			c Shift wind speed using old coefficient
475	heimbach	1.3	ccc rd = rdn/(exf_one + rdn/karman*
476			ccc & (log(hu/zref) - psimh) )
477	heimbach	1.1	rd = rdn/(exf_one - rdn/karman*psimh )
478	heimbach	1.3	shn = sh*rd/rdn
479			uzn = max(shn, umin)
480	heimbach	1.1
481			c Update the transfer coefficients at 10 meters
482			c and neutral stability.
483	heimbach	1.3
484	heimbach	1.1	rdn = sqrt(exf_BulkCdn(uzn))
485
486			c Shift all coefficients to the measurement height
487			c and stability.
488			c rd = rdn/(exf_one + rdn/karman*(log(hu/zref) - psimh))
489			rd = rdn/(exf_one - rdn/karman*psimh)
490			rh = exf_BulkRhn(stable)/(exf_one +
491			& exf_BulkRhn(stable)/
492			& karman*(aln - psixh))
493			re = cdalton/(exf_one + cdalton/karman*(aln - psixh))
494
495			c Update ustar, tstar, qstar using updated, shifted
496			c coefficients.
497			ustar = rd*sh
498			qstar = re*delq
499			tstar = rh*deltap
500	heimbach	1.3	tau = atmrhoustar*2
501			tau = tau*us/sh
502	heimbach	1.1
503			enddo
504
505	heimbach	1.3	#ifdef ALLOW_AUTODIFF_TAMC
506			CADJ STORE ustar = comlev1_exf_1, key = ikey_1
507			CADJ STORE qstar = comlev1_exf_1, key = ikey_1
508			CADJ STORE tstar = comlev1_exf_1, key = ikey_1
509			CADJ STORE tau = comlev1_exf_1, key = ikey_1
510			CADJ STORE cw = comlev1_exf_1, key = ikey_1
511			CADJ STORE sw = comlev1_exf_1, key = ikey_1
512			#endif
513
514	heimbach	1.1	hs(i,j,bi,bj) = atmcptautstar/ustar
515			hl(i,j,bi,bj) = flambtauqstar/ustar
516
517			evap(i,j,bi,bj) = tau*qstar/ustar
518
519			ustress(i,j,bi,bj) = tau*cw
520			vstress(i,j,bi,bj) = tau*sw
521			ce hflux(i,j,bi,bj) = atmcptautstar/ustar +
522			ce & flambtauqstar/ustar
523			else
524			ustress(i,j,bi,bj) = 0. _d 0
525			vstress(i,j,bi,bj) = 0. _d 0
526			hflux (i,j,bi,bj) = 0. _d 0
527			hs(i,j,bi,bj) = 0. _d 0
528			hl(i,j,bi,bj) = 0. _d 0
529			endif
530
531			#else
532			#ifdef ALLOW_ATM_WIND
533			ustress(i,j,bi,bj) = atmrhoexf_BulkCdn(sh)us*
534			& uwind(i,j,bi,bj)
535			vstress(i,j,bi,bj) = atmrhoexf_BulkCdn(sh)us*
536			& vwind(i,j,bi,bj)
537			#endif /* ALLOW_ATM_WIND */
538			#endif /* ALLOW_ATM_TEMP */
539			enddo
540			enddo
541			enddo
542			enddo
543
544			c Add all contributions.
545			do bj = mybylo(mythid),mybyhi(mythid)
546			do bi = mybxlo(mythid),mybxhi(mythid)
547			do j = 1,sny
548			do i = 1,snx
549			c Net surface heat flux.
550			#ifdef ALLOW_ATM_TEMP
551			hfl = 0. _d 0
552			hfl = hfl - hs(i,j,bi,bj)
553			hfl = hfl - hl(i,j,bi,bj)
554			hfl = hfl + lwflux(i,j,bi,bj)
555			#ifndef ALLOW_KPP
556			hfl = hfl + swflux(i,j,bi,bj)
557			#endif /* ALLOW_KPP undef */
558			c Heat flux:
559	heimbach	1.3	hflux(i,j,bi,bj) = hfl*maskc(i,j,1,bi,bj)
560	heimbach	1.1	c Salt flux from Precipitation and Evaporation.
561			sflux(i,j,bi,bj) = precip(i,j,bi,bj) - evap(i,j,bi,bj)
562			#endif /* ALLOW_ATM_TEMP */
563
564			#else
565	heimbach	1.3	hflux(i,j,bi,bj) = hflux(i,j,bi,bj)*maskc(i,j,1,bi,bj)
566			sflux(i,j,bi,bj) = sflux(i,j,bi,bj)*maskc(i,j,1,bi,bj)
567	heimbach	1.1	#endif /* ALLOW_BULKFORMULAE */
568			enddo
569			enddo
570			enddo
571			enddo
572
573			c Update the tile edges.
574			_EXCH_XY_R8(hflux, mythid)
575			_EXCH_XY_R8(sflux, mythid)
576			_EXCH_XY_R8(ustress, mythid)
577			_EXCH_XY_R8(vstress, mythid)
578
579			#ifdef ALLOW_KPP
580			_EXCH_XY_R8(swflux, mythid)
581			#endif /* ALLOW_KPP */
582
583			end