model/src/ini_linear_phisurf.F

C $Header: /u/gcmpack/MITgcm/model/src/ini_linear_phisurf.F,v 1.5 2002/09/25 19:36:50 mlosch 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     | Presently: Initialise -Boyancy at surface level (Bo_surf) 
C     |  to setup the Linear relation: Phi_surf(eta)=Bo_surf*eta  
C     | Futur: might add other things for Non-Linear FreeSurface  
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
        CALL WRITE_FLD_XY_RL( 'Bo_surf',' ',Bo_surf,0,myThid) 
      ENDIF

      RETURN
      END
1	heimbach	1.6	C $Header: /u/gcmpack/MITgcm/model/src/ini_linear_phisurf.F,v 1.5 2002/09/25 19:36:50 mlosch 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			C \| Presently: Initialise -Boyancy at surface level (Bo_surf)
17			C \| to setup the Linear relation: Phi_surf(eta)=Bo_surf*eta
18			C \| Futur: might add other things for Non-Linear FreeSurface
19			C ==========================================================
20			C \ev
21
22			C !USES:
23	jmc	1.1	IMPLICIT NONE
24			C === Global variables ===
25			#include "SIZE.h"
26			#include "EEPARAMS.h"
27			#include "PARAMS.h"
28			#include "GRID.h"
29			#include "SURFACE.h"
30
31	cnh	1.3	C !INPUT/OUTPUT PARAMETERS:
32	jmc	1.1	C === Routine arguments ===
33			C myThid - Thread no. that called this routine.
34			INTEGER myThid
35
36	cnh	1.3	C !LOCAL VARIABLES:
37	jmc	1.1	C === Local variables ===
38			C bi,bj - Loop counters
39			C I,J,K
40			CHARACTER*(MAX_LEN_MBUF) msgBuf
41			INTEGER bi, bj
42			INTEGER I, J, K
43	mlosch	1.4	_RL rhoLoc
44	jmc	1.1	_RL dPIdp
45	cnh	1.3	CEOP
46	heimbach	1.6
47			#ifdef ALLOW_AUTODIFF_TAMC
48			DO bj=myByLo(myThid),myByHi(myThid)
49			DO bi=myBxLo(myThid),myBxHi(myThid)
50			DO J=1-Oly,sNy+Oly
51			DO I=1-Olx,sNx+Olx
52			Bo_surf(I,J,bi,bj) = 0. _d 0
53			recip_Bo(I,J,bi,bj) = 0. _d 0
54			ENDDO
55			ENDDO
56			ENDDO
57			ENDDO
58			#endif
59	jmc	1.1
60			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
61
62			C-- Initialise -Boyancy at surface level : Bo_surf
63			C Bo_surf is defined as d/dr(Phi_surf) and set to g/rtoz (linear free surface)
64			C with rtoz = conversion factor from r-unit to z-unit (=horiVertRatio)
65			C an accurate formulation includes P_surf and T,S_surf effects on rho_surf:
66			C (setting uniformLin_PhiSurf=.FALSE.):
67			C z-ocean (rtoz=1) : Bo_surf = - Boyancy = gravity * rho_surf/rho_0
68			C p-atmos (rtoz=rho_c*g) : Bo_surf = (1/rho)_surf
69			C Note on Phi_surf splitting : Non-linear Time-dependent effects on b_surf
70			C [through eta & (T-tRef)_surf] are included in PhiHyd rather than in Bo_surf
71			C--
72			IF ( buoyancyRelation .eq. 'OCEANIC' ) THEN
73			C- gBaro = gravity (except for External mode test with reduced gravity)
74			DO bj=myByLo(myThid),myByHi(myThid)
75			DO bi=myBxLo(myThid),myBxHi(myThid)
76			DO J=1-Oly,sNy+Oly
77			DO I=1-Olx,sNx+Olx
78			Bo_surf(I,J,bi,bj) = gBaro
79			recip_Bo(I,J,bi,bj) = 1. _d 0 / gBaro
80			ENDDO
81			ENDDO
82			ENDDO
83			ENDDO
84			ELSEIF ( uniformLin_PhiSurf ) THEN
85			C- use a linear (in ps) uniform relation : Phi'_surf = 1/rhoConst * ps'_surf
86			DO bj=myByLo(myThid),myByHi(myThid)
87			DO bi=myBxLo(myThid),myBxHi(myThid)
88			DO J=1-Oly,sNy+Oly
89			DO I=1-Olx,sNx+Olx
90			Bo_surf(I,J,bi,bj) = recip_rhoConst
91			recip_Bo(I,J,bi,bj) = rhoConst
92			ENDDO
93			ENDDO
94			ENDDO
95			ENDDO
96	mlosch	1.4	ELSEIF ( buoyancyRelation .eq. 'OCEANICP' ) THEN
97			DO bj=myByLo(myThid),myByHi(myThid)
98			DO bi=myBxLo(myThid),myBxHi(myThid)
99			DO J=1-Oly,sNy+Oly
100			DO I=1-Olx,sNx+Olx
101			IF ( Ro_surf(I,J,bi,bj).GT.0. _d 0
102			& .AND. ksurfC(I,J,bi,bj).LE.Nr ) THEN
103			k = ksurfC(I,J,bi,bj)
104			CALL FIND_RHO_SCALAR(
105			& tRef(k), sRef(k), Ro_surf(I,J,bi,bj),
106			& rhoLoc, myThid )
107	mlosch	1.5	rhoLoc = rhoLoc + rhoConst
108			if ( rhoLoc .eq. 0. _d 0 ) then
109	mlosch	1.4	Bo_surf(I,J,bi,bj) = 0. _d 0
110			else
111	mlosch	1.5	Bo_surf(I,J,bi,bj) = 1./rhoLoc
112	mlosch	1.4	endif
113	mlosch	1.5	recip_Bo(I,J,bi,bj) = rhoLoc
114	mlosch	1.4	ELSE
115			Bo_surf(I,J,bi,bj) = 0. _d 0
116			recip_Bo(I,J,bi,bj) = 0. _d 0
117			ENDIF
118			ENDDO
119			ENDDO
120			ENDDO
121			ENDDO
122			ELSEIF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN
123	jmc	1.1	C- use a linearized (in ps) Non-uniform relation : Bo_surf(Po_surf,tRef_surf)
124			C--- Bo = d/d_p(Phi_surf) = tRef_surfd/d_p(PI) ; PI = Cp(p/Po)^kappa
125			DO bj=myByLo(myThid),myByHi(myThid)
126			DO bi=myBxLo(myThid),myBxHi(myThid)
127			DO J=1-Oly,sNy+Oly
128			DO I=1-Olx,sNx+Olx
129	jmc	1.2	IF ( Ro_surf(I,J,bi,bj).GT.0. _d 0
130			& .AND. ksurfC(I,J,bi,bj).LE.Nr ) THEN
131	jmc	1.1	dPIdp = (atm_cpatm_kappa/atm_po)
132			& (Ro_surf(I,J,bi,bj)/atm_po)**(atm_kappa-1. _d 0)
133	jmc	1.2	Bo_surf(I,J,bi,bj) = dPIdp*tRef(ksurfC(I,J,bi,bj))
134	jmc	1.1	recip_Bo(I,J,bi,bj) = 1. _d 0 / Bo_surf(I,J,bi,bj)
135			ELSE
136			Bo_surf(I,J,bi,bj) = 0.
137			recip_Bo(I,J,bi,bj) = 0.
138			ENDIF
139			ENDDO
140			ENDDO
141			ENDDO
142			ENDDO
143	mlosch	1.4	ELSE
144			STOP 'INI_LINEAR_PHISURF: We should never reach this point!'
145	jmc	1.1	ENDIF
146
147			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
148
149			C-- Update overlap regions
150			_EXCH_XY_R8(Bo_surf, myThid)
151			_EXCH_XY_R8(recip_Bo, myThid)
152
153	mlosch	1.4	IF ( ( buoyancyRelation .eq. 'ATMOSPHERIC' .OR.
154			& buoyancyRelation .eq. 'OCEANICP' )
155			& .AND. .NOT.uniformLin_PhiSurf ) THEN
156	jmc	1.1	CALL WRITE_FLD_XY_RL( 'Bo_surf',' ',Bo_surf,0,myThid)
157			ENDIF
158
159			RETURN
160			END