/[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.4.1 by heimbach, Mon Apr 8 20:10:38 2002 UTC revision 1.67 by gforget, Fri Jan 20 20:45:58 2017 UTC
# Line 1  Line 1 
1  c $Header$  C $Header$
2    C $Name$
3    
4  #include "EXF_CPPOPTIONS.h"  #include "EXF_OPTIONS.h"
5    #ifdef ALLOW_CTRL
6    # include "CTRL_OPTIONS.h"
7    #endif
8    #ifdef ALLOW_ECCO
9    # include "ECCO_OPTIONS.h"
10    #endif
11    
12          SUBROUTINE EXF_GETFFIELDS( myTime, myIter, myThid )
13    
14    C     ==================================================================
15    C     SUBROUTINE exf_getffields
16    C     ==================================================================
17    C
18    C     o Read-in atmospheric state and/or surface fluxes from files.
19    C
20    C       heimbach@mit.edu, 23-May-2003 totally re-structured
21    C       5-Aug-2003: added USE_EXF_INTERPOLATION for arbitrary input grid
22    C
23    C     ==================================================================
24    C     SUBROUTINE exf_getffields
25    C     ==================================================================
26    
27        subroutine exf_GetFFields( mycurrenttime, mycurrentiter, mythid )        IMPLICIT NONE
28    
29  c     ==================================================================  C     == global variables ==
 c     SUBROUTINE exf_GetFFields  
 c     ==================================================================  
 c  
 c     o Get the surface fluxes either from file or as derived from bulk  
 c       formulae that use the atmospheric state.  
 c  
 c       The calculation of the bulk surface fluxes has been adapted from  
 c       the NCOM model which uses the formulae given in Large and Pond  
 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  
 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  
 c              - Changed Routine names (package prefix: exf_)  
 c  
 c              Patrick Heimbach, heimbach@mit.edu  04-May-2000  
 c  
 c              - changed the handling of precip and sflux with respect  
 c                to CPP options ALLOW_BULKFORMULAE and ALLOW_ATM_TEMP  
 c  
 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  
 c              - statement functions discarded  
 c  
 c              Ralf.Giering@FastOpt.de 25-Mai-2000  
 c  
 c              - total rewrite using new subroutines  
 c  
 c     ==================================================================  
 c     SUBROUTINE exf_GetFFields  
 c     ==================================================================  
   
       implicit none  
   
 c     == global variables ==  
30    
31  #include "EEPARAMS.h"  #include "EEPARAMS.h"
32  #include "SIZE.h"  #include "SIZE.h"
# Line 130  c     == global variables == Line 34  c     == global variables ==
34  #include "DYNVARS.h"  #include "DYNVARS.h"
35  #include "GRID.h"  #include "GRID.h"
36    
37  #include "exf_fields.h"  #include "EXF_PARAM.h"
38  #include "exf_constants.h"  #include "EXF_FIELDS.h"
39    #include "EXF_CONSTANTS.h"
40  #ifdef ALLOW_AUTODIFF_TAMC  
41  #include "tamc.h"  #ifdef ALLOW_CTRL
42    # include "CTRL_SIZE.h"
43    # include "ctrl.h"
44    # include "ctrl_dummy.h"
45    # ifdef ALLOW_GENTIM2D_CONTROL
46    #  include "CTRL_GENARR.h"
47    # endif
48    #endif
49    #if (defined (ALLOW_ECCO) && defined (ECCO_CTRL_DEPRECATED))
50    #  include "ecco_cost.h"
51    #endif
52    
53    C     == routine arguments ==
54          _RL     myTime
55          INTEGER myIter
56          INTEGER myThid
57    
58    C     == local variables ==
59          INTEGER i, j, bi, bj
60    #ifdef ALLOW_ROTATE_UV_CONTROLS
61          _RL     tmpUE(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nsx,nsy)
62          _RL     tmpVN(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nsx,nsy)
63          _RL     tmpUX(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nsx,nsy)
64          _RL     tmpVY(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nsx,nsy)
65    #endif
66    #if (defined (ALLOW_CTRL) && \
67         defined (ALLOW_GENTIM2D_CONTROL))
68          INTEGER iarr
69    #endif
70    
71    C     == end of interface ==
72    
73    C--   read forcing fields from files and temporal interpolation
74    
75    C     Zonal and meridional wind stress.
76          IF ( .NOT.useAtmWind ) THEN
77          CALL EXF_SET_UV(
78         I     ustressfile, ustressstartdate, ustressperiod,
79         I     exf_inscal_ustress,
80         I     ustress_exfremo_intercept, ustress_exfremo_slope,
81         U     ustress, ustress0, ustress1, ustressmask,
82         I     vstressfile, vstressstartdate, vstressperiod,
83         I     exf_inscal_vstress,
84         I     vstress_exfremo_intercept, vstress_exfremo_slope,
85         U     vstress, vstress0, vstress1, vstressmask,
86    #ifdef USE_EXF_INTERPOLATION
87         I     ustress_lon0, ustress_lon_inc, ustress_lat0, ustress_lat_inc,
88         I     ustress_nlon, ustress_nlat, ustress_interpMethod,
89         I     vstress_lon0, vstress_lon_inc, vstress_lat0, vstress_lat_inc,
90         I     vstress_nlon, vstress_nlat, vstress_interpMethod,
91         I     uvInterp_stress,
92    #endif /* USE_EXF_INTERPOLATION */
93         I     myTime, myIter, myThid )
94          ELSE
95           DO bj = myByLo(myThid),myByHi(myThid)
96            DO bi = myBxLo(myThid),mybxhi(myThid)
97             DO j=1-OLy,sNy+OLy
98              DO i=1-OLx,sNx+OLx
99               ustress(i,j,bi,bj) = 0. _d 0
100               vstress(i,j,bi,bj) = 0. _d 0
101              ENDDO
102             ENDDO
103            ENDDO
104           ENDDO
105          ENDIF
106    
107    C--   wind speed
108          CALL EXF_SET_GEN(
109         &     wspeedfile, wspeedstartdate, wspeedperiod,
110         &     exf_inscal_wspeed,
111         &     wspeed_exfremo_intercept, wspeed_exfremo_slope,
112         &     wspeed, wspeed0, wspeed1, wspeedmask,
113    #ifdef USE_EXF_INTERPOLATION
114         &     wspeed_lon0, wspeed_lon_inc,
115         &     wspeed_lat0, wspeed_lat_inc,
116         &     wspeed_nlon, wspeed_nlat, xC, yC, wspeed_interpMethod,
117    #endif
118         &     myTime, myIter, myThid )
119    
120    C     Zonal and meridional wind.
121          IF ( useAtmWind ) THEN
122          CALL EXF_SET_UV(
123         I     uwindfile, uwindstartdate, uwindperiod,
124         I     exf_inscal_uwind,
125         I     uwind_exfremo_intercept, uwind_exfremo_slope,
126         U     uwind, uwind0, uwind1, uwindmask,
127         I     vwindfile, vwindstartdate, vwindperiod,
128         I     exf_inscal_vwind,
129         I     vwind_exfremo_intercept, vwind_exfremo_slope,
130         U     vwind, vwind0, vwind1, vwindmask,
131    #ifdef USE_EXF_INTERPOLATION
132         I     uwind_lon0, uwind_lon_inc, uwind_lat0, uwind_lat_inc,
133         I     uwind_nlon, uwind_nlat, uwind_interpMethod,
134         I     vwind_lon0, vwind_lon_inc, vwind_lat0, vwind_lat_inc,
135         I     vwind_nlon, vwind_nlat, vwind_interpMethod, uvInterp_wind,
136    #endif /* USE_EXF_INTERPOLATION */
137         I     myTime, myIter, myThid )
138    
139          if (useRelativeWind) then
140    C     Subtract UVEL and VVEL from UWIND and VWIND.
141           DO bj = myByLo(myThid),myByHi(myThid)
142            DO bi = myBxLo(myThid),mybxhi(myThid)
143             DO j = 1,sNy
144              DO i = 1,sNx
145               uwind(i,j,bi,bj) = uwind(i,j,bi,bj) - 0.5 _d 0 *
146         &          (uVel(i,j,1,bi,bj)+uVel(i+1,j,1,bi,bj))
147               vwind(i,j,bi,bj) = vwind(i,j,bi,bj) - 0.5 _d 0 *
148         &          (vVel(i,j,1,bi,bj)+vVel(i,j+1,1,bi,bj))
149              ENDDO
150             ENDDO
151            ENDDO
152           ENDDO
153          endif
154    
155          ELSE
156           DO bj = myByLo(myThid),myByHi(myThid)
157            DO bi = myBxLo(myThid),mybxhi(myThid)
158             DO j=1-OLy,sNy+OLy
159              DO i=1-OLx,sNx+OLx
160               uwind(i,j,bi,bj) = 0. _d 0
161               vwind(i,j,bi,bj) = 0. _d 0
162              ENDDO
163             ENDDO
164            ENDDO
165           ENDDO
166          ENDIF
167    
168    C     Atmospheric heat flux.
169          CALL EXF_SET_GEN  (
170         &     hfluxfile, hfluxstartdate, hfluxperiod,
171         &     exf_inscal_hflux,
172         &     hflux_exfremo_intercept, hflux_exfremo_slope,
173         &     hflux, hflux0, hflux1, hfluxmask,
174    #ifdef USE_EXF_INTERPOLATION
175         &     hflux_lon0, hflux_lon_inc, hflux_lat0, hflux_lat_inc,
176         &     hflux_nlon, hflux_nlat, xC, yC, hflux_interpMethod,
177    #endif
178         &     myTime, myIter, myThid )
179    
180    C     Freshwater flux.
181          CALL EXF_SET_GEN  (
182         &     sfluxfile, sfluxstartdate, sfluxperiod,
183         &     exf_inscal_sflux,
184         &     sflux_exfremo_intercept, sflux_exfremo_slope,
185         &     sflux, sflux0, sflux1, sfluxmask,
186    #ifdef USE_EXF_INTERPOLATION
187         &     sflux_lon0, sflux_lon_inc, sflux_lat0, sflux_lat_inc,
188         &     sflux_nlon, sflux_nlat, xC, yC, sflux_interpMethod,
189  #endif  #endif
190         &     myTime, myIter, myThid )
 c     == routine arguments ==  
   
       integer mythid  
       integer mycurrentiter  
       _RL     mycurrenttime  
   
 c     == local variables ==  
   
       integer bi,bj  
       integer i,j,k  
   
 #ifdef ALLOW_BULKFORMULAE  
