source: branches/GRISLIv3/SOURCES/calving_frange_ISMIP_fracture.f90 @ 443

Last change on this file since 443 was 385, checked in by dumas, 16 months ago

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