forked from darwinproject/darwin3
-
Notifications
You must be signed in to change notification settings - Fork 1
/
momentum_correction_step.F
131 lines (115 loc) · 3.7 KB
/
momentum_correction_step.F
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
#include "PACKAGES_CONFIG.h"
#include "CPP_OPTIONS.h"
CBOP
C !ROUTINE: MOMENTUM_CORRECTION_STEP
C !INTERFACE:
SUBROUTINE MOMENTUM_CORRECTION_STEP( myTime, myIter, myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | SUBROUTINE MOMENTUM_CORRECTION_STEP
C *==========================================================*
C |1rst Part : Update U,V.
C |
C | The arrays used for time stepping are cycled.
C | Momentum:
C | V(n) = Gv(n) - dt * grad Eta
C |
C |part1: update U,V
C | U*,V* (contained in gU,gV) have the surface
C | pressure gradient term added and the result stored
C | in U,V (contained in uVel, vVel)
C |
C |part2: Adjustments
C | o Filter U,V (Shapiro Filter, Zonal_Filter)
C *==========================================================*
C \ev
C !USES:
IMPLICIT NONE
C == Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "DYNVARS.h"
#ifdef ALLOW_SHAP_FILT
#include "SHAP_FILT.h"
#endif
#ifdef ALLOW_ZONAL_FILT
#include "ZONAL_FILT.h"
#endif
C !INPUT/OUTPUT PARAMETERS:
C myTime :: Current time in simulation
C myIter :: Current iteration number in simulation
C myThid :: my Thread Id. number
_RL myTime
INTEGER myIter
INTEGER myThid
C !LOCAL VARIABLES:
_RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
INTEGER iMin, iMax
INTEGER jMin, jMax
INTEGER bi, bj
INTEGER i, j
CEOP
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
IF ( momStepping ) THEN
C-- Set up work arrays that need valid initial values
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
phiSurfX(i,j) = 0.
phiSurfY(i,j) = 0.
ENDDO
ENDDO
C Loop range: Gradients of Eta are evaluated so valid
C range is all but first row and column in overlaps.
iMin = 1-OLx+1
iMax = sNx+OLx
jMin = 1-OLy+1
jMax = sNy+OLy
C- Calculate gradient of surface Potentiel
CALL CALC_GRAD_PHI_SURF(
I bi, bj, iMin, iMax, jMin, jMax,
I etaN,
O phiSurfX, phiSurfY,
I myThid )
C- Update velocity fields: V(n) = V** - dt * grad Eta
CALL CORRECTION_STEP(
I bi, bj, iMin, iMax, jMin, jMax,
I phiSurfX, phiSurfY,
I myTime, myIter, myThid )
ENDIF
#ifdef ALLOW_OBCS
IF ( useOBCS ) THEN
CALL OBCS_APPLY_UV( bi, bj, 0, uVel, vVel, myThid )
ENDIF
#endif /* ALLOW_OBCS */
C-- End of 1rst bi,bj loop
ENDDO
ENDDO
C--- 2nd Part : Adjustment.
C-- Filter (and exchange)
#ifdef ALLOW_SHAP_FILT
IF ( useSHAP_FILT ) THEN
IF ( .NOT.shap_filt_uvStar ) THEN
CALL TIMER_START('SHAP_FILT_UV [MOM_CORR_STEP]',myThid)
CALL SHAP_FILT_APPLY_UV( uVel, vVel, myTime, myIter, myThid )
CALL TIMER_STOP ('SHAP_FILT_UV [MOM_CORR_STEP]',myThid)
ENDIF
ENDIF
#endif
#ifdef ALLOW_ZONAL_FILT
IF ( useZONAL_FILT ) THEN
IF ( .NOT.zonal_filt_uvStar ) THEN
CALL TIMER_START('ZONAL_FILT_UV [MOM_CORR_STEP]',myThid)
CALL ZONAL_FILT_APPLY_UV( uVel, vVel, myThid )
CALL TIMER_STOP ('ZONAL_FILT_UV [MOM_CORR_STEP]',myThid)
ENDIF
ENDIF
#endif
C-- Try to fix restart (Pb at machine trucation level accumulating in wVel):
C Apply EXCH to horiz velocity before integrating continuity (-> wVel)
IF ( applyExchUV_early )
& CALL EXCH_UV_3D_RL( uVel, vVel, .TRUE., Nr, myThid )
RETURN
END