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

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

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

revision 1.2 by heimbach, Mon Jul 30 20:41:20 2001 UTC revision 1.2.6.4 by dimitri, Fri Feb 14 23:10:35 2003 UTC
# Line 2  c $Header$ Line 2  c $Header$
2    
3  #include "EXF_CPPOPTIONS.h"  #include "EXF_CPPOPTIONS.h"
4    
5        subroutine exf_GetFFields(        subroutine exf_GetFFields( mycurrenttime, mycurrentiter, mythid )
      I                           mycurrenttime,  
      I                           mycurrentiter,  
      I                           mythid  
      &                         )  
6    
7  c     ==================================================================  c     ==================================================================
8  c     SUBROUTINE exf_GetFFields  c     SUBROUTINE exf_GetFFields
9  c     ==================================================================  c     ==================================================================
10  c  c
11  c     o Get the surface fluxes either from file or as derived from bulk  c     o Read-in atmospheric state and/or surface fluxes from files.
12  c       formulae that use the atmospheric state.  c
13    c     o Use bulk formulae to estimate turbulent and/or radiative
14    c       fluxes at the surface.
15    c
16    c     NOTES:
17    c     ======
18    c
19    c     See EXF_CPPOPTIONS.h for a description of the various possible
20    c     ocean-model forcing configurations.
21    c
22    c     The bulk formulae of pkg/exf are not valid for sea-ice covered
23    c     oceans but they can be used in combination with a sea-ice model,
24    c     for example, pkg/seaice, to specify open water flux contributions.
25    c
26    c     ==================================================================
27  c  c
28  c       The calculation of the bulk surface fluxes has been adapted from  c       The calculation of the bulk surface fluxes has been adapted from
29  c       the NCOM model which uses the formulae given in Large and Pond  c       the NCOM model which uses the formulae given in Large and Pond
# Line 66  c              4) this version is for an Line 76  c              4) this version is for an
76  c                 ht = hq.  c                 ht = hq.
77  c              5) sst enters in Celsius.  c              5) sst enters in Celsius.
78  c  c
79  c       ------------------------------------  c     ==================================================================
 c  
 c       By setting CPP options in the header file *EXF_CPPOPTIONS.h* it  
 c       is possible to combine data sets in four different ways:  
 c  
 c       The following options are available:  
 c  
 c         ALLOW_ATM_TEMP (UAT)  
 c         ALLOW_ATM_WIND (UAW)  
 c  
 c  
 c             UAT    |    UAW      |         action  
 c         ----------------------------------------------------  
 c         undefined  |  undefined  |  Use surface fluxes.  
 c         undefined  |    defined  |  Assume cdn(u) given to  
 c                    |             |  infer the wind stress.  
 c           defined  |  undefined  |  Compute wind field from  
 c                    |             |  given stress assuming a  
 c                    |             |  linear relation.  
 c           defined  |    defined  |  Use the bulk formulae.  
 c         ----------------------------------------------------  
 c  
 c       Implementations of the bulk formulae exist for the follwing  
 c       versions of the MITgcm:  
 c  
 c       MITgcm  : Patrick Heimbach  
 c       MITgcmUV: Christian Eckert  
80  c  c
81  c       started: Christian Eckert eckert@mit.edu 27-Aug-1999  c       started: Christian Eckert eckert@mit.edu 27-Aug-1999
82  c  c
83  c       changed: Christian Eckert eckert@mit.edu 14-Jan-2000  c       changed: Christian Eckert eckert@mit.edu 14-Jan-2000
 c  
84  c              - restructured the original version in order to have a  c              - restructured the original version in order to have a
85  c                better interface to the MITgcmUV.  c                better interface to the MITgcmUV.
86  c  c
87  c              Christian Eckert eckert@mit.edu  12-Feb-2000  c              Christian Eckert eckert@mit.edu  12-Feb-2000
 c  
88  c              - Changed Routine names (package prefix: exf_)  c              - Changed Routine names (package prefix: exf_)
89  c  c
90  c              Patrick Heimbach, heimbach@mit.edu  04-May-2000  c              Patrick Heimbach, heimbach@mit.edu  04-May-2000
 c  
