OpenAD/code_ad_moc/ini_linear_phisurf.F

C $Header: /u/gcmpack/MITgcm_contrib/heimbach/OpenAD/code_ad_moc/ini_linear_phisurf.F,v 1.1 2008/09/15 22:16:00 utke Exp $
C $Name:  $

#include "CPP_OPTIONS.h"

CBOP
C     !ROUTINE: INI_LINEAR_PHISURF
C     !INTERFACE:
      SUBROUTINE INI_LINEAR_PHISURF( myThid )

C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE INI_LINEAR_PHISURF
C     | o Initialise the Linear Relation Phi_surf(eta)
C     *==========================================================*
C     | Initialise -Boyancy at surface level (Bo_surf)
C     |  to setup the Linear relation: Phi_surf(eta)=Bo_surf*eta
C     | Initialise phi0surf = starting point for integrating
C     |                       phiHyd (= phiHyd at r=RoSurf)
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE
C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "SURFACE.h"

C     !INPUT/OUTPUT PARAMETERS:
C     === Routine arguments ===
C     myThid - Thread no. that called this routine.
      INTEGER myThid

C     == Local variables in common ==
C     Hloc  - Temporary array used to write surface topography
C             has to be in common for multi threading
      COMMON / LOCAL_INI_PHISURF / topoHloc
      _RS topoHloc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)

C     !LOCAL VARIABLES:
C     === Local variables ===
C     bi,bj  - Loop counters
C     I,J,K
      INTEGER bi, bj
      INTEGER I, J, K
      _RL     pLoc, rhoLoc
      _RL     dPIdp
CEOP

#ifdef ALLOW_AUTODIFF_TAMC
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO J=1-Oly,sNy+Oly
         DO I=1-Olx,sNx+Olx
           Bo_surf(I,J,bi,bj)  = 0. _d 0
           recip_Bo(I,J,bi,bj) = 0. _d 0
         ENDDO
        ENDDO
       ENDDO
      ENDDO
#endif

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

C-- Initialise -Boyancy at surface level : Bo_surf
C   Bo_surf is defined as d/dr(Phi_surf) and set to g/z2rUnit with
C     z2rUnit = conversion factor from z-unit to r-unit: [z] * z2rUnit = [r]
C   an accurate formulation includes P_surf and T,S_ref effects on rho_surf:
C    (setting uniformLin_PhiSurf=.FALSE.):
C    z-coord (z2rUnit=1): Bo_surf = - Boyancy
C                                 = g * rho_surf(Tref,Sref,pSurf_0)/rho_0
C    p-coord (z2rUnit=rho*g): Bo_surf = 1/rho(Tref(ksurf),pSurf_0)
C   Note: on Phi_surf splitting : Non-linear Time-dependent effects on B_surf
C   [through eta & (T-tRef)_surf] are included in PhiHyd rather than in Bo_surf
C--
      IF ( buoyancyRelation .eq. 'OCEANIC' ) THEN
C-  gBaro = gravity (except for External mode test with reduced gravity)
        DO bj=myByLo(myThid),myByHi(myThid)
         DO bi=myBxLo(myThid),myBxHi(myThid)
          DO J=1-Oly,sNy+Oly
           DO I=1-Olx,sNx+Olx
             Bo_surf(I,J,bi,bj) = gBaro
             recip_Bo(I,J,bi,bj) = 1. _d 0 / gBaro
           ENDDO
          ENDDO
         ENDDO
        ENDDO
      ELSEIF ( uniformLin_PhiSurf ) THEN
C-  use a linear (in ps) uniform relation : Phi'_surf = 1/rhoConst * ps'_surf
        DO bj=myByLo(myThid),myByHi(myThid)
         DO bi=myBxLo(myThid),myBxHi(myThid)
          DO J=1-Oly,sNy+Oly
           DO I=1-Olx,sNx+Olx
