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

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

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

revision 1.20 by jmc, Mon May 14 19:34:57 2007 UTC revision 1.30 by jmc, Sat Oct 5 19:36:12 2013 UTC
# Line 3  C $Name$ Line 3  C $Name$
3    
4  #include "EXF_OPTIONS.h"  #include "EXF_OPTIONS.h"
5    
6        subroutine exf_mapfields( mytime, myiter, mythid )  CBOP 0
7    C     !ROUTINE: EXF_MAPFIELDS
8    
9  c     ==================================================================  C     !INTERFACE:
10  c     SUBROUTINE exf_mapfields        SUBROUTINE EXF_MAPFIELDS( myTime, myIter, myThid )
 c     ==================================================================  
 c  
 c     o Map external forcing fields (ustress, vstress, hflux, sflux,  
 c       swflux, apressure, climsss, climsst, etc.) onto ocean model  
 c       arrays (fu, fv, Qnet, EmPmR, Qsw, pload, sss, sst, etc.).  
 c       This routine is included to separate the ocean state estimation  
 c       tool as much as possible from the ocean model.  Unit and sign  
 c       conventions can be customized using variables exf_outscal_*,  
 c       which are set in exf_readparms.F.  See the header files  
 c       EXF_FIELDS.h and FFIELDS.h for definitions of the various input  
 c       and output fields and for default unit and sign convetions.  
 c  
 c     started: Christian Eckert eckert@mit.edu  09-Aug-1999  
 c  
 c     changed: Christian Eckert eckert@mit.edu  11-Jan-2000  
 c              - Restructured the code in order to create a package  
 c                for 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  06-May-2000  
 c              - added and changed CPP flag structure for  
 c                ALLOW_BULKFORMULAE, ALLOW_ATM_TEMP  
 c  
 c              Patrick Heimbach, heimbach@mit.edu  23-May-2000  
 c              - sign change of ustress/vstress incorporated into  
 c                scaling factors exf_outscal_ust, exf_outscal_vst  
 c  
 c     mods for pkg/seaice: menemenlis@jpl.nasa.gov 20-Dec-2002  
 c  
 c     ==================================================================  
 c     SUBROUTINE exf_mapfields  
 c     ==================================================================  
11    
12        implicit none  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    
50  c     == global variables ==  C     !USES:
51          IMPLICIT NONE
52    
53    C     == global variables ==
54  #include "EEPARAMS.h"  #include "EEPARAMS.h"
55  #include "SIZE.h"  #include "SIZE.h"
56  #include "PARAMS.h"  #include "PARAMS.h"
57  #include "FFIELDS.h"  #include "FFIELDS.h"
58  #include "GRID.h"  #include "GRID.h"
59    #include "DYNVARS.h"
60    
61  #include "EXF_PARAM.h"  #include "EXF_PARAM.h"
62  #include "EXF_CONSTANTS.h"  #include "EXF_CONSTANTS.h"
# Line 59  c     == global variables == Line 65  c     == global variables ==
65  # include "tamc.h"  # include "tamc.h"
66  # include "tamc_keys.h"  # include "tamc_keys.h"
67  #endif  #endif
 c     == routine arguments ==  
