pkg/exf/exf_bulkformulae.F

C $Header: /u/gcmpack/MITgcm/pkg/exf/exf_bulkformulae.F,v 1.15 2007/04/16 23:27:20 jmc Exp $
C $Name:  $

#include "EXF_OPTIONS.h"

      subroutine exf_bulkformulae(mytime, myiter, mythid)

c     ==================================================================
c     SUBROUTINE exf_bulkformulae
c     ==================================================================
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_OPTIONS.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_bulkformulae
c     ==================================================================

      implicit none

c     == global variables ==

#include "EEPARAMS.h"
#include "SIZE.h"
#include "PARAMS.h"
#include "DYNVARS.h"
c #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 myiter
      _RL     mytime

#ifdef ALLOW_BULKFORMULAE

c     == local variables ==

      integer bi,bj
      integer i,j,k

      _RL aln
      _RL czol

      integer iter
      _RL     tmpbulk
      _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

c     == external functions ==

      integer  ilnblnk
      external ilnblnk

      _RL       exf_BulkqSat
      external  exf_BulkqSat
      _RL       exf_BulkCdn
      external  exf_BulkCdn
      _RL       exf_BulkRhn
      external  exf_BulkRhn

c     == end of interface ==

      aln = log(ht/zref)
      czol = hu*karman*gravity_mks

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_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.

#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)
             tmpbulk= exf_BulkqSat(ssttmp + cen2kel)
             ssq    = saltsat*tmpbulk/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(i,j,bi,bj)     = comlev1_exf_1, key = ikey_1
#endif
#ifdef ALLOW_ATM_WIND
             tmpbulk= exf_BulkCdn(sh(i,j,bi,bj))
             rdn    = sqrt(tmpbulk)
             ustar  = rdn*sh(i,j,bi,bj)
#else /* ifndef ALLOW_ATM_WIND */
C     in this case ustress and vstress are defined a u and v points
C     respectively, and we need to average to mass points to avoid
C     tau = 0 near coasts or other boundaries
             tau   = sqrt(0.5*
     &                   (ustress(i,  j,bi,bj)*ustress(i  ,j,bi,bj)
     &                   +ustress(i+1,j,bi,bj)*ustress(i+1,j,bi,bj)
     &                   +vstress(i,j,  bi,bj)*vstress(i,j  ,bi,bj)
     &                   +vstress(i,j+1,bi,bj)*vstress(i,j+1,bi,bj))
     &                   )
             ustar = sqrt(tau/atmrho)
#endif /* ifndef ALLOW_ATM_WIND */
             tmpbulk= exf_BulkRhn(stable)
             tstar  = tmpbulk*deltap
             qstar  = cdalton*delq

             do iter = 1,niter_bulk

#ifdef ALLOW_AUTODIFF_TAMC
                   ikey_2 = iter
     &                  + niter_bulk*(i-1)
     &                  + niter_bulk*sNx*(j-1)
     &                  + niter_bulk*sNx*sNy*act1
     &                  + niter_bulk*sNx*sNy*max1*act2
     &                  + niter_bulk*sNx*sNy*max1*max2*act3
     &                  + niter_bulk*sNx*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(i,j,bi,bj)     = comlev1_exf_2, key = ikey_2
CADJ STORE us(i,j,bi,bj)     = 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.
#ifdef ALLOW_ATM_WIND
                  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)
#endif /* ALLOW_ATM_WIND */
                  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)

#ifdef ALLOW_ATM_WIND
c                 Shift wind speed using old coefficient
                  rd     = rdn/(exf_one - rdn/karman*psimh )
                  shn    = sh(i,j,bi,bj)*rd/rdn
                  uzn    = max(shn, umin)

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

                  tmpbulk= exf_BulkCdn(uzn)
                  rdn    = sqrt(tmpbulk)

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

c                 Update ustar, tstar, qstar using updated, shifted
c                 coefficients.
                  qstar  = re*delq
                  tstar  = rh*deltap
#ifdef ALLOW_ATM_WIND
                  ustar  = rd*sh(i,j,bi,bj)
                  tau    = atmrho*ustar**2
                  tau    = tau*us(i,j,bi,bj)/sh(i,j,bi,bj)
#endif /* ALLOW_ATM_WIND */

                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(i,j,bi,bj)     = comlev1_exf_1, key = ikey_1
CADJ STORE sw(i,j,bi,bj)     = 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
#ifdef ALLOW_ATM_WIND
                ustress(i,j,bi,bj) = tau*cw(i,j,bi,bj)
                vstress(i,j,bi,bj) = tau*sw(i,j,bi,bj)
#endif /* ALLOW_ATM_WIND */
              else
