pkg/exf/exf_bulkformulae.F

6d54cf9ca1 Ed H*0001 #include "EXF_OPTIONS.h"
3a255f48df Gael*0002 #ifdef ALLOW_AUTODIFF
                0003 # include "AUTODIFF_OPTIONS.h"
                0004 #endif
bdec91d862 Patr*0005
d216b1690c Mart*0006 CBOP
                0007 C     !ROUTINE: EXF_BULKFORMULAE
                0008 C     !INTERFACE:
                0009       SUBROUTINE EXF_BULKFORMULAE( myTime, myIter, myThid )
                0010
                0011 C     !DESCRIPTION: \bv
                0012 C     *==========================================================*
                0013 C     | SUBROUTINE EXF_BULKFORMULAE
                0014 C     | o Calculate bulk formula fluxes over open ocean
                0015 C     |   either following
e1fb02e8f0 Jean*0016 C     |   Large and Pond, 1981 & 1982
d216b1690c Mart*0017 C     |   or (if defined ALLOW_BULK_LARGEYEAGER04)
                0018 C     |   Large and Yeager, 2004, NCAR/TN-460+STR.
                0019 C     *==========================================================*
                0020 C     \ev
                0021 C
                0022 C
                0023 C     NOTES:
                0024 C     ======
                0025 C
                0026 C     See EXF_OPTIONS.h for a description of the various possible
                0027 C     ocean-model forcing configurations.
                0028 C
                0029 C     The bulk formulae of pkg/exf are not valid for sea-ice covered
                0030 C     oceans but they can be used in combination with a sea-ice model,
                0031 C     for example, pkg/seaice, to specify open water flux contributions.
                0032 C
                0033 C     ==================================================================
                0034 C
                0035 C     for Large and Pond, 1981 & 1982
                0036 C
                0037 C     The calculation of the bulk surface fluxes has been adapted from
                0038 C     the NCOM model which uses the formulae given in Large and Pond
                0039 C     (1981 & 1982 )
                0040 C
                0041 C
                0042 C     Header taken from NCOM version: ncom1.4.1
                0043 C     -----------------------------------------
                0044 C
                0045 C     Following procedures and coefficients in Large and Pond
                0046 C     (1981 ; 1982)
                0047 C
                0048 C     Output: Bulk estimates of the turbulent surface fluxes.
                0049 C     -------
                0050 C
                0051 C     hs  - sensible heat flux  (W/m^2), into ocean
                0052 C     hl  - latent   heat flux  (W/m^2), into ocean
                0053 C
                0054 C     Input:
                0055 C     ------
                0056 C
                0057 C     us  - mean wind speed (m/s)     at height hu (m)
                0058 C     th  - mean air temperature (K)  at height ht (m)
                0059 C     qh  - mean air humidity (kg/kg) at height hq (m)
                0060 C     sst - sea surface temperature (K)
                0061 C     tk0 - Kelvin temperature at 0 Celsius (K)
                0062 C
                0063 C     Assume 1) a neutral 10m drag coefficient =
                0064 C
                0065 C               cdn = .0027/u10 + .000142 + .0000764 u10
                0066 C
                0067 C            2) a neutral 10m stanton number =
                0068 C
                0069 C               ctn = .0327 sqrt(cdn), unstable
                0070 C               ctn = .0180 sqrt(cdn), stable
                0071 C
                0072 C            3) a neutral 10m dalton number =
                0073 C
                0074 C               cen = .0346 sqrt(cdn)
                0075 C
                0076 C            4) the saturation humidity of air at
                0077 C
                0078 C               t(k) = exf_BulkqSat(t)  (kg/m^3)
                0079 C
                0080 C     Note:  1) here, tstar = <wt>/u*, and qstar = <wq>/u*.
                0081 C            2) wind speeds should all be above a minimum speed,
                0082 C               say 0.5 m/s.
                0083 C            3) with optional iteration loop, niter=3, should suffice.
                0084 C            4) this version is for analyses inputs with hu = 10m and
                0085 C               ht = hq.
                0086 C            5) sst enters in Celsius.
                0087 C
                0088 C     ==================================================================
                0089 C
                0090 C       started: Christian Eckert eckert@mit.edu 27-Aug-1999
                0091 C
                0092 C     changed: Christian Eckert eckert@mit.edu 14-Jan-2000
                0093 C            - restructured the original version in order to have a
                0094 C              better interface to the MITgcmUV.
                0095 C
                0096 C            Christian Eckert eckert@mit.edu  12-Feb-2000
                0097 C            - Changed Routine names (package prefix: exf_)
                0098 C
                0099 C            Patrick Heimbach, heimbach@mit.edu  04-May-2000
                0100 C            - changed the handling of precip and sflux with respect
                0101 C              to CPP options ALLOW_BULKFORMULAE and ALLOW_ATM_TEMP
                0102 C            - included some CPP flags ALLOW_BULKFORMULAE to make
                0103 C              sure ALLOW_ATM_TEMP, ALLOW_ATM_WIND are used only in
                0104 C              conjunction with defined ALLOW_BULKFORMULAE
                0105 C            - statement functions discarded
                0106 C
                0107 C            Ralf.Giering@FastOpt.de 25-Mai-2000
                0108 C            - total rewrite using new subroutines
                0109 C
                0110 C            Detlef Stammer: include river run-off. Nov. 21, 2001
                0111 C
                0112 C            heimbach@mit.edu, 10-Jan-2002
                0113 C            - changes to enable field swapping
                0114 C
                0115 C     mods for pkg/seaice: menemenlis@jpl.nasa.gov 20-Dec-2002
                0116 C
                0117 C     martin.losch@awi.de: merged with exf_bulk_largeyeager04, 21-May-2010
                0118 C
                0119 C     ==================================================================
                0120 C
                0121 C     for Large and Yeager, 2004
                0122 C
                0123 C === Turbulent Fluxes :
                0124 C  * use the approach "B": shift coeff to height & stability of the
                0125 C    atmosphere state (instead of "C": shift temp & humid to the height
                0126 C    of wind, then shift the coeff to this height & stability of the atmos).
                0127 C  * similar to EXF (except over sea-ice) ; default parameter values
                0128 C    taken from Large & Yeager.
                0129 C  * assume that Qair & Tair inputs are from the same height (zq=zt)
                0130 C  * formulae in short:
                0131 C     wind stress = (ust,vst) = rhoA * Cd * Ws * (del.u,del.v)
                0132 C     Sensib Heat flux = fsha = rhoA * Ch * Ws * del.T * CpAir
                0133 C     Latent Heat flux = flha = rhoA * Ce * Ws * del.Q * Lvap
                0134 C                      = -Evap * Lvap
                0135 C   with Ws = wind speed = sqrt(del.u^2 +del.v^2) ;
                0136 C        del.T = Tair - Tsurf ; del.Q = Qair - Qsurf ;
                0137 C        Cd,Ch,Ce = drag coefficient, Stanton number and Dalton number
                0138 C              respectively [no-units], function of height & stability
                0139
                0140 C     !USES:
                0141        IMPLICIT NONE
                0142 C     === Global variables ===