91  c              - changed the handling of precip and sflux with respect  c              - changed the handling of precip and sflux with respect
92  c                to CPP options ALLOW_BULKFORMULAE and ALLOW_ATM_TEMP  c                to CPP options ALLOW_BULKFORMULAE and ALLOW_ATM_TEMP
 c  
93  c              - included some CPP flags ALLOW_BULKFORMULAE to make  c              - included some CPP flags ALLOW_BULKFORMULAE to make
94  c                sure ALLOW_ATM_TEMP, ALLOW_ATM_WIND are used only in  c                sure ALLOW_ATM_TEMP, ALLOW_ATM_WIND are used only in
95  c                conjunction with defined ALLOW_BULKFORMULAE  c                conjunction with defined ALLOW_BULKFORMULAE
 c  
96  c              - statement functions discarded  c              - statement functions discarded
97  c  c
98  c              Ralf.Giering@FastOpt.de 25-Mai-2000  c              Ralf.Giering@FastOpt.de 25-Mai-2000
 c  
99  c              - total rewrite using new subroutines  c              - total rewrite using new subroutines
100  c  c
101    c              Detlef Stammer: include river run-off. Nov. 21, 2001
102    c
103    c              heimbach@mit.edu, 10-Jan-2002
104    c              - changes to enable field swapping
105    c
106    c       mods for pkg/seaice: menemenlis@jpl.nasa.gov 20-Dec-2002
107    c
108  c     ==================================================================  c     ==================================================================
109  c     SUBROUTINE exf_GetFFields  c     SUBROUTINE exf_GetFFields
110  c     ==================================================================  c     ==================================================================
# Line 137  c     == global variables == Line 122  c     == global variables ==
122  #include "exf_fields.h"  #include "exf_fields.h"
123  #include "exf_constants.h"  #include "exf_constants.h"
124    
125    #ifdef ALLOW_AUTODIFF_TAMC
126    #include "tamc.h"
127    #endif
128    
129  c     == routine arguments ==  c     == routine arguments ==
130    
131        integer mythid        integer mythid
# Line 165  c     == local variables == Line 154  c     == local variables ==
154        _RL     re        _RL     re
155        _RL     rdn        _RL     rdn
156        _RL     rh        _RL     rh
157          _RL     ssttmp
158        _RL     ssq        _RL     ssq
159        _RL     stable        _RL     stable
       _RL     tau  
160        _RL     tstar        _RL     tstar
161        _RL     t0        _RL     t0
162        _RL     ustar        _RL     ustar
163        _RL     uzn        _RL     uzn
164          _RL     shn
165        _RL     xsq        _RL     xsq
166        _RL     x        _RL     x
167        _RL     evap(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)        _RL     tau
168    #ifdef ALLOW_AUTODIFF_TAMC
169          integer ikey_1
170          integer ikey_2
171    #endif
172  #endif /* ALLOW_ATM_TEMP */  #endif /* ALLOW_ATM_TEMP */
173    
174          _RL     ustmp
175        _RL     us        _RL     us
176        _RL     cw        _RL     cw
177        _RL     sw        _RL     sw
# Line 201  c     == external functions == Line 196  c     == external functions ==
196        external  exf_BulkRhn        external  exf_BulkRhn
197  #endif /* ALLOW_BULKFORMULAE */  #endif /* ALLOW_BULKFORMULAE */
198    
199  c     == end of interface ==  #ifndef ALLOW_ATM_WIND
200          _RL       TMP1
201          _RL       TMP2
202          _RL       TMP3
203          _RL       TMP4
204          _RL       TMP5
205    #endif
206    
207  c     determine forcing field records  c     == end of interface ==
208    
209  #ifdef ALLOW_BULKFORMULAE  #ifdef ALLOW_BULKFORMULAE
210    cph   This statement cannot be a PARAMETER statement in the header,
211    cph   but must come here; it's not fortran77 standard
212          aln = log(ht/zref)
213    #endif
214    
215  c     Determine where we are in time and set counters, flags and  c--   read forcing fields from files and temporal interpolation
 c     the linear interpolation factors accordingly.  
 #ifdef ALLOW_ATM_TEMP  
 c     Atmospheric temperature.  
       call exf_set_atemp( atemp  
      &                  , mycurrenttime, mycurrentiter, mythid )  
