/[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.8 by dimitri, Tue Feb 18 05:33:54 2003 UTC revision 1.18 by edhill, Thu Oct 23 04:41:40 2003 UTC
# Line 1  Line 1 
1  c $Header$  C
2    C $Header$
3    C $Name$
4    
5  #include "EXF_CPPOPTIONS.h"  #include "AD_CONFIG.h"
6    #include "EXF_OPTIONS.h"
7    
8        subroutine exf_GetFFields( mycurrenttime, mycurrentiter, mythid )        subroutine exf_getffields( mycurrenttime, mycurrentiter, mythid )
9    
10  c     ==================================================================  c     ==================================================================
11  c     SUBROUTINE exf_GetFFields  c     SUBROUTINE exf_getffields
12  c     ==================================================================  c     ==================================================================
13  c  c
14  c     o Read-in atmospheric state and/or surface fluxes from files.  c     o Read-in atmospheric state and/or surface fluxes from files.
15  c  c
16  c     o Use bulk formulae to estimate turbulent and/or radiative  c       heimbach@mit.edu, 23-May-2003 totally re-structured
17  c       fluxes at the surface.  c       5-Aug-2003: added USE_EXF_INTERPOLATION for arbitrary input grid
 c  
 c     NOTES:  
 c     ======  
 c  
 c     See EXF_CPPOPTIONS.h for a description of the various possible  
 c     ocean-model forcing configurations.  
 c  
 c     The bulk formulae of pkg/exf are not valid for sea-ice covered  
 c     oceans but they can be used in combination with a sea-ice model,  
 c     for example, pkg/seaice, to specify open water flux contributions.  
 c  
 c     ==================================================================  
 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       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       mods for pkg/seaice: menemenlis@jpl.nasa.gov 20-Dec-2002  
18  c  c
19  c     ==================================================================  c     ==================================================================
20  c     SUBROUTINE exf_GetFFields  c     SUBROUTINE exf_getffields
21  c     ==================================================================  c     ==================================================================
22    
23        implicit none        implicit none
# Line 119  c     == global variables == Line 30  c     == global variables ==
30  #include "DYNVARS.h"  #include "DYNVARS.h"
31  #include "GRID.h"  #include "GRID.h"
32    
33    #include "exf_param.h"
34  #include "exf_fields.h"  #include "exf_fields.h"
35  #include "exf_constants.h"  #include "exf_constants.h"
36    
37  #ifdef ALLOW_AUTODIFF_TAMC  #if (defined (ALLOW_ADJOINT_RUN) || \
38  #include "tamc.h"       defined (ALLOW_TANGENTLINEAR_RUN) || \
39         defined (ALLOW_ECCO_OPTIMIZATION))
40    # include "ctrl.h"
41    # include "ctrl_dummy.h"
42  #endif  #endif
43    
44  c     == routine arguments ==  c     == routine arguments ==
# Line 134  c     == routine arguments == Line 49  c     == routine arguments ==
49    
50  c     == local variables ==  c     == local variables ==
51    
       integer bi,bj  
       integer i,j,k  
   
 #ifdef ALLOW_BULKFORMULAE  
   
 #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  
 #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 ==  
   
       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  
   
52  c     == end of interface ==  c     == end of interface ==
53    
 #ifdef ALLOW_BULKFORMULAE  
 cph   This statement cannot be a PARAMETER statement in the header,  
 cph   but must come here; it's not fortran77 standard  
       aln = log(ht/zref)  
 #endif  
   