68    
69  c     mythid - thread number for this instance of the routine.  C     !INPUT/OUTPUT PARAMETERS:
70    C     myTime  :: Current time in simulation
71        _RL     mytime  C     myIter  :: Current iteration number
72        integer myiter  C     myThid  :: my Thread Id number
73        integer mythid        _RL     myTime
74          INTEGER myIter
75  c     == local variables ==        INTEGER myThid
76    
77        integer bi,bj  C     !LOCAL VARIABLES:
78        integer i,j,k        INTEGER bi,bj
79          INTEGER i,j,ks
80        INTEGER imin, imax        INTEGER imin, imax
81        INTEGER jmin, jmax        INTEGER jmin, jmax
82        PARAMETER ( imin = 1-OLx , imax = sNx+OLx )        PARAMETER ( imin = 1-OLx , imax = sNx+OLx )
83        PARAMETER ( jmin = 1-OLy , jmax = sNy+OLy )        PARAMETER ( jmin = 1-OLy , jmax = sNy+OLy )
84    CEOP
85    
86  c     == end of interface ==  C--   set surface level index:
87          ks = 1
88    
89        DO bj = myByLo(myThid),myByHi(myThid)        DO bj = myByLo(myThid),myByHi(myThid)
90          DO bi = myBxLo(myThid),myBxHi(myThid)         DO bi = myBxLo(myThid),myBxHi(myThid)
91    
92  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
93            act1 = bi - myBxLo(myThid)            act1 = bi - myBxLo(myThid)
# Line 94  c     == end of interface == Line 102  c     == end of interface ==
102       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
103  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
104    
105            do j = jmin,jmax  C     Heat flux.
106              do i = imin,imax            DO j = jmin,jmax
107  c            Heat flux.              DO i = imin,imax
108               qnet(i,j,bi,bj) = exf_outscal_hflux*hflux(i,j,bi,bj)               Qnet(i,j,bi,bj) = exf_outscal_hflux*hflux(i,j,bi,bj)
109               if ( hfluxfile .EQ. ' ' )              ENDDO
110       &            qnet(i,j,bi,bj) = qnet(i,j,bi,bj) -            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       &            exf_outscal_hflux * ( hflux_exfremo_intercept +       &            exf_outscal_hflux * ( hflux_exfremo_intercept +
116       &            hflux_exfremo_slope*(mytime-starttime) )       &            hflux_exfremo_slope*(myTime-startTime) )
117              enddo              ENDDO
118            enddo             ENDDO
119              ENDIF
120    
121            do j = jmin,jmax  C     Salt flux.
122              do i = imin,imax            DO j = jmin,jmax
123  c            Salt flux.              DO i = imin,imax
124               empmr(i,j,bi,bj)= exf_outscal_sflux*sflux(i,j,bi,bj)               EmPmR(i,j,bi,bj)= exf_outscal_sflux*sflux(i,j,bi,bj)
125               if ( sfluxfile .EQ. ' ' )       &                                          *rhoConstFresh
126       &            empmr(i,j,bi,bj) = empmr(i,j,bi,bj) -              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       &            exf_outscal_sflux * ( sflux_exfremo_intercept +       &            exf_outscal_sflux * ( sflux_exfremo_intercept +
133       &            sflux_exfremo_slope*(mytime-starttime) )       &            sflux_exfremo_slope*(myTime-startTime) )
134              enddo              ENDDO
135            enddo             ENDDO
136              ENDIF
137    
138    #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) 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    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    #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    C--   Account for energy content of Evap:
193               DO j = 1, sNy
194                DO i = 1, sNx
195                  Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)
196         &              + HeatCapacity_Cp
197         &               *( theta(i,j,ks,bi,bj) - temp_EvPrRn )
198         &               *evap(i,j,bi,bj)*rhoConstFresh
199                  Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)*maskC(i,j,ks,bi,bj)
200                ENDDO
201               ENDDO
202              ENDIF
203    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               DO j = 1, sNy
209                DO i = 1, sNx
210                   Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)
211         &              + HeatCapacity_Cp
212         &              *( theta(i,j,ks,bi,bj) - runoftemp(i,j,bi,bj) )
213         &              *runoff(i,j,bi,bj)*rhoConstFresh
214                ENDDO
215               ENDDO
216              ENDIF
217    #endif
218    
219  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
220  CADJ STORE ustress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE ustress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
221  #endif  #endif
222            do j = jmin,jmax            DO j = jmin,jmax
223              do i = imin,imax              DO i = imin,imax
224  c             Zonal wind stress.  C             Zonal wind stress.
225                if (ustress(i,j,bi,bj).gt.windstressmax) then                IF (ustress(i,j,bi,bj).GT.windstressmax) THEN
226                  ustress(i,j,bi,bj)=windstressmax                  ustress(i,j,bi,bj)=windstressmax
227                endif                ENDIF
228              enddo              ENDDO
229            enddo            ENDDO
230  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
231  CADJ STORE ustress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE ustress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
232  #endif  #endif
233            do j = jmin,jmax            DO j = jmin,jmax
234              do i = imin,imax              DO i = imin,imax
235                if (ustress(i,j,bi,bj).lt.-windstressmax) then                IF (ustress(i,j,bi,bj).LT.-windstressmax) THEN
236                  ustress(i,j,bi,bj)=-windstressmax                  ustress(i,j,bi,bj)=-windstressmax
237                endif                ENDIF
238              enddo              ENDDO
239            enddo            ENDDO
240            IF ( stressIsOnCgrid ) THEN            IF ( stressIsOnCgrid ) THEN
241             do j = jmin,jmax             DO j = jmin,jmax
242              do i = imin+1,imax              DO i = imin+1,imax
243                fu(i,j,bi,bj) = exf_outscal_ustress*ustress(i,j,bi,bj)                fu(i,j,bi,bj) = exf_outscal_ustress*ustress(i,j,bi,bj)
244              enddo              ENDDO
245             enddo             ENDDO
246            ELSE            ELSE
247             do j = jmin,jmax             DO j = jmin,jmax
248              do i = imin+1,imax              DO i = imin+1,imax
249  c     Shift wind stresses calculated at Grid-center to W/S points  C     Shift wind stresses calculated at Grid-center to W/S points
250                fu(i,j,bi,bj) = exf_outscal_ustress*                fu(i,j,bi,bj) = exf_outscal_ustress*
251       &              (ustress(i,j,bi,bj)+ustress(i-1,j,bi,bj))       &              (ustress(i,j,bi,bj)+ustress(i-1,j,bi,bj))
252       &              *exf_half*maskW(i,j,1,bi,bj)       &              *exf_half*maskW(i,j,ks,bi,bj)
253              enddo              ENDDO
254             enddo             ENDDO
255            ENDIF            ENDIF
256    
257  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
258  CADJ STORE vstress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE vstress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
259  #endif  #endif
260            do j = jmin,jmax            DO j = jmin,jmax
261              do i = imin,imax              DO i = imin,imax
262  c             Meridional wind stress.  C             Meridional wind stress.
263                if (vstress(i,j,bi,bj).gt.windstressmax) then                IF (vstress(i,j,bi,bj).GT.windstressmax) THEN
264                  vstress(i,j,bi,bj)=windstressmax                  vstress(i,j,bi,bj)=windstressmax
265                endif                ENDIF
266              enddo              ENDDO
267            enddo            ENDDO
268  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
269  CADJ STORE vstress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE vstress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
270  #endif  #endif
271            do j = jmin,jmax            DO j = jmin,jmax
272              do i = imin,imax              DO i = imin,imax
273                if (vstress(i,j,bi,bj).lt.-windstressmax) then                IF (vstress(i,j,bi,bj).LT.-windstressmax) THEN
274                  vstress(i,j,bi,bj)=-windstressmax                  vstress(i,j,bi,bj)=-windstressmax
275                endif                ENDIF
276              enddo              ENDDO
277            enddo            ENDDO
278            IF ( stressIsOnCgrid ) THEN            IF ( stressIsOnCgrid ) THEN
279             do j = jmin+1,jmax             DO j = jmin+1,jmax
280              do i = imin,imax              DO i = imin,imax
281                fv(i,j,bi,bj) = exf_outscal_vstress*vstress(i,j,bi,bj)                fv(i,j,bi,bj) = exf_outscal_vstress*vstress(i,j,bi,bj)
282              enddo              ENDDO
283             enddo             ENDDO
284            ELSE            ELSE
285             do j = jmin+1,jmax             DO j = jmin+1,jmax
286              do i = imin,imax              DO i = imin,imax
287  c     Shift wind stresses calculated at C-points to W/S points  C     Shift wind stresses calculated at C-points to W/S points
288                fv(i,j,bi,bj) = exf_outscal_vstress*                fv(i,j,bi,bj) = exf_outscal_vstress*
289       &              (vstress(i,j,bi,bj)+vstress(i,j-1,bi,bj))       &              (vstress(i,j,bi,bj)+vstress(i,j-1,bi,bj))
290       &              *exf_half*maskS(i,j,1,bi,bj)       &              *exf_half*maskS(i,j,ks,bi,bj)
291              enddo              ENDDO
292             enddo             ENDDO
293            ENDIF            ENDIF
294    
295  #ifdef SHORTWAVE_HEATING  #if defined(ALLOW_ATM_TEMP) || defined(SHORTWAVE_HEATING)
296  c             Short wave radiative flux.  C             Short wave radiative flux.
297            do j = jmin,jmax            DO j = jmin,jmax
298              do i = imin,imax              DO i = imin,imax
299               qsw(i,j,bi,bj)  = exf_outscal_swflux*swflux(i,j,bi,bj)               Qsw(i,j,bi,bj)  = exf_outscal_swflux*swflux(i,j,bi,bj)
300              enddo              ENDDO
301            enddo            ENDDO
302  #endif  #endif
303    
304  #ifdef ALLOW_CLIMSST_RELAXATION  #ifdef ALLOW_CLIMSST_RELAXATION
305            do j = jmin,jmax            DO j = jmin,jmax
306              do i = imin,imax              DO i = imin,imax
307               sst(i,j,bi,bj)  = exf_outscal_sst*climsst(i,j,bi,bj)               SST(i,j,bi,bj)  = exf_outscal_sst*climsst(i,j,bi,bj)
308              enddo              ENDDO
309            enddo            ENDDO
310  #endif  #endif
311    
312  #ifdef ALLOW_CLIMSSS_RELAXATION  #ifdef ALLOW_CLIMSSS_RELAXATION
313            do j = jmin,jmax            DO j = jmin,jmax
314              do i = imin,imax              DO i = imin,imax
315               sss(i,j,bi,bj)  = exf_outscal_sss*climsss(i,j,bi,bj)               SSS(i,j,bi,bj)  = exf_outscal_sss*climsss(i,j,bi,bj)
316              enddo              ENDDO
317            enddo            ENDDO
318  #endif  #endif
319    
320  #ifdef ATMOSPHERIC_LOADING  #ifdef ATMOSPHERIC_LOADING
321            do j = jmin,jmax            DO j = jmin,jmax
322              do i = imin,imax              DO i = imin,imax
323               pload(i,j,bi,bj)=exf_outscal_apressure*apressure(i,j,bi,bj)               pLoad(i,j,bi,bj)=exf_outscal_apressure*apressure(i,j,bi,bj)
324              enddo              ENDDO
325            enddo            ENDDO
326    #endif
327    
328    #ifdef EXF_SEAICE_FRACTION
329              DO j = jmin,jmax
330                DO i = imin,imax
331                  exf_iceFraction(i,j,bi,bj) =
332         &           exf_outscal_areamask*areamask(i,j,bi,bj)
333                  exf_iceFraction(I,J,bi,bj) =
334         &           MIN( MAX(exf_iceFraction(I,J,bi,bj),zeroRS), oneRS )
335                ENDDO
336              ENDDO
337  #endif  #endif
338    
339          ENDDO         ENDDO
340        ENDDO        ENDDO
341    
342  c     Update the tile edges.  C--   Update the tile edges.
343          _EXCH_XY_RS(  Qnet, myThid )
344        _EXCH_XY_R4(  qnet, mythid )        _EXCH_XY_RS( EmPmR, myThid )
       _EXCH_XY_R4( empmr, mythid )  
