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:
701e10a905 Mart*0009       SUBROUTINE EXF_BULKFORMULAE( exf_Tsf, myTime, myIter, myThid )
d216b1690c Mart*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
701e10a905 Mart*0162 C     exf_Tsf :: local copy of global field gcmSST or extrapolated
                0163 C                surface temperature (in deg Celsius)
d216b1690c Mart*0164       _RL     myTime
                0165       INTEGER myIter
                0166       INTEGER myThid
701e10a905 Mart*0167       _RL exf_Tsf(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
d216b1690c Mart*0168 C     output:
                0169 CEOP
bdec91d862 Patr*0170
7109a141b2 Patr*0171 #ifdef ALLOW_BULKFORMULAE
423768d890 Jean*0172 #ifdef ALLOW_ATM_TEMP
d216b1690c Mart*0173 C     == external Functions
7109a141b2 Patr*0174
d216b1690c Mart*0175 C     == Local variables ==
                0176 C     i,j      :: grid point indices
                0177 C     bi,bj    :: tile indices
                0178 C     ssq      :: saturation specific humidity [kg/kg] in eq. with Sea-Surface water
                0179       INTEGER i,j,bi,bj
bdec91d862 Patr*0180
20cf30c902 Patr*0181       _RL czol
d216b1690c Mart*0182       _RL Tsf                   ! surface temperature [K]
                0183       _RL wsm                   ! limited wind speed [m/s] (> umin)
                0184       _RL t0                    ! virtual temperature [K]
                0185 C     these need to be 2D-arrays for vectorizing code
                0186       _RL tstar (1:sNx,1:sNy)   ! turbulent temperature scale [K]
                0187       _RL qstar (1:sNx,1:sNy)   ! turbulent humidity scale  [kg/kg]
                0188       _RL ustar (1:sNx,1:sNy)   ! friction velocity [m/s]
                0189       _RL tau   (1:sNx,1:sNy)   ! surface stress coef = rhoA * Ws * sqrt(Cd)
                0190       _RL rdn   (1:sNx,1:sNy)   ! neutral, zref (=10m) values of rd
                0191       _RL rd    (1:sNx,1:sNy)   ! = sqrt(Cd)          [-]
                0192       _RL delq  (1:sNx,1:sNy)   ! specific humidity difference [kg/kg]
0b902863ee Mart*0193       _RL deltap(1:sNx,1:sNy)
8053a926ae Jean*0194 #ifdef EXF_CALC_ATMRHO
                0195       _RL atmrho_loc(1:sNx,1:sNy) ! local atmospheric density [kg/m^3]
                0196 #endif
                0197
e1fb02e8f0 Jean*0198 #ifdef ALLOW_BULK_LARGEYEAGER04
d216b1690c Mart*0199       _RL dzTmp
e1fb02e8f0 Jean*0200 #endif
d216b1690c Mart*0201       _RL ssq
                0202       _RL re                    ! = Ce / sqrt(Cd)     [-]
                0203       _RL rh                    ! = Ch / sqrt(Cd)     [-]
                0204       _RL ren, rhn              ! neutral, zref (=10m) values of re, rh
                0205       _RL usn, usm
                0206       _RL stable                ! = 1 if stable ; = 0 if unstable
                0207       _RL huol                  ! stability parameter at zwd [-] (=z/Monin-Obuklov length)
                0208       _RL htol                  ! stability parameter at zth [-]
                0209       _RL hqol
                0210       _RL x                     ! stability function  [-]
                0211       _RL xsq                   ! = x^2               [-]
                0212       _RL psimh                 ! momentum stability function
                0213       _RL psixh                 ! latent & sensib. stability function
                0214       _RL zwln                  ! = log(zwd/zref)
                0215       _RL ztln                  ! = log(zth/zref)
                0216       _RL tmpbulk
                0217       _RL recip_rhoConstFresh
0320e25227 Mart*0218       INTEGER ks, kl
d216b1690c Mart*0219       INTEGER iter
                0220 C     solve4Stress :: if F, by-pass momentum turb. flux (wind-stress) calculation
                0221       LOGICAL solve4Stress
e1fb02e8f0 Jean*0222       _RL windStress            ! surface wind-stress (@ grid cell center)
d216b1690c Mart*0223
bdec91d862 Patr*0224 #ifdef ALLOW_AUTODIFF_TAMC
20dee61641 Mart*0225       INTEGER ikey, ikey_1, ikey_2
bdec91d862 Patr*0226 #endif
                0227
e1fb02e8f0 Jean*0228 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0229
423768d890 Jean*0230       IF ( useAtmWind ) THEN
358649780a Gael*0231         solve4Stress = .TRUE.
423768d890 Jean*0232       ELSE
d216b1690c Mart*0233 #ifdef ALLOW_BULK_LARGEYEAGER04
358649780a Gael*0234         solve4Stress = wspeedfile .NE. ' '
d216b1690c Mart*0235 #else
358649780a Gael*0236         solve4Stress = .FALSE.
d216b1690c Mart*0237 #endif
423768d890 Jean*0238       ENDIF
bdec91d862 Patr*0239
0320e25227 Mart*0240       IF ( usingPCoords ) THEN
                0241        ks = Nr
                0242        kl = Nr-1
                0243       ELSE
                0244        ks = 1
                0245        kl = 2
                0246       ENDIF
                0247
                0248 C--   Set surface parameters :
666d41f610 Mart*0249       zwln = LOG(hu/zref)
                0250       ztln = LOG(ht/zref)
20cf30c902 Patr*0251       czol = hu*karman*gravity_mks
d216b1690c Mart*0252 C--   abbreviation
666d41f610 Mart*0253       recip_rhoConstFresh = 1. _d 0/rhoConstFresh
e1fb02e8f0 Jean*0254
8053a926ae Jean*0255 C     Loop over tiles.
bdec91d862 Patr*0256 #ifdef ALLOW_AUTODIFF_TAMC
b7b61e618a Mart*0257 C--   HPF directive to help TAF
bdec91d862 Patr*0258 CHPF$ INDEPENDENT
                0259 #endif
d216b1690c Mart*0260       DO bj = myByLo(myThid),myByHi(myThid)
bdec91d862 Patr*0261 #ifdef ALLOW_AUTODIFF_TAMC
b7b61e618a Mart*0262 C--    HPF directive to help TAF
bdec91d862 Patr*0263 CHPF$  INDEPENDENT
                0264 #endif
d216b1690c Mart*0265        DO bi = myBxLo(myThid),myBxHi(myThid)
0320e25227 Mart*0266
20dee61641 Mart*0267 #ifdef ALLOW_AUTODIFF_TAMC
                0268         ikey = bi + (bj-1)*nSx + (ikey_dynamics-1)*nSx*nSy
                0269 #endif
d216b1690c Mart*0270         DO j = 1,sNy
                0271          DO i = 1,sNx
bdec91d862 Patr*0272
                0273 #ifdef ALLOW_AUTODIFF_TAMC
20dee61641 Mart*0274           ikey_1 = i + (j-1)*sNx + (ikey-1)*sNx*sNy
bdec91d862 Patr*0275 #endif
                0276
423768d890 Jean*0277 #ifdef ALLOW_AUTODIFF
f1644eaf0d Patr*0278           deltap(i,j) = 0. _d 0
                0279           delq(i,j)   = 0. _d 0
                0280 #endif
d216b1690c Mart*0281 C--- Compute turbulent surface fluxes
                0282 C-   Pot. Temp and saturated specific humidity
                0283           IF ( atemp(i,j,bi,bj) .NE. 0. _d 0 ) THEN
                0284 C-   Surface Temp.
701e10a905 Mart*0285            Tsf = exf_Tsf(i,j,bi,bj)
d216b1690c Mart*0286            tmpbulk = cvapor_fac*exp(-cvapor_exp/Tsf)
8053a926ae Jean*0287 #ifdef EXF_CALC_ATMRHO
                0288            atmrho_loc(i,j) = apressure(i,j,bi,bj) /
                0289      &                  (287.04 _d 0*atemp(i,j,bi,bj)
                0290      &                  *(1. _d 0 + humid_fac*aqh(i,j,bi,bj)))
                0291            ssq = saltsat*tmpbulk/atmrho_loc(i,j)
                0292 #else
d216b1690c Mart*0293            ssq = saltsat*tmpbulk/atmrho
8053a926ae Jean*0294 #endif
d216b1690c Mart*0295            deltap(i,j) = atemp(i,j,bi,bj) + gamma_blk*ht - Tsf
0b902863ee Mart*0296            delq(i,j)   = aqh(i,j,bi,bj) - ssq
d216b1690c Mart*0297 C--  no negative evaporation over ocean:
                0298            IF ( noNegativeEvap ) delq(i,j) = MIN( 0. _d 0, delq(i,j) )
                0299
                0300 C--  initial guess for exchange coefficients:
                0301 C    take U_N = del.U ; stability from del.Theta ;
                0302            stable = exf_half + SIGN(exf_half, deltap(i,j))
                0303
                0304            IF ( solve4Stress ) THEN
bdec91d862 Patr*0305 #ifdef ALLOW_AUTODIFF_TAMC
3752238fd8 Patr*0306 CADJ STORE sh(i,j,bi,bj)     = comlev1_exf_1, key = ikey_1
d216b1690c Mart*0307 #endif /* ALLOW_AUTODIFF_TAMC */
                0308 C--   Wind speed
                0309             wsm        = sh(i,j,bi,bj)
cc60455fbb Mart*0310 #ifdef  ALLOW_DRAG_LARGEYEAGER09
                0311 C     Large and Yeager (2009), Climate Dynamics, equation 11a/b
                0312             tmpbulk =  cdrag_1/wsm + cdrag_2 + cdrag_3*wsm
                0313      &           + cdrag_8 * wsm**6
                0314             tmpbulk = exf_scal_BulkCdn * (
                0315      &             ( halfRL - SIGN(halfRL, wsm-umax) )*tmpbulk
                0316      &           + ( halfRL + SIGN(halfRL, wsm-umax) )*cdragMax
                0317      &           )
                0318 #else
d216b1690c Mart*0319             tmpbulk    = exf_scal_BulkCdn *
                0320      &           ( cdrag_1/wsm + cdrag_2 + cdrag_3*wsm )
cc60455fbb Mart*0321 #endif
d216b1690c Mart*0322             rdn(i,j)   = SQRT(tmpbulk)
                0323             ustar(i,j) = rdn(i,j)*wsm
                0324            ELSE
ddea17ebdb Gael*0325             rdn(i,j)   = 0. _d 0
d216b1690c Mart*0326             windStress = wStress(i,j,bi,bj)
8053a926ae Jean*0327 #ifdef EXF_CALC_ATMRHO
423768d890 Jean*0328             ustar(i,j) = SQRT(windStress/atmrho_loc(i,j))
                0329             tau(i,j)   = SQRT(windStress*atmrho_loc(i,j))
8053a926ae Jean*0330 #else
423768d890 Jean*0331             ustar(i,j) = SQRT(windStress/atmrho)
d216b1690c Mart*0332 c           tau(i,j)   = windStress/ustar(i,j)
423768d890 Jean*0333             tau(i,j)   = SQRT(windStress*atmrho)
8053a926ae Jean*0334 #endif
d216b1690c Mart*0335            ENDIF
7311a000e1 Mart*0336
                0337 C--  initial guess for exchange other coefficients:
666d41f610 Mart*0338            rhn = (exf_one-stable)*cstanton_1 + stable*cstanton_2
d216b1690c Mart*0339            ren = cDalton
7311a000e1 Mart*0340 C--  calculate turbulent scales
d216b1690c Mart*0341            tstar(i,j)=rhn*deltap(i,j)
                0342            qstar(i,j)=ren*delq(i,j)
0b902863ee Mart*0343
d216b1690c Mart*0344           ELSE
e1fb02e8f0 Jean*0345 C     atemp(i,j,bi,bj) .EQ. 0.
0b902863ee Mart*0346            tstar (i,j) = 0. _d 0
                0347            qstar (i,j) = 0. _d 0
                0348            ustar (i,j) = 0. _d 0
                0349            tau   (i,j) = 0. _d 0
                0350            rdn   (i,j) = 0. _d 0
d216b1690c Mart*0351           ENDIF
                0352          ENDDO
                0353         ENDDO
                0354         DO iter=1,niter_bulk
                0355          DO j = 1,sNy
                0356           DO i = 1,sNx
                0357            IF ( atemp(i,j,bi,bj) .NE. 0. _d 0 ) THEN
                0358 C--- iterate with psi-functions to find transfer coefficients
                0359
bdec91d862 Patr*0360 #ifdef ALLOW_AUTODIFF_TAMC
20dee61641 Mart*0361             ikey_2 = i + (j-1)*sNx + (iter-1)*sNx*sNy
                0362      &           + (ikey-1)*sNx*sNy*niter_bulk
d216b1690c Mart*0363 CADJ STORE rdn   (i,j)       = comlev1_exf_2, key = ikey_2
                0364 CADJ STORE ustar (i,j)       = comlev1_exf_2, key = ikey_2
                0365 CADJ STORE qstar (i,j)       = comlev1_exf_2, key = ikey_2
                0366 CADJ STORE tstar (i,j)       = comlev1_exf_2, key = ikey_2
                0367 CADJ STORE sh    (i,j,bi,bj) = comlev1_exf_2, key = ikey_2
                0368 CADJ STORE wspeed(i,j,bi,bj) = comlev1_exf_2, key = ikey_2
bdec91d862 Patr*0369 #endif
                0370
d216b1690c Mart*0371             t0   = atemp(i,j,bi,bj)*
0b902863ee Mart*0372      &           (exf_one + humid_fac*aqh(i,j,bi,bj))
d216b1690c Mart*0373             huol = ( tstar(i,j)/t0 +
                0374      &               qstar(i,j)/(exf_one/humid_fac+aqh(i,j,bi,bj))
                0375      &              )*czol/(ustar(i,j)*ustar(i,j))
                0376 #ifdef ALLOW_BULK_LARGEYEAGER04
                0377             tmpbulk = MIN(abs(huol),10. _d 0)
                0378             huol   = SIGN(tmpbulk , huol)
                0379 #else
                0380 C--   Large&Pond1981:
0b902863ee Mart*0381             huol   = max(huol,zolmin)
d216b1690c Mart*0382 #endif /* ALLOW_BULK_LARGEYEAGER04 */
0b902863ee Mart*0383             htol   = huol*ht/hu
                0384             hqol   = huol*hq/hu
666d41f610 Mart*0385             stable = exf_half + sign(exf_half, huol)
d216b1690c Mart*0386
8053a926ae Jean*0387 C     Evaluate all stability functions assuming hq = ht.
d216b1690c Mart*0388             IF ( solve4Stress ) THEN
                0389 #ifdef ALLOW_BULK_LARGEYEAGER04
                0390 C--   Large&Yeager04:
                0391              xsq    = SQRT( ABS(exf_one - huol*16. _d 0) )
                0392 #else
                0393 C--   Large&Pond1981:
423768d890 Jean*0394              xsq    = MAX(SQRT(ABS(exf_one - 16.*huol)),exf_one)
d216b1690c Mart*0395 #endif /* ALLOW_BULK_LARGEYEAGER04 */
                0396              x      = SQRT(xsq)
                0397              psimh  = -psim_fac*huol*stable
                0398      &               +(exf_one-stable)*
                0399      &                ( LOG( (exf_one + exf_two*x + xsq)*
                0400      &                       (exf_one+xsq)*.125 _d 0 )
                0401      &                  -exf_two*ATAN(x) + exf_half*pi )
ddea17ebdb Gael*0402             ELSE
                0403              psimh  = 0. _d 0
d216b1690c Mart*0404             ENDIF
                0405 #ifdef ALLOW_BULK_LARGEYEAGER04
                0406 C--   Large&Yeager04:
                0407             xsq   = SQRT( ABS(exf_one - htol*16. _d 0) )
                0408 #else
                0409 C--   Large&Pond1981:
423768d890 Jean*0410             xsq   = MAX(SQRT(ABS(exf_one - 16.*htol)),exf_one)
d216b1690c Mart*0411 #endif /* ALLOW_BULK_LARGEYEAGER04 */
666d41f610 Mart*0412             psixh = -psim_fac*htol*stable + (exf_one-stable)*
                0413      &               ( exf_two*LOG(exf_half*(exf_one+xsq)) )
d216b1690c Mart*0414
                0415             IF ( solve4Stress ) THEN
7448700841 Mart*0416 #ifdef ALLOW_AUTODIFF_TAMC
8ed910e28e Patr*0417 CADJ STORE rdn(i,j)       = comlev1_exf_2, key = ikey_2
7448700841 Mart*0418 #endif
d216b1690c Mart*0419 C-   Shift wind speed using old coefficient
                0420 #ifdef ALLOW_BULK_LARGEYEAGER04
                0421 C--   Large&Yeager04:
                0422              dzTmp = (zwln-psimh)/karman
                0423              usn   = wspeed(i,j,bi,bj)/(exf_one + rdn(i,j)*dzTmp )
                0424 #else
                0425 C--   Large&Pond1981:
666d41f610 Mart*0426 c           rd(i,j)= rdn(i,j)/(exf_one - rdn(i,j)/karman*psimh )
d216b1690c Mart*0427 c           usn    = sh(i,j,bi,bj)*rd(i,j)/rdn(i,j)
8053a926ae Jean*0428 C     ML: the original formulation above is replaced below to be
                0429 C     similar to largeyeager04, but does not give the same results, strange
d216b1690c Mart*0430              usn   = sh(i,j,bi,bj)/(exf_one - rdn(i,j)/karman*psimh)
                0431 #endif /* ALLOW_BULK_LARGEYEAGER04 */
                0432              usm   = MAX(usn, umin)
                0433
                0434 C-   Update the 10m, neutral stability transfer coefficients (momentum)
cc60455fbb Mart*0435 #ifdef  ALLOW_DRAG_LARGEYEAGER09
                0436 C     Large and Yeager (2009), J.Clim equation 11a/b
                0437              tmpbulk =  cdrag_1/usm + cdrag_2 + cdrag_3*usm
                0438      &            + cdrag_8 * usm**6
                0439              tmpbulk = exf_scal_BulkCdn * (
                0440      &              ( halfRL - SIGN(halfRL, usm-umax) )*tmpbulk
                0441      &            + ( halfRL + SIGN(halfRL, usm-umax) )*cdragMax
                0442      &            )
                0443 #else
d216b1690c Mart*0444              tmpbulk  = exf_scal_BulkCdn *
                0445      &            ( cdrag_1/usm + cdrag_2 + cdrag_3*usm )
cc60455fbb Mart*0446 #endif
d216b1690c Mart*0447              rdn(i,j)   = SQRT(tmpbulk)
                0448 #ifdef ALLOW_BULK_LARGEYEAGER04
                0449              rd(i,j)    = rdn(i,j)/(exf_one + rdn(i,j)*dzTmp)
                0450 #else
                0451              rd(i,j)    = rdn(i,j)/(exf_one - rdn(i,j)/karman*psimh)
                0452 #endif /* ALLOW_BULK_LARGEYEAGER04 */
                0453              ustar(i,j) = rd(i,j)*sh(i,j,bi,bj)
                0454 C-   Coeff:
8053a926ae Jean*0455 #ifdef EXF_CALC_ATMRHO
                0456              tau(i,j)   = atmrho_loc(i,j)*rd(i,j)*wspeed(i,j,bi,bj)
                0457 #else
d216b1690c Mart*0458              tau(i,j)   = atmrho*rd(i,j)*wspeed(i,j,bi,bj)
8053a926ae Jean*0459 #endif
d216b1690c Mart*0460             ENDIF
                0461
                0462 C-   Update the 10m, neutral stability transfer coefficients (sens&evap)
666d41f610 Mart*0463             rhn = (exf_one-stable)*cstanton_1 + stable*cstanton_2
d216b1690c Mart*0464             ren = cDalton
                0465
                0466 C-   Shift all coefficients to the measurement height and stability.
                0467             rh = rhn/(exf_one + rhn*(ztln-psixh)/karman)
                0468             re = ren/(exf_one + ren*(ztln-psixh)/karman)
bdec91d862 Patr*0469
d216b1690c Mart*0470 C--  Update ustar, tstar, qstar using updated, shifted coefficients.
0b902863ee Mart*0471             qstar(i,j) = re*delq(i,j)
                0472             tstar(i,j) = rh*deltap(i,j)
f1644eaf0d Patr*0473
d216b1690c Mart*0474            ENDIF
                0475           ENDDO
                0476          ENDDO
8053a926ae Jean*0477 C end of iteration loop
d216b1690c Mart*0478         ENDDO
                0479         DO j = 1,sNy
                0480          DO i = 1,sNx
                0481           IF ( atemp(i,j,bi,bj) .NE. 0. _d 0 ) THEN
bdec91d862 Patr*0482
                0483 #ifdef ALLOW_AUTODIFF_TAMC
20dee61641 Mart*0484            ikey_1 = i + (j-1)*sNx + (ikey-1)*sNx*sNy
7311a000e1 Mart*0485 CADJ STORE qstar(i,j)    = comlev1_exf_1, key = ikey_1
                0486 CADJ STORE tstar(i,j)    = comlev1_exf_1, key = ikey_1
d216b1690c Mart*0487 CADJ STORE tau  (i,j)    = comlev1_exf_1, key = ikey_1
                0488 CADJ STORE rd   (i,j)    = comlev1_exf_1, key = ikey_1
                0489 #endif /* ALLOW_AUTODIFF_TAMC */
                0490
                0491 C-   Turbulent Fluxes
666d41f610 Mart*0492            hs(i,j,bi,bj) = atmcp*tau(i,j)*tstar(i,j)
                0493            hl(i,j,bi,bj) = flamb*tau(i,j)*qstar(i,j)
bdec91d862 Patr*0494 #ifndef EXF_READ_EVAP
d216b1690c Mart*0495 C   change sign and convert from kg/m^2/s to m/s via rhoConstFresh
                0496 c          evap(i,j,bi,bj) = -recip_rhonil*tau(i,j)*qstar(i,j)
666d41f610 Mart*0497            evap(i,j,bi,bj) = -recip_rhoConstFresh*tau(i,j)*qstar(i,j)
d216b1690c Mart*0498 C   but older version was using rhonil instead:
                0499 c           evap(i,j,bi,bj) = -recip_rhonil*tau(i,j)*qstar(i,j)
bdec91d862 Patr*0500 #endif
56d13a40ed Mart*0501            IF ( useAtmWind .AND. useRelativeWind ) THEN
                0502             tmpbulk =  tau(i,j)*rd(i,j)
                0503             ustress(i,j,bi,bj) = tmpbulk* ( uwind(i,j,bi,bj) -
                0504      &           0.5 _d 0 * (uVel(i,j,ks,bi,bj)+uVel(i+1,j,ks,bi,bj)) )
                0505             vstress(i,j,bi,bj) = tmpbulk* ( vwind(i,j,bi,bj) -
                0506      &           0.5 _d 0 * (vVel(i,j,ks,bi,bj)+vVel(i,j+1,ks,bi,bj)) )
                0507            ELSEIF ( useAtmWind ) THEN
423768d890 Jean*0508 c           ustress(i,j,bi,bj) = tau(i,j)*rd(i,j)*ws*cw(i,j,bi,bj)
                0509 c           vstress(i,j,bi,bj) = tau(i,j)*rd(i,j)*ws*sw(i,j,bi,bj)
d216b1690c Mart*0510 C- jmc: below is how it should be written ; different from above when
                0511 C       both wind-speed & 2 compon. of the wind are specified, and in
                0512 C-      this case, formula below is better.
423768d890 Jean*0513             tmpbulk =  tau(i,j)*rd(i,j)
                0514             ustress(i,j,bi,bj) = tmpbulk*uwind(i,j,bi,bj)
                0515             vstress(i,j,bi,bj) = tmpbulk*vwind(i,j,bi,bj)
                0516            ENDIF
7311a000e1 Mart*0517
                0518 c IF ( ATEMP(i,j,bi,bj) .NE. 0. )
423768d890 Jean*0519           ELSE
                0520            IF ( useAtmWind ) ustress(i,j,bi,bj) = 0. _d 0
                0521            IF ( useAtmWind ) vstress(i,j,bi,bj) = 0. _d 0
0b902863ee Mart*0522            hflux  (i,j,bi,bj) = 0. _d 0
d216b1690c Mart*0523            evap   (i,j,bi,bj) = 0. _d 0
                0524            hs     (i,j,bi,bj) = 0. _d 0
                0525            hl     (i,j,bi,bj) = 0. _d 0
7311a000e1 Mart*0526
                0527 c IF ( ATEMP(i,j,bi,bj) .NE. 0. )
423768d890 Jean*0528           ENDIF
                0529
d216b1690c Mart*0530          ENDDO
                0531         ENDDO
                0532        ENDDO
                0533       ENDDO
bdec91d862 Patr*0534
423768d890 Jean*0535 #endif /* ALLOW_ATM_TEMP */
7109a141b2 Patr*0536 #endif /* ALLOW_BULKFORMULAE */
                0537
d216b1690c Mart*0538       RETURN
                0539       END