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