source: trunk/SOURCES/calving_frange_abuk.f90 @ 334

Last change on this file since 334 was 237, checked in by aquiquet, 5 years ago

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