/[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.17 by heimbach, Thu Oct 21 02:10:34 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  #endif
148    
149  C     == external functions ==  C     == external functions ==
# Line 132  c     external  exf_BulkRhn Line 158  c     external  exf_BulkRhn
158  C     == end of interface ==  C     == end of interface ==
159    
160  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
161            act1 = bi - myBxLo(myThid)  c         act1 = bi - myBxLo(myThid)
162            max1 = myBxHi(myThid) - myBxLo(myThid) + 1  c         max1 = myBxHi(myThid) - myBxLo(myThid) + 1
163            act2 = bj - myByLo(myThid)  c         act2 = bj - myByLo(myThid)
164            max2 = myByHi(myThid) - myByLo(myThid) + 1  c         max2 = myByHi(myThid) - myByLo(myThid) + 1
165            act3 = myThid - 1  c         act3 = myThid - 1
166            max3 = nTx*nTy  c         max3 = nTx*nTy
167            act4 = ikey_dynamics - 1  c         act4 = ikey_dynamics - 1
168    c         ikey_1 = i
169            ikey_1 = i  c    &         + sNx*(j-1)
170       &         + sNx*(j-1)  c    &         + sNx*sNy*act1
171       &         + sNx*sNy*act1  c    &         + sNx*sNy*max1*act2
172       &         + sNx*sNy*max1*act2  c    &         + sNx*sNy*max1*max2*act3
173       &         + sNx*sNy*max1*max2*act3  c    &         + sNx*sNy*max1*max2*max3*act4
      &         + sNx*sNy*max1*max2*max3*act4  
174  #endif  #endif
175    
176    C--   Set surface parameters :
177          zwln = LOG(hu/zref)
178          ztln = LOG(ht/zref)
179          czol = hu*karman*gravity_mks
180          ren  = cDalton
181    C     more abbreviations
182          lath     = flamb+flami
183          qsat_fac = cvapor_fac_ice
184          qsat_exp = cvapor_exp_ice
185    
186    C     initialisation of local arrays
187          DO j = 1-Oly,sNy+Oly
188           DO i = 1-Olx,sNx+Olx
189            tstar(i,j)  = 0. _d 0
190            qstar(i,j)  = 0. _d 0
191            ustar(i,j)  = 0. _d 0
192            rdn(i,j)    = 0. _d 0
193            rd(i,j)     = 0. _d 0
194            rh(i,j)     = 0. _d 0
195            re(i,j)     = 0. _d 0
196            delq(i,j)   = 0. _d 0
197            deltap(i,j) = 0. _d 0
198            ssq(i,j)    = 0. _d 0
199           ENDDO
200          ENDDO
201    C
202          DO j=jMin,jMax
203           DO i=iMin,iMax
204    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7---+----
205    
206    C--   Use atmospheric state to compute surface fluxes.
207            IF ( iceFlag(i,j) .AND. atemp(i,j,bi,bj).NE.0. _d 0 ) THEN
208            
209             IF ( hSnow(i,j).GT.3. _d -1 ) THEN
210              emiss = snow_emissivity
211             ELSE
212              emiss = ice_emissivity
213             ENDIF
214  C     copy a few variables to names used in bulkf_formula_lay  C     copy a few variables to names used in bulkf_formula_lay
215        Tsf    =  tsfCel+cen2kel           Tsf         = tsfCel(i,j)+cen2kel
216        Ts2    = Tsf*Tsf           Ts2         = Tsf*Tsf
217  C     wind speed  C     wind speed
       ws     = wspeed(i,j,bi,bj)  
218  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
219  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
220  #endif  #endif
221        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  
   
222  C--   air - surface difference of temperature & humidity  C--   air - surface difference of temperature & humidity
223  c         tmpbulk= exf_BulkqSat(Tsf)  c        tmpbulk= exf_BulkqSat(Tsf)
224  c         ssq    = saltsat*tmpbulk/atmrho  c        ssq(i,j)    = saltsat*tmpbulk/atmrho
225            tmpbulk= qsat_fac*EXP(-qsat_exp/Tsf)           tmpbulk     = qsat_fac*EXP(-qsat_exp/Tsf)
226            ssq    = tmpbulk/atmrho           ssq(i,j)    = tmpbulk/atmrho
227            deltap = atemp(i,j,bi,bj) + gamma_blk*ht - Tsf           deltap(i,j) = atemp(i,j,bi,bj) + gamma_blk*ht - Tsf
228            delq   = aqh(i,j,bi,bj) - ssq           delq(i,j)   = aqh(i,j,bi,bj) - ssq(i,j)
229    C     Do the part of the output variables that do not depend
230    C     on the ice here to save a few re-computations
231    C     This is not yet dEvdT, but just a cheap way to save a 2D-field
232    C     for ssq and recomputing Ts2 lateron
233             dEvdT(i,j)  = ssq(i,j)*qsat_exp/Ts2
234             flwup       = emiss*stefanBoltzmann*Ts2*Ts2
235             dflwupdT    = emiss*stefanBoltzmann*Ts2*Tsf * 4. _d 0
236    #ifdef ALLOW_DOWNWARD_RADIATION
237    c        flwNet_dwn  =       lwdown(i,j,bi,bj) - flwup
238    C-    assume long-wave albedo = 1 - emissivity :
239             flwNet_dwn  = emiss*lwdown(i,j,bi,bj) - flwup
240    #else
241             STOP
242         &     'ABNORMAL END: S/R THSICE_GET_EXF: DOWNWARD_RADIATION undef'
243    #endif
244    C--   This is not yet the total derivative with respect to surface temperature
245             dFlxdT(i,j)   = -dflwupdT
246    C--   This is not yet the Net downward radiation excluding shortwave
247             flxExcSw(i,j) = flwNet_dwn
248            ENDIF
249           ENDDO
250          ENDDO
251    
252            IF ( useStabilityFct_overIce ) THEN        IF ( useStabilityFct_overIce ) THEN
253           DO j=jMin,jMax
254            DO i=iMin,iMax
255             IF ( iceFlag(i,j) .AND. atemp(i,j,bi,bj).NE.0. _d 0 ) THEN
256  C--   Compute the turbulent surface fluxes (function of stability).  C--   Compute the turbulent surface fluxes (function of stability).
257    
 C--   Set surface parameters :  
            zwln = LOG(hu/zref)  
            ztln = LOG(ht/zref)  
            czol = hu*karman*gravity_mks  
258    
259  C             Initial guess: z/l=0.0; hu=ht=hq=z  C             Initial guess: z/l=0.0; hu=ht=hq=z
260  C             Iterations:    converge on z/l and hence the fluxes.  C             Iterations:    converge on z/l and hence the fluxes.
261    
262             t0     = atemp(i,j,bi,bj)*            t0         = atemp(i,j,bi,bj)*
263       &            (exf_one + humid_fac*aqh(i,j,bi,bj))       &         (exf_one + humid_fac*aqh(i,j,bi,bj))
264             stable = exf_half + SIGN(exf_half, deltap)            stable     = exf_half + SIGN(exf_half, deltap(i,j))
265  c          tmpbulk= exf_BulkCdn(sh(i,j,bi,bj))  c         tmpbulk    = exf_BulkCdn(sh(i,j,bi,bj))
266             tmpbulk= cdrag_1/wsm + cdrag_2 + cdrag_3*wsm            wsm        = sh(i,j,bi,bj)
267             rdn    = SQRT(tmpbulk)            tmpbulk    = cdrag_1/wsm + cdrag_2 + cdrag_3*wsm
268              IF (tmpbulk.NE.0.) THEN
269               rdn(i,j)   = SQRT(tmpbulk)
270              ELSE
271               rdn(i,j)   = 0. _d 0
272              ENDIF
273  C--  initial guess for exchange other coefficients:  C--  initial guess for exchange other coefficients:
274  c          rhn    = exf_BulkRhn(stable)  c         rhn        = exf_BulkRhn(stable)
275             rhn    = (exf_one-stable)*cstanton_1 + stable*cstanton_2            rhn        = (exf_one-stable)*cstanton_1 + stable*cstanton_2
            ren    = cDalton  
276  C--  calculate turbulent scales  C--  calculate turbulent scales
277             ustar  = rdn*wsm            ustar(i,j) = rdn(i,j)*wsm
278             tstar  = rhn*deltap            tstar(i,j) = rhn*deltap(i,j)
279             qstar  = ren*delq            qstar(i,j) = ren*delq(i,j)
280             ENDIF
281             DO iter = 1,niter_bulk          ENDDO
282           ENDDO
283    C     start iteration
284           DO iter = 1,niter_bulk
285            DO j=jMin,jMax
286             DO i=iMin,iMax
287              IF ( iceFlag(i,j) .AND. atemp(i,j,bi,bj).NE.0. _d 0 ) THEN
288    
289  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
290                     ikey_2 = iter  c                  ikey_2 = iter
291       &                  + niter_bulk*(i-1)  c    &                  + niter_bulk*(i-1)
292       &                  + niter_bulk*sNx*(j-1)  c    &                  + niter_bulk*sNx*(j-1)
293       &                  + niter_bulk*sNx*sNy*act1  c    &                  + niter_bulk*sNx*sNy*act1
294       &                  + niter_bulk*sNx*sNy*max1*act2  c    &                  + niter_bulk*sNx*sNy*max1*act2
295       &                  + niter_bulk*sNx*sNy*max1*max2*act3  c    &                  + niter_bulk*sNx*sNy*max1*max2*act3
296       &                  + niter_bulk*sNx*sNy*max1*max2*max3*act4  c    &                  + niter_bulk*sNx*sNy*max1*max2*max3*act4
297    cCADJ STORE rdn    = comlev1_exf_2, key = ikey_2
298  CADJ STORE rdn    = comlev1_exf_2, key = ikey_2  cCADJ STORE ustar  = comlev1_exf_2, key = ikey_2
299  CADJ STORE ustar  = comlev1_exf_2, key = ikey_2  cCADJ STORE qstar  = comlev1_exf_2, key = ikey_2
300  CADJ STORE qstar  = comlev1_exf_2, key = ikey_2  cCADJ STORE tstar  = comlev1_exf_2, key = ikey_2
301  CADJ STORE tstar  = comlev1_exf_2, key = ikey_2  cCADJ STORE sh(i,j,bi,bj)     = comlev1_exf_2, key = ikey_2
 CADJ STORE sh(i,j,bi,bj)     = comlev1_exf_2, key = ikey_2  
302  #endif  #endif
303    
304              huol   = (tstar/t0 +             t0     = atemp(i,j,bi,bj)*
305       &                qstar/(exf_one/humid_fac+aqh(i,j,bi,bj))       &          (exf_one + humid_fac*aqh(i,j,bi,bj))
306       &               )*czol/(ustar*ustar)             huol   = (tstar(i,j)/t0 +
307  C-    Large&Pond_1981 code (zolmin default = -100):       &               qstar(i,j)/(exf_one/humid_fac+aqh(i,j,bi,bj))
308  c           huol   = MAX(huol,zolmin)       &              )*czol/(ustar(i,j)*ustar(i,j))
309    #ifdef ALLOW_BULK_LARGEYEAGER04
310  C-    Large&Yeager_2004 code:  C-    Large&Yeager_2004 code:
311              huol   = MIN( MAX(-10. _d 0,huol), 10. _d 0 )             huol   = MIN( MAX(-10. _d 0,huol), 10. _d 0 )
312              htol   = huol*ht/hu  #else
313              hqol   = huol*hq/hu  C-    Large&Pond_1981 code (zolmin default = -100):
314              stable = exf_half + SIGN(exf_half, huol)             huol   = MAX(huol,zolmin)
315    #endif /* ALLOW_BULK_LARGEYEAGER04 */
316               htol   = huol*ht/hu
317               hqol   = huol*hq/hu
318               stable = exf_half + SIGN(exf_half, huol)
319    
320  C     Evaluate all stability functions assuming hq = ht.  C     Evaluate all stability functions assuming hq = ht.
321  C-    Large&Pond_1981 code:  #ifdef ALLOW_BULK_LARGEYEAGER04
 c           xsq    = MAX(SQRT(ABS(exf_one - huol*16. _d 0)),exf_one)  
322  C-    Large&Yeager_2004 code:  C-    Large&Yeager_2004 code:
323              xsq    = SQRT( ABS(exf_one - huol*16. _d 0) )             xsq    = SQRT( ABS(exf_one - huol*16. _d 0) )
324              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 )  
325  C-    Large&Pond_1981 code:  C-    Large&Pond_1981 code:
326  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)
327  C-    Large&Yeager_2004 code:  #endif /* ALLOW_BULK_LARGEYEAGER04 */
328              xsq    = SQRT( ABS(exf_one - htol*16. _d 0) )             x      = SQRT(xsq)
329              psixh  = -psim_fac*htol*stable             psimh  = -psim_fac*huol*stable
330       &             + (exf_one-stable)       &             + (exf_one-stable)
331         &             *( LOG( (exf_one + exf_two*x + xsq)
332         &                    *(exf_one+xsq)*0.125 _d 0 )
333         &                -exf_two*ATAN(x) + exf_half*pi )
334    #ifdef ALLOW_BULK_LARGEYEAGER04
335    C-    Large&Yeager_2004 code:
336               xsq    = SQRT( ABS(exf_one - htol*16. _d 0) )
337    #else
338    C-    Large&Pond_1981 code:
339               xsq    = MAX(SQRT(ABS(exf_one - htol*16. _d 0)),exf_one)
340    #endif /* ALLOW_BULK_LARGEYEAGER04 */
341               psixh  = -psim_fac*htol*stable
342         &            + (exf_one-stable)
343       &              *exf_two*LOG( exf_half*(exf_one+xsq) )       &              *exf_two*LOG( exf_half*(exf_one+xsq) )
344    
345  C     Shift wind speed using old coefficient  C     Shift wind speed using old coefficient
346              usn    = ws/( exf_one + rdn*(zwln-psimh)/karman )  #ifdef ALLOW_BULK_LARGEYEAGER04
347              usm    = MAX(usn, umin)  C--   Large&Yeager04:
348               usn    = wspeed(i,j,bi,bj)
349         &           /( exf_one + rdn(i,j)*(zwln-psimh)/karman )
350    #else
351    C--   Large&Pond1981:
352               usn   = sh(i,j,bi,bj)/(exf_one - rdn(i,j)/karman*psimh)
353    #endif /* ALLOW_BULK_LARGEYEAGER04 */
354               usm    = MAX(usn, umin)
355    
356  C-    Update the 10m, neutral stability transfer coefficients  C-    Update the 10m, neutral stability transfer coefficients
357  c           tmpbulk= exf_BulkCdn(usm)  c          tmpbulk= exf_BulkCdn(usm)
358              tmpbulk= cdrag_1/usm + cdrag_2 + cdrag_3*usm             tmpbulk= cdrag_1/usm + cdrag_2 + cdrag_3*usm
359              rdn    = SQRT(tmpbulk)             rdn(i,j) = SQRT(tmpbulk)
360  c           rhn    = exf_BulkRhn(stable)  c          rhn    = exf_BulkRhn(stable)
361              rhn    = (exf_one-stable)*cstanton_1 + stable*cstanton_2             rhn    = (exf_one-stable)*cstanton_1 + stable*cstanton_2
362    
363  C     Shift all coefficients to the measurement height and stability.  C     Shift all coefficients to the measurement height and stability.
364              rd     = rdn/( exf_one + rdn*(zwln-psimh)/karman )  #ifdef ALLOW_BULK_LARGEYEAGER04
365              rh     = rhn/( exf_one + rhn*(ztln-psixh)/karman )             rd(i,j)= rdn(i,j)/( exf_one + rdn(i,j)*(zwln-psimh)/karman )
366              re     = ren/( exf_one + ren*(ztln-psixh)/karman )  #else
367               rd(i,j)= rdn(i,j)/( exf_one - rdn(i,j)/karman*psimh )
368    #endif /* ALLOW_BULK_LARGEYEAGER04 */
369               rh(i,j)= rhn/( exf_one + rhn*(ztln-psixh)/karman )
370               re(i,j)= ren/( exf_one + ren*(ztln-psixh)/karman )
371    
372  C     Update ustar, tstar, qstar using updated, shifted coefficients.  C     Update ustar, tstar, qstar using updated, shifted coefficients.
373              ustar  = rd*wsm             ustar(i,j)  = rd(i,j)*sh(i,j,bi,bj)
374              qstar  = re*delq             qstar(i,j)  = re(i,j)*delq(i,j)
375              tstar  = rh*deltap             tstar(i,j)  = rh(i,j)*deltap(i,j)
376             ENDDO            ENDIF
377    C     end i/j-loops
378             tau     = atmrho*rd*ws           ENDDO
379            ENDDO
380             evapLoc = -tau*qstar  C     end iteration loop
381             hlLocal = -lath*evapLoc         ENDDO
382             hsLocal = atmcp*tau*tstar         DO j=jMin,jMax
383  c          ustress = tau*rd*UwindSpeed          DO i=iMin,iMax
384  c          vstress = tau*rd*VwindSpeed           IF ( iceFlag(i,j) .AND. atemp(i,j,bi,bj).NE.0. _d 0 ) THEN
385              
386              tau     = atmrho*rd(i,j)*wspeed(i,j,bi,bj)
387              
388              evapLoc(i,j)  = -tau*qstar(i,j)
389              hlLocal       = -lath*evapLoc(i,j)
390              hsLocal       = atmcp*tau*tstar(i,j)
391    c         ustress = tau*rd(i,j)*UwindSpeed
392    c         vstress = tau*rd(i,j)*VwindSpeed
393    
394  C---  surf.Temp derivative of turbulent Fluxes  C---  surf.Temp derivative of turbulent Fluxes
395             dEvdT  =  (tau*re)*ssq*qsat_exp/Ts2  C     complete computation of dEvdT
396             dflhdT = -lath*dEvdT            dEvdT(i,j)    = (tau*re(i,j))*dEvdT(i,j)
397             dfshdT = -atmcp*tau*rh            dflhdT        = -lath*dEvdT(i,j)
398              dfshdT        = -atmcp*tau*rh(i,j)
399            ELSE  C--   Update total derivative with respect to surface temperature
400              dFlxdT(i,j)   = dFlxdT(i,j)   + dfshdT  + dflhdT
401    C--   Update net downward radiation excluding shortwave
402              flxExcSw(i,j) = flxExcSw(i,j) + hsLocal + hlLocal
403    
404             ENDIF
405            ENDDO
406           ENDDO
407          ELSE
408  C--   Compute the turbulent surface fluxes using fixed transfert Coeffs  C--   Compute the turbulent surface fluxes using fixed transfert Coeffs
409  C      with no stability dependence ( useStabilityFct_overIce = false )  C     with no stability dependence ( useStabilityFct_overIce = false )
410           DO j=jMin,jMax
411             evapLoc =      -atmrho*exf_iceCe*wsm*delq          DO i=iMin,iMax
412             hlLocal = -lath*evapLoc           IF ( iceFlag(i,j) .AND. atemp(i,j,bi,bj).NE.0. _d 0 ) THEN
413             hsLocal = atmcp*atmrho*exf_iceCh*wsm*deltap            
414              wsm           = sh(i,j,bi,bj)
415              tau           = atmrho*exf_iceCe*wsm
416              evapLoc(i,j)  = -tau*delq(i,j)
417              hlLocal       = -lath*evapLoc(i,j)
418              hsLocal       = atmcp*atmrho*exf_iceCh*wsm*deltap(i,j)
419  C---  surf.Temp derivative of turbulent Fluxes  C---  surf.Temp derivative of turbulent Fluxes
420             dEvdT  = (atmrho*exf_iceCe*wsm)*(ssq*qsat_exp/Ts2)  C     complete computation of dEvdT
421             dflhdT = -lath*dEvdT            dEvdT(i,j)    = tau*dEvdT(i,j)
422             dfshdT = -atmcp*atmrho*exf_iceCh*wsm            dflhdT        = -lath*dEvdT(i,j)
423              dfshdT        = -atmcp*atmrho*exf_iceCh*wsm
424            ENDIF  C--   Update total derivative with respect to surface temperature
425  C--- Upward long wave radiation            dFlxdT(i,j)   = dFlxdT(i,j)   + dfshdT  + dflhdT
426            IF ( iceornot.EQ.0 ) THEN  C--   Update net downward radiation excluding shortwave
427              emiss = ocean_emissivity            flxExcSw(i,j) = flxExcSw(i,j) + hsLocal + hlLocal
428            ELSEIF (iceornot.EQ.2) THEN            
429              emiss = snow_emissivity           ENDIF
430            ELSE          ENDDO
431              emiss = ice_emissivity         ENDDO
432            ENDIF  C     endif useStabilityFct_overIce
433            flwup    = emiss*stefanBoltzmann*Ts2*Ts2        ENDIF
434            dflwupdT = emiss*stefanBoltzmann*Ts2*Tsf * 4. _d 0        DO j=jMin,jMax
435           DO i=iMin,iMax
436  C--   Total derivative with respect to surface temperature          IF ( iceFlag(i,j) .AND. atemp(i,j,bi,bj).LE.0. _d 0 ) THEN
437            df0dT = -dflwupdT+dfshdT+dflhdT  C--   in case atemp is zero:
438             flxExcSw(i,j) = 0. _d 0
439  #ifdef ALLOW_DOWNWARD_RADIATION           dFlxdT  (i,j) = 0. _d 0
440  c         flwNet_dwn  =       lwdown(i,j,bi,bj) - flwup           evapLoc (i,j) = 0. _d 0
441  C-    assume long-wave albedo = 1 - emissivity :           dEvdT   (i,j) = 0. _d 0
           flwNet_dwn  = emiss*lwdown(i,j,bi,bj) - flwup  
 #else  
       STOP 'ABNORMAL END: S/R THSICE_GET_EXF: DOWNWARD_RADIATION undef'  
 #endif  
           flxExceptSw = flwNet_dwn + hsLocal + hlLocal  
   
         ELSE  
           flxExceptSw = 0. _d 0  
           df0dT   = 0. _d 0  
           evapLoc = 0. _d 0  
           dEvdT   = 0. _d 0  
442          ENDIF          ENDIF
443    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7---+----
444  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|         ENDDO
445          ENDDO
446    
447  #else /* ALLOW_ATM_TEMP */  #else /* ALLOW_ATM_TEMP */
448        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.17
  ViewVC Help
Powered by ViewVC 1.1.22