model/src/ini_linear_phisurf.F

C $Header: /u/gcmpack/MITgcm/model/src/ini_linear_phisurf.F,v 1.10 2005/11/04 01:19:24 jmc 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/rtoz (linear free surface)
C     with rtoz = conversion factor from r-unit to z-unit  (=horiVertRatio)
C   an accurate formulation includes P_surf and T,S_surf effects on rho_surf:
C    (setting uniformLin_PhiSurf=.FALSE.):
C    z-ocean (rtoz=1) : Bo_surf = - Boyancy = gravity * rho_surf/rho_0 
C    p-atmos (rtoz=rho_c*g) : Bo_surf = (1/rho)_surf  
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
             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 )
             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./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
      _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

        _BEGIN_MASTER( myThid ) 
        CALL WRITE_FLD_XY_RL( 'Bo_surf',' ',Bo_surf,0,myThid) 
        _END_MASTER( 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 )

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

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

        _BEGIN_MASTER( myThid ) 
        CALL WRITE_FLD_XY_RS( 'phi0surf',' ',phi0surf,0,myThid) 
        _END_MASTER( myThid )
       ENDIF

#endif /* ALLOW_AUTODIFF_TAMC */

      ENDIF

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
      RETURN
      END
1	jmc	1.11	C $Header: /u/gcmpack/MITgcm/model/src/ini_linear_phisurf.F,v 1.10 2005/11/04 01:19:24 jmc Exp $
2	jmc	1.1	C $Name: $
3
4			#include "CPP_OPTIONS.h"
5
6	cnh	1.3	CBOP
7			C !ROUTINE: INI_LINEAR_PHISURF
8			C !INTERFACE:
9	jmc	1.1	SUBROUTINE INI_LINEAR_PHISURF( myThid )
10	cnh	1.3
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	jmc	1.7	C \| Initialise -Boyancy at surface level (Bo_surf)
17	cnh	1.3	C \| to setup the Linear relation: Phi_surf(eta)=Bo_surf*eta
18	jmc	1.7	C \| Initialise phi0surf = starting point for integrating
19			C \| phiHyd (= phiHyd at r=RoSurf)
20	cnh	1.3	C ==========================================================
21			C \ev
22
23			C !USES:
24	jmc	1.1	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	cnh	1.3	C !INPUT/OUTPUT PARAMETERS:
33	jmc	1.1	C === Routine arguments ===
34			C myThid - Thread no. that called this routine.
35			INTEGER myThid
36
37	jmc	1.9	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	cnh	1.3	C !LOCAL VARIABLES:
44	jmc	1.1	C === Local variables ===
45			C bi,bj - Loop counters
46			C I,J,K
47			INTEGER bi, bj
48			INTEGER I, J, K
49	jmc	1.11	_RL pLoc, rhoLoc
50	jmc	1.1	_RL dPIdp
51	cnh	1.3	CEOP
52	heimbach	1.6
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	jmc	1.1
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/rtoz (linear free surface)
70			C with rtoz = conversion factor from r-unit to z-unit (=horiVertRatio)
71			C an accurate formulation includes P_surf and T,S_surf effects on rho_surf:
72			C (setting uniformLin_PhiSurf=.FALSE.):
73			C z-ocean (rtoz=1) : Bo_surf = - Boyancy = gravity * rho_surf/rho_0
74			C p-atmos (rtoz=rho_c*g) : Bo_surf = (1/rho)_surf
75			C Note on Phi_surf splitting : Non-linear Time-dependent effects on b_surf
76			C [through eta & (T-tRef)_surf] are included in PhiHyd rather than in Bo_surf
77			C--
78			IF ( buoyancyRelation .eq. 'OCEANIC' ) THEN
79			C- gBaro = gravity (except for External mode test with reduced gravity)
80			DO bj=myByLo(myThid),myByHi(myThid)
81			DO bi=myBxLo(myThid),myBxHi(myThid)
82			DO J=1-Oly,sNy+Oly
83			DO I=1-Olx,sNx+Olx
84			Bo_surf(I,J,bi,bj) = gBaro
85			recip_Bo(I,J,bi,bj) = 1. _d 0 / gBaro
86			ENDDO
87			ENDDO
88			ENDDO
89			ENDDO
90			ELSEIF ( uniformLin_PhiSurf ) THEN
91			C- use a linear (in ps) uniform relation : Phi'_surf = 1/rhoConst * ps'_surf
92			DO bj=myByLo(myThid),myByHi(myThid)
93			DO bi=myBxLo(myThid),myBxHi(myThid)
94			DO J=1-Oly,sNy+Oly
95			DO I=1-Olx,sNx+Olx
96			Bo_surf(I,J,bi,bj) = recip_rhoConst
97			recip_Bo(I,J,bi,bj) = rhoConst
98			ENDDO
99			ENDDO
100			ENDDO
101			ENDDO
102	mlosch	1.4	ELSEIF ( buoyancyRelation .eq. 'OCEANICP' ) THEN
103			DO bj=myByLo(myThid),myByHi(myThid)
104			DO bi=myBxLo(myThid),myBxHi(myThid)
105			DO J=1-Oly,sNy+Oly
106			DO I=1-Olx,sNx+Olx
107			IF ( Ro_surf(I,J,bi,bj).GT.0. _d 0
108			& .AND. ksurfC(I,J,bi,bj).LE.Nr ) THEN
109			k = ksurfC(I,J,bi,bj)
110	jmc	1.11	pLoc = Ro_surf(I,J,bi,bj)
111	mlosch	1.4	CALL FIND_RHO_SCALAR(
112	jmc	1.11	I tRef(k), sRef(k), pLoc,
113			O rhoLoc, myThid )
114	mlosch	1.5	rhoLoc = rhoLoc + rhoConst
115			if ( rhoLoc .eq. 0. _d 0 ) then
116	mlosch	1.4	Bo_surf(I,J,bi,bj) = 0. _d 0
117			else
118	mlosch	1.5	Bo_surf(I,J,bi,bj) = 1./rhoLoc
119	mlosch	1.4	endif
120	mlosch	1.5	recip_Bo(I,J,bi,bj) = rhoLoc
121	mlosch	1.4	ELSE
122			Bo_surf(I,J,bi,bj) = 0. _d 0
123			recip_Bo(I,J,bi,bj) = 0. _d 0
124			ENDIF
125			ENDDO
126			ENDDO
127			ENDDO
128			ENDDO
129			ELSEIF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN
130	jmc	1.1	C- use a linearized (in ps) Non-uniform relation : Bo_surf(Po_surf,tRef_surf)
131			C--- Bo = d/d_p(Phi_surf) = tRef_surfd/d_p(PI) ; PI = Cp(p/Po)^kappa
132			DO bj=myByLo(myThid),myByHi(myThid)
133			DO bi=myBxLo(myThid),myBxHi(myThid)
134			DO J=1-Oly,sNy+Oly
135			DO I=1-Olx,sNx+Olx
136	jmc	1.2	IF ( Ro_surf(I,J,bi,bj).GT.0. _d 0
137			& .AND. ksurfC(I,J,bi,bj).LE.Nr ) THEN
138	jmc	1.7	dPIdp = (atm_Cpatm_kappa/atm_Po)
139			& (Ro_surf(I,J,bi,bj)/atm_Po)**(atm_kappa-1. _d 0)
140	jmc	1.2	Bo_surf(I,J,bi,bj) = dPIdp*tRef(ksurfC(I,J,bi,bj))
141	jmc	1.1	recip_Bo(I,J,bi,bj) = 1. _d 0 / Bo_surf(I,J,bi,bj)
142			ELSE
143			Bo_surf(I,J,bi,bj) = 0.
144			recip_Bo(I,J,bi,bj) = 0.
145			ENDIF
146			ENDDO
147			ENDDO
148			ENDDO
149			ENDDO
150	mlosch	1.4	ELSE
151			STOP 'INI_LINEAR_PHISURF: We should never reach this point!'
152	jmc	1.1	ENDIF
153
154			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
155
156			C-- Update overlap regions
157			_EXCH_XY_R8(Bo_surf, myThid)
158			_EXCH_XY_R8(recip_Bo, myThid)
159
160	mlosch	1.4	IF ( ( buoyancyRelation .eq. 'ATMOSPHERIC' .OR.
161			& buoyancyRelation .eq. 'OCEANICP' )
162			& .AND. .NOT.uniformLin_PhiSurf ) THEN
163	jmc	1.7
164			_BEGIN_MASTER( myThid )
165	jmc	1.1	CALL WRITE_FLD_XY_RL( 'Bo_surf',' ',Bo_surf,0,myThid)
166	jmc	1.7	_END_MASTER( myThid )
167
168	jmc	1.1	ENDIF
169
170	jmc	1.7	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
171
172			C-- Initialise phi0surf: used for atmos. surf. P-loading (ocean, z-coord)
173			C or topographic geopotential anom. (p-coord)
174
175			DO bj=myByLo(myThid),myByHi(myThid)
176			DO bi=myBxLo(myThid),myBxHi(myThid)
177			DO J=1-Oly,sNy+Oly
178			DO I=1-Olx,sNx+Olx
179			phi0surf(I,J,bi,bj) = 0.
180			ENDDO
181			ENDDO
182			ENDDO
183			ENDDO
184
185			IF ( buoyancyRelation .eq. 'ATMOSPHERIC'
186	jmc	1.9	& .AND. topoFile.NE.' ' ) THEN
187	jmc	1.7
188	heimbach	1.8	#ifdef ALLOW_AUTODIFF_TAMC
189			STOP 'CANNOT PRESENTLY USE THIS OPTION WITH ADJOINT'
190			#else
191
192	jmc	1.9	C-- Compute topoH = PhiRef(Po_surf)/g ; is different from original
193			C topoZ(read from file) because of truncation of Po_surf.
194			C NOTE: not clear for now which topoZ needs to be saved in common block
195			C-- AND set phi0surf = starting point for integrating Geopotential;
196
197			CALL INI_P_GROUND( -2,
198			O topoHloc,
199	jmc	1.7	I Ro_surf, myThid )
200
201	jmc	1.9	_BEGIN_MASTER( myThid )
202			CALL WRITE_FLD_XY_RS( 'topo_H',' ',topoHloc,0,myThid)
203			_END_MASTER( myThid )
204
205			IF (selectFindRoSurf.NE.0) THEN
206	jmc	1.7	_EXCH_XY_RS(phi0surf, myThid)
207
208			_BEGIN_MASTER( myThid )
209			CALL WRITE_FLD_XY_RS( 'phi0surf',' ',phi0surf,0,myThid)
210			_END_MASTER( myThid )
211	jmc	1.9	ENDIF
212	heimbach	1.8
213			#endif /* ALLOW_AUTODIFF_TAMC */
214	jmc	1.7
215			ENDIF
216
217			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
218	jmc	1.1	RETURN
219			END