/[MITgcm]/MITgcm/pkg/exf/exf_bulkformulae.F
ViewVC logotype

Diff of /MITgcm/pkg/exf/exf_bulkformulae.F

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

revision 1.22 by mlosch, Fri May 21 07:50:51 2010 UTC revision 1.23 by mlosch, Fri May 21 10:08:44 2010 UTC
# Line 3  C $Name$ Line 3  C $Name$
3    
4  #include "EXF_OPTIONS.h"  #include "EXF_OPTIONS.h"
5    
6        subroutine exf_bulkformulae(mytime, myiter, mythid)  CBOP
7    C     !ROUTINE: EXF_BULKFORMULAE
8  c     ==================================================================  C     !INTERFACE:
9  c     SUBROUTINE exf_bulkformulae        SUBROUTINE EXF_BULKFORMULAE( myTime, myIter, myThid )
10  c     ==================================================================  
11  c  C     !DESCRIPTION: \bv
12  c     o Use bulk formulae to estimate turbulent and/or radiative  C     *==========================================================*
13  c       fluxes at the surface.  C     | SUBROUTINE EXF_BULKFORMULAE
14  c  C     | o Calculate bulk formula fluxes over open ocean
15  c     NOTES:  C     |   either following
16  c     ======  C     |   Large and Pond, 1981 & 1982
17  c  C     |   or (if defined ALLOW_BULK_LARGEYEAGER04)
18  c     See EXF_OPTIONS.h for a description of the various possible  C     |   Large and Yeager, 2004, NCAR/TN-460+STR.
19  c     ocean-model forcing configurations.  C     *==========================================================*
20  c  C     \ev
21  c     The bulk formulae of pkg/exf are not valid for sea-ice covered  C
22  c     oceans but they can be used in combination with a sea-ice model,  C
23  c     for example, pkg/seaice, to specify open water flux contributions.  C     NOTES:
24  c  C     ======
25  c     ==================================================================  C
26  c  C     See EXF_OPTIONS.h for a description of the various possible
27  c       The calculation of the bulk surface fluxes has been adapted from  C     ocean-model forcing configurations.
28  c       the NCOM model which uses the formulae given in Large and Pond  C
29  c       (1981 & 1982 )  C     The bulk formulae of pkg/exf are not valid for sea-ice covered
30  c  C     oceans but they can be used in combination with a sea-ice model,
31  c  C     for example, pkg/seaice, to specify open water flux contributions.
32  c       Header taken from NCOM version: ncom1.4.1  C
33  c       -----------------------------------------  C     ==================================================================
34  c  C
35  c       Following procedures and coefficients in Large and Pond  C     for Large and Pond, 1981 & 1982
36  c       (1981 ; 1982)  C
37  c  C     The calculation of the bulk surface fluxes has been adapted from
38  c       Output: Bulk estimates of the turbulent surface fluxes.  C     the NCOM model which uses the formulae given in Large and Pond
39  c       -------  C     (1981 & 1982 )
40  c  C
41  c       hs  - sensible heat flux  (W/m^2), into ocean  C
42  c       hl  - latent   heat flux  (W/m^2), into ocean  C     Header taken from NCOM version: ncom1.4.1
43  c  C     -----------------------------------------
44  c       Input:  C
45  c       ------  C     Following procedures and coefficients in Large and Pond
46  c  C     (1981 ; 1982)
47  c       us  - mean wind speed (m/s)     at height hu (m)  C
48  c       th  - mean air temperature (K)  at height ht (m)  C     Output: Bulk estimates of the turbulent surface fluxes.
49  c       qh  - mean air humidity (kg/kg) at height hq (m)  C     -------
50  c       sst - sea surface temperature (K)  C
51  c       tk0 - Kelvin temperature at 0 Celsius (K)  C     hs  - sensible heat flux  (W/m^2), into ocean
52  c  C     hl  - latent   heat flux  (W/m^2), into ocean
53  c       Assume 1) a neutral 10m drag coefficient =  C
54  c  C     Input:
55  c                 cdn = .0027/u10 + .000142 + .0000764 u10  C     ------
56  c  C
57  c              2) a neutral 10m stanton number =  C     us  - mean wind speed (m/s)     at height hu (m)
58  c  C     th  - mean air temperature (K)  at height ht (m)
59  c                 ctn = .0327 sqrt(cdn), unstable  C     qh  - mean air humidity (kg/kg) at height hq (m)
60  c                 ctn = .0180 sqrt(cdn), stable  C     sst - sea surface temperature (K)
61  c  C     tk0 - Kelvin temperature at 0 Celsius (K)
62  c              3) a neutral 10m dalton number =  C
63  c  C     Assume 1) a neutral 10m drag coefficient =
64  c                 cen = .0346 sqrt(cdn)  C
65  c  C               cdn = .0027/u10 + .000142 + .0000764 u10
66  c              4) the saturation humidity of air at  C
67  c  C            2) a neutral 10m stanton number =
68  c                 t(k) = exf_BulkqSat(t)  (kg/m^3)  C
69  c  C               ctn = .0327 sqrt(cdn), unstable
70  c       Note:  1) here, tstar = <wt>/u*, and qstar = <wq>/u*.  C               ctn = .0180 sqrt(cdn), stable
71  c              2) wind speeds should all be above a minimum speed,  C
72  c                 say 0.5 m/s.  C            3) a neutral 10m dalton number =
73  c              3) with optional iteration loop, niter=3, should suffice.  C
74  c              4) this version is for analyses inputs with hu = 10m and  C               cen = .0346 sqrt(cdn)
75  c                 ht = hq.  C
76  c              5) sst enters in Celsius.  C            4) the saturation humidity of air at
77  c  C
78  c     ==================================================================  C               t(k) = exf_BulkqSat(t)  (kg/m^3)
79  c  C
80  c       started: Christian Eckert eckert@mit.edu 27-Aug-1999  C     Note:  1) here, tstar = <wt>/u*, and qstar = <wq>/u*.
81  c  C            2) wind speeds should all be above a minimum speed,
82  c       changed: Christian Eckert eckert@mit.edu 14-Jan-2000  C               say 0.5 m/s.
83  c              - restructured the original version in order to have a  C            3) with optional iteration loop, niter=3, should suffice.
84  c                better interface to the MITgcmUV.  C            4) this version is for analyses inputs with hu = 10m and
85  c  C               ht = hq.
86  c              Christian Eckert eckert@mit.edu  12-Feb-2000  C            5) sst enters in Celsius.
87  c              - Changed Routine names (package prefix: exf_)  C
88  c  C     ==================================================================
89  c              Patrick Heimbach, heimbach@mit.edu  04-May-2000  C
90  c              - changed the handling of precip and sflux with respect  C       started: Christian Eckert eckert@mit.edu 27-Aug-1999
91  c                to CPP options ALLOW_BULKFORMULAE and ALLOW_ATM_TEMP  C
92  c              - included some CPP flags ALLOW_BULKFORMULAE to make  C     changed: Christian Eckert eckert@mit.edu 14-Jan-2000
93  c                sure ALLOW_ATM_TEMP, ALLOW_ATM_WIND are used only in  C            - restructured the original version in order to have a
94  c                conjunction with defined ALLOW_BULKFORMULAE  C              better interface to the MITgcmUV.
95  c              - statement functions discarded  C
96  c  C            Christian Eckert eckert@mit.edu  12-Feb-2000
97  c              Ralf.Giering@FastOpt.de 25-Mai-2000  C            - Changed Routine names (package prefix: exf_)
98  c              - total rewrite using new subroutines  C
99  c  C            Patrick Heimbach, heimbach@mit.edu  04-May-2000
100  c              Detlef Stammer: include river run-off. Nov. 21, 2001  C            - changed the handling of precip and sflux with respect
101  c  C              to CPP options ALLOW_BULKFORMULAE and ALLOW_ATM_TEMP
102  c              heimbach@mit.edu, 10-Jan-2002  C            - included some CPP flags ALLOW_BULKFORMULAE to make
103  c              - changes to enable field swapping  C              sure ALLOW_ATM_TEMP, ALLOW_ATM_WIND are used only in
104  c  C              conjunction with defined ALLOW_BULKFORMULAE
105  c       mods for pkg/seaice: menemenlis@jpl.nasa.gov 20-Dec-2002  C            - statement functions discarded
106  c  C
107  c     ==================================================================  C            Ralf.Giering@FastOpt.de 25-Mai-2000
108  c     SUBROUTINE exf_bulkformulae  C            - total rewrite using new subroutines
109  c     ==================================================================  C
110    C            Detlef Stammer: include river run-off. Nov. 21, 2001
111        implicit none  C
112    C            heimbach@mit.edu, 10-Jan-2002
113  c     == global variables ==  C            - changes to enable field swapping
114    C
115    C     mods for pkg/seaice: menemenlis@jpl.nasa.gov 20-Dec-2002
116    C
117    C     martin.losch@awi.de: merged with exf_bulk_largeyeager04, 21-May-2010
118    C
119    C     ==================================================================
120    C
121    C     for Large and Yeager, 2004
122    C
123    C === Turbulent Fluxes :
124    C  * use the approach "B": shift coeff to height & stability of the
125    C    atmosphere state (instead of "C": shift temp & humid to the height
126    C    of wind, then shift the coeff to this height & stability of the atmos).
127    C  * similar to EXF (except over sea-ice) ; default parameter values
128    C    taken from Large & Yeager.
129    C  * assume that Qair & Tair inputs are from the same height (zq=zt)
130    C  * formulae in short:
131    C     wind stress = (ust,vst) = rhoA * Cd * Ws * (del.u,del.v)
132    C     Sensib Heat flux = fsha = rhoA * Ch * Ws * del.T * CpAir
133    C     Latent Heat flux = flha = rhoA * Ce * Ws * del.Q * Lvap
134    C                      = -Evap * Lvap
135    C   with Ws = wind speed = sqrt(del.u^2 +del.v^2) ;
136    C        del.T = Tair - Tsurf ; del.Q = Qair - Qsurf ;
137    C        Cd,Ch,Ce = drag coefficient, Stanton number and Dalton number
138    C              respectively [no-units], function of height & stability
139    
140    C     !USES:
141           IMPLICIT NONE
142    C     === Global variables ===
143  #include "EEPARAMS.h"  #include "EEPARAMS.h"
144  #include "SIZE.h"  #include "SIZE.h"
145  #include "PARAMS.h"  #include "PARAMS.h"
146  #include "DYNVARS.h"  #include "DYNVARS.h"
147  c #include "GRID.h"  #include "GRID.h"
148    
149  #include "EXF_PARAM.h"  #include "EXF_PARAM.h"
150  #include "EXF_FIELDS.h"  #include "EXF_FIELDS.h"
# Line 126  c #include "GRID.h" Line 154  c #include "GRID.h"
154  #include "tamc.h"  #include "tamc.h"
155  #endif  #endif
156    
157  c     == routine arguments ==  C     !INPUT/OUTPUT PARAMETERS:
158    C     input:
159        integer mythid  C     myTime  :: Current time in simulation
160        integer myiter  C     myIter  :: Current iteration number in simulation
161        _RL     mytime  C     myThid  :: My Thread Id number
162          _RL     myTime
163          INTEGER myIter
164          INTEGER myThid
165    C     output:
166    CEOP
167    
168  #ifdef ALLOW_BULKFORMULAE  #ifdef ALLOW_BULKFORMULAE
169    C     == external Functions
170    
171  c     == local variables ==  C     == Local variables ==
172    C     i,j      :: grid point indices
173        integer bi,bj  C     bi,bj    :: tile indices
174        integer i,j,k  C     ssq      :: saturation specific humidity [kg/kg] in eq. with Sea-Surface water
175          INTEGER i,j,bi,bj
176    
177        _RL czol        _RL czol
178        _RL zwln               ! = log(zwd/zref)        _RL Tsf                   ! surface temperature [K]
179        _RL ztln               ! = log(zth/zref)        _RL wsm                   ! limited wind speed [m/s] (> umin)
180        _RL windStress         ! surface wind-stress (@ grid cell center)        _RL t0                    ! virtual temperature [K]
181    C     these need to be 2D-arrays for vectorizing code
182        integer iter        _RL tstar (1:sNx,1:sNy)   ! turbulent temperature scale [K]
183        _RL     tmpbulk        _RL qstar (1:sNx,1:sNy)   ! turbulent humidity scale  [kg/kg]
184        _RL     hqol        _RL ustar (1:sNx,1:sNy)   ! friction velocity [m/s]
185        _RL     htol        _RL tau   (1:sNx,1:sNy)   ! surface stress coef = rhoA * Ws * sqrt(Cd)
186        _RL     huol        _RL rdn   (1:sNx,1:sNy)   ! neutral, zref (=10m) values of rd
187        _RL     psimh        _RL rd    (1:sNx,1:sNy)   ! = sqrt(Cd)          [-]
188        _RL     psixh        _RL delq  (1:sNx,1:sNy)   ! specific humidity difference [kg/kg]
       _RL     re                ! = Ce / sqrt(Cd)     [-]  
       _RL     rh                ! = Ch / sqrt(Cd)     [-]  
       _RL     ren, rhn          ! neutral, zref (=10m) values of re, rh  
       _RL     tsf  
       _RL     ssq  
       _RL     stable  
       _RL     t0  
       _RL     wsm                ! limited wind speed [m/s] (> umin)  
       _RL     uzn  
       _RL     shn  
       _RL     xsq  
       _RL     x  
       _RL     recip_rhoConstFresh  
 C     these need to 2D-arrays for vectorizing code  
       _RL tstar (1:sNx,1:sNy)  
       _RL qstar (1:sNx,1:sNy)  
       _RL ustar (1:sNx,1:sNy)  
       _RL tau   (1:sNx,1:sNy)  
       _RL rd    (1:sNx,1:sNy)  
       _RL rdn   (1:sNx,1:sNy)  
       _RL delq  (1:sNx,1:sNy)  
