/[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.7 by dimitri, Sat Dec 28 10:11:11 2002 UTC revision 1.28 by heimbach, Thu Jul 28 19:52:33 2005 UTC
# Line 1  Line 1 
1  c $Header$  C
2    C $Header$
3    C $Name$
4    
5  #include "EXF_CPPOPTIONS.h"  #include "EXF_OPTIONS.h"
6    
7        subroutine exf_GetFFields( mycurrenttime, mycurrentiter, mythid )        subroutine exf_getffields( mycurrenttime, mycurrentiter, mythid )
8    
9  c     ==================================================================  c     ==================================================================
10  c     SUBROUTINE exf_GetFFields  c     SUBROUTINE exf_getffields
11  c     ==================================================================  c     ==================================================================
12  c  c
13  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.
 c       formulae that use the atmospheric state.  
14  c  c
15  c       The calculation of the bulk surface fluxes has been adapted from  c       heimbach@mit.edu, 23-May-2003 totally re-structured
16  c       the NCOM model which uses the formulae given in Large and Pond  c       5-Aug-2003: added USE_EXF_INTERPOLATION for arbitrary input grid
 c       (1981 & 1982 )  
 c  
 c  
 c       Header taken from NCOM version: ncom1.4.1  
 c       -----------------------------------------  
 c  
 c       Following procedures and coefficients in Large and Pond  
 c       (1981 ; 1982)  
 c  
 c       Output: Bulk estimates of the turbulent surface fluxes.  
 c       -------  
 c  
 c       hs  - sensible heat flux  (W/m^2), into ocean  
 c       hl  - latent   heat flux  (W/m^2), into ocean  
 c  
 c       Input:  
 c       ------  
 c  
 c       us  - mean wind speed (m/s)     at height hu (m)  
 c       th  - mean air temperature (K)  at height ht (m)  
 c       qh  - mean air humidity (kg/kg) at height hq (m)  
 c       sst - sea surface temperature (K)  
 c       tk0 - Kelvin temperature at 0 Celsius (K)  
 c  
 c       Assume 1) a neutral 10m drag coefficient =  
 c  
 c                 cdn = .0027/u10 + .000142 + .0000764 u10  
 c  
 c              2) a neutral 10m stanton number =  
 c  
 c                 ctn = .0327 sqrt(cdn), unstable  
 c                 ctn = .0180 sqrt(cdn), stable  
 c  
 c              3) a neutral 10m dalton number =  
 c  
 c                 cen = .0346 sqrt(cdn)  
 c  
 c              4) the saturation humidity of air at  
 c  
 c                 t(k) = exf_BulkqSat(t)  (kg/m^3)  
 c  
 c       Note:  1) here, tstar = <wt>/u*, and qstar = <wq>/u*.  
 c              2) wind speeds should all be above a minimum speed,  
 c                 say 0.5 m/s.  
 c              3) with optional iteration loop, niter=3, should suffice.  
 c              4) this version is for analyses inputs with hu = 10m and  
 c                 ht = hq.  
 c              5) sst enters in Celsius.  
 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  
 c  
 c       started: Christian Eckert eckert@mit.edu 27-Aug-1999  
 c  
 c       changed: Christian Eckert eckert@mit.edu 14-Jan-2000  
 c              - restructured the original version in order to have a  
 c                better interface to the MITgcmUV.  
 c  
 c              Christian Eckert eckert@mit.edu  12-Feb-2000  
 c              - Changed Routine names (package prefix: exf_)  
 c  
 c              Patrick Heimbach, heimbach@mit.edu  04-May-2000  
 c              - changed the handling of precip and sflux with respect  
 c                to CPP options ALLOW_BULKFORMULAE and ALLOW_ATM_TEMP  
 c              - included some CPP flags ALLOW_BULKFORMULAE to make  
 c                sure ALLOW_ATM_TEMP, ALLOW_ATM_WIND are used only in  
 c                conjunction with defined ALLOW_BULKFORMULAE  
 c              - statement functions discarded  
 c  
 c              Ralf.Giering@FastOpt.de 25-Mai-2000  
 c              - total rewrite using new subroutines  
 c  
 c              Detlef Stammer: include river run-off. Nov. 21, 2001  
 c  
 c              heimbach@mit.edu, 10-Jan-2002  
 c              - changes to enable field swapping  
 c  
 c              menemenlis@jpl.nasa.gov, 20-Dec-2002  
 c              - Added EXF_READ_EVAP and EXF_NO_BULK_COMPUTATIONS.  