216    
217  c     Atmospheric humidity.  #ifdef ALLOW_ATM_WIND
       call exf_set_aqh( aqh  
      &                , mycurrenttime, mycurrentiter, mythid )  
218    
219  c     Net long wave radiative flux.  c     Zonal wind.
220        call exf_set_lwflux( lwflux        call exf_set_uwind  ( mycurrenttime, mycurrentiter, mythid )
      &                   , mycurrenttime, mycurrentiter, mythid )  
221    
222  c     Net short wave radiative flux.  c     Meridional wind.
223        call exf_set_swflux( swflux        call exf_set_vwind  ( mycurrenttime, mycurrentiter, mythid )
      &                   , mycurrenttime, mycurrentiter, mythid )  
224    
225  c     Precipitation.  #ifdef ALLOW_UWIND_CONTROL
226        call exf_set_precip( precip        call ctrl_getuwind  ( mycurrenttime, mycurrentiter, mythid )
227       &                   , mycurrenttime, mycurrentiter, mythid )  #endif
228    
229        aln = log(ht/zref)  #ifdef ALLOW_VWIND_CONTROL
230          call ctrl_getvwind  ( mycurrenttime, mycurrentiter, mythid )
231    #endif
232    
233  #else  #else  /* ifndef ALLOW_ATM_WIND */
234    
235  c     Atmospheric heat flux.  c     Zonal wind stress.
236        call exf_set_hflux( hflux        call exf_set_ustress( mycurrenttime, mycurrentiter, mythid )
      &                  , mycurrenttime, mycurrentiter, mythid )  
237    
238  c     Salt flux.  c     Meridional wind stress.
239        call exf_set_sflux( sflux        call exf_set_vstress( mycurrenttime, mycurrentiter, mythid )
      &                  , mycurrenttime, mycurrentiter, mythid )  
240    
241  #ifdef ALLOW_KPP  #endif /* ifndef ALLOW_ATM_WIND */
 c     Net short wave radiative flux.  
       call exf_set_swflux( swflux  
      &                   , mycurrenttime, mycurrentiter, mythid )  
242    
243  #endif /* ALLOW_KPP */  #ifdef ALLOW_ATM_TEMP
244    
245  #endif /* ALLOW_ATM_TEMP */  c     Atmospheric temperature.
246          call exf_set_atemp  ( mycurrenttime, mycurrentiter, mythid )
247    
248  #ifdef ALLOW_ATM_WIND  c     Atmospheric humidity.
249  c     Zonal wind.        call exf_set_aqh    ( mycurrenttime, mycurrentiter, mythid )
       call exf_set_uwind( uwind  
      &                  , mycurrenttime, mycurrentiter, mythid )  
250    
251  c     Meridional wind.  c     Net long wave radiative flux.
252        call exf_set_vwind( vwind        call exf_set_lwflux ( mycurrenttime, mycurrentiter, mythid )
      &                  , mycurrenttime, mycurrentiter, mythid )  
253    
254  #else  c     Precipitation.
255  c     Zonal wind stress.        call exf_set_precip ( mycurrenttime, mycurrentiter, mythid )
       call exf_set_ustress( ustress  
      &                    , mycurrenttime, mycurrentiter, mythid )  
256    
257  c     Meridional wind stress.  #ifdef ALLOW_ATEMP_CONTROL
258        call exf_set_vstress( vstress        call ctrl_getatemp  ( mycurrenttime, mycurrentiter, mythid )
259       &                    , mycurrenttime, mycurrentiter, mythid )  #endif
260    
261  #endif /* ALLOW_ATM_WIND */  #ifdef ALLOW_AQH_CONTROL
262          call ctrl_getaqh    ( mycurrenttime, mycurrentiter, mythid )
263    #endif
264    
265    #else  /* ifndef ALLOW_ATM_TEMP */
266    
 #else /* ALLOW_BULKFORMULAE undefined */  
