forked from altMITgcm/MITgcm
-
Notifications
You must be signed in to change notification settings - Fork 0
/
seaice_calc_residual.F
134 lines (118 loc) · 4.48 KB
/
seaice_calc_residual.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
132
133
134
C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_calc_residual.F,v 1.7 2015/09/10 15:57:57 mlosch Exp $
C $Name: $
#include "SEAICE_OPTIONS.h"
CBOP
C !ROUTINE: SEAICE_CALC_RESIDUAL
C !INTERFACE:
SUBROUTINE SEAICE_CALC_RESIDUAL(
I uIceLoc, vIceLoc,
O uIceRes, vIceRes,
I newtonIter, krylovIter, myTime, myIter, myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | SUBROUTINE SEAICE_CALC_RESIDUAL
C | o For Jacobian-free Newton-Krylov solver compute
C | the residual of the momentum equations
C *==========================================================*
C | written by Martin Losch, Oct 2012
C *==========================================================*
C \ev
C !USES:
IMPLICIT NONE
C === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "SEAICE_SIZE.h"
#include "SEAICE.h"
C !INPUT/OUTPUT PARAMETERS:
C === Routine arguments ===
C myTime :: Simulation time
C myIter :: Simulation timestep number
C myThid :: my Thread Id. number
C newtonIter :: current iterate of Newton iteration
C krylovIter :: current iterate of Krylov iteration
_RL myTime
INTEGER myIter
INTEGER myThid
INTEGER newtonIter
INTEGER krylovIter
C u/vIceLoc :: local copies of the current ice velocity
_RL uIceLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
_RL vIceLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
C u/vIceRes :: residual of sea-ice momentum equations
_RL uIceRes(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
_RL vIceRes(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
#ifdef SEAICE_ALLOW_JFNK
C u/vIceLHS :: left hand side of momentum equations
_RL uIceLHS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
_RL vIceLHS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
C u/vIceRHS :: righ hand side of momentum equations
_RL uIceRHS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
_RL vIceRHS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
C i,j,bi,bj :: loop indices
INTEGER i,j,bi,bj
CEOP
C Initialise
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
DO J=1-Oly,sNy+Oly
DO I=1-Olx,sNx+Olx
uIceLHS(I,J,bi,bj) = 0. _d 0
vIceLHS(I,J,bi,bj) = 0. _d 0
uIceRHS(I,J,bi,bj) = 0. _d 0
vIceRHS(I,J,bi,bj) = 0. _d 0
ENDDO
ENDDO
ENDDO
ENDDO
C u/vIceLoc have changed so that new drag coefficients and
C viscosities are required
CALL SEAICE_OCEANDRAG_COEFFS(
I uIceLoc, vIceLoc,
O DWATN,
I krylovIter, myTime, myIter, myThid )
CALL SEAICE_CALC_STRAINRATES(
I uIceLoc, vIceLoc,
O e11, e22, e12,
I krylovIter, myTime, myIter, myThid )
CALL SEAICE_CALC_VISCOSITIES(
I e11, e22, e12, zMin, zMax, hEffM, press0, tensileStrength,
O eta, etaZ, zeta, zetaZ, press, deltaC,
I krylovIter, myTime, myIter, myThid )
C The scheme is backward Euler in time, i.e. the rhs-vector contains
C only terms that are independent of u/vIce, except for the time
C derivative part mass*(u/vIce-u/vIceNm1)/deltaT
C compute new right hand side (depends to DWATN=Cdrag)
C sea-surface tilt and wind stress: FORCEX0, FORCEY0
C + mass*(u/vIceNm1)/deltaT
C + Cdrag*(uVel*cosWat - vVel*sinWat)/(vVel*cosWat + uVel*sinWat)
CALL SEAICE_CALC_RHS(
O uIceRHS, vIceRHS,
I newtonIter, krylovIter, myTime, myIter, myThid )
C Left-hand side contributions:
C + mass*(u/vIce)/deltaT
C + Cdrag*(uIce*cosWat - vIce*sinWat)/(vIce*cosWat + uIce*sinWat)
C - mass*f*vIce/+mass*f*uIce
C - dsigma/dx / -dsigma/dy, eta and zeta are only computed once per
C Newton iterate
CALL SEAICE_CALC_LHS(
I uIceLoc, vIceLoc,
O uIceLHS, vIceLHS,
I newtonIter, myTime, myIter, myThid )
C Right-hand side contributions only need to be computed once per
C time step, therefore we will put them into a separate routine
C and call them elsewhere to save floating point operations
C Calculate the residual
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
DO J=1,sNy
DO I=1,sNx
uIceRes(I,J,bi,bj) = uIceLHS(I,J,bi,bj) - uIceRHS(I,J,bi,bj)
vIceRes(I,J,bi,bj) = vIceLHS(I,J,bi,bj) - vIceRHS(I,J,bi,bj)
ENDDO
ENDDO
ENDDO
ENDDO
#endif /* SEAICE_ALLOW_JFNK */
RETURN
END