/[MITgcm]/MITgcm_contrib/heimbach/OpenAD/code_shallow_openad3/ini_linear_phisurf.F
ViewVC logotype

Annotation of /MITgcm_contrib/heimbach/OpenAD/code_shallow_openad3/ini_linear_phisurf.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.1 - (hide annotations) (download)
Sun Jan 14 22:42:54 2007 UTC (19 years, 9 months ago) by utke
Branch: MAIN
replace temporary for index expressions because it doesn't fit the checks for
the simple loop pattern

1 utke 1.1 C $Header: /u/gcmpack/MITgcm/model/src/ini_linear_phisurf.F,v 1.12 2006/10/17 18:52:34 jmc Exp $
2     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/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     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     pLoc = Ro_surf(I,J,bi,bj)
111     CALL FIND_RHO_SCALAR(
112     I tRef(ksurfC(I,J,bi,bj)),
113     &sRef(ksurfC(I,J,bi,bj)), pLoc,
114     O rhoLoc, myThid )
115     rhoLoc = rhoLoc + rhoConst
116     if ( rhoLoc .eq. 0. _d 0 ) then
117     Bo_surf(I,J,bi,bj) = 0. _d 0
118     else
119     Bo_surf(I,J,bi,bj) = 1./rhoLoc
120     endif
121     recip_Bo(I,J,bi,bj) = rhoLoc
122     ELSE
123     Bo_surf(I,J,bi,bj) = 0. _d 0
124     recip_Bo(I,J,bi,bj) = 0. _d 0
125     ENDIF
126     ENDDO
127     ENDDO
128     ENDDO
129     ENDDO
130     ELSEIF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN
131     C- use a linearized (in ps) Non-uniform relation : Bo_surf(Po_surf,tRef_surf)
132     C--- Bo = d/d_p(Phi_surf) = tRef_surf*d/d_p(PI) ; PI = Cp*(p/Po)^kappa
133     DO bj=myByLo(myThid),myByHi(myThid)
134     DO bi=myBxLo(myThid),myBxHi(myThid)
135     DO J=1-Oly,sNy+Oly
136     DO I=1-Olx,sNx+Olx
137     IF ( Ro_surf(I,J,bi,bj).GT.0. _d 0
138     & .AND. ksurfC(I,J,bi,bj).LE.Nr ) THEN
139     dPIdp = (atm_Cp*atm_kappa/atm_Po)*
140     & (Ro_surf(I,J,bi,bj)/atm_Po)**(atm_kappa-1. _d 0)
141     Bo_surf(I,J,bi,bj) = dPIdp*tRef(ksurfC(I,J,bi,bj))
142     recip_Bo(I,J,bi,bj) = 1. _d 0 / Bo_surf(I,J,bi,bj)
143     ELSE
144     Bo_surf(I,J,bi,bj) = 0.
145     recip_Bo(I,J,bi,bj) = 0.
146     ENDIF
147     ENDDO
148     ENDDO
149     ENDDO
150     ENDDO
151     ELSE
152     STOP 'INI_LINEAR_PHISURF: We should never reach this point!'
153     ENDIF
154    
155     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
156    
157     C-- Update overlap regions
158     _EXCH_XY_R8(Bo_surf, myThid)
159     _EXCH_XY_R8(recip_Bo, myThid)
160    
161     IF ( ( buoyancyRelation .eq. 'ATMOSPHERIC' .OR.
162     & buoyancyRelation .eq. 'OCEANICP' )
163     & .AND. .NOT.uniformLin_PhiSurf ) THEN
164    
165     CALL WRITE_FLD_XY_RL( 'Bo_surf',' ',Bo_surf,0,myThid)
166    
167     ENDIF
168    
169     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
170    
171     C-- Initialise phi0surf: used for atmos. surf. P-loading (ocean, z-coord)
172     C or topographic geopotential anom. (p-coord)
173    
174     DO bj=myByLo(myThid),myByHi(myThid)
175     DO bi=myBxLo(myThid),myBxHi(myThid)
176     DO J=1-Oly,sNy+Oly
177     DO I=1-Olx,sNx+Olx
178     phi0surf(I,J,bi,bj) = 0.
179     ENDDO
180     ENDDO
181     ENDDO
182     ENDDO
183    
184     IF ( buoyancyRelation .eq. 'ATMOSPHERIC'
185     & .AND. topoFile.NE.' ' ) THEN
186    
187     #ifdef ALLOW_AUTODIFF_TAMC
188     STOP 'CANNOT PRESENTLY USE THIS OPTION WITH ADJOINT'
189     #else
190    
191     C-- Compute topoH = PhiRef(Po_surf)/g ; is different from original
192     C topoZ(read from file) because of truncation of Po_surf.
193     C NOTE: not clear for now which topoZ needs to be saved in common block
194     C-- AND set phi0surf = starting point for integrating Geopotential;
195    
196     CALL INI_P_GROUND( -2,
197     O topoHloc,
198     I Ro_surf, myThid )
199    
200    
201     IF (selectFindRoSurf.NE.0) THEN
202     _EXCH_XY_RS(phi0surf, myThid)
203    
204     CALL WRITE_FLD_XY_RS( 'phi0surf',' ',phi0surf,0,myThid)
205     ELSE
206     _BARRIER
207     ENDIF
208    
209     CALL WRITE_FLD_XY_RS( 'topo_H',' ',topoHloc,0,myThid)
210    
211     #endif /* ALLOW_AUTODIFF_TAMC */
212    
213     ENDIF
214    
215     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
216     RETURN
217     END

  ViewVC Help
Powered by ViewVC 1.1.22