191    
192  #ifdef ALLOW_ATM_TEMP  #ifdef ALLOW_ATM_TEMP
       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  
       _RL     evap(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)  
 #ifdef ALLOW_AUTODIFF_TAMC  
       integer ikey_1  
       integer ikey_2  
 #endif  
 #endif /* ALLOW_ATM_TEMP */  
   
       _RL     ustmp  
       _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 ==  
193    
194        integer  ilnblnk  C     Atmospheric temperature.
195        external ilnblnk        CALL EXF_SET_GEN(
196         &     atempfile, atempstartdate, atempperiod,
197         &     exf_inscal_atemp,
198         &     atemp_exfremo_intercept, atemp_exfremo_slope,
199         &     atemp, atemp0, atemp1, atempmask,
200    #ifdef USE_EXF_INTERPOLATION
201         &     atemp_lon0, atemp_lon_inc, atemp_lat0, atemp_lat_inc,
202         &     atemp_nlon, atemp_nlat, xC, yC, atemp_interpMethod,
203    #endif
204         &     myTime, myIter, myThid )
205          DO bj = myByLo(myThid),myByHi(myThid)
206           DO bi = myBxLo(myThid),mybxhi(myThid)
207            DO j = 1,sNy
208             DO i = 1,sNx
209              atemp(i,j,bi,bj) = atemp(i,j,bi,bj) + exf_offset_atemp
210             ENDDO
211            ENDDO
212           ENDDO
213          ENDDO
214    
215    C     Atmospheric humidity.
216          CALL EXF_SET_GEN(
217         &     aqhfile, aqhstartdate, aqhperiod,
218         &     exf_inscal_aqh,
219         &     aqh_exfremo_intercept, aqh_exfremo_slope,
220         &     aqh, aqh0, aqh1, aqhmask,
221    #ifdef USE_EXF_INTERPOLATION
222         &     aqh_lon0, aqh_lon_inc, aqh_lat0, aqh_lat_inc,
223         &     aqh_nlon, aqh_nlat, xC, yC, aqh_interpMethod,
224    #endif
225         &     myTime, myIter, myThid )
226    
227    C     Net long wave radiative flux.
228          CALL EXF_SET_GEN(
229         &     lwfluxfile, lwfluxstartdate, lwfluxperiod,
230         &     exf_inscal_lwflux,
231         &     lwflux_exfremo_intercept, lwflux_exfremo_slope,
232         &     lwflux, lwflux0, lwflux1, lwfluxmask,
233    #ifdef USE_EXF_INTERPOLATION
234         &     lwflux_lon0, lwflux_lon_inc, lwflux_lat0, lwflux_lat_inc,
235         &     lwflux_nlon, lwflux_nlat, xC, yC, lwflux_interpMethod,
236    #endif
237         &     myTime, myIter, myThid )
238    
239    C     Precipitation.
240          CALL EXF_SET_GEN(
241         &     precipfile, precipstartdate, precipperiod,
242         &     exf_inscal_precip,
243         &     precip_exfremo_intercept, precip_exfremo_slope,
244         &     precip, precip0, precip1, precipmask,
245    #ifdef USE_EXF_INTERPOLATION
246         &     precip_lon0, precip_lon_inc, precip_lat0, precip_lat_inc,
247         &     precip_nlon, precip_nlat, xC, yC, precip_interpMethod,
248    #endif
249         &     myTime, myIter, myThid )
250    
251    C     Snow.
252          CALL EXF_SET_GEN(
253         &     snowprecipfile, snowprecipstartdate, snowprecipperiod,
254         &     exf_inscal_snowprecip,
255         &     snowprecip_exfremo_intercept, snowprecip_exfremo_slope,
256         &     snowprecip, snowprecip0, snowprecip1, snowprecipmask,
257    #ifdef USE_EXF_INTERPOLATION
258         &     snowprecip_lon0, snowprecip_lon_inc,
259         &     snowprecip_lat0, snowprecip_lat_inc,
260         &     snowprecip_nlon, snowprecip_nlat, xC, yC,
261         &     snowprecip_interpMethod,
262    #endif
263         &     myTime, myIter, myThid )
264    C     Take care of case where total precip is not defined
265          IF ( snowPrecipFile .NE. ' ' ) THEN
266           DO bj = myByLo(myThid),myByHi(myThid)
267            DO bi = myBxLo(myThid),mybxhi(myThid)
268             DO j = 1,sNy
269              DO i = 1,sNx
270               precip(i,j,bi,bj) =
271         &          max( precip(i,j,bi,bj), snowPrecip(i,j,bi,bj) )
272              ENDDO
273             ENDDO
274            ENDDO
275           ENDDO
276          ENDIF
277    
278  #ifdef ALLOW_BULKFORMULAE  #endif /* ALLOW_ATM_TEMP */
       _RL       exf_BulkqSat  
       external  exf_BulkqSat  
       _RL       exf_BulkCdn  
       external  exf_BulkCdn  
       _RL       exf_BulkRhn  
       external  exf_BulkRhn  
 #endif /* ALLOW_BULKFORMULAE */  
   
 c     == end of interface ==  
   
 c     determine forcing field records  
   
 #ifdef ALLOW_BULKFORMULAE  
   
 #if (defined (ALLOW_ATM_TEMP) || defined (ALLOW_ATM_WIND))  
       aln = log(ht/zref)  
 #endif  
   
 c     Determine where we are in time and set counters, flags and  
 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 )  
   
 c     Net long wave radiative flux.  
       call exf_set_lwflux( mycurrenttime, mycurrentiter, mythid )  
   
 c     Net short wave radiative flux.  
       call exf_set_swflux( mycurrenttime, mycurrentiter, mythid )  
