/[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.2 by heimbach, Tue Nov 12 20:34:41 2002 UTC revision 1.32 by gforget, Wed Jan 11 03:45:34 2017 UTC
# Line 1  Line 1 
1  c $Header$  C $Header$
2    C $Name$
3    
4  #include "EXF_CPPOPTIONS.h"  #include "EXF_OPTIONS.h"
5    #ifdef ALLOW_AUTODIFF
6    # include "AUTODIFF_OPTIONS.h"
7    #endif
8    
9    CBOP 0
10    C     !ROUTINE: EXF_MAPFIELDS
11    
12    C     !INTERFACE:
13          SUBROUTINE EXF_MAPFIELDS( myTime, myIter, myThid )
14    
15    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    
53    C     !USES:
54          IMPLICIT NONE
55    
56        subroutine exf_mapfields( mythid )  C     == global variables ==
   
 c     ==================================================================  
 c     SUBROUTINE exf_mapfields  
 c     ==================================================================  
 c  
 c     o Map the external forcing fields on the ocean model arrays. This  
 c       routine is included to separate the ocean state estimation tool  
 c       as much as possible from the ocean model. Unit conversion factors  
 c       are to be set by the user.  
 c  
 c       The units have to be such that the individual forcing record has  
 c       units equal to [quantity/s]. See the header file *FFIELDS.h* of  
 c       the MITgcmuv.  
 c  
 c       Required units such that no scaling has to be applied:  
 c  
 c       heat flux:          input file W/m^2  
 c       salt flux:          input file m/s  
 c       zonal wind stress:  input file N/m^2  
 c       merid. wind stress: input file N/m^2  
 c  
 c       To allow for such unit conversions this routine contains scaling  
 c       factors scal_quantity.  
 c  
 c     started: Christian Eckert eckert@mit.edu  09-Aug-1999  
 c  
 c     changed: Christian Eckert eckert@mit.edu  11-Jan-2000  
 c  
 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  
 c              - Changed Routine names (package prefix: exf_)  
 c  
 c              Patrick Heimbach, heimbach@mit.edu  06-May-2000  
 c  
 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  
 c              - sign change of ustress/vstress incorporated into  
 c                scaling factors scal_ust, scal_vst  
 c  
 c     ==================================================================  
 c     SUBROUTINE exf_mapfields  
 c     ==================================================================  
   
       implicit none  
   
 c     == global variables ==  
   
57  #include "EEPARAMS.h"  #include "EEPARAMS.h"
58  #include "SIZE.h"  #include "SIZE.h"
59    #include "PARAMS.h"
60  #include "FFIELDS.h"  #include "FFIELDS.h"
61  #include "exf_param.h"  #include "GRID.h"
62  #include "exf_constants.h"  #include "DYNVARS.h"
63  #include "exf_fields.h"  
64  #include "exf_clim_fields.h"  #include "EXF_PARAM.h"
65    #include "EXF_CONSTANTS.h"
66    #include "EXF_FIELDS.h"
67  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
68  # include "tamc.h"  # include "tamc.h"
69  # include "tamc_keys.h"  # include "tamc_keys.h"
70  #endif  #endif
 c     == routine arguments ==  
   
 c     mythid - thread number for this instance of the routine.  
   
       integer mythid  
71    
72  c     == local variables ==  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    
80    C     !LOCAL VARIABLES:
81          INTEGER bi,bj
82          INTEGER i,j,ks
83          INTEGER imin, imax
84          INTEGER jmin, jmax
85          PARAMETER ( imin = 1-OLx , imax = sNx+OLx )
86          PARAMETER ( jmin = 1-OLy , jmax = sNy+OLy )
87    CEOP
88    
89        integer bi,bj  C--   set surface level index:
90        integer i,j        ks = 1
       integer jtlo  
       integer jthi  
       integer itlo  
       integer ithi  
       integer jmin  
       integer jmax  
       integer imin  
       integer imax  
   
 c     == end of interface ==  
   
       jtlo = mybylo(mythid)  
       jthi = mybyhi(mythid)  
       itlo = mybxlo(mythid)  
       ithi = mybxhi(mythid)  
       jmin = 1-oly  
       jmax = sny+oly  
       imin = 1-olx  
       imax = snx+olx  
91    
92        do bj = jtlo,jthi        DO bj = myByLo(myThid),myByHi(myThid)
93          do bi = itlo,ithi         DO bi = myBxLo(myThid),myBxHi(myThid)
94    
95  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
96            act1 = bi - myBxLo(myThid)            act1 = bi - myBxLo(myThid)
# Line 114  c     == end of interface == Line 105  c     == end of interface ==
105       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
106  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
107    
108            do j = jmin,jmax  C     Heat flux.
109              do i = imin,imax            DO j = jmin,jmax
110  c             Heat flux.              DO i = imin,imax
111                qnet(i,j,bi,bj)  = scal_hfl*hflux(i,j,bi,bj)               Qnet(i,j,bi,bj) = exf_outscal_hflux*hflux(i,j,bi,bj)
112              enddo              ENDDO
113            enddo            ENDDO
114              IF ( hfluxfile .EQ. ' ' ) THEN
115               DO j = jmin,jmax
116            do j = jmin,jmax              DO i = imin,imax
117              do i = imin,imax               Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) -
118  c             Salt flux.       &            exf_outscal_hflux * ( hflux_exfremo_intercept +
119  #if (defined (ALLOW_BULKFORMULAE)  && defined (ALLOW_ATM_TEMP))       &            hflux_exfremo_slope*(myTime-startTime) )
120                empmr(i,j,bi,bj) = scal_prc*precip(i,j,bi,bj)              ENDDO
121  #else             ENDDO
122                empmr(i,j,bi,bj) = scal_sfl*sflux(i,j,bi,bj)            ENDIF
123    
124    C     Freshwater flux.
125              DO j = jmin,jmax
126                DO i = imin,imax
127                 EmPmR(i,j,bi,bj)= exf_outscal_sflux*sflux(i,j,bi,bj)
128         &                                          *rhoConstFresh
129                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         &            exf_outscal_sflux * ( sflux_exfremo_intercept +
136         &            sflux_exfremo_slope*(myTime-startTime) )
137                ENDDO
138               ENDDO
139              ENDIF
140    
141    #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    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    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    #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    C--   Account for energy content of Evap:
196               DO j = 1, sNy
197                DO i = 1, sNx
198                  Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)
199         &              + HeatCapacity_Cp
200         &               *( theta(i,j,ks,bi,bj) - temp_EvPrRn )
201         &               *evap(i,j,bi,bj)*rhoConstFresh
202                  Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)*maskC(i,j,ks,bi,bj)
203                ENDDO
204               ENDDO
205              ENDIF
206    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               DO j = 1, sNy
212                DO i = 1, sNx
213                   Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)
214         &              + HeatCapacity_Cp
215         &              *( theta(i,j,ks,bi,bj) - runoftemp(i,j,bi,bj) )
216         &              *runoff(i,j,bi,bj)*rhoConstFresh
217                ENDDO
218               ENDDO
219              ENDIF
220  #endif  #endif
             enddo  
           enddo  
