model/src/ini_linear_phisurf.F

C $Header: /u/gcmpack/MITgcm/model/src/ini_linear_phisurf.F,v 1.16 2009/06/14 21:50:26 jmc Exp $
C $Name:  $

#include "PACKAGES_CONFIG.h"
#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     | Compute global domain Area ;
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"
#ifdef ALLOW_EXCH2
# include "W2_EXCH2_SIZE.h"
# include "W2_EXCH2_TOPOLOGY.h"
#endif /* ALLOW_EXCH2 */


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

C     == Local variables in common ==
C     topoHloc had to be in common for multi threading but no longer
C     needed since MDSIO now allows (2009/06/07) to write local arrays
      _RL tileArea(nSx,nSy), threadArea
C     put tileArea in (local) common block to print from master-thread:
      COMMON / LOCAL_INI_PHISURF / tileArea

C     !LOCAL VARIABLES:
C     === Local variables ===
C     topoHloc  :: Temporary array used to write surface topography
C     bi,bj  :: tile indices
C     I,J,K  :: Loop counters
      _RS topoHloc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      INTEGER bi, bj
      INTEGER I, J, K
      _RL     pLoc, rhoLoc
      _RL     dPIdp
      CHARACTER*(MAX_LEN_MBUF) msgBuf
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     Calculate global domain area:
C     use to be in ini_masks_etc.F but has been move after packages_init_fixed
C     in case 1 pkg (e.g., OBCS) modifies the domain size.
      threadArea = 0. _d 0
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
C-      Compute the domain Area:
        tileArea(bi,bj) = 0. _d 0
        DO j=1,sNy
         DO i=1,sNx
          tileArea(bi,bj) = tileArea(bi,bj)
     &                    + rA(i,j,bi,bj)*maskInC(i,j,bi,bj)
c    &                    + rA(i,j,bi,bj)*maskH(i,j,bi,bj)
         ENDDO
        ENDDO
c       threadArea = threadArea + tileArea(bi,bj)
       ENDDO
      ENDDO
c#ifdef ALLOW_AUTODIFF_TAMC
C_jmc: apply GLOBAL_SUM to thread-local variable (not in common block)
c      _GLOBAL_SUM_RL( threadArea, myThid )
c#else
      CALL GLOBAL_SUM_TILE_RL( tileArea, threadArea, myThid )
c#endif
      _BEGIN_MASTER( myThid )
      globalArea = threadArea
C-    list empty tiles:
      msgBuf(1:1) = ' '
      DO bj = 1,nSy
       DO bi = 1,nSx
        IF ( tileArea(bi,bj).EQ.0. _d 0 ) THEN
#ifdef ALLOW_EXCH2
         WRITE(msgBuf,'(A,I6,A,2I4,A)')
     &    'Empty tile: #', W2_myTileList(bi,bj), ' (bi,bj=',bi,bj,' )'
#else
         WRITE(msgBuf,'(A,I6,I6)') 'Empty tile bi,bj=', bi, bj
#endif
         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                       SQUEEZE_RIGHT, myThid )
        ENDIF
       ENDDO
      ENDDO
      IF ( msgBuf(1:1).NE.' ' ) THEN
         WRITE(msgBuf,'(A)') ' '
         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                       SQUEEZE_RIGHT, myThid )
      ENDIF
      _END_MASTER( myThid )

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
             k = ksurfC(I,J,bi,bj)
             pLoc = Ro_surf(I,J,bi,bj)
             CALL FIND_RHO_SCALAR(
     I            tRef(k), sRef(k), pLoc,
     O            rhoLoc, myThid )
             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_RL(Bo_surf, myThid)
      _EXCH_XY_RL(recip_Bo, myThid)

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

        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)
       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	jmc	1.17	C $Header: /u/gcmpack/MITgcm/model/src/ini_linear_phisurf.F,v 1.16 2009/06/14 21:50:26 jmc Exp $
