/[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.4 by heimbach, Tue Nov 12 20:34:41 2002 UTC revision 1.10 by heimbach, Thu Mar 6 00:47:33 2003 UTC
# Line 8  c     ================================== Line 8  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 62  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
# Line 115  c Line 103  c
103  c              heimbach@mit.edu, 10-Jan-2002  c              heimbach@mit.edu, 10-Jan-2002
104  c              - changes to enable field swapping  c              - changes to enable field swapping
105  c  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 129  c     == global variables == Line 119  c     == global variables ==
119  #include "DYNVARS.h"  #include "DYNVARS.h"
120  #include "GRID.h"  #include "GRID.h"
121    
122    #include "exf_param.h"
123  #include "exf_fields.h"  #include "exf_fields.h"
124  #include "exf_constants.h"  #include "exf_constants.h"
125    
# Line 175  c     == local variables == Line 166  c     == local variables ==
166        _RL     xsq        _RL     xsq
167        _RL     x        _RL     x
168        _RL     tau        _RL     tau
       _RL     evap(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)  
169  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
170        integer ikey_1        integer ikey_1
171        integer ikey_2        integer ikey_2
# Line 207  c     == external functions == Line 197  c     == external functions ==
197        external  exf_BulkRhn        external  exf_BulkRhn
198  #endif /* ALLOW_BULKFORMULAE */  #endif /* ALLOW_BULKFORMULAE */
199    
200  c     == end of interface ==  #ifndef ALLOW_ATM_WIND
201          _RL       TMP1
202          _RL       TMP2
203          _RL       TMP3
204          _RL       TMP4
205          _RL       TMP5
206    #endif
207    
208  c     determine forcing field records  c     == end of interface ==
209    
210  #ifdef ALLOW_BULKFORMULAE  #ifdef ALLOW_BULKFORMULAE
   
 #if (defined (ALLOW_ATM_TEMP) || defined (ALLOW_ATM_WIND))  
211  cph   This statement cannot be a PARAMETER statement in the header,  cph   This statement cannot be a PARAMETER statement in the header,
212  cph   but must come here; it's not fortran77 standard  cph   but must come here; it's not fortran77 standard
213        aln = log(ht/zref)        aln = log(ht/zref)
214  #endif  #endif
215    
216  c     Determine where we are in time and set counters, flags and  c--   read forcing fields from files and temporal interpolation
217  c     the linear interpolation factors accordingly.  
218    #ifdef ALLOW_ATM_WIND
219    
220    c     Zonal wind.
221          call exf_set_uwind  ( mycurrenttime, mycurrentiter, mythid )
222    
223    c     Meridional wind.
224          call exf_set_vwind  ( mycurrenttime, mycurrentiter, mythid )
225    
226    #ifdef ALLOW_UWIND_CONTROL
227          call ctrl_getuwind  ( mycurrenttime, mycurrentiter, mythid )
228    #endif
229    
230    #ifdef ALLOW_VWIND_CONTROL
231          call ctrl_getvwind  ( mycurrenttime, mycurrentiter, mythid )
232    #endif
233    
234    #else  /* ifndef ALLOW_ATM_WIND */
235    
236    c     Zonal wind stress.
237          call exf_set_ustress( mycurrenttime, mycurrentiter, mythid )
238    
239    c     Meridional wind stress.
240          call exf_set_vstress( mycurrenttime, mycurrentiter, mythid )
241    
242    #endif /* ifndef ALLOW_ATM_WIND */
243    
244  #ifdef ALLOW_ATM_TEMP  #ifdef ALLOW_ATM_TEMP
245    
246  c     Atmospheric temperature.  c     Atmospheric temperature.
247        call exf_set_atemp( mycurrenttime, mycurrentiter, mythid )        call exf_set_atemp  ( mycurrenttime, mycurrentiter, mythid )
248    
249  c     Atmospheric humidity.  c     Atmospheric humidity.
250        call exf_set_aqh( mycurrenttime, mycurrentiter, mythid )        call exf_set_aqh    ( mycurrenttime, mycurrentiter, mythid )
251    
252  c     Net long wave radiative flux.  c     Net long wave radiative flux.
253        call exf_set_lwflux( mycurrenttime, mycurrentiter, mythid )        call exf_set_lwflux ( mycurrenttime, mycurrentiter, mythid )
   
 c     Net short wave radiative flux.  
       call exf_set_swflux( mycurrenttime, mycurrentiter, mythid )  
254    
255  c     Precipitation.  c     Precipitation.
256        call exf_set_precip( mycurrenttime, mycurrentiter, mythid )        call exf_set_precip ( mycurrenttime, mycurrentiter, mythid )
257    
258  #ifdef ALLOW_ATEMP_CONTROL  #ifdef ALLOW_ATEMP_CONTROL
259        call ctrl_getatemp ( mycurrenttime, mycurrentiter, mythid )        call ctrl_getatemp  ( mycurrenttime, mycurrentiter, mythid )
260  #endif  #endif
261    
262  #ifdef ALLOW_AQH_CONTROL  #ifdef ALLOW_AQH_CONTROL
263        call ctrl_getaqh ( mycurrenttime, mycurrentiter, mythid )        call ctrl_getaqh    ( mycurrenttime, mycurrentiter, mythid )
264  #endif  #endif
265    
266  #else  #else  /* ifndef ALLOW_ATM_TEMP */
267    
268  c     Atmospheric heat flux.  c     Atmospheric heat flux.
269        call exf_set_hflux( mycurrenttime, mycurrentiter, mythid )        call exf_set_hflux  ( mycurrenttime, mycurrentiter, mythid )
270    
271  c     Salt flux.  c     Salt flux.
272        call exf_set_sflux( mycurrenttime, mycurrentiter, mythid )        call exf_set_sflux  ( mycurrenttime, mycurrentiter, mythid )
   
 #ifdef ALLOW_KPP  
 c     Net short wave radiative flux.  
       call exf_set_swflux( mycurrenttime, mycurrentiter, mythid )  
   
 #endif /* ALLOW_KPP */  
273    
274  #endif /* ALLOW_ATM_TEMP */  #endif /* ifndef ALLOW_ATM_TEMP */
   
 #ifdef ALLOW_ATM_WIND  
 c     Zonal wind.  
       call exf_set_uwind( mycurrenttime, mycurrentiter, mythid )  
   
 c     Meridional wind.  
       call exf_set_vwind( mycurrenttime, mycurrentiter, mythid )  
275    
276  #ifdef ALLOW_UWIND_CONTROL  #if defined(ALLOW_ATM_TEMP) || defined(SHORTWAVE_HEATING)
277        call ctrl_getuwind ( mycurrenttime, mycurrentiter, mythid )  c     Net short wave radiative flux.
278          call exf_set_swflux ( mycurrenttime, mycurrentiter, mythid )
279  #endif  #endif
280    
281  #ifdef ALLOW_VWIND_CONTROL  #ifdef EXF_READ_EVAP
282        call ctrl_getvwind ( mycurrenttime, mycurrentiter, mythid )  c     Evaporation
283          call exf_set_evap   ( mycurrenttime, mycurrentiter, mythid )
284  #endif  #endif
285    
286  #else  #ifdef ALLOW_DOWNWARD_RADIATION
 c     Zonal wind stress.  
       call exf_set_ustress( mycurrenttime, mycurrentiter, mythid )  
287    
288  c     Meridional wind stress.  c     Downward shortwave radiation.
289        call exf_set_vstress( mycurrenttime, mycurrentiter, mythid )        call exf_set_swdown ( mycurrenttime, mycurrentiter, mythid )
   
 #endif /* ALLOW_ATM_WIND */  
   
 #else /* ALLOW_BULKFORMULAE undefined */  
 c     Atmospheric heat flux.  
       call exf_set_hflux( mycurrenttime, mycurrentiter, mythid )  
   
 c     Salt flux.  
       call exf_set_sflux( mycurrenttime, mycurrentiter, mythid )  
290    
291  c     Zonal wind stress.  c     Downward longwave radiation.
292        call exf_set_ustress( mycurrenttime, mycurrentiter, mythid )        call exf_set_lwdown ( mycurrenttime, mycurrentiter, mythid )
   
 c     Meridional wind stress.  
       call exf_set_vstress( mycurrenttime, mycurrentiter, mythid )  
293    
 #ifdef ALLOW_KPP  
 c     Net short wave radiative flux.  
       call exf_set_swflux( mycurrenttime, mycurrentiter, mythid )  
294  #endif  #endif
295    
296  #endif /* ALLOW_BULKFORMULAE */  c--   Use atmospheric state to compute surface fluxes.
297    
298  c     Loop over tiles.  c     Loop over tiles.
299  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
300  C--   HPF directive to help TAMC  C--   HPF directive to help TAMC
301  CHPF$ INDEPENDENT  CHPF$ INDEPENDENT
302  #endif /* ALLOW_AUTODIFF_TAMC */  #endif
303        do bj = mybylo(mythid),mybyhi(mythid)        do bj = mybylo(mythid),mybyhi(mythid)
304  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
305  C--    HPF directive to help TAMC  C--    HPF directive to help TAMC
# Line 319  CHPF$  INDEPENDENT Line 309  CHPF$  INDEPENDENT
309    
310            k = 1            k = 1
311    
312            do j = 1-oly,sny+oly            do j = 1,sny
313              do i = 1-olx,snx+olx              do i = 1,snx
314    
315  #ifdef ALLOW_BULKFORMULAE  #ifdef ALLOW_BULKFORMULAE
316    
# Line 341  CHPF$  INDEPENDENT Line 331  CHPF$  INDEPENDENT
331       &              + sNx*sNy*max1*max2*max3*act4       &              + sNx*sNy*max1*max2*max3*act4
332  #endif  #endif
333    
334  c             Compute the turbulent surface fluxes.  #ifdef ALLOW_DOWNWARD_RADIATION
335  c                  (Bulk formulae estimates)  c--   Compute net from downward and downward from net longwave and
336    c     shortwave radiation, if needed.
337    c     lwflux = Stefan-Boltzman constant * emissivity * SST - lwdown
338    c     swflux = - ( 1 - albedo ) * swdown
339    
340    #ifdef ALLOW_ATM_TEMP
341          if ( lwfluxfile .EQ. ' ' .AND. lwdownfile .NE. ' ' )
342         &     lwflux(i,j,bi,bj) = 5.5 _d -08 *
343         &              ((theta(i,j,k,bi,bj)+cen2kel)**4)
344         &              - lwdown(i,j,bi,bj)
345          if ( lwfluxfile .NE. ' ' .AND. lwdownfile .EQ. ' ' )
346         &     lwdown(i,j,bi,bj) = 5.5 _d -08 *
347         &              ((theta(i,j,k,bi,bj)+cen2kel)**4)
348         &              - lwflux(i,j,bi,bj)
349    #endif
350    
351    #if defined(ALLOW_ATM_TEMP) || defined(SHORTWAVE_HEATING)
352          if ( swfluxfile .EQ. ' ' .AND. swdownfile .NE. ' ' )
353         &     swflux(i,j,bi,bj) = -0.9 _d 0 * swdown(i,j,bi,bj)
354          if ( swfluxfile .NE. ' ' .AND. swdownfile .EQ. ' ' )
355         &     swdown(i,j,bi,bj) = -1.111111 _d 0 * swflux(i,j,bi,bj)
356    #endif
357    
358    #endif /* ALLOW_DOWNWARD_RADIATION */
359    
360    c--   Compute the turbulent surface fluxes.
361    
362  #ifdef ALLOW_ATM_WIND  #ifdef ALLOW_ATM_WIND
363  c             Wind speed and direction.  c             Wind speed and direction.
# Line 358  c             Wind speed and direction. Line 373  c             Wind speed and direction.
373                  sw = 0. _d 0                  sw = 0. _d 0
374                endif                endif
375                sh = max(us,umin)                sh = max(us,umin)
376  #else  #else  /* ifndef ALLOW_ATM_WIND */
377  #ifdef ALLOW_ATM_TEMP  #ifdef ALLOW_ATM_TEMP
378    
379  c             The variables us, sh and rdn have to be computed from  c             The variables us, sh and rdn have to be computed from
# Line 402  c             cdn(umps); ustar can be di Line 417  c             cdn(umps); ustar can be di
417    
418                sh    = max(us,umin)                sh    = max(us,umin)
419  #endif /* ALLOW_ATM_TEMP */  #endif /* ALLOW_ATM_TEMP */
420  #endif /* ALLOW_ATM_WIND */  #endif /* ifndef ALLOW_ATM_WIND */
421    
422  #ifdef ALLOW_ATM_TEMP  #ifdef ALLOW_ATM_TEMP
423    
# Line 441  CADJ STORE sh     = comlev1_exf_1, key = Line 456  CADJ STORE sh     = comlev1_exf_1, key =
456       &                  + sNx*niter_bulk*sNy*max1*act2       &                  + sNx*niter_bulk*sNy*max1*act2
457       &                  + sNx*niter_bulk*sNy*max1*max2*act3       &                  + sNx*niter_bulk*sNy*max1*max2*act3
458       &                  + sNx*niter_bulk*sNy*max1*max2*max3*act4       &                  + sNx*niter_bulk*sNy*max1*max2*max3*act4
 #endif  
459    
 #ifdef ALLOW_AUTODIFF_TAMC  
460  CADJ STORE rdn    = comlev1_exf_2, key = ikey_2  CADJ STORE rdn    = comlev1_exf_2, key = ikey_2
461  CADJ STORE ustar  = comlev1_exf_2, key = ikey_2  CADJ STORE ustar  = comlev1_exf_2, key = ikey_2
462  CADJ STORE qstar  = comlev1_exf_2, key = ikey_2  CADJ STORE qstar  = comlev1_exf_2, key = ikey_2
# Line 514  CADJ STORE sw     = comlev1_exf_1, key = Line 527  CADJ STORE sw     = comlev1_exf_1, key =
527    
528                  hs(i,j,bi,bj)      = atmcp*tau*tstar/ustar                  hs(i,j,bi,bj)      = atmcp*tau*tstar/ustar
529                  hl(i,j,bi,bj)      = flamb*tau*qstar/ustar                  hl(i,j,bi,bj)      = flamb*tau*qstar/ustar
530                  evap(i,j,bi,bj)    = tau*qstar/ustar  #ifndef EXF_READ_EVAP
531    cdm             evap(i,j,bi,bj)    = tau*qstar/ustar
532    cdm !!! need to change sign and to convert from kg/m^2/s to m/s !!!
533                    evap(i,j,bi,bj)    = -recip_rhonil*tau*qstar/ustar
534    #endif
535                  ustress(i,j,bi,bj) = tau*cw                  ustress(i,j,bi,bj) = tau*cw
536                  vstress(i,j,bi,bj) = tau*sw                  vstress(i,j,bi,bj) = tau*sw
537                else                else
# Line 525  CADJ STORE sw     = comlev1_exf_1, key = Line 542  CADJ STORE sw     = comlev1_exf_1, key =
542                  hl(i,j,bi,bj)      = 0. _d 0                  hl(i,j,bi,bj)      = 0. _d 0
543                endif                endif
544    
545  #else  #else  /* ifndef ALLOW_ATM_TEMP */
546  #ifdef ALLOW_ATM_WIND  #ifdef ALLOW_ATM_WIND
547                ustress(i,j,bi,bj) = atmrho*exf_BulkCdn(sh)*us*                ustress(i,j,bi,bj) = atmrho*exf_BulkCdn(sh)*us*
548       &                             uwind(i,j,bi,bj)       &                             uwind(i,j,bi,bj)
549                vstress(i,j,bi,bj) = atmrho*exf_BulkCdn(sh)*us*                vstress(i,j,bi,bj) = atmrho*exf_BulkCdn(sh)*us*
550       &                             vwind(i,j,bi,bj)       &                             vwind(i,j,bi,bj)
551  #endif /* ALLOW_ATM_WIND */  #endif
552  #endif /* ALLOW_ATM_TEMP */  #endif /* ifndef ALLOW_ATM_TEMP */
553              enddo              enddo
554            enddo            enddo
555          enddo          enddo
# Line 549  c             Net surface heat flux. Line 566  c             Net surface heat flux.
566                hfl = hfl - hs(i,j,bi,bj)                hfl = hfl - hs(i,j,bi,bj)
567                hfl = hfl - hl(i,j,bi,bj)                hfl = hfl - hl(i,j,bi,bj)
568                hfl = hfl + lwflux(i,j,bi,bj)                hfl = hfl + lwflux(i,j,bi,bj)
569  #ifndef ALLOW_KPP  #ifndef SHORTWAVE_HEATING
570                hfl = hfl + swflux(i,j,bi,bj)                hfl = hfl + swflux(i,j,bi,bj)
571  #endif /* ALLOW_KPP undef */  #endif
572  c             Heat flux:  c             Heat flux:
573                hflux(i,j,bi,bj) = hfl*maskc(i,j,1,bi,bj)                hflux(i,j,bi,bj) = hfl
574  c             Salt flux from Precipitation and Evaporation.  c             Salt flux from Precipitation and Evaporation.
575                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)
576  #endif /* ALLOW_ATM_TEMP */  #endif /* ALLOW_ATM_TEMP */
577    
 #else  
               hflux(i,j,bi,bj) = hflux(i,j,bi,bj)*maskc(i,j,1,bi,bj)  
               sflux(i,j,bi,bj) = sflux(i,j,bi,bj)*maskc(i,j,1,bi,bj)  
