Changeset 179


Ignore:
Timestamp:
02/01/18 10:18:05 (6 years ago)
Author:
dumas
Message:

Calving update : loop to allow calving of several points from the edge

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/SOURCES/calving_frange.f90

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