/[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.14 by heimbach, Mon Jan 2 21:17:02 2006 UTC revision 1.27 by jmc, Thu Apr 4 16:30:50 2013 UTC
# Line 1  Line 1 
1  c $Header$  C $Header$
2    C $Name$
3    
4  #include "EXF_OPTIONS.h"  #include "EXF_OPTIONS.h"
5    
6        subroutine exf_mapfields( 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"
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"
63  #include "exf_fields.h"  #include "EXF_FIELDS.h"
 #include "exf_clim_param.h"  
 #include "exf_clim_fields.h"  
64  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
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    C     myIter  :: Current iteration number
72    C     myThid  :: my Thread Id number
73          _RL     myTime
74          INTEGER myIter
75          INTEGER myThid
76    
77    C     !LOCAL VARIABLES:
78          INTEGER bi,bj
79          INTEGER i,j,ks
80          INTEGER imin, imax
81          INTEGER jmin, jmax
82          PARAMETER ( imin = 1-OLx , imax = sNx+OLx )
83          PARAMETER ( jmin = 1-OLy , jmax = sNy+OLy )
84    CEOP
85    
86        integer mythid  C--   set surface level index:
87          ks = 1
88    
89  c     == local variables ==        DO bj = myByLo(myThid),myByHi(myThid)
90           DO bi = myBxLo(myThid),myBxHi(myThid)
       integer bi,bj  
       integer i,j,k  
       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  
   
       do bj = jtlo,jthi  
         do bi = itlo,ithi  
91    
92  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
93            act1 = bi - myBxLo(myThid)            act1 = bi - myBxLo(myThid)
# Line 105  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              enddo              ENDDO
110            enddo            ENDDO
111              IF ( hfluxfile .EQ. ' ' ) THEN
112               DO j = jmin,jmax
113            do j = jmin,jmax              DO i = imin,imax
114              do i = imin,imax               Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) -
115  c            Salt flux.       &            exf_outscal_hflux * ( hflux_exfremo_intercept +
116               empmr(i,j,bi,bj)= exf_outscal_sflux*sflux(i,j,bi,bj)       &            hflux_exfremo_slope*(myTime-startTime) )
117              enddo              ENDDO
118            enddo             ENDDO
119              ENDIF
120    
121    C     Salt flux.
122              DO j = jmin,jmax
123                DO i = imin,imax
124                 EmPmR(i,j,bi,bj)= exf_outscal_sflux*sflux(i,j,bi,bj)
125         &                                          *rhoConstFresh
126                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 +
133         &            sflux_exfremo_slope*(myTime-startTime) )
134                ENDDO
135               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) Run-Off comes at the temp of surface water (with same Cp)
146    C     5) Evap is released to the Atmos @ surf-temp (=SST); should be using
147    C        the water-vapor heat capacity here and consistently in Bulk-Formulae;
148    C        Could also be put directly into Latent Heat flux.
149               IF ( snowPrecipFile .NE. ' ' ) THEN
150    C--   Melt snow (if provided) into the ocean and account for rain-temp
151                DO j = 1, sNy
152                 DO i = 1, sNx
153                  Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)
154         &              + flami*snowPrecip(i,j,bi,bj)*rhoConstFresh
155         &              - HeatCapacity_Cp
156         &               *( atemp(i,j,bi,bj) - cen2kel - temp_EvPrRn )
157         &               *( precip(i,j,bi,bj)- snowPrecip(i,j,bi,bj) )
158         &               *rhoConstFresh
159                 ENDDO
160                ENDDO
161               ELSE
162    C--   Make snow (according to Air Temp) and melt it in the ocean
163    C     note: here we just use the same criteria as over seaice but would be
164    C           better to consider a higher altitude air temp, e.g., 850.mb
165                DO j = 1, sNy
166                 DO i = 1, sNx
167                  IF ( atemp(i,j,bi,bj).LT.cen2kel ) THEN
168                   Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)
169         &              + flami*precip(i,j,bi,bj)*rhoConstFresh
170                  ELSE
171    C--   Account for rain-temp
172                   Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)
173         &              - HeatCapacity_Cp
174         &               *( atemp(i,j,bi,bj) - cen2kel - temp_EvPrRn )
175         &               *precip(i,j,bi,bj)*rhoConstFresh
176                  ENDIF
177                 ENDDO
178                ENDDO
179               ENDIF
180    #ifdef ALLOW_RUNOFF
181    C--   Account for energy content of RunOff:
182               DO j = 1, sNy
183                DO i = 1, sNx
184                  Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)
185         &              - HeatCapacity_Cp
186         &               *( theta(i,j,ks,bi,bj) - temp_EvPrRn )
187         &               *runoff(i,j,bi,bj)*rhoConstFresh
188                ENDDO
189               ENDDO
190    #endif
191    C--   Account for energy content of Evap:
192               DO j = 1, sNy
193                DO i = 1, sNx
194                  Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)
195         &              + HeatCapacity_Cp
196         &               *( theta(i,j,ks,bi,bj) - temp_EvPrRn )
197         &               *evap(i,j,bi,bj)*rhoConstFresh
198                  Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)*maskC(i,j,ks,bi,bj)
199                ENDDO
200               ENDDO
201              ENDIF
202    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
203    #endif /* ALLOW_ATM_TEMP */
204    
205  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
206  CADJ STORE ustress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE ustress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
207  #endif  #endif
208            do j = jmin,jmax            DO j = jmin,jmax
209              do i = imin,imax              DO i = imin,imax
210  c             Zonal wind stress.  C             Zonal wind stress.
211                if (ustress(i,j,bi,bj).gt.windstressmax) then                IF (ustress(i,j,bi,bj).GT.windstressmax) THEN
212                  ustress(i,j,bi,bj)=windstressmax                  ustress(i,j,bi,bj)=windstressmax
213                endif                ENDIF
214              enddo              ENDDO
215            enddo            ENDDO
216  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
217  CADJ STORE ustress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE ustress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
218  #endif  #endif
219            do j = jmin,jmax            DO j = jmin,jmax
220              do i = imin,imax              DO i = imin,imax
221                if (ustress(i,j,bi,bj).lt.-windstressmax) then                IF (ustress(i,j,bi,bj).LT.-windstressmax) THEN
222                  ustress(i,j,bi,bj)=-windstressmax                  ustress(i,j,bi,bj)=-windstressmax
223                endif                ENDIF
224              enddo              ENDDO
225            enddo            ENDDO
226            do j = jmin,jmax            IF ( stressIsOnCgrid ) THEN
227              do i = imin+1,imax             DO j = jmin,jmax
228  #if (defined (ALLOW_BULKFORMULAE) || defined (USE_EXF_INTERPOLATION))              DO i = imin+1,imax
 c     Shift wind stresses calculated at C-points to W/S points  
               fu(i,j,bi,bj) = exf_outscal_ustress*  
      &              (ustress(i,j,bi,bj)+ustress(i-1,j,bi,bj))/2.*  
      &              maskW(i,j,1,bi,bj)  
 #else  
