/[MITgcm]/MITgcm_contrib/PRM/multi_comp_setup/cg/code/set_ddtvars.F
ViewVC logotype

Annotation of /MITgcm_contrib/PRM/multi_comp_setup/cg/code/set_ddtvars.F

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


Revision 1.5 - (hide annotations) (download)
Sat Mar 29 19:43:14 2008 UTC (17 years, 3 months ago) by jmc
Branch: MAIN
Changes since 1.4: +137 -30 lines
put mypackage momentum tencencies (myPa_TendVelU,V) on to C-grid ;
add option for smoothing momentum tendencies (using horizontal viscosity)

1 jmc 1.5 C $Header: /u/gcmpack/MITgcm_contrib/PRM/multi_comp_setup/cg/code/set_ddtvars.F,v 1.4 2008/02/25 19:01:55 jmc Exp $
2 jmc 1.1 C $Name: $
3    
4     #include "PACKAGES_CONFIG.h"
5 jmc 1.5 #ifdef ALLOW_MYPACKAGE
6 jmc 1.3 # include "MYPACKAGE_OPTIONS.h"
7     #else
8     # include "CPP_OPTIONS.h"
9     #endif
10 jmc 1.1
11     CBOP
12     C !ROUTINE: SET_DDTVARS
13     C !INTERFACE:
14     SUBROUTINE SET_DDTVARS( gUVelC, gVVelC, gThetaC,
15 jmc 1.5 I cnx, cny, cnr,
16 jmc 1.1 I myTime, myIter, myThid )
17     C !DESCRIPTION: \bv
18     C *==========================================================*
19     C | S/R SET_DDTVARS
20 jmc 1.5 C | o Set tendencies to be used as additional forcing
21 jmc 1.1 C *==========================================================*
22     C *==========================================================*
23     C \ev
24    
25     C !USES:
26     IMPLICIT NONE
27    
28     C == Global data ==
29     #include "SIZE.h"
30     #include "EEPARAMS.h"
31     #include "PARAMS.h"
32 jmc 1.5 #include "GRID.h"
33 jmc 1.1 #include "ORIENTATION.h"
34 jmc 1.5 #include "EXT_MOM_TEND.h"
35     #ifdef ALLOW_MYPACKAGE
36 jmc 1.3 # include "MYPACKAGE.h"
37     #endif
38 jmc 1.1
39     C !INPUT/OUTPUT PARAMETERS:
40     C == Routine arguments ==
41     C myTime :: Current time in simulation
42     C myIter :: Current iteration number
43     C myThid :: Thread Id number
44     INTEGER cnx, cny, cnr
45     REAL*8 gUVelC( cnx, cny, cnr)
46     REAL*8 gVVelC( cnx, cny, cnr)
47     REAL*8 gThetaC(cnx, cny, cnr)
48     _RL myTime
49     INTEGER myIter
50     INTEGER myThid
51    
52     C !FUNCTIONS:
53     LOGICAL DIFFERENT_MULTIPLE
54     EXTERNAL DIFFERENT_MULTIPLE
55    
56     C !LOCAL VARIABLES:
57     C == Local variables ==
58     C i,j,k :: Loop counters
59     INTEGER i,j,k
60 jmc 1.5 INTEGER bi,bj
61 jmc 1.1 C Variables used for I/O
62     CHARACTER*(10) suff
63 jmc 1.5 C Variables used for smoothing momentum tendencies:
64     _RL viscAhDivDt, viscAhVortDt
65     _RL locDiv (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
66     _RL locVort(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
67     _RL locFx (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
68     _RL locFy (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
69 jmc 1.1 CEOP
70    
71 jmc 1.5 viscAhDivDt = myPa_param1 * deltaTmom
72     viscAhVortDt = myPa_param2 * deltaTmom
73    
74 jmc 1.1 IF ( cnx .NE. sNx ) THEN
75     STOP 'SET_DTTVARS cnx NE sNx'
76     ENDIF
77     IF ( cny .NE. sNy ) THEN
78     STOP 'SET_DTTVARS cny NE sNy'
79     ENDIF
80     IF ( cnr .NE. Nr ) THEN
81     STOP 'SET_DTTVARS cnr NE Nr'
82     ENDIF
83    
84 jmc 1.5 DO bj = myByLo(myThid), myByHi(myThid)
85     DO bi = myBxLo(myThid), myBxHi(myThid)
86    
87     C--- Temp. Tendency from CG : just do a copy
88     DO k=1,Nr
89     DO j=1,sNy
90     DO i=1,sNx
91     myPa_TendScal1(i,j,k,bi,bj) = gThetaC(i,j,k)
92     ENDDO
93     ENDDO
94     ENDDO
95    
96     C--- Momentum Tendency from CG
97     C-- 1) rotate to get u,v aline with CG axes
98     DO k=1,Nr
99     DO j=1-Oly,sNy+Oly
100     DO i=1-Olx,sNx+Olx
101     ext_gu(i,j,k,bi,bj) = 0. _d 0
102     ext_gv(i,j,k,bi,bj) = 0. _d 0
103     ENDDO
104     ENDDO
105     DO j=1,sNy
106     DO i=1,sNx
107     IF ( velDvRMS(i,j,bi,bj).GT.0. _d 0 ) THEN
108     ext_gu(i,j,k,bi,bj) = gUVelC(i,j,k)*cAngleFG(i,j,bi,bj)
109     & - gVVelC(i,j,k)*sAngleFG(i,j,bi,bj)
110     ext_gv(i,j,k,bi,bj) = gVVelC(i,j,k)*cAngleFG(i,j,bi,bj)
111     & + gUVelC(i,j,k)*sAngleFG(i,j,bi,bj)
112     ENDIF
113     ENDDO
114     ENDDO
115 jmc 1.1 ENDDO
116 jmc 1.5
117     C- end bi,bj loops
118 jmc 1.1 ENDDO
119     ENDDO
120    
121 jmc 1.5 C-- 2) Fill the overlap : ext_gu & ext_gv are both on A-grid
122     CALL EXCH_UV_AGRID_3D_RL( ext_gu, ext_gv,
123     & .TRUE., Nr, myThid )
124    
125     DO bj = myByLo(myThid), myByHi(myThid)
126     DO bi = myBxLo(myThid), myBxHi(myThid)
127     DO k=1,Nr
128     C-- 3) average to C-grid
129     DO j=1-OLy,sNy+OLy
130     DO i=2-OLx,sNx+OLx
131     myPa_TendVelU(i,j,k,bi,bj) = maskW(i,j,k,bi,bj)
132     & *( ext_gu(i-1,j,k,bi,bj)
133     & +ext_gu( i ,j,k,bi,bj) )*0.5 _d 0
134     ENDDO
135     ENDDO
136     DO j=2-OLy,sNy+OLy
137     DO i=1-OLx,sNx+OLx
138     myPa_TendVelV(i,j,k,bi,bj) = maskS(i,j,k,bi,bj)
139     & *( ext_gv(i,j-1,k,bi,bj)
140     & +ext_gv(i, j ,k,bi,bj) )*0.5 _d 0
141     ENDDO
142     ENDDO
143     C-- 4) Smooth Momentum tendency :
144     C apply horizontal viscosity (*time-step = viscAhDivDt) on Divergence
145     C tendency and on vorticity tend. (viscosity*time-step = viscAhVortDt)
146     IF ( viscAhDivDt.GT.0. .OR. viscAhVortDt.GT.0. ) THEN
147     C Fluxes
148     DO j=2-OLy,sNy+OLy
149     DO i=2-OLx,sNx+OLx
150     locFx(i,j) = myPa_TendVelU(i,j,k,bi,bj)
151     & *dyG(i,j,bi,bj)*hFacW(i,j,k,bi,bj)
152     locFy(i,j) = myPa_TendVelV(i,j,k,bi,bj)
153     & *dxG(i,j,bi,bj)*hFacS(i,j,k,bi,bj)
154     ENDDO
155     ENDDO
156     C Divergence
157     DO j=2-OLy,sNy+OLy-1
158     DO i=2-OLx,sNx+OLx-1
159     locDiv(i,j) =
160     & ( ( locFx(i+1,j)-locFx(i,j) )
161     & +( locFy(i,j+1)-locFy(i,j) )
162     & )*recip_rA(i,j,bi,bj)*recip_hFacC(i,j,k,bi,bj)
163     ENDDO
164     ENDDO
165     C Vorticity
166     DO j=2-OLy,sNy+OLy
167     DO i=2-OLx,sNx+OLx
168     locVort(i,j) =
169     & ( ( myPa_TendVelV(i,j,k,bi,bj)*dyC(i,j,bi,bj)
170     & -myPa_TendVelV(i-1,j,k,bi,bj)*dyC(i-1,j,bi,bj) )
171     & -( myPa_TendVelU(i,j,k,bi,bj)*dxC(i,j,bi,bj)
172     & -myPa_TendVelU(i,j-1,k,bi,bj)*dxC(i,j-1,bi,bj) )
173     & )*recip_rAz(I,J,bi,bj)
174     ENDDO
175     ENDDO
176     C Apply to C-grid tendencies:
177     DO j=3-OLy,sNy+OLy-1
178     DO i=3-OLx,sNx+OLx-1
179     myPa_TendVelU(i,j,k,bi,bj) = myPa_TendVelU(i,j,k,bi,bj)
180     & + ( viscAhDivDt * ( locDiv(i,j) - locDiv(i-1,j) )
181     & * recip_dxC(i,j,bi,bj)
182     & - viscAhVortDt* ( locVort(i,j+1)-locVort(i,j) )
183     & * recip_dyG(i,j,bi,bj)
184     & )*maskW(i,j,k,bi,bj)
185     myPa_TendVelV(i,j,k,bi,bj) = myPa_TendVelV(i,j,k,bi,bj)
186     & + ( viscAhDivDt * ( locDiv(i,j) - locDiv(i,j-1) )
187     & * recip_dyC(i,j,bi,bj)
188     & + viscAhVortDt* ( locVort(i+1,j)-locVort(i,j) )
189     & * recip_dxG(i,j,bi,bj)
190     & )*maskS(i,j,k,bi,bj)
191     ENDDO
192     ENDDO
193     C- end smoothing & end of k loop
194     ENDIF
195 jmc 1.1 ENDDO
196 jmc 1.5
197     C- end bi,bj loops
198 jmc 1.1 ENDDO
199     ENDDO
200    
201 jmc 1.5 C-- Diagnostics
202 jmc 1.4 #ifdef ALLOW_DIAGNOSTICS
203     IF ( useDiagnostics ) THEN
204 jmc 1.5 CALL DIAGNOSTICS_FILL ( myPa_TendScal1,'MYPadTdt',
205     & 0, Nr, 0, 1, 1, myThid )
206 jmc 1.4 CALL DIAGNOSTICS_FILL ( myPa_TendVelU, 'MYPadUdt',
207     & 0, Nr, 0, 1, 1, myThid )
208     CALL DIAGNOSTICS_FILL ( myPa_TendVelV, 'MYPadVdt',
209     & 0, Nr, 0, 1, 1, myThid )
210 jmc 1.5 C- This is a hack (to avoid changing mypackage_diagnostics_init.F)
211     CALL DIAGNOSTICS_FILL ( ext_gu, 'MYPaStaU',
212     & 0, Nr, 0, 1, 1, myThid )
213     CALL DIAGNOSTICS_FILL ( ext_gv, 'MYPaStaV',
214 jmc 1.4 & 0, Nr, 0, 1, 1, myThid )
215     ENDIF
216     #endif /* ALLOW_DIAGNOSTICS */
217    
218 jmc 1.5 C-- Write snap-shot output:
219 jmc 1.1 IF ( DIFFERENT_MULTIPLE(dumpFreq,myTime,deltaTClock)
220     & .OR. dumpInitAndLast.AND.( myTime.EQ.endTime .OR.
221     & myTime.EQ.startTime )
222     & ) THEN
223     WRITE(suff,'(I10.10)') myIter
224 jmc 1.5 CALL WRITE_FLD_XYZ_RL( 'ext_gT.', suff, myPa_TendScal1,
225 jmc 1.3 & myIter, myThid )
226 jmc 1.5 CALL WRITE_FLD_XYZ_RL( 'ext_gU.', suff, ext_gu,
227 jmc 1.3 & myIter, myThid )
228 jmc 1.5 CALL WRITE_FLD_XYZ_RL( 'ext_gV.', suff, ext_gv,
229 jmc 1.3 & myIter, myThid )
230 jmc 1.1 ENDIF
231    
232     RETURN
233     END

  ViewVC Help
Powered by ViewVC 1.1.22