source: branches/iLoveclim/SOURCES/calving_frange.f90 @ 90

Last change on this file since 90 was 90, checked in by aquiquet, 8 years ago

GRISLI-loveclim branch, GRISLI freshwater flux needed by iLOVECLIM
File size: 18.9 KB
Line 
1!> \file calving_frange.f90
2!! Module et routines qui calculent le calving avec la methode de Vincent
3!<
4
5
6!> \namespace module3D_phy
7!! Module et routines qui calculent le calving avec la methode de Vincent
8!! \author ...
9!! \date ...
10!! @note recopie direct de icethick5-ant. A revoir en particulier, accroche avec la glace posee
11!! @note Used modules
12!! @note   - use module3D_phy
13!! @todo il faudrait que le calving soit inactif dans les zones ou l epaisseur est imposee
14!<
15
16module calving_frange
17
18use module3D_phy
19implicit none
20
21real, dimension (nx,ny) :: hmhc  ! hauteur au dessus de la coupure
22real, dimension (nx,ny) :: bil_tot ! bilan surface et fond (bm-bmelt)
23real :: hcoup ! epaiseur de coupure au temps time
24real :: hcoup_max  ! epaisseur max
25real :: hcoup_min=100.  ! epaisseur min
26real :: mass_total
27integer :: meth_hcoup  ! pour avoir hcoup dépendant du coefbmshelf
28integer :: ifrange        !  0 pas de traitement particulier pres du bord, 1 -> franges
29integer ::  iin2,jin2,iin3,jin3             ! pour la detection polynies
30logical :: testmij,testpij,testimj,testipj
31logical :: bilan_surf_fond                  ! vrai si bm-bmelt est positif
32logical :: avalw,avale,avals,avaln,interieur
33
34integer :: testH ! somme des flux entrants en un point ij
35
36contains
37
38!---------------------------------------------------------------------------------------
39subroutine init_calving
40
41namelist/calving/Hcoup,ifrange,meth_hcoup
42
43
44
45mass_total = 0d0
46calvingGRIS=0d0
47calv(:,:)=0. 
48
49! formats pour les ecritures dans 42
50428 format(A)
51
52! lecture des parametres du run                      block eaubasale1
53!--------------------------------------------------------------------
54rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
55read(num_param,calving)
56
57
58write(num_rep_42,428)'!___________________________________________________________' 
59write(num_rep_42,428) '&calving              ! nom du bloc calving méthode Vincent '
60write(num_rep_42,*)
61write(num_rep_42,*)   'Hcoup      = ', Hcoup
62write(num_rep_42,*)   'ifrange    = ',ifrange
63write(num_rep_42,*)   'meth_hcoup = ', meth_hcoup
64write(num_rep_42,*)'/'
65                           
66write(num_rep_42,428) '! Hcoup epaisseur de coupure, le max si attache a coefbmshelf '
67write(num_rep_42,428) '! ifrange=0 -> pas de traitement particulier sur les bords'
68write(num_rep_42,428) '! ifrange=1 -> traitement de Vincent avec ice shelves frangeants partout'
69write(num_rep_42,428) '! ifrange=2 -> ice shelves frangeants seulement si bm-bmelt positif'
70write(num_rep_42,*)   '! meth_hcoup pour faire eventuellement varier Hcoup avec le climat'
71write(num_rep_42,*)   '! Hcoup_min=',Hcoup_min
72write(num_rep_42,*)
73
74Hcoup_max=Hcoup
75
76end subroutine init_calving
77!---------------------------------------------------------------------------------------
78subroutine calving
79 
80  real :: calvtest
81 
82! initialisation calving
83  calv(:,:)=0. 
84  mass_total=0d0
85  calvtest=0.
86  mass_total=0d0
87
88! calcul du dhdt lagrangien
89dhdt(:,:)=bm(:,:)-bmelt(:,:)-H(:,:)*(epsxx(:,:)+epsyy(:,:))
90
91! calcul du bilan surface et fond : divise par 2 car utilise dans des moyennes
92bil_tot(:,:)=0.5*(bm(:,:)-bmelt(:,:))
93
94! coupure
95if (meth_hcoup.eq.0) then
96   Hcoup=Hcoup_max
97else if (meth_hcoup.eq.1) then
98   Hcoup=coefbmshelf*Hcoup_max
99   Hcoup=min(Hcoup,Hcoup_max)
100   Hcoup=max(Hcoup,Hcoup_min)
101
102else if (meth_hcoup.eq.2) then
103   Hcoup=coefbmshelf*Hcoup_max
104   Hcoup=max(Hcoup,Hcoup_min)
105endif
106
107! hauteur au dessus de la coupure
108hmhc(:,:)=H(:,:)-Hcoup
109
110! coupure de l'ice shelf
111!---------------------------
112! on divise le test en 2 'ifext' et 'ifint' pour reduire les tests redondants 
113
114
115do j=2,ny-1
116   do i=2,nx-1
117
118
119ifext:  if((flot(i,j)).and.(h(i,j).le.hcoup))  then
120! ifext: pour les noeuds flottants avec h < hcoup
121
122ifint:    if(front(i,j).lt.4) then       
123! ifint: le point doit avoir au - 1 voisin mais ne pas etre entouré de glace
124! ce qui evite la formation des polynies dans les shelfs
125
126
127! on regarde dhdt minimum pour voir si 1 des voisins peut alimenter le point
128
129!                     hmhc est l'épaisseur en plus de l'épaisseur de coupure
130!                     ux/dx=1/(deltat) ou deltat est le temps nécessaire
131!                     à la glace pour passer d'un noeud à l'autre
132!                     la deuxieme partie du test revient à dh/dt*deltat > hmhc 
133
134!                     Rajout vince : si le point amont est pose -> test vrai.
135
136! rappel des differents ifrange
137!-------------------------------
138!Attention ifrange 1,2  ont sans doute un bug (il faudrait abs(uxbar))
139
140! ifrange=1               Rajout vince : si un point voisin est pose -> test vrai.
141!                         comme dans article fenno2007
142
143! ifrange=2               pour les points ayant des voisins poses, garde le point
144!                         seulement si le bilan  surface fond est positif.
145!                         bm(i,j)-bmelt(i,j) > 0.
146
147! ifrange=3               pas de test sur l'epaisseur du point amont.
148!                         test dhdt > -hmhc/deltat= -hmhc abs(U)/dx (correction du bug)
149!                         + test sur bm(i,j)-bmelt(i,j) > 0.
150
151! ifrange = 4             idem 3 mais avec test sur l'epaisseur du point amont
152
153! ifrange = 0             pas de traitement specifique près du continent
154
155! ifrange = -1            ancienne version MIS11
156
157! ifrange = 5             semble idem 0
158
159! ifrange = 6             ??
160
161! ifrange = 7             calving drastique sauf si coefbmshelf < 0.5
162
163!cdc
164! ifrange=44              idem 4 mais avec test sur i et j couple // afq: probably useless
165! ifrange=8               C. Dumas : coupure si alimentation ne permet pas H > Hcoup
166
167if (ifrange.eq.1) then           !   Rajout vince : si un point voisin est pose -> test vrai.
168                                 !   comme dans article fenno2007
169
170   testmij=( ((hmhc(i-1,j).gt.0.).and.(uxbar(i,j).ge.0.)  &  ! voisin (i-1,j) amont et > hcoup
171        .and.(dhdt(i-1,j).gt.(hmhc(i-1,j)*uxbar(i-1,j)/dx)))& 
172        .or.(.not.flot(i-1,j)) ) !
173
174   testpij=( ((hmhc(i+1,j).gt.0.).and.(uxbar(i+1,j).le.0.) & ! voisin (i+1,j) amont et > hcoup
175        .and.(dhdt(i+1,j).gt.(hmhc(i+1,j)*uxbar(i+1,j)/dx))) &
176        .or.(.not.flot(i+1,j)) ) !
177
178   testimj=( ((hmhc(i,j-1).gt.0.).and.(uybar(i,j).ge.0.)  &  ! voisin (i,j-1) amont et > hcoup
179        .and.(dhdt(i,j-1).gt.(hmhc(i,j-1)*uybar(i,j-1)/dx)))&
180        .or.(.not.flot(i,j-1)) ) !
181
182   testipj=( ((hmhc(i,j+1).gt.0.).and.(uybar(i,j+1).le.0.) & ! voisin (i,j+1) amont et > hcoup
183        .and.(dhdt(i,j+1).gt.(hmhc(i,j+1)*uybar(i,j+1)/dx)))&
184        .or.(.not.flot(i,j+1)) ) !
185
186else if (ifrange.eq.2) then       ! pour les points ayant des voisins posés, seulement si le bilan
187                                  ! surface fond est positif. bm(i,j)-bmelt(i,j) > 0.
188
189   bilan_surf_fond=(bm(i,j)-bmelt(i,j).gt.0.)
190
191   testmij=( ((hmhc(i-1,j).gt.0.).and.(uxbar(i,j).ge.0.)  &  ! voisin (i-1,j) amont et > hcoup
192        .and.  (dhdt(i-1,j).gt.(hmhc(i-1,j)*uxbar(i-1,j)/dx))) & 
193        .or.((.not.flot(i-1,j)).and.bilan_surf_fond )) !
194
195   testpij=( ((hmhc(i+1,j).gt.0.).and.(uxbar(i+1,j).le.0.) & ! voisin (i+1,j) amont et > hcoup
196        .and.(dhdt(i+1,j).gt.(hmhc(i+1,j)*uxbar(i+1,j)/dx))) &
197        .or.((.not.flot(i+1,j)).and.bilan_surf_fond) ) !
198
199   testimj=( ((hmhc(i,j-1).gt.0.).and.(uybar(i,j).ge.0.)  &  ! voisin (i,j-1) amont et > hcoup
200        .and.(dhdt(i,j-1).gt.(hmhc(i,j-1)*uybar(i,j-1)/dx)))&
201        .or.((.not.flot(i,j-1)).and.bilan_surf_fond ) ) !
202
203   testipj=( ((hmhc(i,j+1).gt.0.).and.(uybar(i,j+1).le.0.) & ! voisin (i,j+1) amont et > hcoup
204        .and.(dhdt(i,j+1).gt.(hmhc(i,j+1)*uybar(i,j+1)/dx)))&
205        .or.((.not.flot(i,j+1)).and.bilan_surf_fond ) ) !
206
207else if (ifrange.eq.3) then       ! nouvelle formulation Cat mars 08
208
209! pas de test sur l'epaisseur du point amont.
210! test dhdt > -hmhc/deltat= -hmhc abs(U)/dx
211
212! pour les points ayant des voisins posés, seulement si le bilan
213! surface fond est positif. bm(i,j)-bmelt(i,j) > 0.
214
215
216   bilan_surf_fond=(bm(i,j)-bmelt(i,j).gt.0.)
217
218   testmij=( (uxbar(i,j).ge.0.) &  ! voisin (i-1,j) amont et > hcoup
219        .and.(dhdt(i-1,j).gt.(-hmhc(i-1,j)*abs(uxbar(i-1,j))/dx))) & 
220        .or.((.not.flot(i-1,j)).and.bilan_surf_fond ) !
221
222   testpij=( (uxbar(i+1,j).le.0.) & ! voisin (i+1,j) amont et > hcoup
223        .and.(dhdt(i+1,j).gt.(-hmhc(i+1,j)*abs(uxbar(i+1,j))/dx))) &
224        .or.((.not.flot(i+1,j)).and.bilan_surf_fond )  !
225
226   testimj=(( uybar(i,j).ge.0.)  &  ! voisin (i,j-1) amont et > hcoup
227        .and.(dhdt(i,j-1).gt.(-hmhc(i,j-1)*abs(uybar(i,j-1))/dx))) &
228        .or.((.not.flot(i,j-1)).and.bilan_surf_fond )  !
229
230   testipj=((uybar(i,j+1).le.0.) & ! voisin (i,j+1) amont et > hcoup
231        .and.(dhdt(i,j+1).gt.(-hmhc(i,j+1)*abs(uybar(i,j+1))/dx))) &
232        .or.((.not.flot(i,j+1)).and.bilan_surf_fond )  !
233
234else if (ifrange.eq.4) then       ! idem 3 mais avec test sur l'épaisseur du point amont
235
236   bilan_surf_fond=(bm(i,j)-bmelt(i,j).gt.0.)
237
238   testmij=( ((hmhc(i-1,j).gt.0.).and.(uxbar(i,j).ge.0.)  &  ! voisin (i-1,j) amont et > hcoup
239        .and.  (dhdt(i-1,j).gt.(-hmhc(i-1,j)*abs(uxbar(i-1,j)/dx)))) & 
240        .or.((.not.flot(i-1,j)).and.bilan_surf_fond )) !
241
242   testpij=( ((hmhc(i+1,j).gt.0.).and.(uxbar(i+1,j).le.0.) & ! voisin (i+1,j) amont et > hcoup
243        .and.(dhdt(i+1,j).gt.(-hmhc(i+1,j)*abs(uxbar(i+1,j)/dx)))) &
244        .or.((.not.flot(i+1,j)).and.bilan_surf_fond) ) !
245
246   testimj=( ((hmhc(i,j-1).gt.0.).and.(uybar(i,j).ge.0.)  &  ! voisin (i,j-1) amont et > hcoup
247        .and.(dhdt(i,j-1).gt.(-hmhc(i,j-1)*abs(uybar(i,j-1)/dx))))&
248        .or.((.not.flot(i,j-1)).and.bilan_surf_fond ) ) !
249
250   testipj=( ((hmhc(i,j+1).gt.0.).and.(uybar(i,j+1).le.0.) & ! voisin (i,j+1) amont et > hcoup
251        .and.(dhdt(i,j+1).gt.(-hmhc(i,j+1)*abs(uybar(i,j+1)/dx))))&
252        .or.((.not.flot(i,j+1)).and.bilan_surf_fond ) ) !
253
254else if (ifrange.eq.0) then           ! pas de traitement special pres du continent
255
256   testmij= ((hmhc(i-1,j).gt.0.).and.(uxbar(i,j).ge.0.)  &  ! voisin (i-1,j) amont et > hcoup
257        .and.(dhdt(i-1,j).gt.(-hmhc(i-1,j)*abs(uxbar(i-1,j)/dx))))
258
259   testpij= ((hmhc(i+1,j).gt.0.).and.(uxbar(i+1,j).le.0.) & ! voisin (i+1,j) amont et > hcoup
260        .and.(dhdt(i+1,j).gt.(-hmhc(i+1,j)*abs(uxbar(i+1,j))/dx)))
261
262   testimj= ((hmhc(i,j-1).gt.0.).and.(uybar(i,j).ge.0.)  &  ! voisin (i,j-1) amont et > hcoup
263        .and.(dhdt(i,j-1).gt.(-hmhc(i,j-1)*abs(uybar(i,j-1))/dx)))
264
265   testipj= ((hmhc(i,j+1).gt.0.).and.(uybar(i,j+1).le.0.) & ! voisin (i,j+1) amont et > hcoup
266        .and.(dhdt(i,j+1).gt.(-hmhc(i,j+1)*abs(uybar(i,j+1))/dx)))
267
268else if (ifrange.eq.-1) then           ! ancienne version MIS11
269
270testmij=((hmhc(i-1,j).gt.0.).and.(uxbar(i,j).ge.0.)  &  ! voisin (i-1,j) amont et > hcoup
271   .and.(dhdt(i-1,j).gt.(hmhc(i-1,j)*uxbar(i-1,j)/dx)))   
272
273testpij=((hmhc(i+1,j).gt.0.).and.(uxbar(i+1,j).le.0.) & ! voisin (i+1,j) amont et > hcoup
274   .and.(dhdt(i+1,j).gt.(hmhc(i+1,j)*uxbar(i+1,j)/dx)))
275
276testimj=((hmhc(i,j-1).gt.0.).and.(uybar(i,j).ge.0.)  &  ! voisin (i,j-1) amont et > hcoup
277   .and.(dhdt(i,j-1).gt.(hmhc(i,j-1)*uybar(i,j-1)/dx)))
278
279testipj=((hmhc(i,j+1).gt.0.).and.(uybar(i,j+1).le.0.) & ! voisin (i,j+1) amont et > hcoup
280    .and.(dhdt(i,j+1).gt.(hmhc(i,j+1)*uybar(i,j+1)/dx)))
281
282
283else if (ifrange.eq.5) then           ! pas de traitement special pres du continent
284
285   testmij= ((hmhc(i-1,j).gt.0.).and.(uxbar(i,j).ge.0.)  &  ! voisin (i-1,j) amont et > hcoup
286        .and.((bil_tot(i-1,j)+bil_tot(i,j)) &
287        .gt.(-hmhc(i-1,j)*abs(uxbar(i-1,j)/dx))))
288
289   testpij= ((hmhc(i+1,j).gt.0.).and.(uxbar(i+1,j).le.0.) & ! voisin (i+1,j) amont et > hcoup
290        .and.((bil_tot(i+1,j)+bil_tot(i,j)) &
291        .gt.(-hmhc(i+1,j)*abs(uxbar(i+1,j))/dx)))
292
293   testimj= ((hmhc(i,j-1).gt.0.).and.(uybar(i,j).ge.0.)  &  ! voisin (i,j-1) amont et > hcoup
294        .and.((bil_tot(i,j-1)+bil_tot(i,j)) &
295        .gt.(-hmhc(i,j-1)*abs(uybar(i,j-1))/dx)))
296
297   testipj= ((hmhc(i,j+1).gt.0.).and.(uybar(i,j+1).le.0.) & ! voisin (i,j+1) amont et > hcoup
298        .and.((bil_tot(i,j+1)+bil_tot(i,j)) &
299        .gt.(-hmhc(i,j+1)*abs(uybar(i,j+1))/dx)))
300
301else if (ifrange.eq.6) then           ! pas de traitement special pres du continent
302
303   testmij= ((hmhc(i-1,j).gt.0.).and.(uxbar(i,j).ge.0.)  &  ! voisin (i-1,j) amont et > hcoup
304        .and.(2.*bil_tot(i,j) &
305        .gt.(-hmhc(i-1,j)*abs(uxbar(i-1,j)/dx))))
306
307   testpij= ((hmhc(i+1,j).gt.0.).and.(uxbar(i+1,j).le.0.) & ! voisin (i+1,j) amont et > hcoup
308        .and.(2.*bil_tot(i,j) &
309        .gt.(-hmhc(i+1,j)*abs(uxbar(i+1,j))/dx)))
310
311   testimj= ((hmhc(i,j-1).gt.0.).and.(uybar(i,j).ge.0.)  &  ! voisin (i,j-1) amont et > hcoup
312        .and.(2.*bil_tot(i,j) &
313        .gt.(-hmhc(i,j-1)*abs(uybar(i,j-1))/dx)))
314
315   testipj= ((hmhc(i,j+1).gt.0.).and.(uybar(i,j+1).le.0.) & ! voisin (i,j+1) amont et > hcoup
316        .and.(2.*bil_tot(i,j) &
317        .gt.(-hmhc(i,j+1)*abs(uybar(i,j+1))/dx)))
318
319else if (ifrange.eq.7) then           ! drastique calving sauf si coefbmelt.lt.0.5
320
321   testmij= ((hmhc(i-1,j).gt.0.).and.(uxbar(i,j).ge.0.)  &  ! voisin (i-1,j) amont et > hcoup
322        .and.(coefbmshelf.lt.1.))
323   testpij= ((hmhc(i+1,j).gt.0.).and.(uxbar(i+1,j).le.0.) & ! voisin (i+1,j) amont et > hcoup
324        .and.(coefbmshelf.lt.1.))
325   testimj= ((hmhc(i,j-1).gt.0.).and.(uybar(i,j).ge.0.)  &  ! voisin (i,j-1) amont et > hcoup
326         .and.(coefbmshelf.lt.1.))
327   testipj= ((hmhc(i,j+1).gt.0.).and.(uybar(i,j+1).le.0.) & ! voisin (i,j+1) amont et > hcoup
328        .and.(coefbmshelf.lt.1.))
329
330else if (ifrange.eq.44) then       ! idem 4 mais avec test sur i et j couple
331
332   bilan_surf_fond=(bm(i,j)-bmelt(i,j).gt.0.)
333
334   testmij=( ((hmhc(i-1,j).gt.0.).and.(uxbar(i,j).ge.0.)  &  ! voisin (i-1,j) amont et > hcoup
335        .and.  (dhdt(i-1,j).gt.(-hmhc(i-1,j)*abs(uxbar(i-1,j)/dx)))) & 
336        .or.((.not.flot(i-1,j)).and.bilan_surf_fond )) !
337
338   testpij=( ((hmhc(i+1,j).gt.0.).and.(uxbar(i+1,j).le.0.) & ! voisin (i+1,j) amont et > hcoup
339        .and.(dhdt(i+1,j).gt.(-hmhc(i+1,j)*abs(uxbar(i+1,j)/dx)))) &
340        .or.((.not.flot(i+1,j)).and.bilan_surf_fond) ) !
341
342   testimj=( ((hmhc(i,j-1).gt.0.).and.(uybar(i,j).ge.0.)  &  ! voisin (i,j-1) amont et > hcoup
343        .and.(dhdt(i,j-1).gt.(-hmhc(i,j-1)*abs(uybar(i,j-1)/dx))))&
344        .or.((.not.flot(i,j-1)).and.bilan_surf_fond ) ) !
345
346   testipj=( ((hmhc(i,j+1).gt.0.).and.(uybar(i,j+1).le.0.) & ! voisin (i,j+1) amont et > hcoup
347        .and.(dhdt(i,j+1).gt.(-hmhc(i,j+1)*abs(uybar(i,j+1)/dx))))&
348        .or.((.not.flot(i,j+1)).and.bilan_surf_fond ) ) !
349
350else if (ifrange.eq.8) then ! test C Dumas : si somme flux entrant permet H > Hcoup : vrai
351                   !   bilan_surf_fond=(bm(i,j)-bmelt(i,j).gt.0.)
352
353   testmij=( ((hmhc(i-1,j).gt.0.).and.(uxbar(i,j).ge.0.))  &  ! voisin (i-1,j) amont et > hcoup
354                        !        .and.  (dhdt(i-1,j).gt.(-hmhc(i-1,j)*abs(uxbar(i-1,j)/dx)))) &
355        .or.((.not.flot(i-1,j)).and.bilan_surf_fond )) !
356
357   testpij=( ((hmhc(i+1,j).gt.0.).and.(uxbar(i+1,j).le.0.)) & ! voisin (i+1,j) amont et > hcoup
358                        !        .and.(dhdt(i+1,j).gt.(-hmhc(i+1,j)*abs(uxbar(i+1,j)/dx)))) &
359        .or.((.not.flot(i+1,j)).and.bilan_surf_fond) ) !
360
361   testimj=( ((hmhc(i,j-1).gt.0.).and.(uybar(i,j).ge.0.))  &  ! voisin (i,j-1) amont et > hcoup
362                        !        .and.(dhdt(i,j-1).gt.(-hmhc(i,j-1)*abs(uybar(i,j-1)/dx))))&
363        .or.((.not.flot(i,j-1)).and.bilan_surf_fond ) ) !
364
365   testipj=( ((hmhc(i,j+1).gt.0.).and.(uybar(i,j+1).le.0.)) & ! voisin (i,j+1) amont et > hcoup
366                        !        .and.(dhdt(i,j+1).gt.(-hmhc(i,j+1)*abs(uybar(i,j+1)/dx))))&
367        .or.((.not.flot(i,j+1)).and.bilan_surf_fond ) ) !
368
369   testH=0.
370   if (testmij) testH=testH + abs(uxbar(i,j))*H(i-1,j)
371   if (testpij) testH=testH + abs(uxbar(i+1,j))*H(i+1,j)
372   if (testimj) testH=testH + abs(uybar(i,j))*H(i,j-1)
373   if (testipj) testH=testH + abs(uybar(i,j+1))*H(i,j+1)
374   
375   if (testH.gt.Hcoup) then
376      testmij=.true.
377   else
378      testmij=.false.
379      testpij=.false.
380      testimj=.false.
381      testipj=.false.
382   endif
383                   
384endif
385
386
387
388! detection des polynies dans les shelfs. on regarde vers l'aval
389
390! 1 des 3 voisins aval > hcoup du cote west (i-1)
391iin2=max(1,i-2)
392iin3=max(1,i-3)
393
394!cdc avalw=((hmhc(i-1,j).gt.0.).or.(hmhc(iin2,j).gt.0.) &
395!cdc      .or.(hmhc(iin3,j).gt.0.)).and.(uxbar(i,j).le.0.)
396
397avalw=((hmhc(i-1,j).gt.0..and.flot(i-1,j)).or.(hmhc(iin2,j).gt.0..and.flot(iin2,j)) &
398     .or.(hmhc(iin3,j).gt.0..and.flot(iin3,j))).and.(uxbar(i,j).le.0.)   
399
400! 1 des 3 voisins aval > hcoup du cote est (i+1)
401iin2=min(i+2,nx)
402iin3=min(i+3,nx)
403
404!cdc avale=((hmhc(i+1,j).gt.0.).or.(hmhc(iin2,j).gt.0.) &   
405!cdc      .or.(hmhc(iin3,j).gt.0.)).and.(uxbar(i+1,j).ge.0)                         
406   
407avale=((hmhc(i+1,j).gt.0..and.flot(i+1,j)).or.(hmhc(iin2,j).gt.0..and.flot(iin2,j)) &   
408     .or.(hmhc(iin3,j).gt.0..and.flot(iin3,j))).and.(uxbar(i+1,j).ge.0)
409               
410! 1 des 3 voisins aval > hcoup du cote sud (j-1)
411jin2=max(1,j-2)
412jin3=max(1,j-3)
413
414!cdc avals=((hmhc(i,j-1).gt.0.).or.(hmhc(i,jin2).gt.0.) &   
415!cdc     .or.(hmhc(i,jin3).gt.0.)).and.(uybar(i,j).le.0.) 
416                         
417avals=((hmhc(i,j-1).gt.0..and.flot(i,j-1)).or.(hmhc(i,jin2).gt.0..and.flot(i,jin2)) &   
418     .or.(hmhc(i,jin3).gt.0..and.flot(i,jin3))).and.(uybar(i,j).le.0.)
419               
420! 1 des 3 voisins aval > hcoup du cote nord (j-1)
421jin2=min(j+2,ny)
422jin3=min(j+3,ny)
423
424!cdc avaln=((hmhc(i,j+1).gt.0.).or.(hmhc(i,jin2).gt.0.) &   
425!cdc     .or.(hmhc(i,jin2).gt.0.)) .and.(uybar(i,j+1).ge.0.)                         
426
427avaln=((hmhc(i,j+1).gt.0..and.flot(i,j+1)).or.(hmhc(i,jin2).gt.0..and.flot(i,jin2)) &   
428     .or.(hmhc(i,jin2).gt.0..and.flot(i,jin2))) .and.(uybar(i,j+1).ge.0.)
429
430interieur=(avalw.or.avale).and.(avals.or.avaln)
431
432
433
434!cdc if ((.not.(testmij.or.testpij.or.testimj.or.testipj))  &  ! pas suffisament alimente
435!cdc                .and.(.not.interieur)) then                ! et pas interieur
436
437!cdc pour ne pas faire de calving si H > 350m :   if ((.not.(testmij.or.testpij.or.testimj.or.testipj)).and.(H(i,j).LT.350)) then
438! mab: if none of the neighbours doesn't provide ice for the considered point,
439! mab: then the ice is cut at this point and given to calv(i,j)
440                if ((.not.(testmij.or.testpij.or.testimj.or.testipj))) then
441! mab: new version so that grid points with h=1m (all points were calving takes place are afterwards
442! mab: set to a height of 1m) are not taken into account,  provides a problem as soon as you start
443! mab: coupling GRISLI to the iceberg modell
444                   if(h(i,j).gt.1) then
445                     calv(i,j)=-h(i,j) 
446! mab: copy the calving in the right format for coupling                     
447! calvin_GRIS =  volume flux per year (m3/a)
448                     calvin_GRIS(i,j)=calvin_GRIS(i,j)+(-calv(i,j)*dx*dy)
449!                     calvin_GRIS(i,j)=-calv(i,j)*40000.*40000.
450!                     write(15,*) 'calvheight:',h(i,j),calv(i,j),i,j,dtt,dt,calvin_GRIS(i,j)         
451                     endif 
452                   H(i,j)=1             
453                endif
454
455                  end if ifint
456              end if ifext
457         end do
458      end do
459         
460    calvtest=sum(-calv(:,:)*dx*dy)
461    mass_total=sum(calvingGRIS(:,:))
462
463    calvtest = 0.
464
465   
466    end subroutine calving
467!------------------------------------------------------------------------------------------
468  end module calving_frange
469 
Note: See TracBrowser for help on using the repository browser.