/[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.29 - (hide annotations) (download)
Sat Jun 29 19:44:18 2013 UTC (11 years, 3 months ago) by dimitri
Branch: MAIN
CVS Tags: checkpoint64k, checkpoint64j, checkpoint64m, checkpoint64l, checkpoint64o, checkpoint64n
Changes since 1.28: +20 -20 lines
fix sign and make compatible with temp_EvPrRn .NE. UNSET_RL

1 dimitri 1.29 C $Header: /u/gcmpack/MITgcm/pkg/exf/exf_mapfields.F,v 1.28 2013/04/23 19:04:33 dimitri 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 dimitri 1.28 C 4) By default, RunOff comes at the temp of surface water (with same Cp);
146     C ifdef ALLOW_RUNOFTEMP, RunOff temp can be specified in runoftempfile.
147 jmc 1.27 C 5) Evap is released to the Atmos @ surf-temp (=SST); should be using
148     C the water-vapor heat capacity here and consistently in Bulk-Formulae;
149     C Could also be put directly into Latent Heat flux.
150     IF ( snowPrecipFile .NE. ' ' ) THEN
151     C-- Melt snow (if provided) into the ocean and account for rain-temp
152     DO j = 1, sNy
153     DO i = 1, sNx
154     Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)
155     & + flami*snowPrecip(i,j,bi,bj)*rhoConstFresh
156     & - HeatCapacity_Cp
157     & *( atemp(i,j,bi,bj) - cen2kel - temp_EvPrRn )
158     & *( precip(i,j,bi,bj)- snowPrecip(i,j,bi,bj) )
159     & *rhoConstFresh
160     ENDDO
161     ENDDO
162     ELSE
163     C-- Make snow (according to Air Temp) and melt it in the ocean
164     C note: here we just use the same criteria as over seaice but would be
165     C better to consider a higher altitude air temp, e.g., 850.mb
166     DO j = 1, sNy
167     DO i = 1, sNx
168     IF ( atemp(i,j,bi,bj).LT.cen2kel ) THEN
169     Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)
170     & + flami*precip(i,j,bi,bj)*rhoConstFresh
171     ELSE
172     C-- Account for rain-temp
173     Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)
174     & - HeatCapacity_Cp
175     & *( atemp(i,j,bi,bj) - cen2kel - temp_EvPrRn )
176     & *precip(i,j,bi,bj)*rhoConstFresh
177     ENDIF
178     ENDDO
179     ENDDO
180     ENDIF
181 dimitri 1.29 #ifdef ALLOW_RUNOFF
182     C-- Account for energy content of RunOff:
183     DO j = 1, sNy
184     DO i = 1, sNx
185     Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)
186     & - HeatCapacity_Cp
187     & *( theta(i,j,ks,bi,bj) - temp_EvPrRn )
188     & *runoff(i,j,bi,bj)*rhoConstFresh
189     ENDDO
190     ENDDO
191     #endif
192 dimitri 1.28 C-- Account for energy content of Evap:
193 jmc 1.27 DO j = 1, sNy
194     DO i = 1, sNx
195     Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)
196 dimitri 1.28 & + HeatCapacity_Cp
197 jmc 1.27 & *( theta(i,j,ks,bi,bj) - temp_EvPrRn )
198 dimitri 1.28 & *evap(i,j,bi,bj)*rhoConstFresh
199     Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)*maskC(i,j,ks,bi,bj)
200 jmc 1.27 ENDDO
201     ENDDO
202 dimitri 1.28 ENDIF
203 dimitri 1.29 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
204     #endif /* ALLOW_ATM_TEMP */
205     #if defined(ALLOW_RUNOFF) && defined(ALLOW_RUNOFTEMP)
206     IF ( runoftempfile .NE. ' ' ) THEN
207     C-- Add energy content of RunOff
208 jmc 1.27 DO j = 1, sNy
209     DO i = 1, sNx
210 dimitri 1.28 Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)
211 dimitri 1.29 & + HeatCapacity_Cp
212     & *( theta(i,j,ks,bi,bj) - runoftemp(i,j,bi,bj) )
213 dimitri 1.28 & *runoff(i,j,bi,bj)*rhoConstFresh
214 jmc 1.27 ENDDO
215     ENDDO
216     ENDIF
217 dimitri 1.29 #endif
218 jmc 1.27
219 heimbach 1.2 #ifdef ALLOW_AUTODIFF_TAMC
220     CADJ STORE ustress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
221     #endif
222 jmc 1.26 DO j = jmin,jmax
223     DO i = imin,imax
224     C Zonal wind stress.
225     IF (ustress(i,j,bi,bj).GT.windstressmax) THEN
226 mlosch 1.10 ustress(i,j,bi,bj)=windstressmax
227 jmc 1.26 ENDIF
228     ENDDO
229     ENDDO
230 heimbach 1.2 #ifdef ALLOW_AUTODIFF_TAMC
231     CADJ STORE ustress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
232     #endif
233 jmc 1.26 DO j = jmin,jmax
234     DO i = imin,imax
235     IF (ustress(i,j,bi,bj).LT.-windstressmax) THEN
236 mlosch 1.10 ustress(i,j,bi,bj)=-windstressmax
237 jmc 1.26 ENDIF
238     ENDDO
239     ENDDO
240 jmc 1.20 IF ( stressIsOnCgrid ) THEN
241 jmc 1.26 DO j = jmin,jmax
242     DO i = imin+1,imax
243 jmc 1.20 fu(i,j,bi,bj) = exf_outscal_ustress*ustress(i,j,bi,bj)
244 jmc 1.26 ENDDO
245     ENDDO
246 jmc 1.20 ELSE
247 jmc 1.26 DO j = jmin,jmax
248     DO i = imin+1,imax
249     C Shift wind stresses calculated at Grid-center to W/S points
250 heimbach 1.8 fu(i,j,bi,bj) = exf_outscal_ustress*
251 jmc 1.20 & (ustress(i,j,bi,bj)+ustress(i-1,j,bi,bj))
252 jmc 1.24 & *exf_half*maskW(i,j,ks,bi,bj)
253 jmc 1.26 ENDDO
254     ENDDO
255 jmc 1.20 ENDIF
256 heimbach 1.1
257 heimbach 1.2 #ifdef ALLOW_AUTODIFF_TAMC
258     CADJ STORE vstress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
259     #endif
260 jmc 1.26 DO j = jmin,jmax
261     DO i = imin,imax
262     C Meridional wind stress.
263     IF (vstress(i,j,bi,bj).GT.windstressmax) THEN
264 mlosch 1.10 vstress(i,j,bi,bj)=windstressmax
265 jmc 1.26 ENDIF
266     ENDDO
267     ENDDO
268 heimbach 1.2 #ifdef ALLOW_AUTODIFF_TAMC
269     CADJ STORE vstress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
270     #endif
271 jmc 1.26 DO j = jmin,jmax
272     DO i = imin,imax
273     IF (vstress(i,j,bi,bj).LT.-windstressmax) THEN
274 mlosch 1.10 vstress(i,j,bi,bj)=-windstressmax
275 jmc 1.26 ENDIF
276     ENDDO
277     ENDDO
278 jmc 1.20 IF ( stressIsOnCgrid ) THEN
279 jmc 1.26 DO j = jmin+1,jmax
280     DO i = imin,imax
281 jmc 1.20 fv(i,j,bi,bj) = exf_outscal_vstress*vstress(i,j,bi,bj)
282 jmc 1.26 ENDDO
283     ENDDO
284 jmc 1.20 ELSE
285 jmc 1.26 DO j = jmin+1,jmax
286     DO i = imin,imax
287     C Shift wind stresses calculated at C-points to W/S points
288 heimbach 1.8 fv(i,j,bi,bj) = exf_outscal_vstress*
289 jmc 1.20 & (vstress(i,j,bi,bj)+vstress(i,j-1,bi,bj))
290 jmc 1.24 & *exf_half*maskS(i,j,ks,bi,bj)
291 jmc 1.26 ENDDO
292     ENDDO
293 jmc 1.20 ENDIF
294 heimbach 1.1
295 jmc 1.26 #if defined(ALLOW_ATM_TEMP) || defined(SHORTWAVE_HEATING)
296     C Short wave radiative flux.
297     DO j = jmin,jmax
298     DO i = imin,imax
299     Qsw(i,j,bi,bj) = exf_outscal_swflux*swflux(i,j,bi,bj)
300     ENDDO
301     ENDDO
302 heimbach 1.1 #endif
303    
304     #ifdef ALLOW_CLIMSST_RELAXATION
305 jmc 1.26 DO j = jmin,jmax
306     DO i = imin,imax
307     SST(i,j,bi,bj) = exf_outscal_sst*climsst(i,j,bi,bj)
308     ENDDO
309     ENDDO
310 heimbach 1.1 #endif
311    
312     #ifdef ALLOW_CLIMSSS_RELAXATION
313 jmc 1.26 DO j = jmin,jmax
314     DO i = imin,imax
315     SSS(i,j,bi,bj) = exf_outscal_sss*climsss(i,j,bi,bj)
316     ENDDO
317     ENDDO
318 heimbach 1.1 #endif
319    
320 heimbach 1.2 #ifdef ATMOSPHERIC_LOADING
321 jmc 1.26 DO j = jmin,jmax
322     DO i = imin,imax
323     pLoad(i,j,bi,bj)=exf_outscal_apressure*apressure(i,j,bi,bj)
324     ENDDO
325     ENDDO
326 heimbach 1.2 #endif
327    
328 heimbach 1.25 #ifdef EXF_ALLOW_SEAICE_RELAX
329 jmc 1.26 DO j = jmin,jmax
330     DO i = imin,imax
331     obsSIce(i,j,bi,bj) =
332 heimbach 1.25 & exf_outscal_areamask*areamask(i,j,bi,bj)
333     obsSIce(I,J,bi,bj) =
334     & MIN(MAX(obsSIce(I,J,bi,bj), 0.d0 ), 1.d0)
335 jmc 1.26 ENDDO
336     ENDDO
337 heimbach 1.25 #endif
338    
339 jmc 1.26 ENDDO
340 jmc 1.20 ENDDO
341 heimbach 1.1
342 jmc 1.26 C-- Update the tile edges.
343     _EXCH_XY_RS( Qnet, myThid )
344     _EXCH_XY_RS( EmPmR, myThid )
345 cheisey 1.3 CALL EXCH_UV_XY_RS(fu, fv, .TRUE., myThid)
346 jmc 1.26 c#if defined(ALLOW_ATM_TEMP) || defined(SHORTWAVE_HEATING)
347 dimitri 1.5 #ifdef SHORTWAVE_HEATING
348 jmc 1.26 C Qsw used in SHORTWAVE_HEATING code & for diagnostics (<- EXCH not needed)
349     _EXCH_XY_RS( Qsw, myThid )
350 heimbach 1.1 #endif
351     #ifdef ALLOW_CLIMSST_RELAXATION
352 jmc 1.26 _EXCH_XY_RS( SST, myThid )
353 heimbach 1.1 #endif
354     #ifdef ALLOW_CLIMSSS_RELAXATION
355 jmc 1.26 _EXCH_XY_RS( SSS, myThid )
356 heimbach 1.2 #endif
357     #ifdef ATMOSPHERIC_LOADING
358 jmc 1.26 _EXCH_XY_RS( pLoad, myThid )
359 heimbach 1.1 #endif
360 heimbach 1.25 #ifdef EXF_ALLOW_SEAICE_RELAX
361 jmc 1.26 _EXCH_XY_RS( obsSIce, myThid )
362 heimbach 1.25 #endif
363 heimbach 1.1
364 jmc 1.20 RETURN
365     END

  ViewVC Help
Powered by ViewVC 1.1.22