bdec91d862 Patr*0143 #include "EEPARAMS.h"
                0144 #include "SIZE.h"
                0145 #include "PARAMS.h"
                0146 #include "DYNVARS.h"
d216b1690c Mart*0147 #include "GRID.h"
bdec91d862 Patr*0148
082e18c36c Jean*0149 #include "EXF_PARAM.h"
                0150 #include "EXF_FIELDS.h"
                0151 #include "EXF_CONSTANTS.h"
bdec91d862 Patr*0152
                0153 #ifdef ALLOW_AUTODIFF_TAMC
                0154 #include "tamc.h"
                0155 #endif
                0156
d216b1690c Mart*0157 C     !INPUT/OUTPUT PARAMETERS:
                0158 C     input:
                0159 C     myTime  :: Current time in simulation
                0160 C     myIter  :: Current iteration number in simulation
                0161 C     myThid  :: My Thread Id number
                0162       _RL     myTime
                0163       INTEGER myIter
                0164       INTEGER myThid
                0165 C     output:
                0166 CEOP
bdec91d862 Patr*0167
7109a141b2 Patr*0168 #ifdef ALLOW_BULKFORMULAE
423768d890 Jean*0169 #ifdef ALLOW_ATM_TEMP
d216b1690c Mart*0170 C     == external Functions
7109a141b2 Patr*0171
d216b1690c Mart*0172 C     == Local variables ==
                0173 C     i,j      :: grid point indices
                0174 C     bi,bj    :: tile indices
                0175 C     ssq      :: saturation specific humidity [kg/kg] in eq. with Sea-Surface water
                0176       INTEGER i,j,bi,bj
