pkg/exf/exf_mapfields.F

C $Header:  $
C $Name:  $

#include "EXF_OPTIONS.h"

      subroutine exf_mapfields( mytime, myiter, mythid )

c     ==================================================================
c     SUBROUTINE exf_mapfields
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     ==================================================================

      implicit none

c     == global variables ==

#include "EEPARAMS.h"
#include "SIZE.h"
#include "PARAMS.h"
#include "FFIELDS.h"
#include "GRID.h"

#include "EXF_PARAM.h"
#include "EXF_CONSTANTS.h"
#include "EXF_FIELDS.h"
#include "EXF_CLIM_PARAM.h"
#include "EXF_CLIM_FIELDS.h"
#ifdef ALLOW_AUTODIFF_TAMC
# include "tamc.h"
# include "tamc_keys.h"
#endif
c     == routine arguments ==

c     mythid - thread number for this instance of the routine.

      integer mythid
      integer myiter
      _RL     mytime

c     == local variables ==

      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

#ifdef ALLOW_AUTODIFF_TAMC
          act1 = bi - myBxLo(myThid)
          max1 = myBxHi(myThid) - myBxLo(myThid) + 1
          act2 = bj - myByLo(myThid)
          max2 = myByHi(myThid) - myByLo(myThid) + 1
          act3 = myThid - 1
          max3 = nTx*nTy
          act4 = ikey_dynamics - 1
          ikey = (act1 + 1) + act2*max1
     &                      + act3*max1*max2
     &                      + act4*max1*max2*max3
#endif /* ALLOW_AUTODIFF_TAMC */

          do j = jmin,jmax
            do i = imin,imax
c            Heat flux.
             qnet(i,j,bi,bj) = exf_outscal_hflux*hflux(i,j,bi,bj)
             if ( hfluxfile .EQ. ' ' ) 
     &            qnet(i,j,bi,bj) = qnet(i,j,bi,bj) -
     &            exf_outscal_hflux * ( hflux_exfremo_intercept +
     &            hflux_exfremo_slope*(mytime-starttime) )
            enddo
          enddo


          do j = jmin,jmax
            do i = imin,imax
c            Salt flux.
             empmr(i,j,bi,bj)= exf_outscal_sflux*sflux(i,j,bi,bj)
             if ( sfluxfile .EQ. ' ' )
     &            empmr(i,j,bi,bj) = empmr(i,j,bi,bj) - 
     &            exf_outscal_sflux * ( sflux_exfremo_intercept +
     &            sflux_exfremo_slope*(mytime-starttime) )
            enddo
          enddo

#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE ustress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
#endif
          do j = jmin,jmax
            do i = imin,imax
c             Zonal wind stress.
              if (ustress(i,j,bi,bj).gt.windstressmax) then
                ustress(i,j,bi,bj)=windstressmax
              endif
            enddo
          enddo
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE ustress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
#endif
          do j = jmin,jmax
            do i = imin,imax
              if (ustress(i,j,bi,bj).lt.-windstressmax) then
                ustress(i,j,bi,bj)=-windstressmax
              endif
            enddo
          enddo
          do j = jmin,jmax
            do i = imin+1,imax
#if (defined (ALLOW_BULKFORMULAE) || defined (USE_EXF_INTERPOLATION))
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
              fu(i,j,bi,bj) = exf_outscal_ustress*ustress(i,j,bi,bj)
#endif
            enddo
          enddo

#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE vstress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
#endif
          do j = jmin,jmax
            do i = imin,imax
c             Meridional wind stress.
              if (vstress(i,j,bi,bj).gt.windstressmax) then
                vstress(i,j,bi,bj)=windstressmax
              endif
            enddo
          enddo
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE vstress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
#endif
          do j = jmin,jmax
            do i = imin,imax
              if (vstress(i,j,bi,bj).lt.-windstressmax) then
                vstress(i,j,bi,bj)=-windstressmax
              endif
            enddo
          enddo
          do j = jmin+1,jmax
            do i = imin,imax
#if (defined (ALLOW_BULKFORMULAE) || defined (USE_EXF_INTERPOLATION))
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
              fv(i,j,bi,bj) = exf_outscal_vstress*vstress(i,j,bi,bj)
#endif
            enddo
          enddo

#ifdef SHORTWAVE_HEATING
c             Short wave radiative flux.
          do j = jmin,jmax
            do i = imin,imax
             qsw(i,j,bi,bj)  = exf_outscal_swflux*swflux(i,j,bi,bj)
            enddo
          enddo