17  c  c
18  c     ==================================================================  c     ==================================================================
19  c     SUBROUTINE exf_GetFFields  c     SUBROUTINE exf_getffields
20  c     ==================================================================  c     ==================================================================
21    
22        implicit none        implicit none
# Line 132  c     == global variables == Line 29  c     == global variables ==
29  #include "DYNVARS.h"  #include "DYNVARS.h"
30  #include "GRID.h"  #include "GRID.h"
31    
32    #include "exf_param.h"
33  #include "exf_fields.h"  #include "exf_fields.h"
34  #include "exf_constants.h"  #include "exf_constants.h"
35    
36  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF
37  #include "tamc.h"  # include "ctrl.h"
38    # include "ctrl_dummy.h"
39  #endif  #endif
40    
41  c     == routine arguments ==  c     == routine arguments ==
# Line 147  c     == routine arguments == Line 46  c     == routine arguments ==
46    
47  c     == local variables ==  c     == local variables ==
48    
49        integer bi,bj        integer i, j, bi, bj, interp_method
50        integer i,j,k        parameter(interp_method=1)
51    
52  #ifdef ALLOW_BULKFORMULAE  c     == end of interface ==
53    
54  #ifdef ALLOW_ATM_TEMP  c--   read forcing fields from files and temporal interpolation
       integer iter  
       _RL     aln  
       _RL     delq  
       _RL     deltap  
       _RL     hqol  
       _RL     htol  
       _RL     huol  
       _RL     psimh  
       _RL     psixh  
       _RL     qstar  
       _RL     rd  
       _RL     re  
       _RL     rdn  
       _RL     rh  
       _RL     ssttmp  
       _RL     ssq  
       _RL     stable  
       _RL     tstar  
       _RL     t0  
       _RL     ustar  
       _RL     uzn  
       _RL     shn  
       _RL     xsq  
       _RL     x  
       _RL     tau  
 #ifdef ALLOW_AUTODIFF_TAMC  
       integer ikey_1  
       integer ikey_2  
 #endif  
 #endif /* ALLOW_ATM_TEMP */  
55    
56        _RL     ustmp  cph-exf-print      print *, 'ph-exf --------- ----------------------------------'
       _RL     us  
       _RL     cw  
       _RL     sw  
       _RL     sh  
       _RL     hs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)  
       _RL     hl(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)  
       _RL     hfl  
   
 #endif /* ALLOW_BULKFORMULAE */  
   
 c     == external functions ==  
   
       integer  ilnblnk  
       external ilnblnk  
   
 #ifdef ALLOW_BULKFORMULAE  
       _RL       exf_BulkqSat  
       external  exf_BulkqSat  
       _RL       exf_BulkCdn  
       external  exf_BulkCdn  
       _RL       exf_BulkRhn  
       external  exf_BulkRhn  
 #endif /* ALLOW_BULKFORMULAE */  
   
 #ifndef ALLOW_ATM_WIND  
       _RL       TMP1  
       _RL       TMP2  
       _RL       TMP3  
       _RL       TMP4  
       _RL       TMP5  
 #endif  