229                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)
230  #endif              ENDDO
231              enddo             ENDDO
232            enddo            ELSE
233               DO j = jmin,jmax
234                DO i = imin+1,imax
235    C     Shift wind stresses calculated at Grid-center to W/S points
236                  fu(i,j,bi,bj) = exf_outscal_ustress*
237         &              (ustress(i,j,bi,bj)+ustress(i-1,j,bi,bj))
238         &              *exf_half*maskW(i,j,ks,bi,bj)
239                ENDDO
240               ENDDO
241              ENDIF
242    
243  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
244  CADJ STORE vstress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE vstress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
245  #endif  #endif
246            do j = jmin,jmax            DO j = jmin,jmax
247              do i = imin,imax              DO i = imin,imax
248  c             Meridional wind stress.  C             Meridional wind stress.
249                if (vstress(i,j,bi,bj).gt.windstressmax) then                IF (vstress(i,j,bi,bj).GT.windstressmax) THEN
250                  vstress(i,j,bi,bj)=windstressmax                  vstress(i,j,bi,bj)=windstressmax
251                endif                ENDIF
252              enddo              ENDDO
253            enddo            ENDDO
254  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
255  CADJ STORE vstress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE vstress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
256  #endif  #endif
257            do j = jmin,jmax            DO j = jmin,jmax
258              do i = imin,imax              DO i = imin,imax
259                if (vstress(i,j,bi,bj).lt.-windstressmax) then                IF (vstress(i,j,bi,bj).LT.-windstressmax) THEN
260                  vstress(i,j,bi,bj)=-windstressmax                  vstress(i,j,bi,bj)=-windstressmax
261                endif                ENDIF
262              enddo              ENDDO
263            enddo            ENDDO
264            do j = jmin+1,jmax            IF ( stressIsOnCgrid ) THEN
265              do i = imin,imax             DO j = jmin+1,jmax
266  #if (defined (ALLOW_BULKFORMULAE) || defined (USE_EXF_INTERPOLATION))              DO i = imin,imax
 c     Shift wind stresses calculated at C-points to W/S points  
               fv(i,j,bi,bj) = exf_outscal_vstress*  
      &              (vstress(i,j,bi,bj)+vstress(i,j-1,bi,bj))/2.*  
      &              maskS(i,j,1,bi,bj)  
 #else  