54  c--   read forcing fields from files and temporal interpolation  c--   read forcing fields from files and temporal interpolation
55    
56  #ifdef ALLOW_ATM_WIND  #ifdef ALLOW_ATM_WIND
57    
58  c     Zonal wind.  c     Zonal wind.
59        call exf_set_uwind  ( mycurrenttime, mycurrentiter, mythid )        call exf_set_gen(
60         &     uwindfile, uwindstartdate, uwindperiod,
61         &     exf_inscal_uwind,
62         &     uwind, uwind0, uwind1, uwindmask,
63    #ifdef USE_EXF_INTERPOLATION
64         &     uwind_lon0, uwind_lon_inc, uwind_lat0, uwind_lat_inc,
65         &     uwind_nlon, uwind_nlat, xC, yC,
66    #endif
67         &     mycurrenttime, mycurrentiter, mythid )
68    
69  c     Meridional wind.  c     Meridional wind.
70        call exf_set_vwind  ( mycurrenttime, mycurrentiter, mythid )        call exf_set_gen(
71         &     vwindfile, vwindstartdate, vwindperiod,
72         &     exf_inscal_vwind,
73         &     vwind, vwind0, vwind1, vwindmask,
74    #ifdef USE_EXF_INTERPOLATION
75         &     vwind_lon0, vwind_lon_inc, vwind_lat0, vwind_lat_inc,
76         &     vwind_nlon, vwind_nlat, xC, yC,
77    #endif
78         &     mycurrenttime, mycurrentiter, mythid )
79    
80  #ifdef ALLOW_UWIND_CONTROL  #ifdef ALLOW_UWIND_CONTROL
81        call ctrl_getuwind  ( mycurrenttime, mycurrentiter, mythid )        call ctrl_get_gen (
82         &     xx_uwind_file, xx_uwindstartdate, xx_uwindperiod,
83         &     maskw, uwind, xx_uwind0, xx_uwind1, xx_uwind_dummy,
84         &     mycurrenttime, mycurrentiter, mythid )
85  #endif  #endif
86    
87  #ifdef ALLOW_VWIND_CONTROL  #ifdef ALLOW_VWIND_CONTROL
88        call ctrl_getvwind  ( mycurrenttime, mycurrentiter, mythid )        call ctrl_get_gen (
89         &     xx_vwind_file, xx_vwindstartdate, xx_vwindperiod,
90         &     masks, vwind, xx_vwind0, xx_vwind1, xx_vwind_dummy,
91         &     mycurrenttime, mycurrentiter, mythid )
92  #endif  #endif
93    
94  #else  /* ifndef ALLOW_ATM_WIND */  #else  /* ifndef ALLOW_ATM_WIND */
95    
96  c     Zonal wind stress.  c     Zonal wind stress.
97        call exf_set_ustress( mycurrenttime, mycurrentiter, mythid )        call exf_set_gen(
98         &     ustressfile, ustressstartdate, ustressperiod,
99         &     exf_inscal_ustress,
100         &     ustress, ustress0, ustress1, ustressmask,
101    #ifdef USE_EXF_INTERPOLATION
102         &     ustress_lon0, ustress_lon_inc, ustress_lat0, ustress_lat_inc,
103         &     ustress_nlon, ustress_nlat, xG, yC,
104    #endif
105         &     mycurrenttime, mycurrentiter, mythid )
106    
107  c     Meridional wind stress.  c     Meridional wind stress.
108        call exf_set_vstress( mycurrenttime, mycurrentiter, mythid )        call exf_set_gen(
109         &     vstressfile, vstressstartdate, vstressperiod,
110         &     exf_inscal_vstress,
111         &     vstress, vstress0, vstress1, vstressmask,
112    #ifdef USE_EXF_INTERPOLATION
113         &     vstress_lon0, vstress_lon_inc, vstress_lat0, vstress_lat_inc,
114         &     vstress_nlon, vstress_nlat, xC, yG,
115    #endif
116         &     mycurrenttime, mycurrentiter, mythid )
117    
118  #endif /* ifndef ALLOW_ATM_WIND */  #endif /* ifndef ALLOW_ATM_WIND */
119    
120  #ifdef ALLOW_ATM_TEMP  #ifdef ALLOW_ATM_TEMP
121    
122  c     Atmospheric temperature.  c     Atmospheric temperature.
123        call exf_set_atemp  ( mycurrenttime, mycurrentiter, mythid )        call exf_set_gen(
124         &     atempfile, atempstartdate, atempperiod,
125         &     exf_inscal_atemp,
126         &     atemp, atemp0, atemp1, atempmask,
127    #ifdef USE_EXF_INTERPOLATION
128         &     atemp_lon0, atemp_lon_inc, atemp_lat0, atemp_lat_inc,
129         &     atemp_nlon, atemp_nlat, xC, yC,
130    #endif
131         &     mycurrenttime, mycurrentiter, mythid )
132    
133  c     Atmospheric humidity.  c     Atmospheric humidity.
134        call exf_set_aqh    ( mycurrenttime, mycurrentiter, mythid )        call exf_set_gen(
135         &     aqhfile, aqhstartdate, aqhperiod,
136         &     exf_inscal_aqh,
137         &     aqh, aqh0, aqh1, aqhmask,
138    #ifdef USE_EXF_INTERPOLATION
139         &     aqh_lon0, aqh_lon_inc, aqh_lat0, aqh_lat_inc,
140         &     aqh_nlon, aqh_nlat, xC, yC,
141    #endif
142         &     mycurrenttime, mycurrentiter, mythid )
143    
144  c     Net long wave radiative flux.  c     Net long wave radiative flux.
145        call exf_set_lwflux ( mycurrenttime, mycurrentiter, mythid )        call exf_set_gen(
146         &     lwfluxfile, lwfluxstartdate, lwfluxperiod,
147         &     exf_inscal_lwflux,
148         &     lwflux, lwflux0, lwflux1, lwfluxmask,
149    #ifdef USE_EXF_INTERPOLATION
150         &     lwflux_lon0, lwflux_lon_inc, lwflux_lat0, lwflux_lat_inc,
151         &     lwflux_nlon, lwflux_nlat, xC, yC,
152    #endif
153         &     mycurrenttime, mycurrentiter, mythid )
154    
155  c     Precipitation.  c     Precipitation.
156        call exf_set_precip ( mycurrenttime, mycurrentiter, mythid )        call exf_set_gen(
157         &     precipfile, precipstartdate, precipperiod,
158         &     exf_inscal_precip,
159         &     precip, precip0, precip1, precipmask,
160    #ifdef USE_EXF_INTERPOLATION
161         &     precip_lon0, precip_lon_inc, precip_lat0, precip_lat_inc,
162         &     precip_nlon, precip_nlat, xC, yC,
163    #endif
164         &     mycurrenttime, mycurrentiter, mythid )
165    
166  #ifdef ALLOW_ATEMP_CONTROL  #ifdef ALLOW_ATEMP_CONTROL
167        call ctrl_getatemp  ( mycurrenttime, mycurrentiter, mythid )        call ctrl_get_gen (
168         &     xx_atemp_file, xx_atempstartdate, xx_atempperiod,
169         &     maskc, atemp, xx_atemp0, xx_atemp1, xx_atemp_dummy,
170         &     mycurrenttime, mycurrentiter, mythid )
171  #endif  #endif
172    
173  #ifdef ALLOW_AQH_CONTROL  #ifdef ALLOW_AQH_CONTROL
174        call ctrl_getaqh    ( mycurrenttime, mycurrentiter, mythid )        call ctrl_get_gen (
175         &     xx_aqh_file, xx_aqhstartdate, xx_aqhperiod,
176         &     maskc, aqh, xx_aqh0, xx_aqh1, xx_aqh_dummy,
177         &     mycurrenttime, mycurrentiter, mythid )
178  #endif  #endif
179    
180  #else  /* ifndef ALLOW_ATM_TEMP */  #else  /* ifndef ALLOW_ATM_TEMP */
181    
182  c     Atmospheric heat flux.  c     Atmospheric heat flux.
183        call exf_set_hflux  ( mycurrenttime, mycurrentiter, mythid )        call exf_set_gen  (
184         &     hfluxfile, hfluxstartdate, hfluxperiod, exf_inscal_hflux,
185         &     hflux, hflux0, hflux1, hfluxmask,
186    #ifdef USE_EXF_INTERPOLATION
187         &     hflux_lon0, hflux_lon_inc, hflux_lat0, hflux_lat_inc,
188         &     hflux_nlon, hflux_nlat, xC, yC,
189    #endif
190         &     mycurrenttime, mycurrentiter, mythid )
191    
192  c     Salt flux.  c     Salt flux.
193        call exf_set_sflux  ( mycurrenttime, mycurrentiter, mythid )        call exf_set_gen  (
194         &     sfluxfile, sfluxstartdate, sfluxperiod, exf_inscal_sflux,
195         &     sflux, sflux0, sflux1, sfluxmask,
196    #ifdef USE_EXF_INTERPOLATION
197         &     sflux_lon0, sflux_lon_inc, sflux_lat0, sflux_lat_inc,
198         &     sflux_nlon, sflux_nlat, xC, yC,
199    #endif
200         &     mycurrenttime, mycurrentiter, mythid )
201    
202  #endif /* ifndef ALLOW_ATM_TEMP */  #endif /* ifndef ALLOW_ATM_TEMP */
203    
204  #if defined(ALLOW_ATM_TEMP) || defined(SHORTWAVE_HEATING)  #if defined(ALLOW_ATM_TEMP) || defined(SHORTWAVE_HEATING)
205  c     Net short wave radiative flux.  c     Net short wave radiative flux.
206        call exf_set_swflux ( mycurrenttime, mycurrentiter, mythid )        call exf_set_gen  (
207         &     swfluxfile, swfluxstartdate, swfluxperiod, exf_inscal_swflux,
208         &     swflux, swflux0, swflux1, swfluxmask,
209    #ifdef USE_EXF_INTERPOLATION
210         &     swflux_lon0, swflux_lon_inc, swflux_lat0, swflux_lat_inc,
211         &     swflux_nlon, swflux_nlat, xC, yC,
212    #endif
213         &     mycurrenttime, mycurrentiter, mythid )
214  #endif  #endif
215    
216  #ifdef EXF_READ_EVAP  #ifdef EXF_READ_EVAP
217  c     Evaporation  c     Evaporation
218        call exf_set_evap   ( mycurrenttime, mycurrentiter, mythid )        call exf_set_gen  (
219         &     evapfile, evapstartdate, evapperiod, exf_inscal_evap,
220         &     evap, evap0, evap1, evapmask,
221    #ifdef USE_EXF_INTERPOLATION
222         &     evap_lon0, evap_lon_inc, evap_lat0, evap_lat_inc,
223         &     evap_nlon, evap_nlat, xC, yC,
224    #endif
225         &     mycurrenttime, mycurrentiter, mythid )
226  #endif  #endif
227    
228  #ifdef ALLOW_DOWNWARD_RADIATION  #ifdef ALLOW_DOWNWARD_RADIATION
229    
230  c     Downward shortwave radiation.  c     Downward shortwave radiation.
231        call exf_set_swdown ( mycurrenttime, mycurrentiter, mythid )        call exf_set_gen  (
232         &     swdownfile, swdownstartdate, swdownperiod, exf_inscal_swdown,
233  c     Downward longwave radiation.       &     swdown, swdown0, swdown1, swdownmask,
234        call exf_set_lwdown ( mycurrenttime, mycurrentiter, mythid )  #ifdef USE_EXF_INTERPOLATION
235         &     swdown_lon0, swdown_lon_inc, swdown_lat0, swdown_lat_inc,
236  #endif       &     swdown_nlon, swdown_nlat, xC, yC,
   
 c--   Use atmospheric state to compute surface fluxes.  
   
 c     Loop over tiles.  
 #ifdef ALLOW_AUTODIFF_TAMC  
 C--   HPF directive to help TAMC  
 CHPF$ INDEPENDENT  
 #endif  
       do bj = mybylo(mythid),mybyhi(mythid)  
 #ifdef ALLOW_AUTODIFF_TAMC  
 C--    HPF directive to help TAMC  
 CHPF$  INDEPENDENT  
237  #endif  #endif
238          do bi = mybxlo(mythid),mybxhi(mythid)       &     mycurrenttime, mycurrentiter, mythid )
   
           k = 1  
   
 cdm? can olx, oly be eliminated?  
           do j = 1-oly,sny+oly  
             do i = 1-olx,snx+olx  