57    
58  c     == end of interface ==  c     Zonal and meridional wind stress.
59    #ifdef USE_EXF_INTERPOLATION
60          call exf_set_uv(
61         &     ustressfile, ustressstartdate, ustressperiod,
62         &     ustressstartdate1, ustressstartdate2,
63         &     exf_inscal_ustress, ustress, ustress0, ustress1, ustressmask,
64         &     ustress_lon0, ustress_lon_inc, ustress_lat0, ustress_lat_inc,
65         &     ustress_nlon, ustress_nlat,
66         &     vstressfile, vstressstartdate, vstressperiod,
67         &     vstressstartdate1, vstressstartdate2,
68         &     exf_inscal_vstress, vstress, vstress0, vstress1, vstressmask,
69         &     vstress_lon0, vstress_lon_inc, vstress_lat0, vstress_lat_inc,
70         &     vstress_nlon, vstress_nlat,
71         &     mycurrenttime, mycurrentiter, mythid )
72    #else /* ifndef USE_EXF_INTERPOLATION */
73          call exf_set_gen(
74         &     ustressfile, ustressstartdate, ustressperiod,
75         &     ustressstartdate1, ustressstartdate2,
76         &     exf_inscal_ustress,
77         &     ustress, ustress0, ustress1, ustressmask,
78         &     mycurrenttime, mycurrentiter, mythid )
79          call exf_set_gen(
80         &     vstressfile, vstressstartdate, vstressperiod,
81         &     ustressstartdate1, ustressstartdate2,
82         &     exf_inscal_vstress,
83         &     vstress, vstress0, vstress1, vstressmask,
84         &     mycurrenttime, mycurrentiter, mythid )
85    #endif /* USE_EXF_INTERPOLATION */
86    
87  c     determine forcing field records  #ifdef ALLOW_ATM_WIND
88    
89  #ifdef EXF_READ_EVAP  c     Zonal and meridional wind.
90  c     Evaporation  #ifdef USE_EXF_INTERPOLATION
91        call exf_set_evap( mycurrenttime, mycurrentiter, mythid )        call exf_set_uv(
92  #endif EXF_READ_EVAP       &     uwindfile, uwindstartdate, uwindperiod,
93         &     uwindstartdate1, uwindstartdate2,
94         &     exf_inscal_uwind, uwind, uwind0, uwind1, uwindmask,
95         &     uwind_lon0, uwind_lon_inc, uwind_lat0, uwind_lat_inc,
96         &     uwind_nlon, uwind_nlat,
97         &     vwindfile, vwindstartdate, vwindperiod,
98         &     vwindstartdate1, vwindstartdate2,
99         &     exf_inscal_vwind, vwind, vwind0, vwind1, vwindmask,
100         &     vwind_lon0, vwind_lon_inc, vwind_lat0, vwind_lat_inc,
101         &     vwind_nlon, vwind_nlat,
102         &     mycurrenttime, mycurrentiter, mythid )
103    #else /* ifndef USE_EXF_INTERPOLATION */
104          call exf_set_gen(
105         &     uwindfile, uwindstartdate, uwindperiod,
106         &     uwindstartdate1, uwindstartdate2,
107         &     exf_inscal_uwind,
108         &     uwind, uwind0, uwind1, uwindmask,
109         &     mycurrenttime, mycurrentiter, mythid )
110          call exf_set_gen(
111         &     vwindfile, vwindstartdate, vwindperiod,
112         &     vwindstartdate1, vwindstartdate2,
113         &     exf_inscal_vwind,
114         &     vwind, vwind0, vwind1, vwindmask,
115         &     mycurrenttime, mycurrentiter, mythid )
116    #endif /* USE_EXF_INTERPOLATION */
117    
118  #ifdef ALLOW_BULKFORMULAE  #endif /* ALLOW_ATM_WIND */
119    
120    c     Atmospheric heat flux.
121          call exf_set_gen  (
122         &     hfluxfile, hfluxstartdate, hfluxperiod,
123         &     hfluxstartdate1, hfluxstartdate2,
124         &     exf_inscal_hflux,
125         &     hflux, hflux0, hflux1, hfluxmask,
126    #ifdef USE_EXF_INTERPOLATION
127         &     hflux_lon0, hflux_lon_inc, hflux_lat0, hflux_lat_inc,
128         &     hflux_nlon, hflux_nlat, xC, yC, interp_method,
129    #endif
130         &     mycurrenttime, mycurrentiter, mythid )
131    
132  #if (defined (ALLOW_ATM_TEMP) || defined (ALLOW_ATM_WIND))  c     Salt flux.
133  cph   This statement cannot be a PARAMETER statement in the header,        call exf_set_gen  (
134  cph   but must come here; it's not fortran77 standard       &     sfluxfile, sfluxstartdate, sfluxperiod,
135        aln = log(ht/zref)       &     sfluxstartdate1, sfluxstartdate2,
136         &     exf_inscal_sflux,
137         &     sflux, sflux0, sflux1, sfluxmask,
138    #ifdef USE_EXF_INTERPOLATION
139         &     sflux_lon0, sflux_lon_inc, sflux_lat0, sflux_lat_inc,
140         &     sflux_nlon, sflux_nlat, xC, yC, interp_method,
141  #endif  #endif
142         &     mycurrenttime, mycurrentiter, mythid )
143    
 c     Determine where we are in time and set counters, flags and  
 c     the linear interpolation factors accordingly.  