#endif

#ifdef ALLOW_CLIMSST_RELAXATION
          do j = jmin,jmax
            do i = imin,imax
             sst(i,j,bi,bj)  = exf_outscal_sst*climsst(i,j,bi,bj)
            enddo
          enddo
#endif

#ifdef ALLOW_CLIMSSS_RELAXATION
          do j = jmin,jmax
            do i = imin,imax
             sss(i,j,bi,bj)  = exf_outscal_sss*climsss(i,j,bi,bj)
            enddo
          enddo
#endif

#ifdef ATMOSPHERIC_LOADING
          do j = jmin,jmax
            do i = imin,imax
             pload(i,j,bi,bj)=exf_outscal_apressure*apressure(i,j,bi,bj)
            enddo
          enddo
#endif

        enddo
      enddo

c     Update the tile edges.

      _EXCH_XY_R4(  qnet, mythid )
      _EXCH_XY_R4( empmr, mythid )
c      _EXCH_XY_R4(    fu, mythid )
c      _EXCH_XY_R4(    fv, mythid )
       CALL EXCH_UV_XY_RS(fu, fv, .TRUE., myThid)
#ifdef SHORTWAVE_HEATING
      _EXCH_XY_R4(   qsw, mythid )
#endif
#ifdef ALLOW_CLIMSST_RELAXATION
      _EXCH_XY_R4(   sst, mythid )
#endif
#ifdef ALLOW_CLIMSSS_RELAXATION
      _EXCH_XY_R4(   sss, mythid )
#endif
#ifdef ATMOSPHERIC_LOADING
      _EXCH_XY_R4( pload, mythid )
#endif

      end