267  c     Atmospheric heat flux.  c     Atmospheric heat flux.
268        call exf_set_hflux( hflux        call exf_set_hflux  ( mycurrenttime, mycurrentiter, mythid )
      &                  , mycurrenttime, mycurrentiter, mythid )  
269    
270  c     Salt flux.  c     Salt flux.
271        call exf_set_sflux( sflux        call exf_set_sflux  ( mycurrenttime, mycurrentiter, mythid )
      &                  , mycurrenttime, mycurrentiter, mythid )  
272    
273  c     Zonal wind stress.  #endif /* ifndef ALLOW_ATM_TEMP */
       call exf_set_ustress( ustress  
      &                    , mycurrenttime, mycurrentiter, mythid )  
274    
275  c     Meridional wind stress.  #if defined(ALLOW_ATM_TEMP) || defined(SHORTWAVE_HEATING)
       call exf_set_vstress( vstress  
      &                    , mycurrenttime, mycurrentiter, mythid )  
   
 #ifdef ALLOW_KPP  
276  c     Net short wave radiative flux.  c     Net short wave radiative flux.
277        call exf_set_swflux( swflux        call exf_set_swflux ( mycurrenttime, mycurrentiter, mythid )
278       &                   , mycurrenttime, mycurrentiter, mythid )  #endif
 #endif /* ALLOW_KPP */  