267                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)
268  #endif              ENDDO
269              enddo             ENDDO
270            enddo            ELSE
271               DO j = jmin+1,jmax
272  #ifdef SHORTWAVE_HEATING              DO i = imin,imax
273  c             Short wave radiative flux.  C     Shift wind stresses calculated at C-points to W/S points
274            do j = jmin,jmax                fv(i,j,bi,bj) = exf_outscal_vstress*
275              do i = imin,imax       &              (vstress(i,j,bi,bj)+vstress(i,j-1,bi,bj))
276               qsw(i,j,bi,bj)  = exf_outscal_swflux*swflux(i,j,bi,bj)       &              *exf_half*maskS(i,j,ks,bi,bj)
277              enddo              ENDDO
278            enddo             ENDDO
279              ENDIF
280    
281    #if defined(ALLOW_ATM_TEMP) || defined(SHORTWAVE_HEATING)
282    C             Short wave radiative flux.
283              DO j = jmin,jmax
284                DO i = imin,imax
285                 Qsw(i,j,bi,bj)  = exf_outscal_swflux*swflux(i,j,bi,bj)
286                ENDDO
287              ENDDO
288  #endif  #endif
289    
290  #ifdef ALLOW_CLIMSST_RELAXATION  #ifdef ALLOW_CLIMSST_RELAXATION
291            do j = jmin,jmax            DO j = jmin,jmax
292              do i = imin,imax              DO i = imin,imax
293               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)
294              enddo              ENDDO
295            enddo            ENDDO
296  #endif  #endif
297    
298  #ifdef ALLOW_CLIMSSS_RELAXATION  #ifdef ALLOW_CLIMSSS_RELAXATION
299            do j = jmin,jmax            DO j = jmin,jmax
300              do i = imin,imax              DO i = imin,imax
301               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)
302              enddo              ENDDO
303            enddo            ENDDO
304  #endif  #endif
305    
306  #ifdef ATMOSPHERIC_LOADING  #ifdef ATMOSPHERIC_LOADING
307            do j = jmin,jmax            DO j = jmin,jmax
308              do i = imin,imax              DO i = imin,imax
309               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)
310              enddo              ENDDO
311            enddo            ENDDO
312  #endif  #endif
313    
314          enddo  #ifdef EXF_ALLOW_SEAICE_RELAX
315        enddo            DO j = jmin,jmax
316                DO i = imin,imax
317  c     Update the tile edges.                obsSIce(i,j,bi,bj) =
318         &           exf_outscal_areamask*areamask(i,j,bi,bj)
319        _EXCH_XY_R4(  qnet, mythid )                obsSIce(I,J,bi,bj) =
320        _EXCH_XY_R4( empmr, mythid )       &           MIN(MAX(obsSIce(I,J,bi,bj), 0.d0 ), 1.d0)
321  c      _EXCH_XY_R4(    fu, mythid )              ENDDO
322  c      _EXCH_XY_R4(    fv, mythid )            ENDDO
323    #endif
324    
325           ENDDO
326          ENDDO
327    
328    C--   Update the tile edges.
329          _EXCH_XY_RS(  Qnet, myThid )
330          _EXCH_XY_RS( EmPmR, myThid )
331         CALL EXCH_UV_XY_RS(fu, fv, .TRUE., myThid)         CALL EXCH_UV_XY_RS(fu, fv, .TRUE., myThid)
332    c#if defined(ALLOW_ATM_TEMP) || defined(SHORTWAVE_HEATING)
333  #ifdef SHORTWAVE_HEATING  #ifdef SHORTWAVE_HEATING
334        _EXCH_XY_R4(   qsw, mythid )  C     Qsw used in SHORTWAVE_HEATING code & for diagnostics (<- EXCH not needed)
335          _EXCH_XY_RS(   Qsw, myThid )
336  #endif  #endif
337  #ifdef ALLOW_CLIMSST_RELAXATION  #ifdef ALLOW_CLIMSST_RELAXATION
338        _EXCH_XY_R4(   sst, mythid )        _EXCH_XY_RS(   SST, myThid )
339  #endif  #endif
340  #ifdef ALLOW_CLIMSSS_RELAXATION  #ifdef ALLOW_CLIMSSS_RELAXATION
341        _EXCH_XY_R4(   sss, mythid )        _EXCH_XY_RS(   SSS, myThid )
 #endif  
 #ifdef ALLOW_CLIMTEMP_RELAXATION  
       _EXCH_XYZ_R4( thetaStar, mythid )  
 #endif  
 #ifdef ALLOW_CLIMSALT_RELAXATION  
       _EXCH_XYZ_R4( saltStar, mythid )  
342  #endif  #endif
343  #ifdef ATMOSPHERIC_LOADING  #ifdef ATMOSPHERIC_LOADING
344        _EXCH_XY_R4( pload, mythid )        _EXCH_XY_RS( pLoad, myThid )
345    #endif
346    #ifdef EXF_ALLOW_SEAICE_RELAX
347          _EXCH_XY_RS( obsSIce, myThid )
348  #endif  #endif
349    
350        end        RETURN
351          END
Legend:
Removed from v.1.14  
changed lines
  Added in v.1.27
  ViewVC Help
Powered by ViewVC 1.1.22