189        _RL deltap(1:sNx,1:sNy)        _RL deltap(1:sNx,1:sNy)
190    C
191          _RL dzTmp
192          _RL SSTtmp
193          _RL ssq
194          _RL re                    ! = Ce / sqrt(Cd)     [-]
195          _RL rh                    ! = Ch / sqrt(Cd)     [-]
196          _RL ren, rhn              ! neutral, zref (=10m) values of re, rh
197          _RL usn, usm
198          _RL stable                ! = 1 if stable ; = 0 if unstable
199          _RL huol                  ! stability parameter at zwd [-] (=z/Monin-Obuklov length)
200          _RL htol                  ! stability parameter at zth [-]
201          _RL hqol
202          _RL x                     ! stability function  [-]
203          _RL xsq                   ! = x^2               [-]
204          _RL psimh                 ! momentum stability function
205          _RL psixh                 ! latent & sensib. stability function
206          _RL zwln                  ! = log(zwd/zref)
207          _RL ztln                  ! = log(zth/zref)
208          _RL windStress            ! surface wind-stress (@ grid cell center)
209          _RL tmpbulk
210          _RL recip_rhoConstFresh
211          INTEGER ksrf, ksrfp1
212          INTEGER iter
213    C     solve4Stress :: if F, by-pass momentum turb. flux (wind-stress) calculation
214          LOGICAL solve4Stress
215    #ifdef ALLOW_ATM_WIND
216          PARAMETER ( solve4Stress = .TRUE. )
217    #endif
218    
219  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
220        integer ikey_1        integer ikey_1
221        integer ikey_2        integer ikey_2
222  #endif  #endif
223    
224  c     == external functions ==  #ifndef ALLOW_ATM_WIND
225    #ifdef ALLOW_BULK_LARGEYEAGER04
226        integer  ilnblnk        solve4Stress = wspeedfile .NE. ' '
227        external ilnblnk  #else
228          solve4Stress = .FALSE.
229  c     == end of interface ==  #endif
230    #endif
231    
232    C-- Set surface parameters :
233        zwln = LOG(hu/zref)        zwln = LOG(hu/zref)
234        ztln = LOG(ht/zref)        ztln = LOG(ht/zref)
235        czol = hu*karman*gravity_mks        czol = hu*karman*gravity_mks
236    C--   abbreviation
237        recip_rhoConstFresh = 1. _d 0/rhoConstFresh        recip_rhoConstFresh = 1. _d 0/rhoConstFresh
238  c--   Use atmospheric state to compute surface fluxes.        
   