bdec91d862 Patr*0177
20cf30c902 Patr*0178       _RL czol
d216b1690c Mart*0179       _RL Tsf                   ! surface temperature [K]
                0180       _RL wsm                   ! limited wind speed [m/s] (> umin)
                0181       _RL t0                    ! virtual temperature [K]
                0182 C     these need to be 2D-arrays for vectorizing code
                0183       _RL tstar (1:sNx,1:sNy)   ! turbulent temperature scale [K]
                0184       _RL qstar (1:sNx,1:sNy)   ! turbulent humidity scale  [kg/kg]
                0185       _RL ustar (1:sNx,1:sNy)   ! friction velocity [m/s]
                0186       _RL tau   (1:sNx,1:sNy)   ! surface stress coef = rhoA * Ws * sqrt(Cd)
                0187       _RL rdn   (1:sNx,1:sNy)   ! neutral, zref (=10m) values of rd
                0188       _RL rd    (1:sNx,1:sNy)   ! = sqrt(Cd)          [-]
                0189       _RL delq  (1:sNx,1:sNy)   ! specific humidity difference [kg/kg]
0b902863ee Mart*0190       _RL deltap(1:sNx,1:sNy)
8053a926ae Jean*0191 #ifdef EXF_CALC_ATMRHO
                0192       _RL atmrho_loc(1:sNx,1:sNy) ! local atmospheric density [kg/m^3]
                0193 #endif
                0194
e1fb02e8f0 Jean*0195 #ifdef ALLOW_BULK_LARGEYEAGER04
d216b1690c Mart*0196       _RL dzTmp
e1fb02e8f0 Jean*0197 #endif
d216b1690c Mart*0198       _RL SSTtmp
                0199       _RL ssq
                0200       _RL re                    ! = Ce / sqrt(Cd)     [-]
                0201       _RL rh                    ! = Ch / sqrt(Cd)     [-]
                0202       _RL ren, rhn              ! neutral, zref (=10m) values of re, rh
                0203       _RL usn, usm
                0204       _RL stable                ! = 1 if stable ; = 0 if unstable
                0205       _RL huol                  ! stability parameter at zwd [-] (=z/Monin-Obuklov length)
                0206       _RL htol                  ! stability parameter at zth [-]
                0207       _RL hqol
                0208       _RL x                     ! stability function  [-]
                0209       _RL xsq                   ! = x^2               [-]
                0210       _RL psimh                 ! momentum stability function
                0211       _RL psixh                 ! latent & sensib. stability function
                0212       _RL zwln                  ! = log(zwd/zref)
                0213       _RL ztln                  ! = log(zth/zref)
                0214       _RL tmpbulk
                0215       _RL recip_rhoConstFresh
0320e25227 Mart*0216       INTEGER ks, kl
d216b1690c Mart*0217       INTEGER iter
                0218 C     solve4Stress :: if F, by-pass momentum turb. flux (wind-stress) calculation
                0219       LOGICAL solve4Stress
e1fb02e8f0 Jean*0220       _RL windStress            ! surface wind-stress (@ grid cell center)
d216b1690c Mart*0221
bdec91d862 Patr*0222 #ifdef ALLOW_AUTODIFF_TAMC
20dee61641 Mart*0223       INTEGER ikey, ikey_1, ikey_2
bdec91d862 Patr*0224 #endif
                0225
e1fb02e8f0 Jean*0226 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0227
423768d890 Jean*0228       IF ( useAtmWind ) THEN
358649780a Gael*0229         solve4Stress = .TRUE.
423768d890 Jean*0230       ELSE
d216b1690c Mart*0231 #ifdef ALLOW_BULK_LARGEYEAGER04
358649780a Gael*0232         solve4Stress = wspeedfile .NE. ' '
d216b1690c Mart*0233 #else
358649780a Gael*0234         solve4Stress = .FALSE.
d216b1690c Mart*0235 #endif
423768d890 Jean*0236       ENDIF
bdec91d862 Patr*0237
0320e25227 Mart*0238       IF ( usingPCoords ) THEN
                0239        ks = Nr
                0240        kl = Nr-1
                0241       ELSE
                0242        ks = 1
                0243        kl = 2
                0244       ENDIF
                0245
                0246 C--   Set surface parameters :
