/[MITgcm]/MITgcm_contrib/dgoldberg/streamice/streamice_adv_front.F
ViewVC logotype

Annotation of /MITgcm_contrib/dgoldberg/streamice/streamice_adv_front.F

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


Revision 1.3 - (hide annotations) (download)
Thu May 3 15:52:06 2012 UTC (13 years, 2 months ago) by dgoldberg
Branch: MAIN
Changes since 1.2: +7 -5 lines
add stdout statement to indicate number of iterations taken

1 dgoldberg 1.3 C $Header: /u/gcmpack/MITgcm_contrib/dgoldberg/streamice/streamice_adv_front.F,v 1.2 2012/05/03 15:39:22 heimbach Exp $
2 heimbach 1.1 C $Name: $
3    
4     #include "STREAMICE_OPTIONS.h"
5    
6     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
7    
8     CBOP
9     SUBROUTINE STREAMICE_ADV_FRONT ( myThid, time_step )
10    
11     C /============================================================\
12     C | SUBROUTINE |
13     C | o |
14     C |============================================================|
15     C | |
16     C \============================================================/
17     IMPLICIT NONE
18    
19     C === Global variables ===
20     #include "SIZE.h"
21     #include "GRID.h"
22     #include "EEPARAMS.h"
23     #include "PARAMS.h"
24     #include "STREAMICE.h"
25     #include "STREAMICE_ADV.h"
26 heimbach 1.2 #ifdef ALLOW_AUTODIFF_TAMC
27     # include "tamc.h"
28     #endif
29 heimbach 1.1
30     INTEGER myThid
31     _RL time_step
32    
33     #ifdef ALLOW_STREAMICE
34    
35 dgoldberg 1.3 INTEGER i, j, bi, bj, k, iter_count, iter_rpt
36 heimbach 1.1 INTEGER Gi, Gj
37     INTEGER new_partial(4)
38 heimbach 1.2 INTEGER ikey_front, ikey_1
39     _RL iter_flag
40     _RL n_flux_1, n_flux_2
41 heimbach 1.1 _RL href, rho, partial_vol, tot_flux, hpot
42    
43     rho = streamice_density
44 heimbach 1.2 cph iter_count = 0
45     iter_flag = 1. _d 0
46 dgoldberg 1.3 iter_rpt = 0
47 heimbach 1.2
48     DO iter_count = 0, 3
49    
50     #ifdef ALLOW_AUTODIFF_TAMC
51     ikey_front = (ikey_dynamics-1)*4 + iter_count + 1
52     CADJ STORE area_shelf_streamice
53     CADJ & = comlev1_stream_front, key = ikey_front
54     CADJ STORE h_streamice
55     CADJ & = comlev1_stream_front, key = ikey_front
56     CADJ STORE hflux_x_si
57     CADJ & = comlev1_stream_front, key = ikey_front
58     CADJ STORE hflux_x_si2
59     CADJ & = comlev1_stream_front, key = ikey_front
60     CADJ STORE hflux_y_si
61     CADJ & = comlev1_stream_front, key = ikey_front
62     CADJ STORE hflux_y_si2
63     CADJ & = comlev1_stream_front, key = ikey_front
64     CADJ STORE streamice_hmask
65     CADJ & = comlev1_stream_front, key = ikey_front
66     CADJ STORE iter_flag
67     CADJ & = comlev1_stream_front, key = ikey_front
68     #endif
69 heimbach 1.1
70 heimbach 1.2 IF ( iter_flag .GT. 0. ) THEN
71    
72     iter_flag = 0. _d 0
73 heimbach 1.1
74     IF (iter_count .gt. 0) then
75     DO bj=myByLo(myThid),myByHi(myThid)
76     DO bi=myBxLo(myThid),myBxHi(myThid)
77     DO j=1-OLy,sNy+OLy
78     DO i=1-OLx,sNx+OLx
79     hflux_x_SI(i,j,bi,bj)=hflux_x_SI2(i,j,bi,bj)
80     hflux_y_SI(i,j,bi,bj)=hflux_y_SI2(i,j,bi,bj)
81     hflux_x_SI2(i,j,bi,bj) = 0. _d 0
82     hflux_y_SI2(i,j,bi,bj) = 0. _d 0
83     ENDDO
84     ENDDO
85     ENDDO
86     ENDDO
87     ENDIF
88    
89 heimbach 1.2 ! iter_count = iter_count + 1
90 dgoldberg 1.3 iter_rpt = iter_rpt + 1
91 heimbach 1.1
92     DO bj=myByLo(myThid),myByHi(myThid)
93     DO bi=myBxLo(myThid),myBxHi(myThid)
94 heimbach 1.2
95 heimbach 1.1 DO j=1-1,sNy+1
96     Gj = (myYGlobalLo-1)+(bj-1)*sNy+j
97 heimbach 1.2 cph IF ((Gj .ge. 1) .and. (Gj .le. Ny)) THEN
98 heimbach 1.1 DO i=1-1,sNx+1
99     Gi = (myXGlobalLo-1)+(bi-1)*sNx+i
100 heimbach 1.2
101     #ifdef ALLOW_AUTODIFF_TAMC
102     act1 = bi - myBxLo(myThid)
103     max1 = myBxHi(myThid) - myBxLo(myThid) + 1
104     act2 = bj - myByLo(myThid)
105     max2 = myByHi(myThid) - myByLo(myThid) + 1
106     act3 = myThid - 1
107     max3 = nTx*nTy
108     act4 = ikey_front - 1
109     ikey_1 = i
110     & + sNx*(j-1)
111     & + sNx*sNy*act1
112     & + sNx*sNy*max1*act2
113     & + sNx*sNy*max1*max2*act3
114     & + sNx*sNy*max1*max2*max3*act4
115     CADJ STORE area_shelf_streamice(i,j,bi,bj)
116     CADJ & = comlev1_stream_ij, key = ikey_1
117     CADJ STORE h_streamice(i,j,bi,bj)
118     CADJ & = comlev1_stream_ij, key = ikey_1
119     CADJ STORE hflux_x_si(i,j,bi,bj)
120     CADJ & = comlev1_stream_ij, key = ikey_1
121     CADJ STORE hflux_y_si(i,j,bi,bj)
122     CADJ & = comlev1_stream_ij, key = ikey_1
123     CADJ STORE streamice_hmask(i,j,bi,bj)
124     CADJ & = comlev1_stream_ij, key = ikey_1
125     #endif
126    
127     cph IF ((Gi .ge. 1) .and. (Gi .le. Nx) .and.
128     IF ((STREAMICE_Hmask(i,j,bi,bj).eq.0.0 .or.
129 heimbach 1.1 & STREAMICE_Hmask(i,j,bi,bj).eq.2.0)) THEN
130 heimbach 1.2 n_flux_1 = 0. _d 0
131 heimbach 1.1 href = 0. _d 0
132     tot_flux = 0. _d 0
133    
134 heimbach 1.2 #ifdef ALLOW_AUTODIFF_TAMC
135     CADJ STORE hflux_x_SI(i,j,bi,bj)
136     CADJ & = comlev1_stream_ij, key = ikey_1
137     #endif
138 heimbach 1.1 IF (hflux_x_SI(i,j,bi,bj).gt. 0. _d 0) THEN
139 heimbach 1.2 n_flux_1 = n_flux_1 + 1. _d 0
140 heimbach 1.1 href = href + H_streamice(i-1,j,bi,bj)
141     tot_flux = tot_flux + hflux_x_SI(i,j,bi,bj) *
142     & dxG(i,j,bi,bj) * time_step
143     hflux_x_SI(i,j,bi,bj) = 0. _d 0
144     ENDIF
145    
146 heimbach 1.2 #ifdef ALLOW_AUTODIFF_TAMC
147     CADJ STORE hflux_x_SI(i,j,bi,bj)
148     CADJ & = comlev1_stream_ij, key = ikey_1
149     #endif
150 heimbach 1.1 IF (hflux_x_SI(i+1,j,bi,bj).lt. 0. _d 0) THEN
151 heimbach 1.2 n_flux_1 = n_flux_1 + 1. _d 0
152 heimbach 1.1 href = href + H_streamice(i+1,j,bi,bj)
153     tot_flux = tot_flux - hflux_x_SI(i+1,j,bi,bj) *
154     & dxG(i+1,j,bi,bj) * time_step
155     hflux_x_SI(i+1,j,bi,bj) = 0. _d 0
156     ENDIF
157    
158 heimbach 1.2 #ifdef ALLOW_AUTODIFF_TAMC
159     CADJ STORE hflux_y_SI(i,j,bi,bj)
160     CADJ & = comlev1_stream_ij, key = ikey_1
161     #endif
162 heimbach 1.1 IF (hflux_y_SI(i,j,bi,bj).gt. 0. _d 0) THEN
163 heimbach 1.2 n_flux_1 = n_flux_1 + 1. _d 0
164 heimbach 1.1 href = href + H_streamice(i,j-1,bi,bj)
165     tot_flux = tot_flux + hflux_y_SI(i,j,bi,bj) *
166     & dyG(i,j,bi,bj) * time_step
167     hflux_y_SI(i,j,bi,bj) = 0. _d 0
168     ENDIF
169    
170 heimbach 1.2 #ifdef ALLOW_AUTODIFF_TAMC
171     CADJ STORE hflux_y_SI(i,j,bi,bj)
172     CADJ & = comlev1_stream_ij, key = ikey_1
173     #endif
174 heimbach 1.1 IF (hflux_y_SI(i,j+1,bi,bj).lt. 0. _d 0) THEN
175 heimbach 1.2 n_flux_1 = n_flux_1 + 1. _d 0
176 heimbach 1.1 href = href + H_streamice(i,j+1,bi,bj)
177     tot_flux = tot_flux - hflux_y_SI(i,j+1,bi,bj) *
178     & dyG(i,j+1,bi,bj) * time_step
179     hflux_y_SI(i,j+1,bi,bj) = 0. _d 0
180     ENDIF
181    
182 heimbach 1.2 IF (n_flux_1 .gt. 0.) THEN
183 heimbach 1.1
184 heimbach 1.2 href = href / n_flux_1
185 heimbach 1.1 partial_vol = H_streamice (i,j,bi,bj) *
186     & area_shelf_streamice (i,j,bi,bj) + tot_flux
187     hpot = partial_vol * recip_rA(i,j,bi,bj)
188    
189     IF (hpot .eq. href) THEN ! cell is exactly covered, no overflow
190     STREAMICE_hmask (i,j,bi,bj) = 1.0
191     H_streamice (i,j,bi,bj) = href
192     area_shelf_streamice(i,j,bi,bj) =
193     & rA(i,j,bi,bj)
194     ELSEIF (hpot .lt. href) THEN ! cell still unfilled
195    
196     ! PRINT *, "PARTIAL CELL INCREASED", tot_flux, i,
197     ! & area_shelf_streamice (i,j,bi,bj),
198     ! & H_streamice (i,j,bi,bj)
199    
200     STREAMICE_hmask (i,j,bi,bj) = 2.0
201     area_shelf_streamice (i,j,bi,bj) = partial_vol / href
202     H_streamice (i,j,bi,bj) = href
203     ELSE ! cell is filled - do overflow
204    
205     ! PRINT *, "CELL FILLED"
206    
207     STREAMICE_hmask (i,j,bi,bj) = 1.0
208     area_shelf_streamice(i,j,bi,bj) =
209     & rA(i,j,bi,bj)
210    
211    
212     partial_vol = partial_vol - href * rA(i,j,bi,bj)
213    
214 heimbach 1.2 iter_flag = 1. _d 0
215 heimbach 1.1
216 heimbach 1.2 n_flux_2 = 0. _d 0 ;
217 heimbach 1.1 DO k=1,4
218     new_partial (:) = 0
219     ENDDO
220    
221     DO k=1,2
222     IF (STREAMICE_ufacemask(i-1+k,j,bi,bj).eq.2.0) THEN ! at a permanent calving boundary - no advance allowed
223 heimbach 1.2 n_flux_2 = n_flux_2 + 1. _d 0
224 heimbach 1.1 ELSEIF (STREAMICE_hmask(i+2*k-3,j,bi,bj).eq.0 _d 0) THEN ! adjacent cell is completely ice free
225 heimbach 1.2 n_flux_2 = n_flux_2 + 1. _d 0
226 heimbach 1.1 new_partial (k) = 1
227     ENDIF
228     ENDDO
229     DO k=1,2
230     IF (STREAMICE_vfacemask (i,j-1+k,bi,bj).eq.2.0) THEN
231 heimbach 1.2 n_flux_2 = n_flux_2 + 1. _d 0
232 heimbach 1.1 ELSEIF (STREAMICE_hmask(i,j+2*k-3,bi,bj).eq.0 _d 0) THEN
233 heimbach 1.2 n_flux_2 = n_flux_2 + 1. _d 0
234 heimbach 1.1 new_partial (k+2) = 1
235     ENDIF
236     ENDDO
237    
238 heimbach 1.2 IF (n_flux_2 .eq. 0.) THEN ! there is nowhere to put the extra ice!
239 heimbach 1.1 H_streamice(i,j,bi,bj) = href + partial_vol *
240     & recip_rA(i,j,bi,bj)
241     ELSE
242     H_streamice(i,j,bi,bj) = href
243    
244     DO k=1,2
245     IF (new_partial(k) .eq. 1) THEN
246     hflux_x_SI2(i-1+k,j,bi,bj) =
247 heimbach 1.2 & partial_vol/time_step/n_flux_2/
248 heimbach 1.1 & dxG(i-1+k,j,bi,bj)
249     ENDIF
250     ENDDO
251    
252     DO k=1,2
253     IF (new_partial(k+2) .eq. 1) THEN
254     hflux_y_SI2(i,j-1+k,bi,bj) =
255 heimbach 1.2 & partial_vol/time_step/n_flux_2/
256 heimbach 1.1 & dxG(i,j-1+k,bi,bj)
257     ENDIF
258     ENDDO
259    
260     ENDIF
261     ENDIF
262     ENDIF
263 heimbach 1.2
264 heimbach 1.1 ENDIF
265     ENDDO
266 heimbach 1.2 cph ENDIF
267 heimbach 1.1 ENDDO
268 heimbach 1.2 c
269 heimbach 1.1 ENDDO
270     ENDDO
271 heimbach 1.2 c
272     ENDIF
273 heimbach 1.1 ENDDO
274    
275 dgoldberg 1.3 IF (iter_rpt.gt.1) THEN
276     PRINT *, "FRONT ADVANCE: ", iter_rpt, " ITERATIONS"
277     ENDIF
278 heimbach 1.1
279    
280    
281     #endif
282     RETURN
283     END

  ViewVC Help
Powered by ViewVC 1.1.22