279    
280  c     Precipitation.  #if defined(ALLOW_ATM_TEMP) || defined(SHORTWAVE_HEATING)
281        call exf_set_precip( mycurrenttime, mycurrentiter, mythid )  C     Net short wave radiative flux.
282          CALL EXF_SET_GEN  (
283         &     swfluxfile, swfluxstartdate, swfluxperiod,
284         &     exf_inscal_swflux,
285         &     swflux_exfremo_intercept, swflux_exfremo_slope,
286         &     swflux, swflux0, swflux1, swfluxmask,
287    #ifdef USE_EXF_INTERPOLATION
288         &     swflux_lon0, swflux_lon_inc, swflux_lat0, swflux_lat_inc,
289         &     swflux_nlon, swflux_nlat, xC, yC, swflux_interpMethod,
290    #endif
291         &     myTime, myIter, myThid )
292    #endif
293    
294    #ifdef EXF_READ_EVAP
295    C     Evaporation
296          CALL EXF_SET_GEN  (
297         &     evapfile, evapstartdate, evapperiod,
298         &     exf_inscal_evap,
299         &     evap_exfremo_intercept, evap_exfremo_slope,
300         &     evap, evap0, evap1, evapmask,
301    #ifdef USE_EXF_INTERPOLATION
302         &     evap_lon0, evap_lon_inc, evap_lat0, evap_lat_inc,
303         &     evap_nlon, evap_nlat, xC, yC, evap_interpMethod,
304    #endif
305         &     myTime, myIter, myThid )
306    #endif
307    
308    #ifdef ALLOW_DOWNWARD_RADIATION
309    
310    C     Downward shortwave radiation.
311          CALL EXF_SET_GEN  (
312         &     swdownfile, swdownstartdate, swdownperiod,
313         &     exf_inscal_swdown,
314         &     swdown_exfremo_intercept, swdown_exfremo_slope,
315         &     swdown, swdown0, swdown1, swdownmask,
316    #ifdef USE_EXF_INTERPOLATION
317         &     swdown_lon0, swdown_lon_inc, swdown_lat0, swdown_lat_inc,
318         &     swdown_nlon, swdown_nlat, xC, yC, swdown_interpMethod,
319    #endif
320         &     myTime, myIter, myThid )
321    
322    C     Downward longwave radiation.
323          CALL EXF_SET_GEN  (
324         &     lwdownfile, lwdownstartdate, lwdownperiod,
325         &     exf_inscal_lwdown,
326         &     lwdown_exfremo_intercept, lwdown_exfremo_slope,
327         &     lwdown, lwdown0, lwdown1, lwdownmask,
328    #ifdef USE_EXF_INTERPOLATION
329         &     lwdown_lon0, lwdown_lon_inc, lwdown_lat0, lwdown_lat_inc,
330         &     lwdown_nlon, lwdown_nlat, xC, yC, lwdown_interpMethod,
331    #endif
332         &     myTime, myIter, myThid )
333    
334    #endif
335    
336    #ifdef ATMOSPHERIC_LOADING
337    C     Atmos. pressure forcing
338          CALL EXF_SET_GEN  (
339         &     apressurefile, apressurestartdate, apressureperiod,
340         &     exf_inscal_apressure,
341         &     apressure_exfremo_intercept, apressure_exfremo_slope,
342         &     apressure, apressure0, apressure1, apressuremask,
343    #ifdef USE_EXF_INTERPOLATION
344         &     apressure_lon0, apressure_lon_inc,
345         &     apressure_lat0, apressure_lat_inc,
346         &     apressure_nlon,apressure_nlat,xC,yC, apressure_interpMethod,
347    #endif
348         &     myTime, myIter, myThid )
349    #endif
350    
351    #ifdef EXF_SEAICE_FRACTION
352    C     fractional ice-covered area MASK
353          CALL EXF_SET_GEN  (
354         &     areamaskfile, areamaskstartdate, areamaskperiod,
355         &     exf_inscal_areamask,
356         &     areamask_exfremo_intercept, areamask_exfremo_slope,
357         &     areamask, areamask0, areamask1, areamaskmask,
358    #ifdef USE_EXF_INTERPOLATION
359         &     areamask_lon0, areamask_lon_inc,
360         &     areamask_lat0, areamask_lat_inc,
361         &     areamask_nlon, areamask_nlat, xC, yC, areamask_interpMethod,
362    #endif
363         &     myTime, myIter, myThid )
364    #endif
365    
366    #ifdef ALLOW_RUNOFF
367    C     Runoff
368          CALL EXF_SET_GEN  (
369         &     runofffile, runoffstartdate, runoffperiod,
370         &     exf_inscal_runoff,
371         &     runoff_exfremo_intercept, runoff_exfremo_slope,
372         &     runoff, runoff0, runoff1, runoffmask,
373    #ifdef USE_EXF_INTERPOLATION
374         &     runoff_lon0, runoff_lon_inc, runoff_lat0, runoff_lat_inc,
375         &     runoff_nlon, runoff_nlat, xC, yC, runoff_interpMethod,
376    #endif
377         &     myTime, myIter, myThid )
378    #endif /* ALLOW_RUNOFF */
379    
380    #ifdef ALLOW_RUNOFTEMP
381    C     Runoff temperature
382          CALL EXF_SET_GEN  (
383         &     runoftempfile, runoffstartdate, runoffperiod,
384         &     exf_inscal_runoftemp,
385         &     runoftemp_exfremo_intercept, runoftemp_exfremo_slope,
386         &     runoftemp, runoftemp0, runoftemp1, runoffmask,
387    #ifdef USE_EXF_INTERPOLATION
388         &     runoff_lon0, runoff_lon_inc, runoff_lat0, runoff_lat_inc,
389         &     runoff_nlon, runoff_nlat, xC, yC, runoff_interpMethod,
390    #endif
391         &     myTime, myIter, myThid )
392    #endif /* ALLOW_RUNOFTEMP */
393    
394    #ifdef ALLOW_SALTFLX
395    C     Salt flux
396          CALL EXF_SET_GEN  (
397         &     saltflxfile, saltflxstartdate, saltflxperiod,
398         &     exf_inscal_saltflx,
399         &     saltflx_exfremo_intercept, saltflx_exfremo_slope,
400         &     saltflx, saltflx0, saltflx1, saltflxmask,
401    #ifdef USE_EXF_INTERPOLATION
402         &     saltflx_lon0, saltflx_lon_inc,
403         &     saltflx_lat0, saltflx_lat_inc,
404         &     saltflx_nlon, saltflx_nlat, xC, yC, saltflx_interpMethod,
405    #endif
406         &     myTime, myIter, myThid )
407    #endif
408    
409    #ifdef ALLOW_ROTATE_UV_CONTROLS
410          if ( useCTRL ) then
411            DO bj = myByLo(myThid),myByHi(myThid)
412             DO bi = myBxLo(myThid),mybxhi(myThid)
413              DO j = 1-OLy,sNy+OLy
414               DO i = 1-OLx,sNx+OLx
415                 tmpUE(i,j,bi,bj) = 0. _d 0
416                 tmpVN(i,j,bi,bj) = 0. _d 0
417                 tmpUX(i,j,bi,bj) = 0. _d 0
418                 tmpVY(i,j,bi,bj) = 0. _d 0
419               ENDDO
420              ENDDO
421             ENDDO
422            ENDDO
423          endif
424    #endif
425    
426    # if (!defined (ALLOW_ECCO) || defined (ECCO_CTRL_DEPRECATED))
427    
428    C-- Control variables for atmos. state
429    #ifdef ALLOW_CTRL
430          if (.NOT.ctrlUseGen) then
431    
432  #ifdef ALLOW_ATEMP_CONTROL  #ifdef ALLOW_ATEMP_CONTROL
433        call ctrl_getatemp ( mycurrenttime, mycurrentiter, mythid )        CALL CTRL_GET_GEN (
434         &     xx_atemp_file, xx_atempstartdate, xx_atempperiod,
435         &     maskc, atemp, xx_atemp0, xx_atemp1, xx_atemp_dummy,
436         &     xx_atemp_remo_intercept, xx_atemp_remo_slope,
437         &     watemp, myTime, myIter, myThid )
438  #endif  #endif
439    
440  #ifdef ALLOW_AQH_CONTROL  #ifdef ALLOW_AQH_CONTROL
441        call ctrl_getaqh ( mycurrenttime, mycurrentiter, mythid )        CALL CTRL_GET_GEN (
442         &     xx_aqh_file, xx_aqhstartdate, xx_aqhperiod,
443         &     maskc, aqh, xx_aqh0, xx_aqh1, xx_aqh_dummy,
444         &     xx_aqh_remo_intercept, xx_aqh_remo_slope,
445         &     waqh, myTime, myIter, myThid )
446    #endif
447    
448    #ifdef ALLOW_PRECIP_CONTROL
449          CALL CTRL_GET_GEN (
450         &     xx_precip_file, xx_precipstartdate, xx_precipperiod,
451         &     maskc, precip, xx_precip0, xx_precip1, xx_precip_dummy,
452         &     xx_precip_remo_intercept, xx_precip_remo_slope,
453         &     wprecip, myTime, myIter, myThid )
454    #endif
455    
456    #ifdef ALLOW_SWDOWN_CONTROL
457          CALL CTRL_GET_GEN (
458         &     xx_swdown_file, xx_swdownstartdate, xx_swdownperiod,
459         &     maskc, swdown, xx_swdown0, xx_swdown1, xx_swdown_dummy,
460         &     xx_swdown_remo_intercept, xx_swdown_remo_slope,
461         &     wswdown, myTime, myIter, myThid )
462    #endif
463    
464    #ifdef ALLOW_LWDOWN_CONTROL
465          CALL CTRL_GET_GEN (
466         &     xx_lwdown_file, xx_lwdownstartdate, xx_lwdownperiod,
467         &     maskc, lwdown, xx_lwdown0, xx_lwdown1, xx_lwdown_dummy,
468         &     xx_lwdown_remo_intercept, xx_lwdown_remo_slope,
469         &     wlwdown, myTime, myIter, myThid )
470    #endif
471    
472          endif !if (.NOT.ctrlUseGen) then
473    
474    #ifdef ALLOW_SWFLUX_CONTROL
475          CALL CTRL_GET_GEN (
476         &     xx_swflux_file, xx_swfluxstartdate, xx_swfluxperiod,
477         &     maskc, swflux, xx_swflux0, xx_swflux1, xx_swflux_dummy,
478         &     xx_swflux_remo_intercept, xx_swflux_remo_slope,
479         &     wswflux, myTime, myIter, myThid )
480    #endif
481    
482    #ifdef ALLOW_LWFLUX_CONTROL
483          CALL CTRL_GET_GEN (
484         &     xx_lwflux_file, xx_lwfluxstartdate, xx_lwfluxperiod,
485         &     maskc, lwflux, xx_lwflux0, xx_lwflux1, xx_lwflux_dummy,
486         &     xx_lwflux_remo_intercept, xx_lwflux_remo_slope,
487         &     wswflux, myTime, myIter, myThid )
488    #endif
489    
490    #ifdef ALLOW_EVAP_CONTROL
491          CALL CTRL_GET_GEN (
492         &     xx_evap_file, xx_evapstartdate, xx_evapperiod,
493         &     maskc, evap, xx_evap0, xx_evap1, xx_evap_dummy,
494         &     xx_evap_remo_intercept, xx_evap_remo_slope,
495         &     wevap, myTime, myIter, myThid )
496    #endif
497    
498    #ifdef ALLOW_SNOWPRECIP_CONTROL
499          CALL CTRL_GET_GEN (
500         &     xx_snowprecip_file, xx_snowprecipstartdate,
501         &     xx_snowprecipperiod,
502         &     maskc, snowprecip, xx_snowprecip0, xx_snowprecip1,
503         &     xx_snowprecip_dummy,
504         &     xx_snowprecip_remo_intercept, xx_snowprecip_remo_slope,
505         &     wsnowprecip, myTime, myIter, myThid )
506    #endif
507    
508    #ifdef ALLOW_APRESSURE_CONTROL
509          CALL CTRL_GET_GEN (
510         &     xx_apressure_file, xx_apressurestartdate,
511         &     xx_apressureperiod,
512         &     maskc, apressure, xx_apressure0, xx_apressure1,
513         &     xx_apressure_dummy,
514         &     xx_apressure_remo_intercept, xx_apressure_remo_slope,
515         &     wapressure, myTime, myIter, myThid )
516  #endif  #endif
517    
518  #else        IF ( useAtmWind ) THEN
519    #ifndef ALLOW_ROTATE_UV_CONTROLS
 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 */  
   
 #endif /* 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 )  