666d41f610 Mart*0247       zwln = LOG(hu/zref)
                0248       ztln = LOG(ht/zref)
20cf30c902 Patr*0249       czol = hu*karman*gravity_mks
d216b1690c Mart*0250 C--   abbreviation
666d41f610 Mart*0251       recip_rhoConstFresh = 1. _d 0/rhoConstFresh
e1fb02e8f0 Jean*0252
8053a926ae Jean*0253 C     Loop over tiles.
bdec91d862 Patr*0254 #ifdef ALLOW_AUTODIFF_TAMC
                0255 C--   HPF directive to help TAMC
                0256 CHPF$ INDEPENDENT
                0257 #endif
d216b1690c Mart*0258       DO bj = myByLo(myThid),myByHi(myThid)
bdec91d862 Patr*0259 #ifdef ALLOW_AUTODIFF_TAMC
                0260 C--    HPF directive to help TAMC
                0261 CHPF$  INDEPENDENT
                0262 #endif
d216b1690c Mart*0263        DO bi = myBxLo(myThid),myBxHi(myThid)
0320e25227 Mart*0264
20dee61641 Mart*0265 #ifdef ALLOW_AUTODIFF_TAMC
                0266         ikey = bi + (bj-1)*nSx + (ikey_dynamics-1)*nSx*nSy
                0267 #endif
d216b1690c Mart*0268         DO j = 1,sNy
                0269          DO i = 1,sNx
bdec91d862 Patr*0270
                0271 #ifdef ALLOW_AUTODIFF_TAMC
20dee61641 Mart*0272           ikey_1 = i + (j-1)*sNx + (ikey-1)*sNx*sNy
bdec91d862 Patr*0273 #endif
                0274
423768d890 Jean*0275 #ifdef ALLOW_AUTODIFF
f1644eaf0d Patr*0276           deltap(i,j) = 0. _d 0
                0277           delq(i,j)   = 0. _d 0
                0278 #endif
d216b1690c Mart*0279 C--- Compute turbulent surface fluxes
                0280 C-   Pot. Temp and saturated specific humidity
                0281           IF ( atemp(i,j,bi,bj) .NE. 0. _d 0 ) THEN
                0282 C-   Surface Temp.
0320e25227 Mart*0283            Tsf = theta(i,j,ks,bi,bj) + cen2kel
c2c9eed210 Jean*0284            IF ( Nr.GE.2 .AND. sstExtrapol.GT.0. _d 0 ) THEN
d216b1690c Mart*0285             SSTtmp = sstExtrapol
0320e25227 Mart*0286      &              *( theta(i,j,ks,bi,bj)-theta(i,j,kl,bi,bj) )
                0287      &              *  maskC(i,j,kl,bi,bj)
d216b1690c Mart*0288             Tsf = Tsf + MAX( SSTtmp, 0. _d 0 )
                0289            ENDIF
8053a926ae Jean*0290
d216b1690c Mart*0291            tmpbulk = cvapor_fac*exp(-cvapor_exp/Tsf)
8053a926ae Jean*0292 #ifdef EXF_CALC_ATMRHO
                0293            atmrho_loc(i,j) = apressure(i,j,bi,bj) /
                0294      &                  (287.04 _d 0*atemp(i,j,bi,bj)
                0295      &                  *(1. _d 0 + humid_fac*aqh(i,j,bi,bj)))
                0296            ssq = saltsat*tmpbulk/atmrho_loc(i,j)
                0297 #else
d216b1690c Mart*0298            ssq = saltsat*tmpbulk/atmrho
8053a926ae Jean*0299 #endif
d216b1690c Mart*0300            deltap(i,j) = atemp(i,j,bi,bj) + gamma_blk*ht - Tsf
0b902863ee Mart*0301            delq(i,j)   = aqh(i,j,bi,bj) - ssq
d216b1690c Mart*0302 C--  no negative evaporation over ocean:
                0303            IF ( noNegativeEvap ) delq(i,j) = MIN( 0. _d 0, delq(i,j) )
                0304
                0305 C--  initial guess for exchange coefficients:
                0306 C    take U_N = del.U ; stability from del.Theta ;
                0307            stable = exf_half + SIGN(exf_half, deltap(i,j))
                0308
                0309            IF ( solve4Stress ) THEN
