pkg/exf/exf_mapfields.F

C $Header: /u/gcmpack/MITgcm/pkg/exf/exf_mapfields.F,v 1.18 2007/04/17 18:53:32 jmc Exp $
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"
#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 (ALLOW_ATM_WIND) ) \
     || 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 (ALLOW_ATM_WIND) ) \
     || 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	heimbach	1.19	C $Header: /u/gcmpack/MITgcm/pkg/exf/exf_mapfields.F,v 1.18 2007/04/17 18:53:32 jmc Exp $
2	jmc	1.17	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	heimbach	1.2	#ifdef ALLOW_AUTODIFF_TAMC
59			# include "tamc.h"
60			# include "tamc_keys.h"
61			#endif
62	heimbach	1.1	c == routine arguments ==
63
64			c mythid - thread number for this instance of the routine.
65
66			integer mythid
67	heimbach	1.16	integer myiter
68			_RL mytime
69	heimbach	1.1
70			c == local variables ==
71
72			integer bi,bj
73	heimbach	1.13	integer i,j,k
74	heimbach	1.1	integer jtlo
75			integer jthi
76			integer itlo
77			integer ithi
78			integer jmin
79			integer jmax
80			integer imin
81			integer imax
82
83			c == end of interface ==
84
85			jtlo = mybylo(mythid)
86			jthi = mybyhi(mythid)
87			itlo = mybxlo(mythid)
88			ithi = mybxhi(mythid)
89			jmin = 1-oly
90			jmax = sny+oly
91			imin = 1-olx
92			imax = snx+olx
93
94			do bj = jtlo,jthi
95			do bi = itlo,ithi
96	heimbach	1.2
97			#ifdef ALLOW_AUTODIFF_TAMC
98			act1 = bi - myBxLo(myThid)
99			max1 = myBxHi(myThid) - myBxLo(myThid) + 1
100			act2 = bj - myByLo(myThid)
101			max2 = myByHi(myThid) - myByLo(myThid) + 1
102			act3 = myThid - 1
103			max3 = nTx*nTy
104			act4 = ikey_dynamics - 1
105			ikey = (act1 + 1) + act2*max1
106			& + act3max1max2
107			& + act4max1max2*max3
108			#endif /* ALLOW_AUTODIFF_TAMC */
109
110	heimbach	1.1	do j = jmin,jmax
111			do i = imin,imax
112	dimitri	1.5	c Heat flux.
113	heimbach	1.6	qnet(i,j,bi,bj) = exf_outscal_hflux*hflux(i,j,bi,bj)
114	heimbach	1.16	if ( hfluxfile .EQ. ' ' )
115			& qnet(i,j,bi,bj) = qnet(i,j,bi,bj) -
116			& exf_outscal_hflux * ( hflux_exfremo_intercept +
117			& hflux_exfremo_slope*(mytime-starttime) )
118	heimbach	1.2	enddo
119			enddo
120	heimbach	1.1
121	heimbach	1.2
122			do j = jmin,jmax
123			do i = imin,imax
124	dimitri	1.5	c Salt flux.
125	heimbach	1.6	empmr(i,j,bi,bj)= exf_outscal_sflux*sflux(i,j,bi,bj)
126	heimbach	1.16	if ( sfluxfile .EQ. ' ' )
127			& empmr(i,j,bi,bj) = empmr(i,j,bi,bj) -
128			& exf_outscal_sflux * ( sflux_exfremo_intercept +
129			& sflux_exfremo_slope*(mytime-starttime) )
130	heimbach	1.2	enddo
131			enddo
132	heimbach	1.1
133	heimbach	1.2	#ifdef ALLOW_AUTODIFF_TAMC
134			CADJ STORE ustress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
135			#endif
136			do j = jmin,jmax
137			do i = imin,imax
138	heimbach	1.1	c Zonal wind stress.
139	mlosch	1.10	if (ustress(i,j,bi,bj).gt.windstressmax) then
140			ustress(i,j,bi,bj)=windstressmax
141	heimbach	1.2	endif
142			enddo
143			enddo
144			#ifdef ALLOW_AUTODIFF_TAMC
145			CADJ STORE ustress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
146			#endif
147			do j = jmin,jmax
148			do i = imin,imax
149	mlosch	1.10	if (ustress(i,j,bi,bj).lt.-windstressmax) then
150			ustress(i,j,bi,bj)=-windstressmax
151	heimbach	1.2	endif
152			enddo
153			enddo
154			do j = jmin,jmax
155	heimbach	1.8	do i = imin+1,imax
156	jmc	1.18	#if ( ( defined (ALLOW_BULKFORMULAE) && defined (ALLOW_ATM_WIND) ) \
157			\|\| defined (USE_EXF_INTERPOLATION) )
158	heimbach	1.8	c Shift wind stresses calculated at C-points to W/S points
159			fu(i,j,bi,bj) = exf_outscal_ustress*
160			& (ustress(i,j,bi,bj)+ustress(i-1,j,bi,bj))/2.*
161			& maskW(i,j,1,bi,bj)
162			#else
163			fu(i,j,bi,bj) = exf_outscal_ustress*ustress(i,j,bi,bj)
164			#endif
165	heimbach	1.2	enddo
166			enddo
167	heimbach	1.1
168	heimbach	1.2	#ifdef ALLOW_AUTODIFF_TAMC
169			CADJ STORE vstress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
170			#endif
171			do j = jmin,jmax
172			do i = imin,imax
173	heimbach	1.1	c Meridional wind stress.
174	mlosch	1.10	if (vstress(i,j,bi,bj).gt.windstressmax) then
175			vstress(i,j,bi,bj)=windstressmax
176	heimbach	1.2	endif
177			enddo
178			enddo
179			#ifdef ALLOW_AUTODIFF_TAMC
180			CADJ STORE vstress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
181			#endif
182			do j = jmin,jmax
183			do i = imin,imax
184	mlosch	1.10	if (vstress(i,j,bi,bj).lt.-windstressmax) then
185			vstress(i,j,bi,bj)=-windstressmax
186	heimbach	1.2	endif
187			enddo
188			enddo
189	heimbach	1.8	do j = jmin+1,jmax
190	heimbach	1.2	do i = imin,imax
191	jmc	1.18	#if ( ( defined (ALLOW_BULKFORMULAE) && defined (ALLOW_ATM_WIND) ) \
192			\|\| 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