model/src/ini_linear_phisurf.F

C $Header: /u/gcmpack/MITgcm/model/src/ini_linear_phisurf.F,v 1.6 2002/12/02 21:53:29 heimbach 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:
C     === Local variables ===
C     bi,bj  - Loop counters
C     I,J,K
      CHARACTER*(MAX_LEN_MBUF) msgBuf
      INTEGER bi, bj
      INTEGER I, J, K
      _RL     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)
             CALL FIND_RHO_SCALAR( 
     &            tRef(k), sRef(k), Ro_surf(I,J,bi,bj),
     &            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. selectFindRoSurf.EQ.1 ) THEN

C-    set phi0surf = starting point for integrating Geopotential:
        CALL INI_P_GROUND( -selectFindRoSurf, 
     O                     phi0surf, 
     I                     Ro_surf, myThid )

        _EXCH_XY_RS(phi0surf, myThid)

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

      ENDIF

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
      RETURN
      END
1	jmc	1.7	C $Header: /u/gcmpack/MITgcm/model/src/ini_linear_phisurf.F,v 1.6 2002/12/02 21:53:29 heimbach 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	cnh	1.3	C !LOCAL VARIABLES:
38	jmc	1.1	C === Local variables ===
39			C bi,bj - Loop counters
40			C I,J,K
41			CHARACTER*(MAX_LEN_MBUF) msgBuf
42			INTEGER bi, bj
43			INTEGER I, J, K
44	mlosch	1.4	_RL rhoLoc
45	jmc	1.1	_RL dPIdp
46	cnh	1.3	CEOP
47	heimbach	1.6
48			#ifdef ALLOW_AUTODIFF_TAMC
49			DO bj=myByLo(myThid),myByHi(myThid)
50			DO bi=myBxLo(myThid),myBxHi(myThid)
51			DO J=1-Oly,sNy+Oly
52			DO I=1-Olx,sNx+Olx
53			Bo_surf(I,J,bi,bj) = 0. _d 0
54			recip_Bo(I,J,bi,bj) = 0. _d 0
55			ENDDO
56			ENDDO
57			ENDDO
58			ENDDO
59			#endif
60	jmc	1.1
61			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
62
63			C-- Initialise -Boyancy at surface level : Bo_surf
64			C Bo_surf is defined as d/dr(Phi_surf) and set to g/rtoz (linear free surface)
65			C with rtoz = conversion factor from r-unit to z-unit (=horiVertRatio)
66			C an accurate formulation includes P_surf and T,S_surf effects on rho_surf:
67			C (setting uniformLin_PhiSurf=.FALSE.):
68			C z-ocean (rtoz=1) : Bo_surf = - Boyancy = gravity * rho_surf/rho_0
69			C p-atmos (rtoz=rho_c*g) : Bo_surf = (1/rho)_surf
70			C Note on Phi_surf splitting : Non-linear Time-dependent effects on b_surf
71			C [through eta & (T-tRef)_surf] are included in PhiHyd rather than in Bo_surf
72			C--
73			IF ( buoyancyRelation .eq. 'OCEANIC' ) THEN
74			C- gBaro = gravity (except for External mode test with reduced gravity)
75			DO bj=myByLo(myThid),myByHi(myThid)
76			DO bi=myBxLo(myThid),myBxHi(myThid)
77			DO J=1-Oly,sNy+Oly
78			DO I=1-Olx,sNx+Olx
79			Bo_surf(I,J,bi,bj) = gBaro
80			recip_Bo(I,J,bi,bj) = 1. _d 0 / gBaro
81			ENDDO
82			ENDDO
83			ENDDO
84			ENDDO
85			ELSEIF ( uniformLin_PhiSurf ) THEN
86			C- use a linear (in ps) uniform relation : Phi'_surf = 1/rhoConst * ps'_surf
87			DO bj=myByLo(myThid),myByHi(myThid)
88			DO bi=myBxLo(myThid),myBxHi(myThid)
89			DO J=1-Oly,sNy+Oly
90			DO I=1-Olx,sNx+Olx
91			Bo_surf(I,J,bi,bj) = recip_rhoConst
92			recip_Bo(I,J,bi,bj) = rhoConst
93			ENDDO
94			ENDDO
95			ENDDO
96			ENDDO
97	mlosch	1.4	ELSEIF ( buoyancyRelation .eq. 'OCEANICP' ) THEN
98			DO bj=myByLo(myThid),myByHi(myThid)
99			DO bi=myBxLo(myThid),myBxHi(myThid)
100			DO J=1-Oly,sNy+Oly
101			DO I=1-Olx,sNx+Olx
102			IF ( Ro_surf(I,J,bi,bj).GT.0. _d 0
103			& .AND. ksurfC(I,J,bi,bj).LE.Nr ) THEN
104			k = ksurfC(I,J,bi,bj)
105			CALL FIND_RHO_SCALAR(
106			& tRef(k), sRef(k), Ro_surf(I,J,bi,bj),
107			& rhoLoc, myThid )
108	mlosch	1.5	rhoLoc = rhoLoc + rhoConst
109			if ( rhoLoc .eq. 0. _d 0 ) then
110	mlosch	1.4	Bo_surf(I,J,bi,bj) = 0. _d 0
111			else
112	mlosch	1.5	Bo_surf(I,J,bi,bj) = 1./rhoLoc
113	mlosch	1.4	endif
114	mlosch	1.5	recip_Bo(I,J,bi,bj) = rhoLoc
115	mlosch	1.4	ELSE
116			Bo_surf(I,J,bi,bj) = 0. _d 0
117			recip_Bo(I,J,bi,bj) = 0. _d 0
118			ENDIF
119			ENDDO
120			ENDDO
121			ENDDO
122			ENDDO
123			ELSEIF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN
124	jmc	1.1	C- use a linearized (in ps) Non-uniform relation : Bo_surf(Po_surf,tRef_surf)
125			C--- Bo = d/d_p(Phi_surf) = tRef_surfd/d_p(PI) ; PI = Cp(p/Po)^kappa
126			DO bj=myByLo(myThid),myByHi(myThid)
127			DO bi=myBxLo(myThid),myBxHi(myThid)
128			DO J=1-Oly,sNy+Oly
129			DO I=1-Olx,sNx+Olx
130	jmc	1.2	IF ( Ro_surf(I,J,bi,bj).GT.0. _d 0
131			& .AND. ksurfC(I,J,bi,bj).LE.Nr ) THEN
132	jmc	1.7	dPIdp = (atm_Cpatm_kappa/atm_Po)
133			& (Ro_surf(I,J,bi,bj)/atm_Po)**(atm_kappa-1. _d 0)
134	jmc	1.2	Bo_surf(I,J,bi,bj) = dPIdp*tRef(ksurfC(I,J,bi,bj))
135	jmc	1.1	recip_Bo(I,J,bi,bj) = 1. _d 0 / Bo_surf(I,J,bi,bj)
136			ELSE
137			Bo_surf(I,J,bi,bj) = 0.
138			recip_Bo(I,J,bi,bj) = 0.
139			ENDIF
140			ENDDO
141			ENDDO
142			ENDDO
143			ENDDO
144	mlosch	1.4	ELSE
145			STOP 'INI_LINEAR_PHISURF: We should never reach this point!'
146	jmc	1.1	ENDIF
147
148			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
149
150			C-- Update overlap regions
151			_EXCH_XY_R8(Bo_surf, myThid)
152			_EXCH_XY_R8(recip_Bo, myThid)
153
154	mlosch	1.4	IF ( ( buoyancyRelation .eq. 'ATMOSPHERIC' .OR.
155			& buoyancyRelation .eq. 'OCEANICP' )
156			& .AND. .NOT.uniformLin_PhiSurf ) THEN
157	jmc	1.7
158			_BEGIN_MASTER( myThid )
159	jmc	1.1	CALL WRITE_FLD_XY_RL( 'Bo_surf',' ',Bo_surf,0,myThid)
160	jmc	1.7	_END_MASTER( myThid )
161
162	jmc	1.1	ENDIF
163
164	jmc	1.7	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
165
166			C-- Initialise phi0surf: used for atmos. surf. P-loading (ocean, z-coord)
167			C or topographic geopotential anom. (p-coord)
168
169			DO bj=myByLo(myThid),myByHi(myThid)
170			DO bi=myBxLo(myThid),myBxHi(myThid)
171			DO J=1-Oly,sNy+Oly
172			DO I=1-Olx,sNx+Olx
173			phi0surf(I,J,bi,bj) = 0.
174			ENDDO
175			ENDDO
176			ENDDO
177			ENDDO
178
179			IF ( buoyancyRelation .eq. 'ATMOSPHERIC'
180			& .AND. selectFindRoSurf.EQ.1 ) THEN
181
182			C- set phi0surf = starting point for integrating Geopotential:
183			CALL INI_P_GROUND( -selectFindRoSurf,
184			O phi0surf,
185			I Ro_surf, myThid )
186
187			_EXCH_XY_RS(phi0surf, myThid)
188
189			_BEGIN_MASTER( myThid )
190			CALL WRITE_FLD_XY_RS( 'phi0surf',' ',phi0surf,0,myThid)
191			_END_MASTER( myThid )
192
193			ENDIF
194
195			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
196	jmc	1.1	RETURN
197			END