239  c     Loop over tiles.  c     Loop over tiles.
240  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
241  C--   HPF directive to help TAMC  C--   HPF directive to help TAMC
242  CHPF$ INDEPENDENT  CHPF$ INDEPENDENT
243  #endif  #endif
244        do bj = mybylo(mythid),mybyhi(mythid)        DO bj = myByLo(myThid),myByHi(myThid)
245  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
246  C--    HPF directive to help TAMC  C--    HPF directive to help TAMC
247  CHPF$  INDEPENDENT  CHPF$  INDEPENDENT
248  #endif  #endif
249         do bi = mybxlo(mythid),mybxhi(mythid)         DO bi = myBxLo(myThid),myBxHi(myThid)
250          k = 1          ksrf   = 1
251          do j = 1,sny          ksrfp1 = 2
252           do i = 1,snx          DO j = 1,sNy
253             DO i = 1,sNx
254    
255  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
256            act1 = bi - myBxLo(myThid)            act1 = bi - myBxLo(myThid)
# Line 223  CHPF$  INDEPENDENT Line 269  CHPF$  INDEPENDENT
269       &         + sNx*sNy*max1*max2*max3*act4       &         + sNx*sNy*max1*max2*max3*act4
270  #endif  #endif
271    
 c--   Compute the turbulent surface fluxes.  
   