144  #ifdef ALLOW_ATM_TEMP  #ifdef ALLOW_ATM_TEMP
145    
146  c     Atmospheric temperature.  c     Atmospheric temperature.
147        call exf_set_atemp( mycurrenttime, mycurrentiter, mythid )        call exf_set_gen(
148         &     atempfile, atempstartdate, atempperiod,
149         &     atempstartdate1, atempstartdate2,
150         &     exf_inscal_atemp,
151         &     atemp, atemp0, atemp1, atempmask,
152    #ifdef USE_EXF_INTERPOLATION
153         &     atemp_lon0, atemp_lon_inc, atemp_lat0, atemp_lat_inc,
154         &     atemp_nlon, atemp_nlat, xC, yC, interp_method,
155    #endif
156         &     mycurrenttime, mycurrentiter, mythid )
157          do bj = mybylo(mythid),mybyhi(mythid)
158             do bi = mybxlo(mythid),mybxhi(mythid)
159                do j = 1,sny
160                   do i = 1,snx
161                      atemp(i,j,bi,bj) = atemp(i,j,bi,bj) + exf_offset_atemp
162                   enddo
163                enddo
164             enddo
165          enddo
166    
167  c     Atmospheric humidity.  c     Atmospheric humidity.
168        call exf_set_aqh( mycurrenttime, mycurrentiter, mythid )        call exf_set_gen(
169         &     aqhfile, aqhstartdate, aqhperiod,
170         &     aqhstartdate1, aqhstartdate2,
171         &     exf_inscal_aqh,
172         &     aqh, aqh0, aqh1, aqhmask,
173    #ifdef USE_EXF_INTERPOLATION
174         &     aqh_lon0, aqh_lon_inc, aqh_lat0, aqh_lat_inc,
175         &     aqh_nlon, aqh_nlat, xC, yC, interp_method,
176    #endif
177         &     mycurrenttime, mycurrentiter, mythid )
178    
179  c     Net long wave radiative flux.  c     Net long wave radiative flux.
180        call exf_set_lwflux( mycurrenttime, mycurrentiter, mythid )        call exf_set_gen(
181         &     lwfluxfile, lwfluxstartdate, lwfluxperiod,
182  c     Net short wave radiative flux.       &     lwfluxstartdate1, lwfluxstartdate2,
183        call exf_set_swflux( mycurrenttime, mycurrentiter, mythid )       &     exf_inscal_lwflux,
184         &     lwflux, lwflux0, lwflux1, lwfluxmask,
185  c     Precipitation.  #ifdef USE_EXF_INTERPOLATION
186        call exf_set_precip( mycurrenttime, mycurrentiter, mythid )       &     lwflux_lon0, lwflux_lon_inc, lwflux_lat0, lwflux_lat_inc,
187         &     lwflux_nlon, lwflux_nlat, xC, yC, interp_method,
 #ifdef ALLOW_ATEMP_CONTROL  
       call ctrl_getatemp ( mycurrenttime, mycurrentiter, mythid )  
188  #endif  #endif
189         &     mycurrenttime, mycurrentiter, mythid )
190    
191  #ifdef ALLOW_AQH_CONTROL  c     Precipitation.
192        call ctrl_getaqh ( mycurrenttime, mycurrentiter, mythid )        call exf_set_gen(
193         &     precipfile, precipstartdate, precipperiod,
194         &     precipstartdate1, precipstartdate2,
195         &     exf_inscal_precip,
196         &     precip, precip0, precip1, precipmask,
197    #ifdef USE_EXF_INTERPOLATION
198         &     precip_lon0, precip_lon_inc, precip_lat0, precip_lat_inc,
199         &     precip_nlon, precip_nlat, xC, yC, interp_method,
200  #endif  #endif
201         &     mycurrenttime, mycurrentiter, mythid )
 #else  
   
 c     Atmospheric heat flux.  
       call exf_set_hflux( mycurrenttime, mycurrentiter, mythid )  
   
 c     Salt flux.  
       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 */  
202    
203  #endif /* ALLOW_ATM_TEMP */  #endif /* ALLOW_ATM_TEMP */
204    
205  #ifdef ALLOW_ATM_WIND  #if defined(ALLOW_ATM_TEMP) || defined(SHORTWAVE_HEATING)
206  c     Zonal wind.  c     Net short wave radiative flux.
207        call exf_set_uwind( mycurrenttime, mycurrentiter, mythid )        call exf_set_gen  (
208         &     swfluxfile, swfluxstartdate, swfluxperiod,
209  c     Meridional wind.       &     swfluxstartdate1, swfluxstartdate2,
210        call exf_set_vwind( mycurrenttime, mycurrentiter, mythid )       &     exf_inscal_swflux,
211         &     swflux, swflux0, swflux1, swfluxmask,
212  #ifdef ALLOW_UWIND_CONTROL  #ifdef USE_EXF_INTERPOLATION
213        call ctrl_getuwind ( mycurrenttime, mycurrentiter, mythid )       &     swflux_lon0, swflux_lon_inc, swflux_lat0, swflux_lat_inc,
214         &     swflux_nlon, swflux_nlat, xC, yC, interp_method,
215  #endif  #endif
216         &     mycurrenttime, mycurrentiter, mythid )
 #ifdef ALLOW_VWIND_CONTROL  
       call ctrl_getvwind ( mycurrenttime, mycurrentiter, mythid )  
