/[MITgcm]/MITgcm/pkg/exf/exf_mapfields.F
ViewVC logotype

Annotation of /MITgcm/pkg/exf/exf_mapfields.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.27 - (hide annotations) (download)
Thu Apr 4 16:30:50 2013 UTC (11 years, 5 months ago) by jmc
Branch: MAIN
Changes since 1.26: +69 -1 lines
In case Energy-Reference-Level is used (temp_EvPrRn=0), account
for energy content of Precip + RunOff & Evap. Assumes:
 1) Rain has same temp as Air
 2) Snow has no heat capacity (consistent with seaice & thsice pkgs)
 3) No distinction between sea-water Cp and fresh-water Cp
 4) Run-Off comes at the temp of surface water (with same Cp)
 5) Evap is released to the Atmos @ surf-temp (=SST); should be using
    the water-vapor heat capacity here and consistently in Bulk-Formulae;
    and could also be put directly into Latent Heat flux.

1 jmc 1.27 C $Header: /u/gcmpack/MITgcm/pkg/exf/exf_mapfields.F,v 1.26 2013/01/10 22:59:24 jmc Exp $
2 jmc 1.17 C $Name: $
3 heimbach 1.1
4 edhill 1.7 #include "EXF_OPTIONS.h"
5 heimbach 1.1
6 jmc 1.26 CBOP 0
7     C !ROUTINE: EXF_MAPFIELDS
8 heimbach 1.1
9 jmc 1.26 C !INTERFACE:
10     SUBROUTINE EXF_MAPFIELDS( myTime, myIter, myThid )
11 heimbach 1.1
12 jmc 1.26 C !DESCRIPTION:
13     C ==================================================================
14     C SUBROUTINE EXF_MAPFIELDS
15     C ==================================================================
16     C
17     C o Map external forcing fields (ustress, vstress, hflux, sflux,
18     C swflux, apressure, climsss, climsst, etc.) onto ocean model
19     C arrays (fu, fv, Qnet, EmPmR, Qsw, pLoad, SSS, SST, etc.).
20     C This routine is included to separate the ocean state estimation
21     C tool as much as possible from the ocean model. Unit and sign
22     C conventions can be customized using variables exf_outscal_*,
23     C which are set in exf_readparms.F. See the header files
24     C EXF_FIELDS.h and FFIELDS.h for definitions of the various input
25     C and output fields and for default unit and sign convetions.
26     C
27     C started: Christian Eckert eckert@mit.edu 09-Aug-1999
28     C
29     C changed: Christian Eckert eckert@mit.edu 11-Jan-2000
30     C - Restructured the code in order to create a package
31     C for the MITgcmUV.
32     C
33     C Christian Eckert eckert@mit.edu 12-Feb-2000
34     C - Changed Routine names (package prefix: exf_)
35     C
36     C Patrick Heimbach, heimbach@mit.edu 06-May-2000
37     C - added and changed CPP flag structure for
38     C ALLOW_BULKFORMULAE, ALLOW_ATM_TEMP
39     C
40     C Patrick Heimbach, heimbach@mit.edu 23-May-2000
41     C - sign change of ustress/vstress incorporated into
42     C scaling factors exf_outscal_ust, exf_outscal_vst
43     C
44     C mods for pkg/seaice: menemenlis@jpl.nasa.gov 20-Dec-2002
45     C
46     C ==================================================================
47     C SUBROUTINE EXF_MAPFIELDS
48     C ==================================================================
49 heimbach 1.1
50 jmc 1.26 C !USES:
51     IMPLICIT NONE
52 heimbach 1.1
53 jmc 1.26 C == global variables ==
54 heimbach 1.1 #include "EEPARAMS.h"
55     #include "SIZE.h"
56 heimbach 1.16 #include "PARAMS.h"
57 heimbach 1.1 #include "FFIELDS.h"
58 heimbach 1.8 #include "GRID.h"
59 jmc 1.27 #include "DYNVARS.h"
60 heimbach 1.8
61 jmc 1.17 #include "EXF_PARAM.h"
62     #include "EXF_CONSTANTS.h"
63     #include "EXF_FIELDS.h"
64 heimbach 1.2 #ifdef ALLOW_AUTODIFF_TAMC
65     # include "tamc.h"
66     # include "tamc_keys.h"
67     #endif
68 heimbach 1.1
69 jmc 1.26 C !INPUT/OUTPUT PARAMETERS:
70     C myTime :: Current time in simulation
71     C myIter :: Current iteration number
72     C myThid :: my Thread Id number
73     _RL myTime
74     INTEGER myIter
75     INTEGER myThid
76 heimbach 1.1
77 jmc 1.26 C !LOCAL VARIABLES:
78 jmc 1.24 INTEGER bi,bj
79     INTEGER i,j,ks
80 jmc 1.20 INTEGER imin, imax
81     INTEGER jmin, jmax
82     PARAMETER ( imin = 1-OLx , imax = sNx+OLx )
83     PARAMETER ( jmin = 1-OLy , jmax = sNy+OLy )
84 jmc 1.26 CEOP
85 heimbach 1.1
86 jmc 1.24 C-- set surface level index:
87     ks = 1
88    
89 jmc 1.20 DO bj = myByLo(myThid),myByHi(myThid)
90 jmc 1.26 DO bi = myBxLo(myThid),myBxHi(myThid)
91 heimbach 1.2
92     #ifdef ALLOW_AUTODIFF_TAMC
93     act1 = bi - myBxLo(myThid)
94     max1 = myBxHi(myThid) - myBxLo(myThid) + 1
95     act2 = bj - myByLo(myThid)
96     max2 = myByHi(myThid) - myByLo(myThid) + 1
97     act3 = myThid - 1
98     max3 = nTx*nTy
99     act4 = ikey_dynamics - 1
100     ikey = (act1 + 1) + act2*max1
101     & + act3*max1*max2
102     & + act4*max1*max2*max3
103     #endif /* ALLOW_AUTODIFF_TAMC */
104    
105 jmc 1.26 C Heat flux.
106     DO j = jmin,jmax
107     DO i = imin,imax
108     Qnet(i,j,bi,bj) = exf_outscal_hflux*hflux(i,j,bi,bj)
109     ENDDO
110     ENDDO
111     IF ( hfluxfile .EQ. ' ' ) THEN
112     DO j = jmin,jmax
113     DO i = imin,imax
114     Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) -
115 heimbach 1.16 & exf_outscal_hflux * ( hflux_exfremo_intercept +
116 jmc 1.26 & hflux_exfremo_slope*(myTime-startTime) )
117     ENDDO
118     ENDDO
119     ENDIF
120    
121     C Salt flux.
122     DO j = jmin,jmax
123     DO i = imin,imax
124 jmc 1.22 EmPmR(i,j,bi,bj)= exf_outscal_sflux*sflux(i,j,bi,bj)
125     & *rhoConstFresh
126 jmc 1.26 ENDDO
127     ENDDO
128     IF ( sfluxfile .EQ. ' ' ) THEN
129     DO j = jmin,jmax
130     DO i = imin,imax
131     EmPmR(i,j,bi,bj) = EmPmR(i,j,bi,bj) - rhoConstFresh*
132 heimbach 1.16 & exf_outscal_sflux * ( sflux_exfremo_intercept +
133 jmc 1.26 & sflux_exfremo_slope*(myTime-startTime) )
134     ENDDO
135     ENDDO
136     ENDIF
137 heimbach 1.1
138 jmc 1.27 #ifdef ALLOW_ATM_TEMP
139     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
140     IF ( temp_EvPrRn .NE. UNSET_RL ) THEN
141     C-- Account for energy content of Precip + RunOff & Evap. Assumes:
142     C 1) Rain has same temp as Air
143     C 2) Snow has no heat capacity (consistent with seaice & thsice pkgs)
144     C 3) No distinction between sea-water Cp and fresh-water Cp
145     C 4) Run-Off comes at the temp of surface water (with same Cp)
146     C 5) Evap is released to the Atmos @ surf-temp (=SST); should be using
147     C the water-vapor heat capacity here and consistently in Bulk-Formulae;
148     C Could also be put directly into Latent Heat flux.
149     IF ( snowPrecipFile .NE. ' ' ) THEN
150     C-- Melt snow (if provided) into the ocean and account for rain-temp
151     DO j = 1, sNy
152     DO i = 1, sNx
153     Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)
154     & + flami*snowPrecip(i,j,bi,bj)*rhoConstFresh
155     & - HeatCapacity_Cp
156     & *( atemp(i,j,bi,bj) - cen2kel - temp_EvPrRn )
157     & *( precip(i,j,bi,bj)- snowPrecip(i,j,bi,bj) )
158     & *rhoConstFresh
159     ENDDO
160     ENDDO
161     ELSE
162     C-- Make snow (according to Air Temp) and melt it in the ocean
163     C note: here we just use the same criteria as over seaice but would be
164     C better to consider a higher altitude air temp, e.g., 850.mb
165     DO j = 1, sNy
166     DO i = 1, sNx
167     IF ( atemp(i,j,bi,bj).LT.cen2kel ) THEN
168     Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)
169     & + flami*precip(i,j,bi,bj)*rhoConstFresh
170     ELSE
171     C-- Account for rain-temp
172     Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)
173     & - HeatCapacity_Cp
174     & *( atemp(i,j,bi,bj) - cen2kel - temp_EvPrRn )
175     & *precip(i,j,bi,bj)*rhoConstFresh
176     ENDIF
177     ENDDO
178     ENDDO
179     ENDIF
180     #ifdef ALLOW_RUNOFF
181     C-- Account for energy content of RunOff:
182     DO j = 1, sNy
183     DO i = 1, sNx
184     Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)
185     & - HeatCapacity_Cp
186     & *( theta(i,j,ks,bi,bj) - temp_EvPrRn )
187     & *runoff(i,j,bi,bj)*rhoConstFresh
188     ENDDO
189     ENDDO
190     #endif
191     C-- Account for energy content of Evap:
192     DO j = 1, sNy
193     DO i = 1, sNx
194     Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)
195     & + HeatCapacity_Cp
196     & *( theta(i,j,ks,bi,bj) - temp_EvPrRn )
197     & *evap(i,j,bi,bj)*rhoConstFresh
198     Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)*maskC(i,j,ks,bi,bj)
199     ENDDO
200     ENDDO
201     ENDIF
202     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
203     #endif /* ALLOW_ATM_TEMP */
204    
205 heimbach 1.2 #ifdef ALLOW_AUTODIFF_TAMC
206     CADJ STORE ustress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
207     #endif
208 jmc 1.26 DO j = jmin,jmax
209     DO i = imin,imax
210     C Zonal wind stress.
211     IF (ustress(i,j,bi,bj).GT.windstressmax) THEN
212 mlosch 1.10 ustress(i,j,bi,bj)=windstressmax
213 jmc 1.26 ENDIF
214     ENDDO
215     ENDDO
216 heimbach 1.2 #ifdef ALLOW_AUTODIFF_TAMC
217     CADJ STORE ustress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
218     #endif
219 jmc 1.26 DO j = jmin,jmax
220     DO i = imin,imax
221     IF (ustress(i,j,bi,bj).LT.-windstressmax) THEN
222 mlosch 1.10 ustress(i,j,bi,bj)=-windstressmax
223 jmc 1.26 ENDIF
224     ENDDO
225     ENDDO
226 jmc 1.20 IF ( stressIsOnCgrid ) THEN
227 jmc 1.26 DO j = jmin,jmax
228     DO i = imin+1,imax
229 jmc 1.20 fu(i,j,bi,bj) = exf_outscal_ustress*ustress(i,j,bi,bj)
230 jmc 1.26 ENDDO
231     ENDDO
232 jmc 1.20 ELSE
233 jmc 1.26 DO j = jmin,jmax
234     DO i = imin+1,imax
235     C Shift wind stresses calculated at Grid-center to W/S points
236 heimbach 1.8 fu(i,j,bi,bj) = exf_outscal_ustress*
237 jmc 1.20 & (ustress(i,j,bi,bj)+ustress(i-1,j,bi,bj))
238 jmc 1.24 & *exf_half*maskW(i,j,ks,bi,bj)
239 jmc 1.26 ENDDO
240     ENDDO
241 jmc 1.20 ENDIF
242 heimbach 1.1
243 heimbach 1.2 #ifdef ALLOW_AUTODIFF_TAMC
244     CADJ STORE vstress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
245     #endif
246 jmc 1.26 DO j = jmin,jmax
247     DO i = imin,imax
248     C Meridional wind stress.
249     IF (vstress(i,j,bi,bj).GT.windstressmax) THEN
250 mlosch 1.10 vstress(i,j,bi,bj)=windstressmax
251 jmc 1.26 ENDIF
252     ENDDO
253     ENDDO
254 heimbach 1.2 #ifdef ALLOW_AUTODIFF_TAMC
255     CADJ STORE vstress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
256     #endif
257 jmc 1.26 DO j = jmin,jmax
258     DO i = imin,imax
259     IF (vstress(i,j,bi,bj).LT.-windstressmax) THEN
260 mlosch 1.10 vstress(i,j,bi,bj)=-windstressmax
261 jmc 1.26 ENDIF
262     ENDDO
263     ENDDO
264 jmc 1.20 IF ( stressIsOnCgrid ) THEN
265 jmc 1.26 DO j = jmin+1,jmax
266     DO i = imin,imax
267 jmc 1.20 fv(i,j,bi,bj) = exf_outscal_vstress*vstress(i,j,bi,bj)
268 jmc 1.26 ENDDO
269     ENDDO
270 jmc 1.20 ELSE
271 jmc 1.26 DO j = jmin+1,jmax
272     DO i = imin,imax
273     C Shift wind stresses calculated at C-points to W/S points
274 heimbach 1.8 fv(i,j,bi,bj) = exf_outscal_vstress*
275 jmc 1.20 & (vstress(i,j,bi,bj)+vstress(i,j-1,bi,bj))
276 jmc 1.24 & *exf_half*maskS(i,j,ks,bi,bj)
277 jmc 1.26 ENDDO
278     ENDDO
279 jmc 1.20 ENDIF
280 heimbach 1.1
281 jmc 1.26 #if defined(ALLOW_ATM_TEMP) || defined(SHORTWAVE_HEATING)
282     C Short wave radiative flux.
283     DO j = jmin,jmax
284     DO i = imin,imax
285     Qsw(i,j,bi,bj) = exf_outscal_swflux*swflux(i,j,bi,bj)
286     ENDDO
287     ENDDO
288 heimbach 1.1 #endif
289    
290     #ifdef ALLOW_CLIMSST_RELAXATION
291 jmc 1.26 DO j = jmin,jmax
292     DO i = imin,imax
293     SST(i,j,bi,bj) = exf_outscal_sst*climsst(i,j,bi,bj)
294     ENDDO
295     ENDDO
296 heimbach 1.1 #endif
297    
298     #ifdef ALLOW_CLIMSSS_RELAXATION
299 jmc 1.26 DO j = jmin,jmax
300     DO i = imin,imax
301     SSS(i,j,bi,bj) = exf_outscal_sss*climsss(i,j,bi,bj)
302     ENDDO
303     ENDDO
304 heimbach 1.1 #endif
305    
306 heimbach 1.2 #ifdef ATMOSPHERIC_LOADING
307 jmc 1.26 DO j = jmin,jmax
308     DO i = imin,imax
309     pLoad(i,j,bi,bj)=exf_outscal_apressure*apressure(i,j,bi,bj)
310     ENDDO
311     ENDDO
312 heimbach 1.2 #endif
313    
314 heimbach 1.25 #ifdef EXF_ALLOW_SEAICE_RELAX
315 jmc 1.26 DO j = jmin,jmax
316     DO i = imin,imax
317     obsSIce(i,j,bi,bj) =
318 heimbach 1.25 & exf_outscal_areamask*areamask(i,j,bi,bj)
319     obsSIce(I,J,bi,bj) =
320     & MIN(MAX(obsSIce(I,J,bi,bj), 0.d0 ), 1.d0)
321 jmc 1.26 ENDDO
322     ENDDO
323 heimbach 1.25 #endif
324    
325 jmc 1.26 ENDDO
326 jmc 1.20 ENDDO
327 heimbach 1.1
328 jmc 1.26 C-- Update the tile edges.
329     _EXCH_XY_RS( Qnet, myThid )
330     _EXCH_XY_RS( EmPmR, myThid )
331 cheisey 1.3 CALL EXCH_UV_XY_RS(fu, fv, .TRUE., myThid)
332 jmc 1.26 c#if defined(ALLOW_ATM_TEMP) || defined(SHORTWAVE_HEATING)
333 dimitri 1.5 #ifdef SHORTWAVE_HEATING
334 jmc 1.26 C Qsw used in SHORTWAVE_HEATING code & for diagnostics (<- EXCH not needed)
335     _EXCH_XY_RS( Qsw, myThid )
336 heimbach 1.1 #endif
337     #ifdef ALLOW_CLIMSST_RELAXATION
338 jmc 1.26 _EXCH_XY_RS( SST, myThid )
339 heimbach 1.1 #endif
340     #ifdef ALLOW_CLIMSSS_RELAXATION
341 jmc 1.26 _EXCH_XY_RS( SSS, myThid )
342 heimbach 1.2 #endif
343     #ifdef ATMOSPHERIC_LOADING
344 jmc 1.26 _EXCH_XY_RS( pLoad, myThid )
345 heimbach 1.1 #endif
346 heimbach 1.25 #ifdef EXF_ALLOW_SEAICE_RELAX
347 jmc 1.26 _EXCH_XY_RS( obsSIce, myThid )
348 heimbach 1.25 #endif
349 heimbach 1.1
350 jmc 1.20 RETURN
351     END

  ViewVC Help
Powered by ViewVC 1.1.22