/[MITgcm]/MITgcm/pkg/thsice/thsice_get_exf.F
ViewVC logotype

Diff of /MITgcm/pkg/thsice/thsice_get_exf.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.11 by jmc, Mon May 14 20:48:26 2007 UTC revision 1.18 by jmc, Fri Dec 24 00:53:42 2010 UTC
# Line 10  CBOP Line 10  CBOP
10  C     !ROUTINE: THSICE_GET_EXF  C     !ROUTINE: THSICE_GET_EXF
11  C     !INTERFACE:  C     !INTERFACE:
12        SUBROUTINE THSICE_GET_EXF(        SUBROUTINE THSICE_GET_EXF(
13       I                         iceornot, tsfCel,       I                  bi, bj,
14       O                         flxExceptSw, df0dT, evapLoc, dEvdT,       I                  iMin,iMax, jMin,jMax,
15       I                         i,j,bi,bj,myThid )       I                  iceFlag, hSnow, tsfCel,
16         O                  flxExcSw, dFlxdT, evapLoc, dEvdT,
17         I                  myTime, myIter, myThid )
18    
19  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
20  C     *==========================================================*  C     *==========================================================*
21  C     | S/R  THSICE_GET_EXF  C     | S/R  THSICE_GET_EXF
# Line 40  C     == Global data == Line 43  C     == Global data ==
43    
44  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
45  C     === Routine arguments ===  C     === Routine arguments ===
46  C     iceornot    :: 0=open water, 1=ice cover  C     bi,bj       :: tile indices
47    C     iMin,iMax   :: computation domain: 1rst index range
48    C     jMin,jMax   :: computation domain: 2nd  index range
49    C     iceFlag     :: True= get fluxes at this location ; False= do nothing
50    C     hSnow       :: snow height [m]
51  C     tsfCel      :: surface (ice or snow) temperature (oC)  C     tsfCel      :: surface (ice or snow) temperature (oC)
52  C     flxExceptSw :: net (downward) surface heat flux, except short-wave [W/m2]  C     flxExcSw    :: net (downward) surface heat flux, except short-wave [W/m2]
53  C     df0dT       :: deriv of flx with respect to Tsf    [W/m/K]  C     dFlxdT      :: deriv of flx with respect to Tsf    [W/m/K]
54  C     evapLoc     :: surface evaporation (>0 if evaporate) [kg/m2/s]  C     evapLoc     :: surface evaporation (>0 if evaporate) [kg/m2/s]
55  C     dEvdT       :: deriv of evap. with respect to Tsf  [kg/m2/s/K]  C     dEvdT       :: deriv of evap. with respect to Tsf  [kg/m2/s/K]
56  C     i,j, bi,bj  :: current grid point indices  C     myTime      :: current Time of simulation [s]
57  C     myThid      :: My Thread Id number  C     myIter      :: current Iteration number in simulation
58        INTEGER i,j, bi,bj  C     myThid      :: my Thread Id number
59          INTEGER bi, bj
60          INTEGER iMin, iMax
61          INTEGER jMin, jMax
62          LOGICAL iceFlag (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
63          _RL     hSnow   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
64          _RL     tsfCel  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
65          _RL     flxExcSw(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
66          _RL     dFlxdT  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
67          _RL     evapLoc (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
68          _RL     dEvdT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
69          _RL     myTime
70          INTEGER myIter
71        INTEGER myThid        INTEGER myThid
       INTEGER iceornot  
       _RL  tsfCel  
       _RL  flxExceptSw  
       _RL  df0dT  
       _RL  evapLoc  
       _RL  dEvdT  
72  CEOP  CEOP
73    
74  #ifdef ALLOW_EXF  #ifdef ALLOW_EXF
# Line 67  C     hsLocal, hlLocal :: sensible & lat Line 80  C     hsLocal, hlLocal :: sensible & lat
80  C             t0       :: virtual temperature (K)  C             t0       :: virtual temperature (K)
81  C             ssq      :: saturation specific humidity (kg/kg)  C             ssq      :: saturation specific humidity (kg/kg)
82  C             deltap   :: potential temperature diff (K)  C             deltap   :: potential temperature diff (K)
83          _RL hsLocal
84        _RL     hsLocal, hlLocal        _RL hlLocal
85        INTEGER iter        INTEGER iter
86        _RL     delq        INTEGER i, j
       _RL     deltap  
87        _RL czol        _RL czol
       _RL ws                 ! wind speed [m/s] (unlimited)  
88        _RL wsm                ! limited wind speed [m/s] (> umin)        _RL wsm                ! limited wind speed [m/s] (> umin)
89        _RL t0                 ! virtual temperature [K]        _RL t0                 ! virtual temperature [K]
90        _RL ustar              ! friction velocity [m/s]  C     copied from exf_bulkformulae:
91        _RL tstar              ! turbulent temperature scale [K]  C     these need to be 2D-arrays for vectorizing code
92        _RL qstar              ! turbulent humidity scale  [kg/kg]  C     turbulent temperature scale [K]
93        _RL ssq        _RL tstar (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
94        _RL rd                 ! = sqrt(Cd)          [-]  C     turbulent humidity scale  [kg/kg]
95        _RL re                 ! = Ce / sqrt(Cd)     [-]        _RL qstar (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
96        _RL rh                 ! = Ch / sqrt(Cd)     [-]  C     friction velocity [m/s]
97        _RL rdn, ren, rhn      ! neutral, zref (=10m) values of rd, re, rh        _RL ustar (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
98    C     neutral, zref (=10m) values of rd
99          _RL rdn   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
100          _RL rd    (1-OLx:sNx+OLx,1-OLy:sNy+OLy) ! = sqrt(Cd)          [-]
101          _RL rh    (1-OLx:sNx+OLx,1-OLy:sNy+OLy) ! = Ch / sqrt(Cd)     [-]
102          _RL re    (1-OLx:sNx+OLx,1-OLy:sNy+OLy) ! = Ce / sqrt(Cd)     [-]
103    C     specific humidity difference [kg/kg]
104          _RL delq  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
105          _RL deltap(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
106    C
107          _RL ssq   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
108          _RL ren, rhn           ! neutral, zref (=10m) values of re, rh
109        _RL usn, usm           ! neutral, zref (=10m) wind-speed (+limited)        _RL usn, usm           ! neutral, zref (=10m) wind-speed (+limited)
110        _RL stable             ! = 1 if stable ; = 0 if unstable        _RL stable             ! = 1 if stable ; = 0 if unstable
111        _RL huol               ! stability parameter at zwd [-] (=z/Monin-Obuklov length)  C     stability parameter at zwd [-] (=z/Monin-Obuklov length)
112          _RL huol
113        _RL htol               ! stability parameter at zth [-]        _RL htol               ! stability parameter at zth [-]
114        _RL hqol        _RL hqol
115        _RL x                  ! stability function  [-]        _RL x                  ! stability function  [-]
# Line 105  C     net (downward) LW at surface (W m- Line 128  C     net (downward) LW at surface (W m-
128        _RL  flwNet_dwn        _RL  flwNet_dwn
129  C     gradients of latent/sensible net upward heat flux  C     gradients of latent/sensible net upward heat flux
130  C     w/ respect to temperature  C     w/ respect to temperature
131        _RL dflhdT, dfshdT, dflwupdT        _RL dflhdT
132          _RL dfshdT
133          _RL dflwupdT
134  C     emissivities, called emittance in exf  C     emissivities, called emittance in exf
135        _RL     emiss        _RL     emiss
136  C     Tsf    :: surface temperature [K]  C     Tsf    :: surface temperature [K]
# Line 114  C     Ts2    :: surface temperature squa Line 139  C     Ts2    :: surface temperature squa
139        _RL Ts2        _RL Ts2
140  C     latent heat of evaporation or sublimation [J/kg]  C     latent heat of evaporation or sublimation [J/kg]
141        _RL lath        _RL lath
142        _RL qsat_fac, qsat_exp        _RL qsat_fac
143          _RL qsat_exp
144  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
145        INTEGER ikey_1  c     INTEGER ikey_1
146        INTEGER ikey_2  c     INTEGER ikey_2
147    #endif
148    #ifdef ALLOW_DBUG_THSICE
149          LOGICAL dBugFlag
150          INTEGER stdUnit
151  #endif  #endif
152    
153  C     == external functions ==  C     == external functions ==
# Line 131  c     external  exf_BulkRhn Line 161  c     external  exf_BulkRhn
161    
162  C     == end of interface ==  C     == end of interface ==
163    
164    C-    Define grid-point location where to print debugging values
165    #include "THSICE_DEBUG.h"
166    
167    #ifdef ALLOW_DBUG_THSICE
168          dBugFlag = debugLevel.GE.debLevB
169          stdUnit = standardMessageUnit
170    #endif
171    
172  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
173            act1 = bi - myBxLo(myThid)  c         act1 = bi - myBxLo(myThid)
174            max1 = myBxHi(myThid) - myBxLo(myThid) + 1  c         max1 = myBxHi(myThid) - myBxLo(myThid) + 1
175            act2 = bj - myByLo(myThid)  c         act2 = bj - myByLo(myThid)
176            max2 = myByHi(myThid) - myByLo(myThid) + 1  c         max2 = myByHi(myThid) - myByLo(myThid) + 1
177            act3 = myThid - 1  c         act3 = myThid - 1
178            max3 = nTx*nTy  c         max3 = nTx*nTy
179            act4 = ikey_dynamics - 1  c         act4 = ikey_dynamics - 1
180    c         ikey_1 = i
181            ikey_1 = i  c    &         + sNx*(j-1)
182       &         + sNx*(j-1)  c    &         + sNx*sNy*act1
183       &         + sNx*sNy*act1  c    &         + sNx*sNy*max1*act2
184       &         + sNx*sNy*max1*act2  c    &         + sNx*sNy*max1*max2*act3
185       &         + sNx*sNy*max1*max2*act3  c    &         + sNx*sNy*max1*max2*max3*act4
      &         + sNx*sNy*max1*max2*max3*act4  
186  #endif  #endif
187    
188    C--   Set surface parameters :
189          zwln = LOG(hu/zref)
190          ztln = LOG(ht/zref)
191          czol = hu*karman*gravity_mks
192          ren  = cDalton
193    C     more abbreviations
194          lath     = flamb+flami
195          qsat_fac = cvapor_fac_ice
196          qsat_exp = cvapor_exp_ice
197    
198    C     initialisation of local arrays
199          DO j = 1-Oly,sNy+Oly
200           DO i = 1-Olx,sNx+Olx
201            tstar(i,j)  = 0. _d 0
202            qstar(i,j)  = 0. _d 0
203            ustar(i,j)  = 0. _d 0
204            rdn(i,j)    = 0. _d 0
205            rd(i,j)     = 0. _d 0
206            rh(i,j)     = 0. _d 0
207            re(i,j)     = 0. _d 0
208            delq(i,j)   = 0. _d 0
209            deltap(i,j) = 0. _d 0
210            ssq(i,j)    = 0. _d 0
211           ENDDO
212          ENDDO
213    C
214          DO j=jMin,jMax
215           DO i=iMin,iMax
216    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
217    #ifdef ALLOW_DBUG_THSICE
218            IF ( dBug(i,j,bi,bj) .AND. iceFlag(i,j) ) WRITE(stdUnit,
219         &    '(A,2I4,2I2,2F12.6)') 'ThSI_GET_EXF: i,j,atemp,lwd=',
220         &           i,j,bi,bj, atemp(i,j,bi,bj),lwdown(i,j,bi,bj)
221    #endif
222    
223    C--   Use atmospheric state to compute surface fluxes.
224            IF ( iceFlag(i,j) .AND. atemp(i,j,bi,bj).NE.0. _d 0 ) THEN
225             IF ( hSnow(i,j).GT.3. _d -1 ) THEN
226              emiss = snow_emissivity
227             ELSE
228              emiss = ice_emissivity
229             ENDIF
230  C     copy a few variables to names used in bulkf_formula_lay  C     copy a few variables to names used in bulkf_formula_lay
231        Tsf    =  tsfCel+cen2kel           Tsf         = tsfCel(i,j)+cen2kel
232        Ts2    = Tsf*Tsf           Ts2         = Tsf*Tsf
233  C     wind speed  C     wind speed
       ws     = wspeed(i,j,bi,bj)  
234  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
235  CADJ STORE sh(i,j,bi,bj)     = comlev1_exf_1, key = ikey_1  cCADJ STORE sh(i,j,bi,bj)     = comlev1_exf_1, key = ikey_1
236  #endif  #endif
237        wsm    = sh(i,j,bi,bj)           wsm         = sh(i,j,bi,bj)
       IF ( iceornot.EQ.0 ) THEN  
         lath  = flamb  
         qsat_fac = cvapor_fac  
         qsat_exp = cvapor_exp  
       ELSE  
         lath  = flamb+flami  
         qsat_fac = cvapor_fac_ice  
         qsat_exp = cvapor_exp_ice  
       ENDIF  
   
 C--   Use atmospheric state to compute surface fluxes.  
         IF ( atemp(i,j,bi,bj) .NE. 0. _d 0 ) THEN  
   
238  C--   air - surface difference of temperature & humidity  C--   air - surface difference of temperature & humidity
239  c         tmpbulk= exf_BulkqSat(Tsf)  c        tmpbulk= exf_BulkqSat(Tsf)
240  c         ssq    = saltsat*tmpbulk/atmrho  c        ssq(i,j)    = saltsat*tmpbulk/atmrho
241            tmpbulk= qsat_fac*EXP(-qsat_exp/Tsf)           tmpbulk     = qsat_fac*EXP(-qsat_exp/Tsf)
242            ssq    = tmpbulk/atmrho           ssq(i,j)    = tmpbulk/atmrho
243            deltap = atemp(i,j,bi,bj) + gamma_blk*ht - Tsf           deltap(i,j) = atemp(i,j,bi,bj) + gamma_blk*ht - Tsf
244            delq   = aqh(i,j,bi,bj) - ssq           delq(i,j)   = aqh(i,j,bi,bj) - ssq(i,j)
245    C     Do the part of the output variables that do not depend
246    C     on the ice here to save a few re-computations
247    C     This is not yet dEvdT, but just a cheap way to save a 2D-field
248    C     for ssq and recomputing Ts2 lateron
249             dEvdT(i,j)  = ssq(i,j)*qsat_exp/Ts2
250             flwup       = emiss*stefanBoltzmann*Ts2*Ts2
251             dflwupdT    = emiss*stefanBoltzmann*Ts2*Tsf * 4. _d 0
252    #ifdef ALLOW_DOWNWARD_RADIATION
253    c        flwNet_dwn  =       lwdown(i,j,bi,bj) - flwup
254    C-    assume long-wave albedo = 1 - emissivity :
255             flwNet_dwn  = emiss*lwdown(i,j,bi,bj) - flwup
256    #else
257             STOP
258         &     'ABNORMAL END: S/R THSICE_GET_EXF: DOWNWARD_RADIATION undef'
259    #endif
260    C--   This is not yet the total derivative with respect to surface temperature
261             dFlxdT(i,j)   = -dflwupdT
262    C--   This is not yet the Net downward radiation excluding shortwave
263             flxExcSw(i,j) = flwNet_dwn
264            ENDIF
265           ENDDO
266          ENDDO
267    
268            IF ( useStabilityFct_overIce ) THEN  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
269    
270          IF ( useStabilityFct_overIce ) THEN
271           DO j=jMin,jMax
272            DO i=iMin,iMax
273             IF ( iceFlag(i,j) .AND. atemp(i,j,bi,bj).NE.0. _d 0 ) THEN
274  C--   Compute the turbulent surface fluxes (function of stability).  C--   Compute the turbulent surface fluxes (function of stability).
275    
 C--   Set surface parameters :  
            zwln = LOG(hu/zref)  
            ztln = LOG(ht/zref)  
            czol = hu*karman*gravity_mks  
276    
277  C             Initial guess: z/l=0.0; hu=ht=hq=z  C             Initial guess: z/l=0.0; hu=ht=hq=z
278  C             Iterations:    converge on z/l and hence the fluxes.  C             Iterations:    converge on z/l and hence the fluxes.
279    
280             t0     = atemp(i,j,bi,bj)*            t0         = atemp(i,j,bi,bj)*
281       &            (exf_one + humid_fac*aqh(i,j,bi,bj))       &         (exf_one + humid_fac*aqh(i,j,bi,bj))
282             stable = exf_half + SIGN(exf_half, deltap)            stable     = exf_half + SIGN(exf_half, deltap(i,j))
283  c          tmpbulk= exf_BulkCdn(sh(i,j,bi,bj))  c         tmpbulk    = exf_BulkCdn(sh(i,j,bi,bj))
284             tmpbulk= cdrag_1/wsm + cdrag_2 + cdrag_3*wsm            wsm        = sh(i,j,bi,bj)
285             rdn    = SQRT(tmpbulk)            tmpbulk    = cdrag_1/wsm + cdrag_2 + cdrag_3*wsm
286              IF (tmpbulk.NE.0.) THEN
287               rdn(i,j)   = SQRT(tmpbulk)
288              ELSE
289               rdn(i,j)   = 0. _d 0
290              ENDIF
291  C--  initial guess for exchange other coefficients:  C--  initial guess for exchange other coefficients:
292  c          rhn    = exf_BulkRhn(stable)  c         rhn        = exf_BulkRhn(stable)
293             rhn    = (exf_one-stable)*cstanton_1 + stable*cstanton_2            rhn        = (exf_one-stable)*cstanton_1 + stable*cstanton_2
            ren    = cDalton  
294  C--  calculate turbulent scales  C--  calculate turbulent scales
295             ustar  = rdn*wsm            ustar(i,j) = rdn(i,j)*wsm
296             tstar  = rhn*deltap            tstar(i,j) = rhn*deltap(i,j)
297             qstar  = ren*delq            qstar(i,j) = ren*delq(i,j)
298             ENDIF
299             DO iter = 1,niter_bulk          ENDDO
300           ENDDO
301    C     start iteration
302           DO iter = 1,niter_bulk
303            DO j=jMin,jMax
304             DO i=iMin,iMax
305              IF ( iceFlag(i,j) .AND. atemp(i,j,bi,bj).NE.0. _d 0 ) THEN
306    
307  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
308                     ikey_2 = iter  c                  ikey_2 = iter
309       &                  + niter_bulk*(i-1)  c    &                  + niter_bulk*(i-1)
310       &                  + niter_bulk*sNx*(j-1)  c    &                  + niter_bulk*sNx*(j-1)
311       &                  + niter_bulk*sNx*sNy*act1  c    &                  + niter_bulk*sNx*sNy*act1
312       &                  + niter_bulk*sNx*sNy*max1*act2  c    &                  + niter_bulk*sNx*sNy*max1*act2
313       &                  + niter_bulk*sNx*sNy*max1*max2*act3  c    &                  + niter_bulk*sNx*sNy*max1*max2*act3
314       &                  + niter_bulk*sNx*sNy*max1*max2*max3*act4  c    &                  + niter_bulk*sNx*sNy*max1*max2*max3*act4
315    cCADJ STORE rdn    = comlev1_exf_2, key = ikey_2
316  CADJ STORE rdn    = comlev1_exf_2, key = ikey_2  cCADJ STORE ustar  = comlev1_exf_2, key = ikey_2
317  CADJ STORE ustar  = comlev1_exf_2, key = ikey_2  cCADJ STORE qstar  = comlev1_exf_2, key = ikey_2
318  CADJ STORE qstar  = comlev1_exf_2, key = ikey_2  cCADJ STORE tstar  = comlev1_exf_2, key = ikey_2
319  CADJ STORE tstar  = comlev1_exf_2, key = ikey_2  cCADJ STORE sh(i,j,bi,bj)     = comlev1_exf_2, key = ikey_2
320  CADJ STORE sh(i,j,bi,bj)     = comlev1_exf_2, key = ikey_2  #endif
321  #endif  
322               t0     = atemp(i,j,bi,bj)*
323              huol   = (tstar/t0 +       &          (exf_one + humid_fac*aqh(i,j,bi,bj))
324       &                qstar/(exf_one/humid_fac+aqh(i,j,bi,bj))             huol   = (tstar(i,j)/t0 +
325       &               )*czol/(ustar*ustar)       &               qstar(i,j)/(exf_one/humid_fac+aqh(i,j,bi,bj))
326  C-    Large&Pond_1981 code (zolmin default = -100):       &              )*czol/(ustar(i,j)*ustar(i,j))
327  c           huol   = MAX(huol,zolmin)  #ifdef ALLOW_BULK_LARGEYEAGER04
328  C-    Large&Yeager_2004 code:  C-    Large&Yeager_2004 code:
329              huol   = MIN( MAX(-10. _d 0,huol), 10. _d 0 )             huol   = MIN( MAX(-10. _d 0,huol), 10. _d 0 )
330              htol   = huol*ht/hu  #else
331              hqol   = huol*hq/hu  C-    Large&Pond_1981 code (zolmin default = -100):
332              stable = exf_half + SIGN(exf_half, huol)             huol   = MAX(huol,zolmin)
333    #endif /* ALLOW_BULK_LARGEYEAGER04 */
334               htol   = huol*ht/hu
335               hqol   = huol*hq/hu
336               stable = exf_half + SIGN(exf_half, huol)
337    
338  C     Evaluate all stability functions assuming hq = ht.  C     Evaluate all stability functions assuming hq = ht.
339  C-    Large&Pond_1981 code:  #ifdef ALLOW_BULK_LARGEYEAGER04
 c           xsq    = MAX(SQRT(ABS(exf_one - huol*16. _d 0)),exf_one)  
340  C-    Large&Yeager_2004 code:  C-    Large&Yeager_2004 code:
341              xsq    = SQRT( ABS(exf_one - huol*16. _d 0) )             xsq    = SQRT( ABS(exf_one - huol*16. _d 0) )
342              x      = SQRT(xsq)  #else
             psimh  = -psim_fac*huol*stable  
      &             + (exf_one-stable)  
      &              *( LOG( (exf_one + exf_two*x + xsq)  
      &                     *(exf_one+xsq)*0.125 _d 0 )  
      &                 -exf_two*ATAN(x) + exf_half*pi )  
343  C-    Large&Pond_1981 code:  C-    Large&Pond_1981 code:
344  c           xsq    = MAX(SQRT(ABS(exf_one - htol*16. _d 0)),exf_one)             xsq    = MAX(SQRT(ABS(exf_one - huol*16. _d 0)),exf_one)
345  C-    Large&Yeager_2004 code:  #endif /* ALLOW_BULK_LARGEYEAGER04 */
346              xsq    = SQRT( ABS(exf_one - htol*16. _d 0) )             x      = SQRT(xsq)
347              psixh  = -psim_fac*htol*stable             psimh  = -psim_fac*huol*stable
348       &             + (exf_one-stable)       &             + (exf_one-stable)
349         &             *( LOG( (exf_one + exf_two*x + xsq)
350         &                    *(exf_one+xsq)*0.125 _d 0 )
351         &                -exf_two*ATAN(x) + exf_half*pi )
352    #ifdef ALLOW_BULK_LARGEYEAGER04
353    C-    Large&Yeager_2004 code:
354               xsq    = SQRT( ABS(exf_one - htol*16. _d 0) )
355    #else
356    C-    Large&Pond_1981 code:
357               xsq    = MAX(SQRT(ABS(exf_one - htol*16. _d 0)),exf_one)
358    #endif /* ALLOW_BULK_LARGEYEAGER04 */
359               psixh  = -psim_fac*htol*stable
360         &            + (exf_one-stable)
361       &              *exf_two*LOG( exf_half*(exf_one+xsq) )       &              *exf_two*LOG( exf_half*(exf_one+xsq) )
362    
363  C     Shift wind speed using old coefficient  C     Shift wind speed using old coefficient
364              usn    = ws/( exf_one + rdn*(zwln-psimh)/karman )  #ifdef ALLOW_BULK_LARGEYEAGER04
365              usm    = MAX(usn, umin)  C--   Large&Yeager04:
366               usn    = wspeed(i,j,bi,bj)
367         &           /( exf_one + rdn(i,j)*(zwln-psimh)/karman )
368    #else
369    C--   Large&Pond1981:
370               usn   = sh(i,j,bi,bj)/(exf_one - rdn(i,j)/karman*psimh)
371    #endif /* ALLOW_BULK_LARGEYEAGER04 */
372               usm    = MAX(usn, umin)
373    
374  C-    Update the 10m, neutral stability transfer coefficients  C-    Update the 10m, neutral stability transfer coefficients
375  c           tmpbulk= exf_BulkCdn(usm)  c          tmpbulk= exf_BulkCdn(usm)
376              tmpbulk= cdrag_1/usm + cdrag_2 + cdrag_3*usm             tmpbulk= cdrag_1/usm + cdrag_2 + cdrag_3*usm
377              rdn    = SQRT(tmpbulk)             rdn(i,j) = SQRT(tmpbulk)
378  c           rhn    = exf_BulkRhn(stable)  c          rhn    = exf_BulkRhn(stable)
379              rhn    = (exf_one-stable)*cstanton_1 + stable*cstanton_2             rhn    = (exf_one-stable)*cstanton_1 + stable*cstanton_2
380    
381  C     Shift all coefficients to the measurement height and stability.  C     Shift all coefficients to the measurement height and stability.
382              rd     = rdn/( exf_one + rdn*(zwln-psimh)/karman )  #ifdef ALLOW_BULK_LARGEYEAGER04
383              rh     = rhn/( exf_one + rhn*(ztln-psixh)/karman )             rd(i,j)= rdn(i,j)/( exf_one + rdn(i,j)*(zwln-psimh)/karman )
384              re     = ren/( exf_one + ren*(ztln-psixh)/karman )  #else
385               rd(i,j)= rdn(i,j)/( exf_one - rdn(i,j)/karman*psimh )
386    #endif /* ALLOW_BULK_LARGEYEAGER04 */
387               rh(i,j)= rhn/( exf_one + rhn*(ztln-psixh)/karman )
388               re(i,j)= ren/( exf_one + ren*(ztln-psixh)/karman )
389    
390  C     Update ustar, tstar, qstar using updated, shifted coefficients.  C     Update ustar, tstar, qstar using updated, shifted coefficients.
391              ustar  = rd*wsm             ustar(i,j)  = rd(i,j)*sh(i,j,bi,bj)
392              qstar  = re*delq             qstar(i,j)  = re(i,j)*delq(i,j)
393              tstar  = rh*deltap             tstar(i,j)  = rh(i,j)*deltap(i,j)
394             ENDDO            ENDIF
395    C     end i/j-loops
396             tau     = atmrho*rd*ws           ENDDO
397            ENDDO
398             evapLoc = -tau*qstar  C     end iteration loop
399             hlLocal = -lath*evapLoc         ENDDO
400             hsLocal = atmcp*tau*tstar         DO j=jMin,jMax
401  c          ustress = tau*rd*UwindSpeed          DO i=iMin,iMax
402  c          vstress = tau*rd*VwindSpeed           IF ( iceFlag(i,j) .AND. atemp(i,j,bi,bj).NE.0. _d 0 ) THEN
403              tau     = atmrho*rd(i,j)*wspeed(i,j,bi,bj)
404              evapLoc(i,j)  = -tau*qstar(i,j)
405              hlLocal       = -lath*evapLoc(i,j)
406              hsLocal       = atmcp*tau*tstar(i,j)
407    c         ustress = tau*rd(i,j)*UwindSpeed
408    c         vstress = tau*rd(i,j)*VwindSpeed
409    
410  C---  surf.Temp derivative of turbulent Fluxes  C---  surf.Temp derivative of turbulent Fluxes
411             dEvdT  =  (tau*re)*ssq*qsat_exp/Ts2  C     complete computation of dEvdT
412             dflhdT = -lath*dEvdT            dEvdT(i,j)    = (tau*re(i,j))*dEvdT(i,j)
413             dfshdT = -atmcp*tau*rh            dflhdT        = -lath*dEvdT(i,j)
414              dfshdT        = -atmcp*tau*rh(i,j)
415            ELSE  C--   Update total derivative with respect to surface temperature
416              dFlxdT(i,j)   = dFlxdT(i,j)   + dfshdT  + dflhdT
417    C--   Update net downward radiation excluding shortwave
418              flxExcSw(i,j) = flxExcSw(i,j) + hsLocal + hlLocal
419    
420             ENDIF
421            ENDDO
422           ENDDO
423          ELSE
424    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
425  C--   Compute the turbulent surface fluxes using fixed transfert Coeffs  C--   Compute the turbulent surface fluxes using fixed transfert Coeffs
426  C      with no stability dependence ( useStabilityFct_overIce = false )  C     with no stability dependence ( useStabilityFct_overIce = false )
427           DO j=jMin,jMax
428             evapLoc =      -atmrho*exf_iceCe*wsm*delq          DO i=iMin,iMax
429             hlLocal = -lath*evapLoc           IF ( iceFlag(i,j) .AND. atemp(i,j,bi,bj).NE.0. _d 0 ) THEN
430             hsLocal = atmcp*atmrho*exf_iceCh*wsm*deltap            wsm           = sh(i,j,bi,bj)
431              tau           = atmrho*exf_iceCe*wsm
432              evapLoc(i,j)  = -tau*delq(i,j)
433              hlLocal       = -lath*evapLoc(i,j)
434              hsLocal       = atmcp*atmrho*exf_iceCh*wsm*deltap(i,j)
435    #ifdef ALLOW_DBUG_THSICE
436              IF ( dBug(i,j,bi,bj) ) WRITE(stdUnit,'(A,4F12.6)')
437         &      'ThSI_GET_EXF: wsm,hl,hs,Lw=',
438         &                     wsm,hlLocal,hsLocal,flxExcSw(i,j)
439    #endif
440  C---  surf.Temp derivative of turbulent Fluxes  C---  surf.Temp derivative of turbulent Fluxes
441             dEvdT  = (atmrho*exf_iceCe*wsm)*(ssq*qsat_exp/Ts2)  C     complete computation of dEvdT
442             dflhdT = -lath*dEvdT            dEvdT(i,j)    = tau*dEvdT(i,j)
443             dfshdT = -atmcp*atmrho*exf_iceCh*wsm            dflhdT        = -lath*dEvdT(i,j)
444              dfshdT        = -atmcp*atmrho*exf_iceCh*wsm
445            ENDIF  C--   Update total derivative with respect to surface temperature
446  C--- Upward long wave radiation            dFlxdT(i,j)   = dFlxdT(i,j)   + dfshdT  + dflhdT
447            IF ( iceornot.EQ.0 ) THEN  C--   Update net downward radiation excluding shortwave
448              emiss = ocean_emissivity            flxExcSw(i,j) = flxExcSw(i,j) + hsLocal + hlLocal
449            ELSEIF (iceornot.EQ.2) THEN  #ifdef ALLOW_DBUG_THSICE
450              emiss = snow_emissivity            IF ( dBug(i,j,bi,bj) ) WRITE(stdUnit,'(A,4F12.6)')
451            ELSE       &      'ThSI_GET_EXF: flx,dFlxdT,evap,dEvdT',
452              emiss = ice_emissivity       &       flxExcSw(i,j), dFlxdT(i,j), evapLoc(i,j),dEvdT(i,j)
           ENDIF  
           flwup    = emiss*stefanBoltzmann*Ts2*Ts2  
           dflwupdT = emiss*stefanBoltzmann*Ts2*Tsf * 4. _d 0  
   
 C--   Total derivative with respect to surface temperature  
           df0dT = -dflwupdT+dfshdT+dflhdT  
   
 #ifdef ALLOW_DOWNWARD_RADIATION  
 c         flwNet_dwn  =       lwdown(i,j,bi,bj) - flwup  
 C-    assume long-wave albedo = 1 - emissivity :  
           flwNet_dwn  = emiss*lwdown(i,j,bi,bj) - flwup  
 #else  
       STOP 'ABNORMAL END: S/R THSICE_GET_EXF: DOWNWARD_RADIATION undef'  
453  #endif  #endif
454            flxExceptSw = flwNet_dwn + hsLocal + hlLocal           ENDIF
455            ENDDO
456          ELSE         ENDDO
457            flxExceptSw = 0. _d 0  C     endif useStabilityFct_overIce
458            df0dT   = 0. _d 0        ENDIF
           evapLoc = 0. _d 0  
           dEvdT   = 0. _d 0  
         ENDIF  
   
459  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
460          DO j=jMin,jMax
461           DO i=iMin,iMax
462            IF ( iceFlag(i,j) .AND. atemp(i,j,bi,bj).LE.0. _d 0 ) THEN
463    C--   in case atemp is zero:
464             flxExcSw(i,j) = 0. _d 0
465             dFlxdT  (i,j) = 0. _d 0
466             evapLoc (i,j) = 0. _d 0
467             dEvdT   (i,j) = 0. _d 0
468            ENDIF
469           ENDDO
470          ENDDO
471    
472  #else /* ALLOW_ATM_TEMP */  #else /* ALLOW_ATM_TEMP */
473        STOP 'ABNORMAL END: S/R THSICE_GET_EXF: ATM_TEMP undef'        STOP 'ABNORMAL END: S/R THSICE_GET_EXF: ATM_TEMP undef'
Legend:
Removed from v.1.11  
changed lines
  Added in v.1.18
  ViewVC Help
Powered by ViewVC 1.1.22