578  #endif /* ALLOW_BULKFORMULAE */  #endif /* ALLOW_BULKFORMULAE */
579    
580  #ifdef ALLOW_RUNOFF  #ifdef ALLOW_RUNOFF
581                sflux(i,j,bi,bj) = sflux(i,j,bi,bj) +                sflux(i,j,bi,bj) = sflux(i,j,bi,bj) - runoff(i,j,bi,bj)
582       &                           runoff(i,j,bi,bj)*maskc(i,j,1,bi,bj)  #endif
583  #endif /* ALLOW_RUNOFF */  
584                  hflux(i,j,bi,bj) = hflux(i,j,bi,bj)*maskc(i,j,1,bi,bj)
585                  sflux(i,j,bi,bj) = sflux(i,j,bi,bj)*maskc(i,j,1,bi,bj)
586    
587              enddo              enddo
588            enddo            enddo
# Line 576  c             Salt flux from Precipitati Line 592  c             Salt flux from Precipitati
592  c     Update the tile edges.  c     Update the tile edges.
593        _EXCH_XY_R8(hflux,   mythid)        _EXCH_XY_R8(hflux,   mythid)
594        _EXCH_XY_R8(sflux,   mythid)        _EXCH_XY_R8(sflux,   mythid)
595        _EXCH_XY_R8(ustress, mythid)  c      _EXCH_XY_R8(ustress, mythid)
596        _EXCH_XY_R8(vstress, mythid)  c      _EXCH_XY_R8(vstress, mythid)
597           CALL EXCH_UV_XY_RL(ustress, vstress, .TRUE., myThid)
598    
599  #ifdef ALLOW_KPP  #ifdef SHORTWAVE_HEATING
600        _EXCH_XY_R8(swflux, mythid)        _EXCH_XY_R8(swflux, mythid)
601  #endif /* ALLOW_KPP */  #endif
602    
603        end        end
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.10
  ViewVC Help
Powered by ViewVC 1.1.22