pkg/exf/exf_getffields.F

c $Header: /u/gcmpack/MITgcm/pkg/exf/exf_getffields.F,v 1.2.4.9 2003/03/03 06:58:38 dimitri Exp $

#include "EXF_CPPOPTIONS.h"

      subroutine exf_GetFFields( mycurrenttime, mycurrentiter, mythid )

c     ==================================================================
c     SUBROUTINE exf_GetFFields
c     ==================================================================
c
c     o Read-in atmospheric state and/or surface fluxes from files.
c
c     o Use bulk formulae to estimate turbulent and/or radiative
c       fluxes at the surface.
c
c     NOTES:
c     ======
c
c     See EXF_CPPOPTIONS.h for a description of the various possible
c     ocean-model forcing configurations.
c
c     The bulk formulae of pkg/exf are not valid for sea-ice covered
c     oceans but they can be used in combination with a sea-ice model,
c     for example, pkg/seaice, to specify open water flux contributions.
c
c     ==================================================================
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       started: Christian Eckert eckert@mit.edu 27-Aug-1999
c
c       changed: Christian Eckert eckert@mit.edu 14-Jan-2000
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              - Changed Routine names (package prefix: exf_)
c
c              Patrick Heimbach, heimbach@mit.edu  04-May-2000
c              - changed the handling of precip and sflux with respect
c                to CPP options ALLOW_BULKFORMULAE and ALLOW_ATM_TEMP
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              - statement functions discarded
c
c              Ralf.Giering@FastOpt.de 25-Mai-2000
c              - total rewrite using new subroutines
c
c              Detlef Stammer: include river run-off. Nov. 21, 2001
c
c              heimbach@mit.edu, 10-Jan-2002
c              - changes to enable field swapping
c
c       mods for pkg/seaice: menemenlis@jpl.nasa.gov 20-Dec-2002
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_param.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

      _RL     aln

#ifdef ALLOW_ATM_TEMP
      integer iter
      _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
#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 */

#ifndef ALLOW_ATM_WIND
      _RL       TMP1
      _RL       TMP2
      _RL       TMP3
      _RL       TMP4
      _RL       TMP5
#endif

c     == end of interface ==

#ifdef ALLOW_BULKFORMULAE
cph   This statement cannot be a PARAMETER statement in the header,
cph   but must come here; it's not fortran77 standard
      aln = log(ht/zref)
#endif

c--   read forcing fields from files and temporal interpolation

#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  /* ifndef ALLOW_ATM_WIND */

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

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

#endif /* ifndef ALLOW_ATM_WIND */

#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     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  /* ifndef ALLOW_ATM_TEMP */

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

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

#endif /* ifndef ALLOW_ATM_TEMP */

#if defined(ALLOW_ATM_TEMP) || defined(SHORTWAVE_HEATING)
c     Net short wave radiative flux.
      call exf_set_swflux ( mycurrenttime, mycurrentiter, mythid )
#endif

#ifdef EXF_READ_EVAP
c     Evaporation
      call exf_set_evap   ( mycurrenttime, mycurrentiter, mythid )
#endif

#ifdef ALLOW_DOWNWARD_RADIATION

c     Downward shortwave radiation.
      call exf_set_swdown ( mycurrenttime, mycurrentiter, mythid )

c     Downward longwave radiation.
      call exf_set_lwdown ( mycurrenttime, mycurrentiter, mythid )

#endif

#ifdef ATMOSPHERIC_LOADING
c     Atmos. pressure forcing
      call exf_set_apressure ( mycurrenttime, mycurrentiter, mythid )
#endif

c--   Use atmospheric state to compute surface fluxes.

c     Loop over tiles.
#ifdef ALLOW_AUTODIFF_TAMC
C--   HPF directive to help TAMC
CHPF$ INDEPENDENT
#endif
      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,sny
            do i = 1,snx

#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

#ifdef ALLOW_DOWNWARD_RADIATION
c--   Compute net from downward and downward from net longwave and
c     shortwave radiation, if needed.
c     lwflux = Stefan-Boltzman constant * emissivity * SST - lwdown
c     swflux = - ( 1 - albedo ) * swdown

#ifdef ALLOW_ATM_TEMP
      if ( lwfluxfile .EQ. ' ' .AND. lwdownfile .NE. ' ' )
     &     lwflux(i,j,bi,bj) = 5.5 _d -08 *
     &              ((theta(i,j,k,bi,bj)+cen2kel)**4)
     &              - lwdown(i,j,bi,bj)
      if ( lwfluxfile .NE. ' ' .AND. lwdownfile .EQ. ' ' )
     &     lwdown(i,j,bi,bj) = 5.5 _d -08 *
     &              ((theta(i,j,k,bi,bj)+cen2kel)**4)
     &              - lwflux(i,j,bi,bj)
#endif

#if defined(ALLOW_ATM_TEMP) || defined(SHORTWAVE_HEATING)
      if ( swfluxfile .EQ. ' ' .AND. swdownfile .NE. ' ' )
     &     swflux(i,j,bi,bj) = -0.9 _d 0 * swdown(i,j,bi,bj)
      if ( swfluxfile .NE. ' ' .AND. swdownfile .EQ. ' ' )
     &     swdown(i,j,bi,bj) = -1.111111 _d 0 * swflux(i,j,bi,bj)
#endif

#endif /* ALLOW_DOWNWARD_RADIATION */

c--   Compute the turbulent surface fluxes.

#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  /* ifndef ALLOW_ATM_WIND */
#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 /* ifndef 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

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
#ifndef EXF_READ_EVAP
cdm             evap(i,j,bi,bj)    = tau*qstar/ustar
cdm !!! need to change sign and to convert from kg/m^2/s to m/s !!!
                evap(i,j,bi,bj)    = -recip_rhonil*tau*qstar/ustar
#endif
                ustress(i,j,bi,bj) = tau*cw
                vstress(i,j,bi,bj) = tau*sw
              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  /* ifndef ALLOW_ATM_TEMP */
#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
#endif /* ifndef 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 SHORTWAVE_HEATING
              hfl = hfl + swflux(i,j,bi,bj)
#endif
c             Heat flux:
              hflux(i,j,bi,bj) = hfl
c             Salt flux from Precipitation and Evaporation.
              sflux(i,j,bi,bj) = evap(i,j,bi,bj) - precip(i,j,bi,bj)
#endif /* ALLOW_ATM_TEMP */

#endif /* ALLOW_BULKFORMULAE */

#ifdef ALLOW_RUNOFF
              sflux(i,j,bi,bj) = sflux(i,j,bi,bj) - runoff(i,j,bi,bj)
#endif

              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)

            enddo
          enddo
        enddo
      enddo

c     Update the tile edges.
      _EXCH_XY_R8(hflux,   mythid)
      _EXCH_XY_R8(sflux,   mythid)
c      _EXCH_XY_R8(ustress, mythid)
c      _EXCH_XY_R8(vstress, mythid)
       CALL EXCH_UV_XY_RL(ustress, vstress, .TRUE., myThid)

#ifdef SHORTWAVE_HEATING
      _EXCH_XY_R8(swflux, mythid)
#endif

#ifdef ATMOSPHERIC_LOADING
      _EXCH_XY_R8(apressure, mythid)
#endif

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