pkg/thsice/thsice_get_exf.F

2a9474d935 Mart*0001 #include "THSICE_OPTIONS.h"
6b47d550f4 Mart*0002 #ifdef ALLOW_AUTODIFF
                0003 # include "AUTODIFF_OPTIONS.h"
                0004 #endif
2a9474d935 Mart*0005 #ifdef ALLOW_EXF
                0006 #include "EXF_OPTIONS.h"
                0007 #endif
                0008
                0009 CBOP
                0010 C     !ROUTINE: THSICE_GET_EXF
                0011 C     !INTERFACE:
                0012       SUBROUTINE THSICE_GET_EXF(
a163d11794 Patr*0013      I                  bi, bj, it2,
9dcf02c6ac Jean*0014      I                  iMin,iMax, jMin,jMax,
c1c3d0f9d7 Patr*0015      I                  icFlag, hSnow1, tsfCel,
9dcf02c6ac Jean*0016      O                  flxExcSw, dFlxdT, evapLoc, dEvdT,
                0017      I                  myTime, myIter, myThid )
                0018
2a9474d935 Mart*0019 C     !DESCRIPTION: \bv
                0020 C     *==========================================================*
                0021 C     | S/R  THSICE_GET_EXF
                0022 C     *==========================================================*
                0023 C     | Interface S/R : get Surface Fluxes from pkg EXF
                0024 C     *==========================================================*
                0025 C     \ev
                0026
                0027 C     !USES:
                0028       IMPLICIT NONE
                0029
                0030 C     == Global data ==
1fbbacd21d Jean*0031 #include "SIZE.h"
                0032 #include "EEPARAMS.h"
                0033 #include "PARAMS.h"
2a9474d935 Mart*0034 #ifdef ALLOW_EXF
681d55b157 Jean*0035 # include "EXF_CONSTANTS.h"
                0036 # include "EXF_PARAM.h"
                0037 # include "EXF_FIELDS.h"
2a9474d935 Mart*0038 #endif
77008a74ba Patr*0039 #ifdef ALLOW_AUTODIFF_TAMC
                0040 # include "tamc.h"
a163d11794 Patr*0041 # include "THSICE_SIZE.h"
77008a74ba Patr*0042 #endif
2a9474d935 Mart*0043
                0044 C     !INPUT/OUTPUT PARAMETERS:
                0045 C     === Routine arguments ===
9dcf02c6ac Jean*0046 C     bi,bj       :: tile indices
a163d11794 Patr*0047 C     it          :: solv4temp iteration
9dcf02c6ac Jean*0048 C     iMin,iMax   :: computation domain: 1rst index range
                0049 C     jMin,jMax   :: computation domain: 2nd  index range
c1c3d0f9d7 Patr*0050 C     icFlag     :: True= get fluxes at this location ; False= do nothing
                0051 C     hSnow1       :: snow height [m]
b20c2a3388 Jean*0052 C     tsfCel      :: surface (ice or snow) temperature (oC)
9dcf02c6ac Jean*0053 C     flxExcSw    :: net (downward) surface heat flux, except short-wave [W/m2]
                0054 C     dFlxdT      :: deriv of flx with respect to Tsf    [W/m/K]
2a9474d935 Mart*0055 C     evapLoc     :: surface evaporation (>0 if evaporate) [kg/m2/s]
                0056 C     dEvdT       :: deriv of evap. with respect to Tsf  [kg/m2/s/K]
9dcf02c6ac Jean*0057 C     myTime      :: current Time of simulation [s]
                0058 C     myIter      :: current Iteration number in simulation
                0059 C     myThid      :: my Thread Id number
                0060       INTEGER bi, bj
a163d11794 Patr*0061       INTEGER it2
9dcf02c6ac Jean*0062       INTEGER iMin, iMax
                0063       INTEGER jMin, jMax