217  #endif  #endif
218    
219  #else  #ifdef EXF_READ_EVAP
220  c     Zonal wind stress.  c     Evaporation
221        call exf_set_ustress( mycurrenttime, mycurrentiter, mythid )        call exf_set_gen  (
222         &     evapfile, evapstartdate, evapperiod,
223  c     Meridional wind stress.       &     evapstartdate1, evapstartdate2,
224        call exf_set_vstress( mycurrenttime, mycurrentiter, mythid )       &     exf_inscal_evap,
225         &     evap, evap0, evap1, evapmask,
226  #endif /* ALLOW_ATM_WIND */  #ifdef USE_EXF_INTERPOLATION
227         &     evap_lon0, evap_lon_inc, evap_lat0, evap_lat_inc,
228  #else /* ALLOW_BULKFORMULAE undefined */       &     evap_nlon, evap_nlat, xC, yC, interp_method,
229  c     Atmospheric heat flux.  #endif
230        call exf_set_hflux( mycurrenttime, mycurrentiter, mythid )       &     mycurrenttime, mycurrentiter, mythid )
231    #endif
232  c     Salt flux.  
233        call exf_set_sflux( mycurrenttime, mycurrentiter, mythid )  #ifdef ALLOW_DOWNWARD_RADIATION
234    
235  c     Zonal wind stress.  c     Downward shortwave radiation.
236        call exf_set_ustress( mycurrenttime, mycurrentiter, mythid )        call exf_set_gen  (
237         &     swdownfile, swdownstartdate, swdownperiod,
238  c     Meridional wind stress.       &     swdownstartdate1, swdownstartdate2,
239        call exf_set_vstress( mycurrenttime, mycurrentiter, mythid )       &     exf_inscal_swdown,
240         &     swdown, swdown0, swdown1, swdownmask,
241    #ifdef USE_EXF_INTERPOLATION
242         &     swdown_lon0, swdown_lon_inc, swdown_lat0, swdown_lat_inc,
243         &     swdown_nlon, swdown_nlat, xC, yC, interp_method,
244    #endif
245         &     mycurrenttime, mycurrentiter, mythid )
246    
247    c     Downward longwave radiation.
248          call exf_set_gen  (
249         &     lwdownfile, lwdownstartdate, lwdownperiod,
250         &     lwdownstartdate1, lwdownstartdate2,
251         &     exf_inscal_lwdown,
252         &     lwdown, lwdown0, lwdown1, lwdownmask,
253    #ifdef USE_EXF_INTERPOLATION
254         &     lwdown_lon0, lwdown_lon_inc, lwdown_lat0, lwdown_lat_inc,
255         &     lwdown_nlon, lwdown_nlat, xC, yC, interp_method,
256    #endif
257         &     mycurrenttime, mycurrentiter, mythid )
258    
 #ifdef ALLOW_KPP  
 c     Net short wave radiative flux.  
       call exf_set_swflux( mycurrenttime, mycurrentiter, mythid )  
259  #endif  #endif
260    
261  #endif /* ALLOW_BULKFORMULAE */  #ifdef ATMOSPHERIC_LOADING
262    c     Atmos. pressure forcing
263  #if ~defined(EXF_NO_BULK_COMPUTATIONS) || ~defined(EXF_READ_EVAP)        call exf_set_gen  (
264  C--   Use atmospheric state to compute surace fluxes.       &     apressurefile, apressurestartdate, apressureperiod,
265         &     apressurestartdate1, apressurestartdate2,
266  c     Loop over tiles.       &     exf_inscal_apressure,
267  #ifdef ALLOW_AUTODIFF_TAMC       &     apressure, apressure0, apressure1, apressuremask,
268  C--   HPF directive to help TAMC  #ifdef USE_EXF_INTERPOLATION
269  CHPF$ INDEPENDENT       &     apressure_lon0, apressure_lon_inc, apressure_lat0,
270  #endif /* ALLOW_AUTODIFF_TAMC */       &     apressure_lat_inc, apressure_nlon, apressure_nlat, xC, yC,
271        do bj = mybylo(mythid),mybyhi(mythid)       &     interp_method,
272  #ifdef ALLOW_AUTODIFF_TAMC  #endif
273  C--    HPF directive to help TAMC       &     mycurrenttime, mycurrentiter, mythid )
 CHPF$  INDEPENDENT  
 #endif  
         do bi = mybxlo(mythid),mybxhi(mythid)  
   
           k = 1  
   
           do j = 1-oly,sny+oly  
             do i = 1-olx,snx+olx  
   
 #ifdef ALLOW_BULKFORMULAE  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
                act1 = bi - myBxLo(myThid)  
                max1 = myBxHi(myThid) - myBxLo(myThid) + 1  
                act2 = bj - myByLo(myThid)  
                max2 = myByHi(myThid) - myByLo(myThid) + 1  
                act3 = myThid - 1  
                max3 = nTx*nTy  
                act4 = ikey_dynamics - 1  
   
                ikey_1 = i  
      &              + sNx*(j-1)  
      &              + sNx*sNy*act1  
      &              + sNx*sNy*max1*act2  
      &              + sNx*sNy*max1*max2*act3  
      &              + sNx*sNy*max1*max2*max3*act4  