279    
280  #endif /* ALLOW_BULKFORMULAE */  #ifdef EXF_READ_EVAP
281    c     Evaporation
282          call exf_set_evap   ( mycurrenttime, mycurrentiter, mythid )
283    #endif
284    
285    #ifdef ALLOW_DOWNWARD_RADIATION
286    
287    c     Downward shortwave radiation.
288          call exf_set_swdown ( mycurrenttime, mycurrentiter, mythid )
289    
290    c     Downward longwave radiation.
291          call exf_set_lwdown ( mycurrenttime, mycurrentiter, mythid )
292    
293    #endif
294    
295    c--   Use atmospheric state to compute surface fluxes.
296    
297  c     Loop over tiles.  c     Loop over tiles.
298    #ifdef ALLOW_AUTODIFF_TAMC
299    C--   HPF directive to help TAMC
300    CHPF$ INDEPENDENT
301    #endif
302        do bj = mybylo(mythid),mybyhi(mythid)        do bj = mybylo(mythid),mybyhi(mythid)
303    #ifdef ALLOW_AUTODIFF_TAMC
304    C--    HPF directive to help TAMC
305    CHPF$  INDEPENDENT
306    #endif
307          do bi = mybxlo(mythid),mybxhi(mythid)          do bi = mybxlo(mythid),mybxhi(mythid)
308    
309            k = 1            k = 1
310    
311    cdm? can olx, oly be eliminated?
312            do j = 1-oly,sny+oly            do j = 1-oly,sny+oly
313              do i = 1-olx,snx+olx              do i = 1-olx,snx+olx
314    
315  #ifdef ALLOW_BULKFORMULAE  #ifdef ALLOW_BULKFORMULAE
316    
317  c             Compute the turbulent surface fluxes.  #ifdef ALLOW_AUTODIFF_TAMC
318  c                  (Bulk formulae estimates)                 act1 = bi - myBxLo(myThid)
319                   max1 = myBxHi(myThid) - myBxLo(myThid) + 1
320                   act2 = bj - myByLo(myThid)
321                   max2 = myByHi(myThid) - myByLo(myThid) + 1
322                   act3 = myThid - 1
323                   max3 = nTx*nTy
324                   act4 = ikey_dynamics - 1
325    
326                   ikey_1 = i
327         &              + sNx*(j-1)
328         &              + sNx*sNy*act1
329         &              + sNx*sNy*max1*act2
330         &              + sNx*sNy*max1*max2*act3
331         &              + sNx*sNy*max1*max2*max3*act4
332    #endif
333    
334    #ifdef ALLOW_DOWNWARD_RADIATION
335    c--   Compute net longwave and shortwave radiation:
336    c     lwflux = Stefan-Boltzman constant * emissivity * SST - lwdown
337    c     swflux = - ( 1 - albedo ) * swdown
338                   lwflux(i,j,bi,bj) = 5.5 _d -08 *
339         &              ((theta(i,j,k,bi,bj)+cen2kel)**4)
340         &              - lwdown(i,j,bi,bj)
341                   swflux(i,j,bi,bj) = -0.9 _d 0 * swdown(i,j,bi,bj)
342    #endif
343    
344    c--   Compute the turbulent surface fluxes.
345    
346  #ifdef ALLOW_ATM_WIND  #ifdef ALLOW_ATM_WIND
347  c             Wind speed and direction.  c             Wind speed and direction.
348                us = sqrt(uwind(i,j,bi,bj)*uwind(i,j,bi,bj) +                ustmp = uwind(i,j,bi,bj)*uwind(i,j,bi,bj) +
349       &                  vwind(i,j,bi,bj)*vwind(i,j,bi,bj))       &                vwind(i,j,bi,bj)*vwind(i,j,bi,bj)
350                if ( us .ne. 0. _d 0 ) then                if ( ustmp .ne. 0. _d 0 ) then
351                    us = sqrt(ustmp)
352                  cw = uwind(i,j,bi,bj)/us                  cw = uwind(i,j,bi,bj)/us
353                  sw = vwind(i,j,bi,bj)/us                  sw = vwind(i,j,bi,bj)/us
354                else                else
355                    us = 0. _d 0
356                  cw = 0. _d 0                  cw = 0. _d 0
357                  sw = 0. _d 0                  sw = 0. _d 0
358                endif                endif
359                sh = max(us,umin)                sh = max(us,umin)
360  #else  #else  /* ifndef ALLOW_ATM_WIND */
361  #ifdef ALLOW_ATM_TEMP  #ifdef ALLOW_ATM_TEMP
362    
363  c             The variables us, sh and rdn have to be computed from given  c             The variables us, sh and rdn have to be computed from
364  c             wind stresses inverting relationship for neutral drag coeff.  c             given wind stresses inverting relationship for neutral
365  c             cdn.  c             drag coeff. cdn.
366  c             The inversion is based on linear and quadratic form of  c             The inversion is based on linear and quadratic form of
367  c             cdn(umps); ustar can be directly computed from stress;  c             cdn(umps); ustar can be directly computed from stress;
368    
369                ustar = sqrt(ustress(i,j,bi,bj)*ustress(i,j,bi,bj) +                ustmp = ustress(i,j,bi,bj)*ustress(i,j,bi,bj) +
370       &                     vstress(i,j,bi,bj)*vstress(i,j,bi,bj))/       &                vstress(i,j,bi,bj)*vstress(i,j,bi,bj)
371       &                    atmrho                if ( ustmp .ne. 0. _d 0 ) then
372                cw = ustress(i,j,bi,bj)/ustar                  ustar = sqrt(ustmp/atmrho)
373                sw = ustress(i,j,bi,bj)/ustar                  cw = ustress(i,j,bi,bj)/sqrt(ustmp)
374                    sw = vstress(i,j,bi,bj)/sqrt(ustmp)
375                  else
376                     ustar = 0. _d 0
377                     cw    = 0. _d 0
378                     sw    = 0. _d 0
379                  endif
380    
381                if ( ustar .eq. 0. _d 0 ) then                if ( ustar .eq. 0. _d 0 ) then
382                  us = 0. _d 0                  us = 0. _d 0
# Line 359  c             cdn(umps); ustar can be di Line 401  c             cdn(umps); ustar can be di
401    
402                sh    = max(us,umin)                sh    = max(us,umin)
403  #endif /* ALLOW_ATM_TEMP */  #endif /* ALLOW_ATM_TEMP */
404  #endif /* ALLOW_ATM_WIND */  #endif /* ifndef ALLOW_ATM_WIND */
405    
406  #ifdef ALLOW_ATM_TEMP  #ifdef ALLOW_ATM_TEMP
407    
408  c             Initial guess: z/l=0.0; hu=ht=hq=z  c             Initial guess: z/l=0.0; hu=ht=hq=z
409  c             Iterations:    converge on z/l and hence the fluxes.  c             Iterations:    converge on z/l and hence the fluxes.
410  c             t0     : virtual temperature (K)  c             t0     : virtual temperature (K)
# Line 371  c             deltap : potential tempera Line 414  c             deltap : potential tempera
414                if ( atemp(i,j,bi,bj) .ne. 0. _d 0 ) then                if ( atemp(i,j,bi,bj) .ne. 0. _d 0 ) then
415                  t0     = atemp(i,j,bi,bj)*                  t0     = atemp(i,j,bi,bj)*
416       &                   (exf_one + humid_fac*aqh(i,j,bi,bj))       &                   (exf_one + humid_fac*aqh(i,j,bi,bj))
417                    ssttmp = theta(i,j,k,bi,bj)
418                  ssq    = saltsat*                  ssq    = saltsat*
419       &                   exf_BulkqSat(theta(i,j,1,bi,bj) + cen2kel)/       &                   exf_BulkqSat(ssttmp + cen2kel)/
420       &                   atmrho       &                   atmrho
421                  deltap = atemp(i,j,bi,bj)   + gamma_blk*ht -                  deltap = atemp(i,j,bi,bj)   + gamma_blk*ht -
422       &                   theta(i,j,1,bi,bj) - cen2kel       &                   ssttmp - cen2kel
423                  delq   = aqh(i,j,bi,bj) - ssq                  delq   = aqh(i,j,bi,bj) - ssq
424                  stable = exf_half + sign(exf_half, deltap)                  stable = exf_half + sign(exf_half, deltap)
425    #ifdef ALLOW_AUTODIFF_TAMC
426    CADJ STORE sh     = comlev1_exf_1, key = ikey_1
427    #endif
428                  rdn    = sqrt(exf_BulkCdn(sh))                  rdn    = sqrt(exf_BulkCdn(sh))
429                  ustar  = rdn*sh                  ustar  = rdn*sh
430                  tstar  = exf_BulkRhn(stable)*deltap                  tstar  = exf_BulkRhn(stable)*deltap
# Line 385  c             deltap : potential tempera Line 432  c             deltap : potential tempera
432    
433                  do iter = 1,niter_bulk                  do iter = 1,niter_bulk
434    
435    #ifdef ALLOW_AUTODIFF_TAMC
436                       ikey_2 = iter
437         &                  + niter_bulk*(i-1)
438         &                  + sNx*niter_bulk*(j-1)
439         &                  + sNx*niter_bulk*sNy*act1
440         &                  + sNx*niter_bulk*sNy*max1*act2
441         &                  + sNx*niter_bulk*sNy*max1*max2*act3
442         &                  + sNx*niter_bulk*sNy*max1*max2*max3*act4
443    
444    CADJ STORE rdn    = comlev1_exf_2, key = ikey_2
445    CADJ STORE ustar  = comlev1_exf_2, key = ikey_2
446    CADJ STORE qstar  = comlev1_exf_2, key = ikey_2
447    CADJ STORE tstar  = comlev1_exf_2, key = ikey_2
448    CADJ STORE sh     = comlev1_exf_2, key = ikey_2
449    CADJ STORE us     = comlev1_exf_2, key = ikey_2
450    #endif
451    
452                    huol   = czol*(tstar/t0 +                    huol   = czol*(tstar/t0 +
453       &                     qstar/(exf_one/humid_fac+aqh(i,j,bi,bj)))/       &                     qstar/(exf_one/humid_fac+aqh(i,j,bi,bj)))/
454       &                           ustar**2       &                           ustar**2
# Line 406  c                 Evaluate all stability Line 470  c                 Evaluate all stability
470       &                     exf_two*log((exf_one + xsq)/exf_two)       &                     exf_two*log((exf_one + xsq)/exf_two)
471    
472  c                 Shift wind speed using old coefficient  c                 Shift wind speed using old coefficient
473  c                 rd   = rdn/(exf_one + rdn/karman*(log(hu/zref) - psimh) )  ccc                  rd   = rdn/(exf_one + rdn/karman*
474    ccc     &                 (log(hu/zref) - psimh) )
475                    rd   = rdn/(exf_one - rdn/karman*psimh )                    rd   = rdn/(exf_one - rdn/karman*psimh )
476                    uzn  = max(sh*rd/rdn, umin)                    shn  = sh*rd/rdn
477                      uzn  = max(shn, umin)
478    
479  c                 Update the transfer coefficients at 10 meters  c                 Update the transfer coefficients at 10 meters
480  c                 and neutral stability.  c                 and neutral stability.
481    
482                    rdn = sqrt(exf_BulkCdn(uzn))                    rdn = sqrt(exf_BulkCdn(uzn))
483    
484  c                 Shift all coefficients to the measurement height  c                 Shift all coefficients to the measurement height
# Line 428  c                 coefficients. Line 495  c                 coefficients.
495                    ustar = rd*sh                      ustar = rd*sh  
496                    qstar = re*delq                    qstar = re*delq
497                    tstar = rh*deltap                    tstar = rh*deltap
498                      tau   = atmrho*ustar**2
499                      tau   = tau*us/sh
500    
501                  enddo                  enddo
502    
503                  tau                = atmrho*ustar**2    #ifdef ALLOW_AUTODIFF_TAMC
504                  tau                = tau*us/sh  CADJ STORE ustar  = comlev1_exf_1, key = ikey_1
505    CADJ STORE qstar  = comlev1_exf_1, key = ikey_1
506    CADJ STORE tstar  = comlev1_exf_1, key = ikey_1
507    CADJ STORE tau    = comlev1_exf_1, key = ikey_1
508    CADJ STORE cw     = comlev1_exf_1, key = ikey_1
509    CADJ STORE sw     = comlev1_exf_1, key = ikey_1
510    #endif
511    
512                  hs(i,j,bi,bj)      = atmcp*tau*tstar/ustar                  hs(i,j,bi,bj)      = atmcp*tau*tstar/ustar
513                  hl(i,j,bi,bj)      = flamb*tau*qstar/ustar                  hl(i,j,bi,bj)      = flamb*tau*qstar/ustar
514    #ifndef EXF_READ_EVAP
515                  evap(i,j,bi,bj)    = tau*qstar/ustar  cdm             evap(i,j,bi,bj)    = tau*qstar/ustar
516    cdm !!! need to change sign and to convert from kg/m^2/s to m/s !!!
517                    evap(i,j,bi,bj)    = -recip_rhonil*tau*qstar/ustar
518    #endif
519                  ustress(i,j,bi,bj) = tau*cw                  ustress(i,j,bi,bj) = tau*cw
520                  vstress(i,j,bi,bj) = tau*sw                  vstress(i,j,bi,bj) = tau*sw
 ce              hflux(i,j,bi,bj)   = atmcp*tau*tstar/ustar +  
 ce     &                             flamb*tau*qstar/ustar  