c            Bo_surf(I,J,bi,bj)  = rVel2wUnit(1)*gravity
c            recip_Bo(I,J,bi,bj) = wUnit2rVel(1)*recip_gravity
             Bo_surf(I,J,bi,bj)  = recip_rhoConst
             recip_Bo(I,J,bi,bj) = rhoConst
           ENDDO
          ENDDO
         ENDDO
        ENDDO
      ELSEIF ( buoyancyRelation .eq. 'OCEANICP' ) THEN
        DO bj=myByLo(myThid),myByHi(myThid)
         DO bi=myBxLo(myThid),myBxHi(myThid)
          DO J=1-Oly,sNy+Oly
           DO I=1-Olx,sNx+Olx
            IF ( Ro_surf(I,J,bi,bj).GT.0. _d 0
     &          .AND. ksurfC(I,J,bi,bj).LE.Nr ) THEN
             pLoc = Ro_surf(I,J,bi,bj)
             CALL FIND_RHO_SCALAR(
     I            tRef(ksurfC(I,J,bi,bj)), 
     I            sRef(ksurfC(I,J,bi,bj)), 
     I            pLoc,
     O            rhoLoc, myThid )
             rhoLoc = rhoLoc + rhoConst
             IF ( rhoLoc .EQ. 0. _d 0 ) THEN
              Bo_surf(I,J,bi,bj) = 0. _d 0
             ELSE
              Bo_surf(I,J,bi,bj) = 1. _d 0/rhoLoc
             ENDIF
             recip_Bo(I,J,bi,bj) =  rhoLoc
            ELSE
              Bo_surf(I,J,bi,bj)  = 0. _d 0
              recip_Bo(I,J,bi,bj) = 0. _d 0
            ENDIF
           ENDDO
          ENDDO
         ENDDO
        ENDDO
      ELSEIF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN
C-  use a linearized (in ps) Non-uniform relation : Bo_surf(Po_surf,tRef_surf)
C--- Bo = d/d_p(Phi_surf) = tRef_surf*d/d_p(PI) ; PI = Cp*(p/Po)^kappa
        DO bj=myByLo(myThid),myByHi(myThid)
         DO bi=myBxLo(myThid),myBxHi(myThid)
          DO J=1-Oly,sNy+Oly
           DO I=1-Olx,sNx+Olx
            IF ( Ro_surf(I,J,bi,bj).GT.0. _d 0
     &          .AND. ksurfC(I,J,bi,bj).LE.Nr ) THEN
              dPIdp = (atm_Cp*atm_kappa/atm_Po)*
     &         (Ro_surf(I,J,bi,bj)/atm_Po)**(atm_kappa-1. _d 0)
              Bo_surf(I,J,bi,bj) = dPIdp*tRef(ksurfC(I,J,bi,bj))
              recip_Bo(I,J,bi,bj) = 1. _d 0 / Bo_surf(I,J,bi,bj)
            ELSE
              Bo_surf(I,J,bi,bj) = 0.
              recip_Bo(I,J,bi,bj) = 0.
            ENDIF
           ENDDO
          ENDDO
         ENDDO
        ENDDO
      ELSE
        STOP 'INI_LINEAR_PHISURF: We should never reach this point!'
      ENDIF

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

C--   Update overlap regions (jmc: is it really needed ?)
      _EXCH_XY_R8(Bo_surf, myThid)
      _EXCH_XY_R8(recip_Bo, myThid)

      IF ( ( buoyancyRelation .eq. 'ATMOSPHERIC' .OR.
     &       buoyancyRelation .eq. 'OCEANICP'         )
     &       .AND. .NOT.uniformLin_PhiSurf              ) THEN

C--     EXCH (above) contains barrier calls
c       _BARRIER
        CALL WRITE_FLD_XY_RL( 'Bo_surf',' ',Bo_surf,0,myThid)

      ENDIF

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