1	jmc	1.17	C $Header: $
2			C $Name: $
3	heimbach	1.1
4	edhill	1.7	#include "EXF_OPTIONS.h"
5	heimbach	1.1
6	heimbach	1.16	subroutine exf_mapfields( mytime, myiter, mythid )
7	heimbach	1.1
8			c ==================================================================
9	heimbach	1.2	c SUBROUTINE exf_mapfields
10	heimbach	1.1	c ==================================================================
11			c
12	dimitri	1.5	c o Map external forcing fields (ustress, vstress, hflux, sflux,
13			c swflux, apressure, climsss, climsst, etc.) onto ocean model
14			c arrays (fu, fv, Qnet, EmPmR, Qsw, pload, sss, sst, etc.).
15			c This routine is included to separate the ocean state estimation
16			c tool as much as possible from the ocean model. Unit and sign
17			c conventions can be customized using variables exf_outscal_*,
18			c which are set in exf_readparms.F. See the header files
19	jmc	1.17	c EXF_FIELDS.h and FFIELDS.h for definitions of the various input
20	dimitri	1.5	c and output fields and for default unit and sign convetions.
21	heimbach	1.1	c
22			c started: Christian Eckert eckert@mit.edu 09-Aug-1999
23			c
24			c changed: Christian Eckert eckert@mit.edu 11-Jan-2000
25			c - Restructured the code in order to create a package
26			c for the MITgcmUV.
27			c
28			c Christian Eckert eckert@mit.edu 12-Feb-2000
29			c - Changed Routine names (package prefix: exf_)
30			c
31			c Patrick Heimbach, heimbach@mit.edu 06-May-2000
32			c - added and changed CPP flag structure for
33			c ALLOW_BULKFORMULAE, ALLOW_ATM_TEMP
34			c
35			c Patrick Heimbach, heimbach@mit.edu 23-May-2000
36			c - sign change of ustress/vstress incorporated into
37	dimitri	1.5	c scaling factors exf_outscal_ust, exf_outscal_vst
38	dimitri	1.4	c
39	dimitri	1.5	c mods for pkg/seaice: menemenlis@jpl.nasa.gov 20-Dec-2002
40	dimitri	1.4	c
41	heimbach	1.1	c ==================================================================
42	heimbach	1.2	c SUBROUTINE exf_mapfields
43	heimbach	1.1	c ==================================================================
44
45			implicit none
46
47			c == global variables ==
48
49			#include "EEPARAMS.h"
50			#include "SIZE.h"
51	heimbach	1.16	#include "PARAMS.h"
52	heimbach	1.1	#include "FFIELDS.h"
53	heimbach	1.8	#include "GRID.h"
54
55	jmc	1.17	#include "EXF_PARAM.h"
56			#include "EXF_CONSTANTS.h"
57			#include "EXF_FIELDS.h"
58			#include "EXF_CLIM_PARAM.h"
59			#include "EXF_CLIM_FIELDS.h"
60	heimbach	1.2	#ifdef ALLOW_AUTODIFF_TAMC
61			# include "tamc.h"
62			# include "tamc_keys.h"
63			#endif
64	heimbach	1.1	c == routine arguments ==
65
66			c mythid - thread number for this instance of the routine.
67
68			integer mythid
69	heimbach	1.16	integer myiter
70			_RL mytime
71	heimbach	1.1
72			c == local variables ==
73
74			integer bi,bj
75	heimbach	1.13	integer i,j,k
76	heimbach	1.1	integer jtlo
77			integer jthi
78			integer itlo
79			integer ithi
80			integer jmin
81			integer jmax
82			integer imin
83			integer imax
84
85			c == end of interface ==
86
87			jtlo = mybylo(mythid)
88			jthi = mybyhi(mythid)
89			itlo = mybxlo(mythid)
90			ithi = mybxhi(mythid)
91			jmin = 1-oly
92			jmax = sny+oly
93			imin = 1-olx
94			imax = snx+olx
95
96			do bj = jtlo,jthi
97			do bi = itlo,ithi
98	heimbach	1.2
99			#ifdef ALLOW_AUTODIFF_TAMC
100			act1 = bi - myBxLo(myThid)
101			max1 = myBxHi(myThid) - myBxLo(myThid) + 1
102			act2 = bj - myByLo(myThid)
103			max2 = myByHi(myThid) - myByLo(myThid) + 1
104			act3 = myThid - 1
105			max3 = nTx*nTy
106			act4 = ikey_dynamics - 1
107			ikey = (act1 + 1) + act2*max1
108			& + act3max1max2
109			& + act4max1max2*max3
110			#endif /* ALLOW_AUTODIFF_TAMC */
111
112	heimbach	1.1	do j = jmin,jmax
113			do i = imin,imax
114	dimitri	1.5	c Heat flux.
115	heimbach	1.6	qnet(i,j,bi,bj) = exf_outscal_hflux*hflux(i,j,bi,bj)
116	heimbach	1.16	if ( hfluxfile .EQ. ' ' )
117			& qnet(i,j,bi,bj) = qnet(i,j,bi,bj) -
118			& exf_outscal_hflux * ( hflux_exfremo_intercept +
119			& hflux_exfremo_slope*(mytime-starttime) )
120	heimbach	1.2	enddo
121			enddo
122	heimbach	1.1
123	heimbach	1.2
124			do j = jmin,jmax
125			do i = imin,imax
126	dimitri	1.5	c Salt flux.
127	heimbach	1.6	empmr(i,j,bi,bj)= exf_outscal_sflux*sflux(i,j,bi,bj)
128	heimbach	1.16	if ( sfluxfile .EQ. ' ' )
129			& empmr(i,j,bi,bj) = empmr(i,j,bi,bj) -
130			& exf_outscal_sflux * ( sflux_exfremo_intercept +
131			& sflux_exfremo_slope*(mytime-starttime) )
132	heimbach	1.2	enddo
133			enddo
134	heimbach	1.1
135	heimbach	1.2	#ifdef ALLOW_AUTODIFF_TAMC
136			CADJ STORE ustress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
137			#endif
138			do j = jmin,jmax
139			do i = imin,imax
140	heimbach	1.1	c Zonal wind stress.
141	mlosch	1.10	if (ustress(i,j,bi,bj).gt.windstressmax) then
142			ustress(i,j,bi,bj)=windstressmax
143	heimbach	1.2	endif
144			enddo
145			enddo
146			#ifdef ALLOW_AUTODIFF_TAMC
147			CADJ STORE ustress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
148			#endif
149			do j = jmin,jmax
150			do i = imin,imax
151	mlosch	1.10	if (ustress(i,j,bi,bj).lt.-windstressmax) then
152			ustress(i,j,bi,bj)=-windstressmax
153	heimbach	1.2	endif
154			enddo
155			enddo
156			do j = jmin,jmax
157	heimbach	1.8	do i = imin+1,imax
158	dimitri	1.9	#if (defined (ALLOW_BULKFORMULAE) \|\| defined (USE_EXF_INTERPOLATION))
159	heimbach	1.8	c Shift wind stresses calculated at C-points to W/S points
160			fu(i,j,bi,bj) = exf_outscal_ustress*
161			& (ustress(i,j,bi,bj)+ustress(i-1,j,bi,bj))/2.*
162			& maskW(i,j,1,bi,bj)
163			#else
164			fu(i,j,bi,bj) = exf_outscal_ustress*ustress(i,j,bi,bj)
165			#endif
166	heimbach	1.2	enddo
167			enddo
168	heimbach	1.1
169	heimbach	1.2	#ifdef ALLOW_AUTODIFF_TAMC
170			CADJ STORE vstress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
171			#endif
172			do j = jmin,jmax
173			do i = imin,imax
174	heimbach	1.1	c Meridional wind stress.
175	mlosch	1.10	if (vstress(i,j,bi,bj).gt.windstressmax) then
176			vstress(i,j,bi,bj)=windstressmax
177	heimbach	1.2	endif
178			enddo
179			enddo
180			#ifdef ALLOW_AUTODIFF_TAMC
181			CADJ STORE vstress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
182			#endif
183			do j = jmin,jmax
184			do i = imin,imax
185	mlosch	1.10	if (vstress(i,j,bi,bj).lt.-windstressmax) then
186			vstress(i,j,bi,bj)=-windstressmax
187	heimbach	1.2	endif
188			enddo
189			enddo
190	heimbach	1.8	do j = jmin+1,jmax
191	heimbach	1.2	do i = imin,imax
192	dimitri	1.9	#if (defined (ALLOW_BULKFORMULAE) \|\| defined (USE_EXF_INTERPOLATION))
193	heimbach	1.8	c Shift wind stresses calculated at C-points to W/S points
194			fv(i,j,bi,bj) = exf_outscal_vstress*
195			& (vstress(i,j,bi,bj)+vstress(i,j-1,bi,bj))/2.*
196			& maskS(i,j,1,bi,bj)
197			#else
198			fv(i,j,bi,bj) = exf_outscal_vstress*vstress(i,j,bi,bj)
199			#endif
200	heimbach	1.2	enddo
201			enddo
202	heimbach	1.1
203	dimitri	1.5	#ifdef SHORTWAVE_HEATING
204	heimbach	1.1	c Short wave radiative flux.
205	heimbach	1.2	do j = jmin,jmax
206			do i = imin,imax
207	heimbach	1.6	qsw(i,j,bi,bj) = exf_outscal_swflux*swflux(i,j,bi,bj)
208	heimbach	1.2	enddo
209			enddo
210	heimbach	1.1	#endif
211
212			#ifdef ALLOW_CLIMSST_RELAXATION
213	heimbach	1.2	do j = jmin,jmax
214			do i = imin,imax
215	dimitri	1.5	sst(i,j,bi,bj) = exf_outscal_sst*climsst(i,j,bi,bj)
216	heimbach	1.2	enddo
217			enddo
218	heimbach	1.1	#endif
219
220			#ifdef ALLOW_CLIMSSS_RELAXATION
221	heimbach	1.2	do j = jmin,jmax
222			do i = imin,imax
223	dimitri	1.5	sss(i,j,bi,bj) = exf_outscal_sss*climsss(i,j,bi,bj)
224	heimbach	1.2	enddo
225			enddo
226	heimbach	1.1	#endif
227
228	heimbach	1.2	#ifdef ATMOSPHERIC_LOADING
229			do j = jmin,jmax
230			do i = imin,imax
231	dimitri	1.5	pload(i,j,bi,bj)=exf_outscal_apressure*apressure(i,j,bi,bj)
232	heimbach	1.1	enddo
233			enddo
234	heimbach	1.2	#endif
235
236	heimbach	1.1	enddo
237			enddo
238
239			c Update the tile edges.
240
241			_EXCH_XY_R4( qnet, mythid )
242			_EXCH_XY_R4( empmr, mythid )
243	cheisey	1.3	c _EXCH_XY_R4( fu, mythid )
244			c _EXCH_XY_R4( fv, mythid )
245			CALL EXCH_UV_XY_RS(fu, fv, .TRUE., myThid)
246	dimitri	1.5	#ifdef SHORTWAVE_HEATING
247	heimbach	1.1	_EXCH_XY_R4( qsw, mythid )
248			#endif
249			#ifdef ALLOW_CLIMSST_RELAXATION
250			_EXCH_XY_R4( sst, mythid )
251			#endif
252			#ifdef ALLOW_CLIMSSS_RELAXATION
253			_EXCH_XY_R4( sss, mythid )
254	heimbach	1.2	#endif
255			#ifdef ATMOSPHERIC_LOADING
256			_EXCH_XY_R4( pload, mythid )
257	heimbach	1.1	#endif
258
259			end