/[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.6.2 by heimbach, Thu Apr 4 11:08:03 2002 UTC revision 1.2.6.3 by dimitri, Thu Feb 13 19:28:38 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 175  c     == local variables == Line 165  c     == local variables ==
165        _RL     xsq        _RL     xsq
166        _RL     x        _RL     x
167        _RL     tau        _RL     tau
       _RL     evap(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)  
168  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
169        integer ikey_1        integer ikey_1
170        integer ikey_2        integer ikey_2
# Line 207  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
   
 #if (defined (ALLOW_ATM_TEMP) || defined (ALLOW_ATM_WIND))  
210  cph   This statement cannot be a PARAMETER statement in the header,  cph   This statement cannot be a PARAMETER statement in the header,
211  cph   but must come here; it's not fortran77 standard  cph   but must come here; it's not fortran77 standard
212        aln = log(ht/zref)        aln = log(ht/zref)
213  #endif  #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( mycurrenttime, mycurrentiter, mythid )  
   
 c     Atmospheric humidity.  
       call exf_set_aqh( mycurrenttime, mycurrentiter, mythid )  
216    
217  c     Net long wave radiative flux.  #ifdef ALLOW_ATM_WIND
       call exf_set_lwflux( mycurrenttime, mycurrentiter, mythid )  
218    
219  c     Net short wave radiative flux.  c     Zonal wind.
220        call exf_set_swflux( mycurrenttime, mycurrentiter, mythid )        call exf_set_uwind  ( mycurrenttime, mycurrentiter, mythid )
221    
222  c     Precipitation.  c     Meridional wind.
223        call exf_set_precip( mycurrenttime, mycurrentiter, mythid )        call exf_set_vwind  ( mycurrenttime, mycurrentiter, mythid )
224    
225  #ifdef ALLOW_ATEMP_CONTROL  #ifdef ALLOW_UWIND_CONTROL
226        call ctrl_getatemp ( mycurrenttime, mycurrentiter, mythid )        call ctrl_getuwind  ( mycurrenttime, mycurrentiter, mythid )
227  #endif  #endif
228    
229  #ifdef ALLOW_AQH_CONTROL  #ifdef ALLOW_VWIND_CONTROL
230        call ctrl_getaqh ( mycurrenttime, mycurrentiter, mythid )        call ctrl_getvwind  ( mycurrenttime, mycurrentiter, mythid )
231  #endif  #endif
232    
233  #else  #else  /* ifndef ALLOW_ATM_WIND */
234    
235  c     Atmospheric heat flux.  c     Zonal wind stress.
236        call exf_set_hflux( mycurrenttime, mycurrentiter, mythid )        call exf_set_ustress( mycurrenttime, mycurrentiter, mythid )
237    
238  c     Salt flux.  c     Meridional wind stress.
239        call exf_set_sflux( mycurrenttime, mycurrentiter, mythid )        call exf_set_vstress( mycurrenttime, mycurrentiter, mythid )
240    
241  #ifdef ALLOW_KPP  #endif /* ifndef ALLOW_ATM_WIND */
 c     Net short wave radiative flux.  
       call exf_set_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( mycurrenttime, mycurrentiter, mythid )  
250    
251  c     Meridional wind.  c     Net long wave radiative flux.
252        call exf_set_vwind( mycurrenttime, mycurrentiter, mythid )        call exf_set_lwflux ( mycurrenttime, mycurrentiter, mythid )
253    
254  #ifdef ALLOW_UWIND_CONTROL  c     Precipitation.
255        call ctrl_getuwind ( mycurrenttime, mycurrentiter, mythid )        call exf_set_precip ( mycurrenttime, mycurrentiter, mythid )
 #endif  
256    
257  #ifdef ALLOW_VWIND_CONTROL  #ifdef ALLOW_ATEMP_CONTROL
258        call ctrl_getvwind ( mycurrenttime, mycurrentiter, mythid )        call ctrl_getatemp  ( mycurrenttime, mycurrentiter, mythid )
259  #endif  #endif
260    
261  #else  #ifdef ALLOW_AQH_CONTROL
262  c     Zonal wind stress.        call ctrl_getaqh    ( mycurrenttime, mycurrentiter, mythid )
263        call exf_set_ustress( mycurrenttime, mycurrentiter, mythid )  #endif
   
 c     Meridional wind stress.  
       call exf_set_vstress( mycurrenttime, mycurrentiter, mythid )  
264    
265  #endif /* ALLOW_ATM_WIND */  #else  /* ifndef ALLOW_ATM_TEMP */
266    
 #else /* ALLOW_BULKFORMULAE undefined */  
267  c     Atmospheric heat flux.  c     Atmospheric heat flux.
268        call exf_set_hflux( mycurrenttime, mycurrentiter, mythid )        call exf_set_hflux  ( mycurrenttime, mycurrentiter, mythid )
269    
270  c     Salt flux.  c     Salt flux.
271        call exf_set_sflux( mycurrenttime, mycurrentiter, mythid )        call exf_set_sflux  ( mycurrenttime, mycurrentiter, mythid )
272    
273  c     Zonal wind stress.  #endif /* ifndef ALLOW_ATM_TEMP */
       call exf_set_ustress( mycurrenttime, mycurrentiter, mythid )  
274    
275  c     Meridional wind stress.  #if defined(ALLOW_ATM_TEMP) || defined(SHORTWAVE_HEATING)
       call exf_set_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( mycurrenttime, mycurrentiter, mythid )        call exf_set_swflux ( mycurrenttime, mycurrentiter, mythid )
278  #endif  #endif
279    
280  #endif /* ALLOW_BULKFORMULAE */  #ifdef EXF_READ_EVAP
281    c     Evaporation
282          call exf_set_evap   ( mycurrenttime, mycurrentiter, mythid )
283    #endif EXF_READ_EVAP
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  #ifdef ALLOW_AUTODIFF_TAMC
299  C--   HPF directive to help TAMC  C--   HPF directive to help TAMC
300  CHPF$ INDEPENDENT  CHPF$ INDEPENDENT
301  #endif /* ALLOW_AUTODIFF_TAMC */  #endif
302        do bj = mybylo(mythid),mybyhi(mythid)        do bj = mybylo(mythid),mybyhi(mythid)
303  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
304  C--    HPF directive to help TAMC  C--    HPF directive to help TAMC
# Line 319  CHPF$  INDEPENDENT Line 308  CHPF$  INDEPENDENT
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    
# 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 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.
# Line 358  c             Wind speed and direction. Line 357  c             Wind speed and direction.
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  c             The variables us, sh and rdn have to be computed from
# Line 402  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    
# Line 441  CADJ STORE sh     = comlev1_exf_1, key = Line 440  CADJ STORE sh     = comlev1_exf_1, key =
440       &                  + sNx*niter_bulk*sNy*max1*act2       &                  + sNx*niter_bulk*sNy*max1*act2
441       &                  + sNx*niter_bulk*sNy*max1*max2*act3       &                  + sNx*niter_bulk*sNy*max1*max2*act3
442       &                  + sNx*niter_bulk*sNy*max1*max2*max3*act4       &                  + sNx*niter_bulk*sNy*max1*max2*max3*act4
 #endif  
443    
 #ifdef ALLOW_AUTODIFF_TAMC  
444  CADJ STORE rdn    = comlev1_exf_2, key = ikey_2  CADJ STORE rdn    = comlev1_exf_2, key = ikey_2
445  CADJ STORE ustar  = comlev1_exf_2, key = ikey_2  CADJ STORE ustar  = comlev1_exf_2, key = ikey_2
446  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 511  CADJ STORE sw     = comlev1_exf_1, key =
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                  evap(i,j,bi,bj)    = tau*qstar/ustar  #ifndef EXF_READ_EVAP
515    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
521                else                else
# Line 525  CADJ STORE sw     = comlev1_exf_1, key = Line 526  CADJ STORE sw     = comlev1_exf_1, key =
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
# Line 549  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,1,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,1,bi,bj)  
               sflux(i,j,bi,bj) = sflux(i,j,bi,bj)*maskc(i,j,1,bi,bj)  
562  #endif /* ALLOW_BULKFORMULAE */  #endif /* ALLOW_BULKFORMULAE */
563    
564  #ifdef ALLOW_RUNOFF  #ifdef ALLOW_RUNOFF
565                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)
566       &                           runoff(i,j,bi,bj)*maskc(i,j,1,bi,bj)  #endif
567  #endif /* ALLOW_RUNOFF */  
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
# Line 579  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.6.2  
changed lines
  Added in v.1.2.6.3
  ViewVC Help
Powered by ViewVC 1.1.22