274  #endif  #endif
275    
276  c             Compute the turbulent surface fluxes.  c-- Control variables for atmos. state
 c                  (Bulk formulae estimates)  
   
 #ifdef ALLOW_ATM_WIND  
 c             Wind speed and direction.  
               ustmp = uwind(i,j,bi,bj)*uwind(i,j,bi,bj) +  
      &                vwind(i,j,bi,bj)*vwind(i,j,bi,bj)  
               if ( ustmp .ne. 0. _d 0 ) then  
                 us = sqrt(ustmp)  
                 cw = uwind(i,j,bi,bj)/us  
                 sw = vwind(i,j,bi,bj)/us  
               else  
                 us = 0. _d 0  
                 cw = 0. _d 0  
                 sw = 0. _d 0  
               endif  
               sh = max(us,umin)  
 #else  
 #ifdef ALLOW_ATM_TEMP  
   
 c             The variables us, sh and rdn have to be computed from  
 c             given wind stresses inverting relationship for neutral  
 c             drag coeff. cdn.  
 c             The inversion is based on linear and quadratic form of  
 c             cdn(umps); ustar can be directly computed from stress;  
   
               ustmp = ustress(i,j,bi,bj)*ustress(i,j,bi,bj) +  
      &                vstress(i,j,bi,bj)*vstress(i,j,bi,bj)  
               if ( ustmp .ne. 0. _d 0 ) then  
                 ustar = sqrt(ustmp/atmrho)  
                 cw = ustress(i,j,bi,bj)/sqrt(ustmp)  
                 sw = vstress(i,j,bi,bj)/sqrt(ustmp)  
               else  
                  ustar = 0. _d 0  
                  cw    = 0. _d 0  
                  sw    = 0. _d 0  
               endif  
   
               if ( ustar .eq. 0. _d 0 ) then  
                 us = 0. _d 0  
               else if ( ustar .lt. ustofu11 ) then  
                 tmp1 = -cquadrag_2/cquadrag_1/2  
                 tmp2 = sqrt(tmp1*tmp1 + ustar*ustar/cquadrag_1)  
                 us   = sqrt(tmp1 + tmp2)  
               else  
                 tmp3 = clindrag_2/clindrag_1/3  
                 tmp4 = ustar*ustar/clindrag_1/2 - tmp3**3  
                 tmp5 = sqrt(ustar*ustar/clindrag_1*  
      &                      (ustar*ustar/clindrag_1/4 - tmp3**3))  
                 us   = (tmp4 + tmp5)**(1/3) +  
      &                 tmp3**2 * (tmp4 + tmp5)**(-1/3) - tmp3  
               endif  
   
               if ( us .ne. 0 ) then  
                  rdn = ustar/us  
               else  
                  rdn = 0. _d 0  
               end if  
   
               sh    = max(us,umin)  
 #endif /* ALLOW_ATM_TEMP */  
 #endif /* ALLOW_ATM_WIND */  
   
 #ifdef ALLOW_ATM_TEMP  
   
 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)  
   
               if ( atemp(i,j,bi,bj) .ne. 0. _d 0 ) then  
                 t0     = atemp(i,j,bi,bj)*  
      &                   (exf_one + humid_fac*aqh(i,j,bi,bj))  
                 ssttmp = theta(i,j,k,bi,bj)  
                 ssq    = saltsat*  
      &                   exf_BulkqSat(ssttmp + cen2kel)/  
      &                   atmrho  
                 deltap = atemp(i,j,bi,bj)   + gamma_blk*ht -  
      &                   ssttmp - cen2kel  
                 delq   = aqh(i,j,bi,bj) - ssq  
                 stable = exf_half + sign(exf_half, deltap)  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE sh     = comlev1_exf_1, key = ikey_1  
 #endif  
                 rdn    = sqrt(exf_BulkCdn(sh))  
                 ustar  = rdn*sh  
                 tstar  = exf_BulkRhn(stable)*deltap  
                 qstar  = cdalton*delq  
   
                 do iter = 1,niter_bulk  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
                    ikey_2 = iter  
      &                  + niter_bulk*(i-1)  
      &                  + sNx*niter_bulk*(j-1)  
      &                  + sNx*niter_bulk*sNy*act1  
      &                  + sNx*niter_bulk*sNy*max1*act2  
      &                  + sNx*niter_bulk*sNy*max1*max2*act3  
      &                  + sNx*niter_bulk*sNy*max1*max2*max3*act4  
 #endif  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE rdn    = comlev1_exf_2, key = ikey_2  
 CADJ STORE ustar  = comlev1_exf_2, key = ikey_2  
 CADJ STORE qstar  = comlev1_exf_2, key = ikey_2  
 CADJ STORE tstar  = comlev1_exf_2, key = ikey_2  
 CADJ STORE sh     = comlev1_exf_2, key = ikey_2  
 CADJ STORE us     = comlev1_exf_2, key = ikey_2  
 #endif  
   
                   huol   = czol*(tstar/t0 +  
      &                     qstar/(exf_one/humid_fac+aqh(i,j,bi,bj)))/  
      &                           ustar**2  
                   huol   = max(huol,zolmin)  
                   stable = exf_half + sign(exf_half, huol)  
                   htol   = huol*ht/hu  
                   hqol   = huol*hq/hu  
   
 c                 Evaluate all stability functions assuming hq = ht.  
                   xsq    = max(sqrt(abs(exf_one - 16.*huol)),exf_one)  
                    x     = sqrt(xsq)  
                   psimh  = -psim_fac*huol*stable +  
      &                     (exf_one - stable)*  
      &                     log((exf_one + x*(exf_two + x))*  
      &                     (exf_one + xsq)/8.) - exf_two*atan(x) +  
      &                     pi*exf_half  
                   xsq    = max(sqrt(abs(exf_one - 16.*htol)),exf_one)  
                   psixh  = -psim_fac*htol*stable + (exf_one - stable)*  
      &                     exf_two*log((exf_one + xsq)/exf_two)  
   
 c                 Shift wind speed using old coefficient  
 ccc                  rd   = rdn/(exf_one + rdn/karman*  
 ccc     &                 (log(hu/zref) - psimh) )  
                   rd   = rdn/(exf_one - rdn/karman*psimh )  
                   shn  = sh*rd/rdn  
                   uzn  = max(shn, umin)  
   
 c                 Update the transfer coefficients at 10 meters  
 c                 and neutral stability.  
   
                   rdn = sqrt(exf_BulkCdn(uzn))  
   
 c                 Shift all coefficients to the measurement height  
 c                 and stability.  
 c                 rd = rdn/(exf_one + rdn/karman*(log(hu/zref) - psimh))  
                   rd = rdn/(exf_one - rdn/karman*psimh)  
                   rh = exf_BulkRhn(stable)/(exf_one +  
      &                                      exf_BulkRhn(stable)/  
      &                                     karman*(aln - psixh))  
                   re = cdalton/(exf_one + cdalton/karman*(aln - psixh))  
   
 c                 Update ustar, tstar, qstar using updated, shifted  
 c                 coefficients.  
                   ustar = rd*sh    
                   qstar = re*delq  
                   tstar = rh*deltap  
                   tau   = atmrho*ustar**2  
                   tau   = tau*us/sh  
   
                 enddo  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE ustar  = comlev1_exf_1, key = ikey_1  
 CADJ STORE qstar  = comlev1_exf_1, key = ikey_1  
 CADJ STORE tstar  = comlev1_exf_1, key = ikey_1  
 CADJ STORE tau    = comlev1_exf_1, key = ikey_1  
 CADJ STORE cw     = comlev1_exf_1, key = ikey_1  
 CADJ STORE sw     = comlev1_exf_1, key = ikey_1  
 #endif  
   
                 hs(i,j,bi,bj)      = atmcp*tau*tstar/ustar  
                 hl(i,j,bi,bj)      = flamb*tau*qstar/ustar  
 #ifndef EXF_READ_EVAP  
                 evap(i,j,bi,bj)    = tau*qstar/ustar  
 #endif EXF_READ_EVAP  
                 ustress(i,j,bi,bj) = tau*cw  
                 vstress(i,j,bi,bj) = tau*sw  
               else  
                 ustress(i,j,bi,bj) = 0. _d 0  
                 vstress(i,j,bi,bj) = 0. _d 0  
                 hflux  (i,j,bi,bj) = 0. _d 0  
                 hs(i,j,bi,bj)      = 0. _d 0  
                 hl(i,j,bi,bj)      = 0. _d 0  
               endif  
