/[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.33 - (hide annotations) (download)
Fri Oct 6 00:03:56 2017 UTC (6 years, 7 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, HEAD
Changes since 1.32: +12 -1 lines
- add specific forcing field for tides to feed model new geopotential anomaly
  forcing, for now within #ifdef EXF_ALLOW_TIDES.

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

  ViewVC Help
Powered by ViewVC 1.1.22