bdec91d862 Patr*0310 #ifdef ALLOW_AUTODIFF_TAMC
3752238fd8 Patr*0311 CADJ STORE sh(i,j,bi,bj)     = comlev1_exf_1, key = ikey_1
d216b1690c Mart*0312 #endif /* ALLOW_AUTODIFF_TAMC */
                0313 C--   Wind speed
                0314             wsm        = sh(i,j,bi,bj)
cc60455fbb Mart*0315 #ifdef  ALLOW_DRAG_LARGEYEAGER09
                0316 C     Large and Yeager (2009), Climate Dynamics, equation 11a/b
                0317             tmpbulk =  cdrag_1/wsm + cdrag_2 + cdrag_3*wsm
                0318      &           + cdrag_8 * wsm**6
                0319             tmpbulk = exf_scal_BulkCdn * (
                0320      &             ( halfRL - SIGN(halfRL, wsm-umax) )*tmpbulk
                0321      &           + ( halfRL + SIGN(halfRL, wsm-umax) )*cdragMax
                0322      &           )
                0323 #else
d216b1690c Mart*0324             tmpbulk    = exf_scal_BulkCdn *
                0325      &           ( cdrag_1/wsm + cdrag_2 + cdrag_3*wsm )
cc60455fbb Mart*0326 #endif
d216b1690c Mart*0327             rdn(i,j)   = SQRT(tmpbulk)
                0328             ustar(i,j) = rdn(i,j)*wsm
                0329            ELSE
ddea17ebdb Gael*0330             rdn(i,j)   = 0. _d 0
d216b1690c Mart*0331             windStress = wStress(i,j,bi,bj)
8053a926ae Jean*0332 #ifdef EXF_CALC_ATMRHO
423768d890 Jean*0333             ustar(i,j) = SQRT(windStress/atmrho_loc(i,j))
                0334             tau(i,j)   = SQRT(windStress*atmrho_loc(i,j))
8053a926ae Jean*0335 #else
423768d890 Jean*0336             ustar(i,j) = SQRT(windStress/atmrho)
d216b1690c Mart*0337 c           tau(i,j)   = windStress/ustar(i,j)
423768d890 Jean*0338             tau(i,j)   = SQRT(windStress*atmrho)
8053a926ae Jean*0339 #endif
d216b1690c Mart*0340            ENDIF
7311a000e1 Mart*0341
                0342 C--  initial guess for exchange other coefficients:
666d41f610 Mart*0343            rhn = (exf_one-stable)*cstanton_1 + stable*cstanton_2
d216b1690c Mart*0344            ren = cDalton
7311a000e1 Mart*0345 C--  calculate turbulent scales
d216b1690c Mart*0346            tstar(i,j)=rhn*deltap(i,j)
                0347            qstar(i,j)=ren*delq(i,j)
0b902863ee Mart*0348
d216b1690c Mart*0349           ELSE
e1fb02e8f0 Jean*0350 C     atemp(i,j,bi,bj) .EQ. 0.
0b902863ee Mart*0351            tstar (i,j) = 0. _d 0
                0352            qstar (i,j) = 0. _d 0
                0353            ustar (i,j) = 0. _d 0
                0354            tau   (i,j) = 0. _d 0
                0355            rdn   (i,j) = 0. _d 0
d216b1690c Mart*0356           ENDIF
                0357          ENDDO
                0358         ENDDO
                0359         DO iter=1,niter_bulk
                0360          DO j = 1,sNy
                0361           DO i = 1,sNx
                0362            IF ( atemp(i,j,bi,bj) .NE. 0. _d 0 ) THEN
                0363 C--- iterate with psi-functions to find transfer coefficients
                0364
bdec91d862 Patr*0365 #ifdef ALLOW_AUTODIFF_TAMC
20dee61641 Mart*0366             ikey_2 = i + (j-1)*sNx + (iter-1)*sNx*sNy
                0367      &           + (ikey-1)*sNx*sNy*niter_bulk
