source: branches/GRISLIv3/SOURCES/calving_frange_abuk.f90 @ 455

Last change on this file since 455 was 446, checked in by aquiquet, 9 months ago

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