C--   Initialise phi0surf: used for atmos. surf. P-loading (ocean, z-coord)
C                               or topographic geopotential anom. (p-coord)

        DO bj=myByLo(myThid),myByHi(myThid)
         DO bi=myBxLo(myThid),myBxHi(myThid)
          DO J=1-Oly,sNy+Oly
           DO I=1-Olx,sNx+Olx
             phi0surf(I,J,bi,bj) = 0.
           ENDDO
          ENDDO
         ENDDO
        ENDDO

      IF ( buoyancyRelation .eq. 'ATMOSPHERIC'
     &                   .AND. topoFile.NE.' ' ) THEN

#ifdef ALLOW_AUTODIFF_TAMC
         STOP 'CANNOT PRESENTLY USE THIS OPTION WITH ADJOINT'
#else

C--   Compute topoH = PhiRef(Po_surf)/g ; is different from original
C      topoZ(read from file) because of truncation of Po_surf.
C     NOTE: not clear for now which topoZ needs to be saved in common block
C--   AND set phi0surf = starting point for integrating Geopotential;

        CALL INI_P_GROUND( -2,
     O                     topoHloc,
     I                     Ro_surf, myThid )


       IF (selectFindRoSurf.NE.0) THEN
        _EXCH_XY_RS(phi0surf, myThid)

        CALL WRITE_FLD_XY_RS( 'phi0surf',' ',phi0surf,0,myThid)
       ELSE
        _BARRIER
       ENDIF

        CALL WRITE_FLD_XY_RS( 'topo_H',' ',topoHloc,0,myThid)

#endif /* ALLOW_AUTODIFF_TAMC */

      ENDIF

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
      RETURN
      END