d216b1690c Mart*0368 CADJ STORE rdn   (i,j)       = comlev1_exf_2, key = ikey_2
                0369 CADJ STORE ustar (i,j)       = comlev1_exf_2, key = ikey_2
                0370 CADJ STORE qstar (i,j)       = comlev1_exf_2, key = ikey_2
                0371 CADJ STORE tstar (i,j)       = comlev1_exf_2, key = ikey_2
                0372 CADJ STORE sh    (i,j,bi,bj) = comlev1_exf_2, key = ikey_2
                0373 CADJ STORE wspeed(i,j,bi,bj) = comlev1_exf_2, key = ikey_2
bdec91d862 Patr*0374 #endif
                0375
d216b1690c Mart*0376             t0   = atemp(i,j,bi,bj)*
0b902863ee Mart*0377      &           (exf_one + humid_fac*aqh(i,j,bi,bj))
d216b1690c Mart*0378             huol = ( tstar(i,j)/t0 +
                0379      &               qstar(i,j)/(exf_one/humid_fac+aqh(i,j,bi,bj))
                0380      &              )*czol/(ustar(i,j)*ustar(i,j))
                0381 #ifdef ALLOW_BULK_LARGEYEAGER04
                0382             tmpbulk = MIN(abs(huol),10. _d 0)
                0383             huol   = SIGN(tmpbulk , huol)
                0384 #else
                0385 C--   Large&Pond1981:
0b902863ee Mart*0386             huol   = max(huol,zolmin)
d216b1690c Mart*0387 #endif /* ALLOW_BULK_LARGEYEAGER04 */
0b902863ee Mart*0388             htol   = huol*ht/hu
                0389             hqol   = huol*hq/hu
666d41f610 Mart*0390             stable = exf_half + sign(exf_half, huol)
d216b1690c Mart*0391
8053a926ae Jean*0392 C     Evaluate all stability functions assuming hq = ht.
d216b1690c Mart*0393             IF ( solve4Stress ) THEN
                0394 #ifdef ALLOW_BULK_LARGEYEAGER04
                0395 C--   Large&Yeager04:
                0396              xsq    = SQRT( ABS(exf_one - huol*16. _d 0) )
                0397 #else
                0398 C--   Large&Pond1981:
423768d890 Jean*0399              xsq    = MAX(SQRT(ABS(exf_one - 16.*huol)),exf_one)
d216b1690c Mart*0400 #endif /* ALLOW_BULK_LARGEYEAGER04 */
                0401              x      = SQRT(xsq)
                0402              psimh  = -psim_fac*huol*stable
                0403      &               +(exf_one-stable)*
                0404      &                ( LOG( (exf_one + exf_two*x + xsq)*
                0405      &                       (exf_one+xsq)*.125 _d 0 )
                0406      &                  -exf_two*ATAN(x) + exf_half*pi )
ddea17ebdb Gael*0407             ELSE
                0408              psimh  = 0. _d 0
d216b1690c Mart*0409             ENDIF
                0410 #ifdef ALLOW_BULK_LARGEYEAGER04
                0411 C--   Large&Yeager04:
                0412             xsq   = SQRT( ABS(exf_one - htol*16. _d 0) )
                0413 #else
                0414 C--   Large&Pond1981:
423768d890 Jean*0415             xsq   = MAX(SQRT(ABS(exf_one - 16.*htol)),exf_one)
d216b1690c Mart*0416 #endif /* ALLOW_BULK_LARGEYEAGER04 */
666d41f610 Mart*0417             psixh = -psim_fac*htol*stable + (exf_one-stable)*
                0418      &               ( exf_two*LOG(exf_half*(exf_one+xsq)) )
d216b1690c Mart*0419
                0420             IF ( solve4Stress ) THEN
7448700841 Mart*0421 #ifdef ALLOW_AUTODIFF_TAMC
8ed910e28e Patr*0422 CADJ STORE rdn(i,j)       = comlev1_exf_2, key = ikey_2
7448700841 Mart*0423 #endif
d216b1690c Mart*0424 C-   Shift wind speed using old coefficient
                0425 #ifdef ALLOW_BULK_LARGEYEAGER04
                0426 C--   Large&Yeager04:
                0427              dzTmp = (zwln-psimh)/karman
                0428              usn   = wspeed(i,j,bi,bj)/(exf_one + rdn(i,j)*dzTmp )
                0429 #else
                0430 C--   Large&Pond1981:
666d41f610 Mart*0431 c           rd(i,j)= rdn(i,j)/(exf_one - rdn(i,j)/karman*psimh )
d216b1690c Mart*0432 c           usn    = sh(i,j,bi,bj)*rd(i,j)/rdn(i,j)
8053a926ae Jean*0433 C     ML: the original formulation above is replaced below to be
                0434 C     similar to largeyeager04, but does not give the same results, strange
d216b1690c Mart*0435              usn   = sh(i,j,bi,bj)/(exf_one - rdn(i,j)/karman*psimh)
                0436 #endif /* ALLOW_BULK_LARGEYEAGER04 */
                0437              usm   = MAX(usn, umin)
                0438
                0439 C-   Update the 10m, neutral stability transfer coefficients (momentum)
cc60455fbb Mart*0440 #ifdef  ALLOW_DRAG_LARGEYEAGER09
                0441 C     Large and Yeager (2009), J.Clim equation 11a/b
                0442              tmpbulk =  cdrag_1/usm + cdrag_2 + cdrag_3*usm
                0443      &            + cdrag_8 * usm**6
                0444              tmpbulk = exf_scal_BulkCdn * (
                0445      &              ( halfRL - SIGN(halfRL, usm-umax) )*tmpbulk
                0446      &            + ( halfRL + SIGN(halfRL, usm-umax) )*cdragMax
                0447      &            )
                0448 #else
d216b1690c Mart*0449              tmpbulk  = exf_scal_BulkCdn *
                0450      &            ( cdrag_1/usm + cdrag_2 + cdrag_3*usm )
cc60455fbb Mart*0451 #endif
d216b1690c Mart*0452              rdn(i,j)   = SQRT(tmpbulk)
                0453 #ifdef ALLOW_BULK_LARGEYEAGER04
                0454              rd(i,j)    = rdn(i,j)/(exf_one + rdn(i,j)*dzTmp)
                0455 #else
                0456              rd(i,j)    = rdn(i,j)/(exf_one - rdn(i,j)/karman*psimh)
                0457 #endif /* ALLOW_BULK_LARGEYEAGER04 */
                0458              ustar(i,j) = rd(i,j)*sh(i,j,bi,bj)
                0459 C-   Coeff:
8053a926ae Jean*0460 #ifdef EXF_CALC_ATMRHO
                0461              tau(i,j)   = atmrho_loc(i,j)*rd(i,j)*wspeed(i,j,bi,bj)
                0462 #else
d216b1690c Mart*0463              tau(i,j)   = atmrho*rd(i,j)*wspeed(i,j,bi,bj)
8053a926ae Jean*0464 #endif
d216b1690c Mart*0465             ENDIF
                0466
                0467 C-   Update the 10m, neutral stability transfer coefficients (sens&evap)
666d41f610 Mart*0468             rhn = (exf_one-stable)*cstanton_1 + stable*cstanton_2
d216b1690c Mart*0469             ren = cDalton
                0470
                0471 C-   Shift all coefficients to the measurement height and stability.
                0472             rh = rhn/(exf_one + rhn*(ztln-psixh)/karman)
                0473             re = ren/(exf_one + ren*(ztln-psixh)/karman)
bdec91d862 Patr*0474
d216b1690c Mart*0475 C--  Update ustar, tstar, qstar using updated, shifted coefficients.
0b902863ee Mart*0476             qstar(i,j) = re*delq(i,j)
                0477             tstar(i,j) = rh*deltap(i,j)
f1644eaf0d Patr*0478
d216b1690c Mart*0479            ENDIF
                0480           ENDDO
                0481          ENDDO
8053a926ae Jean*0482 C end of iteration loop
d216b1690c Mart*0483         ENDDO
                0484         DO j = 1,sNy
                0485          DO i = 1,sNx
                0486           IF ( atemp(i,j,bi,bj) .NE. 0. _d 0 ) THEN
bdec91d862 Patr*0487
                0488 #ifdef ALLOW_AUTODIFF_TAMC