277    
278  #else  #ifdef ALLOW_ATEMP_CONTROL
279  #ifdef ALLOW_ATM_WIND        call ctrl_get_gen (
280                ustress(i,j,bi,bj) = atmrho*exf_BulkCdn(sh)*us*       &     xx_atemp_file, xx_atempstartdate, xx_atempperiod,
281       &                             uwind(i,j,bi,bj)       &     maskc, atemp, xx_atemp0, xx_atemp1, xx_atemp_dummy,
282                vstress(i,j,bi,bj) = atmrho*exf_BulkCdn(sh)*us*       &     mycurrenttime, mycurrentiter, mythid )
283       &                             vwind(i,j,bi,bj)  #endif
 #endif /* ALLOW_ATM_WIND */  
 #endif /* ALLOW_ATM_TEMP */  
             enddo  
           enddo  
         enddo  
       enddo  
   
 c     Add all contributions.  
       do bj = mybylo(mythid),mybyhi(mythid)  
         do bi = mybxlo(mythid),mybxhi(mythid)  
           do j = 1,sny  
             do i = 1,snx  
 c             Net surface heat flux.  
 #ifdef ALLOW_ATM_TEMP  
               hfl = 0. _d 0  
               hfl = hfl - hs(i,j,bi,bj)  
               hfl = hfl - hl(i,j,bi,bj)  
               hfl = hfl + lwflux(i,j,bi,bj)  
 #ifndef ALLOW_KPP  
               hfl = hfl + swflux(i,j,bi,bj)  
 #endif /* ALLOW_KPP undef */  
 c             Heat flux:  
               hflux(i,j,bi,bj) = hfl*maskc(i,j,1,bi,bj)  
 c             Salt flux from Precipitation and Evaporation.  
               sflux(i,j,bi,bj) = precip(i,j,bi,bj) - evap(i,j,bi,bj)  
 #endif /* ALLOW_ATM_TEMP */  
