/[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.1 by heimbach, Mon May 14 22:08:41 2001 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        subroutine exf_MapFields(  C     !INTERFACE:
13       I                          mythid        SUBROUTINE EXF_MAPFIELDS( myTime, myIter, myThid )
      &                        )  
   
 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     ==================================================================  
14    
15        implicit none  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     == global variables ==  C     !USES:
54          IMPLICIT NONE
55    
56    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_constants.h"  #include "GRID.h"
62  #include "exf_fields.h"  #include "DYNVARS.h"
63  #include "exf_clim_fields.h"  
64    #include "EXF_PARAM.h"
65  c     == routine arguments ==  #include "EXF_CONSTANTS.h"
66    #include "EXF_FIELDS.h"
67  c     mythid - thread number for this instance of the routine.  #ifdef ALLOW_AUTODIFF_TAMC
68    # include "tamc.h"
69        integer mythid  # include "tamc_keys.h"
70    #endif
71  c     == local variables ==  
72    C     !INPUT/OUTPUT PARAMETERS:
73        integer bi,bj  C     myTime  :: Current time in simulation
74        integer i,j  C     myIter  :: Current iteration number
75        integer jtlo  C     myThid  :: my Thread Id number
76        integer jthi        _RL     myTime
77        integer itlo        INTEGER myIter
78        integer ithi        INTEGER myThid
79        integer jmin  
80        integer jmax  C     !LOCAL VARIABLES:
81        integer imin        INTEGER bi,bj
82        integer imax        INTEGER i,j,ks
83        _RL     scal_hfl        INTEGER imin, imax
84        _RL     scal_ust        INTEGER jmin, jmax
85        _RL     scal_vst        PARAMETER ( imin = 1-OLx , imax = sNx+OLx )
86        _RL     scal_swf        PARAMETER ( jmin = 1-OLy , jmax = sNy+OLy )
87        _RL     scal_sst  CEOP
88        _RL     scal_sss  
89  #if (defined (ALLOW_BULKFORMULAE) && defined (ALLOW_ATM_TEMP))  C--   set surface level index:
90        _RL     scal_prc        ks = 1
91  #else  
92        _RL     scal_sfl        DO bj = myByLo(myThid),myByHi(myThid)
93  #endif         DO bi = myBxLo(myThid),myBxHi(myThid)
94    
95  c     == end of interface ==  #ifdef ALLOW_AUTODIFF_TAMC
96              act1 = bi - myBxLo(myThid)
97        jtlo = mybylo(mythid)            max1 = myBxHi(myThid) - myBxLo(myThid) + 1
98        jthi = mybyhi(mythid)            act2 = bj - myByLo(myThid)
99        itlo = mybxlo(mythid)            max2 = myByHi(myThid) - myByLo(myThid) + 1
100        ithi = mybxhi(mythid)            act3 = myThid - 1
101        jmin = 1-oly            max3 = nTx*nTy
102        jmax = sny+oly            act4 = ikey_dynamics - 1
103        imin = 1-olx            ikey = (act1 + 1) + act2*max1
104        imax = snx+olx       &                      + act3*max1*max2
105         &                      + act4*max1*max2*max3
106        scal_hfl =  1. _d 0  #endif /* ALLOW_AUTODIFF_TAMC */
107        scal_ust =  -1. _d 0  
108        scal_vst =  -1. _d 0  C     Heat flux.
109        scal_swf =  1. _d 0            DO j = jmin,jmax
110        scal_sst =  1. _d 0              DO i = imin,imax
111        scal_sss =  1. _d 0               Qnet(i,j,bi,bj) = exf_outscal_hflux*hflux(i,j,bi,bj)
112  #if (defined (ALLOW_BULKFORMULAE) && defined (ALLOW_ATM_TEMP))              ENDDO
113        scal_prc =  1. _d 0            ENDDO
114  #else            IF ( hfluxfile .EQ. ' ' ) THEN
115        scal_sfl =  1. _d 0                 DO j = jmin,jmax
116  #endif              DO i = imin,imax
117                 Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) -
118        do bj = jtlo,jthi       &            exf_outscal_hflux * ( hflux_exfremo_intercept +
119          do bi = itlo,ithi       &            hflux_exfremo_slope*(myTime-startTime) )
120            do j = jmin,jmax              ENDDO
121              do i = imin,imax             ENDDO
122              ENDIF
123  c             Heat flux.  
124                qnet(i,j,bi,bj)  = scal_hfl*hflux(i,j,bi,bj)  C     Freshwater flux.
125              DO j = jmin,jmax
126  c             Salt flux.              DO i = imin,imax
127  #if (defined (ALLOW_BULKFORMULAE)  && defined (ALLOW_ATM_TEMP))               EmPmR(i,j,bi,bj)= exf_outscal_sflux*sflux(i,j,bi,bj)
128                empmr(i,j,bi,bj) = scal_prc*precip(i,j,bi,bj)       &                                          *rhoConstFresh
129  #else              ENDDO
130                empmr(i,j,bi,bj) = scal_sfl*sflux(i,j,bi,bj)            ENDDO
131  #endif            IF ( sfluxfile .EQ. ' ' ) THEN
132               DO j = jmin,jmax
133  c             Zonal wind stress.              DO i = imin,imax
134                fu(i,j,bi,bj)    = scal_ust*ustress(i,j,bi,bj)               EmPmR(i,j,bi,bj) = EmPmR(i,j,bi,bj) - rhoConstFresh*
135         &            exf_outscal_sflux * ( sflux_exfremo_intercept +
136  c             Meridional wind stress.       &            sflux_exfremo_slope*(myTime-startTime) )
137                fv(i,j,bi,bj)    = scal_vst*vstress(i,j,bi,bj)              ENDDO
138               ENDDO
139  #ifdef ALLOW_KPP || (defined (ALLOW_BULKFORMULAE) && defined (ALLOW_ATM_TEMP)))            ENDIF
140  c             Short wave radiative flux.  
141                qsw(i,j,bi,bj)   = scal_swf*swflux(i,j,bi,bj)  #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
221    
222    #ifdef ALLOW_AUTODIFF_TAMC
223    CADJ STORE ustress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
224    #endif
225              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                    ustress(i,j,bi,bj)=windstressmax
230                  ENDIF
231                ENDDO
232              ENDDO
233    #ifdef ALLOW_AUTODIFF_TAMC
234    CADJ STORE ustress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
235    #endif
236              DO j = jmin,jmax
237                DO i = imin,imax
238                  IF (ustress(i,j,bi,bj).LT.-windstressmax) THEN
239                    ustress(i,j,bi,bj)=-windstressmax
240                  ENDIF
241                ENDDO
242              ENDDO
243              IF ( stressIsOnCgrid ) THEN
244               DO j = jmin,jmax
245                DO i = imin+1,imax
246                  fu(i,j,bi,bj) = exf_outscal_ustress*ustress(i,j,bi,bj)
247                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
261    CADJ STORE vstress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
262    #endif
263              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                    vstress(i,j,bi,bj)=windstressmax
268                  ENDIF
269                ENDDO
270              ENDDO
271    #ifdef ALLOW_AUTODIFF_TAMC
272    CADJ STORE vstress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
273    #endif
274              DO j = jmin,jmax
275                DO i = imin,imax
276                  IF (vstress(i,j,bi,bj).LT.-windstressmax) THEN
277                    vstress(i,j,bi,bj)=-windstressmax
278                  ENDIF
279                ENDDO
280              ENDDO
281              IF ( stressIsOnCgrid ) THEN
282               DO j = jmin+1,jmax
283                DO i = imin,imax
284                  fv(i,j,bi,bj) = exf_outscal_vstress*vstress(i,j,bi,bj)
285                ENDDO
286               ENDDO
287              ELSE
288               DO j = jmin+1,jmax
289                DO i = imin,imax
290    C     Shift wind stresses calculated at C-points to W/S points
291                  fv(i,j,bi,bj) = exf_outscal_vstress*
292         &              (vstress(i,j,bi,bj)+vstress(i,j-1,bi,bj))
293         &              *exf_half*maskS(i,j,ks,bi,bj)
294                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                sst(i,j,bi,bj)   = scal_sst*climsst(i,j,bi,bj)            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  #endif  #endif
314    
315  #ifdef ALLOW_CLIMSSS_RELAXATION  #ifdef ALLOW_CLIMSSS_RELAXATION
316                sss(i,j,bi,bj)   = scal_sss*climsss(i,j,bi,bj)            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    #endif
322    
323    #ifdef ATMOSPHERIC_LOADING
324              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    #endif
330    
331    #ifdef ALLOW_SALTFLX
332              DO j = jmin,jmax
333                DO i = imin,imax
334                  saltFlux(I,J,bi,bj) = saltflx(I,J,bi,bj)
335                ENDDO
336              ENDDO
337  #endif  #endif
338    
339              enddo  #ifdef EXF_SEAICE_FRACTION
340            enddo            DO j = jmin,jmax
341          enddo              DO i = imin,imax
342        enddo                exf_iceFraction(i,j,bi,bj) =
343         &           exf_outscal_areamask*areamask(i,j,bi,bj)
344  c     Update the tile edges.                exf_iceFraction(I,J,bi,bj) =
345         &           MIN( MAX(exf_iceFraction(I,J,bi,bj),zeroRS), oneRS )
346        _EXCH_XY_R4(  qnet, mythid )              ENDDO
347        _EXCH_XY_R4( empmr, mythid )            ENDDO
348        _EXCH_XY_R4(    fu, mythid )  #endif
349        _EXCH_XY_R4(    fv, mythid )  
350  #ifdef ALLOW_KPP         ENDDO
351        _EXCH_XY_R4(   qsw, mythid )        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
368    #ifdef ATMOSPHERIC_LOADING
369          _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.1  
changed lines
  Added in v.1.32
  ViewVC Help
Powered by ViewVC 1.1.22