pkg/exf/exf_mapfields.F

C $Header: /u/gcmpack/MITgcm/pkg/exf/exf_mapfields.F,v 1.17 2007/04/16 23:27:21 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"
#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 (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	jmc	1.18	C $Header: /u/gcmpack/MITgcm/pkg/exf/exf_mapfields.F,v 1.17 2007/04/16 23:27:21 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			#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	jmc	1.18	#if ( ( defined (ALLOW_BULKFORMULAE) && defined (ALLOW_ATM_WIND) ) \
159			\|\| defined (USE_EXF_INTERPOLATION) )
160	heimbach	1.8	c Shift wind stresses calculated at C-points to W/S points
161			fu(i,j,bi,bj) = exf_outscal_ustress*
162			& (ustress(i,j,bi,bj)+ustress(i-1,j,bi,bj))/2.*
163			& maskW(i,j,1,bi,bj)
164			#else
165			fu(i,j,bi,bj) = exf_outscal_ustress*ustress(i,j,bi,bj)
166			#endif
167	heimbach	1.2	enddo
168			enddo
169	heimbach	1.1
170	heimbach	1.2	#ifdef ALLOW_AUTODIFF_TAMC
171			CADJ STORE vstress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
172			#endif
173			do j = jmin,jmax
174			do i = imin,imax
175	heimbach	1.1	c Meridional wind stress.
176	mlosch	1.10	if (vstress(i,j,bi,bj).gt.windstressmax) then
177			vstress(i,j,bi,bj)=windstressmax
178	heimbach	1.2	endif
179			enddo
180			enddo
181			#ifdef ALLOW_AUTODIFF_TAMC
182			CADJ STORE vstress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
183			#endif
184			do j = jmin,jmax
185			do i = imin,imax
186	mlosch	1.10	if (vstress(i,j,bi,bj).lt.-windstressmax) then
187			vstress(i,j,bi,bj)=-windstressmax
188	heimbach	1.2	endif
189			enddo
190			enddo
191	heimbach	1.8	do j = jmin+1,jmax
192	heimbach	1.2	do i = imin,imax
193	jmc	1.18	#if ( ( defined (ALLOW_BULKFORMULAE) && defined (ALLOW_ATM_WIND) ) \
194			\|\| defined (USE_EXF_INTERPOLATION) )
195	heimbach	1.8	c Shift wind stresses calculated at C-points to W/S points
196			fv(i,j,bi,bj) = exf_outscal_vstress*
197			& (vstress(i,j,bi,bj)+vstress(i,j-1,bi,bj))/2.*
198			& maskS(i,j,1,bi,bj)
199			#else
200			fv(i,j,bi,bj) = exf_outscal_vstress*vstress(i,j,bi,bj)
201			#endif
202	heimbach	1.2	enddo
203			enddo
204	heimbach	1.1
205	dimitri	1.5	#ifdef SHORTWAVE_HEATING
206	heimbach	1.1	c Short wave radiative flux.
207	heimbach	1.2	do j = jmin,jmax
208			do i = imin,imax
209	heimbach	1.6	qsw(i,j,bi,bj) = exf_outscal_swflux*swflux(i,j,bi,bj)
210	heimbach	1.2	enddo
211			enddo
212	heimbach	1.1	#endif
213
214			#ifdef ALLOW_CLIMSST_RELAXATION
215	heimbach	1.2	do j = jmin,jmax
216			do i = imin,imax
217	dimitri	1.5	sst(i,j,bi,bj) = exf_outscal_sst*climsst(i,j,bi,bj)
218	heimbach	1.2	enddo
219			enddo
220	heimbach	1.1	#endif
221
222			#ifdef ALLOW_CLIMSSS_RELAXATION
223	heimbach	1.2	do j = jmin,jmax
224			do i = imin,imax
225	dimitri	1.5	sss(i,j,bi,bj) = exf_outscal_sss*climsss(i,j,bi,bj)
226	heimbach	1.2	enddo
227			enddo
228	heimbach	1.1	#endif
229
230	heimbach	1.2	#ifdef ATMOSPHERIC_LOADING
231			do j = jmin,jmax
232			do i = imin,imax
233	dimitri	1.5	pload(i,j,bi,bj)=exf_outscal_apressure*apressure(i,j,bi,bj)
234	heimbach	1.1	enddo
235			enddo
236	heimbach	1.2	#endif
237
238	heimbach	1.1	enddo
239			enddo
240
241			c Update the tile edges.
242
243			_EXCH_XY_R4( qnet, mythid )
244			_EXCH_XY_R4( empmr, mythid )
245	cheisey	1.3	c _EXCH_XY_R4( fu, mythid )
246			c _EXCH_XY_R4( fv, mythid )
247			CALL EXCH_UV_XY_RS(fu, fv, .TRUE., myThid)
248	dimitri	1.5	#ifdef SHORTWAVE_HEATING
249	heimbach	1.1	_EXCH_XY_R4( qsw, mythid )
250			#endif
251			#ifdef ALLOW_CLIMSST_RELAXATION
252			_EXCH_XY_R4( sst, mythid )
253			#endif
254			#ifdef ALLOW_CLIMSSS_RELAXATION
255			_EXCH_XY_R4( sss, mythid )
256	heimbach	1.2	#endif
257			#ifdef ATMOSPHERIC_LOADING
258			_EXCH_XY_R4( pload, mythid )
259	heimbach	1.1	#endif
260
261			end