c1c3d0f9d7 Patr*0064       _RL     icFlag  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
a85293d087 Mart*0065       _RL     hSnow1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
9dcf02c6ac Jean*0066       _RL     tsfCel  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0067       _RL     flxExcSw(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0068       _RL     dFlxdT  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0069       _RL     evapLoc (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0070       _RL     dEvdT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0071       _RL     myTime
                0072       INTEGER myIter
170766e9fd Jean*0073       INTEGER myThid
2a9474d935 Mart*0074 CEOP
                0075
                0076 #ifdef ALLOW_EXF
5465695fc7 Jean*0077 #ifdef ALLOW_ATM_TEMP
bd97d4af55 Jean*0078 #ifdef ALLOW_DOWNWARD_RADIATION
2a9474d935 Mart*0079
                0080 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0081 C     === Local variables ===
a3f1d55259 Jean*0082 C     hsLocal, hlLocal :: sensible & latent heat flux over sea-ice
5de01d9c5b Jean*0083 C             t0       :: virtual temperature (K)
                0084 C             ssq      :: saturation specific humidity (kg/kg)
                0085 C             deltap   :: potential temperature diff (K)
5f374c5868 Mart*0086       _RL hsLocal
                0087       _RL hlLocal
5de01d9c5b Jean*0088       INTEGER iter
9dcf02c6ac Jean*0089       INTEGER i, j
5de01d9c5b Jean*0090       _RL czol
                0091       _RL wsm                ! limited wind speed [m/s] (> umin)
                0092       _RL t0                 ! virtual temperature [K]
6de40c6370 Mart*0093 C     copied from exf_bulkformulae:
                0094 C     these need to be 2D-arrays for vectorizing code
                0095 C     turbulent temperature scale [K]
6575a2505e Jean*0096       _RL tstar (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
6de40c6370 Mart*0097 C     turbulent humidity scale  [kg/kg]
6575a2505e Jean*0098       _RL qstar (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
6de40c6370 Mart*0099 C     friction velocity [m/s]
                0100       _RL ustar (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0101 C     neutral, zref (=10m) values of rd
                0102       _RL rdn   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0103       _RL rd    (1-OLx:sNx+OLx,1-OLy:sNy+OLy) ! = sqrt(Cd)          [-]
                0104       _RL rh    (1-OLx:sNx+OLx,1-OLy:sNy+OLy) ! = Ch / sqrt(Cd)     [-]
                0105       _RL re    (1-OLx:sNx+OLx,1-OLy:sNy+OLy) ! = Ce / sqrt(Cd)     [-]
                0106 C     specific humidity difference [kg/kg]
6575a2505e Jean*0107       _RL delq  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
6de40c6370 Mart*0108       _RL deltap(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
8053a926ae Jean*0109 #ifdef EXF_CALC_ATMRHO
                0110 C     local atmospheric density [kg/m^3]
                0111       _RL atmrho_loc(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0112 #endif
6de40c6370 Mart*0113 C
a85293d087 Mart*0114       _RL ssq
6de40c6370 Mart*0115       _RL ren, rhn           ! neutral, zref (=10m) values of re, rh
5de01d9c5b Jean*0116       _RL usn, usm           ! neutral, zref (=10m) wind-speed (+limited)
                0117       _RL stable             ! = 1 if stable ; = 0 if unstable
6de40c6370 Mart*0118 C     stability parameter at zwd [-] (=z/Monin-Obuklov length)
6575a2505e Jean*0119       _RL huol
5de01d9c5b Jean*0120       _RL htol               ! stability parameter at zth [-]
                0121       _RL hqol
                0122       _RL x                  ! stability function  [-]
                0123       _RL xsq                ! = x^2               [-]
                0124       _RL psimh              ! momentum stability function
                0125       _RL psixh              ! latent & sensib. stability function
                0126       _RL zwln               ! = log(zwd/zref)
                0127       _RL ztln               ! = log(zth/zref)
                0128       _RL tau                ! surface stress coef = rhoA * Ws * sqrt(Cd)
                0129       _RL tmpbulk
2a9474d935 Mart*0130
                0131 C     additional variables that are copied from bulkf_formula_lay:
                0132 C     upward LW at surface (W m-2)
                0133       _RL  flwup
                0134 C     net (downward) LW at surface (W m-2)
                0135       _RL  flwNet_dwn
b20c2a3388 Jean*0136 C     gradients of latent/sensible net upward heat flux
2a9474d935 Mart*0137 C     w/ respect to temperature
5f374c5868 Mart*0138       _RL dflhdT
                0139       _RL dfshdT
6de40c6370 Mart*0140       _RL dflwupdT
2a9474d935 Mart*0141 C     emissivities, called emittance in exf
                0142       _RL     emiss
b20c2a3388 Jean*0143 C     Tsf    :: surface temperature [K]
                0144 C     Ts2    :: surface temperature square [K^2]
                0145       _RL Tsf
2a9474d935 Mart*0146       _RL Ts2
b20c2a3388 Jean*0147 C     latent heat of evaporation or sublimation [J/kg]
                0148       _RL lath
6de40c6370 Mart*0149       _RL qsat_fac
                0150       _RL qsat_exp
7c50f07931 Mart*0151 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0152 C     ikey  :: tape key (depends on tiles and external iteration)
                0153 C     nbkey :: tape key (depends on ikey and bulk formulae iteration)
                0154 C     ijkey :: tape key (depends on nbkey,  and i,j indices)
                0155       INTEGER ikey, nbkey, ijkey
7c50f07931 Mart*0156 #endif
6575a2505e Jean*0157 #ifdef ALLOW_DBUG_THSICE
                0158       LOGICAL dBugFlag
                0159       INTEGER stdUnit
                0160 #endif
2a9474d935 Mart*0161
b20c2a3388 Jean*0162 C     == external functions ==
2a9474d935 Mart*0163
b20c2a3388 Jean*0164 c     _RL       exf_BulkqSat
                0165 c     external  exf_BulkqSat
5de01d9c5b Jean*0166 c     _RL       exf_BulkCdn
                0167 c     external  exf_BulkCdn
                0168 c     _RL       exf_BulkRhn
                0169 c     external  exf_BulkRhn
2a9474d935 Mart*0170
b20c2a3388 Jean*0171 C     == end of interface ==
2a9474d935 Mart*0172
6575a2505e Jean*0173 C-    Define grid-point location where to print debugging values
                0174 #include "THSICE_DEBUG.h"
                0175
                0176 #ifdef ALLOW_DBUG_THSICE
ae605e558b Jean*0177       dBugFlag = debugLevel.GE.debLevC
6575a2505e Jean*0178       stdUnit = standardMessageUnit
                0179 #endif
                0180
6de40c6370 Mart*0181 C--   Set surface parameters :
                0182       zwln = LOG(hu/zref)
                0183       ztln = LOG(ht/zref)
                0184       czol = hu*karman*gravity_mks
                0185       ren  = cDalton
                0186 C     more abbreviations
                0187       lath     = flamb+flami
                0188       qsat_fac = cvapor_fac_ice
                0189       qsat_exp = cvapor_exp_ice
                0190
                0191 C     initialisation of local arrays
bd97d4af55 Jean*0192       DO j = 1-OLy,sNy+OLy
                0193        DO i = 1-OLx,sNx+OLx
6de40c6370 Mart*0194         tstar(i,j)  = 0. _d 0
                0195         qstar(i,j)  = 0. _d 0
                0196         ustar(i,j)  = 0. _d 0
                0197         rdn(i,j)    = 0. _d 0
                0198         rd(i,j)     = 0. _d 0
                0199         rh(i,j)     = 0. _d 0
                0200         re(i,j)     = 0. _d 0
                0201         delq(i,j)   = 0. _d 0
                0202         deltap(i,j) = 0. _d 0
                0203        ENDDO
                0204       ENDDO
                0205 C
9dcf02c6ac Jean*0206       DO j=jMin,jMax
                0207        DO i=iMin,iMax
6575a2505e Jean*0208 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0209 #ifdef ALLOW_DBUG_THSICE
8053a926ae Jean*0210         IF ( dBug(i,j,bi,bj) .AND. (icFlag(i,j).GT.0. _d 0) )
c1c3d0f9d7 Patr*0211      &    WRITE(stdUnit,'(A,2I4,2I2,2F12.6)')
                0212      &    'ThSI_GET_EXF: i,j,atemp,lwd=',
6575a2505e Jean*0213      &           i,j,bi,bj, atemp(i,j,bi,bj),lwdown(i,j,bi,bj)
                0214 #endif
9dcf02c6ac Jean*0215
                0216 C--   Use atmospheric state to compute surface fluxes.
8053a926ae Jean*0217         IF ( (icFlag(i,j).GT.0. _d 0) .AND.
c1c3d0f9d7 Patr*0218      &       (atemp(i,j,bi,bj).NE.0. _d 0) ) THEN
                0219          IF ( hSnow1(i,j).GT.3. _d -1 ) THEN
6de40c6370 Mart*0220           emiss = snow_emissivity
                0221          ELSE
                0222           emiss = ice_emissivity
                0223          ENDIF
2a9474d935 Mart*0224 C     copy a few variables to names used in bulkf_formula_lay
6de40c6370 Mart*0225          Tsf         = tsfCel(i,j)+cen2kel
                0226          Ts2         = Tsf*Tsf
5de01d9c5b Jean*0227 C     wind speed
6de40c6370 Mart*0228          wsm         = sh(i,j,bi,bj)
1fbbacd21d Jean*0229 C--   air - surface difference of temperature & humidity
6de40c6370 Mart*0230 c        tmpbulk= exf_BulkqSat(Tsf)
a85293d087 Mart*0231 c        ssq         = saltsat*tmpbulk/atmrho
6de40c6370 Mart*0232          tmpbulk     = qsat_fac*EXP(-qsat_exp/Tsf)
8053a926ae Jean*0233 #ifdef EXF_CALC_ATMRHO
                0234          atmrho_loc(i,j) = apressure(i,j,bi,bj) /
                0235      &                  (287.04 _d 0*atemp(i,j,bi,bj)
                0236      &                  *(1. _d 0 + humid_fac*aqh(i,j,bi,bj)))
a85293d087 Mart*0237          ssq         = tmpbulk/atmrho_loc(i,j)
8053a926ae Jean*0238 #else
a85293d087 Mart*0239          ssq         = tmpbulk/atmrho
8053a926ae Jean*0240 #endif
6de40c6370 Mart*0241          deltap(i,j) = atemp(i,j,bi,bj) + gamma_blk*ht - Tsf
a85293d087 Mart*0242          delq(i,j)   = aqh(i,j,bi,bj) - ssq
6575a2505e Jean*0243 C     Do the part of the output variables that do not depend
6de40c6370 Mart*0244 C     on the ice here to save a few re-computations
                0245 C     This is not yet dEvdT, but just a cheap way to save a 2D-field
                0246 C     for ssq and recomputing Ts2 lateron
a85293d087 Mart*0247          dEvdT(i,j)  = ssq*qsat_exp/Ts2
6de40c6370 Mart*0248          flwup       = emiss*stefanBoltzmann*Ts2*Ts2
                0249          dflwupdT    = emiss*stefanBoltzmann*Ts2*Tsf * 4. _d 0
                0250 c        flwNet_dwn  =       lwdown(i,j,bi,bj) - flwup
                0251 C-    assume long-wave albedo = 1 - emissivity :
                0252          flwNet_dwn  = emiss*lwdown(i,j,bi,bj) - flwup
                0253 C--   This is not yet the total derivative with respect to surface temperature
                0254          dFlxdT(i,j)   = -dflwupdT
                0255 C--   This is not yet the Net downward radiation excluding shortwave
                0256          flxExcSw(i,j) = flwNet_dwn
                0257         ENDIF
                0258        ENDDO
                0259       ENDDO
2a9474d935 Mart*0260
6575a2505e Jean*0261 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0262
6de40c6370 Mart*0263       IF ( useStabilityFct_overIce ) THEN
edb6656069 Mart*0264 #ifdef ALLOW_AUTODIFF_TAMC
                0265 C     Storing is only necessary if useStabilityFct_overIce = T, so we
                0266 C     compute the keys here, just before we need them.
                0267        ikey = bi + (bj-1)*nSx + (ikey_dynamics-1)*nSx*nSy
                0268        ikey = (ikey-1)*(MaxTsf+1) + it2
                0269 #endif
6de40c6370 Mart*0270        DO j=jMin,jMax
                0271         DO i=iMin,iMax
c1c3d0f9d7 Patr*0272 C--
8053a926ae Jean*0273          IF ( (icFlag(i,j).GT.0. _d 0) .AND.
c1c3d0f9d7 Patr*0274      &        (atemp(i,j,bi,bj).NE.0. _d 0) ) THEN
1fbbacd21d Jean*0275 C--   Compute the turbulent surface fluxes (function of stability).
                0276
b20c2a3388 Jean*0277 C             Initial guess: z/l=0.0; hu=ht=hq=z
                0278 C             Iterations:    converge on z/l and hence the fluxes.
2a9474d935 Mart*0279
6de40c6370 Mart*0280           t0         = atemp(i,j,bi,bj)*
                0281      &         (exf_one + humid_fac*aqh(i,j,bi,bj))
                0282           stable     = exf_half + SIGN(exf_half, deltap(i,j))
                0283 c         tmpbulk    = exf_BulkCdn(sh(i,j,bi,bj))
                0284           wsm        = sh(i,j,bi,bj)
                0285           tmpbulk    = cdrag_1/wsm + cdrag_2 + cdrag_3*wsm
d4b7a0d2bc Patr*0286           IF (tmpbulk.NE.0.) THEN
                0287            rdn(i,j)   = SQRT(tmpbulk)
                0288           ELSE
                0289            rdn(i,j)   = 0. _d 0
                0290           ENDIF
5de01d9c5b Jean*0291 C--  initial guess for exchange other coefficients:
6de40c6370 Mart*0292 c         rhn        = exf_BulkRhn(stable)
                0293           rhn        = (exf_one-stable)*cstanton_1 + stable*cstanton_2
5de01d9c5b Jean*0294 C--  calculate turbulent scales
6de40c6370 Mart*0295           ustar(i,j) = rdn(i,j)*wsm
                0296           tstar(i,j) = rhn*deltap(i,j)
                0297           qstar(i,j) = ren*delq(i,j)
                0298          ENDIF
                0299         ENDDO
                0300        ENDDO
c1c3d0f9d7 Patr*0301
6de40c6370 Mart*0302 C     start iteration
                0303        DO iter = 1,niter_bulk
a85293d087 Mart*0304 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0305         nbkey = (ikey-1)*niter_bulk + iter
                0306 CADJ STORE ustar = comlev1_thsice_exf_niter, key = nbkey, kind = isbyte
                0307 CADJ STORE qstar = comlev1_thsice_exf_niter, key = nbkey, kind = isbyte
                0308 CADJ STORE tstar = comlev1_thsice_exf_niter, key = nbkey, kind = isbyte
a85293d087 Mart*0309 #endif
6de40c6370 Mart*0310         DO j=jMin,jMax
                0311          DO i=iMin,iMax
8053a926ae Jean*0312           IF ( (icFlag(i,j).GT.0. _d 0) .AND.
c1c3d0f9d7 Patr*0313      &         (atemp(i,j,bi,bj).NE.0. _d 0) ) THEN
2a9474d935 Mart*0314
77008a74ba Patr*0315 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0316            ijkey = i + sNx*(j-1) + (nbkey-1)*sNx*sNy
a85293d087 Mart*0317 C     Not sure why this has to be done here and not outside the ij-loop.
edb6656069 Mart*0318 CADJ STORE rdn(i,j) = comlev1_thsice_exf_ij, key = ijkey, kind = isbyte
77008a74ba Patr*0319 #endif
                0320
6de40c6370 Mart*0321            t0     = atemp(i,j,bi,bj)*
                0322      &          (exf_one + humid_fac*aqh(i,j,bi,bj))
                0323            huol   = (tstar(i,j)/t0 +
                0324      &               qstar(i,j)/(exf_one/humid_fac+aqh(i,j,bi,bj))
                0325      &              )*czol/(ustar(i,j)*ustar(i,j))
                0326 #ifdef ALLOW_BULK_LARGEYEAGER04
5de01d9c5b Jean*0327 C-    Large&Yeager_2004 code:
6de40c6370 Mart*0328            huol   = MIN( MAX(-10. _d 0,huol), 10. _d 0 )
                0329 #else
                0330 C-    Large&Pond_1981 code (zolmin default = -100):
                0331            huol   = MAX(huol,zolmin)
                0332 #endif /* ALLOW_BULK_LARGEYEAGER04 */
                0333            htol   = huol*ht/hu
                0334            hqol   = huol*hq/hu
                0335            stable = exf_half + SIGN(exf_half, huol)
b20c2a3388 Jean*0336
                0337 C     Evaluate all stability functions assuming hq = ht.
6de40c6370 Mart*0338 #ifdef ALLOW_BULK_LARGEYEAGER04
5de01d9c5b Jean*0339 C-    Large&Yeager_2004 code:
6de40c6370 Mart*0340            xsq    = SQRT( ABS(exf_one - huol*16. _d 0) )
                0341 #else
5de01d9c5b Jean*0342 C-    Large&Pond_1981 code:
6de40c6370 Mart*0343            xsq    = MAX(SQRT(ABS(exf_one - huol*16. _d 0)),exf_one)
                0344 #endif /* ALLOW_BULK_LARGEYEAGER04 */
                0345            x      = SQRT(xsq)
                0346            psimh  = -psim_fac*huol*stable
5de01d9c5b Jean*0347      &             + (exf_one-stable)
6de40c6370 Mart*0348      &             *( LOG( (exf_one + exf_two*x + xsq)
                0349      &                    *(exf_one+xsq)*0.125 _d 0 )
                0350      &                -exf_two*ATAN(x) + exf_half*pi )
                0351 #ifdef ALLOW_BULK_LARGEYEAGER04
                0352 C-    Large&Yeager_2004 code:
                0353            xsq    = SQRT( ABS(exf_one - htol*16. _d 0) )
                0354 #else
                0355 C-    Large&Pond_1981 code:
                0356            xsq    = MAX(SQRT(ABS(exf_one - htol*16. _d 0)),exf_one)
                0357 #endif /* ALLOW_BULK_LARGEYEAGER04 */
                0358            psixh  = -psim_fac*htol*stable
                0359      &            + (exf_one-stable)
5de01d9c5b Jean*0360      &              *exf_two*LOG( exf_half*(exf_one+xsq) )
b20c2a3388 Jean*0361
                0362 C     Shift wind speed using old coefficient
6de40c6370 Mart*0363 #ifdef ALLOW_BULK_LARGEYEAGER04
                0364 C--   Large&Yeager04:
                0365            usn    = wspeed(i,j,bi,bj)
                0366      &           /( exf_one + rdn(i,j)*(zwln-psimh)/karman )
                0367 #else
                0368 C--   Large&Pond1981:
                0369            usn   = sh(i,j,bi,bj)/(exf_one - rdn(i,j)/karman*psimh)
                0370 #endif /* ALLOW_BULK_LARGEYEAGER04 */
                0371            usm    = MAX(usn, umin)
5de01d9c5b Jean*0372
                0373 C-    Update the 10m, neutral stability transfer coefficients
6de40c6370 Mart*0374 c          tmpbulk= exf_BulkCdn(usm)
                0375            tmpbulk= cdrag_1/usm + cdrag_2 + cdrag_3*usm
                0376            rdn(i,j) = SQRT(tmpbulk)
                0377 c          rhn    = exf_BulkRhn(stable)
                0378            rhn    = (exf_one-stable)*cstanton_1 + stable*cstanton_2
5de01d9c5b Jean*0379
                0380 C     Shift all coefficients to the measurement height and stability.
6de40c6370 Mart*0381 #ifdef ALLOW_BULK_LARGEYEAGER04
                0382            rd(i,j)= rdn(i,j)/( exf_one + rdn(i,j)*(zwln-psimh)/karman )
                0383 #else
                0384            rd(i,j)= rdn(i,j)/( exf_one - rdn(i,j)/karman*psimh )
                0385 #endif /* ALLOW_BULK_LARGEYEAGER04 */
                0386            rh(i,j)= rhn/( exf_one + rhn*(ztln-psixh)/karman )
                0387            re(i,j)= ren/( exf_one + ren*(ztln-psixh)/karman )
5de01d9c5b Jean*0388
                0389 C     Update ustar, tstar, qstar using updated, shifted coefficients.
6de40c6370 Mart*0390            ustar(i,j)  = rd(i,j)*sh(i,j,bi,bj)
                0391            qstar(i,j)  = re(i,j)*delq(i,j)
                0392            tstar(i,j)  = rh(i,j)*deltap(i,j)
                0393           ENDIF
                0394 C     end i/j-loops
                0395          ENDDO
                0396         ENDDO
                0397 C     end iteration loop
                0398        ENDDO
a85293d087 Mart*0399 #ifdef ALLOW_AUTODIFF_TAMC
                0400 C     these stores avoid recomputing this routine in the ad-code (not
                0401 C     strictly necessary)
edb6656069 Mart*0402 CADJ STORE rd    = comlev1_thsice_exf, key = ikey, kind = isbyte
                0403 CADJ STORE rh    = comlev1_thsice_exf, key = ikey, kind = isbyte
                0404 CADJ STORE re    = comlev1_thsice_exf, key = ikey, kind = isbyte
                0405 CADJ STORE qstar = comlev1_thsice_exf, key = ikey, kind = isbyte
                0406 CADJ STORE tstar = comlev1_thsice_exf, key = ikey, kind = isbyte
a85293d087 Mart*0407 #endif
6de40c6370 Mart*0408        DO j=jMin,jMax
                0409         DO i=iMin,iMax
8053a926ae Jean*0410          IF ( (icFlag(i,j).GT.0. _d 0) .AND.
c1c3d0f9d7 Patr*0411      &        (atemp(i,j,bi,bj).NE.0. _d 0) ) THEN
8053a926ae Jean*0412 #ifdef EXF_CALC_ATMRHO
                0413           tau     = atmrho_loc(i,j)*rd(i,j)*wspeed(i,j,bi,bj)
                0414 #else
6de40c6370 Mart*0415           tau     = atmrho*rd(i,j)*wspeed(i,j,bi,bj)
8053a926ae Jean*0416 #endif
5f374c5868 Mart*0417           evapLoc(i,j)  = -tau*qstar(i,j)
                0418           hlLocal       = -lath*evapLoc(i,j)
                0419           hsLocal       = atmcp*tau*tstar(i,j)
6de40c6370 Mart*0420 c         ustress = tau*rd(i,j)*UwindSpeed
                0421 c         vstress = tau*rd(i,j)*VwindSpeed
2a9474d935 Mart*0422
                0423 C---  surf.Temp derivative of turbulent Fluxes
6de40c6370 Mart*0424 C     complete computation of dEvdT
5f374c5868 Mart*0425           dEvdT(i,j)    = (tau*re(i,j))*dEvdT(i,j)
                0426           dflhdT        = -lath*dEvdT(i,j)
                0427           dfshdT        = -atmcp*tau*rh(i,j)
                0428 C--   Update total derivative with respect to surface temperature
                0429           dFlxdT(i,j)   = dFlxdT(i,j)   + dfshdT  + dflhdT
                0430 C--   Update net downward radiation excluding shortwave
                0431           flxExcSw(i,j) = flxExcSw(i,j) + hsLocal + hlLocal
2a9474d935 Mart*0432
6de40c6370 Mart*0433          ENDIF
                0434         ENDDO
                0435        ENDDO
                0436       ELSE
6575a2505e Jean*0437 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
1fbbacd21d Jean*0438 C--   Compute the turbulent surface fluxes using fixed transfert Coeffs
6de40c6370 Mart*0439 C     with no stability dependence ( useStabilityFct_overIce = false )
                0440        DO j=jMin,jMax
                0441         DO i=iMin,iMax
8053a926ae Jean*0442          IF ( (icFlag(i,j).GT.0. _d 0) .AND.
c1c3d0f9d7 Patr*0443      &        (atemp(i,j,bi,bj).NE.0. _d 0) ) THEN
5f374c5868 Mart*0444           wsm           = sh(i,j,bi,bj)
8053a926ae Jean*0445 #ifdef EXF_CALC_ATMRHO
                0446           tau           = atmrho_loc(i,j)*exf_iceCe*wsm
                0447 #else
5f374c5868 Mart*0448           tau           = atmrho*exf_iceCe*wsm
8053a926ae Jean*0449 #endif
5f374c5868 Mart*0450           evapLoc(i,j)  = -tau*delq(i,j)
                0451           hlLocal       = -lath*evapLoc(i,j)
8053a926ae Jean*0452 #ifdef EXF_CALC_ATMRHO
                0453           hsLocal       = atmcp*atmrho_loc(i,j)
                0454      &                                *exf_iceCh*wsm*deltap(i,j)
                0455 #else
5f374c5868 Mart*0456           hsLocal       = atmcp*atmrho*exf_iceCh*wsm*deltap(i,j)
8053a926ae Jean*0457 #endif
6575a2505e Jean*0458 #ifdef ALLOW_DBUG_THSICE
                0459           IF ( dBug(i,j,bi,bj) ) WRITE(stdUnit,'(A,4F12.6)')
                0460      &      'ThSI_GET_EXF: wsm,hl,hs,Lw=',
                0461      &                     wsm,hlLocal,hsLocal,flxExcSw(i,j)
                0462 #endif
1fbbacd21d Jean*0463 C---  surf.Temp derivative of turbulent Fluxes
6de40c6370 Mart*0464 C     complete computation of dEvdT
5f374c5868 Mart*0465           dEvdT(i,j)    = tau*dEvdT(i,j)
                0466           dflhdT        = -lath*dEvdT(i,j)
8053a926ae Jean*0467 #ifdef EXF_CALC_ATMRHO
                0468           dfshdT        = -atmcp*atmrho_loc(i,j)*exf_iceCh*wsm
                0469 #else
5f374c5868 Mart*0470           dfshdT        = -atmcp*atmrho*exf_iceCh*wsm
8053a926ae Jean*0471 #endif
5f374c5868 Mart*0472 C--   Update total derivative with respect to surface temperature
                0473           dFlxdT(i,j)   = dFlxdT(i,j)   + dfshdT  + dflhdT
                0474 C--   Update net downward radiation excluding shortwave
                0475           flxExcSw(i,j) = flxExcSw(i,j) + hsLocal + hlLocal
6575a2505e Jean*0476 #ifdef ALLOW_DBUG_THSICE
                0477           IF ( dBug(i,j,bi,bj) ) WRITE(stdUnit,'(A,4F12.6)')
                0478      &      'ThSI_GET_EXF: flx,dFlxdT,evap,dEvdT',
                0479      &       flxExcSw(i,j), dFlxdT(i,j), evapLoc(i,j),dEvdT(i,j)
                0480 #endif
6de40c6370 Mart*0481          ENDIF
                0482         ENDDO
                0483        ENDDO
                0484 C     endif useStabilityFct_overIce
                0485       ENDIF
6575a2505e Jean*0486 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
6de40c6370 Mart*0487       DO j=jMin,jMax
                0488        DO i=iMin,iMax
8053a926ae Jean*0489         IF ( (icFlag(i,j).GT.0. _d 0) .AND.
c1c3d0f9d7 Patr*0490      &       (atemp(i,j,bi,bj).LE.0. _d 0) ) THEN
6de40c6370 Mart*0491 C--   in case atemp is zero:
                0492          flxExcSw(i,j) = 0. _d 0
                0493          dFlxdT  (i,j) = 0. _d 0
                0494          evapLoc (i,j) = 0. _d 0
                0495          dEvdT   (i,j) = 0. _d 0
1fbbacd21d Jean*0496         ENDIF
9dcf02c6ac Jean*0497        ENDDO
                0498       ENDDO
2a9474d935 Mart*0499
bd97d4af55 Jean*0500 #else /* ALLOW_DOWNWARD_RADIATION */
                0501       STOP 'ABNORMAL END: S/R THSICE_GET_EXF: DOWNWARD_RADIATION undef'
                0502 #endif /* ALLOW_DOWNWARD_RADIATION */
5465695fc7 Jean*0503 #else /* ALLOW_ATM_TEMP */
                0504       STOP 'ABNORMAL END: S/R THSICE_GET_EXF: ATM_TEMP undef'
                0505 #endif /* ALLOW_ATM_TEMP */
5de01d9c5b Jean*0506 #ifdef EXF_READ_EVAP
                0507       STOP 'ABNORMAL END: S/R THSICE_GET_EXF: EXF_READ_EVAP defined'
                0508 #endif /* EXF_READ_EVAP */
2a9474d935 Mart*0509 #endif /* ALLOW_EXF */
                0510
                0511       RETURN
                0512       END