221    
222  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
223  CADJ STORE ustress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE ustress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
224  #endif  #endif
225            do j = jmin,jmax            DO j = jmin,jmax
226              do i = imin,imax              DO i = imin,imax
227  c             Zonal wind stress.  C             Zonal wind stress.
228                if (ustress(i,j,bi,bj).gt.2.0D0) then                IF (ustress(i,j,bi,bj).GT.windstressmax) THEN
229                  ustress(i,j,bi,bj)=2.0D0                  ustress(i,j,bi,bj)=windstressmax
230                endif                ENDIF
231              enddo              ENDDO
232            enddo            ENDDO
233  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
234  CADJ STORE ustress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE ustress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
235  #endif  #endif
236            do j = jmin,jmax            DO j = jmin,jmax
237              do i = imin,imax              DO i = imin,imax
238                if (ustress(i,j,bi,bj).lt.-2.0D0) then                IF (ustress(i,j,bi,bj).LT.-windstressmax) THEN
239                  ustress(i,j,bi,bj)=-2.0D0                  ustress(i,j,bi,bj)=-windstressmax
240                endif                ENDIF
241              enddo              ENDDO
242            enddo            ENDDO
243            do j = jmin,jmax            IF ( stressIsOnCgrid ) THEN
244              do i = imin,imax             DO j = jmin,jmax
245                fu(i,j,bi,bj)    = scal_ust*ustress(i,j,bi,bj)              DO i = imin+1,imax
246              enddo                fu(i,j,bi,bj) = exf_outscal_ustress*ustress(i,j,bi,bj)
247            enddo              ENDDO
248               ENDDO
249              ELSE
250               DO j = jmin,jmax
251                DO i = imin+1,imax
252    C     Shift wind stresses calculated at Grid-center to W/S points
253                  fu(i,j,bi,bj) = exf_outscal_ustress*
254         &              (ustress(i,j,bi,bj)+ustress(i-1,j,bi,bj))
255         &              *exf_half*maskW(i,j,ks,bi,bj)
256                ENDDO
257               ENDDO
258              ENDIF
259    
260  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
261  CADJ STORE vstress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE vstress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
262  #endif  #endif
263            do j = jmin,jmax            DO j = jmin,jmax
264              do i = imin,imax              DO i = imin,imax
265  c             Meridional wind stress.  C             Meridional wind stress.
266                if (vstress(i,j,bi,bj).gt.2.0D0) then                IF (vstress(i,j,bi,bj).GT.windstressmax) THEN
267                  vstress(i,j,bi,bj)=2.0D0                  vstress(i,j,bi,bj)=windstressmax
268                endif                ENDIF
269              enddo              ENDDO
270            enddo            ENDDO
271  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
272  CADJ STORE vstress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE vstress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
273  #endif  #endif
274            do j = jmin,jmax            DO j = jmin,jmax
275              do i = imin,imax              DO i = imin,imax
276                if (vstress(i,j,bi,bj).lt.-2.0D0) then                IF (vstress(i,j,bi,bj).LT.-windstressmax) THEN
277                  vstress(i,j,bi,bj)=-2.0D0                  vstress(i,j,bi,bj)=-windstressmax
278                endif                ENDIF
279              enddo              ENDDO
280            enddo            ENDDO
281            do j = jmin,jmax            IF ( stressIsOnCgrid ) THEN
282              do i = imin,imax             DO j = jmin+1,jmax
283                fv(i,j,bi,bj)    = scal_vst*vstress(i,j,bi,bj)              DO i = imin,imax
284              enddo                fv(i,j,bi,bj) = exf_outscal_vstress*vstress(i,j,bi,bj)
285            enddo              ENDDO
286               ENDDO
287  #ifdef ALLOW_KPP || \            ELSE
288   (defined (ALLOW_BULKFORMULAE) && defined (ALLOW_ATM_TEMP)))             DO j = jmin+1,jmax
289  c             Short wave radiative flux.              DO i = imin,imax
290            do j = jmin,jmax  C     Shift wind stresses calculated at C-points to W/S points
291              do i = imin,imax                fv(i,j,bi,bj) = exf_outscal_vstress*
292                qsw(i,j,bi,bj)   = scal_swf*swflux(i,j,bi,bj)       &              (vstress(i,j,bi,bj)+vstress(i,j-1,bi,bj))
293              enddo       &              *exf_half*maskS(i,j,ks,bi,bj)
294            enddo              ENDDO
295               ENDDO
296              ENDIF
297    
298    #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  #endif  #endif
306    
307  #ifdef ALLOW_CLIMSST_RELAXATION  #ifdef ALLOW_CLIMSST_RELAXATION
308            do j = jmin,jmax            DO j = jmin,jmax
309              do i = imin,imax              DO i = imin,imax
310                sst(i,j,bi,bj)   = scal_sst*climsst(i,j,bi,bj)               SST(i,j,bi,bj)  = exf_outscal_sst*climsst(i,j,bi,bj)
311              enddo              ENDDO
312            enddo            ENDDO
313  #endif  #endif
314    
315  #ifdef ALLOW_CLIMSSS_RELAXATION  #ifdef ALLOW_CLIMSSS_RELAXATION
316            do j = jmin,jmax            DO j = jmin,jmax
317              do i = imin,imax              DO i = imin,imax
318                sss(i,j,bi,bj)   = scal_sss*climsss(i,j,bi,bj)               SSS(i,j,bi,bj)  = exf_outscal_sss*climsss(i,j,bi,bj)
319              enddo              ENDDO
320            enddo            ENDDO
321  #endif  #endif
322    
323  #ifdef ATMOSPHERIC_LOADING  #ifdef ATMOSPHERIC_LOADING
324            do j = jmin,jmax            DO j = jmin,jmax
325              do i = imin,imax              DO i = imin,imax
326                pload(i,j,bi,bj) = scal_apressure*apressure(i,j,bi,bj)               pLoad(i,j,bi,bj)=exf_outscal_apressure*apressure(i,j,bi,bj)
327              enddo              ENDDO
328            enddo            ENDDO
329  #endif  #endif
330    
331    #ifdef ALLOW_SALTFLX
332          enddo            DO j = jmin,jmax
333        enddo              DO i = imin,imax
334                  saltFlux(I,J,bi,bj) = saltflx(I,J,bi,bj)
335  c     Update the tile edges.              ENDDO
336              ENDDO
337        _EXCH_XY_R4(  qnet, mythid )  #endif
338        _EXCH_XY_R4( empmr, mythid )  
339        _EXCH_XY_R4(    fu, mythid )  #ifdef EXF_SEAICE_FRACTION
340        _EXCH_XY_R4(    fv, mythid )            DO j = jmin,jmax
341  #ifdef ALLOW_KPP              DO i = imin,imax
342        _EXCH_XY_R4(   qsw, mythid )                exf_iceFraction(i,j,bi,bj) =
343         &           exf_outscal_areamask*areamask(i,j,bi,bj)
344                  exf_iceFraction(I,J,bi,bj) =
345         &           MIN( MAX(exf_iceFraction(I,J,bi,bj),zeroRS), oneRS )
346                ENDDO
347              ENDDO
348    #endif
349    
350           ENDDO
351          ENDDO
352    
353    C--   Update the tile edges.
354          _EXCH_XY_RS(  Qnet, myThid )
355          _EXCH_XY_RS( EmPmR, myThid )
356           CALL EXCH_UV_XY_RS(fu, fv, .TRUE., myThid)
357    c#if defined(ALLOW_ATM_TEMP) || defined(SHORTWAVE_HEATING)
358    #ifdef SHORTWAVE_HEATING
359    C     Qsw used in SHORTWAVE_HEATING code & for diagnostics (<- EXCH not needed)
360          _EXCH_XY_RS(   Qsw, myThid )
361  #endif  #endif
362  #ifdef ALLOW_CLIMSST_RELAXATION  #ifdef ALLOW_CLIMSST_RELAXATION
363        _EXCH_XY_R4(   sst, mythid )        _EXCH_XY_RS(   SST, myThid )
364  #endif  #endif
365  #ifdef ALLOW_CLIMSSS_RELAXATION  #ifdef ALLOW_CLIMSSS_RELAXATION
366        _EXCH_XY_R4(   sss, mythid )        _EXCH_XY_RS(   SSS, myThid )
367  #endif  #endif
368  #ifdef ATMOSPHERIC_LOADING  #ifdef ATMOSPHERIC_LOADING
369        _EXCH_XY_R4( pload, mythid )        _EXCH_XY_RS( pLoad, myThid )
370    #endif
371    #ifdef EXF_SEAICE_FRACTION
372          _EXCH_XY_RS( exf_iceFraction, myThid )
373  #endif  #endif
374    
375        end        RETURN
376          END
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.32
  ViewVC Help
Powered by ViewVC 1.1.22