521                else                else
522                  ustress(i,j,bi,bj) = 0. _d 0                  ustress(i,j,bi,bj) = 0. _d 0
523                  vstress(i,j,bi,bj) = 0. _d 0                  vstress(i,j,bi,bj) = 0. _d 0
# Line 450  ce     &                             fla Line 526  ce     &                             fla
526                  hl(i,j,bi,bj)      = 0. _d 0                  hl(i,j,bi,bj)      = 0. _d 0
527                endif                endif
528    
529  #else  #else  /* ifndef ALLOW_ATM_TEMP */
530  #ifdef ALLOW_ATM_WIND  #ifdef ALLOW_ATM_WIND
531                ustress(i,j,bi,bj) = atmrho*exf_BulkCdn(sh)*us*                ustress(i,j,bi,bj) = atmrho*exf_BulkCdn(sh)*us*
532       &                             uwind(i,j,bi,bj)       &                             uwind(i,j,bi,bj)
533                vstress(i,j,bi,bj) = atmrho*exf_BulkCdn(sh)*us*                vstress(i,j,bi,bj) = atmrho*exf_BulkCdn(sh)*us*
534       &                             vwind(i,j,bi,bj)       &                             vwind(i,j,bi,bj)
535  #endif /* ALLOW_ATM_WIND */  #endif
536  #endif /* ALLOW_ATM_TEMP */  #endif /* ifndef ALLOW_ATM_TEMP */
537              enddo              enddo
538            enddo            enddo
539          enddo          enddo
540        enddo        enddo
541    
542  c     Add all contributions.  c     Add all contributions.
       k = 1  
