1 |
C $Header: $ |
2 |
C $Name: $ |
3 |
|
4 |
#include "SEAICE_OPTIONS.h" |
5 |
|
6 |
CBOP |
7 |
C !ROUTINE: SEAICE_ADVDIFF |
8 |
|
9 |
C !INTERFACE: ========================================================== |
10 |
SUBROUTINE SEAICE_ADVDIFF( |
11 |
I myTime, myIter, myThid ) |
12 |
|
13 |
C !DESCRIPTION: \bv |
14 |
C /===========================================================\ |
15 |
C | SUBROUTINE SEAICE_ADVDIFF | |
16 |
C | o driver for different advection routines | |
17 |
C | calls an adaption of gad_advection to call different | |
18 |
C | advection routines of pkg/generic_advdiff | |
19 |
C \===========================================================/ |
20 |
IMPLICIT NONE |
21 |
c \ev |
22 |
|
23 |
C !USES: =============================================================== |
24 |
C === Global variables === |
25 |
#include "SIZE.h" |
26 |
#include "EEPARAMS.h" |
27 |
#include "PARAMS.h" |
28 |
#include "GRID.h" |
29 |
#include "GAD.h" |
30 |
#include "SEAICE_PARAMS.h" |
31 |
#include "SEAICE.h" |
32 |
|
33 |
#ifdef ALLOW_AUTODIFF_TAMC |
34 |
# include "tamc.h" |
35 |
#endif |
36 |
|
37 |
C !INPUT PARAMETERS: =================================================== |
38 |
C === Routine arguments === |
39 |
C UICE/VICE - ice velocity |
40 |
C HEFF - scalar field to be advected |
41 |
C HEFFM - mask for scalar field |
42 |
C myTime - current time |
43 |
C myIter - iteration number |
44 |
C myThid - Thread no. that called this routine. |
45 |
CML _RL UICE (1-OLx:sNx+OLx,1-OLy:sNy+OLy,3,nSx,nSy) |
46 |
CML _RL VICE (1-OLx:sNx+OLx,1-OLy:sNy+OLy,3,nSx,nSy) |
47 |
CML _RL HEFF (1-OLx:sNx+OLx,1-OLy:sNy+OLy,3,nSx,nSy) |
48 |
CML _RL HEFFM (1-OLx:sNx+OLx,1-OLy:sNy+OLy, nSx,nSy) |
49 |
_RL myTime |
50 |
INTEGER myIter |
51 |
INTEGER myThid |
52 |
CEndOfInterface |
53 |
|
54 |
#ifdef ALLOW_SEAICE |
55 |
C !LOCAL VARIABLES: ==================================================== |
56 |
C === Local variables === |
57 |
C i,j,bi,bj - Loop counters |
58 |
C uc/vc - current ice velocity on C-grid |
59 |
C fld - copy of scalar field |
60 |
C gfld - tendency of scalar field |
61 |
INTEGER i, j, bi, bj, k3 |
62 |
LOGICAL SEAICEmultiDimAdvection |
63 |
|
64 |
_RL uc (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) |
65 |
_RL vc (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) |
66 |
_RL fld (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) |
67 |
_RL gfld (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy) |
68 |
CEOP |
69 |
|
70 |
SEAICEmultidimadvection = .TRUE. |
71 |
IF ( SEAICEadvScheme.EQ.ENUM_CENTERED_2ND |
72 |
& .OR.SEAICEadvScheme.EQ.ENUM_UPWIND_3RD |
73 |
& .OR.SEAICEadvScheme.EQ.ENUM_CENTERED_4TH ) THEN |
74 |
SEAICEmultiDimAdvection = .FALSE. |
75 |
ENDIF |
76 |
|
77 |
|
78 |
IF ( SEAICEmultiDimAdvection ) THEN |
79 |
C This has to be done to comply with the time stepping in advect.F: |
80 |
C Making sure that the following routines see the different |
81 |
C time levels correctly |
82 |
C At the end of the routine ADVECT, |
83 |
C timelevel 1 is updated with advection contribution |
84 |
C and diffusion contribution |
85 |
C (which was computed in DIFFUS on timelevel 3) |
86 |
C timelevel 2 is the previous timelevel 1 |
87 |
C timelevel 3 is the total diffusion tendency * deltaT |
88 |
C (empty if no diffusion) |
89 |
|
90 |
DO bj=myByLo(myThid),myByHi(myThid) |
91 |
DO bi=myBxLo(myThid),myBxHi(myThid) |
92 |
DO j=1-OLy,sNy+OLy |
93 |
DO i=1-OLx,sNx+OLx |
94 |
HEFF(I,J,3,bi,bj) = 0. _d 0 !HEFF(I,J,2,bi,bj) |
95 |
HEFF(I,J,2,bi,bj) = HEFF(I,J,1,bi,bj) |
96 |
AREA(I,J,3,bi,bj) = 0. _d 0 !AREA(I,J,2,bi,bj) |
97 |
AREA(I,J,2,bi,bj) = AREA(I,J,1,bi,bj) |
98 |
ENDDO |
99 |
ENDDO |
100 |
ENDDO |
101 |
ENDDO |
102 |
|
103 |
DO bj=myByLo(myThid),myByHi(myThid) |
104 |
DO bi=myBxLo(myThid),myBxHi(myThid) |
105 |
C average seaice velocity to c-grid |
106 |
DO j=1-Oly,sNy+Oly-1 |
107 |
DO i=1-Olx,sNx+Olx-1 |
108 |
uc(I,J,bi,bj)=.5 _d 0*(UICE(I,J,1,bi,bj)+UICE(I,J+1,1,bi,bj)) |
109 |
vc(I,J,bi,bj)=.5 _d 0*(VICE(I,J,1,bi,bj)+VICE(I+1,J,1,bi,bj)) |
110 |
ENDDO |
111 |
ENDDO |
112 |
ENDDO |
113 |
ENDDO |
114 |
C Do we need this? I am afraid so. |
115 |
CALL EXCH_UV_XY_RL(uc,vc,.TRUE.,myThid) |
116 |
|
117 |
C Thickness (Volume) |
118 |
DO bj=myByLo(myThid),myByHi(myThid) |
119 |
DO bi=myBxLo(myThid),myBxHi(myThid) |
120 |
DO j=1-Oly,sNy+Oly |
121 |
DO i=1-Olx,sNx+Olx |
122 |
fld (I,J,bi,bj) = HEFF(I,J,1,bi,bj) |
123 |
gfld(I,J,bi,bj) = 0. _d 0 |
124 |
ENDDO |
125 |
ENDDO |
126 |
ENDDO |
127 |
ENDDO |
128 |
|
129 |
DO bj=myByLo(myThid),myByHi(myThid) |
130 |
DO bi=myBxLo(myThid),myBxHi(myThid) |
131 |
C advection |
132 |
CALL SEAICE_ADVECTION( |
133 |
I SEAICEadvScheme, GAD_HEFF, |
134 |
I uc, vc, fld, |
135 |
O gfld, |
136 |
I bi, bj, myTime, myIter, myThid) |
137 |
C now do the "explicit" time step |
138 |
DO j=1,sNy |
139 |
DO i=1,sNx |
140 |
HEFF(I,J,1,bi,bj) = HEFFM(I,J,bi,bj) * ( |
141 |
& HEFF(I,J,1,bi,bj) + SEAICE_deltaTtherm * gFld(I,J,bi,bj) |
142 |
& ) |
143 |
ENDDO |
144 |
ENDDO |
145 |
ENDDO |
146 |
ENDDO |
147 |
IF ( DIFF1 .GT. 0. _d 0 ) THEN |
148 |
C Do we need this? |
149 |
CALL SEAICE_EXCH( HEFF, myThid ) |
150 |
C diffusion |
151 |
CALL SEAICE_DIFFUSION( |
152 |
U HEFF, |
153 |
I HEFFM, SEAICE_deltaTtherm, myTime, myIter, myThid ) |
154 |
ENDIF |
155 |
|
156 |
C Fractional Area |
157 |
DO bj=myByLo(myThid),myByHi(myThid) |
158 |
DO bi=myBxLo(myThid),myBxHi(myThid) |
159 |
DO j=1-Oly,sNy+Oly |
160 |
DO i=1-Olx,sNx+Olx |
161 |
fld (I,J,bi,bj) = AREA(I,J,1,bi,bj) |
162 |
gfld(I,J,bi,bj) = 0. _d 0 |
163 |
ENDDO |
164 |
ENDDO |
165 |
ENDDO |
166 |
ENDDO |
167 |
|
168 |
DO bj=myByLo(myThid),myByHi(myThid) |
169 |
DO bi=myBxLo(myThid),myBxHi(myThid) |
170 |
C advection |
171 |
CALL SEAICE_ADVECTION( |
172 |
I SEAICEadvScheme, GAD_AREA, |
173 |
I uc, vc, fld, |
174 |
O gfld, |
175 |
I bi, bj, myTime, myIter, myThid) |
176 |
C now do the "explicit" time step |
177 |
DO j=1,sNy |
178 |
DO i=1,sNx |
179 |
AREA(I,J,1,bi,bj) = HEFFM(I,J,bi,bj) * ( |
180 |
& AREA(I,J,1,bi,bj) + SEAICE_deltaTtherm * gFld(I,J,bi,bj) |
181 |
& ) |
182 |
ENDDO |
183 |
ENDDO |
184 |
ENDDO |
185 |
ENDDO |
186 |
|
187 |
IF ( DIFF1 .GT. 0. _d 0 ) THEN |
188 |
C Do we need this? Probably not, but it is done in ADVECT |
189 |
CALL SEAICE_EXCH( AREA, myThid ) |
190 |
C diffusion |
191 |
CALL SEAICE_DIFFUSION( |
192 |
U AREA, |
193 |
I HEFFM, SEAICE_deltaTtherm, myTime, myIter, myThid ) |
194 |
ENDIF |
195 |
|
196 |
ELSE |
197 |
C if not multiDimAdvection |
198 |
CALL ADVECT( UICE, VICE, HEFF, HEFFM, myThid ) |
199 |
CALL ADVECT( UICE, VICE, AREA, HEFFM, myThid ) |
200 |
ENDIF |
201 |
|
202 |
#endif /* ALLOW_SEAICE */ |
203 |
|
204 |
RETURN |
205 |
END |