239    
240  #ifdef ALLOW_BULKFORMULAE  c     Downward longwave radiation.
241          call exf_set_gen  (
242  #ifdef ALLOW_AUTODIFF_TAMC       &     lwdownfile, lwdownstartdate, lwdownperiod, exf_inscal_lwdown,
243                 act1 = bi - myBxLo(myThid)       &     lwdown, lwdown0, lwdown1, lwdownmask,
244                 max1 = myBxHi(myThid) - myBxLo(myThid) + 1  #ifdef USE_EXF_INTERPOLATION
245                 act2 = bj - myByLo(myThid)       &     lwdown_lon0, lwdown_lon_inc, lwdown_lat0, lwdown_lat_inc,
246                 max2 = myByHi(myThid) - myByLo(myThid) + 1       &     lwdown_nlon, lwdown_nlat, xC, yC,
247                 act3 = myThid - 1  #endif
248                 max3 = nTx*nTy       &     mycurrenttime, mycurrentiter, mythid )
                act4 = ikey_dynamics - 1  
249    
                ikey_1 = i  
      &              + sNx*(j-1)  
      &              + sNx*sNy*act1  
      &              + sNx*sNy*max1*act2  
      &              + sNx*sNy*max1*max2*act3  
      &              + sNx*sNy*max1*max2*max3*act4  