543        do bj = mybylo(mythid),mybyhi(mythid)        do bj = mybylo(mythid),mybyhi(mythid)
544          do bi = mybxlo(mythid),mybxhi(mythid)          do bi = mybxlo(mythid),mybxhi(mythid)
545            do j = 1,sny            do j = 1,sny
# Line 475  c             Net surface heat flux. Line 550  c             Net surface heat flux.
550                hfl = hfl - hs(i,j,bi,bj)                hfl = hfl - hs(i,j,bi,bj)
551                hfl = hfl - hl(i,j,bi,bj)                hfl = hfl - hl(i,j,bi,bj)
552                hfl = hfl + lwflux(i,j,bi,bj)                hfl = hfl + lwflux(i,j,bi,bj)
553  #ifndef ALLOW_KPP  #ifndef SHORTWAVE_HEATING
554                hfl = hfl + swflux(i,j,bi,bj)                hfl = hfl + swflux(i,j,bi,bj)
555  #endif /* ALLOW_KPP undef */  #endif
556  c             Heat flux:  c             Heat flux:
557                hflux(i,j,bi,bj) = hfl*maskc(i,j,k,bi,bj)                hflux(i,j,bi,bj) = hfl
558  c             Salt flux from Precipitation and Evaporation.  c             Salt flux from Precipitation and Evaporation.
559                sflux(i,j,bi,bj) = precip(i,j,bi,bj) - evap(i,j,bi,bj)                sflux(i,j,bi,bj) = evap(i,j,bi,bj) - precip(i,j,bi,bj)
560  #endif /* ALLOW_ATM_TEMP */  #endif /* ALLOW_ATM_TEMP */
561    
 #else  
               hflux(i,j,bi,bj) = hflux(i,j,bi,bj)*maskc(i,j,k,bi,bj)  
               sflux(i,j,bi,bj) = sflux(i,j,bi,bj)*maskc(i,j,k,bi,bj)  
