Changeset 90 for branches


Ignore:
Timestamp:
10/14/16 14:45:46 (8 years ago)
Author:
aquiquet
Message:

GRISLI-loveclim branch, GRISLI freshwater flux needed by iLOVECLIM

Location:
branches/iLoveclim/SOURCES
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • branches/iLoveclim/SOURCES/3D-physique-gen_mod.f90

    r89 r90  
    251251  !hassine 
    252252  real, dimension(nx,ny) :: beta_centre !< beta on major node (average) 
     253  real (kind = 8),dimension(nx,ny) :: calvingGRIS 
    253254  real (kind = 8),dimension(nx,ny),target :: calvin_GRIS ! calving flux[m3/s] given to CLIO 
    254255  real, dimension(nx,ny) :: coef_drag   !< coefficient de la loi de friction non lineaire : depend de la valeur de alpha_drag 

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

    r77 r90  
    2727!real ::  pdfmean                      ! moyenne PDF 
    2828logical,dimension(nx,ny,13) :: mask_cal ! masque regions calotte  
     29REAL, dimension(nx,ny) :: old_H_dtt, old_H_new    ! Epaisseur de glace au pas de temps precedent 
     30real, dimension(nx,ny) :: runof_oc_sngl 
    2931integer, dimension(nx,ny) :: write_mask 
    3032CONTAINS 
    … …  
    3638mask_cal(:,:,:)=.false. 
    3739mask_cal(:,:,1)=.true. 
     40old_H_dtt(:,:)=H(:,:) ! afq missing in mab? 
    3841! calcul d'un masque pour les regions des calottes 
    3942do j=1,ny 
    … …  
    178181!      REAL ABLA(NX,NY) 
    179182      REAL DELTAVOL 
     183      REAL, dimension(nx,ny,13) :: delta_H_dtt 
     184      REAL, dimension(nx,ny) :: delta_H_test 
    180185 
    181186 
    … …  
    207212        TBM(kk) = 0. 
    208213        ITJJA(kk) = 0. 
     214        delta_H_dtt(:,:,kk) = 0. 
     215        runof_oc(:,:) = 0. 
    209216         
    210217! nouveau tof mai 2009 
    … …  
    231238 
    232239 
     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 
     244     endwhere 
     245     old_H_new(:,:) = delta_H_dtt(:,:,1) 
     246     delta_H_dtt(:,:,1) = old_H_dtt(:,:) - delta_H_dtt(:,:,1) 
     247      
     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      
    233258!!$ 
    234259!!$ 

  • branches/iLoveclim/SOURCES/calving_frange.f90

    r4 r90  
    2424real :: hcoup_max  ! epaisseur max 
    2525real :: hcoup_min=100.  ! epaisseur min 
     26real :: mass_total 
    2627integer :: meth_hcoup  ! pour avoir hcoup dépendant du coefbmshelf 
    2728integer :: ifrange        !  0 pas de traitement particulier pres du bord, 1 -> franges 
    … …  
    3132logical :: avalw,avale,avals,avaln,interieur 
    3233 
     34integer :: testH ! somme des flux entrants en un point ij 
     35 
    3336contains 
    3437 
    … …  
    4043 
    4144 
     45mass_total = 0d0 
     46calvingGRIS=0d0 
    4247calv(:,:)=0.  
    4348 
    … …  
    7277!--------------------------------------------------------------------------------------- 
    7378subroutine calving 
    74  
     79   
     80  real :: calvtest 
     81   
    7582! initialisation calving 
    7683  calv(:,:)=0.  
     84  mass_total=0d0 
     85  calvtest=0. 
     86  mass_total=0d0 
    7787 
    7888! calcul du dhdt lagrangien 
    … …  
    151161! ifrange = 7             calving drastique sauf si coefbmshelf < 0.5  
    152162 
     163!cdc  
     164! ifrange=44              idem 4 mais avec test sur i et j couple // afq: probably useless 
     165! ifrange=8               C. Dumas : coupure si alimentation ne permet pas H > Hcoup 
    153166 
    154167if (ifrange.eq.1) then           !   Rajout vince : si un point voisin est pose -> test vrai. 
    … …  
    315328        .and.(coefbmshelf.lt.1.)) 
    316329 
    317  
     330else if (ifrange.eq.44) then       ! idem 4 mais avec test sur i et j couple 
     331 
     332   bilan_surf_fond=(bm(i,j)-bmelt(i,j).gt.0.) 
     333 
     334   testmij=( ((hmhc(i-1,j).gt.0.).and.(uxbar(i,j).ge.0.)  &  ! voisin (i-1,j) amont et > hcoup 
     335        .and.  (dhdt(i-1,j).gt.(-hmhc(i-1,j)*abs(uxbar(i-1,j)/dx)))) &  
     336        .or.((.not.flot(i-1,j)).and.bilan_surf_fond )) ! 
     337 
     338   testpij=( ((hmhc(i+1,j).gt.0.).and.(uxbar(i+1,j).le.0.) & ! voisin (i+1,j) amont et > hcoup 
     339        .and.(dhdt(i+1,j).gt.(-hmhc(i+1,j)*abs(uxbar(i+1,j)/dx)))) & 
     340        .or.((.not.flot(i+1,j)).and.bilan_surf_fond) ) ! 
     341 
     342   testimj=( ((hmhc(i,j-1).gt.0.).and.(uybar(i,j).ge.0.)  &  ! voisin (i,j-1) amont et > hcoup 
     343        .and.(dhdt(i,j-1).gt.(-hmhc(i,j-1)*abs(uybar(i,j-1)/dx))))& 
     344        .or.((.not.flot(i,j-1)).and.bilan_surf_fond ) ) ! 
     345 
     346   testipj=( ((hmhc(i,j+1).gt.0.).and.(uybar(i,j+1).le.0.) & ! voisin (i,j+1) amont et > hcoup 
     347        .and.(dhdt(i,j+1).gt.(-hmhc(i,j+1)*abs(uybar(i,j+1)/dx))))& 
     348        .or.((.not.flot(i,j+1)).and.bilan_surf_fond ) ) ! 
     349 
     350else if (ifrange.eq.8) then ! test C Dumas : si somme flux entrant permet H > Hcoup : vrai  
     351                   !   bilan_surf_fond=(bm(i,j)-bmelt(i,j).gt.0.) 
     352 
     353   testmij=( ((hmhc(i-1,j).gt.0.).and.(uxbar(i,j).ge.0.))  &  ! voisin (i-1,j) amont et > hcoup 
     354                        !        .and.  (dhdt(i-1,j).gt.(-hmhc(i-1,j)*abs(uxbar(i-1,j)/dx)))) &  
     355        .or.((.not.flot(i-1,j)).and.bilan_surf_fond )) ! 
     356 
     357   testpij=( ((hmhc(i+1,j).gt.0.).and.(uxbar(i+1,j).le.0.)) & ! voisin (i+1,j) amont et > hcoup 
     358                        !        .and.(dhdt(i+1,j).gt.(-hmhc(i+1,j)*abs(uxbar(i+1,j)/dx)))) & 
     359        .or.((.not.flot(i+1,j)).and.bilan_surf_fond) ) ! 
     360 
     361   testimj=( ((hmhc(i,j-1).gt.0.).and.(uybar(i,j).ge.0.))  &  ! voisin (i,j-1) amont et > hcoup 
     362                        !        .and.(dhdt(i,j-1).gt.(-hmhc(i,j-1)*abs(uybar(i,j-1)/dx))))& 
     363        .or.((.not.flot(i,j-1)).and.bilan_surf_fond ) ) ! 
     364 
     365   testipj=( ((hmhc(i,j+1).gt.0.).and.(uybar(i,j+1).le.0.)) & ! voisin (i,j+1) amont et > hcoup 
     366                        !        .and.(dhdt(i,j+1).gt.(-hmhc(i,j+1)*abs(uybar(i,j+1)/dx))))& 
     367        .or.((.not.flot(i,j+1)).and.bilan_surf_fond ) ) ! 
     368 
     369   testH=0. 
     370   if (testmij) testH=testH + abs(uxbar(i,j))*H(i-1,j) 
     371   if (testpij) testH=testH + abs(uxbar(i+1,j))*H(i+1,j) 
     372   if (testimj) testH=testH + abs(uybar(i,j))*H(i,j-1) 
     373   if (testipj) testH=testH + abs(uybar(i,j+1))*H(i,j+1) 
     374    
     375   if (testH.gt.Hcoup) then 
     376      testmij=.true. 
     377   else 
     378      testmij=.false. 
     379      testpij=.false. 
     380      testimj=.false. 
     381      testipj=.false. 
     382   endif 
     383                    
    318384endif 
    319385 
    … …  
    326392iin3=max(1,i-3) 
    327393 
    328 avalw=((hmhc(i-1,j).gt.0.).or.(hmhc(iin2,j).gt.0.) & 
    329      .or.(hmhc(iin3,j).gt.0.)).and.(uxbar(i,j).le.0.)                          
    330  
     394!cdc avalw=((hmhc(i-1,j).gt.0.).or.(hmhc(iin2,j).gt.0.) & 
     395!cdc      .or.(hmhc(iin3,j).gt.0.)).and.(uxbar(i,j).le.0.) 
     396 
     397avalw=((hmhc(i-1,j).gt.0..and.flot(i-1,j)).or.(hmhc(iin2,j).gt.0..and.flot(iin2,j)) & 
     398     .or.(hmhc(iin3,j).gt.0..and.flot(iin3,j))).and.(uxbar(i,j).le.0.)    
    331399 
    332400! 1 des 3 voisins aval > hcoup du cote est (i+1) 
    … …  
    334402iin3=min(i+3,nx) 
    335403 
    336 avale=((hmhc(i+1,j).gt.0.).or.(hmhc(iin2,j).gt.0.) &     
    337      .or.(hmhc(iin3,j).gt.0.)).and.(uxbar(i+1,j).ge.0)                           
    338  
     404!cdc avale=((hmhc(i+1,j).gt.0.).or.(hmhc(iin2,j).gt.0.) &     
     405!cdc      .or.(hmhc(iin3,j).gt.0.)).and.(uxbar(i+1,j).ge.0)                           
     406    
     407avale=((hmhc(i+1,j).gt.0..and.flot(i+1,j)).or.(hmhc(iin2,j).gt.0..and.flot(iin2,j)) &     
     408     .or.(hmhc(iin3,j).gt.0..and.flot(iin3,j))).and.(uxbar(i+1,j).ge.0) 
     409                 
    339410! 1 des 3 voisins aval > hcoup du cote sud (j-1) 
    340411jin2=max(1,j-2) 
    341412jin3=max(1,j-3) 
    342413 
    343 avals=((hmhc(i,j-1).gt.0.).or.(hmhc(i,jin2).gt.0.) &     
    344      .or.(hmhc(i,jin3).gt.0.)).and.(uybar(i,j).le.0.)   
     414!cdc avals=((hmhc(i,j-1).gt.0.).or.(hmhc(i,jin2).gt.0.) &     
     415!cdc     .or.(hmhc(i,jin3).gt.0.)).and.(uybar(i,j).le.0.)   
    345416                           
     417avals=((hmhc(i,j-1).gt.0..and.flot(i,j-1)).or.(hmhc(i,jin2).gt.0..and.flot(i,jin2)) &     
     418     .or.(hmhc(i,jin3).gt.0..and.flot(i,jin3))).and.(uybar(i,j).le.0.) 
     419                 
    346420! 1 des 3 voisins aval > hcoup du cote nord (j-1) 
    347421jin2=min(j+2,ny) 
    348422jin3=min(j+3,ny) 
    349423 
    350 avaln=((hmhc(i,j+1).gt.0.).or.(hmhc(i,jin2).gt.0.) &     
    351      .or.(hmhc(i,jin2).gt.0.)) .and.(uybar(i,j+1).ge.0.)                          
     424!cdc avaln=((hmhc(i,j+1).gt.0.).or.(hmhc(i,jin2).gt.0.) &     
     425!cdc     .or.(hmhc(i,jin2).gt.0.)) .and.(uybar(i,j+1).ge.0.)                          
     426 
     427avaln=((hmhc(i,j+1).gt.0..and.flot(i,j+1)).or.(hmhc(i,jin2).gt.0..and.flot(i,jin2)) &     
     428     .or.(hmhc(i,jin2).gt.0..and.flot(i,jin2))) .and.(uybar(i,j+1).ge.0.) 
    352429 
    353430interieur=(avalw.or.avale).and.(avals.or.avaln) 
    … …  
    355432 
    356433 
    357  if ((.not.(testmij.or.testpij.or.testimj.or.testipj))  &  ! pas suffisament alimente 
    358                 .and.(.not.interieur)) then                ! et pas interieur 
     434!cdc if ((.not.(testmij.or.testpij.or.testimj.or.testipj))  &  ! pas suffisament alimente 
     435!cdc                .and.(.not.interieur)) then                ! et pas interieur 
     436 
     437!cdc pour ne pas faire de calving si H > 350m :   if ((.not.(testmij.or.testpij.or.testimj.or.testipj)).and.(H(i,j).LT.350)) then 
     438! mab: if none of the neighbours doesn't provide ice for the considered point, 
     439! mab: then the ice is cut at this point and given to calv(i,j) 
     440                if ((.not.(testmij.or.testpij.or.testimj.or.testipj))) then 
     441! mab: new version so that grid points with h=1m (all points were calving takes place are afterwards 
     442! mab: set to a height of 1m) are not taken into account,  provides a problem as soon as you start  
     443! mab: coupling GRISLI to the iceberg modell 
     444                   if(h(i,j).gt.1) then 
    359445                     calv(i,j)=-h(i,j)   
    360                      H(i,j)=1              
     446! mab: copy the calving in the right format for coupling                      
     447! calvin_GRIS =  volume flux per year (m3/a) 
     448                     calvin_GRIS(i,j)=calvin_GRIS(i,j)+(-calv(i,j)*dx*dy) 
     449!                     calvin_GRIS(i,j)=-calv(i,j)*40000.*40000. 
     450!                     write(15,*) 'calvheight:',h(i,j),calv(i,j),i,j,dtt,dt,calvin_GRIS(i,j)          
    361451                     endif   
     452                   H(i,j)=1              
     453                endif 
    362454 
    363455                  end if ifint 
    … …  
    365457         end do 
    366458      end do 
     459          
     460    calvtest=sum(-calv(:,:)*dx*dy) 
     461    mass_total=sum(calvingGRIS(:,:)) 
     462 
     463    calvtest = 0. 
     464 
     465     
    367466    end subroutine calving 
    368467!------------------------------------------------------------------------------------------ 
Note: See TracChangeset for help on using the changeset viewer.