345         CALL EXCH_UV_XY_RS(fu, fv, .TRUE., myThid)         CALL EXCH_UV_XY_RS(fu, fv, .TRUE., myThid)
346    c#if defined(ALLOW_ATM_TEMP) || defined(SHORTWAVE_HEATING)
347  #ifdef SHORTWAVE_HEATING  #ifdef SHORTWAVE_HEATING
348        _EXCH_XY_R4(   qsw, mythid )  C     Qsw used in SHORTWAVE_HEATING code & for diagnostics (<- EXCH not needed)
349          _EXCH_XY_RS(   Qsw, myThid )
350  #endif  #endif
351  #ifdef ALLOW_CLIMSST_RELAXATION  #ifdef ALLOW_CLIMSST_RELAXATION
352        _EXCH_XY_R4(   sst, mythid )        _EXCH_XY_RS(   SST, myThid )
353  #endif  #endif
354  #ifdef ALLOW_CLIMSSS_RELAXATION  #ifdef ALLOW_CLIMSSS_RELAXATION
355        _EXCH_XY_R4(   sss, mythid )        _EXCH_XY_RS(   SSS, myThid )
356  #endif  #endif
357  #ifdef ATMOSPHERIC_LOADING  #ifdef ATMOSPHERIC_LOADING
358        _EXCH_XY_R4( pload, mythid )        _EXCH_XY_RS( pLoad, myThid )
359    #endif
360    #ifdef EXF_SEAICE_FRACTION
361          _EXCH_XY_RS( exf_iceFraction, myThid )
362  #endif  #endif
363    
364        RETURN        RETURN
Legend:
Removed from v.1.20  
changed lines
  Added in v.1.30
  ViewVC Help
Powered by ViewVC 1.1.22