Changeset 161


Ignore:
Timestamp:
11/21/17 17:00:25 (6 years ago)
Author:
dumas
Message:

Attempt to debug the classical crash in diagnoshelf : determin_tache now called by steps_time_loop | isostasy after diagnoshelf | flgzmx & flgzmy .false. if there is no ice

Location:
trunk/SOURCES
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/SOURCES/dragging_param_beta_mod.f90

    r127 r161  
    230230flgzmy(:,:) = flgzmy(:,:) .or. gzmy(:,:) 
    231231 
     232where (hmx(:,:).eq.0.) 
     233  flgzmx(:,:) = .false. 
     234endwhere 
     235where (hmy(:,:).eq.0.) 
     236  flgzmy(:,:) = .false. 
     237endwhere 
     238 
    232239fleuvemx(:,:)=gzmx(:,:) 
    233240fleuvemy(:,:)=gzmy(:,:) 

  • trunk/SOURCES/dragging_param_beta_nolin_mod.f90

    r152 r161  
    188188where (ilemy(:,:)) betamy(:,:) = betamy(:,:) * coef_ile 
    189189 
    190 uslmag(:,:) = (ux(:,:,nz)**2+uy(:,:,nz)**2)**(0.5) 
    191  
    192 where (uslmag(:,:).gt.0) 
    193    betamx(:,:) = betamx(:,:) / uslmag(:,:) 
    194    betamy(:,:) = betamy(:,:) / uslmag(:,:) 
     190!uslmag(:,:) = (ux(:,:,nz)**2+uy(:,:,nz)**2)**(0.5) 
     191!where (uslmag(:,:).gt.0) 
     192!   betamx(:,:) = betamx(:,:) / uslmag(:,:) 
     193!   betamy(:,:) = betamy(:,:) / uslmag(:,:) 
     194!endwhere 
     195 
     196where (abs(ux(:,:,nz)).gt.0) 
     197   betamx(:,:) = betamx(:,:) / min(abs(ux(:,:,nz)),100.) 
     198endwhere 
     199where (abs(uy(:,:,nz)).gt.0) 
     200   betamy(:,:) = betamy(:,:) / min(abs(uy(:,:,nz)),100.) 
    195201endwhere 
    196202 
    … …  
    238244flgzmy(:,:) = flgzmy(:,:) .or. gzmy(:,:) 
    239245 
     246where (hmx(:,:).eq.0.) 
     247  flgzmx(:,:) = .false. 
     248endwhere 
     249where (hmy(:,:).eq.0.) 
     250  flgzmy(:,:) = .false. 
     251endwhere 
     252 
    240253fleuvemx(:,:)=gzmx(:,:) 
    241254fleuvemy(:,:)=gzmy(:,:) 

  • trunk/SOURCES/flottab2-0.7.f90

    r156 r161  
    554554    flgzmx(:,:)=(marine.and.(flotmx(:,:).or.gzmx(:,:).or.ilemx(:,:)))   & 
    555555         .or.(.not.marine.and.flotmx(:,:)) 
     556    where (hmx(:,:).eq.0.) 
     557      flgzmx(:,:) = .false. 
     558    endwhere 
    556559    flgzmy(:,:)=(marine.and.(flotmy(:,:).or.gzmy(:,:).or.ilemy(:,:)))   & 
    557560         .or.(.not.marine.and.flotmy(:,:)) 
     561    where (hmy(:,:).eq.0.) 
     562      flgzmy(:,:) = .false. 
     563    endwhere 
    558564    !$OMP END WORKSHARE 
    559565 
    … …  
    610616    !$OMP END PARALLEL 
    611617!print*,'H(90,179) flottab 3',H(90,179),ice(90,179),flot(90,179),Bsoc(90,179)+H(90,179)*ro/row -sealevel 
    612     call DETERMIN_TACHE  
    613     !~  
    614     !~ synchro :    if (isynchro.eq.1) then 
    615     !~ !!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf) 
    616     !~ !!$if (itestf.gt.0) then 
    617     !~ !!$   write(6,*) 'dans flottab avant DETERMIN_TACHE  asymetrie sur H  pour time=',time 
    618     !~ !!$   stop 
    619     !~ !!$else 
    620     !~ !!$   write(6,*) 'dans flottab apres  DETERMIN_TACHE pas d asymetrie sur H  pour time=',time 
    621     !~ !!$ 
    622     !~ !!$end if 
    623     !~  
    624     !~  
    625     !~ !----------------------------------------------! 
    626     !~ !On determine les differents ice strean/shelf  ! 
    627     !~ !      call DETERMIN_TACHE                      ! 
    628     !~ !----------------------------------------------! 
    629     !~  
    630     !~  
    631     !~ !!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf) 
    632     !~ !!$if (itestf.gt.0) then 
    633     !~ !!$   write(6,*) 'dans flottab apres DETERMIN_TACHE  asymetrie sur H  pour time=',time 
    634     !~ !!$   stop 
    635     !~ !!$else 
    636     !~ !!$   write(6,*) 'dans flottab apres DETERMIN_TACHE pas d asymetrie sur H  pour time=',time 
    637     !~ !!$ 
    638     !~ !!$end if 
    639      
    640     !~ !ice(:,:)=0  ! Attention, voir si ca marche toujours pour l'Antarctique et heminord ! 
    641      
    642     !On compte comme englacé uniquement les calottes dont une partie est posée       
    643     !$OMP PARALLEL PRIVATE(smax_,smax_coord,smax_i,smax_j,mask_tache_ij) 
    644     !$OMP DO 
    645     do j=3,ny-2 
    646        do i=3,nx-2 
    647     test1:  if (.not.iceberg1D(table_out(i,j))) then ! on est pas sur un iceberg  
    648              if (nb_pts_tache(table_out(i,j)).ge.1) then 
    649                  ice(i,j)=1 
    650     !~             if (nb_pts_tache(table_out(i,j)).le.10) then  ! les petites iles sont en sia 
    651     !~ !    write(6,*) 'petite ile ',i,j 
    652     !~               flgzmx(i,j)=.false. 
    653     !~                flgzmx(i+1,j)=.false. 
    654     !~                flgzmy(i,j)=.false. 
    655     !~                flgzmy(i,j+1)=.false. 
    656     !~                gzmx(i:i+1,j)=.false. 
    657     !~                gzmy(i,j:j+1)=.false. 
    658     !~             endif 
    659     !~  
    660     !~  
    661     ! ici on est sur une tache non iceberg >= 5 points                        
    662     ! on teste si la tache n'est pas completement ice stream 
    663      
    664     test2:       if (icetrim(table_out(i,j))) then   ! on a une ile d'ice stream 
    665      
    666        mask_tache_ij(:,:)=.false. 
    667        mask_tache_ij(:,:)=(table_out(:,:).eq.table_out(i,j))       ! pour toute la tache 
    668      
    669          smax_=maxval(S(:,:),MASK=mask_tache_ij(:,:)) 
    670          smax_coord(:)=maxloc(S(:,:),MASK=mask_tache_ij(:,:)) 
    671          smax_i=smax_coord(1) 
    672          smax_j=smax_coord(2) 
    673      
    674     !~ !!$               smax_i=0 ; smax_j=0 ; smax_=sealevel 
    675     !~ !!$               do ii=3,nx-2 
    676     !~ !!$                  do jj=3,ny-2 
    677     !~ !!$                     if (table_out(ii,jj).eq.table_out(i,j)) then 
    678     !~ !!$                        if (s(ii,jj).gt.smax_) then                            
    679     !~ !!$                           smax_ =s(ii,jj) 
    680     !~ !!$                           smax_i=ii 
    681     !~ !!$                           smax_j=jj 
    682     !~ !!$                        endif 
    683     !~ !!$                     endif 
    684     !~ !!$                  end do 
    685     !~ !!$               end do 
    686      
    687      
    688                    gzmx(smax_i,smax_j)=.false. ; gzmx(smax_i+1,smax_j)=.false. 
    689                    gzmy(smax_i,smax_j)=.false. ; gzmx(smax_i,smax_j+1)=.false. 
    690                    flgzmx(smax_i,smax_j)=.false. ; flgzmx(smax_i+1,smax_j)=.false. 
    691                    flgzmy(smax_i,smax_j)=.false. ; flgzmx(smax_i,smax_j+1)=.false. 
    692       
    693                    if (Smax_.le.sealevel) then 
    694                       write(num_tracebug,*)'Attention, une ile avec la surface sous l eau' 
    695                       write(num_tracebug,*)'time=',time,'   coord:',smax_i,smax_j 
    696                    end if 
    697      
    698                 endif test2  
    699              end if                                             ! endif deplace 
    700 !cdc transfere dans calving :              
    701            else ! on est sur un iceberg                          !   test1 
    702                 iceberg(i,j)=iceberg1D(table_out(i,j)) 
    703 !~              ice(i,j)=0 
    704 !~              h(i,j)=0. !1. afq, we should put everything in calving! 
    705 !~              surnet=H(i,j)*(1.-ro/row) 
    706 !~              S(i,j)=surnet+sealevel 
    707 !~              B(i,j)=S(i,j)-H(i,j) 
    708               
    709           endif test1 
    710      
    711      
    712        end do 
    713     end do 
    714      
    715     !$OMP END DO 
    716     !$OMP END PARALLEL 
    717  
    718     !~  
    719     !~ !---------------------------------------------- 
    720     !~ ! On caracterise le front des ice shelfs/streams 
    721     !~  
    722     !~ !      call DETERMIN_FRONT     
    723     !~  
    724     !~ !----------------------------------------------       
    725     !~ !!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf) 
    726     !~ !!$if (itestf.gt.0) then 
    727     !~ !!$   write(6,*) 'dans flottab apres DETERMIN_front  asymetrie sur H  pour time=',time 
    728     !~ !!$   stop 
    729     !~ !!$else 
    730     !~ !!$   write(6,*) 'dans flottab apres DETERMIN_front pas d asymetrie sur H  pour time=',time 
    731     !~ !!$ 
    732     !~ !!$end if 
    733     !~  
    734     !~ endif synchro 
    735  
    736     ! correction momentanée pour symetrie Heino 
    737     !where ((.not.flot(:,:)).and.(ice(:,:).eq.0)) H(:,:)=0.       
    738  
    739     !fin de routine flottab2 
    740     !print*, 'front',front(50,30),ice(50,30),flotmx(i,j),uxbar(i,j) 
    741     !~ print*,'fin flottab',S(132,183),H(132,183),BSOC(132,183),B(132,183),sealevel 
    742     !~ print*,'fin flottab',flot(132,183),ice(132,183),gzmx(132,183),gzmy(132,183),ilemx(132,183),ilemy(132,183) 
    743     !read(5,*) 
    744  
     618  
    745619    !call determin_front ! cette version ne conserve pas la masse !!! 
    746620    call determin_front_tof ! version simplifiee 
    … …  
    796670  do i=2,nx-1 
    797671 
    798  
    799  
    800672     IF (ice(i,j).ge.1) THEN ! on est sur la glace-----------------------------! 
    801673                         
    … …  
    826698 
    827699           ! si 2 taches differentes sont dans le masque, il faut les identifier dans compt 
    828            ! i.e. les plus grands numeros correspondent au plus petit  
    829            ! on lui affecte le numero de la tache fondamentale avec un signe - 
    830            ! pour indiquer le changement 
     700           ! on lui affecte le numero de la tache fondamentale 
    831701 
    832702           do indice=1,mask_nb 
    … …  
    835705              endif 
    836706           enddo 
     707                                                                   ! exemple on est sur le point X : 5  X 
     708          do indice=1,mask_nb                                      !                                   20 
     709            if(mask(indice).gt.label) then                         ! mask(2)=20 > 5              
     710              compt(mask(indice))=label                            ! compt(20)=5 
     711              if (.not.iceberg1D(mask(indice))) iceberg1D(label)=.false. ! si la tache n'etais pas un iceberg => iceberg =.false.  iceberg(5)=.false.                
     712              if (.not.icetrim(mask(indice))) icetrim(label)=.false. 
     713              where (table_out(:,:).eq.mask(indice))               ! where table_out(:,:)=mask(2)=20 
     714                table_out(:,:)=label                        ! table_out(:,:)=label=5                  
     715              endwhere 
     716            endif 
     717          enddo    
    837718 
    838719        else !aucun des voisins est une tache 
    … …  
    851732              
    852733           label_max  = label_max+1 
    853            if (label_max.gt.n_ta_max) print*,'trop de taches=',label_max  
     734           if (label_max.gt.n_ta_max) print*,'ATTENTION trop de taches icebergs=',label_max  
    854735        endif 
    855736 
    … …  
    869750! On indique aussi le nb de point que contient chaque taches (nb_pts_tache) 
    870751 
    871 do indice=1,label_max 
    872    vartemp = compt(indice) 
    873    if (compt(indice).lt.0) then  
    874       compt(indice)= compt(-vartemp) 
    875       if (.not.iceberg1D(indice)) iceberg1D(-vartemp)=.false. 
    876       if (.not.icetrim(indice)) icetrim(-vartemp)=.false. 
    877    endif     
    878 enddo 
    879  
    880 !$OMP PARALLEL 
    881 !$OMP DO REDUCTION(+:nb_pts_tache) 
    882 do j=1,ny 
    883    do i=1,nx 
     752  !$OMP PARALLEL 
     753  !$OMP DO 
     754  do j=1,ny 
     755    do i=1,nx 
    884756      if (table_out(i,j).ne.0) then 
    885          table_out(i,j)=compt(table_out(i,j))   
    886          nb_pts_tache(compt(table_out(i,j)))= nb_pts_tache(compt(table_out(i,j)))+1    
     757        nb_pts_tache(compt(table_out(i,j)))= nb_pts_tache(compt(table_out(i,j)))+1 
    887758      endif 
    888    enddo 
    889 enddo 
    890 !$OMP END DO 
    891 !$OMP END PARALLEL 
    892  
    893  
    894 !!$tablebis(:,:)=table_out(:,:) 
    895 !!$do j=1,ny 
    896 !!$   do i=1,nx 
    897 !!$      if (tablebis(i,j).ne.0) then     ! tache de glace 
    898 !!$         numtache=table_out(i,j) 
    899 !!$         nb_pt=count(table_out(:,:).eq.numtache)  ! compte tous les points de la tache 
    900 !!$         nb_pts_tache(table_out(i,j))=nb_pt       !  
    901 !!$ 
    902 !!$         where (table_out(:,:).eq.numtache)           
    903 !!$            tablebis(:,:)=0                       ! la table de tache est remise a 0 pour eviter de repasser 
    904 !!$         end where 
    905 !!$          write(6,*) i,j,nb_pt,table_out(i,j) 
    906 !!$      endif 
    907 !!$   enddo 
    908 !!$enddo 
    909  
    910 !!$do j=1,ny 
    911 !!$   do i=1,nx 
    912 !!$      debug_3d(i,j,56)=nb_pts_tache(table_out(i,j)) 
    913 !!$   end do 
    914 !!$end do 
    915 !!$ 
     759    enddo 
     760  enddo 
     761  !$OMP END DO 
     762  !$OMP END PARALLEL 
     763 
     764  !On compte comme englacé uniquement les calottes dont une partie est posée       
     765  !$OMP PARALLEL PRIVATE(smax_,smax_coord,smax_i,smax_j,mask_tache_ij) 
     766  !$OMP DO 
     767  do j=3,ny-2 
     768    do i=3,nx-2 
     769      test1: if (.not.iceberg1D(table_out(i,j))) then ! on est pas sur un iceberg  
     770        if (nb_pts_tache(table_out(i,j)).ge.1) then 
     771          ice(i,j)=1 
     772          ! ici on est sur une tache non iceberg >= 5 points                        
     773          ! on teste si la tache n'est pas completement ice stream 
     774     
     775          test2: if (icetrim(table_out(i,j))) then   ! on a une ile d'ice stream 
     776     
     777            mask_tache_ij(:,:)=.false. 
     778            mask_tache_ij(:,:)=(table_out(:,:).eq.table_out(i,j))       ! pour toute la tache 
     779     
     780            smax_=maxval(S(:,:),MASK=mask_tache_ij(:,:)) 
     781            smax_coord(:)=maxloc(S(:,:),MASK=mask_tache_ij(:,:)) 
     782            smax_i=smax_coord(1) 
     783            smax_j=smax_coord(2) 
     784   
     785            gzmx(smax_i,smax_j)=.false. ; gzmx(smax_i+1,smax_j)=.false. 
     786            gzmy(smax_i,smax_j)=.false. ; gzmx(smax_i,smax_j+1)=.false. 
     787            flgzmx(smax_i,smax_j)=.false. ; flgzmx(smax_i+1,smax_j)=.false. 
     788            flgzmy(smax_i,smax_j)=.false. ; flgzmx(smax_i,smax_j+1)=.false. 
     789      
     790            if (Smax_.le.sealevel) then 
     791              write(num_tracebug,*)'Attention, une ile avec la surface sous l eau' 
     792              write(num_tracebug,*)'time=',time,'   coord:',smax_i,smax_j 
     793            end if     
     794          endif test2  
     795        end if                                             ! endif deplace 
     796!cdc transfere dans calving :              
     797      else ! on est sur un iceberg                          !   test1 
     798        iceberg(i,j)=iceberg1D(table_out(i,j)) 
     799!~        ice(i,j)=0 
     800!~        h(i,j)=0. !1. afq, we should put everything in calving! 
     801!~        surnet=H(i,j)*(1.-ro/row) 
     802!~        S(i,j)=surnet+sealevel 
     803!~        B(i,j)=S(i,j)-H(i,j)              
     804      endif test1     
     805    end do 
     806  end do     
     807  !$OMP END DO 
     808  !$OMP END PARALLEL 
    916809 
    917810debug_3D(:,:,121)=real(table_out(:,:)) 

  • trunk/SOURCES/steps_time_loop.f90

    r157 r161  
    5858     call icethick3 
    5959     call flottab 
     60     call determin_tache 
    6061     call calving 
    6162     call ablation_bord 
    … …  
    304305     end if                                      ! Dans le shelf flgz = Ubar = Vbil et udef=0 
    305306 
    306  
    307      ! isostasy 
    308      !================= 
    309  
    310  
    311      spinup_3_bed: if (ispinup.eq.0.or.nt.lt.nt_init_tm) then 
    312         !--------------------------------------------------------------------------------- 
    313         ! spinup_3_bed 
    314         !---------------- 
    315         ! Ne passe dans ce bloc que quand on veut l'isostasie, donc si variations epaisseur 
    316         ! run standard  (ispinup=0) et pour les initialisations                   
    317         !--------------------------------------------------------------------------------- 
    318  
    319         if (itracebug.eq.1)  call tracebug(' Dans spinup_3_bed') 
    320  
    321         if ((nbed.eq.1).and.nt.ne.1.and.isynchro.eq.1.and.(mod(abs(TIME),50.).lt.dtmin)) then 
    322            call bedrock                                              !  bedrock adjustment 
    323         endif 
    324  
    325         call lake_to_slv                                             ! pour traiter dynamiquement les shelves sur lac 
    326  
    327      end if spinup_3_bed 
    328  
    329307     call tracer      ! aurel neem (attention, position a verifier)    
    330308 
    … …  
    430408 
    431409  end if spinup_4_vitdyn 
     410  
     411   
     412  ! isostasy 
     413  !================= 
     414 
     415  spinup_3_bed: if ((ISYNCHRO.eq.1).and.(ispinup.eq.0.or.nt.lt.nt_init_tm)) then 
     416  !--------------------------------------------------------------------------------- 
     417  ! spinup_3_bed 
     418  !---------------- 
     419  ! Ne passe dans ce bloc que quand on veut l'isostasie, donc si variations epaisseur 
     420  ! run standard  (ispinup=0) et pour les initialisations                   
     421  !--------------------------------------------------------------------------------- 
     422 
     423  if (itracebug.eq.1)  call tracebug(' Dans spinup_3_bed') 
     424 
     425  if ((nbed.eq.1).and.nt.ne.1.and.isynchro.eq.1.and.(mod(abs(TIME),50.).lt.dtmin)) then 
     426    call bedrock                 !  bedrock adjustment 
     427  endif 
     428 
     429  call lake_to_slv               ! pour traiter dynamiquement les shelves sur lac 
     430 
     431  end if spinup_3_bed 
     432   
    432433 
    433434end subroutine step_thermomeca 

  • trunk/SOURCES/steps_time_loop_avec_iterbeta.f90

    r157 r161  
    5858     call icethick3 
    5959     call flottab 
     60     call determin_tache 
    6061     call calving 
    6162     call ablation_bord 
    … …  
    305306     end if                                      ! Dans le shelf flgz = Ubar = Vbil et udef=0 
    306307 
    307  
    308      ! isostasy 
    309      !================= 
    310  
    311  
    312      spinup_3_bed: if (ispinup.eq.0.or.nt.lt.nt_init_tm) then 
    313         !--------------------------------------------------------------------------------- 
    314         ! spinup_3_bed 
    315         !---------------- 
    316         ! Ne passe dans ce bloc que quand on veut l'isostasie, donc si variations epaisseur 
    317         ! run standard  (ispinup=0) et pour les initialisations                   
    318         !--------------------------------------------------------------------------------- 
    319  
    320         if (itracebug.eq.1)  call tracebug(' Dans spinup_3_bed') 
    321  
    322         if ((nbed.eq.1).and.nt.ne.1.and.isynchro.eq.1.and.(mod(abs(TIME),50.).lt.dtmin)) then            
    323            call bedrock()                                            !  bedrock adjustment 
    324         endif 
    325  
    326         call lake_to_slv                                             ! pour traiter dynamiquement les shelves sur lac 
    327  
    328      end if spinup_3_bed 
    329  
    330308     call tracer      ! aurel neem (attention, position a verifier)    
    331309 
    … …  
    441419 
    442420  end if spinup_4_vitdyn 
     421   
     422  ! isostasy 
     423  !================= 
     424 
     425  spinup_3_bed: if ((ISYNCHRO.eq.1).and.(ispinup.eq.0.or.nt.lt.nt_init_tm)) then 
     426  !--------------------------------------------------------------------------------- 
     427  ! spinup_3_bed 
     428  !---------------- 
     429  ! Ne passe dans ce bloc que quand on veut l'isostasie, donc si variations epaisseur 
     430  ! run standard  (ispinup=0) et pour les initialisations                   
     431  !--------------------------------------------------------------------------------- 
     432 
     433  if (itracebug.eq.1)  call tracebug(' Dans spinup_3_bed') 
     434 
     435  if ((nbed.eq.1).and.nt.ne.1.and.isynchro.eq.1.and.(mod(abs(TIME),50.).lt.dtmin)) then 
     436    call bedrock                 !  bedrock adjustment 
     437  endif 
     438 
     439  call lake_to_slv               ! pour traiter dynamiquement les shelves sur lac 
     440 
     441  end if spinup_3_bed 
     442     
    443443 
    444444end subroutine step_thermomeca 
Note: See TracChangeset for help on using the changeset viewer.