520    
521  #ifdef ALLOW_UWIND_CONTROL  #ifdef ALLOW_UWIND_CONTROL
522        call ctrl_getuwind ( mycurrenttime, mycurrentiter, mythid )        CALL CTRL_GET_GEN (
523  #endif       &     xx_uwind_file, xx_uwindstartdate, xx_uwindperiod,
524         &     maskc, uwind, xx_uwind0, xx_uwind1, xx_uwind_dummy,
525         &     xx_uwind_remo_intercept, xx_uwind_remo_slope,
526         &     wuwind, myTime, myIter, myThid )
527    #endif /* ALLOW_UWIND_CONTROL */
528    
529  #ifdef ALLOW_VWIND_CONTROL  #ifdef ALLOW_VWIND_CONTROL
530        call ctrl_getvwind ( mycurrenttime, mycurrentiter, mythid )        CALL CTRL_GET_GEN (
531  #endif       &     xx_vwind_file, xx_vwindstartdate, xx_vwindperiod,
532         &     maskc, vwind, xx_vwind0, xx_vwind1, xx_vwind_dummy,
533         &     xx_vwind_remo_intercept, xx_vwind_remo_slope,
534         &     wvwind, myTime, myIter, myThid )
535    #endif /* ALLOW_VWIND_CONTROL */
536    
537  #else  #else
 c     Zonal wind stress.  
       call exf_set_ustress( mycurrenttime, mycurrentiter, mythid )  