250  #endif  #endif
251    
252  #ifdef ALLOW_DOWNWARD_RADIATION  #ifdef ATMOSPHERIC_LOADING
253  c--   Compute net longwave and shortwave radiation:  c     Atmos. pressure forcing
254  c     lwflux = Stefan-Boltzman constant * emissivity * SST - lwdown        call exf_set_gen  (
255  c     swflux = - ( 1 - albedo ) * swdown       &     apressurefile, apressurestartdate, apressureperiod,
256                 lwflux(i,j,bi,bj) = 5.5 _d -08 *       &     exf_inscal_apressure,
257       &              ((theta(i,j,k,bi,bj)+cen2kel)**4)       &     apressure, apressure0, apressure1, apressuremask,
258       &              - lwdown(i,j,bi,bj)  #ifdef USE_EXF_INTERPOLATION
259                 swflux(i,j,bi,bj) = -0.9 _d 0 * swdown(i,j,bi,bj)       &     apressure_lon0, apressure_lon_inc, apressure_lat0,
260         &     apressure_lat_inc, apressure_nlon, apressure_nlat, xC, yC,
261  #endif  #endif
262         &     mycurrenttime, mycurrentiter, mythid )
 c--   Compute the turbulent surface fluxes.  
   
 #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  /* ifndef ALLOW_ATM_WIND */  
 #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 /* ifndef 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  
   
 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  
 cdm             evap(i,j,bi,bj)    = tau*qstar/ustar  
 cdm !!! need to change sign and to convert from kg/m^2/s to m/s !!!  
                 evap(i,j,bi,bj)    = -recip_rhonil*tau*qstar/ustar  
 #endif  
                 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  
   
 #else  /* ifndef ALLOW_ATM_TEMP */  
 #ifdef ALLOW_ATM_WIND  
               ustress(i,j,bi,bj) = atmrho*exf_BulkCdn(sh)*us*  
      &                             uwind(i,j,bi,bj)  
               vstress(i,j,bi,bj) = atmrho*exf_BulkCdn(sh)*us*  
      &                             vwind(i,j,bi,bj)  
263  #endif  #endif
 #endif /* ifndef 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 SHORTWAVE_HEATING  
               hfl = hfl + swflux(i,j,bi,bj)  
 #endif  
 c             Heat flux:  
               hflux(i,j,bi,bj) = hfl  
 c             Salt flux from Precipitation and Evaporation.  
               sflux(i,j,bi,bj) = evap(i,j,bi,bj) - precip(i,j,bi,bj)  
 #endif /* ALLOW_ATM_TEMP */  
   
 #endif /* ALLOW_BULKFORMULAE */  
   
 #ifdef ALLOW_RUNOFF  
               sflux(i,j,bi,bj) = sflux(i,j,bi,bj) - runoff(i,j,bi,bj)  
 #endif  
   
               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)  
   
             enddo  
           enddo  
         enddo  
       enddo  
   
 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)  
264    
 #ifdef SHORTWAVE_HEATING  
       _EXCH_XY_R8(swflux, mythid)  
 #endif  
265    
266        end        end
Legend:
Removed from v.1.8  
changed lines
  Added in v.1.18
  ViewVC Help
Powered by ViewVC 1.1.22