2	jmc	1.1	C $Name: $
3
4	jmc	1.17	#include "PACKAGES_CONFIG.h"
5	jmc	1.1	#include "CPP_OPTIONS.h"
6
7	cnh	1.3	CBOP
8			C !ROUTINE: INI_LINEAR_PHISURF
9			C !INTERFACE:
10	jmc	1.1	SUBROUTINE INI_LINEAR_PHISURF( myThid )
11	cnh	1.3
12			C !DESCRIPTION: \bv
13			C ==========================================================
14	jmc	1.13	C \| SUBROUTINE INI_LINEAR_PHISURF
15			C \| o Initialise the Linear Relation Phi_surf(eta)
16	cnh	1.3	C ==========================================================
17	jmc	1.17	C \| Compute global domain Area ;
18	jmc	1.13	C \| Initialise -Boyancy at surface level (Bo_surf)
19			C \| to setup the Linear relation: Phi_surf(eta)=Bo_surf*eta
20			C \| Initialise phi0surf = starting point for integrating
21	jmc	1.7	C \| phiHyd (= phiHyd at r=RoSurf)
22	cnh	1.3	C ==========================================================
23			C \ev
24
25			C !USES:
26	jmc	1.1	IMPLICIT NONE
27			C === Global variables ===
28			#include "SIZE.h"
29			#include "EEPARAMS.h"
30			#include "PARAMS.h"
31			#include "GRID.h"
32			#include "SURFACE.h"
33	jmc	1.17	#ifdef ALLOW_EXCH2
34			# include "W2_EXCH2_SIZE.h"
35			# include "W2_EXCH2_TOPOLOGY.h"
36			#endif /* ALLOW_EXCH2 */
37
38	jmc	1.1
39	cnh	1.3	C !INPUT/OUTPUT PARAMETERS:
40	jmc	1.1	C === Routine arguments ===
41			C myThid - Thread no. that called this routine.
42			INTEGER myThid
43
44	jmc	1.9	C == Local variables in common ==
45	jmc	1.17	C topoHloc had to be in common for multi threading but no longer
46			C needed since MDSIO now allows (2009/06/07) to write local arrays
47			_RL tileArea(nSx,nSy), threadArea
48			C put tileArea in (local) common block to print from master-thread:
49			COMMON / LOCAL_INI_PHISURF / tileArea
50	jmc	1.9
51	cnh	1.3	C !LOCAL VARIABLES:
52	jmc	1.1	C === Local variables ===
53	jmc	1.17	C topoHloc :: Temporary array used to write surface topography
54	jmc	1.16	C bi,bj :: tile indices
55			C I,J,K :: Loop counters
56			_RS topoHloc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
57	jmc	1.1	INTEGER bi, bj
58			INTEGER I, J, K
59	jmc	1.11	_RL pLoc, rhoLoc
60	jmc	1.1	_RL dPIdp
61	jmc	1.17	CHARACTER*(MAX_LEN_MBUF) msgBuf
62	cnh	1.3	CEOP
63	heimbach	1.6
64			#ifdef ALLOW_AUTODIFF_TAMC
65			DO bj=myByLo(myThid),myByHi(myThid)
66			DO bi=myBxLo(myThid),myBxHi(myThid)
67			DO J=1-Oly,sNy+Oly
68			DO I=1-Olx,sNx+Olx
69			Bo_surf(I,J,bi,bj) = 0. _d 0
70			recip_Bo(I,J,bi,bj) = 0. _d 0
71			ENDDO
72			ENDDO
73			ENDDO
74			ENDDO
75			#endif
76	jmc	1.1
77			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
78
79	jmc	1.17	C Calculate global domain area:
80			C use to be in ini_masks_etc.F but has been move after packages_init_fixed
81			C in case 1 pkg (e.g., OBCS) modifies the domain size.
82			threadArea = 0. _d 0
83			DO bj = myByLo(myThid), myByHi(myThid)
84			DO bi = myBxLo(myThid), myBxHi(myThid)
85			C- Compute the domain Area:
86			tileArea(bi,bj) = 0. _d 0
87			DO j=1,sNy
88			DO i=1,sNx
89			tileArea(bi,bj) = tileArea(bi,bj)
90			& + rA(i,j,bi,bj)*maskInC(i,j,bi,bj)
91			c & + rA(i,j,bi,bj)*maskH(i,j,bi,bj)
92			ENDDO
93			ENDDO
94			c threadArea = threadArea + tileArea(bi,bj)
95			ENDDO
96			ENDDO
97			c#ifdef ALLOW_AUTODIFF_TAMC
98			C_jmc: apply GLOBAL_SUM to thread-local variable (not in common block)
99			c _GLOBAL_SUM_RL( threadArea, myThid )
100			c#else
101			CALL GLOBAL_SUM_TILE_RL( tileArea, threadArea, myThid )
102			c#endif
103			_BEGIN_MASTER( myThid )
104			globalArea = threadArea
105			C- list empty tiles:
106			msgBuf(1:1) = ' '
107			DO bj = 1,nSy
108			DO bi = 1,nSx
109			IF ( tileArea(bi,bj).EQ.0. _d 0 ) THEN
110			#ifdef ALLOW_EXCH2
111			WRITE(msgBuf,'(A,I6,A,2I4,A)')
112			& 'Empty tile: #', W2_myTileList(bi,bj), ' (bi,bj=',bi,bj,' )'
113			#else
114			WRITE(msgBuf,'(A,I6,I6)') 'Empty tile bi,bj=', bi, bj
115			#endif
116			CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
117			& SQUEEZE_RIGHT, myThid )
118			ENDIF
119			ENDDO
120			ENDDO
121			IF ( msgBuf(1:1).NE.' ' ) THEN
122			WRITE(msgBuf,'(A)') ' '
123			CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
124			& SQUEEZE_RIGHT, myThid )
125			ENDIF
126			_END_MASTER( myThid )
127
128			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
129
130	jmc	1.1	C-- Initialise -Boyancy at surface level : Bo_surf
131	jmc	1.13	C Bo_surf is defined as d/dr(Phi_surf) and set to g/z2rUnit with
132			C z2rUnit = conversion factor from z-unit to r-unit: [z] * z2rUnit = [r]
133			C an accurate formulation includes P_surf and T,S_ref effects on rho_surf:
134	jmc	1.1	C (setting uniformLin_PhiSurf=.FALSE.):
135	jmc	1.13	C z-coord (z2rUnit=1): Bo_surf = - Boyancy
136			C = g * rho_surf(Tref,Sref,pSurf_0)/rho_0
137			C p-coord (z2rUnit=rho*g): Bo_surf = 1/rho(Tref(ksurf),pSurf_0)
138			C Note: on Phi_surf splitting : Non-linear Time-dependent effects on B_surf
139			C [through eta & (T-tRef)_surf] are included in PhiHyd rather than in Bo_surf
140	jmc	1.1	C--
141	jmc	1.14	IF ( buoyancyRelation .EQ. 'OCEANIC' ) THEN
142	jmc	1.1	C- gBaro = gravity (except for External mode test with reduced gravity)
143			DO bj=myByLo(myThid),myByHi(myThid)
144			DO bi=myBxLo(myThid),myBxHi(myThid)
145			DO J=1-Oly,sNy+Oly
146			DO I=1-Olx,sNx+Olx
147			Bo_surf(I,J,bi,bj) = gBaro
148			recip_Bo(I,J,bi,bj) = 1. _d 0 / gBaro
149			ENDDO
150			ENDDO
151			ENDDO
152			ENDDO
153			ELSEIF ( uniformLin_PhiSurf ) THEN
154			C- use a linear (in ps) uniform relation : Phi'_surf = 1/rhoConst * ps'_surf
155			DO bj=myByLo(myThid),myByHi(myThid)
156			DO bi=myBxLo(myThid),myBxHi(myThid)
157			DO J=1-Oly,sNy+Oly
158			DO I=1-Olx,sNx+Olx
159	jmc	1.13	c Bo_surf(I,J,bi,bj) = rVel2wUnit(1)*gravity
160			c recip_Bo(I,J,bi,bj) = wUnit2rVel(1)*recip_gravity
161			Bo_surf(I,J,bi,bj) = recip_rhoConst
162	jmc	1.1	recip_Bo(I,J,bi,bj) = rhoConst
163			ENDDO
164			ENDDO
165			ENDDO
166			ENDDO
167	jmc	1.14	ELSEIF ( buoyancyRelation .EQ. 'OCEANICP' ) THEN
168	mlosch	1.4	DO bj=myByLo(myThid),myByHi(myThid)
169			DO bi=myBxLo(myThid),myBxHi(myThid)
170			DO J=1-Oly,sNy+Oly
171			DO I=1-Olx,sNx+Olx
172	jmc	1.13	IF ( Ro_surf(I,J,bi,bj).GT.0. _d 0
173	mlosch	1.4	& .AND. ksurfC(I,J,bi,bj).LE.Nr ) THEN
174			k = ksurfC(I,J,bi,bj)
175	jmc	1.11	pLoc = Ro_surf(I,J,bi,bj)
176	jmc	1.13	CALL FIND_RHO_SCALAR(
177	jmc	1.11	I tRef(k), sRef(k), pLoc,
178			O rhoLoc, myThid )
179	jmc	1.13	IF ( rhoLoc .EQ. 0. _d 0 ) THEN
180	mlosch	1.4	Bo_surf(I,J,bi,bj) = 0. _d 0
181	jmc	1.13	ELSE
182			Bo_surf(I,J,bi,bj) = 1. _d 0/rhoLoc
183			ENDIF
184	mlosch	1.5	recip_Bo(I,J,bi,bj) = rhoLoc
185	mlosch	1.4	ELSE
186			Bo_surf(I,J,bi,bj) = 0. _d 0
187			recip_Bo(I,J,bi,bj) = 0. _d 0
188			ENDIF
189			ENDDO
190			ENDDO
191			ENDDO
192			ENDDO
193	jmc	1.14	ELSEIF ( buoyancyRelation .EQ. 'ATMOSPHERIC' ) THEN
194	jmc	1.1	C- use a linearized (in ps) Non-uniform relation : Bo_surf(Po_surf,tRef_surf)
195			C--- Bo = d/d_p(Phi_surf) = tRef_surfd/d_p(PI) ; PI = Cp(p/Po)^kappa
196			DO bj=myByLo(myThid),myByHi(myThid)
197			DO bi=myBxLo(myThid),myBxHi(myThid)
198			DO J=1-Oly,sNy+Oly
199			DO I=1-Olx,sNx+Olx
200	jmc	1.13	IF ( Ro_surf(I,J,bi,bj).GT.0. _d 0
201	jmc	1.2	& .AND. ksurfC(I,J,bi,bj).LE.Nr ) THEN
202	jmc	1.7	dPIdp = (atm_Cpatm_kappa/atm_Po)
203			& (Ro_surf(I,J,bi,bj)/atm_Po)**(atm_kappa-1. _d 0)
204	jmc	1.2	Bo_surf(I,J,bi,bj) = dPIdp*tRef(ksurfC(I,J,bi,bj))
205	jmc	1.1	recip_Bo(I,J,bi,bj) = 1. _d 0 / Bo_surf(I,J,bi,bj)
206			ELSE
207			Bo_surf(I,J,bi,bj) = 0.
208			recip_Bo(I,J,bi,bj) = 0.
209			ENDIF
210			ENDDO
211			ENDDO
212			ENDDO
213			ENDDO
214	mlosch	1.4	ELSE
215			STOP 'INI_LINEAR_PHISURF: We should never reach this point!'
216	jmc	1.1	ENDIF
217
218			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
219
220	jmc	1.13	C-- Update overlap regions (jmc: is it really needed ?)
221	jmc	1.15	_EXCH_XY_RL(Bo_surf, myThid)
222			_EXCH_XY_RL(recip_Bo, myThid)
223	jmc	1.1
224	jmc	1.14	IF ( ( buoyancyRelation .EQ. 'ATMOSPHERIC' .OR.
225			& buoyancyRelation .EQ. 'OCEANICP' )
226	mlosch	1.4	& .AND. .NOT.uniformLin_PhiSurf ) THEN
227	jmc	1.7
228	jmc	1.13	CALL WRITE_FLD_XY_RL( 'Bo_surf',' ',Bo_surf,0,myThid)
229	jmc	1.7
230	jmc	1.1	ENDIF
231
232	jmc	1.7	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
233
234	jmc	1.13	C-- Initialise phi0surf: used for atmos. surf. P-loading (ocean, z-coord)
235	jmc	1.7	C or topographic geopotential anom. (p-coord)
236
237			DO bj=myByLo(myThid),myByHi(myThid)
238			DO bi=myBxLo(myThid),myBxHi(myThid)
239			DO J=1-Oly,sNy+Oly
240			DO I=1-Olx,sNx+Olx
241			phi0surf(I,J,bi,bj) = 0.
242			ENDDO
243			ENDDO
244			ENDDO
245			ENDDO
246
247	jmc	1.14	IF ( buoyancyRelation .EQ. 'ATMOSPHERIC'
248	jmc	1.9	& .AND. topoFile.NE.' ' ) THEN
249	jmc	1.7
250	heimbach	1.8	#ifdef ALLOW_AUTODIFF_TAMC
251			STOP 'CANNOT PRESENTLY USE THIS OPTION WITH ADJOINT'
252			#else
253
254	jmc	1.13	C-- Compute topoH = PhiRef(Po_surf)/g ; is different from original
255	jmc	1.9	C topoZ(read from file) because of truncation of Po_surf.
256			C NOTE: not clear for now which topoZ needs to be saved in common block
257			C-- AND set phi0surf = starting point for integrating Geopotential;
258
259	jmc	1.13	CALL INI_P_GROUND( -2,
260			O topoHloc,
261	jmc	1.7	I Ro_surf, myThid )
262
263	jmc	1.9
264			IF (selectFindRoSurf.NE.0) THEN
265	jmc	1.7	_EXCH_XY_RS(phi0surf, myThid)
266	jmc	1.13	CALL WRITE_FLD_XY_RS( 'phi0surf',' ',phi0surf,0,myThid)
267	jmc	1.9	ENDIF
268	heimbach	1.8
269	jmc	1.12	CALL WRITE_FLD_XY_RS( 'topo_H',' ',topoHloc,0,myThid)
270
271	heimbach	1.8	#endif /* ALLOW_AUTODIFF_TAMC */
272	jmc	1.7
273			ENDIF
274
275			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
276	jmc	1.1	RETURN
277			END