/[MITgcm]/MITgcm/model/src/calc_grad_phi_hyd.F
ViewVC logotype

Annotation of /MITgcm/model/src/calc_grad_phi_hyd.F

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

Revision 1.8 - (hide annotations) (download)
Thu Dec 8 15:44:33 2005 UTC (18 years, 5 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint58b_post, checkpoint57y_post, checkpoint58, checkpoint58f_post, checkpoint58d_post, checkpoint58a_post, checkpoint57z_post, checkpoint58m_post, checkpoint58o_post, checkpoint58p_post, checkpoint58q_post, checkpoint58e_post, checkpoint58r_post, checkpoint58n_post, checkpoint58k_post, checkpoint58l_post, checkpoint58g_post, checkpoint58h_post, checkpoint58j_post, checkpoint58i_post, checkpoint58c_post
Changes since 1.7: +3 -1 lines 
First step for a NLFS adjoint
o initially suppress rStar (new flag DISABLE_RSTAR_CODE)
o new init. routines for calc_r_star, calc_surf_dr
o still need to deal with ini_masks_etc
o testreport seemed happy
1 heimbach 1.8 C $Header: /u/gcmpack/MITgcm/model/src/calc_grad_phi_hyd.F,v 1.7 2005/11/05 01:00:57 jmc Exp $
2 jmc 1.1 C $Name: $
3    
4     #include "CPP_OPTIONS.h"
5    
6     CBOP
7     C !ROUTINE: CALC_GRAD_PHI_HYD
8     C !INTERFACE:
9     SUBROUTINE CALC_GRAD_PHI_HYD(
10     I k, bi, bj, iMin,iMax, jMin,jMax,
11 jmc 1.3 I phiHydC, alphRho, tFld, sFld,
12 jmc 1.1 O dPhiHydX, dPhiHydY,
13     I myTime, myIter, myThid)
14     C !DESCRIPTION: \bv
15     C *==========================================================*
16     C | S/R CALC_GRAD_PHI_HYD
17     C | o Calculate the gradient of Hydrostatic potential anomaly
18     C *==========================================================*
19     C \ev
20    
21     C !USES:
22     IMPLICIT NONE
23     C == Global variables ==
24     #include "SIZE.h"
25     #include "EEPARAMS.h"
26     #include "PARAMS.h"
27     #include "GRID.h"
28     #include "SURFACE.h"
29 jmc 1.3 #include "DYNVARS.h"
30 jmc 1.1
31     C !INPUT/OUTPUT PARAMETERS:
32     C == Routine Arguments ==
33     C bi,bj :: tile index
34     C iMin,iMax,jMin,jMax :: Loop counters
35 jmc 1.3 C phiHydC :: Hydrostatic Potential anomaly
36 jmc 1.1 C (atmos: =Geopotential ; ocean-z: =Pressure/rho)
37     C alphRho :: Density (z-coord) or specific volume (p-coord)
38     C tFld :: Potential temp.
39     C sFld :: Salinity
40     C dPhiHydX,Y :: Gradient (X & Y directions) of Hyd. Potential
41     C myTime :: Current time
42     C myIter :: Current iteration number
43     C myThid :: Instance number for this call of the routine.
44     INTEGER k, bi,bj, iMin,iMax, jMin,jMax
45 jmc 1.3 c _RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
46     _RL phiHydC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
47 jmc 1.1 _RL alphRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
48     _RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
49     _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
50     _RL dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
51     _RL dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
52     _RL myTime
53     INTEGER myIter, myThid
54    
55     #ifdef INCLUDE_PHIHYD_CALCULATION_CODE
56    
57     C !LOCAL VARIABLES:
58     C == Local variables ==
59     C i,j :: Loop counters
60     INTEGER i,j
61     _RL varLoc(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
62 jmc 1.7 #ifdef NONLIN_FRSURF
63 jmc 1.3 _RL factorZ, factorP, conv_theta2T
64 jmc 1.5 _RL factPI
65     CHARACTER*(MAX_LEN_MBUF) msgBuf
66 jmc 1.7 #endif
67 jmc 1.1 CEOP
68    
69 jmc 1.3 #ifdef NONLIN_FRSURF
70 jmc 1.5 IF (select_rStar.GE.2 .AND. nonlinFreeSurf.GE.4 ) THEN
71 heimbach 1.8 # ifndef DISABLE_RSTAR_CODE
72 jmc 1.5 C- Integral of b.dr = rStarFac * Integral of b.dr* :
73     C and will add later (select_rStar=2) the contribution of
74     C the slope of the r* coordinate.
75     IF ( buoyancyRelation .EQ. 'ATMOSPHERIC' ) THEN
76     C- Consistent with Phi'= Integr[ theta'.dPi ] :
77     DO j=jMin-1,jMax
78     DO i=iMin-1,iMax
79     varLoc(i,j) = phiHydC(i,j)*rStarFacC(i,j,bi,bj)**atm_kappa
80     & + phi0surf(i,j,bi,bj)
81     ENDDO
82     ENDDO
83     ELSE
84     DO j=jMin-1,jMax
85     DO i=iMin-1,iMax
86     varLoc(i,j) = phiHydC(i,j)*rStarFacC(i,j,bi,bj)
87     & + phi0surf(i,j,bi,bj)
88     ENDDO
89     ENDDO
90     ENDIF
91     ELSEIF (select_rStar.GE.1 .AND. nonlinFreeSurf.GE.4 ) THEN
92     C- Integral of b.dr but scaled to correspond to a fixed r-level (=r*)
93     C no contribution of the slope of the r* coordinate (select_rStar=1)
94     IF ( buoyancyRelation .EQ. 'ATMOSPHERIC' ) THEN
95     C- Consistent with Phi'= Integr[ theta'.dPi ] :
96     DO j=jMin-1,jMax
97     DO i=iMin-1,iMax
98     IF (Ro_surf(i,j,bi,bj).EQ.rC(k)) THEN
99     factPI=atm_Cp*( ((etaH(i,j,bi,bj)+rC(k))/atm_Po)**atm_kappa
100     & -( rC(k) /atm_Po)**atm_kappa
101     & )
102     varLoc(i,j) = factPI*alphRho(i,j)
103     ELSEIF (Ro_surf(i,j,bi,bj).NE.0. _d 0) THEN
104     factPI = (rC(k)/Ro_surf(i,j,bi,bj))**atm_kappa
105     varLoc(i,j) = phiHydC(i,j)
106     & *(rStarFacC(i,j,bi,bj)**atm_kappa - factPI)
107     & /(1. _d 0 -factPI)
108     & + phi0surf(i,j,bi,bj)
109     ENDIF
110     ENDDO
111     ENDDO
112     ELSE
113     DO j=jMin-1,jMax
114     DO i=iMin-1,iMax
115     IF (Ro_surf(i,j,bi,bj).EQ.rC(k)) THEN
116 jmc 1.6 WRITE(msgBuf,'(3A)') 'CALC_GRAD_PHI_HYD: ',
117     & 'Problem when Ro_surf=rC',
118     & ' with select_rStar,integr_GeoPot=1,4'
119 jmc 1.5 CALL PRINT_ERROR( msgBuf , myThid)
120     STOP 'CALC_GRAD_PHI_HYD: Pb in r* options implementation'
121     ELSE
122     varLoc(i,j) = phiHydC(i,j)
123     & *(etaH(i,j,bi,bj)+Ro_surf(i,j,bi,bj)-rC(k))
124     & / (Ro_surf(i,j,bi,bj)-rC(k))
125     & + phi0surf(i,j,bi,bj)
126     ENDIF
127     ENDDO
128 jmc 1.3 ENDDO
129 jmc 1.5 ENDIF
130 heimbach 1.8 # endif /* DISABLE_RSTAR_CODE */
131 jmc 1.3 ELSE
132     #else /* NONLIN_FRSURF */
133     IF (.TRUE.) THEN
134     #endif /* NONLIN_FRSURF */
135     DO j=jMin-1,jMax
136     DO i=iMin-1,iMax
137     varLoc(i,j) = phiHydC(i,j)+phi0surf(i,j,bi,bj)
138     ENDDO
139     ENDDO
140     ENDIF
141 jmc 1.1
142 jmc 1.3 C-- Zonal & Meridional gradient of potential anomaly
143 jmc 1.1 DO j=jMin,jMax
144     DO i=iMin,iMax
145     dPhiHydX(i,j) = _recip_dxC(i,j,bi,bj)
146     & *( varLoc(i,j)-varLoc(i-1,j) )
147     dPhiHydY(i,j) = _recip_dyC(i,j,bi,bj)
148     & *( varLoc(i,j)-varLoc(i,j-1) )
149     ENDDO
150     ENDDO
151 jmc 1.3
152     #ifdef NONLIN_FRSURF
153     IF (select_rStar.GE.2 .AND. nonlinFreeSurf.GE.1 ) THEN
154     IF ( buoyancyRelation .EQ. 'OCEANIC' ) THEN
155     C-- z* coordinate slope term: rho'/rho0 * Grad_r(g.z)
156     factorZ = gravity*recip_rhoConst*0.5 _d 0
157     DO j=jMin-1,jMax
158     DO i=iMin-1,iMax
159     varLoc(i,j) = etaH(i,j,bi,bj)
160     & *(1. _d 0 + rC(k)*recip_Rcol(i,j,bi,bj))
161     ENDDO
162     ENDDO
163     DO j=jMin,jMax
164     DO i=iMin,iMax
165     dPhiHydX(i,j) = dPhiHydX(i,j)
166     & +factorZ*(alphRho(i-1,j)+alphRho(i,j))
167     & *(varLoc(i,j)-varLoc(i-1,j))
168     & *recip_dxC(i,j,bi,bj)
169     dPhiHydY(i,j) = dPhiHydY(i,j)
170     & +factorZ*(alphRho(i,j-1)+alphRho(i,j))
171     & *(varLoc(i,j)-varLoc(i,j-1))
172     & *recip_dyC(i,j,bi,bj)
173     ENDDO
174     ENDDO
175     ELSEIF (buoyancyRelation .EQ. 'OCEANICP' ) THEN
176     C-- p* coordinate slope term: alpha' * Grad_r( p )
177     factorP = 0.5 _d 0
178     DO j=jMin,jMax
179     DO i=iMin,iMax
180     dPhiHydX(i,j) = dPhiHydX(i,j)
181     & +factorP*(alphRho(i-1,j)+alphRho(i,j))
182     & *(rStarFacC(i,j,bi,bj)-rStarFacC(i-1,j,bi,bj))
183     & *rC(k)*recip_dxC(i,j,bi,bj)
184     dPhiHydY(i,j) = dPhiHydY(i,j)
185     & +factorP*(alphRho(i,j-1)+alphRho(i,j))
186     & *(rStarFacC(i,j,bi,bj)-rStarFacC(i,j-1,bi,bj))
187     & *rC(k)*recip_dyC(i,j,bi,bj)
188     ENDDO
189     ENDDO
190     ELSEIF ( buoyancyRelation .EQ. 'ATMOSPHERIC' ) THEN
191     C-- p* coordinate slope term: alpha' * Grad_r( p )
192     conv_theta2T = (rC(k)/atm_Po)**atm_kappa
193     factorP = (atm_Rd/rC(k))*conv_theta2T*0.5 _d 0
194     DO j=jMin,jMax
195     DO i=iMin,iMax
196     dPhiHydX(i,j) = dPhiHydX(i,j)
197 jmc 1.5 & +factorP*(alphRho(i-1,j)+alphRho(i,j))
198 jmc 1.3 & *(rStarFacC(i,j,bi,bj)-rStarFacC(i-1,j,bi,bj))
199     & *rC(k)*recip_dxC(i,j,bi,bj)
200     dPhiHydY(i,j) = dPhiHydY(i,j)
201 jmc 1.5 & +factorP*(alphRho(i,j-1)+alphRho(i,j))
202 jmc 1.3 & *(rStarFacC(i,j,bi,bj)-rStarFacC(i,j-1,bi,bj))
203     & *rC(k)*recip_dyC(i,j,bi,bj)
204     ENDDO
205     ENDDO
206     ENDIF
207     ENDIF
208     #endif /* NONLIN_FRSURF */
209 jmc 1.1
210     #endif /* INCLUDE_PHIHYD_CALCULATION_CODE */
211    
212     RETURN
213     END
  ViewVC Help
Powered by ViewVC 1.1.22