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