20dee61641 Mart*0489            ikey_1 = i + (j-1)*sNx + (ikey-1)*sNx*sNy
7311a000e1 Mart*0490 CADJ STORE qstar(i,j)    = comlev1_exf_1, key = ikey_1
                0491 CADJ STORE tstar(i,j)    = comlev1_exf_1, key = ikey_1
d216b1690c Mart*0492 CADJ STORE tau  (i,j)    = comlev1_exf_1, key = ikey_1
                0493 CADJ STORE rd   (i,j)    = comlev1_exf_1, key = ikey_1
                0494 #endif /* ALLOW_AUTODIFF_TAMC */
                0495
                0496 C-   Turbulent Fluxes
666d41f610 Mart*0497            hs(i,j,bi,bj) = atmcp*tau(i,j)*tstar(i,j)
                0498            hl(i,j,bi,bj) = flamb*tau(i,j)*qstar(i,j)
bdec91d862 Patr*0499 #ifndef EXF_READ_EVAP
d216b1690c Mart*0500 C   change sign and convert from kg/m^2/s to m/s via rhoConstFresh
                0501 c          evap(i,j,bi,bj) = -recip_rhonil*tau(i,j)*qstar(i,j)
666d41f610 Mart*0502            evap(i,j,bi,bj) = -recip_rhoConstFresh*tau(i,j)*qstar(i,j)
d216b1690c Mart*0503 C   but older version was using rhonil instead:
                0504 c           evap(i,j,bi,bj) = -recip_rhonil*tau(i,j)*qstar(i,j)
bdec91d862 Patr*0505 #endif
56d13a40ed Mart*0506            IF ( useAtmWind .AND. useRelativeWind ) THEN
                0507             tmpbulk =  tau(i,j)*rd(i,j)
                0508             ustress(i,j,bi,bj) = tmpbulk* ( uwind(i,j,bi,bj) -
                0509      &           0.5 _d 0 * (uVel(i,j,ks,bi,bj)+uVel(i+1,j,ks,bi,bj)) )
                0510             vstress(i,j,bi,bj) = tmpbulk* ( vwind(i,j,bi,bj) -
                0511      &           0.5 _d 0 * (vVel(i,j,ks,bi,bj)+vVel(i,j+1,ks,bi,bj)) )
                0512            ELSEIF ( useAtmWind ) THEN
423768d890 Jean*0513 c           ustress(i,j,bi,bj) = tau(i,j)*rd(i,j)*ws*cw(i,j,bi,bj)
                0514 c           vstress(i,j,bi,bj) = tau(i,j)*rd(i,j)*ws*sw(i,j,bi,bj)
d216b1690c Mart*0515 C- jmc: below is how it should be written ; different from above when
                0516 C       both wind-speed & 2 compon. of the wind are specified, and in
                0517 C-      this case, formula below is better.
423768d890 Jean*0518             tmpbulk =  tau(i,j)*rd(i,j)
                0519             ustress(i,j,bi,bj) = tmpbulk*uwind(i,j,bi,bj)
                0520             vstress(i,j,bi,bj) = tmpbulk*vwind(i,j,bi,bj)
                0521            ENDIF
7311a000e1 Mart*0522
                0523 c IF ( ATEMP(i,j,bi,bj) .NE. 0. )
423768d890 Jean*0524           ELSE
                0525            IF ( useAtmWind ) ustress(i,j,bi,bj) = 0. _d 0
                0526            IF ( useAtmWind ) vstress(i,j,bi,bj) = 0. _d 0
0b902863ee Mart*0527            hflux  (i,j,bi,bj) = 0. _d 0
d216b1690c Mart*0528            evap   (i,j,bi,bj) = 0. _d 0
                0529            hs     (i,j,bi,bj) = 0. _d 0
                0530            hl     (i,j,bi,bj) = 0. _d 0
7311a000e1 Mart*0531
                0532 c IF ( ATEMP(i,j,bi,bj) .NE. 0. )
423768d890 Jean*0533           ENDIF
                0534
d216b1690c Mart*0535          ENDDO
                0536         ENDDO
                0537        ENDDO
                0538       ENDDO
bdec91d862 Patr*0539
423768d890 Jean*0540 #endif /* ALLOW_ATM_TEMP */
7109a141b2 Patr*0541 #endif /* ALLOW_BULKFORMULAE */
                0542
d216b1690c Mart*0543       RETURN
                0544       END