Changeset 86


Ignore:
Timestamp:
09/15/16 17:52:05 (8 years ago)
Author:
dumas
Message:

Bug correction in flottab : call to determin_tache suppressed | sealevel added in Netcdf output class 1

Location:
trunk/SOURCES
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/SOURCES/Ant40_files/lect-anteis_mod.f90

    r85 r86  
    114114     call Read_Ncdf_var('H',topo_ref,tab) 
    115115     H0(:,:)  = tab(:,:) 
    116 !cdc correction de 3 pts qui posent probleme (pente tres forte): 
    117      S0(38,53)=1400. 
    118      Bsoc0(38,53)=-1000. 
    119      H0(38,53)=S0(38,53)-Bsoc0(38,53) 
    120       
    121      S0(35,53)=1300. 
    122      Bsoc0(35,53)=-1000. 
    123      H0(35,53)=S0(35,53)-Bsoc0(35,53) 
    124       
    125      S0(35,56)=1200. 
    126      Bsoc0(35,56)=-1000. 
    127      H0(35,56)=S0(35,56)-Bsoc0(35,56) 
     116!~ !cdc correction de 3 pts qui posent probleme (pente tres forte): 
     117!~      S0(38,53)=1400. 
     118!~      Bsoc0(38,53)=-1000. 
     119!~      H0(38,53)=S0(38,53)-Bsoc0(38,53) 
     120!~       
     121!~      S0(35,53)=1300. 
     122!~      Bsoc0(35,53)=-1000. 
     123!~      H0(35,53)=S0(35,53)-Bsoc0(35,53) 
     124!~       
     125!~      S0(35,56)=1200. 
     126!~      Bsoc0(35,56)=-1000. 
     127!~      H0(35,56)=S0(35,56)-Bsoc0(35,56) 
    128128 
    129129!cdc correction point pole sud : 
    … …  
    161161     H(:,:)  = tab(:,:) 
    162162      
    163 !cdc correction de 3 pts qui posent probleme (pente tres forte): 
    164 !~     do i=34,39 
    165 !~       do j=52,57 
    166 !~         print*,'i j', i, j 
    167 !~         print*,'SHB',S(i,j),Bsoc(i,j),H(i,j) 
    168 !~       enddo 
    169 !~    enddo 
    170  
    171      S(39,54)=1400. 
    172      Bsoc(39,54)=-1000. 
    173      H(39,54)=S(39,54)-Bsoc(39,54) 
    174       
    175      S(36,54)=1300. 
    176      Bsoc(36,54)=-1000. 
    177      H(36,54)=S(36,54)-Bsoc(36,54) 
    178       
    179      S(36,57)=1200. 
    180      Bsoc(36,57)=-1000. 
    181      H(36,57)=S(36,57)-Bsoc(36,57) 
     163!~ !cdc correction de 3 pts qui posent probleme (pente tres forte): 
     164!~      S(39,54)=1400. 
     165!~      Bsoc(39,54)=-1000. 
     166!~      H(39,54)=S(39,54)-Bsoc(39,54) 
     167!~       
     168!~      S(36,54)=1300. 
     169!~      Bsoc(36,54)=-1000. 
     170!~      H(36,54)=S(36,54)-Bsoc(36,54) 
     171!~       
     172!~      S(36,57)=1200. 
     173!~      Bsoc(36,57)=-1000. 
     174!~      H(36,57)=S(36,57)-Bsoc(36,57) 
    182175     
    183176!    S(71,71)=(S(71,70)+S(71,72)+S(70,71)+S(72,71))/4. 
    … …  
    235228    do J=1,NY 
    236229       do I=1,NX 
    237           if ((BSOC(I,J)+H(I,J)*RO/ROW -SEALEVEL).LT.0.) then 
     230          if (((BSOC(I,J)+H(I,J)*RO/ROW -SEALEVEL).LT.0.).and.(H(I,J).gt.1.E-3)) then 
    238231             FLOT(I,J)=.TRUE. 
    239232          else 

  • trunk/SOURCES/Fichiers-parametres/hemin40_TEMPS-NETCDF.dat

    r60 r86  
    1515------------------------------------------------------------------------------------------- 
    16164       ndtsortie      lecture de dtsortie_hz 
    17 1000. 
     175000. 
    181850000. 
    1919100000. 
    20201.e10 
    2121------------------------------------------------------------------------------------------- 
    22 15      npredeft        lecture de predef_tsort 
     224       npredeft        lecture de predef_tsort 
    232310 
    2424100 
    2525200 
    2626500 
    27 100010 
    28 100100 
    29 200010 
    30 200100 
    31 -20990 
    32 -20900 
    33 -20500 
    34 -19500 
    35 -121990 
    36 -125990 
    37 -125900 
    38 -419990 
    39 -419500 
    40 -419000 
    41  
    42  
    4327------------------------------------------------------------------------------------------- 

  • trunk/SOURCES/GrIce2sea_files/climat_GrIce2sea_years_perturb_mod.f90

    r68 r86  
    425425        coefbmshelf=max(coefbmshelf,0.)   
    426426        coefbmshelf=min(coefbmshelf,2.) 
     427         
     428!        coefbmshelf=1.  !afq for calibration Greenland 
    427429   
    428430     case(1) 

  • trunk/SOURCES/Hemin40_files/module_choix-hemin40-0.4.f90

    r82 r86  
    4747!--------------Module climat --------------- 
    4848!use climat_forcage_mois_mod ! forcage mensuel GCM 1 Snapshot Fev 2015 
    49 use climat_Grice2sea_years_mod ! forcage avec SMB type MAR Fev2015 
    50 !use climat_Grice2sea_years_perturb_mod ! forcage avec SMB type MAR + perturbation temperature 
     49!use climat_Grice2sea_years_mod ! forcage avec SMB type MAR Fev2015 
     50use climat_Grice2sea_years_perturb_mod ! forcage avec SMB type MAR + perturbation temperature 
    5151 
    5252! Anciens modules pas encore valides 

  • trunk/SOURCES/Netcdf-routines/sortie_netcdf_GRISLI_mod.0.2-hassine.f90

    r66 r86  
    111111                                                             !< des variables dans chaque fichier  
    112112  integer,dimension(:),allocatable:: idef_time               !< 1 ou 0 si le temp est defini ou non 
     113  integer,dimension(:),allocatable:: idef_sealevel           !< 1 ou 0 si sealevel est defini ou non 
    113114  integer,dimension(:,:),allocatable:: idef                  !< 1 ou 0 si la varaible est defini ou non 
    114115  integer,dimension(:),allocatable :: num_ncdf_file          !< compteur des fichiers netcdf par class   
    … …  
    504505    if (.not.allocated(nrecs) .and. .not.allocated(nbsnap) .and. .not.allocated(num_ncdf_file) & 
    505506         .and. .not. allocated(idef) .and. .not. allocated(idef_time)) then 
    506        allocate(nrecs(nclassout),nbsnap(nclassout),num_ncdf_file(nclassout),idef(nclassout,ntab),idef_time(nclassout)) 
     507       allocate(nrecs(nclassout),nbsnap(nclassout),num_ncdf_file(nclassout),idef(nclassout,ntab),idef_time(nclassout),idef_sealevel(nclassout)) 
    507508       nrecs=1 
    508509       idef=0 
    509510       idef_time=0 
     511       idef_sealevel=0 
    510512       num_ncdf_file=0 
    511513    end if 
    … …  
    516518          idef(j,:)=0 
    517519          idef_time(j)=0 
     520          idef_sealevel(j)=0 
    518521          nbsnap(j)=0 
    519522          ! numerote le fichier sortie 
    … …  
    558561       idef(posis,:)=0 
    559562       idef_time(posis)=0 
     563       idef_sealevel(posis)=0 
    560564       nbsnap(posis)=0 
    561565       ! numerote le fichier sortie 
    … …  
    654658    character(len=20) :: nametmp                         !< nom intermediaire 
    655659    real*8,pointer,dimension(:) :: liste_time => null()  !< liste des snapshot des variables ecrites en netcdf 
     660    real*8,pointer,dimension(:) :: sealevel_p => null()  !< pointeur vers sealevel pour ecriture netcdf 
    656661    real*8,pointer,dimension(:) :: x,y,x1,y1,z,nzzm    
    657662    real*8,pointer,dimension(:,:):: lat,lon => null() 
    … …  
    683688       liste_time(1)=-1 
    684689    end if 
     690    if (.not.associated(sealevel_p)) then 
     691       allocate(sealevel_p(1)) 
     692       sealevel_p(1)=-1 
     693    end if 
    685694 
    686695liste_times:  if ((liste_time(1) .ne.timetmp) .or.(liste_time(1) .eq. -1) ) then 
    687696       liste_time(1)= timetmp 
     697       sealevel_p(1)= sealevel 
    688698       ! print*,"time outncdf=",liste_time(1) 
    689699       write(charint,'(i0)') floor(timetmp) 
    … …  
    720730! ecrit le temps 
    721731          call write_ncdf_var('time','time',trim(fil_sortie(k)),liste_time,nbsnap(k)+1,idef_time(k),'double') 
     732           
     733          sealevel_p=sealevel 
     734          call write_ncdf_var('sealevel','time',trim(fil_sortie(k)),sealevel_p,nrecs(k),idef_sealevel(k),'double') 
    722735 
    723736          fait = .FALSE. 

  • trunk/SOURCES/flottab2-0.7.f90

    r73 r86  
    7878!  @note    - flottant  flot sur le noeud majeur, flotmx sur le noeud mineur 
    7979!> 
    80        subroutine flottab() 
    81 ! 
    82 !                                       Vince 5 Jan 95 
    83 !                                       Modifie 20 Jan 95 
    84 !                                       Modifie le 30 Novembre 98 
    85 !    Passage f90 + determination des fronts Vincent dec 2003 
    86 !    nettoyage et nouvelle détermination de gzmx et gzmy Cat le 10 juillet 2005 
    87 !    Re-nettoyage par Cat en aout 2006. 
    88 !    Le calcul de gzmx et gzmy pour les points intérieurs passe dans la subroutine dragging 
    89 !    pour les points cotiers, toujours fait dans flottab 
    90 ! 
    91 !                                       ----------------- 
    92 ! 
    93 !      Cette routine determine les endroits ou la glace 
    94 !      flotte , les iles, et la position du front de glace 
    95 ! 
    96 !      Il y a 4 sortes de zone 
    97 !      Pose 
    98 !      Grounding zone et streams    gzmx et gzmy 
    99 !      Iles             ilemx, ilemy 
    100 !      flottant  flot sur le noeud majeur, flotmx sur le noeud mineur 
    101  
    102 !   passage dans flottab tous les pas de temps dt  !  
    103 ! 
    104 ! _________________________________________________________ 
    105  
    106  
    107  
    108 if (itracebug.eq.1)  call tracebug(' Entree dans routine flottab') 
    109  
    110  SHELFY = .FALSE. 
    111  
    112  
    113 ! cas particulier des runs paleo ou on impose un masque grounded 
    114  
    115 !$OMP PARALLEL PRIVATE(archim,surnet) 
    116  if (igrdline.eq.2) then 
     80  subroutine flottab() 
     81    ! 
     82    !                                   Vince 5 Jan 95 
     83    !                                       Modifie 20 Jan 95 
     84    !                                       Modifie le 30 Novembre 98 
     85    !    Passage f90 + determination des fronts Vincent dec 2003 
     86    !    nettoyage et nouvelle détermination de gzmx et gzmy Cat le 10 juillet 2005 
     87    !    Re-nettoyage par Cat en aout 2006. 
     88    !    Le calcul de gzmx et gzmy pour les points intérieurs passe dans la subroutine dragging 
     89    !    pour les points cotiers, toujours fait dans flottab 
     90    ! 
     91    !                                       ----------------- 
     92    ! 
     93    !      Cette routine determine les endroits ou la glace 
     94    !      flotte , les iles, et la position du front de glace 
     95    ! 
     96    !      Il y a 4 sortes de zone 
     97    !      Pose 
     98    !      Grounding zone et streams    gzmx et gzmy 
     99    !      Iles             ilemx, ilemy 
     100    !      flottant  flot sur le noeud majeur, flotmx sur le noeud mineur 
     101 
     102    !   passage dans flottab tous les pas de temps dt  !  
     103    ! 
     104    ! _________________________________________________________ 
     105 
     106    !~ print*,'debut flottab',S(132,183),H(132,183),BSOC(132,183),B(132,183),sealevel 
     107    !~ print*,'debut flottab',flot(132,183),ice(132,183) 
     108 
     109    if (itracebug.eq.1)  call tracebug(' Entree dans routine flottab') 
     110 
     111    SHELFY = .FALSE. 
     112 
     113 
     114    ! cas particulier des runs paleo ou on impose un masque grounded 
     115 
     116    !$OMP PARALLEL PRIVATE(archim,surnet) 
     117    if (igrdline.eq.2) then 
     118       !$OMP WORKSHARE 
     119       where ( mk_init(:,:).eq.1)                 ! pose 
     120          flot(:,:) = .False. 
     121          H(:,:)=max(H(:,:),(10.+sealevel-Bsoc(:,:))*row/ro)  ! pour avoir archim=10 
     122          S(:,:) = Bsoc(:,:) + H(:,:)  
     123       elsewhere 
     124          flot(:,:) = .True. 
     125       end where 
     126       !$OMP END WORKSHARE 
     127    end if 
     128 
     129 
     130    ! 1-INITIALISATION 
     131    ! ---------------- 
     132    ! initialisation des variables pour detecter les points qui se mettent 
     133    ! a flotter entre 2 dtt 
     134 
     135    appel_new_flot=.false. 
     136    !$OMP DO 
     137    do j=1,ny 
     138       do i=1,nx 
     139          new_flot_point(i,j)=.false. 
     140          new_flotmx(i,j)=.false. 
     141          new_flotmy(i,j)=.false. 
     142       enddo 
     143    enddo 
     144    !$OMP END DO 
     145 
     146    ! ICE(:,:)=(H(:,:).gt.1) ! ice=.true. si epaisseur > 1m 
     147 
    117148    !$OMP WORKSHARE 
    118     where ( mk_init(:,:).eq.1)                 ! pose 
    119        flot(:,:) = .False. 
    120        H(:,:)=max(H(:,:),(10.+sealevel-Bsoc(:,:))*row/ro)  ! pour avoir archim=10 
    121        S(:,:) = Bsoc(:,:) + H(:,:)  
    122     elsewhere 
    123        flot(:,:) = .True. 
    124     end where 
     149    ICE(:,:)=0 
     150    front(:,:)=0 
     151    frontfacex(:,:)=0 
     152    frontfacey(:,:)=0 
     153    isolx(:,:)=.FALSE. 
     154    isoly(:,:)=.FALSE. 
     155    cotemx(:,:)=.false. 
     156    cotemy(:,:)=.false. 
     157    boost=.false. 
    125158    !$OMP END WORKSHARE 
    126  end if 
    127  
    128  
    129 ! 1-INITIALISATION 
    130 ! ---------------- 
    131 ! initialisation des variables pour detecter les points qui se mettent 
    132 ! a flotter entre 2 dtt 
    133  
    134        appel_new_flot=.false. 
    135        !$OMP DO 
    136        do j=1,ny 
    137           do i=1,nx 
    138              new_flot_point(i,j)=.false. 
    139              new_flotmx(i,j)=.false. 
    140              new_flotmy(i,j)=.false. 
    141           enddo 
    142        enddo 
    143        !$OMP END DO 
    144         
    145 ! ICE(:,:)=(H(:,:).gt.1) ! ice=.true. si epaisseur > 1m 
    146  
    147       !$OMP WORKSHARE 
    148       ICE(:,:)=0 
    149       front(:,:)=0 
    150       frontfacex(:,:)=0 
    151       frontfacey(:,:)=0 
    152       isolx(:,:)=.FALSE. 
    153       isoly(:,:)=.FALSE. 
    154       cotemx(:,:)=.false. 
    155       cotemy(:,:)=.false. 
    156       boost=.false. 
    157       !$OMP END WORKSHARE 
    158  
    159 ! fin de l'initialisation 
    160 !_____________________________________________________________________ 
    161  
    162 ! 2-TESTE LES NOUVEAUX POINTS FLOTTANTS 
    163 ! ------------------------------------- 
    164  
    165  !$OMP DO 
    166  do j=1,ny 
    167    do i=1,nx 
    168  
    169       uxs1(i,j)=uxbar(i,j) 
    170       uys1(i,j)=uybar(i,j) 
    171  
    172       archim = Bsoc(i,j)+H(i,j)*ro/row -sealevel 
    173  
    174  
    175 arch: if ((ARCHIM.LT.0.).and.(H(I,J).gt.1.E-3)) then    ! le point flotte 
    176    mk(i,j)=1 
    177  
    178  
    179 ex_pose: if ((.not.FLOT(I,J)).and.(isynchro.eq.1)) then  !  il ne flottait pas avant 
    180             FLOT(I,J)=.true. 
    181             BOOST=.false. 
    182  
    183 ! Attention :  avec le bloc dessous il faut faire un calcul de flottaison a la lecture  
    184   
    185              if (igrdline.eq.1) then   ! en cas de grounding line prescrite 
    186                 flot(i,j)=.false. 
    187                 H(i,j)=(10.+sealevel-Bsoc(i,j))*row/ro  ! pour avoir archim=10 
    188                 new_flot_point(i,j)=.false. 
    189              endif 
    190  
    191  
    192  
    193  
    194          else                                       ! isynchro=0 ou il flottait déja 
    195                        
    196             if (.not.FLOT(I,J)) then               ! il ne flottait pas (isynchro=0) 
    197                new_flot_point(i,j)=.true.          ! signale un point qui ne flottait pas 
    198                                                    ! au pas de temps precedent 
    199                        flot(i,j)=.true. 
    200   
     159 
     160    ! fin de l'initialisation 
     161    !_____________________________________________________________________ 
     162 
     163    ! 2-TESTE LES NOUVEAUX POINTS FLOTTANTS 
     164    ! ------------------------------------- 
     165 
     166    !$OMP DO 
     167    do j=1,ny 
     168       do i=1,nx 
     169 
     170          uxs1(i,j)=uxbar(i,j) 
     171          uys1(i,j)=uybar(i,j) 
     172 
     173          archim = Bsoc(i,j)+H(i,j)*ro/row -sealevel 
     174          !      if ((i.eq.132).and.(j.eq.183)) print*,'archim=',archim 
     175 
     176 
     177          arch: if ((ARCHIM.LT.0.).and.(H(I,J).gt.1.e-3)) then    ! le point flotte 
     178             mk(i,j)=1 
     179 
     180 
     181             ex_pose: if ((.not.FLOT(I,J)).and.(isynchro.eq.1)) then  !  il ne flottait pas avant 
     182                FLOT(I,J)=.true. 
     183                BOOST=.false. 
     184 
    201185                if (igrdline.eq.1) then   ! en cas de grounding line prescrite 
    202                     flot(i,j)=.false. 
    203                     H(i,j)=(10.+sealevel-Bsoc(i,j))*row/ro  ! pour avoir archim=10 
     186                   flot(i,j)=.false. 
     187                   H(i,j)=(10.+sealevel-Bsoc(i,j))*row/ro  ! pour avoir archim=10 
    204188                   new_flot_point(i,j)=.false. 
    205                endif 
    206  
    207             endif 
    208          endif ex_pose 
    209  
    210  
    211        else                           !    le point ne flotte pas 
    212 mk(i,j)=0 
    213  
    214            if(FLOT(I,J)) then         !  mais il flottait avant 
     189                endif 
     190 
     191             else                                       ! isynchro=0 ou il flottait déja 
     192 
     193                if (.not.FLOT(I,J)) then               ! il ne flottait pas (isynchro=0) 
     194                   new_flot_point(i,j)=.true.          ! signale un point qui ne flottait pas 
     195                   ! au pas de temps precedent 
     196                   flot(i,j)=.true. 
     197 
     198                   if (igrdline.eq.1) then   ! en cas de grounding line prescrite 
     199                      flot(i,j)=.false. 
     200                      H(i,j)=(10.+sealevel-Bsoc(i,j))*row/ro  ! pour avoir archim=10 
     201                      new_flot_point(i,j)=.false. 
     202                   endif 
     203 
     204                endif 
     205             endif ex_pose 
     206 
     207 
     208          else if ((H(i,j).gt.0.).and.(archim.GE.0.)) then    !    le point ne flotte pas et est englace 
     209             mk(i,j)=0 
     210 
     211             if(FLOT(I,J)) then         !  mais il flottait avant 
    215212                FLOT(I,J)=.false. 
    216213                BOOST=.false. 
    217            endif 
    218        endif arch 
    219  
    220 !   Si la glace flotte -> PRUDENCE !!! 
    221 !   S et B sont alors determines avec la condition de flottabilite 
    222  
    223        if (flot(i,j)) then 
    224           shelfy = .true. 
    225            
    226           surnet=H(i,j)*(1.-ro/row) 
    227           S(i,j)=surnet+sealevel 
    228           B(i,j)=S(i,j)-H(i,j) 
    229        end if 
    230  
     214 
     215             endif 
     216             !cdc correction topo pour suivre  variations sealevel 
     217             !cdd           S(i,j)=Bsoc(i,j)+H(i,j) 
     218             S(i,j)=Bsoc(i,j)+H(i,j) 
     219             B(i,j)=Bsoc(i,j) 
     220 
     221 
     222          else if  ((H(i,j).LE.0.).and.(archim.LT.0.)) then    !    terre deglace qui devient ocean  
     223             ice(i,j)=0 
     224             h(i,j)=1. 
     225             surnet=H(i,j)*(1.-ro/row) 
     226             S(i,j)=surnet+sealevel 
     227             B(i,j)=S(i,j)-H(i,j) 
     228          endif arch 
     229 
     230          !   Si la glace flotte -> PRUDENCE !!! 
     231          !   S et B sont alors determines avec la condition de flottabilite 
     232 
     233          if (flot(i,j)) then 
     234             shelfy = .true. 
     235 
     236             surnet=H(i,j)*(1.-ro/row) 
     237             S(i,j)=surnet+sealevel 
     238             B(i,j)=S(i,j)-H(i,j) 
     239          end if 
     240 
     241       end do 
    231242    end do 
    232  end do 
    233  !$OMP END DO 
     243    !$OMP END DO 
     244    !~ print*,'flottab 2',S(132,183),H(132,183),BSOC(132,183),B(132,183),sealevel 
     245    !~ print*,'flottab 2',flot(132,183),ice(132,183) 
     246    !~ if (S(132,183).LT.sealevel) print*,'BUGGGGGGGGGGGGG !!!!!!!!!!!!!!!!' 
    234247 
    235248!!$ do i=1,nx 
    … …  
    247260 
    248261 
    249 !----------------------------------------------------------------------- 
    250  !$OMP DO 
    251 domain_x: do j=1,ny        
    252          do i=2,nx 
    253  
    254 !    3_x    A- NOUVELLE DEFINITION DE FLOTMX, LES POINTS 
    255 !            AYANT UN DES VOISINS FLOTTANTS SONT FLOTMX ET GZMX 
    256 !         ------------------------------------------- 
    257  
    258 !            fl1 est l'ancienne valeur de flotmx 
    259              gz1mx(i,j)=gzmx(i,j)         !     gz1 est l'ancienne valeur de gzmx 
    260              fl1mx(i,j)=flotmx(i,j)       !     fl1 est l'ancienne valeur de flotmx 
    261  
    262              flotmx(i,j)=flot(i,j).and.flot(i-1,j) 
    263  
    264 ! test pour detecter les nouveaux flotmx entre 2 dtt : 
    265  
    266              if (flotmx(i,j).and.(new_flot_point(i,j).or. & 
    267                 new_flot_point(i-1,j))) then 
    268                 appel_new_flot=.true. 
    269                 new_flotmx(i,j)=.true. 
    270              endif 
    271  
    272  
    273 !   premiere determination de gzmx 
    274 !__________________________________________________________________________ 
    275  
    276 ! gzmx si un des deux voisins est flottant et l'autre posé 
    277 !                                                                   i-1   i 
    278              gzmx(i,j)=((flot(i,j).and..not.flot(i-1,j)) & !         F    P 
    279              .or.(.not.flot(i,j).and.flot(i-1,j)))         !         P    F 
    280  
    281 !                 A condition d'etre assez proche de la flottaison 
    282 !                 sur le demi noeud condition archim < 100 m 
    283  
    284              archim=(Bsoc(i,j)+Bsoc(i-1,j))*0.5-sealevel+ro/row*Hmx(i,j) 
    285              gzmx(i,j)=gzmx(i,j).and.(archim.le.100.)  
    286              cotemx(i,j)=gzmx(i,j) 
    287  
    288           end do 
    289        end do domain_x 
    290  !$OMP END DO 
    291 !if (itracebug.eq.1)  call tracebug('  routine flottab apres domain_x') 
    292  
    293 !     3_y    B- NOUVELLE DEFINITION DE FLOTMY 
    294 !         -------------------------------- 
    295  !$OMP DO 
    296 domain_y: do j=2,ny        
    297          do i=1,nx 
    298  
    299              gz1my(i,j)=gzmy(i,j)           !       gz1 est l'ancienne valeur de gzmy 
    300              fl1my(i,j)=flotmy(i,j)         !       fl1 est l'ancienne valeur de flotmy 
    301  
    302              flotmy(i,j)=flot(i,j).and.flot(i,j-1) 
    303  
    304 ! test pour detecter les nouveaux flotmy entre 2 dtt : 
    305  
    306              if (flotmy(i,j).and.(new_flot_point(i,j).or.     & 
    307              new_flot_point(i,j-1))) then 
    308                 appel_new_flot=.true. 
    309                 new_flotmy(i,j)=.true. 
    310              endif 
    311  
    312 !   premiere determination de gzmy 
    313 !__________________________________________________________________________ 
    314  
    315 ! gzmy si un des deux voisins est flottant et l'autre posé 
    316  
    317              gzmy(i,j)=((flot(i,j).and..not.flot(i,j-1))          & 
    318              .or.(.not.flot(i,j).and.flot(i,j-1)))                  
    319  
    320 !                 A condition d'etre assez proche de la flottaison 
    321 !                 sur le demi noeud condition archim > 100 m 
    322  
    323              archim=(Bsoc(i,j)+Bsoc(i,j-1))*0.5-sealevel+ro/row*Hmy(i,j) 
    324              gzmy(i,j)=gzmy(i,j).and.(archim.le.100.)  
    325              cotemy(i,j)=gzmy(i,j) 
    326  
    327           end do 
    328        end do domain_y 
    329  !$OMP END DO 
    330  
    331  
    332 !------------------------------------------------------------------------------------- 
    333 ! attention : pour expériences Heino 
    334 !             gzmy(i,j)=gzmy_heino(i,j) 
    335 !             shelfy=.true. 
    336 !             shelfy=.false. 
    337 !____________________________________________________________________________________  
     262    !----------------------------------------------------------------------- 
     263    !$OMP DO 
     264    domain_x: do j=1,ny        
     265       do i=2,nx 
     266 
     267          !    3_x    A- NOUVELLE DEFINITION DE FLOTMX, LES POINTS 
     268          !            AYANT UN DES VOISINS FLOTTANTS SONT FLOTMX ET GZMX 
     269          !         ------------------------------------------- 
     270 
     271          !            fl1 est l'ancienne valeur de flotmx 
     272          gz1mx(i,j)=gzmx(i,j)         !     gz1 est l'ancienne valeur de gzmx 
     273          fl1mx(i,j)=flotmx(i,j)       !     fl1 est l'ancienne valeur de flotmx 
     274 
     275          flotmx(i,j)=flot(i,j).and.flot(i-1,j) 
     276 
     277          ! test pour detecter les nouveaux flotmx entre 2 dtt : 
     278 
     279          if (flotmx(i,j).and.(new_flot_point(i,j).or. & 
     280               new_flot_point(i-1,j))) then 
     281             appel_new_flot=.true. 
     282             new_flotmx(i,j)=.true. 
     283          endif 
     284 
     285 
     286          !   premiere determination de gzmx 
     287          !__________________________________________________________________________ 
     288 
     289          ! gzmx si un des deux voisins est flottant et l'autre posé 
     290          !                                                                   i-1   i 
     291          gzmx(i,j)=((flot(i,j).and..not.flot(i-1,j)) & !         F    P 
     292               .or.(.not.flot(i,j).and.flot(i-1,j)))         !         P    F 
     293 
     294          !                 A condition d'etre assez proche de la flottaison 
     295          !                 sur le demi noeud condition archim < 100 m 
     296 
     297          archim=(Bsoc(i,j)+Bsoc(i-1,j))*0.5-sealevel+ro/row*Hmx(i,j) 
     298          gzmx(i,j)=gzmx(i,j).and.(archim.le.100.)  
     299          cotemx(i,j)=gzmx(i,j) 
     300 
     301       end do 
     302    end do domain_x 
     303    !$OMP END DO 
     304    !if (itracebug.eq.1)  call tracebug('  routine flottab apres domain_x') 
     305 
     306    !     3_y    B- NOUVELLE DEFINITION DE FLOTMY 
     307    !         -------------------------------- 
     308    !$OMP DO 
     309    domain_y: do j=2,ny        
     310       do i=1,nx 
     311 
     312          gz1my(i,j)=gzmy(i,j)           !       gz1 est l'ancienne valeur de gzmy 
     313          fl1my(i,j)=flotmy(i,j)         !       fl1 est l'ancienne valeur de flotmy 
     314 
     315          flotmy(i,j)=flot(i,j).and.flot(i,j-1) 
     316 
     317          ! test pour detecter les nouveaux flotmy entre 2 dtt : 
     318 
     319          if (flotmy(i,j).and.(new_flot_point(i,j).or.     & 
     320               new_flot_point(i,j-1))) then 
     321             appel_new_flot=.true. 
     322             new_flotmy(i,j)=.true. 
     323          endif 
     324 
     325          !   premiere determination de gzmy 
     326          !__________________________________________________________________________ 
     327 
     328          ! gzmy si un des deux voisins est flottant et l'autre posé 
     329 
     330          gzmy(i,j)=((flot(i,j).and..not.flot(i,j-1))          & 
     331               .or.(.not.flot(i,j).and.flot(i,j-1)))                  
     332 
     333          !                 A condition d'etre assez proche de la flottaison 
     334          !                 sur le demi noeud condition archim > 100 m 
     335 
     336          archim=(Bsoc(i,j)+Bsoc(i,j-1))*0.5-sealevel+ro/row*Hmy(i,j) 
     337          gzmy(i,j)=gzmy(i,j).and.(archim.le.100.)  
     338          cotemy(i,j)=gzmy(i,j) 
     339 
     340       end do 
     341    end do domain_y 
     342    !$OMP END DO 
     343 
     344 
     345    !------------------------------------------------------------------------------------- 
     346    ! attention : pour expériences Heino 
     347    !             gzmy(i,j)=gzmy_heino(i,j) 
     348    !             shelfy=.true. 
     349    !             shelfy=.false. 
     350    !____________________________________________________________________________________  
    338351 
    339352 
    … …  
    358371 
    359372 
    360 !      4- determination des iles 
    361 !      ------------------------- 
    362           !$OMP WORKSHARE 
    363           ilemx(:,:)=.false. 
    364           ilemy(:,:)=.false. 
    365           !$OMP END WORKSHARE 
    366  
    367 !       selon x 
    368   !$OMP DO 
    369 ilesx:  do j=2,ny-1 
    370            do i=3,nx-2 
    371 !                       F   G   F    
    372 !                           x 
    373 ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m) 
    374             if ((flot(i-1,j).and..not.flot(i,j).and.flot(i+1,j)).and. & 
    375                 (sdx(i,j).LT.1.E-02)) then 
     373    !      4- determination des iles 
     374    !      ------------------------- 
     375    !$OMP WORKSHARE 
     376    ilemx(:,:)=.false. 
     377    ilemy(:,:)=.false. 
     378    !$OMP END WORKSHARE 
     379 
     380    !       selon x 
     381    !$OMP DO 
     382    ilesx:  do j=2,ny-1 
     383       do i=3,nx-2 
     384          !                       F   G   F    
     385          !                           x 
     386          ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m) 
     387          if ((flot(i-1,j).and..not.flot(i,j).and.flot(i+1,j)).and. & 
     388               (sdx(i,j).LT.1.E-02)) then 
    376389             ilemx(i,j)=.true. 
    377390             ilemx(i+1,j)=.true. 
    378391 
    379 !                      F  G   G   F 
    380 !                         x 
    381 ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m) 
    382             else if ((flot(i-1,j).and..not.flot(i,j)                  & 
    383                 .and..not.flot(i+1,j)).and.flot(i+2,j).and.           & 
    384                  (sdx(i,j).LT.1.E-02.and.sdx(i+1,j).LT.1.E-02)) then 
    385               ilemx(i,j)=.true. 
    386               ilemx(i+1,j)=.true. 
    387               ilemx(i+2,j)=.true. 
    388  
    389 !                      F   G   G   F 
    390 !                              x 
    391 ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m) 
    392             else if ((flot(i-2,j).and..not.flot(i-1,j)                & 
    393                 .and..not.flot(i,j)).and.flot(i+1,j).and.             & 
    394                  (sdx(i,j).LT.1.E-02.and.sdx(i-1,j).LT.1.E-02)) then 
    395               ilemx(i-1,j)=.true. 
    396               ilemx(i,j)=.true. 
    397               ilemx(i+1,j)=.true. 
    398  
    399 !                      F   G   G   G    F 
    400 !                              x 
    401 ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m) 
    402             else if ((i.lt.nx-2)                                      & 
    403                 .and.(flot(i-2,j).and..not.flot(i-1,j)                & 
    404                 .and..not.flot(i,j)).and..not.flot(i+1,j)             & 
    405                  .and.flot(i+2,j).and.                                & 
    406                  (sdx(i,j).LT.1.E-02.and.sdx(i-1,j).LT.1.E-02         & 
    407                  .and.sdx(i+1,j).LT.1.E-02)) then 
    408               ilemx(i-1,j)=.true. 
    409               ilemx(i,j)=.true. 
    410               ilemx(i+1,j)=.true. 
    411               ilemx(i+2,j)=.true. 
    412  
    413             endif 
    414  
    415          end do 
    416       end do ilesx 
    417    !$OMP END DO 
    418  
    419 !       selon y 
    420    !$OMP DO 
    421 ilesy: do j=3,ny-2 
    422            do i=2,nx-1 
    423 !                       F   G   F    
    424 !                           x 
    425 ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m) 
    426             if ((flot(i,j-1).and..not.flot(i,j).and.flot(i,j+1)).and. & 
    427                  (sdy(i,j).LT.1.E-02)) then 
     392             !                      F  G   G   F 
     393             !                         x 
     394             ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m) 
     395          else if ((flot(i-1,j).and..not.flot(i,j)                  & 
     396               .and..not.flot(i+1,j)).and.flot(i+2,j).and.           & 
     397               (sdx(i,j).LT.1.E-02.and.sdx(i+1,j).LT.1.E-02)) then 
     398             ilemx(i,j)=.true. 
     399             ilemx(i+1,j)=.true. 
     400             ilemx(i+2,j)=.true. 
     401 
     402             !                      F   G   G   F 
     403             !                              x 
     404             ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m) 
     405          else if ((flot(i-2,j).and..not.flot(i-1,j)                & 
     406               .and..not.flot(i,j)).and.flot(i+1,j).and.             & 
     407               (sdx(i,j).LT.1.E-02.and.sdx(i-1,j).LT.1.E-02)) then 
     408             ilemx(i-1,j)=.true. 
     409             ilemx(i,j)=.true. 
     410             ilemx(i+1,j)=.true. 
     411 
     412             !                      F   G   G   G    F 
     413             !                              x 
     414             ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m) 
     415          else if ((i.lt.nx-2)                                      & 
     416               .and.(flot(i-2,j).and..not.flot(i-1,j)                & 
     417               .and..not.flot(i,j)).and..not.flot(i+1,j)             & 
     418               .and.flot(i+2,j).and.                                & 
     419               (sdx(i,j).LT.1.E-02.and.sdx(i-1,j).LT.1.E-02         & 
     420               .and.sdx(i+1,j).LT.1.E-02)) then 
     421             ilemx(i-1,j)=.true. 
     422             ilemx(i,j)=.true. 
     423             ilemx(i+1,j)=.true. 
     424             ilemx(i+2,j)=.true. 
     425 
     426          endif 
     427 
     428       end do 
     429    end do ilesx 
     430    !$OMP END DO 
     431 
     432    !       selon y 
     433    !$OMP DO 
     434    ilesy: do j=3,ny-2 
     435       do i=2,nx-1 
     436          !                       F   G   F    
     437          !                           x 
     438          ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m) 
     439          if ((flot(i,j-1).and..not.flot(i,j).and.flot(i,j+1)).and. & 
     440               (sdy(i,j).LT.1.E-02)) then 
    428441             ilemy(i,j)=.true. 
    429442             ilemy(i,j+1)=.true. 
    430443 
    431 !                      F  G   G   F 
    432 !                         x 
    433 ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m) 
    434             else if ((flot(i,j-1).and..not.flot(i,j)                  & 
    435                 .and..not.flot(i,j+1)).and.flot(i,j+2).and.           & 
    436                  (sdy(i,j).LT.1.E-02.and.sdy(i,j+1).LT.1.E-02)) then 
    437               ilemy(i,j)=.true. 
    438               ilemy(i,j+1)=.true. 
    439               ilemy(i,j+2)=.true. 
    440  
    441 !                      F   G   G   F 
    442 !                              x 
    443 ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m) 
    444             else if ((flot(i,j-2).and..not.flot(i,j-1)                 & 
    445                 .and..not.flot(i,j)).and.flot(i,j+1).and.              & 
    446                  (sdy(i,j).LT.1.E-02.and.sdy(i,j-1).LT.1.E-02)) then 
    447               ilemy(i,j-1)=.true. 
    448               ilemy(i,j)=.true. 
    449               ilemy(i,j+1)=.true. 
    450  
    451 !                      F   G   G   G    F 
    452 !                              x 
    453 ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m) 
    454             else if ((j.lt.ny-2)                                      & 
    455                 .and.(flot(i,j-2).and..not.flot(i,j-1)                & 
    456                 .and..not.flot(i,j)).and..not.flot(i,j+1)             & 
    457                  .and.flot(i,j+2).and.                                & 
    458                  (sdy(i,j).LT.1.E-02.and.sdy(i,j-1).LT.1.E-02         & 
    459                  .and.sdy(i,j+1).LT.1.E-02)) then  
    460               ilemy(i,j-1)=.true. 
    461               ilemy(i,j)=.true. 
    462               ilemy(i,j+1)=.true. 
    463               ilemy(i,j+2)=.true. 
    464             endif 
    465          end do 
    466       end do ilesy 
     444             !                      F  G   G   F 
     445             !                         x 
     446             ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m) 
     447          else if ((flot(i,j-1).and..not.flot(i,j)                  & 
     448               .and..not.flot(i,j+1)).and.flot(i,j+2).and.           & 
     449               (sdy(i,j).LT.1.E-02.and.sdy(i,j+1).LT.1.E-02)) then 
     450             ilemy(i,j)=.true. 
     451             ilemy(i,j+1)=.true. 
     452             ilemy(i,j+2)=.true. 
     453 
     454             !                      F   G   G   F 
     455             !                              x 
     456             ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m) 
     457          else if ((flot(i,j-2).and..not.flot(i,j-1)                 & 
     458               .and..not.flot(i,j)).and.flot(i,j+1).and.              & 
     459               (sdy(i,j).LT.1.E-02.and.sdy(i,j-1).LT.1.E-02)) then 
     460             ilemy(i,j-1)=.true. 
     461             ilemy(i,j)=.true. 
     462             ilemy(i,j+1)=.true. 
     463 
     464             !                      F   G   G   G    F 
     465             !                              x 
     466             ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m) 
     467          else if ((j.lt.ny-2)                                      & 
     468               .and.(flot(i,j-2).and..not.flot(i,j-1)                & 
     469               .and..not.flot(i,j)).and..not.flot(i,j+1)             & 
     470               .and.flot(i,j+2).and.                                & 
     471               (sdy(i,j).LT.1.E-02.and.sdy(i,j-1).LT.1.E-02         & 
     472               .and.sdy(i,j+1).LT.1.E-02)) then  
     473             ilemy(i,j-1)=.true. 
     474             ilemy(i,j)=.true. 
     475             ilemy(i,j+1)=.true. 
     476             ilemy(i,j+2)=.true. 
     477          endif 
     478       end do 
     479    end do ilesy 
    467480    !$OMP END DO 
    468481    !$OMP END PARALLEL 
    469 ! fin des iles 
     482    ! fin des iles 
    470483 
    471484!!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf) 
    … …  
    478491!!$end if 
    479492 
    480 ! 5- calcule les noeuds qui sont streams a l'interieur et donne le betamx et betamy 
    481 !---------------------------------------------------------------------------------- 
    482 if (iter_beta.eq.0) then 
    483  Call dragging           
    484 endif 
    485 if (itracebug.eq.1)  call tracebug(' routine flottab apres call dragging') 
     493    ! 5- calcule les noeuds qui sont streams a l'interieur et donne le betamx et betamy 
     494    !---------------------------------------------------------------------------------- 
     495    if (iter_beta.eq.0) then 
     496       Call dragging           
     497    endif 
     498    if (itracebug.eq.1)  call tracebug(' routine flottab apres call dragging') 
    486499!!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf) 
    487500!!$if (itestf.gt.0) then 
    … …  
    493506!!$end if 
    494507 
    495 !   6- calcule les vitesses des points qui sont devenus gzm 
    496 !$OMP PARALLEL 
    497 !$OMP DO          
    498 do j=1,ny 
    499    do i=2,nx-1 
    500 !      si le point etait posé (non gz) et devient gzmx 
    501 !      definir la direction de la vitesse (moyenne des points) 
    502  
    503       if ((.not.gz1mx(i,j)).and.(.not.fl1mx(i,j)).and.gzmx(i,j).and. & 
    504        (i.gt.2).and.(i.lt.nx))  then 
    505          uxs1(i,j)=(uxbar(i+1,j)+uxbar(i-1,j))/2. 
    506       endif 
    507  
    508    end do 
    509 end do 
    510 !$OMP END DO 
    511  
    512 !$OMP DO  
    513 do j=2,ny-1 
    514    do i=1,nx 
    515 !      si le point etait posé (non gz) et devient gzmy 
    516 !      definir la direction de la vitesse (moyenne des points) 
    517       if ((.not.gz1my(i,j)).and.(.not.fl1my(i,j)).and.gzmy(i,j).and. & 
    518        (j.gt.2).and.(j.lt.ny))  then 
    519           uys1(i,j)=(uybar(i,j+1)+uybar(i,j-1))/2. 
    520        endif 
    521     
     508    !   6- calcule les vitesses des points qui sont devenus gzm 
     509    !$OMP PARALLEL 
     510    !$OMP DO          
     511    do j=1,ny 
     512       do i=2,nx-1 
     513          !      si le point etait posé (non gz) et devient gzmx 
     514          !      definir la direction de la vitesse (moyenne des points) 
     515 
     516          if ((.not.gz1mx(i,j)).and.(.not.fl1mx(i,j)).and.gzmx(i,j).and. & 
     517               (i.gt.2).and.(i.lt.nx))  then 
     518             uxs1(i,j)=(uxbar(i+1,j)+uxbar(i-1,j))/2. 
     519          endif 
     520 
     521       end do 
    522522    end do 
    523 end do 
    524 !$OMP END DO 
    525  
    526  
    527 !  7-On determine finalement la position des noeuds stream ou shelf 
    528 !  ------------------------------------------------------------- 
    529  
    530 if (nt.ge.2) then   ! pour ne pas faire ce calcul lors du premier passage 
    531    !$OMP WORKSHARE 
    532    uxbar(:,:)=uxs1(:,:) 
    533    uybar(:,:)=uys1(:,:) 
    534    !$OMP END WORKSHARE 
    535 endif 
    536  
    537 !$OMP WORKSHARE 
    538 flgzmx(:,:)=(marine.and.(flotmx(:,:).or.gzmx(:,:).or.ilemx(:,:)))   & 
    539      .or.(.not.marine.and.flotmx(:,:)) 
    540 flgzmy(:,:)=(marine.and.(flotmy(:,:).or.gzmy(:,:).or.ilemy(:,:)))   & 
    541      .or.(.not.marine.and.flotmy(:,:)) 
    542 !$OMP END WORKSHARE 
    543  
    544  
    545 ! 8- Pour la fusion basale sous les ice shelves- region proche de la grounding line  
    546 !--------------------------------------------------------------------------------- 
    547 !       fbm est vrai si le point est flottant mais un des voisins est pose  
    548 !_________________________________________________________________________ 
    549 !$OMP DO 
    550 do j=2,ny-1 
    551   do i=2,nx-1 
    552  
    553 ! if (i.gt.2.AND.i.lt.nx) then    
    554         fbm(i,j)=flot(i,j).and.      &  
    555              ((.not.flot(i+1,j)).or.(.not.flot(i,j+1)) & 
    556              .or.(.not.flot(i-1,j)).or.(.not.flot(i,j-1))) 
    557 !     endif 
    558   end do 
    559 end do 
    560 !$OMP END DO 
    561  
    562  
    563 !  9-On determine maintenant la position du front de glace 
    564 !  ------------------------------------------------------- 
    565 !  C'est ici que l'on determine la position et le type de front 
    566 !       print*,'on est dans flottab pour definir les fronts' 
    567        
     523    !$OMP END DO 
     524 
     525    !$OMP DO  
     526    do j=2,ny-1 
     527       do i=1,nx 
     528          !      si le point etait posé (non gz) et devient gzmy 
     529          !      definir la direction de la vitesse (moyenne des points) 
     530          if ((.not.gz1my(i,j)).and.(.not.fl1my(i,j)).and.gzmy(i,j).and. & 
     531               (j.gt.2).and.(j.lt.ny))  then 
     532             uys1(i,j)=(uybar(i,j+1)+uybar(i,j-1))/2. 
     533          endif 
     534 
     535       end do 
     536    end do 
     537    !$OMP END DO 
     538 
     539 
     540    !  7-On determine finalement la position des noeuds stream ou shelf 
     541    !  ------------------------------------------------------------- 
     542 
     543    if (nt.ge.2) then   ! pour ne pas faire ce calcul lors du premier passage 
     544       !$OMP WORKSHARE 
     545       uxbar(:,:)=uxs1(:,:) 
     546       uybar(:,:)=uys1(:,:) 
     547       !$OMP END WORKSHARE 
     548    endif 
     549 
     550    !$OMP WORKSHARE 
     551    flgzmx(:,:)=(marine.and.(flotmx(:,:).or.gzmx(:,:).or.ilemx(:,:)))   & 
     552         .or.(.not.marine.and.flotmx(:,:)) 
     553    flgzmy(:,:)=(marine.and.(flotmy(:,:).or.gzmy(:,:).or.ilemy(:,:)))   & 
     554         .or.(.not.marine.and.flotmy(:,:)) 
     555    !$OMP END WORKSHARE 
     556 
     557 
     558    ! 8- Pour la fusion basale sous les ice shelves- region proche de la grounding line  
     559    !--------------------------------------------------------------------------------- 
     560    !       fbm est vrai si le point est flottant mais un des voisins est pose  
     561    !_________________________________________________________________________ 
     562    !$OMP DO 
     563    do j=2,ny-1 
     564       do i=2,nx-1 
     565 
     566          ! if (i.gt.2.AND.i.lt.nx) then    
     567          fbm(i,j)=flot(i,j).and.      &  
     568               ((.not.flot(i+1,j)).or.(.not.flot(i,j+1)) & 
     569               .or.(.not.flot(i-1,j)).or.(.not.flot(i,j-1))) 
     570          !     endif 
     571       end do 
     572    end do 
     573    !$OMP END DO 
     574 
     575 
     576    !  9-On determine maintenant la position du front de glace 
     577    !  ------------------------------------------------------- 
     578    !  C'est ici que l'on determine la position et le type de front 
     579    !       print*,'on est dans flottab pour definir les fronts' 
     580 
    568581 
    569582 
    … …  
    574587!!$end do 
    575588 
    576 !$OMP WORKSHARE 
    577 where (flot(:,:)) 
    578    where (H(:,:).gt.(1.1))  
    579       ice(:,:)=1 
    580    elsewhere  
    581       ice(:,:)=0 
    582    end where 
    583 elsewhere 
    584   where (H(:,:).gt.(.1))  
    585       ice(:,:)=1 
    586    elsewhere  
    587       ice(:,:)=0 
    588    end where 
    589 end where 
    590 !$OMP END WORKSHARE 
    591 !$OMP END PARALLEL 
    592  
    593 call DETERMIN_TACHE  
    594  
    595 synchro :    if (isynchro.eq.1) then 
    596 !!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf) 
    597 !!$if (itestf.gt.0) then 
    598 !!$   write(6,*) 'dans flottab avant DETERMIN_TACHE  asymetrie sur H  pour time=',time 
    599 !!$   stop 
    600 !!$else 
    601 !!$   write(6,*) 'dans flottab apres  DETERMIN_TACHE pas d asymetrie sur H  pour time=',time 
    602 !!$ 
    603 !!$end if 
    604  
    605  
    606 !----------------------------------------------! 
    607 !On determine les differents ice strean/shelf  ! 
    608 !      call DETERMIN_TACHE                      ! 
    609 !----------------------------------------------! 
    610  
    611  
    612 !!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf) 
    613 !!$if (itestf.gt.0) then 
    614 !!$   write(6,*) 'dans flottab apres DETERMIN_TACHE  asymetrie sur H  pour time=',time 
    615 !!$   stop 
    616 !!$else 
    617 !!$   write(6,*) 'dans flottab apres DETERMIN_TACHE pas d asymetrie sur H  pour time=',time 
    618 !!$ 
    619 !!$end if 
    620  
    621 !ice(:,:)=0   Attention, voir si ca marche toujours pour l'Antarctique et heminord ! 
    622  
    623 !On compte comme englacé uniquement les calottes dont une partie est posée       
    624 !$OMP PARALLEL PRIVATE(smax_,smax_coord,smax_i,smax_j,mask_tache_ij) 
    625 !$OMP DO 
    626 do j=3,ny-2 
    627    do i=3,nx-2 
    628 test1:  if (.not.iceberg(table_out(i,j))) then ! on est pas sur un iceberg  
    629          if (nb_pts_tache(table_out(i,j)).ge.1) then 
    630             ice(i,j)=1 
    631             if (nb_pts_tache(table_out(i,j)).le.10) then  ! les petites iles sont en sia 
    632 !    write(6,*) 'petite ile ',i,j 
    633               flgzmx(i,j)=.false. 
    634                flgzmx(i+1,j)=.false. 
    635                flgzmy(i,j)=.false. 
    636                flgzmy(i,j+1)=.false. 
    637                gzmx(i:i+1,j)=.false. 
    638                gzmy(i,j:j+1)=.false. 
    639             endif 
    640  
    641  
    642 ! ici on est sur une tache non iceberg >= 5 points                        
    643 ! on teste si la tache n'est pas completement ice stream 
    644  
    645 test2:       if (icetrim(table_out(i,j))) then   ! on a une ile d'ice stream 
    646  
    647    mask_tache_ij(:,:)=.false. 
    648    mask_tache_ij(:,:)=(table_out(:,:).eq.table_out(i,j))       ! pour toute la tache 
    649  
    650      smax_=maxval(S(:,:),MASK=mask_tache_ij(:,:)) 
    651      smax_coord(:)=maxloc(S(:,:),MASK=mask_tache_ij(:,:)) 
    652      smax_i=smax_coord(1) 
    653      smax_j=smax_coord(2) 
    654  
    655 !!$               smax_i=0 ; smax_j=0 ; smax_=sealevel 
    656 !!$               do ii=3,nx-2 
    657 !!$                  do jj=3,ny-2 
    658 !!$                     if (table_out(ii,jj).eq.table_out(i,j)) then 
    659 !!$                        if (s(ii,jj).gt.smax_) then                            
    660 !!$                           smax_ =s(ii,jj) 
    661 !!$                           smax_i=ii 
    662 !!$                           smax_j=jj 
    663 !!$                        endif 
    664 !!$                     endif 
    665 !!$                  end do 
    666 !!$               end do 
    667  
    668  
    669                gzmx(smax_i,smax_j)=.false. ; gzmx(smax_i+1,smax_j)=.false. 
    670                gzmy(smax_i,smax_j)=.false. ; gzmx(smax_i,smax_j+1)=.false. 
    671                flgzmx(smax_i,smax_j)=.false. ; flgzmx(smax_i+1,smax_j)=.false. 
    672                flgzmy(smax_i,smax_j)=.false. ; flgzmx(smax_i,smax_j+1)=.false. 
    673   
    674                if (Smax_.le.sealevel) then 
    675                   write(num_tracebug,*)'Attention, une ile avec la surface sous l eau' 
    676                   write(num_tracebug,*)'time=',time,'   coord:',smax_i,smax_j 
    677                end if 
    678  
    679             endif test2  
    680          end if                                             ! endif deplace 
    681           
    682       else ! on est sur un iceberg                          !   test1 
    683          ice(i,j)=0 
    684          h(i,j)=1. 
    685          surnet=H(i,j)*(1.-ro/row) 
    686          S(i,j)=surnet+sealevel 
    687          B(i,j)=S(i,j)-H(i,j) 
    688           
    689       endif test1 
    690  
    691  
    692    end do 
    693 end do 
    694 !$OMP END DO 
    695 !$OMP END PARALLEL 
    696  
    697 !---------------------------------------------- 
    698 ! On caracterise le front des ice shelfs/streams 
    699  
    700 !      call DETERMIN_FRONT     
    701  
    702 !----------------------------------------------       
    703 !!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf) 
    704 !!$if (itestf.gt.0) then 
    705 !!$   write(6,*) 'dans flottab apres DETERMIN_front  asymetrie sur H  pour time=',time 
    706 !!$   stop 
    707 !!$else 
    708 !!$   write(6,*) 'dans flottab apres DETERMIN_front pas d asymetrie sur H  pour time=',time 
    709 !!$ 
    710 !!$end if 
    711  
    712 endif synchro 
    713  
    714 ! correction momentanée pour symetrie Heino 
    715 !where ((.not.flot(:,:)).and.(ice(:,:).eq.0)) H(:,:)=0.       
    716  
    717 !fin de routine flottab2 
    718 !print*, 'front',front(50,30),ice(50,30),flotmx(i,j),uxbar(i,j) 
    719  
    720 end subroutine flottab 
     589    !$OMP WORKSHARE 
     590    where (flot(:,:)) 
     591       where (H(:,:).gt.(1.1))  
     592          ice(:,:)=1 
     593       elsewhere  
     594          ice(:,:)=0 
     595       end where 
     596    elsewhere 
     597       where (H(:,:).gt.(.1))  
     598          ice(:,:)=1 
     599       elsewhere  
     600          ice(:,:)=0 
     601       end where 
     602    end where 
     603    !$OMP END WORKSHARE 
     604    !$OMP END PARALLEL 
     605 
     606    !~ call DETERMIN_TACHE  
     607    !~  
     608    !~ synchro :    if (isynchro.eq.1) then 
     609    !~ !!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf) 
     610    !~ !!$if (itestf.gt.0) then 
     611    !~ !!$   write(6,*) 'dans flottab avant DETERMIN_TACHE  asymetrie sur H  pour time=',time 
     612    !~ !!$   stop 
     613    !~ !!$else 
     614    !~ !!$   write(6,*) 'dans flottab apres  DETERMIN_TACHE pas d asymetrie sur H  pour time=',time 
     615    !~ !!$ 
     616    !~ !!$end if 
     617    !~  
     618    !~  
     619    !~ !----------------------------------------------! 
     620    !~ !On determine les differents ice strean/shelf  ! 
     621    !~ !      call DETERMIN_TACHE                      ! 
     622    !~ !----------------------------------------------! 
     623    !~  
     624    !~  
     625    !~ !!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf) 
     626    !~ !!$if (itestf.gt.0) then 
     627    !~ !!$   write(6,*) 'dans flottab apres DETERMIN_TACHE  asymetrie sur H  pour time=',time 
     628    !~ !!$   stop 
     629    !~ !!$else 
     630    !~ !!$   write(6,*) 'dans flottab apres DETERMIN_TACHE pas d asymetrie sur H  pour time=',time 
     631    !~ !!$ 
     632    !~ !!$end if 
     633    !~  
     634    !~ !ice(:,:)=0   Attention, voir si ca marche toujours pour l'Antarctique et heminord ! 
     635    !~  
     636    !~ !On compte comme englacé uniquement les calottes dont une partie est posée       
     637    !~ !$OMP PARALLEL PRIVATE(smax_,smax_coord,smax_i,smax_j,mask_tache_ij) 
     638    !~ !$OMP DO 
     639    !~ do j=3,ny-2 
     640    !~    do i=3,nx-2 
     641    !~ test1:  if (.not.iceberg(table_out(i,j))) then ! on est pas sur un iceberg  
     642    !~          if (nb_pts_tache(table_out(i,j)).ge.1) then 
     643    !~             ice(i,j)=1 
     644    !~             if (nb_pts_tache(table_out(i,j)).le.10) then  ! les petites iles sont en sia 
     645    !~ !    write(6,*) 'petite ile ',i,j 
     646    !~               flgzmx(i,j)=.false. 
     647    !~                flgzmx(i+1,j)=.false. 
     648    !~                flgzmy(i,j)=.false. 
     649    !~                flgzmy(i,j+1)=.false. 
     650    !~                gzmx(i:i+1,j)=.false. 
     651    !~                gzmy(i,j:j+1)=.false. 
     652    !~             endif 
     653    !~  
     654    !~  
     655    !~ ! ici on est sur une tache non iceberg >= 5 points                        
     656    !~ ! on teste si la tache n'est pas completement ice stream 
     657    !~  
     658    !~ test2:       if (icetrim(table_out(i,j))) then   ! on a une ile d'ice stream 
     659    !~  
     660    !~    mask_tache_ij(:,:)=.false. 
     661    !~    mask_tache_ij(:,:)=(table_out(:,:).eq.table_out(i,j))       ! pour toute la tache 
     662    !~  
     663    !~      smax_=maxval(S(:,:),MASK=mask_tache_ij(:,:)) 
     664    !~      smax_coord(:)=maxloc(S(:,:),MASK=mask_tache_ij(:,:)) 
     665    !~      smax_i=smax_coord(1) 
     666    !~      smax_j=smax_coord(2) 
     667    !~  
     668    !~ !!$               smax_i=0 ; smax_j=0 ; smax_=sealevel 
     669    !~ !!$               do ii=3,nx-2 
     670    !~ !!$                  do jj=3,ny-2 
     671    !~ !!$                     if (table_out(ii,jj).eq.table_out(i,j)) then 
     672    !~ !!$                        if (s(ii,jj).gt.smax_) then                            
     673    !~ !!$                           smax_ =s(ii,jj) 
     674    !~ !!$                           smax_i=ii 
     675    !~ !!$                           smax_j=jj 
     676    !~ !!$                        endif 
     677    !~ !!$                     endif 
     678    !~ !!$                  end do 
     679    !~ !!$               end do 
     680    !~  
     681    !~  
     682    !~                gzmx(smax_i,smax_j)=.false. ; gzmx(smax_i+1,smax_j)=.false. 
     683    !~                gzmy(smax_i,smax_j)=.false. ; gzmx(smax_i,smax_j+1)=.false. 
     684    !~                flgzmx(smax_i,smax_j)=.false. ; flgzmx(smax_i+1,smax_j)=.false. 
     685    !~                flgzmy(smax_i,smax_j)=.false. ; flgzmx(smax_i,smax_j+1)=.false. 
     686    !~   
     687    !~                if (Smax_.le.sealevel) then 
     688    !~                   write(num_tracebug,*)'Attention, une ile avec la surface sous l eau' 
     689    !~                   write(num_tracebug,*)'time=',time,'   coord:',smax_i,smax_j 
     690    !~                end if 
     691    !~  
     692    !~             endif test2  
     693    !~          end if                                             ! endif deplace 
     694    !~           
     695    !~       else ! on est sur un iceberg                          !   test1 
     696    !~          ice(i,j)=0 
     697    !~          h(i,j)=1. 
     698    !~          surnet=H(i,j)*(1.-ro/row) 
     699    !~          S(i,j)=surnet+sealevel 
     700    !~          B(i,j)=S(i,j)-H(i,j) 
     701    !~           
     702    !~       endif test1 
     703    !~  
     704    !~  
     705    !~    end do 
     706    !~ end do 
     707    !~ !$OMP END DO 
     708    !~ !$OMP END PARALLEL 
     709    !~  
     710    !~ !---------------------------------------------- 
     711    !~ ! On caracterise le front des ice shelfs/streams 
     712    !~  
     713    !~ !      call DETERMIN_FRONT     
     714    !~  
     715    !~ !----------------------------------------------       
     716    !~ !!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf) 
     717    !~ !!$if (itestf.gt.0) then 
     718    !~ !!$   write(6,*) 'dans flottab apres DETERMIN_front  asymetrie sur H  pour time=',time 
     719    !~ !!$   stop 
     720    !~ !!$else 
     721    !~ !!$   write(6,*) 'dans flottab apres DETERMIN_front pas d asymetrie sur H  pour time=',time 
     722    !~ !!$ 
     723    !~ !!$end if 
     724    !~  
     725    !~ endif synchro 
     726 
     727    ! correction momentanée pour symetrie Heino 
     728    !where ((.not.flot(:,:)).and.(ice(:,:).eq.0)) H(:,:)=0.       
     729 
     730    !fin de routine flottab2 
     731    !print*, 'front',front(50,30),ice(50,30),flotmx(i,j),uxbar(i,j) 
     732    !~ print*,'fin flottab',S(132,183),H(132,183),BSOC(132,183),B(132,183),sealevel 
     733    !~ print*,'fin flottab',flot(132,183),ice(132,183),gzmx(132,183),gzmy(132,183),ilemx(132,183),ilemy(132,183) 
     734    !read(5,*) 
     735  end subroutine flottab 
    721736!--------------------------------------------------------------------   
    722737 
Note: See TracChangeset for help on using the changeset viewer.