272  #ifdef ALLOW_ATM_TEMP  #ifdef ALLOW_ATM_TEMP
273    
 c             Initial guess: z/l=0.0; hu=ht=hq=z  
 c             Iterations:    converge on z/l and hence the fluxes.  
 c             t0     : virtual temperature (K)  
 c             ssq    : sea surface humidity (kg/kg)  
 c             deltap : potential temperature diff (K)  
   
274  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
275            deltap(i,j) = 0. _d 0            deltap(i,j) = 0. _d 0
276            delq(i,j)   = 0. _d 0            delq(i,j)   = 0. _d 0
277  #endif  #endif
278            if ( atemp(i,j,bi,bj) .ne. 0. _d 0 ) then  C--- Compute turbulent surface fluxes
279             tsf         = theta(i,j,k,bi,bj) + cen2kel  C-   Pot. Temp and saturated specific humidity
280             tmpbulk     = cvapor_fac * exp( - cvapor_exp/tsf )            IF ( atemp(i,j,bi,bj) .NE. 0. _d 0 ) THEN
281             ssq         = saltsat*tmpbulk/atmrho  C-   Surface Temp.
282             deltap(i,j) = atemp(i,j,bi,bj) + gamma_blk*ht - tsf             Tsf = theta(i,j,ksrf,bi,bj) + cen2kel
283               IF ( sstExtrapol.GT.0. _d 0 ) THEN
284                SSTtmp = sstExtrapol
285         &              *( theta(i,j,ksrf,bi,bj)-theta(i,j,ksrfp1,bi,bj) )
286         &              *  maskC(i,j,ksrfp1,bi,bj)
287                Tsf = Tsf + MAX( SSTtmp, 0. _d 0 )
288               ENDIF
289    c
290               tmpbulk = cvapor_fac*exp(-cvapor_exp/Tsf)
291               ssq = saltsat*tmpbulk/atmrho
292               deltap(i,j) = atemp(i,j,bi,bj) + gamma_blk*ht - Tsf
293             delq(i,j)   = aqh(i,j,bi,bj) - ssq             delq(i,j)   = aqh(i,j,bi,bj) - ssq
294             stable      = exf_half + sign(exf_half, deltap(i,j))  C--  no negative evaporation over ocean:
295               IF ( noNegativeEvap ) delq(i,j) = MIN( 0. _d 0, delq(i,j) )
296    
297    C--  initial guess for exchange coefficients:
298    C    take U_N = del.U ; stability from del.Theta ;
299               stable = exf_half + SIGN(exf_half, deltap(i,j))
300    
301    #ifndef ALLOW_ATM_WIND
302               IF ( solve4Stress ) THEN
303    #endif
304  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
305  CADJ STORE sh(i,j,bi,bj)     = comlev1_exf_1, key = ikey_1  CADJ STORE sh(i,j,bi,bj)     = comlev1_exf_1, key = ikey_1
306  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
307  #ifdef ALLOW_ATM_WIND  C--   Wind speed
308             wsm         = sh(i,j,bi,bj)              wsm        = sh(i,j,bi,bj)
309             tmpbulk     = exf_scal_BulkCdn *              tmpbulk    = exf_scal_BulkCdn *
310       &          ( cdrag_1/wsm + cdrag_2 + cdrag_3*wsm )       &           ( cdrag_1/wsm + cdrag_2 + cdrag_3*wsm )
311             rdn(i,j)    = sqrt(tmpbulk)              rdn(i,j)   = SQRT(tmpbulk)
312             ustar(i,j)  = rdn(i,j)*wsm              ustar(i,j) = rdn(i,j)*wsm
313  #else /* ifndef ALLOW_ATM_WIND */  #ifndef ALLOW_ATM_WIND
314  C     in this case ustress and vstress are defined a u and v points             ELSE
315  C     respectively, and we need to average to mass points to avoid              windStress = wStress(i,j,bi,bj)
316  C     tau = 0 near coasts or other boundaries              ustar(i,j) = sqrt(windStress/atmrho)
317  c            tau(i,j)   = sqrt(0.5*  c           tau(i,j)   = windStress/ustar(i,j)
318  c    &                   (ustress(i,  j,bi,bj)*ustress(i  ,j,bi,bj)              tau(i,j)   = sqrt(windStress*atmrho)
319  c    &                   +ustress(i+1,j,bi,bj)*ustress(i+1,j,bi,bj)             ENDIF
320  c    &                   +vstress(i,j,  bi,bj)*vstress(i,j  ,bi,bj)  #endif /* ALLOW_ATM_WIND */
 c    &                   +vstress(i,j+1,bi,bj)*vstress(i,j+1,bi,bj))  
 c    &                   )  
            windStress  = wStress(i,j,bi,bj)  
            ustar(i,j)  = sqrt(windStress/atmrho)  
 c          tau(i,j)    = wStress/ustar  
            tau(i,j)    = sqrt(windStress*atmrho)  
 #endif /* ifndef ALLOW_ATM_WIND */  