#ifdef ALLOW_ATM_WIND
                ustress(i,j,bi,bj) = 0. _d 0
                vstress(i,j,bi,bj) = 0. _d 0
#endif /* ALLOW_ATM_WIND */
                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
              tmpbulk= exf_BulkCdn(sh(i,j,bi,bj))
              ustress(i,j,bi,bj) = atmrho*tmpbulk*us(i,j,bi,bj)*
     &                             uwind(i,j,bi,bj)
              vstress(i,j,bi,bj) = atmrho*tmpbulk*us(i,j,bi,bj)*
     &                             vwind(i,j,bi,bj)
#endif
#endif /* ifndef ALLOW_ATM_TEMP */
            enddo
          enddo
        enddo
      enddo

#endif /* ALLOW_BULKFORMULAE */

      return
      end
1	heimbach	1.16	C $Header: /u/gcmpack/MITgcm/pkg/exf/exf_bulkformulae.F,v 1.15 2007/04/16 23:27:20 jmc Exp $
2	jmc	1.13	C $Name: $
3	heimbach	1.2
4	edhill	1.4	#include "EXF_OPTIONS.h"
5	heimbach	1.2
6	heimbach	1.12	subroutine exf_bulkformulae(mytime, myiter, mythid)
7	heimbach	1.2
8			c ==================================================================
9			c SUBROUTINE exf_bulkformulae
10			c ==================================================================
11			c
12			c o Use bulk formulae to estimate turbulent and/or radiative
13			c fluxes at the surface.
14			c
15			c NOTES:
16			c ======
17			c
18	edhill	1.4	c See EXF_OPTIONS.h for a description of the various possible
19	heimbach	1.2	c ocean-model forcing configurations.
20			c
21			c The bulk formulae of pkg/exf are not valid for sea-ice covered
22			c oceans but they can be used in combination with a sea-ice model,
23			c for example, pkg/seaice, to specify open water flux contributions.
24			c
25			c ==================================================================
26			c
27			c The calculation of the bulk surface fluxes has been adapted from
28			c the NCOM model which uses the formulae given in Large and Pond
29			c (1981 & 1982 )
30			c
31			c
32			c Header taken from NCOM version: ncom1.4.1
33			c -----------------------------------------
34			c
35			c Following procedures and coefficients in Large and Pond
36			c (1981 ; 1982)
37			c
38			c Output: Bulk estimates of the turbulent surface fluxes.
39			c -------
40			c
41			c hs - sensible heat flux (W/m^2), into ocean
42			c hl - latent heat flux (W/m^2), into ocean
43			c
44			c Input:
45			c ------
46			c
47			c us - mean wind speed (m/s) at height hu (m)
48			c th - mean air temperature (K) at height ht (m)
49			c qh - mean air humidity (kg/kg) at height hq (m)
50			c sst - sea surface temperature (K)
51			c tk0 - Kelvin temperature at 0 Celsius (K)
52			c
53			c Assume 1) a neutral 10m drag coefficient =
54			c
55			c cdn = .0027/u10 + .000142 + .0000764 u10
56			c
57			c 2) a neutral 10m stanton number =
58			c
59			c ctn = .0327 sqrt(cdn), unstable
60			c ctn = .0180 sqrt(cdn), stable
61			c
62			c 3) a neutral 10m dalton number =
63			c
64			c cen = .0346 sqrt(cdn)
65			c
66			c 4) the saturation humidity of air at
67			c
68			c t(k) = exf_BulkqSat(t) (kg/m^3)
69			c
70			c Note: 1) here, tstar = <wt>/u, and qstar = <wq>/u.
71			c 2) wind speeds should all be above a minimum speed,
72			c say 0.5 m/s.
73			c 3) with optional iteration loop, niter=3, should suffice.
74			c 4) this version is for analyses inputs with hu = 10m and
75			c ht = hq.
76			c 5) sst enters in Celsius.
77			c
78			c ==================================================================
79			c
80			c started: Christian Eckert eckert@mit.edu 27-Aug-1999
81			c
82			c changed: Christian Eckert eckert@mit.edu 14-Jan-2000
83			c - restructured the original version in order to have a
84			c better interface to the MITgcmUV.
85			c
86			c Christian Eckert eckert@mit.edu 12-Feb-2000
87			c - Changed Routine names (package prefix: exf_)
88			c
89			c Patrick Heimbach, heimbach@mit.edu 04-May-2000
90			c - changed the handling of precip and sflux with respect
91			c to CPP options ALLOW_BULKFORMULAE and ALLOW_ATM_TEMP
92			c - included some CPP flags ALLOW_BULKFORMULAE to make
93			c sure ALLOW_ATM_TEMP, ALLOW_ATM_WIND are used only in
94			c conjunction with defined ALLOW_BULKFORMULAE
95			c - statement functions discarded
96			c
97			c Ralf.Giering@FastOpt.de 25-Mai-2000
98			c - total rewrite using new subroutines
99			c
100			c Detlef Stammer: include river run-off. Nov. 21, 2001
101			c
102			c heimbach@mit.edu, 10-Jan-2002
103			c - changes to enable field swapping
104			c
105			c mods for pkg/seaice: menemenlis@jpl.nasa.gov 20-Dec-2002
106			c
107			c ==================================================================
108			c SUBROUTINE exf_bulkformulae
109			c ==================================================================
110
111			implicit none
112
113			c == global variables ==
114
115			#include "EEPARAMS.h"
116			#include "SIZE.h"
117			#include "PARAMS.h"
118			#include "DYNVARS.h"
119	jmc	1.13	c #include "GRID.h"
120	heimbach	1.2
121	jmc	1.15	#include "EXF_PARAM.h"
122			#include "EXF_FIELDS.h"
123			#include "EXF_CONSTANTS.h"
124	heimbach	1.2
125			#ifdef ALLOW_AUTODIFF_TAMC
126			#include "tamc.h"
127			#endif
128
129			c == routine arguments ==
130
131			integer mythid
132	heimbach	1.12	integer myiter
133			_RL mytime
134	heimbach	1.2
135	heimbach	1.3	#ifdef ALLOW_BULKFORMULAE
136
137	heimbach	1.2	c == local variables ==
138
139			integer bi,bj
140			integer i,j,k
141
142	heimbach	1.16	_RL aln
143			_RL czol
144	heimbach	1.2
145			integer iter
146	heimbach	1.16	_RL tmpbulk
147	heimbach	1.2	_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			_RL ssttmp
160			_RL ssq
161			_RL stable
162			_RL tstar
163			_RL t0
164			_RL ustar
165			_RL uzn
166			_RL shn
167			_RL xsq
168			_RL x
169			_RL tau
170			#ifdef ALLOW_AUTODIFF_TAMC
171			integer ikey_1
172			integer ikey_2
173			#endif
174
175			c == external functions ==
176
177			integer ilnblnk
178			external ilnblnk
179
180			_RL exf_BulkqSat
181			external exf_BulkqSat
182			_RL exf_BulkCdn
183			external exf_BulkCdn
184			_RL exf_BulkRhn
185			external exf_BulkRhn
186
187			c == end of interface ==
188
189			aln = log(ht/zref)
190	heimbach	1.16	czol = hukarmangravity_mks
191	heimbach	1.2
192			c-- Use atmospheric state to compute surface fluxes.
193
194			c Loop over tiles.
195			#ifdef ALLOW_AUTODIFF_TAMC
196			C-- HPF directive to help TAMC
197			CHPF$ INDEPENDENT
198			#endif
199			do bj = mybylo(mythid),mybyhi(mythid)
200			#ifdef ALLOW_AUTODIFF_TAMC
201			C-- HPF directive to help TAMC
202			CHPF$ INDEPENDENT
203			#endif
204	heimbach	1.12	do bi = mybxlo(mythid),mybxhi(mythid)
205			k = 1
206			do j = 1,sny
207			do i = 1,snx
208	heimbach	1.2
209			#ifdef ALLOW_AUTODIFF_TAMC
210	heimbach	1.12	act1 = bi - myBxLo(myThid)
211			max1 = myBxHi(myThid) - myBxLo(myThid) + 1
212			act2 = bj - myByLo(myThid)
213			max2 = myByHi(myThid) - myByLo(myThid) + 1
214			act3 = myThid - 1
215			max3 = nTx*nTy
216			act4 = ikey_dynamics - 1
217
218			ikey_1 = i
219			& + sNx*(j-1)
220			& + sNxsNyact1
221			& + sNxsNymax1*act2
222			& + sNxsNymax1max2act3
223			& + sNxsNymax1max2max3*act4
224	heimbach	1.2	#endif
225
226			c-- Compute the turbulent surface fluxes.
227
228			#ifdef ALLOW_ATM_TEMP
229
230			c Initial guess: z/l=0.0; hu=ht=hq=z
231			c Iterations: converge on z/l and hence the fluxes.
232			c t0 : virtual temperature (K)
233			c ssq : sea surface humidity (kg/kg)
234			c deltap : potential temperature diff (K)
235
236	heimbach	1.12	if ( atemp(i,j,bi,bj) .ne. 0. _d 0 ) then
237			t0 = atemp(i,j,bi,bj)*
238			& (exf_one + humid_fac*aqh(i,j,bi,bj))
239			ssttmp = theta(i,j,k,bi,bj)
240			tmpbulk= exf_BulkqSat(ssttmp + cen2kel)
241			ssq = saltsat*tmpbulk/atmrho
242			deltap = atemp(i,j,bi,bj) + gamma_blk*ht -
243			& ssttmp - cen2kel
244			delq = aqh(i,j,bi,bj) - ssq
245			stable = exf_half + sign(exf_half, deltap)
246	heimbach	1.2	#ifdef ALLOW_AUTODIFF_TAMC
247	heimbach	1.12	CADJ STORE sh(i,j,bi,bj) = comlev1_exf_1, key = ikey_1
248	heimbach	1.2	#endif
249	jmc	1.13	#ifdef ALLOW_ATM_WIND
250	heimbach	1.12	tmpbulk= exf_BulkCdn(sh(i,j,bi,bj))
251			rdn = sqrt(tmpbulk)
252			ustar = rdn*sh(i,j,bi,bj)
253	jmc	1.13	#else /* ifndef ALLOW_ATM_WIND */
254	mlosch	1.14	C in this case ustress and vstress are defined a u and v points
255			C respectively, and we need to average to mass points to avoid
256			C tau = 0 near coasts or other boundaries
257			tau = sqrt(0.5*
258			& (ustress(i, j,bi,bj)*ustress(i ,j,bi,bj)
259			& +ustress(i+1,j,bi,bj)*ustress(i+1,j,bi,bj)
260			& +vstress(i,j, bi,bj)*vstress(i,j ,bi,bj)
261			& +vstress(i,j+1,bi,bj)*vstress(i,j+1,bi,bj))
262			& )
263	jmc	1.13	ustar = sqrt(tau/atmrho)
264			#endif /* ifndef ALLOW_ATM_WIND */
265	heimbach	1.12	tmpbulk= exf_BulkRhn(stable)
266	jmc	1.13	tstar = tmpbulk*deltap
267			qstar = cdalton*delq
268	heimbach	1.2
269	heimbach	1.12	do iter = 1,niter_bulk
270	heimbach	1.2
271			#ifdef ALLOW_AUTODIFF_TAMC
272			ikey_2 = iter
273			& + niter_bulk*(i-1)
274	heimbach	1.6	& + niter_bulksNx(j-1)
275			& + niter_bulksNxsNy*act1
276			& + niter_bulksNxsNymax1act2
277			& + niter_bulksNxsNymax1max2*act3
278			& + niter_bulksNxsNymax1max2max3act4
279	heimbach	1.2
280			CADJ STORE rdn = comlev1_exf_2, key = ikey_2
281			CADJ STORE ustar = comlev1_exf_2, key = ikey_2
282			CADJ STORE qstar = comlev1_exf_2, key = ikey_2
283			CADJ STORE tstar = comlev1_exf_2, key = ikey_2
284	heimbach	1.12	CADJ STORE sh(i,j,bi,bj) = comlev1_exf_2, key = ikey_2
285			CADJ STORE us(i,j,bi,bj) = comlev1_exf_2, key = ikey_2
286	heimbach	1.2	#endif
287
288			huol = czol*(tstar/t0 +
289			& qstar/(exf_one/humid_fac+aqh(i,j,bi,bj)))/
290			& ustar**2
291			huol = max(huol,zolmin)
292			stable = exf_half + sign(exf_half, huol)
293			htol = huol*ht/hu
294			hqol = huol*hq/hu
295
296			c Evaluate all stability functions assuming hq = ht.
297	jmc	1.13	#ifdef ALLOW_ATM_WIND
298	heimbach	1.2	xsq = max(sqrt(abs(exf_one - 16.*huol)),exf_one)
299			x = sqrt(xsq)
300			psimh = -psim_fachuolstable +
301			& (exf_one - stable)*
302	heimbach	1.9	& (log((exf_one + x(exf_two + x))
303	heimbach	1.2	& (exf_one + xsq)/8.) - exf_two*atan(x) +
304	heimbach	1.9	& pi*exf_half)
305	jmc	1.13	#endif /* ALLOW_ATM_WIND */
306	heimbach	1.2	xsq = max(sqrt(abs(exf_one - 16.*htol)),exf_one)
307			psixh = -psim_fachtolstable + (exf_one - stable)*
308			& exf_two*log((exf_one + xsq)/exf_two)
309
310	jmc	1.13	#ifdef ALLOW_ATM_WIND
311	heimbach	1.2	c Shift wind speed using old coefficient
312	heimbach	1.6	rd = rdn/(exf_one - rdn/karman*psimh )
313	heimbach	1.12	shn = sh(i,j,bi,bj)*rd/rdn
314	heimbach	1.6	uzn = max(shn, umin)
315	heimbach	1.2
316			c Update the transfer coefficients at 10 meters
317			c and neutral stability.
318
319	heimbach	1.6	tmpbulk= exf_BulkCdn(uzn)
320			rdn = sqrt(tmpbulk)
321	heimbach	1.2
322			c Shift all coefficients to the measurement height
323			c and stability.
324	heimbach	1.6	rd = rdn/(exf_one - rdn/karman*psimh)
325	jmc	1.13	#endif /* ALLOW_ATM_WIND */
326	heimbach	1.6	tmpbulk= exf_BulkRhn(stable)
327	jmc	1.13	rh = tmpbulk/( exf_one +
328	heimbach	1.6	& tmpbulk/karman*(aln - psixh) )
329	jmc	1.13	re = cdalton/( exf_one +
330	heimbach	1.6	& cdalton/karman*(aln - psixh) )
331	heimbach	1.2
332			c Update ustar, tstar, qstar using updated, shifted
333			c coefficients.
334	jmc	1.13	qstar = re*delq
335			tstar = rh*deltap
336			#ifdef ALLOW_ATM_WIND
337	heimbach	1.12	ustar = rd*sh(i,j,bi,bj)
338	heimbach	1.6	tau = atmrhoustar*2
339	heimbach	1.12	tau = tau*us(i,j,bi,bj)/sh(i,j,bi,bj)
340	jmc	1.13	#endif /* ALLOW_ATM_WIND */
341	heimbach	1.2
342			enddo
343
344			#ifdef ALLOW_AUTODIFF_TAMC
345			CADJ STORE ustar = comlev1_exf_1, key = ikey_1
346			CADJ STORE qstar = comlev1_exf_1, key = ikey_1
347			CADJ STORE tstar = comlev1_exf_1, key = ikey_1
348			CADJ STORE tau = comlev1_exf_1, key = ikey_1
349	heimbach	1.12	CADJ STORE cw(i,j,bi,bj) = comlev1_exf_1, key = ikey_1
350			CADJ STORE sw(i,j,bi,bj) = comlev1_exf_1, key = ikey_1
351	heimbach	1.2	#endif
352			hs(i,j,bi,bj) = atmcptautstar/ustar
353			hl(i,j,bi,bj) = flambtauqstar/ustar
354			#ifndef EXF_READ_EVAP
355			cdm evap(i,j,bi,bj) = tau*qstar/ustar
356			cdm !!! need to change sign and to convert from kg/m^2/s to m/s !!!
357			evap(i,j,bi,bj) = -recip_rhoniltauqstar/ustar
358			#endif
359	jmc	1.13	#ifdef ALLOW_ATM_WIND
360	heimbach	1.12	ustress(i,j,bi,bj) = tau*cw(i,j,bi,bj)
361			vstress(i,j,bi,bj) = tau*sw(i,j,bi,bj)
362	jmc	1.13	#endif /* ALLOW_ATM_WIND */
363	heimbach	1.2	else
364	jmc	1.13	#ifdef ALLOW_ATM_WIND
365	heimbach	1.2	ustress(i,j,bi,bj) = 0. _d 0
366			vstress(i,j,bi,bj) = 0. _d 0
367	jmc	1.13	#endif /* ALLOW_ATM_WIND */
368	heimbach	1.2	hflux (i,j,bi,bj) = 0. _d 0
369			hs(i,j,bi,bj) = 0. _d 0
370			hl(i,j,bi,bj) = 0. _d 0
371			endif
372
373			#else /* ifndef ALLOW_ATM_TEMP */
374			#ifdef ALLOW_ATM_WIND
375	heimbach	1.12	tmpbulk= exf_BulkCdn(sh(i,j,bi,bj))
376			ustress(i,j,bi,bj) = atmrhotmpbulkus(i,j,bi,bj)*
377	heimbach	1.2	& uwind(i,j,bi,bj)
378	heimbach	1.12	vstress(i,j,bi,bj) = atmrhotmpbulkus(i,j,bi,bj)*
379	heimbach	1.2	& vwind(i,j,bi,bj)
380			#endif
381			#endif /* ifndef ALLOW_ATM_TEMP */
382			enddo
383			enddo
384			enddo
385			enddo
386
387	heimbach	1.3	#endif /* ALLOW_BULKFORMULAE */
388	heimbach	1.2
389	jmc	1.13	return
390	heimbach	1.2	end