538    
539  c     Meridional wind stress.  #if defined(ALLOW_UWIND_CONTROL) && defined(ALLOW_VWIND_CONTROL)
540        call exf_set_vstress( mycurrenttime, mycurrentiter, mythid )        CALL CTRL_GET_GEN (
541         &     xx_uwind_file, xx_uwindstartdate, xx_uwindperiod,
542  #endif /* ALLOW_ATM_WIND */       &     maskc, tmpUE, xx_uwind0, xx_uwind1, xx_uwind_dummy,
543         &     xx_uwind_remo_intercept, xx_uwind_remo_slope,
544  #else /* ALLOW_BULKFORMULAE undefined */       &     wuwind, myTime, myIter, myThid )
545  c     Atmospheric heat flux.  
546        call exf_set_hflux( mycurrenttime, mycurrentiter, mythid )        CALL CTRL_GET_GEN (
547         &     xx_vwind_file, xx_vwindstartdate, xx_vwindperiod,
548  c     Salt flux.       &     maskc, tmpVN, xx_vwind0, xx_vwind1, xx_vwind_dummy,
549        call exf_set_sflux( mycurrenttime, mycurrentiter, mythid )       &     xx_vwind_remo_intercept, xx_vwind_remo_slope,
550         &     wvwind, myTime, myIter, myThid )
551  c     Zonal wind stress.  
552        call exf_set_ustress( mycurrenttime, mycurrentiter, mythid )        CALL ROTATE_UV2EN_RL(tmpUX,tmpVY,tmpUE,tmpVN,
553         &     .FALSE.,.FALSE.,.TRUE.,1,myThid)
554  c     Meridional wind stress.  
555        call exf_set_vstress( mycurrenttime, mycurrentiter, mythid )          DO bj = myByLo(myThid),myByHi(myThid)
556             DO bi = myBxLo(myThid),mybxhi(myThid)
557  #ifdef ALLOW_KPP            DO j = 1,sNy
558  c     Net short wave radiative flux.             DO i = 1,sNx
559        call exf_set_swflux( mycurrenttime, mycurrentiter, mythid )               uwind(i,j,bi,bj)=uwind(i,j,bi,bj)+tmpUX(i,j,bi,bj)
560  #endif               vwind(i,j,bi,bj)=vwind(i,j,bi,bj)+tmpVY(i,j,bi,bj)
561               ENDDO
562  #endif /* ALLOW_BULKFORMULAE */            ENDDO
563             ENDDO
564  c     Loop over tiles.          ENDDO
565  #ifdef ALLOW_AUTODIFF_TAMC  #endif
566  C--   HPF directive to help TAMC  
567  CHPF$ INDEPENDENT  #endif /* ALLOW_ROTATE_UV_CONTROLS */
568  #endif /* ALLOW_AUTODIFF_TAMC */        ENDIF
569        do bj = mybylo(mythid),mybyhi(mythid)  
570  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_ATM_MEAN_CONTROL
571  C--    HPF directive to help TAMC        DO bj = myByLo(myThid),myByHi(myThid)
572  CHPF$  INDEPENDENT         DO bi = myBxLo(myThid),mybxhi(myThid)
573  #endif          DO j = 1,sNy
574          do bi = mybxlo(mythid),mybxhi(mythid)           DO i = 1,sNx
575    # ifdef ALLOW_ATEMP_CONTROL
576            k = 1            atemp(i,j,bi,bj) =atemp(i,j,bi,bj) +xx_atemp_mean(i,j,bi,bj)
577    # endif
578            do j = 1-oly,sny+oly  # ifdef ALLOW_AQH_CONTROL
579              do i = 1-olx,snx+olx            aqh(i,j,bi,bj)   =aqh(i,j,bi,bj)   +xx_aqh_mean(i,j,bi,bj)
580    # endif
581  #ifdef ALLOW_BULKFORMULAE  # ifdef ALLOW_PRECIP_CONTROL
582              precip(i,j,bi,bj)=precip(i,j,bi,bj)+xx_precip_mean(i,j,bi,bj)
583  #ifdef ALLOW_AUTODIFF_TAMC  # endif
584                 act1 = bi - myBxLo(myThid)  # ifdef ALLOW_SWDOWN_CONTROL
585                 max1 = myBxHi(myThid) - myBxLo(myThid) + 1            swdown(i,j,bi,bj)=swdown(i,j,bi,bj)+xx_swdown_mean(i,j,bi,bj)
586                 act2 = bj - myByLo(myThid)  # endif
587                 max2 = myByHi(myThid) - myByLo(myThid) + 1  # ifdef ALLOW_UWIND_CONTROL
588                 act3 = myThid - 1            uwind(i,j,bi,bj) =uwind(i,j,bi,bj) +xx_uwind_mean(i,j,bi,bj)
589                 max3 = nTx*nTy  # endif
590                 act4 = ikey_dynamics - 1  # ifdef ALLOW_VWIND_CONTROL
591              vwind(i,j,bi,bj) =vwind(i,j,bi,bj) +xx_vwind_mean(i,j,bi,bj)
592                 ikey_1 = i  # endif
593       &              + sNx*(j-1)           ENDDO
594       &              + sNx*sNy*act1          ENDDO
595       &              + sNx*sNy*max1*act2         ENDDO
596       &              + sNx*sNy*max1*max2*act3        ENDDO
597       &              + sNx*sNy*max1*max2*max3*act4  #endif /* ALLOW_ATM_MEAN_CONTROL */
598  #endif  
599    cdm transferred from exf_init_runoff.F
600  c             Compute the turbulent surface fluxes.  cdm functionality needs to be checked before turning on
601  c                  (Bulk formulae estimates)  cdm #ifdef ALLOW_RUNOFF_CONTROL
602    cdm       CALL CTRL_GET_GEN (
603  #ifdef ALLOW_ATM_WIND  cdm      &     xx_runoff_file, xx_runoffstartdate, xx_runoffperiod,
604  c             Wind speed and direction.  cdm      &     maskc, runoff, xx_runoff0, xx_runoff1, xx_runoff_dummy,
605                ustmp = uwind(i,j,bi,bj)*uwind(i,j,bi,bj) +  cdm      &     xx_runoff_remo_intercept, xx_runoff_remo_slope,
606       &                vwind(i,j,bi,bj)*vwind(i,j,bi,bj)  cdm      &     wrunoff, 0., 0., myThid )
607                if ( ustmp .ne. 0. _d 0 ) then  cdm #endif
608                  us = sqrt(ustmp)  
609                  cw = uwind(i,j,bi,bj)/us  #endif /* ALLOW_CTRL */
610                  sw = vwind(i,j,bi,bj)/us  
611                else  #endif /* undef ALLOW_ECCO) || def ECCO_CTRL_DEPRECATED */
612                  us = 0. _d 0  
613                  cw = 0. _d 0  #if (defined (ALLOW_CTRL) && defined (ALLOW_GENTIM2D_CONTROL))
614                  sw = 0. _d 0        if ( useCTRL.AND.ctrlUseGen ) then
615                endif        DO bj = myByLo(myThid),myByHi(myThid)
616                sh = max(us,umin)         DO bi = myBxLo(myThid),mybxhi(myThid)
617  #else          DO j = 1,sNy
618             DO i = 1,sNx
619               do iarr = 1, maxCtrlTim2D
620  #ifdef ALLOW_ATM_TEMP  #ifdef ALLOW_ATM_TEMP
621               if (xx_gentim2d_file(iarr)(1:8).EQ.'xx_atemp')
622  c             The variables us, sh and rdn have to be computed from       &       atemp(i,j,bi,bj)=atemp(i,j,bi,bj)+
623  c             given wind stresses inverting relationship for neutral       &                         xx_gentim2d(i,j,bi,bj,iarr)
624  c             drag coeff. cdn.             if (xx_gentim2d_file(iarr)(1:6).EQ.'xx_aqh')
625  c             The inversion is based on linear and quadratic form of       &       aqh(i,j,bi,bj)=aqh(i,j,bi,bj)+
626  c             cdn(umps); ustar can be directly computed from stress;       &                         xx_gentim2d(i,j,bi,bj,iarr)
627               if (xx_gentim2d_file(iarr)(1:9).EQ.'xx_precip')
628                ustmp = ustress(i,j,bi,bj)*ustress(i,j,bi,bj) +       &       precip(i,j,bi,bj)=precip(i,j,bi,bj)+
629       &                vstress(i,j,bi,bj)*vstress(i,j,bi,bj)       &                         xx_gentim2d(i,j,bi,bj,iarr)
630                if ( ustmp .ne. 0. _d 0 ) then             if (xx_gentim2d_file(iarr)(1:9).EQ.'xx_lwflux')
631                  ustar = sqrt(ustmp/atmrho)       &       lwflux(i,j,bi,bj)=lwflux(i,j,bi,bj)+
632                  cw = ustress(i,j,bi,bj)/sqrt(ustmp)       &                         xx_gentim2d(i,j,bi,bj,iarr)
633                  sw = vstress(i,j,bi,bj)/sqrt(ustmp)  #endif
634                else  #if defined(ALLOW_ATM_TEMP) || defined(SHORTWAVE_HEATING)
635                   ustar = 0. _d 0             if (xx_gentim2d_file(iarr)(1:9).EQ.'xx_swflux')
636                   cw    = 0. _d 0       &       swflux(i,j,bi,bj)=swflux(i,j,bi,bj)+
637                   sw    = 0. _d 0       &                         xx_gentim2d(i,j,bi,bj,iarr)
638                endif  #endif
639    #ifdef ALLOW_DOWNWARD_RADIATION
640                if ( ustar .eq. 0. _d 0 ) then             if (xx_gentim2d_file(iarr)(1:9).EQ.'xx_swdown')
641                  us = 0. _d 0       &       swdown(i,j,bi,bj)=swdown(i,j,bi,bj)+
642                else if ( ustar .lt. ustofu11 ) then       &                         xx_gentim2d(i,j,bi,bj,iarr)
643                  tmp1 = -cquadrag_2/cquadrag_1/2             if (xx_gentim2d_file(iarr)(1:9).EQ.'xx_lwdown')
644                  tmp2 = sqrt(tmp1*tmp1 + ustar*ustar/cquadrag_1)       &       lwdown(i,j,bi,bj)=lwdown(i,j,bi,bj)+
645                  us   = sqrt(tmp1 + tmp2)       &                         xx_gentim2d(i,j,bi,bj,iarr)
646                else  #endif
647                  tmp3 = clindrag_2/clindrag_1/3  #ifdef ALLOW_RUNOFF
648                  tmp4 = ustar*ustar/clindrag_1/2 - tmp3**3             if (xx_gentim2d_file(iarr)(1:9).EQ.'xx_runoff')
649                  tmp5 = sqrt(ustar*ustar/clindrag_1*       &       runoff(i,j,bi,bj)=runoff(i,j,bi,bj)+
650       &                      (ustar*ustar/clindrag_1/4 - tmp3**3))       &                         xx_gentim2d(i,j,bi,bj,iarr)
651                  us   = (tmp4 + tmp5)**(1/3) +  #endif
652       &                 tmp3**2 * (tmp4 + tmp5)**(-1/3) - tmp3  #ifdef EXF_READ_EVAP
653                endif             if (xx_gentim2d_file(iarr)(1:7).EQ.'xx_evap')
654         &       evap(i,j,bi,bj)=evap(i,j,bi,bj)+
655                if ( us .ne. 0 ) then       &                         xx_gentim2d(i,j,bi,bj,iarr)
656                   rdn = ustar/us  #endif
657                else  #ifdef ATMOSPHERIC_LOADING
658                   rdn = 0. _d 0             if (xx_gentim2d_file(iarr)(1:12).EQ.'xx_apressure')
659                end if       &       apressure(i,j,bi,bj)=apressure(i,j,bi,bj)+
660         &                         xx_gentim2d(i,j,bi,bj,iarr)
661                sh    = max(us,umin)  #endif
662  #endif /* ALLOW_ATM_TEMP */  #ifdef EXF_SEAICE_FRACTION
663  #endif /* ALLOW_ATM_WIND */             if (xx_gentim2d_file(iarr)(1:11).EQ.'xx_areamask')
664         &       areamask(i,j,bi,bj)=areamask(i,j,bi,bj)+
665  #ifdef ALLOW_ATM_TEMP       &                         xx_gentim2d(i,j,bi,bj,iarr)
666    #endif
667  c             Initial guess: z/l=0.0; hu=ht=hq=z  #ifndef ALLOW_ROTATE_UV_CONTROLS
668  c             Iterations:    converge on z/l and hence the fluxes.             if (xx_gentim2d_file(iarr)(1:8).EQ.'xx_uwind')
669  c             t0     : virtual temperature (K)       &       uwind(i,j,bi,bj)=uwind(i,j,bi,bj)+
670  c             ssq    : sea surface humidity (kg/kg)       &                         xx_gentim2d(i,j,bi,bj,iarr)
671  c             deltap : potential temperature diff (K)             if (xx_gentim2d_file(iarr)(1:8).EQ.'xx_vwind')
672         &       vwind(i,j,bi,bj)=vwind(i,j,bi,bj)+
673                if ( atemp(i,j,bi,bj) .ne. 0. _d 0 ) then       &                         xx_gentim2d(i,j,bi,bj,iarr)
                 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  
   
                 evap(i,j,bi,bj)    = tau*qstar/ustar  
   
                 ustress(i,j,bi,bj) = tau*cw  
                 vstress(i,j,bi,bj) = tau*sw  
 ce              hflux(i,j,bi,bj)   = atmcp*tau*tstar/ustar +  
 ce     &                             flamb*tau*qstar/ustar  
               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  
   
674  #else  #else
675  #ifdef ALLOW_ATM_WIND             if (xx_gentim2d_file(iarr)(1:8).EQ.'xx_uwind')
676                ustress(i,j,bi,bj) = atmrho*exf_BulkCdn(sh)*us*       &       tmpUE(i,j,bi,bj)=tmpUE(i,j,bi,bj)+
677       &                             uwind(i,j,bi,bj)       &                         xx_gentim2d(i,j,bi,bj,iarr)
678                vstress(i,j,bi,bj) = atmrho*exf_BulkCdn(sh)*us*             if (xx_gentim2d_file(iarr)(1:8).EQ.'xx_vwind')
679       &                             vwind(i,j,bi,bj)       &       tmpVN(i,j,bi,bj)=tmpVN(i,j,bi,bj)+
680  #endif /* ALLOW_ATM_WIND */       &                         xx_gentim2d(i,j,bi,bj,iarr)
681  #endif /* ALLOW_ATM_TEMP */  #endif
682              enddo             enddo
683            enddo           ENDDO
684          enddo          ENDDO
685        enddo         ENDDO
686          ENDDO
687  c     Add all contributions.  #ifdef ALLOW_ROTATE_UV_CONTROLS
688        do bj = mybylo(mythid),mybyhi(mythid)        CALL ROTATE_UV2EN_RL(tmpUX,tmpVY,tmpUE,tmpVN,
689          do bi = mybxlo(mythid),mybxhi(mythid)       &     .FALSE.,.FALSE.,.TRUE.,1,myThid)
690            do j = 1,sny  
691              do i = 1,snx         DO bj = myByLo(myThid),myByHi(myThid)
692  c             Net surface heat flux.           DO bi = myBxLo(myThid),mybxhi(myThid)
693  #ifdef ALLOW_ATM_TEMP            DO j = 1,sNy
694                hfl = 0. _d 0             DO i = 1,sNx
695                hfl = hfl - hs(i,j,bi,bj)               uwind(i,j,bi,bj)=uwind(i,j,bi,bj)+tmpUX(i,j,bi,bj)
696                hfl = hfl - hl(i,j,bi,bj)               vwind(i,j,bi,bj)=vwind(i,j,bi,bj)+tmpVY(i,j,bi,bj)
697                hfl = hfl + lwflux(i,j,bi,bj)             ENDDO
698  #ifndef ALLOW_KPP            ENDDO
699                hfl = hfl + swflux(i,j,bi,bj)           ENDDO
700  #endif /* ALLOW_KPP undef */         ENDDO
701  c             Heat flux:  #endif /* ALLOW_ROTATE_UV_CONTROLS */
               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 */  
702    
703  #else        endif !if (ctrlUseGen) then
704                hflux(i,j,bi,bj) = hflux(i,j,bi,bj)*maskc(i,j,1,bi,bj)  #endif
               sflux(i,j,bi,bj) = sflux(i,j,bi,bj)*maskc(i,j,1,bi,bj)  
 #endif /* ALLOW_BULKFORMULAE */  
             enddo  
           enddo  
         enddo  
       enddo  
   
 c     Update the tile edges.  
       _EXCH_XY_R8(hflux,   mythid)  
       _EXCH_XY_R8(sflux,   mythid)  
       _EXCH_XY_R8(ustress, mythid)  
       _EXCH_XY_R8(vstress, mythid)  
   
 #ifdef ALLOW_KPP  
       _EXCH_XY_R8(swflux, mythid)  
 #endif /* ALLOW_KPP */  
705    
706        end        RETURN
707          END
Legend:
Removed from v.1.2.4.1  
changed lines
  Added in v.1.67
  ViewVC Help
Powered by ViewVC 1.1.22