321    
322  C--  initial guess for exchange other coefficients:  C--  initial guess for exchange other coefficients:
323             rhn = (exf_one-stable)*cstanton_1 + stable*cstanton_2             rhn = (exf_one-stable)*cstanton_1 + stable*cstanton_2
324             ren = cdalton             ren = cDalton
325  C--  calculate turbulent scales  C--  calculate turbulent scales
326             tstar(i,j)  = rhn*deltap(i,j)             tstar(i,j)=rhn*deltap(i,j)
327             qstar(i,j)  = ren*delq(i,j)             qstar(i,j)=ren*delq(i,j)
328    
329            else            ELSE
330  C     atemp(i,j,bi,bj) .EQ. 0.  C     atemp(i,j,bi,bj) .EQ. 0.
331             tstar (i,j) = 0. _d 0             tstar (i,j) = 0. _d 0
332             qstar (i,j) = 0. _d 0             qstar (i,j) = 0. _d 0
333             ustar (i,j) = 0. _d 0             ustar (i,j) = 0. _d 0
334             tau   (i,j) = 0. _d 0             tau   (i,j) = 0. _d 0
335             rdn   (i,j) = 0. _d 0             rdn   (i,j) = 0. _d 0
336            end if            ENDIF
337           end do           ENDDO
338          end do          ENDDO
339            DO iter=1,niter_bulk
340          do iter = 1,niter_bulk           DO j = 1,sNy
341           do j = 1,sny            DO i = 1,sNx
342            do i = 1,snx             IF ( atemp(i,j,bi,bj) .NE. 0. _d 0 ) THEN
343             if ( atemp(i,j,bi,bj) .ne. 0. _d 0 ) then  C--- iterate with psi-functions to find transfer coefficients
344                
345  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
346              ikey_2 = i              ikey_2 = i
347       &           + sNx*(j-1)       &           + sNx*(j-1)
# Line 300  C     atemp(i,j,bi,bj) .EQ. 0. Line 350  C     atemp(i,j,bi,bj) .EQ. 0.
350       &           + sNx*sNy*niter_bulk*max1*act2       &           + sNx*sNy*niter_bulk*max1*act2
351       &           + sNx*sNy*niter_bulk*max1*max2*act3       &           + sNx*sNy*niter_bulk*max1*max2*act3
352       &           + sNx*sNy*niter_bulk*max1*max2*max3*act4       &           + sNx*sNy*niter_bulk*max1*max2*max3*act4
353  CADJ STORE rdn(i,j)              = comlev1_exf_2, key = ikey_2  CADJ STORE rdn   (i,j)       = comlev1_exf_2, key = ikey_2
354  CADJ STORE ustar(i,j)            = comlev1_exf_2, key = ikey_2  CADJ STORE ustar (i,j)       = comlev1_exf_2, key = ikey_2
355  CADJ STORE qstar(i,j)            = comlev1_exf_2, key = ikey_2  CADJ STORE qstar (i,j)       = comlev1_exf_2, key = ikey_2
356  CADJ STORE tstar(i,j)            = comlev1_exf_2, key = ikey_2  CADJ STORE tstar (i,j)       = comlev1_exf_2, key = ikey_2
357  CADJ 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
358  CADJ STORE wspeed(i,j,bi,bj)     = comlev1_exf_2, key = ikey_2  CADJ STORE wspeed(i,j,bi,bj) = comlev1_exf_2, key = ikey_2
359  #endif  #endif
360    
361              t0     = atemp(i,j,bi,bj)*              t0   = atemp(i,j,bi,bj)*
362       &           (exf_one + humid_fac*aqh(i,j,bi,bj))       &           (exf_one + humid_fac*aqh(i,j,bi,bj))
363              huol   = ( tstar(i,j)/t0 +              huol = ( tstar(i,j)/t0 +
364       &                 qstar(i,j)/(exf_one/humid_fac+aqh(i,j,bi,bj))       &               qstar(i,j)/(exf_one/humid_fac+aqh(i,j,bi,bj))
365       &                )*czol/(ustar(i,j)*ustar(i,j))       &              )*czol/(ustar(i,j)*ustar(i,j))
366    #ifdef ALLOW_BULK_LARGEYEAGER04
367                tmpbulk = MIN(abs(huol),10. _d 0)
368                huol   = SIGN(tmpbulk , huol)
369    #else
370    C--   Large&Pond1981:
371              huol   = max(huol,zolmin)              huol   = max(huol,zolmin)
372    #endif /* ALLOW_BULK_LARGEYEAGER04 */
373              htol   = huol*ht/hu              htol   = huol*ht/hu
374              hqol   = huol*hq/hu              hqol   = huol*hq/hu
375              stable = exf_half + sign(exf_half, huol)              stable = exf_half + sign(exf_half, huol)
376                
377  c     Evaluate all stability functions assuming hq = ht.  c                 Evaluate all stability functions assuming hq = ht.
378  #ifdef ALLOW_ATM_WIND  #ifndef ALLOW_ATM_WIND
379              xsq    = max(sqrt(abs(exf_one - 16.*huol)),exf_one)              IF ( solve4Stress ) THEN
380              x      = sqrt(xsq)  #endif
381              psimh  = -psim_fac*huol*stable  #ifdef ALLOW_BULK_LARGEYEAGER04
382       &              +(exf_one-stable)*  C--   Large&Yeager04:
383       &               ( LOG( (exf_one + exf_two*x + xsq)*               xsq    = SQRT( ABS(exf_one - huol*16. _d 0) )
384       &                      (exf_one+xsq)*.125 _d 0 )  #else
385       &                 -exf_two*ATAN(x) + exf_half*pi )  C--   Large&Pond1981:
386  #endif /* ALLOW_ATM_WIND */               xsq    = max(sqrt(abs(exf_one - 16.*huol)),exf_one)
387              xsq    = max(sqrt(abs(exf_one - 16.*htol)),exf_one)  #endif /* ALLOW_BULK_LARGEYEAGER04 */
388                 x      = SQRT(xsq)
389                 psimh  = -psim_fac*huol*stable
390         &               +(exf_one-stable)*
391         &                ( LOG( (exf_one + exf_two*x + xsq)*
392         &                       (exf_one+xsq)*.125 _d 0 )
393         &                  -exf_two*ATAN(x) + exf_half*pi )
394    #ifndef ALLOW_ATM_WIND
395                ENDIF
396    #endif
397    #ifdef ALLOW_BULK_LARGEYEAGER04
398    C--   Large&Yeager04:
399                xsq   = SQRT( ABS(exf_one - htol*16. _d 0) )
400    #else
401    C--   Large&Pond1981:
402                xsq   = max(sqrt(abs(exf_one - 16.*htol)),exf_one)
403    #endif /* ALLOW_BULK_LARGEYEAGER04 */
404              psixh = -psim_fac*htol*stable + (exf_one-stable)*              psixh = -psim_fac*htol*stable + (exf_one-stable)*
405       &               ( exf_two*LOG(exf_half*(exf_one+xsq)) )       &               ( exf_two*LOG(exf_half*(exf_one+xsq)) )
406                
407  #ifdef ALLOW_ATM_WIND  #ifndef ALLOW_ATM_WIND
408  c     Shift wind speed using old coefficient              IF ( solve4Stress ) THEN
409    #endif
410    C-   Shift wind speed using old coefficient
411    #ifdef ALLOW_BULK_LARGEYEAGER04
412    C--   Large&Yeager04:
413                 dzTmp = (zwln-psimh)/karman
414                 usn   = wspeed(i,j,bi,bj)/(exf_one + rdn(i,j)*dzTmp )
415    #else
416    C--   Large&Pond1981:
417  c           rd(i,j)= rdn(i,j)/(exf_one - rdn(i,j)/karman*psimh )  c           rd(i,j)= rdn(i,j)/(exf_one - rdn(i,j)/karman*psimh )
418  c           shn    = sh(i,j,bi,bj)*rd(i,j)/rdn(i,j)  c           usn    = sh(i,j,bi,bj)*rd(i,j)/rdn(i,j)
419  c     ML: the original formulation above is replaced below, but does  c     ML: the original formulation above is replaced below to be
420  c     not give the same results, strange  c     similar to largeyeager04, but does not give the same results, strange
421              shn    = sh(i,j,bi,bj)/(exf_one - rdn(i,j)/karman*psimh )               usn   = sh(i,j,bi,bj)/(exf_one - rdn(i,j)/karman*psimh)
422              uzn    = max(shn, umin)  #endif /* ALLOW_BULK_LARGEYEAGER04 */
423                 usm   = MAX(usn, umin)
424  c     Update the transfer coefficients at 10 meters  
425  c     and neutral stability.  C-   Update the 10m, neutral stability transfer coefficients (momentum)
426                 tmpbulk  = exf_scal_BulkCdn *
427              tmpbulk  = exf_scal_BulkCdn *       &            ( cdrag_1/usm + cdrag_2 + cdrag_3*usm )
428       &           ( cdrag_1/uzn + cdrag_2 + cdrag_3*uzn )               rdn(i,j)   = SQRT(tmpbulk)
429              rdn(i,j) = sqrt(tmpbulk)  #ifdef ALLOW_BULK_LARGEYEAGER04
430                 rd(i,j)    = rdn(i,j)/(exf_one + rdn(i,j)*dzTmp)
431  c     Shift all coefficients to the measurement height  #else
432  c     and stability.               rd(i,j)    = rdn(i,j)/(exf_one - rdn(i,j)/karman*psimh)
433              rd(i,j)  = rdn(i,j)/(exf_one - rdn(i,j)/karman*psimh)  #endif /* ALLOW_BULK_LARGEYEAGER04 */
434  #endif /* ALLOW_ATM_WIND */               ustar(i,j) = rd(i,j)*sh(i,j,bi,bj)
435    C-   Coeff:
436                 tau(i,j)   = atmrho*rd(i,j)*wspeed(i,j,bi,bj)
437    #ifndef ALLOW_ATM_WIND
438                ENDIF
439    #endif
440    
441    C-   Update the 10m, neutral stability transfer coefficients (sens&evap)
442              rhn = (exf_one-stable)*cstanton_1 + stable*cstanton_2              rhn = (exf_one-stable)*cstanton_1 + stable*cstanton_2
443              ren = cdalton              ren = cDalton
444              rh     = rhn/( exf_one + rhn/karman*(ztln - psixh) )  
445              re     = ren/( exf_one + ren/karman*(ztln - psixh) )  C-   Shift all coefficients to the measurement height and stability.
446                rh = rhn/(exf_one + rhn*(ztln-psixh)/karman)
447                re = ren/(exf_one + ren*(ztln-psixh)/karman)
448    
449  c                 Update ustar, tstar, qstar using updated, shifted  C--  Update ustar, tstar, qstar using updated, shifted coefficients.
 c                 coefficients.  