1	utke	1.2	C $Header: /u/gcmpack/MITgcm_contrib/heimbach/OpenAD/code_ad_moc/ini_linear_phisurf.F,v 1.1 2008/09/15 22:16:00 utke Exp $
2	utke	1.1	C $Name: $
3
4			#include "CPP_OPTIONS.h"
5
6			CBOP
7			C !ROUTINE: INI_LINEAR_PHISURF
8			C !INTERFACE:
9			SUBROUTINE INI_LINEAR_PHISURF( myThid )
10
11			C !DESCRIPTION: \bv
12			C ==========================================================
13			C \| SUBROUTINE INI_LINEAR_PHISURF
14			C \| o Initialise the Linear Relation Phi_surf(eta)
15			C ==========================================================
16			C \| Initialise -Boyancy at surface level (Bo_surf)
17			C \| to setup the Linear relation: Phi_surf(eta)=Bo_surf*eta
18			C \| Initialise phi0surf = starting point for integrating
19			C \| phiHyd (= phiHyd at r=RoSurf)
20			C ==========================================================
21			C \ev
22
23			C !USES:
24			IMPLICIT NONE
25			C === Global variables ===
26			#include "SIZE.h"
27			#include "EEPARAMS.h"
28			#include "PARAMS.h"
29			#include "GRID.h"
30			#include "SURFACE.h"
31
32			C !INPUT/OUTPUT PARAMETERS:
33			C === Routine arguments ===
34			C myThid - Thread no. that called this routine.
35			INTEGER myThid
36
37			C == Local variables in common ==
38			C Hloc - Temporary array used to write surface topography
39			C has to be in common for multi threading
40			COMMON / LOCAL_INI_PHISURF / topoHloc
41			_RS topoHloc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
42
43			C !LOCAL VARIABLES:
44			C === Local variables ===
45			C bi,bj - Loop counters
46			C I,J,K
47			INTEGER bi, bj
48			INTEGER I, J, K
49			_RL pLoc, rhoLoc
50			_RL dPIdp
51			CEOP
52
53			#ifdef ALLOW_AUTODIFF_TAMC
54			DO bj=myByLo(myThid),myByHi(myThid)
55			DO bi=myBxLo(myThid),myBxHi(myThid)
56			DO J=1-Oly,sNy+Oly
57			DO I=1-Olx,sNx+Olx
58			Bo_surf(I,J,bi,bj) = 0. _d 0
59			recip_Bo(I,J,bi,bj) = 0. _d 0
60			ENDDO
61			ENDDO
62			ENDDO
63			ENDDO
64			#endif
65
66			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
67
68			C-- Initialise -Boyancy at surface level : Bo_surf
69			C Bo_surf is defined as d/dr(Phi_surf) and set to g/z2rUnit with
70			C z2rUnit = conversion factor from z-unit to r-unit: [z] * z2rUnit = [r]
71			C an accurate formulation includes P_surf and T,S_ref effects on rho_surf:
72			C (setting uniformLin_PhiSurf=.FALSE.):
73			C z-coord (z2rUnit=1): Bo_surf = - Boyancy
74			C = g * rho_surf(Tref,Sref,pSurf_0)/rho_0
75			C p-coord (z2rUnit=rho*g): Bo_surf = 1/rho(Tref(ksurf),pSurf_0)
76			C Note: on Phi_surf splitting : Non-linear Time-dependent effects on B_surf
77			C [through eta & (T-tRef)_surf] are included in PhiHyd rather than in Bo_surf
78			C--
79			IF ( buoyancyRelation .eq. 'OCEANIC' ) THEN
80			C- gBaro = gravity (except for External mode test with reduced gravity)
81			DO bj=myByLo(myThid),myByHi(myThid)
82			DO bi=myBxLo(myThid),myBxHi(myThid)
83			DO J=1-Oly,sNy+Oly
84			DO I=1-Olx,sNx+Olx
85			Bo_surf(I,J,bi,bj) = gBaro
86			recip_Bo(I,J,bi,bj) = 1. _d 0 / gBaro
87			ENDDO
88			ENDDO
89			ENDDO
90			ENDDO
91			ELSEIF ( uniformLin_PhiSurf ) THEN
92			C- use a linear (in ps) uniform relation : Phi'_surf = 1/rhoConst * ps'_surf
93			DO bj=myByLo(myThid),myByHi(myThid)
94			DO bi=myBxLo(myThid),myBxHi(myThid)
95			DO J=1-Oly,sNy+Oly
96			DO I=1-Olx,sNx+Olx
97			c Bo_surf(I,J,bi,bj) = rVel2wUnit(1)*gravity
98			c recip_Bo(I,J,bi,bj) = wUnit2rVel(1)*recip_gravity
99			Bo_surf(I,J,bi,bj) = recip_rhoConst
100			recip_Bo(I,J,bi,bj) = rhoConst
101			ENDDO
102			ENDDO
103			ENDDO
104			ENDDO
105			ELSEIF ( buoyancyRelation .eq. 'OCEANICP' ) THEN
106			DO bj=myByLo(myThid),myByHi(myThid)
107			DO bi=myBxLo(myThid),myBxHi(myThid)
108			DO J=1-Oly,sNy+Oly
109			DO I=1-Olx,sNx+Olx
110			IF ( Ro_surf(I,J,bi,bj).GT.0. _d 0
111			& .AND. ksurfC(I,J,bi,bj).LE.Nr ) THEN
112			pLoc = Ro_surf(I,J,bi,bj)
113			CALL FIND_RHO_SCALAR(
114			I tRef(ksurfC(I,J,bi,bj)),
115			I sRef(ksurfC(I,J,bi,bj)),
116			I pLoc,
117			O rhoLoc, myThid )
118			rhoLoc = rhoLoc + rhoConst
119			IF ( rhoLoc .EQ. 0. _d 0 ) THEN
120			Bo_surf(I,J,bi,bj) = 0. _d 0
121			ELSE
122			Bo_surf(I,J,bi,bj) = 1. _d 0/rhoLoc
123			ENDIF
124			recip_Bo(I,J,bi,bj) = rhoLoc
125			ELSE
126			Bo_surf(I,J,bi,bj) = 0. _d 0
127			recip_Bo(I,J,bi,bj) = 0. _d 0
128			ENDIF
129			ENDDO
130			ENDDO
131			ENDDO
132			ENDDO
133			ELSEIF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN
134			C- use a linearized (in ps) Non-uniform relation : Bo_surf(Po_surf,tRef_surf)
135			C--- Bo = d/d_p(Phi_surf) = tRef_surfd/d_p(PI) ; PI = Cp(p/Po)^kappa
136			DO bj=myByLo(myThid),myByHi(myThid)
137			DO bi=myBxLo(myThid),myBxHi(myThid)
138			DO J=1-Oly,sNy+Oly
139			DO I=1-Olx,sNx+Olx
140			IF ( Ro_surf(I,J,bi,bj).GT.0. _d 0
141			& .AND. ksurfC(I,J,bi,bj).LE.Nr ) THEN
142			dPIdp = (atm_Cpatm_kappa/atm_Po)
143			& (Ro_surf(I,J,bi,bj)/atm_Po)**(atm_kappa-1. _d 0)
144			Bo_surf(I,J,bi,bj) = dPIdp*tRef(ksurfC(I,J,bi,bj))
145			recip_Bo(I,J,bi,bj) = 1. _d 0 / Bo_surf(I,J,bi,bj)
146			ELSE
147			Bo_surf(I,J,bi,bj) = 0.
148			recip_Bo(I,J,bi,bj) = 0.
149			ENDIF
150			ENDDO
151			ENDDO
152			ENDDO
153			ENDDO
154			ELSE
155			STOP 'INI_LINEAR_PHISURF: We should never reach this point!'
156			ENDIF
157
158			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
159
160			C-- Update overlap regions (jmc: is it really needed ?)
161			_EXCH_XY_R8(Bo_surf, myThid)
162			_EXCH_XY_R8(recip_Bo, myThid)
163
164			IF ( ( buoyancyRelation .eq. 'ATMOSPHERIC' .OR.
165			& buoyancyRelation .eq. 'OCEANICP' )
166			& .AND. .NOT.uniformLin_PhiSurf ) THEN
167
168			C-- EXCH (above) contains barrier calls
169			c _BARRIER
170			CALL WRITE_FLD_XY_RL( 'Bo_surf',' ',Bo_surf,0,myThid)
171
172			ENDIF
173
174			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
175
176			C-- Initialise phi0surf: used for atmos. surf. P-loading (ocean, z-coord)
177			C or topographic geopotential anom. (p-coord)
178
179			DO bj=myByLo(myThid),myByHi(myThid)
180			DO bi=myBxLo(myThid),myBxHi(myThid)
181			DO J=1-Oly,sNy+Oly
182			DO I=1-Olx,sNx+Olx
183			phi0surf(I,J,bi,bj) = 0.
184			ENDDO
185			ENDDO
186			ENDDO
187			ENDDO
188
189			IF ( buoyancyRelation .eq. 'ATMOSPHERIC'
190			& .AND. topoFile.NE.' ' ) THEN
191
192			#ifdef ALLOW_AUTODIFF_TAMC
193			STOP 'CANNOT PRESENTLY USE THIS OPTION WITH ADJOINT'
194			#else
195
196			C-- Compute topoH = PhiRef(Po_surf)/g ; is different from original
197			C topoZ(read from file) because of truncation of Po_surf.
198			C NOTE: not clear for now which topoZ needs to be saved in common block
199			C-- AND set phi0surf = starting point for integrating Geopotential;
200
201			CALL INI_P_GROUND( -2,
202			O topoHloc,
203			I Ro_surf, myThid )
204
205
206			IF (selectFindRoSurf.NE.0) THEN
207			_EXCH_XY_RS(phi0surf, myThid)
208
209			CALL WRITE_FLD_XY_RS( 'phi0surf',' ',phi0surf,0,myThid)
210			ELSE
211			_BARRIER
212			ENDIF
213
214			CALL WRITE_FLD_XY_RS( 'topo_H',' ',topoHloc,0,myThid)
215
216			#endif /* ALLOW_AUTODIFF_TAMC */
217
218			ENDIF
219
220			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
221			RETURN
222			END