source: trunk/SOURCES/calving_frange_ISMIP_fracture.f90 @ 334

Last change on this file since 334 was 270, checked in by dumas, 5 years ago

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