450              qstar(i,j) = re*delq(i,j)              qstar(i,j) = re*delq(i,j)
451              tstar(i,j) = rh*deltap(i,j)              tstar(i,j) = rh*deltap(i,j)
 #ifdef ALLOW_ATM_WIND  
             ustar(i,j) = rd(i,j)*sh(i,j,bi,bj)  
             tau(i,j)   = atmrho*rd(i,j)*wspeed(i,j,bi,bj)  
 #endif /* ALLOW_ATM_WIND */  
452    
453             endif             ENDIF
454            enddo            ENDDO
455           enddo           ENDDO
456  c        end of bulk_iter loop  c end of iteration loop
457          enddo          ENDDO
458            DO j = 1,sNy
459          do j = 1,sny           DO i = 1,sNx
460           do i = 1,snx            IF ( atemp(i,j,bi,bj) .NE. 0. _d 0 ) THEN
           if ( atemp(i,j,bi,bj) .ne. 0. _d 0 ) then  
461    
462  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
463            ikey_1 = i             ikey_1 = i
464       &         + sNx*(j-1)       &          + sNx*(j-1)
465       &         + sNx*sNy*act1       &          + sNx*sNy*act1
466       &         + sNx*sNy*max1*act2       &          + sNx*sNy*max1*act2
467       &         + sNx*sNy*max1*max2*act3       &          + sNx*sNy*max1*max2*act3
468       &         + sNx*sNy*max1*max2*max3*act4       &          + sNx*sNy*max1*max2*max3*act4
469  CADJ STORE qstar(i,j)    = comlev1_exf_1, key = ikey_1  CADJ STORE qstar(i,j)    = comlev1_exf_1, key = ikey_1
470  CADJ STORE tstar(i,j)    = comlev1_exf_1, key = ikey_1  CADJ STORE tstar(i,j)    = comlev1_exf_1, key = ikey_1
471  CADJ STORE tau(i,j)      = comlev1_exf_1, key = ikey_1  CADJ STORE tau  (i,j)    = comlev1_exf_1, key = ikey_1
472  CADJ STORE rd(i,j)       = comlev1_exf_1, key = ikey_1  CADJ STORE rd   (i,j)    = comlev1_exf_1, key = ikey_1
473  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
474    
475    C-   Turbulent Fluxes
476             hs(i,j,bi,bj) = atmcp*tau(i,j)*tstar(i,j)             hs(i,j,bi,bj) = atmcp*tau(i,j)*tstar(i,j)
477             hl(i,j,bi,bj) = flamb*tau(i,j)*qstar(i,j)             hl(i,j,bi,bj) = flamb*tau(i,j)*qstar(i,j)
478  #ifndef EXF_READ_EVAP  #ifndef EXF_READ_EVAP
479  cdm        evap(i,j,bi,bj) = tau(i,j)*qstar(i,j)  C   change sign and convert from kg/m^2/s to m/s via rhoConstFresh
480  cdm !!! need to change sign and to convert from kg/m^2/s to m/s !!!  c          evap(i,j,bi,bj) = -recip_rhonil*tau(i,j)*qstar(i,j)
481             evap(i,j,bi,bj) = -recip_rhoConstFresh*tau(i,j)*qstar(i,j)             evap(i,j,bi,bj) = -recip_rhoConstFresh*tau(i,j)*qstar(i,j)
482    C   but older version was using rhonil instead:
483    c           evap(i,j,bi,bj) = -recip_rhonil*tau(i,j)*qstar(i,j)
484  #endif  #endif
485  #ifdef ALLOW_ATM_WIND  #ifdef ALLOW_ATM_WIND
486    c       ustress(i,j,bi,bj) = tau(i,j)*rd(i,j)*ws*cw(i,j,bi,bj)
487    c       vstress(i,j,bi,bj) = tau(i,j)*rd(i,j)*ws*sw(i,j,bi,bj)
488    C- jmc: below is how it should be written ; different from above when
489    C       both wind-speed & 2 compon. of the wind are specified, and in
490    C-      this case, formula below is better.
491             tmpbulk =  tau(i,j)*rd(i,j)             tmpbulk =  tau(i,j)*rd(i,j)
492             ustress(i,j,bi,bj) = tmpbulk*uwind(i,j,bi,bj)             ustress(i,j,bi,bj) = tmpbulk*uwind(i,j,bi,bj)
493             vstress(i,j,bi,bj) = tmpbulk*vwind(i,j,bi,bj)             vstress(i,j,bi,bj) = tmpbulk*vwind(i,j,bi,bj)
494  #endif /* ALLOW_ATM_WIND */  #endif
495    
496  c IF ( ATEMP(i,j,bi,bj) .NE. 0. )  c IF ( ATEMP(i,j,bi,bj) .NE. 0. )
497            else            else
# Line 408  c IF ( ATEMP(i,j,bi,bj) .NE. 0. ) Line 500  c IF ( ATEMP(i,j,bi,bj) .NE. 0. )
500             vstress(i,j,bi,bj) = 0. _d 0             vstress(i,j,bi,bj) = 0. _d 0
501  #endif /* ALLOW_ATM_WIND */  #endif /* ALLOW_ATM_WIND */
502             hflux  (i,j,bi,bj) = 0. _d 0             hflux  (i,j,bi,bj) = 0. _d 0
503             hs(i,j,bi,bj)      = 0. _d 0             evap   (i,j,bi,bj) = 0. _d 0
504             hl(i,j,bi,bj)      = 0. _d 0             hs     (i,j,bi,bj) = 0. _d 0
505               hl     (i,j,bi,bj) = 0. _d 0
506    
507  c IF ( ATEMP(i,j,bi,bj) .NE. 0. )  c IF ( ATEMP(i,j,bi,bj) .NE. 0. )
508            endif            endif
509    
510  #else  /* ifndef ALLOW_ATM_TEMP */  #else  /* ifndef ALLOW_ATM_TEMP */
511  #ifdef ALLOW_ATM_WIND  #ifdef ALLOW_ATM_WIND
512            wsm                = sh(i,j,bi,bj)            wsm     = sh(i,j,bi,bj)
513            tmpbulk            = exf_scal_BulkCdn *            tmpbulk = exf_scal_BulkCdn *
514       &         ( cdrag_1/wsm + cdrag_2 + cdrag_3*wsm )       &         ( cdrag_1/wsm + cdrag_2 + cdrag_3*wsm )
515            ustress(i,j,bi,bj) = atmrho*tmpbulk*wspeed(i,j,bi,bj)*            ustress(i,j,bi,bj) = atmrho*tmpbulk*wspeed(i,j,bi,bj)*
516       &         uwind(i,j,bi,bj)       &                         uwind(i,j,bi,bj)
517            vstress(i,j,bi,bj) = atmrho*tmpbulk*wspeed(i,j,bi,bj)*            vstress(i,j,bi,bj) = atmrho*tmpbulk*wspeed(i,j,bi,bj)*
518       &         vwind(i,j,bi,bj)       &                         vwind(i,j,bi,bj)
519  #endif  #endif
520  #endif /* ifndef ALLOW_ATM_TEMP */  #endif /* ifndef ALLOW_ATM_TEMP */
521           enddo           ENDDO
522          enddo          ENDDO
523         enddo         ENDDO
524        enddo        ENDDO
525    
526  #endif /* ALLOW_BULKFORMULAE */  #endif /* ALLOW_BULKFORMULAE */
527    
528        return        RETURN
529        end        END
Legend:
Removed from v.1.22  
changed lines
  Added in v.1.23
  ViewVC Help
Powered by ViewVC 1.1.22