562  #endif /* ALLOW_BULKFORMULAE */  #endif /* ALLOW_BULKFORMULAE */
563    
564    #ifdef ALLOW_RUNOFF
565                  sflux(i,j,bi,bj) = sflux(i,j,bi,bj) - runoff(i,j,bi,bj)
566    #endif
567    
568                  hflux(i,j,bi,bj) = hflux(i,j,bi,bj)*maskc(i,j,1,bi,bj)
569                  sflux(i,j,bi,bj) = sflux(i,j,bi,bj)*maskc(i,j,1,bi,bj)
570    
571              enddo              enddo
572            enddo            enddo
573          enddo          enddo
# Line 499  c     Update the tile edges. Line 579  c     Update the tile edges.
579        _EXCH_XY_R8(ustress, mythid)        _EXCH_XY_R8(ustress, mythid)
580        _EXCH_XY_R8(vstress, mythid)        _EXCH_XY_R8(vstress, mythid)
581    
582  #ifdef ALLOW_KPP  #ifdef SHORTWAVE_HEATING
583        _EXCH_XY_R8(swflux, mythid)        _EXCH_XY_R8(swflux, mythid)
584  #endif /* ALLOW_KPP */  #endif
585    
586        end        end
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.2.6.4
  ViewVC Help
Powered by ViewVC 1.1.22