284    
285  #else  #ifdef ALLOW_AQH_CONTROL
286                hflux(i,j,bi,bj) = hflux(i,j,bi,bj)*maskc(i,j,1,bi,bj)        call ctrl_get_gen (
287                sflux(i,j,bi,bj) = sflux(i,j,bi,bj)*maskc(i,j,1,bi,bj)       &     xx_aqh_file, xx_aqhstartdate, xx_aqhperiod,
288  #endif /* ALLOW_BULKFORMULAE */       &     maskc, aqh, xx_aqh0, xx_aqh1, xx_aqh_dummy,
289         &     mycurrenttime, mycurrentiter, mythid )
290  #ifdef ALLOW_RUNOFF  #endif
291                sflux(i,j,bi,bj) = sflux(i,j,bi,bj) +  
292       &                           runoff(i,j,bi,bj)*maskc(i,j,1,bi,bj)  #ifdef ALLOW_PRECIP_CONTROL
293  #endif /* ALLOW_RUNOFF */        call ctrl_get_gen (
294         &     xx_precip_file, xx_precipstartdate, xx_precipperiod,
295         &     maskc, precip, xx_precip0, xx_precip1, xx_precip_dummy,
296         &     mycurrenttime, mycurrentiter, mythid )
297    #endif
298    
299    #ifdef ALLOW_SWFLUX_CONTROL
300          call ctrl_get_gen (
301         &     xx_swflux_file, xx_swfluxstartdate, xx_swfluxperiod,
302         &     maskc, swflux, xx_swflux0, xx_swflux1, xx_swflux_dummy,
303         &     mycurrenttime, mycurrentiter, mythid )
304    #endif
305    
306              enddo  #ifdef ALLOW_UWIND_CONTROL
307            enddo        call ctrl_get_gen (
308          enddo       &     xx_uwind_file, xx_uwindstartdate, xx_uwindperiod,
309        enddo       &     maskc, uwind, xx_uwind0, xx_uwind1, xx_uwind_dummy,
310         &     mycurrenttime, mycurrentiter, mythid )
311    #endif /* ALLOW_UWIND_CONTROL */
312    
313  #endif EXF_NO_BULK_COMPUTATIONS  #ifdef ALLOW_VWIND_CONTROL
314          call ctrl_get_gen (
315         &     xx_vwind_file, xx_vwindstartdate, xx_vwindperiod,
316         &     maskc, vwind, xx_vwind0, xx_vwind1, xx_vwind_dummy,
317         &     mycurrenttime, mycurrentiter, mythid )
318    #endif /* ALLOW_VWIND_CONTROL */
319    
320    #ifdef ALLOW_LWFLUX_CONTROL
321          call ctrl_get_gen (
322    NOT YET IMPLEMENTED
323         &     mytime, myiter, mythid )
324    #endif
325    
 c     Update the tile edges.  
       _EXCH_XY_R8(hflux,   mythid)  
       _EXCH_XY_R8(sflux,   mythid)  
 c      _EXCH_XY_R8(ustress, mythid)  
 c      _EXCH_XY_R8(vstress, mythid)  
        CALL EXCH_UV_XY_RS(ustress, vstress, .TRUE., myThid)  
   
 #ifdef ALLOW_KPP  
       _EXCH_XY_R8(swflux, mythid)  
 #endif /* ALLOW_KPP */  
326    
327        end        end
Legend:
Removed from v.1.7  
changed lines
  Added in v.1.28
  ViewVC Help
Powered by ViewVC 1.1.22