Changeset 122 for branches


Ignore:
Timestamp:
06/29/17 16:27:48 (7 years ago)
Author:
aquiquet
Message:

Updating iLOVECLIM branch (calving / flottab / outputs) , working version but will be merged to trunk soon for water conservation

Location:
branches/iLoveclim
Files:
1 added
5 edited

Legend:

Unmodified
Added
Removed

  • branches/iLoveclim/SOURCES/Hemin40_files/lect-hemin40_mod.f90

    r88 r122  
    113113 
    114114! calcul des courbures du socle 
    115 !     call courbure(nx,ny,dx,Bsoc,bidon(:,:,1),bidon(:,:,2),bidon(:,:,3), & 
    116 !          bidon(:,:,4),socle_cry,bidon(:,:,5)) 
    117 !     socle_cry(:,:)=socle_cry(:,:)*dx*dx 
     115     call courbure(nx,ny,dx,Bsoc,bidon(:,:,1),bidon(:,:,2),bidon(:,:,3), & 
     116          bidon(:,:,4),socle_cry,bidon(:,:,5)) 
     117     socle_cry(:,:)=socle_cry(:,:)*dx*dx 
    118118 
    119119! lecture des coordonnées geographiques 

  • branches/iLoveclim/SOURCES/Hemin40_files/output_hemin40_mod.f90

    r90 r122  
    22 
    33       USE module3D_phy 
    4  
     4       use calving_frange, only : hcoup  !afq -- needed for the freshwater flux calculation 
     5        
    56implicit none 
    67 
    … …  
    2728!real ::  pdfmean                      ! moyenne PDF 
    2829logical,dimension(nx,ny,13) :: mask_cal ! masque regions calotte  
    29 REAL, dimension(nx,ny) :: old_H_dtt, old_H_new    ! Epaisseur de glace au pas de temps precedent 
    30 real, dimension(nx,ny) :: runof_oc_sngl 
     30REAL, dimension(nx,ny) :: old_H_dtt      ! Epaisseur de glace au pas de temps precedent 
    3131integer, dimension(nx,ny) :: write_mask 
     32 
    3233CONTAINS 
    3334!_________________________________________________________________________ 
    … …  
    181182!      REAL ABLA(NX,NY) 
    182183      REAL DELTAVOL 
    183       REAL, dimension(nx,ny,13) :: delta_H_dtt 
    184       REAL, dimension(nx,ny) :: delta_H_test 
     184      REAL, dimension(nx,ny) :: delta_H_dtt 
    185185 
    186186 
    … …  
    198198      BDOTMEAN=0. 
    199199 
     200      runof_oc(:,:)    = 0. 
     201      bmelt_oc(:,:)    = 0. 
     202      delta_H_dtt(:,:) = 0. 
     203       
    200204      DO kk = 1,13 
    201205        INP(kk) = 0 
    … …  
    212216        TBM(kk) = 0. 
    213217        ITJJA(kk) = 0. 
    214         delta_H_dtt(:,:,kk) = 0. 
    215         runof_oc(:,:) = 0. 
    216218         
    217219! nouveau tof mai 2009 
    … …  
    237239     enddo 
    238240 
    239  
    240      where (mask_cal(:,:,1).and.(H(:,:).gt.2.).and..not.flot(:,:))  
    241           delta_H_dtt(:,:,1) = H(:,:) 
    242      elsewhere 
    243           delta_H_dtt(:,:,1) = 0d0 
     241! afq -- thickness variation  
     242     delta_H_dtt(:,:) = old_H_dtt(:,:) - H(:,:) 
     243! afq -- continental runoff      
     244!        this goes to ECBilt, will be sum up, then routed to CLIO 
     245!     where (mask_cal(:,:,1).and.((H(:,:).gt.1.).or.(old_H_dtt(:,:).gt.1.)).and..not.flot(:,:)) 
     246!     afq I think we do not need the H condition: 
     247     where  (mask_cal(:,:,1).and..not.flot(:,:)) 
     248        runof_oc(:,:) = runof_oc(:,:) + DBLE(delta_H_dtt(:,:))*DX*DY 
    244249     endwhere 
    245      old_H_new(:,:) = delta_H_dtt(:,:,1) 
    246      delta_H_dtt(:,:,1) = old_H_dtt(:,:) - delta_H_dtt(:,:,1) 
     250! afq -- shelf runoff 
     251!        this goes to bilinear interpolation, then put to CLIO 
     252!        we have to make sure that we are not considering a "calved" point 
     253     where (mask_cal(:,:,1).and.(H(:,:).gt.1.).and.flot(:,:).and.(calv(:,:).eq.0.)) 
     254        bmelt_oc(:,:) = bmelt_oc(:,:) + DBLE(delta_H_dtt(:,:))*DX*DY 
     255     endwhere 
     256 
     257! afq -- old_h_dtt has to be reset 
     258     old_H_dtt(:,:)  = H(:,:) 
     259 
     260! afq -- the accumulated calving need to be reset, BUT not here because of coupling freq.     
     261     calvin_GRIS(:,:) = calvin_GRIS(:,:)+calvingGRIS(:,:) 
     262     calvingGRIS(:,:) = 0. 
    247263      
    248      runof_oc(:,:) = DBLE(delta_H_dtt(:,:,1))*DX*DY 
    249       
    250      runof_oc_sngl(:,:) = delta_H_dtt(:,:,1)*DX*DY 
    251  
    252       
    253      old_H_dtt(:,:)  = old_H_new(:,:) 
    254       
    255      calvin_GRIS(:,:) = calvingGRIS(:,:) 
    256      calvingGRIS(:,:) =0. 
    257       
    258 !!$ 
    259 !!$ 
    260 !!$            IF(H(i,j).gt.2.)  THEN  
    261 !!$             if (mk(i,j).eq.1) then 
    262 !!$               INF(2) = INF(2) + 1  
    263 !!$               ISVOLF(2) = ISVOLF(2) + H(I,J) 
    264 !!$             else 
    265 !!$               INP(2) = INP(2) + 1 
    266 !!$               ISVOL(2) = ISVOL(2) + H(I,J) 
    267 !!$               ISACC(2) = ISACC(2) + ACC(I,J) 
    268 !!$               ISABL(2) = ISABL(2) + ABL(I,J) 
    269 !!$               ITJJA(2) = ITJJA(2) + TJULY(I,J) 
    270 !!$               if (H(I,J).gt.HMAX_(2)) HMAX_(2)=H(I,J)  
    271 !!$             endif 
    272 !!$            ENDIF 
    273 !!$            TACC(2) = TACC(2) + ACC(I,J) 
    274 !!$             ISABLBORD(2) = ISABLBORD(2) + ABLBORD(I,J) 
    275 !!$             ABLATOT(2) = ISABL(2) + ISABLBORD(2) 
    276 !!$             ISCALV(2) = ISCALV(2) + CALV(I,J) 
    277 !!$             ISBM(2) = ISACC(2)+ISABL(2)+ISCALV(2)+ISABLBORD(2) 
    278 !!$             TBM(2) = TACC(2)+ISABL(2)+ISCALV(2)+ISABLBORD(2) 
    279 !!$        ENDIF 
    280 !!$ 
    281 !!$ 
    282 !!$      END DO   
    283  
    284264 
    285265   

  • branches/iLoveclim/SOURCES/calving_frange.f90

    r90 r122  
    2121real, dimension (nx,ny) :: hmhc  ! hauteur au dessus de la coupure 
    2222real, dimension (nx,ny) :: bil_tot ! bilan surface et fond (bm-bmelt) 
    23 real :: hcoup ! epaiseur de coupure au temps time 
    24 real :: hcoup_max  ! epaisseur max 
    25 real :: hcoup_min=100.  ! epaisseur min 
     23real, dimension (nx,ny) :: hcoup ! epaiseur de coupure au temps time 
    2624real :: mass_total 
     25real :: hcoup_plateau  ! coupure points peu profonds 
     26real :: hcoup_abysses  ! coupure points ocean profond 
     27real :: prof_plateau  ! profondeur max des points peu profonds 
     28real :: prof_abysses  ! profondeur min des points ocean profond 
    2729integer :: meth_hcoup  ! pour avoir hcoup dépendant du coefbmshelf 
    2830integer :: ifrange        !  0 pas de traitement particulier pres du bord, 1 -> franges 
    … …  
    3941subroutine init_calving 
    4042 
    41 namelist/calving/Hcoup,ifrange,meth_hcoup 
     43namelist/calving/Hcoup_plateau,Hcoup_abysses,prof_plateau,prof_abysses,ifrange,meth_hcoup 
    4244 
    4345 
    … …  
    5961write(num_rep_42,428) '&calving              ! nom du bloc calving méthode Vincent ' 
    6062write(num_rep_42,*) 
    61 write(num_rep_42,*)   'Hcoup      = ', Hcoup 
     63write(num_rep_42,*)   'Hcoup_plateau  = ', Hcoup_plateau 
     64write(num_rep_42,*)   'Hcoup_abysses  = ', Hcoup_abysses 
     65write(num_rep_42,*)   'prof_plateau  = ', prof_plateau 
     66write(num_rep_42,*)   'prof_abysses  = ', prof_abysses 
    6267write(num_rep_42,*)   'ifrange    = ',ifrange 
    6368write(num_rep_42,*)   'meth_hcoup = ', meth_hcoup 
    6469write(num_rep_42,*)'/' 
    65                              
    66 write(num_rep_42,428) '! Hcoup epaisseur de coupure, le max si attache a coefbmshelf ' 
     70 
     71write(num_rep_42,428) '! Hcoup epaisseurs de coupure pour les zones peu prodondes et profondes' 
     72write(num_rep_42,428) '! Hcoup_plateau<Hcoup_abysses && prof_plateau<prof_abysses' 
     73write(num_rep_42,428) '! prof profondeur delimitant les zones peu prodondes et profondes' 
    6774write(num_rep_42,428) '! ifrange=0 -> pas de traitement particulier sur les bords' 
    6875write(num_rep_42,428) '! ifrange=1 -> traitement de Vincent avec ice shelves frangeants partout' 
    6976write(num_rep_42,428) '! ifrange=2 -> ice shelves frangeants seulement si bm-bmelt positif' 
    7077write(num_rep_42,*)   '! meth_hcoup pour faire eventuellement varier Hcoup avec le climat' 
    71 write(num_rep_42,*)   '! Hcoup_min=',Hcoup_min 
     78write(num_rep_42,*)   '! Hcoup_plateau=',Hcoup_plateau 
     79write(num_rep_42,*)   '! Hcoup_abysses=',Hcoup_abysses 
     80write(num_rep_42,*)   '! prof_plateau=',prof_plateau 
     81write(num_rep_42,*)   '! prof_abysses=',prof_abysses 
    7282write(num_rep_42,*) 
    7383 
    74 Hcoup_max=Hcoup 
     84 
     85! afq -- coupure depend de la profondeur: 
     86! 
     87!  hcoup                       prof_abysses 
     88!   ^                               v  
     89!   |                                _______ hcoup_abysses 
     90!   |                               / 
     91!   |                              / 
     92!   |                             / 
     93!   |                            / 
     94!   |      hcoup_plateau _______/ 
     95!                               ^ 
     96!                          prof_plateau      
     97 
     98Hcoup(:,:) = min ( max(                                         & 
     99     (-(Bsoc0(:,:)-sealevel) - prof_plateau)/(prof_abysses-prof_plateau)    & 
     100     *(hcoup_abysses-hcoup_plateau)+hcoup_plateau               & 
     101     , hcoup_plateau), hcoup_abysses ) 
     102 
     103if (meth_hcoup.eq.1) then 
     104   Hcoup(:,:)=coefbmshelf*Hcoup(:,:) 
     105   Hcoup(:,:)=min( max(Hcoup (:,:),Hcoup_plateau),Hcoup_abysses) 
     106else if (meth_hcoup.eq.2) then 
     107   Hcoup(:,:)=coefbmshelf*Hcoup(:,:) 
     108   Hcoup(:,:)=max(Hcoup(:,:),Hcoup_plateau) 
     109endif 
     110 
    75111 
    76112end subroutine init_calving 
    … …  
    93129 
    94130! coupure 
    95 if (meth_hcoup.eq.0) then 
    96    Hcoup=Hcoup_max 
    97 else if (meth_hcoup.eq.1) then 
    98    Hcoup=coefbmshelf*Hcoup_max 
    99    Hcoup=min(Hcoup,Hcoup_max) 
    100    Hcoup=max(Hcoup,Hcoup_min) 
    101  
     131Hcoup(:,:) = min ( max(                                         & 
     132     (-(Bsoc(:,:)-sealevel) - prof_plateau)/(prof_abysses-prof_plateau)    & 
     133     *(hcoup_abysses-hcoup_plateau)+hcoup_plateau               & 
     134     , hcoup_plateau), hcoup_abysses ) 
     135 
     136if (meth_hcoup.eq.1) then 
     137   Hcoup(:,:)=coefbmshelf*Hcoup(:,:) 
     138   Hcoup(:,:)=min( max(Hcoup (:,:),Hcoup_plateau),Hcoup_abysses) 
    102139else if (meth_hcoup.eq.2) then 
    103    Hcoup=coefbmshelf*Hcoup_max 
    104    Hcoup=max(Hcoup,Hcoup_min) 
     140   Hcoup(:,:)=coefbmshelf*Hcoup(:,:) 
     141   Hcoup(:,:)=max(Hcoup(:,:),Hcoup_plateau) 
    105142endif 
    106143 
    107144! hauteur au dessus de la coupure 
    108 hmhc(:,:)=H(:,:)-Hcoup 
     145hmhc(:,:)=H(:,:)-Hcoup(:,:) 
    109146 
    110147! coupure de l'ice shelf 
    … …  
    117154 
    118155 
    119 ifext:  if((flot(i,j)).and.(h(i,j).le.hcoup))  then 
     156ifext:  if((flot(i,j)).and.(h(i,j).le.hcoup(i,j)))  then 
    120157! ifext: pour les noeuds flottants avec h < hcoup  
    121158 
    … …  
    373410   if (testipj) testH=testH + abs(uybar(i,j+1))*H(i,j+1) 
    374411    
    375    if (testH.gt.Hcoup) then 
     412   if (testH.gt.Hcoup(i,j)) then 
    376413      testmij=.true. 
    377414   else 
    … …  
    446483! mab: copy the calving in the right format for coupling                      
    447484! calvin_GRIS =  volume flux per year (m3/a) 
    448                      calvin_GRIS(i,j)=calvin_GRIS(i,j)+(-calv(i,j)*dx*dy) 
     485                     calvingGRIS(i,j)=calvingGRIS(i,j)+(-calv(i,j)*dx*dy) 
    449486!                     calvin_GRIS(i,j)=-calv(i,j)*40000.*40000. 
    450487!                     write(15,*) 'calvheight:',h(i,j),calv(i,j),i,j,dtt,dt,calvin_GRIS(i,j)          

  • branches/iLoveclim/SOURCES/flottab2-0.7.f90

    r88 r122  
    104104    ! _________________________________________________________ 
    105105 
    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) 
     106 
    108107 
    109108    if (itracebug.eq.1)  call tracebug(' Entree dans routine flottab') 
    … …  
    172171 
    173172          archim = Bsoc(i,j)+H(i,j)*ro/row -sealevel 
    174           !      if ((i.eq.132).and.(j.eq.183)) print*,'archim=',archim 
    175173 
    176174 
    … …  
    182180                FLOT(I,J)=.true. 
    183181                BOOST=.false. 
     182 
     183                ! Attention :  avec le bloc dessous il faut faire un calcul de flottaison a la lecture  
    184184 
    185185                if (igrdline.eq.1) then   ! en cas de grounding line prescrite 
    … …  
    206206 
    207207 
     208          !else                           !    le point ne flotte pas 
    208209          else if ((H(i,j).gt.0.).and.(archim.GE.0.)) then    !    le point ne flotte pas et est englace 
    209210             mk(i,j)=0 
    … …  
    212213                FLOT(I,J)=.false. 
    213214                BOOST=.false. 
    214  
    215215             endif 
     216              
    216217             !cdc correction topo pour suivre  variations sealevel 
    217218             !cdd           S(i,j)=Bsoc(i,j)+H(i,j) 
    218219             S(i,j)=Bsoc(i,j)+H(i,j) 
    219220             B(i,j)=Bsoc(i,j) 
    220  
    221221 
    222222          else if  ((H(i,j).LE.0.).and.(archim.LT.0.)) then    !    terre deglace qui devient ocean  
    … …  
    226226             S(i,j)=surnet+sealevel 
    227227             B(i,j)=S(i,j)-H(i,j) 
     228              
    228229          endif arch 
    229230 
    … …  
    242243    end do 
    243244    !$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 !!!!!!!!!!!!!!!!' 
    247245 
    248246!!$ do i=1,nx 
    … …  
    318316 
    319317          if (flotmy(i,j).and.(new_flot_point(i,j).or.     & 
    320                new_flot_point(i,j-1))) then 
     318             new_flot_point(i,j-1))) then 
    321319             appel_new_flot=.true. 
    322320             new_flotmy(i,j)=.true. 
    … …  
    604602    !$OMP END PARALLEL 
    605603 
    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 
     604    call DETERMIN_TACHE  
     605 
     606    synchro :    if (isynchro.eq.1) then 
     607!!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf) 
     608!!$if (itestf.gt.0) then 
     609!!$   write(6,*) 'dans flottab avant DETERMIN_TACHE  asymetrie sur H  pour time=',time 
     610!!$   stop 
     611!!$else 
     612!!$   write(6,*) 'dans flottab apres  DETERMIN_TACHE pas d asymetrie sur H  pour time=',time 
     613!!$ 
     614!!$end if 
     615 
     616 
     617      !----------------------------------------------! 
     618      !On determine les differents ice strean/shelf  ! 
     619      !      call DETERMIN_TACHE                      ! 
     620      !----------------------------------------------! 
     621 
     622 
     623!!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf) 
     624!!$if (itestf.gt.0) then 
     625!!$   write(6,*) 'dans flottab apres DETERMIN_TACHE  asymetrie sur H  pour time=',time 
     626!!$   stop 
     627!!$else 
     628!!$   write(6,*) 'dans flottab apres DETERMIN_TACHE pas d asymetrie sur H  pour time=',time 
     629!!$ 
     630!!$end if 
     631 
     632      !ice(:,:)=0   Attention, voir si ca marche toujours pour l'Antarctique et heminord ! 
     633 
     634      !On compte comme englacé uniquement les calottes dont une partie est posée       
     635      !$OMP PARALLEL PRIVATE(smax_,smax_coord,smax_i,smax_j,mask_tache_ij) 
     636      !$OMP DO 
     637      do j=3,ny-2 
     638          do i=3,nx-2 
     639            test1:  if (.not.iceberg(table_out(i,j))) then ! on est pas sur un iceberg  
     640                if (nb_pts_tache(table_out(i,j)).ge.1) then 
     641                   ice(i,j)=1 
     642                   if (nb_pts_tache(table_out(i,j)).le.10) then  ! les petites iles sont en sia 
     643                      !    write(6,*) 'petite ile ',i,j 
     644                      flgzmx(i,j)=.false. 
     645                      flgzmx(i+1,j)=.false. 
     646                      flgzmy(i,j)=.false. 
     647                      flgzmy(i,j+1)=.false. 
     648                      gzmx(i:i+1,j)=.false. 
     649                      gzmy(i,j:j+1)=.false. 
     650                   endif 
     651 
     652 
     653                  ! ici on est sur une tache non iceberg >= 5 points                        
     654                  ! on teste si la tache n'est pas completement ice stream 
     655 
     656                  test2:       if (icetrim(table_out(i,j))) then   ! on a une ile d'ice stream 
     657 
     658                      mask_tache_ij(:,:)=.false. 
     659                      mask_tache_ij(:,:)=(table_out(:,:).eq.table_out(i,j))       ! pour toute la tache 
     660 
     661                      smax_=maxval(S(:,:),MASK=mask_tache_ij(:,:)) 
     662                      smax_coord(:)=maxloc(S(:,:),MASK=mask_tache_ij(:,:)) 
     663                      smax_i=smax_coord(1) 
     664                      smax_j=smax_coord(2) 
     665 
     666!!$               smax_i=0 ; smax_j=0 ; smax_=sealevel 
     667!!$               do ii=3,nx-2 
     668!!$                  do jj=3,ny-2 
     669!!$                     if (table_out(ii,jj).eq.table_out(i,j)) then 
     670!!$                        if (s(ii,jj).gt.smax_) then                            
     671!!$                           smax_ =s(ii,jj) 
     672!!$                           smax_i=ii 
     673!!$                           smax_j=jj 
     674!!$                        endif 
     675!!$                     endif 
     676!!$                  end do 
     677!!$               end do 
     678 
     679 
     680                      gzmx(smax_i,smax_j)=.false. ; gzmx(smax_i+1,smax_j)=.false. 
     681                      gzmy(smax_i,smax_j)=.false. ; gzmx(smax_i,smax_j+1)=.false. 
     682                      flgzmx(smax_i,smax_j)=.false. ; flgzmx(smax_i+1,smax_j)=.false. 
     683                      flgzmy(smax_i,smax_j)=.false. ; flgzmx(smax_i,smax_j+1)=.false. 
     684 
     685                      if (Smax_.le.sealevel) then 
     686                         write(num_tracebug,*)'Attention, une ile avec la surface sous l eau' 
     687                         write(num_tracebug,*)'time=',time,'   coord:',smax_i,smax_j 
     688                      end if 
     689 
     690                   endif test2 
     691                end if                                             ! endif deplace 
     692 
     693             else ! on est sur un iceberg                          !   test1 
     694                ice(i,j)=0 
     695                h(i,j)=1. 
     696                surnet=H(i,j)*(1.-ro/row) 
     697                S(i,j)=surnet+sealevel 
     698                B(i,j)=S(i,j)-H(i,j) 
     699 
     700             endif test1 
     701 
     702 
     703          end do 
     704      end do 
     705      !$OMP END DO 
     706      !$OMP END PARALLEL 
     707 
     708      !---------------------------------------------- 
     709      ! On caracterise le front des ice shelfs/streams 
     710 
     711      !      call DETERMIN_FRONT     
     712 
     713      !----------------------------------------------       
     714!!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf) 
     715!!$if (itestf.gt.0) then 
     716!!$   write(6,*) 'dans flottab apres DETERMIN_front  asymetrie sur H  pour time=',time 
     717!!$   stop 
     718!!$else 
     719!!$   write(6,*) 'dans flottab apres DETERMIN_front pas d asymetrie sur H  pour time=',time 
     720!!$ 
     721!!$end if 
     722 
     723    endif synchro 
    726724 
    727725    ! correction momentanée pour symetrie Heino 
    … …  
    730728    !fin de routine flottab2 
    731729    !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,*) 
     730 
    735731  end subroutine flottab 
    736732!--------------------------------------------------------------------   

  • branches/iLoveclim/SOURCES/main3D-0.4-40km.f90

    r91 r122  
    131131  PRINT*,'TIME = ',TIME,' TEND = ',TEND 
    132132 
     133!afq -- reset the GRISLI fluxes to CLIO: 
     134     calvin_GRIS(:,:) = 0.   
     135     runof_oc(:,:) = 0.  
     136     bmelt_oc(:,:) = 0. 
     137      
    133138  time_loop: DO WHILE (time.LT.tend)  !____________________________ debut boucle temporelle 
    134139     call step_time_loop 
    … …  
    137142  end do time_loop 
    138143 
     144!afq -- reset the CLIO fluxes to GRISLI: 
    139145  bmshelfCLIO(:,:,:) = 0d0 
     146!afq -- put the right units for GRISLI fluxes to CLIO: 
     147     calvin_GRIS(:,:) = calvin_GRIS(:,:)/timCplGRISyr 
     148     runof_oc(:,:) = runof_oc(:,:)/timCplGRISyr 
     149     bmelt_oc(:,:) = bmelt_oc(:,:)/timCplGRISyr 
    140150 
    141151  if (itracebug.eq.1)  call tracebug('dans main avant call out_recovery ') 
Note: See TracChangeset for help on using the changeset viewer.