Changeset 237 for trunk


Ignore:
Timestamp:
01/09/19 17:09:26 (5 years ago)
Author:
aquiquet
Message:

Sealevel is now treated as a 2D variable (sealevel_2d while sealevel remains the eustatic sea level), results should remain identical as sealevel_2d is equal to sealevel in this revision.

Location:
trunk/SOURCES
Files:
51 edited

Legend:

Unmodified
Added
Removed
  • trunk/SOURCES/3D-physique-gen_mod.f90

    r166 r237  
    141141  real ::  SECYEAR                      !< for relation an/seconds 
    142142  !real ::  S22                          ! ct for PDD calculation 
    143   real ::  sealevel                     !< niveau des mers  
     143  real ::  sealevel                     !< niveau des mers, afq: now used as eustatic sea level 
    144144  ! real ::  SIGMA                        ! variabilite Tday 
    145145  real ::  SURF                         !<  
     
    235235  integer,dimension(nx,ny) :: Mk_init      !< initial mask (with islands, outcrops, ... 
    236236  integer,dimension(nx,ny) :: MK           !< masks (ice sheet, max, above water, below water, 1) 
    237   integer,dimension(nx,ny) :: M_sealev     !< masks for sea level calculation (present ice sheet) 
    238237  integer,dimension(nx,ny) :: MNEG         !< masks (ice sheet, max, above water, below water, 1) 
    239238  integer,dimension(nx,ny) :: IBASE     !< type de base (froide, temperee) 
     
    245244  real,dimension(nx,ny) :: ACQUA        !< Surface des surfaces en eau 
    246245  real,dimension(nx,ny) :: B            !< Altitude de la base de la glace  'o' 
    247   real,dimension(nx,ny) :: B_sealev     !< Altitude de la base de la glace pour le 
    248   !< calcul du nivx de la mer (/present) 
    249246  real,dimension(nx,ny) :: BDOT         !< derivee de B / t 
    250247  real,dimension(nx,ny) :: B1           !<  
     
    310307 
    311308 
    312   real,dimension(nx,ny) :: H_sealev     !< ice thickness for sea level calculation        ! (present ice sheet) 
    313309  real,dimension(nx,ny) :: HDOT         !< ice thickness derivee / t 
    314310  real,dimension(nx,ny) :: HDOTWATER 
     
    334330  real,dimension(nx,ny) :: SW           !< for bedrock isostasy 
    335331  real,dimension(nx,ny) :: S            !< altitude of ice sheet surface 
    336   real,dimension(nx,ny) :: S_sealev     !< altitude of ice sheet surface for sea 
    337   !< level calculation (present ice sheet) 
     332  real,dimension(nx,ny) :: sealevel_2d  !< local sea surface elevation 
    338333  real,dimension(nx,ny) :: SLOPE        !< slope 'o' 
    339334  real,dimension(nx,ny) :: SDX          !< slope derivee / x '>' 
     
    344339  real,dimension(nx,ny) :: SLOPE2my     !< = Sdy**2 + Sdxmy**2 '^' 
    345340  real,dimension(nx,ny) :: S0           !< altitude actuelle de la surface 
    346   real,dimension(nx,ny) :: slv          !< niveau de flottaison (sealevel et lakes) 
     341!afq -- replaced by sealevel_2D:  real,dimension(nx,ny) :: slv          !< niveau de flottaison (sealevel et lakes) 
    347342  real,dimension(nx,ny) :: TJULY        !< Ground air temperature July 
    348343  real,dimension(nx,ny) :: TANN         !< Ground air temperature annual 
  • trunk/SOURCES/Ant15_CISM_files/output_anta_mod-0.4.f90

    r4 r237  
    8383 
    8484!        calcul de la hauteur au dessus de la flottaison 
    85          if (sealevel-B(i,j).le.0.) then    ! socle au dessus du niveau des mers 
     85         if (sealevel_2d(i,j)-B(i,j).le.0.) then    ! socle au dessus du niveau des mers 
    8686               volf=volf+h(i,j) 
    8787         else 
    88             volf=volf+h(i,j)-row/ro*(sealevel-b(i,j)) 
     88            volf=volf+h(i,j)-row/ro*(sealevel_2d(i,j)-b(i,j)) 
    8989         endif 
    9090 
  • trunk/SOURCES/Ant16_files/climat_InitMIP_years_perturb_mod.f90

    r180 r237  
    1616 
    1717 
    18 use module3d_phy,only: nx,ny,S,S0,H,Bsoc,acc,abl,BM,Tann,Tjuly,Ts,time,dt,num_param,num_rep_42,num_forc,dirforcage,dirnameinp,tafor,time,sealevel,coefbmshelf 
     18use module3d_phy,only: nx,ny,S,S0,H,Bsoc,acc,abl,BM,Tann,Tjuly,Ts,time,dt,num_param,num_rep_42,num_forc,dirforcage,dirnameinp,tafor,time,sealevel,sealevel_2d,coefbmshelf 
    1919use netcdf 
    2020use io_netcdf_grisli 
     
    285285100 continue 
    286286 
     287  sealevel_2d(:,:) = sealevel 
     288 
    287289  Tann (:,:) = Ta0 (:,:) + T_lapse_rate * (S(:,:)-S0(:,:)) +Tafor 
    288290  Ts(:,:)    = Tann(:,:) 
  • trunk/SOURCES/Ant40_files/lect-anteis_mod.f90

    r227 r237  
    228228    do J=1,NY 
    229229       do I=1,NX 
    230           if ((BSOC(I,J)+H(I,J)*RO/ROW -SEALEVEL).LT.0.) then 
     230          if ((BSOC(I,J)+H(I,J)*RO/ROW -sealevel_2d(i,j)).LT.0.) then 
    231231             FLOT(I,J)=.TRUE. 
    232232          else 
     
    260260    ymax=ycc(nx,ny)/1000. 
    261261 
    262 !!$    !     lecture du fichier de reference pour le calcul du niveau des mers : etat actuel 
    263 !!$    open(88,file=TRIM(DIRNAMEINP)//'grzone56-k000.pl') 
    264 !!$ 
    265 !!$    read(88,*) 
    266 !!$    read(88,*) 
    267 !!$    read(88,*) 
    268 !!$    do j=1,ny 
    269 !!$       do i=1,nx 
    270 !!$          read(88,*) S_sealev(i,j),H_sealev(i,j),B_sealev(i,j),M_sealev(i,j) 
    271 !!$       enddo 
    272 !!$    enddo 
    273 !!$    close(88) 
    274 !!$ 
    275 !!$    write(num_rep_42,*) 'fichier reference pour le niveau des mers : ', &  
    276 !!$         TRIM(DIRNAMEINP)//'grzone56-k000.pl' 
    277  
    278262    ! lecture du flux geothermique de Shapiro 
    279263    open(88,file=ghf_fich) 
  • trunk/SOURCES/Ant40_files/output_anta40_mod-0.4.f90

    r148 r237  
    191191 
    192192!         calcul de la hauteur au dessus de la flottaison 
    193             if (sealevel-B(i,j).le.0.) then    ! socle au dessus du niveau des mers 
     193            if (sealevel_2d(i,j)-B(i,j).le.0.) then    ! socle au dessus du niveau des mers 
    194194              volf=volf+h(i,j)*corrsurf(i,j)                          ! volume au-dessus de la flottaison 
    195195            else 
    196               volf=volf+(h(i,j)-row/ro*(sealevel-b(i,j)))*corrsurf(i,j) ! volume au-dessus de la flottaison 
     196              volf=volf+(h(i,j)-row/ro*(sealevel_2d(i,j)-b(i,j)))*corrsurf(i,j) ! volume au-dessus de la flottaison 
    197197            endif 
    198198 
  • trunk/SOURCES/Ant45_CISM_files/output_anta_mod-0.4.f90

    r4 r237  
    8484 
    8585!        calcul de la hauteur au dessus de la flottaison 
    86          if (sealevel-B(i,j).le.0.) then    ! socle au dessus du niveau des mers 
     86         if (sealevel_2d(i,j)-B(i,j).le.0.) then    ! socle au dessus du niveau des mers 
    8787               volf=volf+h(i,j) 
    8888         else 
    89             volf=volf+h(i,j)-row/ro*(sealevel-b(i,j)) 
     89            volf=volf+h(i,j)-row/ro*(sealevel_2d(i,j)-b(i,j)) 
    9090         endif 
    9191 
  • trunk/SOURCES/Antarctique_general_files/output_anta_mod-0.4.f90

    r68 r237  
    192192              
    193193             !        calcul de la hauteur au dessus de la flottaison 
    194              if (sealevel-B(i,j).le.0.) then    ! socle au dessus du niveau des mers 
     194             if (sealevel_2d(i,j)-B(i,j).le.0.) then    ! socle au dessus du niveau des mers 
    195195                volf=volf+h(i,j) 
    196196             else 
    197                 volf=volf+h(i,j)-row/ro*(sealevel-b(i,j)) 
     197                volf=volf+h(i,j)-row/ro*(sealevel_2d(i,j)-b(i,j)) 
    198198             endif 
    199199              
     
    289289                scal_np         =  scal_np + 1 
    290290                scal_vol        =  scal_vol + H(i,j) 
    291                 scal_volbuoy    =  scal_volbuoy + H(i,j) - row/ro * max(0.,sealevel-b(i,j)) 
     291                scal_volbuoy    =  scal_volbuoy + H(i,j) - row/ro * max(0.,sealevel_2d(i,j)-b(i,j)) 
    292292                scal_meanhdot   =  scal_meanhdot + Hdot(i,j) 
    293293                scal_sigma_hdot =  scal_sigma_hdot + Hdot(i,j)**2 
  • trunk/SOURCES/Draggings_modules/dragging_prescr_beta_buoyency_mod.f90

    r70 r237  
    500500 
    501501 
    502     where ( sealevel-B1(:,:).le.0.)             ! socle au dessus du niveau des mers 
     502    where ( sealevel_2d(:,:)-B1(:,:).le.0.)             ! socle au dessus du niveau des mers 
    503503       H_buoy(:,:) = H1(:,:) 
    504504 
    505505    elsewhere 
    506        H_buoy(:,:) = H1(:,:)-row/ro*(sealevel-B1(:,:)) 
     506       H_buoy(:,:) = H1(:,:)-row/ro*(sealevel_2d(:,:)-B1(:,:)) 
    507507        
    508508    end where 
  • trunk/SOURCES/Draggings_modules/dragging_prescr_beta_nolin_mod.f90

    r168 r237  
    606606 
    607607 
    608     where ( sealevel-B1(:,:).le.0.)             ! socle au dessus du niveau des mers 
     608    where ( sealevel_2d(:,:)-B1(:,:).le.0.)             ! socle au dessus du niveau des mers 
    609609       H_buoy(:,:) = H1(:,:) 
    610610 
    611611    elsewhere 
    612        H_buoy(:,:) = H1(:,:)-row/ro*(sealevel-B1(:,:)) 
     612       H_buoy(:,:) = H1(:,:)-row/ro*(sealevel_2d(:,:)-B1(:,:)) 
    613613        
    614614    end where 
  • trunk/SOURCES/Eurasie40_files/bmelt-eurasie-depth-lake_mod.f90

    r4 r237  
    8585         do i=1,nx 
    8686 
    87             if ((sealevel-BSOC(i,j)).lt.450) then 
     87            if ((sealevel_2d(i,j)-BSOC(i,j)).lt.450) then 
    8888                    bmshelf(i,j)=2.00 !3.00! 2.55! 0.45!   0.65   
    8989                    bmgrz(i,j)=2.00 !3.00! 2.55! 0.45 !  0.65 
     
    127127  
    128128 
    129                    if (slv(i,j).gt.0) then 
     129                   if (sealevel_2d(i,j).gt.0) then 
    130130                         bmshelf(i,j)=0.05!5.0     
    131131                         bmgrz(i,j)=0.05!5.0 
     
    135135         
    136136                      
    137              if (slv(i,j).gt.sealevel) then 
     137             if (sealevel_2d(i,j).gt.sealevel) then 
    138138                     bmelt(i,j)=bmshelf(i,j) 
    139139                     if (fbm(i,j))  bmelt(i,j)=bmgrz(i,j) 
     
    166166!   en fonction du nombre de points flottants 
    167167 
    168            if (slv(i,j).gt.sealevel) then 
     168           if (sealevel_2d(i,j).gt.sealevel) then 
    169169                 bmelt(i,j)= ngr/4.*bmgrz(i,j)+(1.-ngr/4.)*bmelt(i,j)          
    170170           else 
  • trunk/SOURCES/Eurasie40_files/climat-forcage-eurasie_mod-0.4.f90

    r4 r237  
    493493   nom_table='lvflo' 
    494494   call printtable_r(levelflot,nom_table) 
    495      
     495    
     496   sealevel_2d(:,:) = sealevel 
     497 
    496498!   coefbmshelf(:,:)=(1.+delTatime(:,:)/7.)   
    497499!coefbmshelf=1. ! modif pour test de sensibilite (tof 20 avril 01) 
  • trunk/SOURCES/Eurasie40_files/lakes-prescribed_mod-0.1.f90

    r4 r237  
    5454subroutine lake_to_slv 
    5555 
    56 ! inclut les lacs dans le tableau slv des niveaux de flottaison 
     56! inclut les lacs dans le tableau sealevel_2d des niveaux de flottaison 
    5757 
    58 slv(:,:)=sealevel 
     58sealevel_2d(:,:)=sealevel 
    5959 
    6060! VOIR COMMENT ACTIVER-DESACTIVER LES LACS EN FONCTION DU TEMPS 
     
    6464      do i_lake=1,nb_lakes 
    6565         if (lake(i,j).eq.i_lake) then  
    66             slv(i,j)=max(level_lakes(i_lake),Bsoc(i,j)) 
     66            sealevel_2d(i,j)=max(level_lakes(i_lake),Bsoc(i,j)) 
    6767         end if 
    6868      end do 
  • trunk/SOURCES/Eurasie40_files/output_eurasie40_mod-0.1.f90

    r4 r237  
    8585 
    8686!        calcul de la hauteur au dessus de la flottaison 
    87          if (sealevel-B(i,j).le.0.) then    ! socle au dessus du niveau des mers 
     87         if (sealevel_2d(i,j)-B(i,j).le.0.) then    ! socle au dessus du niveau des mers 
    8888               volf=volf+h(i,j) 
    8989         else 
    90             volf=volf+h(i,j)-row/ro*(sealevel-b(i,j)) 
     90            volf=volf+h(i,j)-row/ro*(sealevel_2d(i,j)-b(i,j)) 
    9191         endif 
    9292 
  • trunk/SOURCES/GrIce2sea_files/climat_GrIce2sea_years_perturb_mod.f90

    r225 r237  
    1616 
    1717 
    18 use module3d_phy,only: nx,ny,S,S0,H,Bsoc,acc,abl,BM,Tann,Tjuly,Ts,time,dt,num_param,num_rep_42,num_forc,dirforcage,dirnameinp,tafor,time,sealevel,coefbmshelf 
     18use module3d_phy,only: nx,ny,S,S0,H,Bsoc,acc,abl,BM,Tann,Tjuly,Ts,time,dt,num_param,num_rep_42,num_forc,dirforcage,dirnameinp,tafor,time,sealevel,sealevel_2d,coefbmshelf 
    1919!use lect_climref_Ice2sea 
    2020use netcdf 
     
    382382        endif 
    383383100     continue 
     384 
     385        sealevel_2d(:,:) = sealevel 
    384386 
    385387        Tann (:,:) = Ta0 (:,:) + T_lapse_rate * (S(:,:)-S0(:,:)) +Tafor 
  • trunk/SOURCES/GrIce2sea_files/lect_GrIce2sea_gen_nc.f90

    r4 r237  
    207207 
    208208 
    209     !     lecture du fichier de reference pour le calcul du niveau des mers : etat actuel 
    210     !     pas de version correcte pour l'instant -> valeurs initiales 
    211     S_sealev(:,:) = S0(:,:) 
    212     H_sealev(:,:) = H0(:,:) 
    213     B_sealev(:,:) = Bsoc(:,:) 
    214     M_sealev(:,:) = Mk0(:,:) 
    215  
    216209    if (itracebug.eq.1)  call tracebug(' apres lect entete') 
    217210   
  • trunk/SOURCES/GrIce2sea_files/output_Grice2sea_mod.f90

    r11 r237  
    196196              
    197197             !        calcul de la hauteur au dessus de la flottaison 
    198              if (sealevel-B(i,j).le.0.) then    ! socle au dessus du niveau des mers 
     198             if (sealevel_2d(i,j)-B(i,j).le.0.) then    ! socle au dessus du niveau des mers 
    199199                volf=volf+h(i,j) 
    200200             else 
    201                 volf=volf+h(i,j)-row/ro*(sealevel-b(i,j)) 
     201                volf=volf+h(i,j)-row/ro*(sealevel_2d(i,j)-b(i,j)) 
    202202             endif 
    203203              
     
    298298                scal_np         =  scal_np + 1 
    299299                scal_vol        =  scal_vol + H(i,j) 
    300                 scal_volbuoy    =  scal_volbuoy + H(i,j) - row/ro * max(0.,sealevel-b(i,j)) 
     300                scal_volbuoy    =  scal_volbuoy + H(i,j) - row/ro * max(0.,sealevel_2d(i,j)-b(i,j)) 
    301301                scal_meanhdot   =  scal_meanhdot + Hdot(i,j) 
    302302                scal_sigma_hdot =  scal_sigma_hdot + Hdot(i,j)**2 
  • trunk/SOURCES/Greeneem_files/lect-clim-act-greeneem_mar_mod.f90

    r4 r237  
    5555   do i=1,nx 
    5656      read (20,*) S0CLIM(i,j) 
    57       S0CLIM(i,j)=max(S0CLIM(i,j),slv(i,j))   ! pour etre au niveau des mers 
     57      S0CLIM(i,j)=max(S0CLIM(i,j),sealevel_2d(i,j))   ! pour etre au niveau des mers 
    5858   end do 
    5959end do 
     
    146146  do j=1,ny 
    147147     do i=1,nx 
    148         alti=max(slv(i,j),S(I,J)) 
     148        alti=max(sealevel_2d(i,j),S(I,J)) 
    149149        Tann(i,j)=TA0(i,j)+lapse_ann*(alti-S0CLIM(i,j)) 
    150150        Tjuly(i,j)=TJ0(i,j)+lapse_july*(alti-S0CLIM(i,j)) 
  • trunk/SOURCES/Greeneem_files/lect-clim-act-greeneem_mois_lapsecouche_mod.f90

    r4 r237  
    114114            Tm_0(i,j,mo)=Tm(i,j,mo)      ! temp sur la calotte du GCM a l'instant initial 
    115115             
    116             zel=max(slv(i,j),S0CLIM(i,j)) 
     116            zel=max(sealevel_2d(i,j),S0CLIM(i,j)) 
    117117! on calcule une temperature au sol a partir de la temperature sur la calotte de ref du GCM 
    118118 
     
    131131            end if 
    132132 
    133             !zel=max(slv(i,j),S(i,j))     
     133            !zel=max(sealevel_2d(i,j),S(i,j))     
    134134            zel=max(0.,S(i,j))     
    135135! puis on calcule la temperature sur la calotte S 
  • trunk/SOURCES/Greeneem_files/lect-clim-act-greeneem_mois_mod.f90

    r4 r237  
    7474      do i=1,nx 
    7575         read (20,*) S0CLIM(i,j) 
    76          S0CLIM(i,j)=max(S0CLIM(i,j),slv(i,j))   ! pour etre au niveau des mers 
     76         S0CLIM(i,j)=max(S0CLIM(i,j),sealevel_2d(i,j))   ! pour etre au niveau des mers 
    7777      end do 
    7878   end do 
     
    120120      do i=1,nx 
    121121! Zs est l'altitude de la surface qu'elle soit mer, glace ou lac 
    122          ZS(i,j)=max(slv(i,j),S0CLIM(I,J)) !slv=sealev pour nolake 
     122         ZS(i,j)=max(sealevel_2d(i,j),S0CLIM(I,J)) !sealevel_2d=sealev pour nolake 
    123123 
    124124         do mo=1,mois 
     
    261261do j=1,ny 
    262262   do i=1,nx 
    263       ZS(i,j)=max(slv(i,j),S(i,j)) 
     263      ZS(i,j)=max(sealevel_2d(i,j),S(i,j)) 
    264264      Tann(i,j)=d_ann+gamma_ann*S(i,j)+c_ann*ylat(i,j)+kappa_ann*(360-xlong(i,j)) 
    265265      Tjuly(i,j)=d_july+gamma_july*S(i,j)+c_july*ylat(i,j)+kappa_july*(360-xlong(i,j)) 
  • trunk/SOURCES/Greeneem_files/lect-greeneem_mod.f90

    r113 r237  
    8686       do i=1,nx 
    8787 
    88           archimed = Bsoc(i,j)+H(i,j)*ro/row -sealevel 
     88          archimed = Bsoc(i,j)+H(i,j)*ro/row -sealevel_2d(i,j) 
    8989 
    9090          if ((archimed.LT.0.).and.(H(I,J).gt.1.E-3)) then    ! le point flotte 
  • trunk/SOURCES/Greeneem_files/output_greeneem_mod-0.4.f90

    r102 r237  
    141141     do j=1,ny 
    142142        do i=1,nx 
    143            if (mask_cal(i,j,k).and..not.flot(i,j).and.H(i,j).GT.1..and.((sealevel-B(i,j)).LE.0.)) then 
     143           if (mask_cal(i,j,k).and..not.flot(i,j).and.H(i,j).GT.1..and.((sealevel_2d(i,j)-B(i,j)).LE.0.)) then 
    144144              isvolf(kk) = isvolf(kk) + H(i,j) 
    145145           else 
    146               isvolf(kk) = isvolf(kk) + H(i,j)-row/ro*(sealevel-B(i,j)) 
     146              isvolf(kk) = isvolf(kk) + H(i,j)-row/ro*(sealevel_2d(i,j)-B(i,j)) 
    147147           endif 
    148148        enddo 
  • trunk/SOURCES/Greenmint40_files/lect-greenmint_mod.f90

    r4 r237  
    112112   do i=1,nx 
    113113 
    114       archim = Bsoc(i,j)+H(i,j)*ro/row -sealevel 
     114      archim = Bsoc(i,j)+H(i,j)*ro/row -sealevel_2d(i,j) 
    115115 
    116116      if ((ARCHIM.LT.0.).and.(H(I,J).gt.1.E-3)) then    ! le point flotte 
  • trunk/SOURCES/Greenmint40_files/output_green_mod-0.4.f90

    r4 r237  
    8989!        calcul de la hauteur au dessus de la flottaison 
    9090 
    91          if (sealevel-B(i,j).le.0.) then    ! socle au dessus du niveau des mers 
     91         if (sealevel_2d(i,j)-B(i,j).le.0.) then    ! socle au dessus du niveau des mers 
    9292               volf=volf+h(i,j) 
    9393         else 
    94             volf=volf+h(i,j)-row/ro*(sealevel-b(i,j)) 
     94            volf=volf+h(i,j)-row/ro*(sealevel_2d(i,j)-b(i,j)) 
    9595         endif 
    9696 
  • trunk/SOURCES/Heino_files/flottab2-0.5-heino.f90

    r4 r237  
    135135      uys1(i,j)=uybar(i,j) 
    136136 
    137       archim = Bsoc(i,j)+H(i,j)*ro/row -sealevel 
     137      archim = Bsoc(i,j)+H(i,j)*ro/row -sealevel_2d(i,j) 
    138138  
    139139arch: if ((ARCHIM.LT.0.).and.(H(I,J).gt.1.E-3)) then    ! le point flotte 
     
    146146             if (igrdline.eq.1) then   ! en cas de grounding line prescrite 
    147147                flot(i,j)=.false. 
    148                 H(i,j)=(10.+sealevel-Bsoc(i,j))*row/ro  ! pour avoir archim=10 
     148                H(i,j)=(10.+sealevel_2d(i,j)-Bsoc(i,j))*row/ro  ! pour avoir archim=10 
    149149                new_flot_point(i,j)=.false. 
    150150             endif 
     
    162162                if (igrdline.eq.1) then   ! en cas de grounding line prescrite 
    163163                    flot(i,j)=.false. 
    164                     H(i,j)=(10.+sealevel-Bsoc(i,j))*row/ro  ! pour avoir archim=10 
     164                    H(i,j)=(10.+sealevel_2d(i,j)-Bsoc(i,j))*row/ro  ! pour avoir archim=10 
    165165                   new_flot_point(i,j)=.false. 
    166166               endif 
     
    184184 
    185185                surnet=H(i,j)*(1.-ro/row) 
    186                 S(i,j)=surnet+sealevel 
     186                S(i,j)=surnet+sealevel_2d(i,j) 
    187187                B(i,j)=S(i,j)-H(i,j) 
    188188             end if 
     
    234234!                 sur le demi noeud condition archim > 100 m 
    235235 
    236              archim=(Bsoc(i,j)+Bsoc(i-1,j))*0.5-sealevel+ro/row*Hmx(i,j) 
     236             archim=(Bsoc(i,j)+Bsoc(i-1,j))*0.5-sealevel_2d(i,j)+ro/row*Hmx(i,j) 
    237237             gzmx(i,j)=gzmx(i,j).and.(archim.le.100.)  
    238238 
     
    301301!                 sur le demi noeud condition archim > 100 m 
    302302 
    303              archim=(Bsoc(i,j)+Bsoc(i,j-1))*0.5-sealevel+ro/row*Hmy(i,j) 
     303             archim=(Bsoc(i,j)+Bsoc(i,j-1))*0.5-sealevel_2d(i,j)+ro/row*Hmy(i,j) 
    304304             gzmy(i,j)=gzmy(i,j).and.(archim.le.100.)  
    305305 
     
    554554!!!!!!!!!                      on teste si la tache n'est pas completement ice stream 
    555555      if (icetrim(table_out(i,j))) THEN ! on a une ile d'ice stream   !!!! TEST2 
    556               smax_i=0 ; smax_j=0 ; smax_=sealevel 
     556              smax_i=0 ; smax_j=0 ; smax_=sealevel !afq -- warning: what about local sealevel? 
    557557              DO II=3,NX-2 
    558558                DO JJ=3,NY-2 
  • trunk/SOURCES/Hemin40_files/bmelt-hemin40-depth_mod.f90

    r4 r237  
    8585         do i=1,nx 
    8686 
    87             if ((sealevel-bsoc(i,j)).lt.600) then 
     87            if ((sealevel_2d(i,j)-bsoc(i,j)).lt.600) then 
    8888                    bmshelf(i,j)=0.2 ! 0.45!   0.65   
    8989                    bmgrz(i,j)=0.2 ! 0.45 !  0.65 
     
    129129shelf:    if (flot(i,j)) then    ! partie flottante 
    130130 
    131             if ((sealevel-bsoc(i,j)).lt.600) then 
     131            if ((sealevel_2d(i,j)-bsoc(i,j)).lt.600) then 
    132132                    if ((j.lt.75).and.(i.lt.145) ) then  !Atlantique  
    133133                    bmelt(i,j)=(.6+.4*coefbmshelf)*bmshelf(i,j) 
  • trunk/SOURCES/Hudson_files/bmelt_hudson_mod.f90

    r4 r237  
    5252!             bmgrz(i,j)= 10. 
    5353 
    54           bmshelf(i,j)=(sealevel-BSOC(i,j))/200  ! jalv: eq lineaire en fonction de la profondeur de la mer 
     54          bmshelf(i,j)=(sealevel_2d(i,j)-BSOC(i,j))/200  ! jalv: eq lineaire en fonction de la profondeur de la mer 
    5555          bmgrz(i,j)=bmshelf(i,j) 
    5656!          endif 
  • trunk/SOURCES/Hudson_files/climat-hudson_mod.f90

    r4 r237  
    102102!jalv: test pour que le bilan de masse soit 0 sur l'ocean plus profond de 600 m: 
    103103!      c'est le bon critere??: 
    104 where (sealevel-BSOC(:,:).gt.600) acc(:,:)=0. 
    105 where (sealevel-BSOC(:,:).gt.600) abl(:,:)=0. 
     104where (sealevel_2d(:,:)-BSOC(:,:).gt.600) acc(:,:)=0. 
     105where (sealevel_2d(:,:)-BSOC(:,:).gt.600) abl(:,:)=0. 
    106106 
    107 where ((sealevel-BSOC(:,:).gt.100).and.(dist_heino(:,:).gt.2500.)) abl(:,:)=-1. 
    108 where ((sealevel-BSOC(:,:).gt.100).and.(dist_heino(:,:).gt.2600.)) abl(:,:)=-5. 
    109 where ((sealevel-BSOC(:,:).gt.100).and.(dist_heino(:,:).gt.2650.)) abl(:,:)=-10. 
    110 where ((sealevel-BSOC(:,:).gt.100).and.(dist_heino(:,:).gt.2700.)) abl(:,:)=-20. 
     107where ((sealevel_2d(:,:)-BSOC(:,:).gt.100).and.(dist_heino(:,:).gt.2500.)) abl(:,:)=-1. 
     108where ((sealevel_2d(:,:)-BSOC(:,:).gt.100).and.(dist_heino(:,:).gt.2600.)) abl(:,:)=-5. 
     109where ((sealevel_2d(:,:)-BSOC(:,:).gt.100).and.(dist_heino(:,:).gt.2650.)) abl(:,:)=-10. 
     110where ((sealevel_2d(:,:)-BSOC(:,:).gt.100).and.(dist_heino(:,:).gt.2700.)) abl(:,:)=-20. 
    111111 
    112112bm(:,:)=acc(:,:)+abl(:,:) 
     
    135135!      c'est le bon critere??: 
    136136!where (BSOC(:,:).lt.B(:,:)) Tann(:,:)=0. 
    137 where (sealevel-BSOC(:,:).gt.600) Tann(:,:)=0.  
     137where (sealevel_2d(:,:)-BSOC(:,:).gt.600) Tann(:,:)=0.  
    138138 
    139139 
  • trunk/SOURCES/Hudson_files/eaubasale-0.5_hudson_mod.f90

    r4 r237  
    139139        if (flot(i,j))then  ! points flottants 
    140140           klimit(i,j)=1 
    141            limit_hw(i,j)=(sealevel-Bsoc(i,j))*rowg/rofreshg 
     141           limit_hw(i,j)=(sealevel_2d(i,j)-Bsoc(i,j))*rowg/rofreshg 
    142142 
    143143        else if (IBASE(I,J).eq.1) then ! base froide 
     
    182182           pot_w(I,J)=pot_w(I,J)+rog*H(I,J) 
    183183 
    184            pot_f(I,J)=rowg*(sealevel-S(i,j)+H(I,J))  ! pression a la base de l'ice shelf 
     184           pot_f(I,J)=rowg*(sealevel_2d(i,j)-S(i,j)+H(I,J))  ! pression a la base de l'ice shelf 
    185185        enddo 
    186186     enddo 
     
    308308 
    309309              if (flot(i,j)) then  ! if flot > hwater=hwater in ocean  
    310                  hwater(i,j)= sealevel - bsoc(i,j) 
     310                 hwater(i,j)= sealevel_2d(i,j) - bsoc(i,j) 
    311311                 !             hwater(i,j)= max(0.,hwater(i,j)) 
    312312                 if (hwater(i,j).lt.0.) hwater(i,j)=0. 
     
    320320 
    321321              ! Attention le bloc suivant est pour le run 20  
    322               !           Zflot=row/ro*(sealevel-Bsoc(i,j))-10. 
     322              !           Zflot=row/ro*(sealevel_2d(i,j)-Bsoc(i,j))-10. 
    323323              Zflot=H(i,j)*rog/rofreshg 
    324324 
  • trunk/SOURCES/MISMIP3D_files/lect-mismip3d_mod.f90

    r4 r237  
    111111  do j=1,ny 
    112112     do i=1,nx 
    113         if ((bsoc(i,j)+h(i,j)*ro/row -sealevel).lt.0.) then 
     113        if ((bsoc(i,j)+h(i,j)*ro/row -sealevel_2d(i,j)).lt.0.) then 
    114114           flot(i,j)=.true. 
    115115        else 
  • trunk/SOURCES/New-remplimat/remplimat-shelves-tabTu.f90

    r180 r237  
    483483 
    484484         case(-2)                             ! bord Est 
    485             call vel_U_East(pvi(i,j),pvi(i-1,j),rog,H(i-1,j),rowg,slv(i-1,j),B(i-1,j),dx,Tu(i,j,:,:),Tv(i,j,:,:),opposx(i,j)) 
     485            call vel_U_East(pvi(i,j),pvi(i-1,j),rog,H(i-1,j),rowg,sealevel_2d(i-1,j),B(i-1,j),dx,Tu(i,j,:,:),Tv(i,j,:,:),opposx(i,j)) 
    486486        
    487487         case(-3)                             ! bord nord 
     
    489489 
    490490         case(-4)                             ! bord West 
    491             call vel_U_West(pvi(i,j),pvi(i-1,j),rog,H(i,j),rowg,slv(i,j),B(i,j),dx,Tu(i,j,:,:),Tv(i,j,:,:),opposx(i,j)) 
     491            call vel_U_West(pvi(i,j),pvi(i-1,j),rog,H(i,j),rowg,sealevel_2d(i,j),B(i,j),dx,Tu(i,j,:,:),Tv(i,j,:,:),opposx(i,j)) 
    492492 
    493493         case(-12)                            ! coin SE 
     
    527527             
    528528         case(-1)                             ! bord sud 
    529             call vel_V_South(pvi(i,j),pvi(i,j-1),rog,H(i,j),rowg,slv(i,j),B(i,j),dx,Sv(i,j,:,:),Su(i,j,:,:),opposy(i,j)) 
     529            call vel_V_South(pvi(i,j),pvi(i,j-1),rog,H(i,j),rowg,sealevel_2d(i,j),B(i,j),dx,Sv(i,j,:,:),Su(i,j,:,:),opposy(i,j)) 
    530530             
    531531         case(-2)                             ! bord Est 
     
    533533             
    534534         case(-3)                             ! bord nord 
    535             call vel_V_North(pvi(i,j),pvi(i,j-1),rog,H(i-1,j),rowg,slv(i-1,j),B(i-1,j),dx,Sv(i,j,:,:),Su(i,j,:,:),opposy(i,j)) 
     535            call vel_V_North(pvi(i,j),pvi(i,j-1),rog,H(i-1,j),rowg,sealevel_2d(i-1,j),B(i-1,j),dx,Sv(i,j,:,:),Su(i,j,:,:),opposy(i,j)) 
    536536             
    537537         case(-4)                             ! bord West 
  • trunk/SOURCES/Snowball_files/bmelt-snowball-depth_mod.f90

    r57 r237  
    6666         do i=1,nx 
    6767 
    68             if ((sealevel-bsoc(i,j)).lt.600) then 
     68            if ((sealevel_2d(i,j)-bsoc(i,j)).lt.600) then 
    6969                    bmshelf(i,j)=0.2 ! 0.45!   0.65   
    7070                    bmgrz(i,j)=0.2 ! 0.45 !  0.65 
     
    106106shelf:    if (flot(i,j)) then    ! partie flottante 
    107107 
    108             if ((sealevel-bsoc(i,j)).lt.600) then 
     108            if ((sealevel_2d(i,j)-bsoc(i,j)).lt.600) then 
    109109                    if ((j.lt.75).and.(i.lt.145) ) then  !Atlantique  
    110110                    bmelt(i,j)=(.6+.4*coefbmshelf)*bmshelf(i,j) 
  • trunk/SOURCES/Snowball_files/output_snowball_mod-0.4.f90

    r57 r237  
    8383                  vol=vol+h(i,j)  
    8484!        calcul de la hauteur au dessus de la flottaison 
    85                   if (sealevel-B(i,j).le.0.) then    ! socle au dessus du niveau des mers 
     85                  if (sealevel_2d(i,j)-B(i,j).le.0.) then    ! socle au dessus du niveau des mers 
    8686                     volf=volf+h(i,j) 
    8787                  else 
    88                      volf=volf+h(i,j)-row/ro*(sealevel-b(i,j)) 
     88                     volf=volf+h(i,j)-row/ro*(sealevel_2d(i,j)-b(i,j)) 
    8989                  endif 
    9090 
  • trunk/SOURCES/bilan_eau_mod.f90

    r192 r237  
    4545subroutine init_bilan_eau 
    4646 ! initialisation des variables 
    47         diff_H=0. 
    48         sum_H_old = sum(H(2:nx-1,2:ny-1),mask=ice(2:nx-1,2:ny-1)==1) 
    49         tot_water(:,:)=0. 
    50         bm_dt(:,:)=0. 
    51         bmelt_dt(:,:)=0. 
     47        diff_H=0. 
     48        sum_H_old = sum(H(2:nx-1,2:ny-1),mask=ice(2:nx-1,2:ny-1)==1) 
     49        tot_water(:,:)=0. 
     50        bm_dt(:,:)=0. 
     51        bmelt_dt(:,:)=0. 
    5252        alpha_flot=ro/row 
    53         archimtab(:,:) = Bsoc(:,:)+H(:,:)*alpha_flot - sealevel 
     53        archimtab(:,:) = Bsoc(:,:)+H(:,:)*alpha_flot - sealevel_2d(:,:) 
    5454        gr_line(:,:)=0 
    5555        do j=1,ny 
    5656           do i=1,nx 
    57               !afq    if ((H(i,j).gt.0.).and.(archimtab(i,j).GE.0.).and.(Bsoc(i,j).LE.sealevel)) then ! grounded with ice 
     57              !afq    if ((H(i,j).gt.0.).and.(archimtab(i,j).GE.0.).and.(Bsoc(i,j).LE.sealevel_2d(i,j))) then ! grounded with ice 
    5858              if ((H(i,j).gt.0.).and.(archimtab(i,j).GE.0.)) then ! grounded with ice 
    5959                 if (archimtab(i-1,j).LT.0..and.Uxbar(i,j).LT.0..and..not.flot_marais(i-1,j)) gr_line(i,j)=1 
     
    6767end subroutine init_bilan_eau 
    6868 
    69          
     69 
    7070 
    7171 
     
    101101 
    102102where (ice(2:nx-1,2:ny-1).eq.1) 
    103         Bm_dtt(2:nx-1,2:ny-1) = Bm_dtt(2:nx-1,2:ny-1) + Bm_dt(2:nx-1,2:ny-1) !* dt ! somme Bm sur dt 
    104         bmelt_dtt(2:nx-1,2:ny-1) = bmelt_dtt(2:nx-1,2:ny-1) + bmelt_dt(2:nx-1,2:ny-1) ! * dt ! somme bmelt sur dt 
     103        Bm_dtt(2:nx-1,2:ny-1) = Bm_dtt(2:nx-1,2:ny-1) + Bm_dt(2:nx-1,2:ny-1) !* dt ! somme Bm sur dt 
     104        bmelt_dtt(2:nx-1,2:ny-1) = bmelt_dtt(2:nx-1,2:ny-1) + bmelt_dt(2:nx-1,2:ny-1) ! * dt ! somme bmelt sur dt 
    105105endwhere 
    106106 
    107 archimtab(:,:) = Bsoc(:,:)+H(:,:)*alpha_flot - sealevel 
     107archimtab(:,:) = Bsoc(:,:)+H(:,:)*alpha_flot - sealevel_2d(:,:) 
    108108gr_line(:,:)=0 
    109109do j=1,ny 
    110110  do i=1,nx 
    111      !afq    if ((H(i,j).gt.0.).and.(archimtab(i,j).GE.0.).and.(Bsoc(i,j).LE.sealevel)) then ! grounded with ice 
     111     !afq    if ((H(i,j).gt.0.).and.(archimtab(i,j).GE.0.).and.(Bsoc(i,j).LE.sealevel_2d(i,j))) then ! grounded with ice 
    112112     if ((H(i,j).gt.0.).and.(archimtab(i,j).GE.0.)) then ! grounded with ice 
    113113      if (archimtab(i-1,j).LT.0..and.Uxbar(i,j).LT.0..and..not.flot_marais(i-1,j)) gr_line(i,j)=1 
     
    141141    
    142142! bilan d'eau sur la grille : 
    143         water_bilan=sum(tot_water(:,:)) 
    144         diff_H = diff_H/dtt 
     143        water_bilan=sum(tot_water(:,:)) 
     144        diff_H = diff_H/dtt 
    145145 
    146146!999 format(f0.2,1x,e15.8,1x,i10,8(1x,e15.8)) 
    147147!       write(6,999),time,sum_H,count(ice(:,:)==1),diff_H,water_bilan,sum(calv_dtt(:,:))/dtt,sum(ablbord_dtt(:,:))/dtt,sum(bmelt_dtt(:,:),mask=ice(:,:)==1)/dtt,sum(bm(:,:),mask=ice(:,:)==1),sum(Bm_dtt(:,:))/dtt,sum(bmelt_dtt(:,:))/dtt 
    148         diff_H_water_bilan(2:nx-1,2:ny-1)=tot_water(2:nx-1,2:ny-1)-diff_H_2D(2:nx-1,2:ny-1) 
     148        diff_H_water_bilan(2:nx-1,2:ny-1)=tot_water(2:nx-1,2:ny-1)-diff_H_2D(2:nx-1,2:ny-1) 
    149149 
    150150endif 
  • trunk/SOURCES/calving_frange.f90

    r224 r237  
    8989 
    9090    Hcoup(:,:) = min ( max(                                         & 
    91        (-(Bsoc0(:,:)-sealevel) - prof_plateau)/(prof_abysses-prof_plateau)    & 
     91       (-(Bsoc0(:,:)-sealevel_2d(:,:)) - prof_plateau)/(prof_abysses-prof_plateau)    & 
    9292       *(hcoup_abysses-hcoup_plateau)+hcoup_plateau               & 
    9393       , hcoup_plateau), hcoup_abysses ) 
     
    121121 
    122122    Hcoup(:,:) = min ( max(                                         & 
    123        (-(Bsoc(:,:)-sealevel) - prof_plateau)/(prof_abysses-prof_plateau)    & 
     123       (-(Bsoc(:,:)-sealevel_2d(:,:)) - prof_plateau)/(prof_abysses-prof_plateau)    & 
    124124       *(hcoup_abysses-hcoup_plateau)+hcoup_plateau               & 
    125125       , hcoup_plateau), hcoup_abysses ) 
     
    409409                      calv(i,j)=-h(i,j)  
    410410                      !cdc                     H(i,j)=1.  
    411                       !cdc 1m                     H(i,j)=min(1.,max(0.,(sealevel - Bsoc(i,j))*row/ro-0.01)) 
     411                      !cdc 1m                     H(i,j)=min(1.,max(0.,(sealevel_2d(i,j) - Bsoc(i,j))*row/ro-0.01)) 
    412412                      H(i,j)=0. 
    413                       S(i,j)=H(i,j)*(1.-ro/row) + sealevel 
     413                      S(i,j)=H(i,j)*(1.-ro/row) + sealevel_2d(i,j) !afq -- WARNING: est-ce qu'on veut vraiment mettre S a la valeur locale du niveau marin? 
    414414                      B(i,j)=S(i,j) - H(i,j) 
    415415                      ! ATTENTION ne pas mettre ice=0 sinon degradation bilan d'eau (bm et bmelt non comptabilises dans ce cas) 
     
    439439        ice(:,:)=0 
    440440        H(:,:)=0. 
    441         S(:,:)=H(:,:)*(1.-ro/row) + sealevel 
     441        S(:,:)=H(:,:)*(1.-ro/row) + sealevel_2d(:,:) !afq -- WARNING: est-ce qu'on veut vraiment mettre S a la valeur locale du niveau marin? 
    442442        B(:,:)=S(:,:) - H(:,:) 
    443443    endwhere 
  • trunk/SOURCES/calving_frange_abuk.f90

    r204 r237  
    9696 
    9797    Hcoup(:,:) = min ( max(                                         & 
    98        (-(Bsoc0(:,:)-sealevel) - prof_plateau)/(prof_abysses-prof_plateau)    & 
     98       (-(Bsoc0(:,:)-sealevel_2d(:,:)) - prof_plateau)/(prof_abysses-prof_plateau)    & 
    9999       *(hcoup_abysses-hcoup_plateau)+hcoup_plateau               & 
    100100       , hcoup_plateau), hcoup_abysses ) 
     
    128128 
    129129    Hcoup(:,:) = min ( max(                                         & 
    130        (-(Bsoc(:,:)-sealevel) - prof_plateau)/(prof_abysses-prof_plateau)    & 
     130       (-(Bsoc(:,:)-sealevel_2d(:,:)) - prof_plateau)/(prof_abysses-prof_plateau)    & 
    131131       *(hcoup_abysses-hcoup_plateau)+hcoup_plateau               & 
    132132       , hcoup_plateau), hcoup_abysses ) 
     
    416416                      calv(i,j)=-h(i,j)  
    417417                      !cdc                     H(i,j)=1.  
    418                       !cdc 1m                     H(i,j)=min(1.,max(0.,(sealevel - Bsoc(i,j))*row/ro-0.01)) 
     418                      !cdc 1m                     H(i,j)=min(1.,max(0.,(sealevel_2d(i,j) - Bsoc(i,j))*row/ro-0.01)) 
    419419                      H(i,j)=0. 
    420                       S(i,j)=H(i,j)*(1.-ro/row) + sealevel 
     420                      S(i,j)=H(i,j)*(1.-ro/row) + sealevel_2d(i,j) !afq -- WARNING: est-ce qu'on veut vraiment mettre S a la valeur locale du niveau marin? 
    421421                      B(i,j)=S(i,j) - H(i,j) 
    422422                      ! ATTENTION ne pas mettre ice=0 sinon degradation bilan d'eau (bm et bmelt non comptabilises dans ce cas) 
     
    446446!~         ice(:,:)=0 
    447447!~         H(:,:)=0. 
    448 !~         S(:,:)=H(:,:)*(1.-ro/row) + sealevel 
     448!~         S(:,:)=H(:,:)*(1.-ro/row) + sealevel_2d(:,:) 
    449449!~         B(:,:)=S(:,:) - H(:,:) 
    450450!~     endwhere 
     
    455455        ice(:,:)=0 
    456456        H(:,:)=0. 
    457         S(:,:)=H(:,:)*(1.-ro/row) + sealevel 
     457        S(:,:)=H(:,:)*(1.-ro/row) + sealevel_2d(:,:) !afq -- WARNING: est-ce qu'on veut vraiment mettre S a la valeur locale du niveau marin? 
    458458        B(:,:)=S(:,:) - H(:,:) 
    459459    endwhere     
  • trunk/SOURCES/climat-forcage-insolation_mod.f90

    r65 r237  
    1212! C. Dumas 07/2015 
    1313 
    14 use module3D_phy,only:nx,ny,S,H,flot,slv,Tann,Tjuly,Tmois,acc,coefbmshelf,ro,num_param,num_rep_42,dirnameinp,time,pi,dx,xlong,ylat,bsoc          
     14use module3D_phy,only:nx,ny,S,H,flot,Tann,Tjuly,Tmois,acc,coefbmshelf,sealevel_2d,ro,num_param,num_rep_42,dirnameinp,time,pi,dx,xlong,ylat,bsoc          
    1515use netcdf 
    1616use io_netcdf_grisli 
     
    637637do j=1,ny 
    638638   do i=1,nx 
    639       Zs(i,j)=max(slv(i,j),S(i,j)) 
     639      Zs(i,j)=max(sealevel_2d(i,j),S(i,j)) 
    640640      do ictr=CO2SnapMin,CO2SnapMax,1 
    641641         do igtr=1,gtr 
  • trunk/SOURCES/climat-forcage-insolation_mod_oneway.f90

    r108 r237  
    1212! C. Dumas 06/2015 
    1313 
    14 USE module3D_phy,only:nx,ny,S,slv,Tann,Tjuly,Tmois,acc,coefbmshelf,ro,num_param,num_rep_42,dirnameinp,time 
     14USE module3D_phy,only:nx,ny,S,sealevel_2d,Tann,Tjuly,Tmois,acc,coefbmshelf,ro,num_param,num_rep_42,dirnameinp,time 
    1515!use interface_input 
    1616use netcdf 
     
    245245do j=1,ny 
    246246   do i=1,nx 
    247       Zs(i,j)=max(slv(i,j),S(i,j)) 
     247      Zs(i,j)=max(sealevel_2d(i,j),S(i,j)) 
    248248      !Il faut mettre S0(i,j) si pas d'interpolation sinon Ssnap(i,j,ictr,igtr) 
    249249      !if (Zs(i,j).ge.Ssnap(i,j,ictr,igtr)) then 
  • trunk/SOURCES/climat-perturb_mod-0.4.f90

    r93 r237  
    1313 
    1414 
    15 use module3d_phy,only:nx,ny,S,S0,Tann,Tjuly,precip,acc,Ylat,num_forc,num_param,num_rep_42,dirforcage,dirnameinp,tafor,time,sealevel,coefbmshelf 
     15use module3d_phy,only:nx,ny,S,S0,Tann,Tjuly,precip,acc,Ylat,num_forc,num_param,num_rep_42,dirforcage,dirnameinp,tafor,time,sealevel,sealevel_2d,coefbmshelf 
    1616use netcdf 
    1717use io_netcdf_grisli 
     
    307307100 continue 
    308308 
     309  sealevel_2d(:,:)=sealevel 
    309310 
    310311  !  coefmshelf est un coefficient qui fait varier bmgrz et bmshelf 
  • trunk/SOURCES/climat_forcage_mod.f90

    r224 r237  
    66! possibilite d'utiliser autant de snapshots que l'on veut (ntr) 
    77! fichiers de forcage, Snapshot et valeurs de lapse rate definis dans fichier param 
    8   use module3d_phy,only:nx,ny,sealevel,slv,S,S0,Tmois,Tann,Tjuly,acc,num_forc,num_param,num_rep_42,coefbmshelf,PI,ro,time,dirforcage,dirnameinp 
     8  use module3d_phy,only:nx,ny,sealevel,sealevel_2d,S,S0,Tmois,Tann,Tjuly,acc,num_forc,num_param,num_rep_42,coefbmshelf,PI,ro,time,dirforcage,dirnameinp 
    99  use netcdf 
    1010  use io_netcdf_grisli 
     
    279279100 continue 
    280280 
    281   slv(:,:)=sealevel 
     281  sealevel_2d(:,:)=sealevel 
    282282 
    283283  if(timeloc.le.ttr(1)) then   ! time en dehors des limites du fichier forcage 
     
    333333!     calcul de Tjuly et et Tann 
    334334  !$OMP WORKSHARE 
    335   Zs(:,:)=max(slv(:,:),S(:,:)) 
     335  Zs(:,:)=max(sealevel_2d(:,:),S(:,:)) 
    336336  deltaZs(:,:)= -lapserate*(Zs(:,:)-S0(:,:)) ! difference de temperature entre Zs et S0 
    337337  !$OMP END WORKSHARE 
  • trunk/SOURCES/eaubasale-0.5_mod.f90

    r196 r237  
    158158        if (flot(i,j))then  ! points flottants 
    159159           klimit(i,j)=1 
    160            limit_hw(i,j)=(sealevel-Bsoc(i,j))*rowg/rofreshg 
     160           limit_hw(i,j)=(sealevel_2d(i,j)-Bsoc(i,j))*rowg/rofreshg 
    161161 
    162162        else if (IBASE(I,J).eq.1) then ! base froide 
     
    200200           pot_w(I,J)=pot_w(I,J)+rog*H(I,J) 
    201201 
    202            pot_f(I,J)=rowg*(sealevel-S(i,j)+H(I,J))  ! pression a la base de l'ice shelf 
     202           pot_f(I,J)=rowg*(sealevel_2d(i,j)-S(i,j)+H(I,J))  ! pression a la base de l'ice shelf 
    203203        enddo 
    204204     enddo 
     
    314314 
    315315              if (flot(i,j)) then  ! if flot > hwater=hwater in ocean  
    316                  hwater(i,j)= sealevel - bsoc(i,j) 
     316                 hwater(i,j)= sealevel_2d(i,j) - bsoc(i,j) 
    317317                 !             hwater(i,j)= max(0.,hwater(i,j)) 
    318318                 if (hwater(i,j).lt.0.) hwater(i,j)=0. 
  • trunk/SOURCES/flottab2-0.7.f90

    r223 r237  
    122122       where ( mk_init(:,:).eq.1)                 ! pose 
    123123          flot(:,:) = .False. 
    124           H(:,:)=max(H(:,:),(10.+sealevel-Bsoc(:,:))*row/ro)  ! pour avoir archim=10 
     124          H(:,:)=max(H(:,:),(10.+sealevel_2d(:,:)-Bsoc(:,:))*row/ro)  ! pour avoir archim=10 
    125125          S(:,:) = Bsoc(:,:) + H(:,:)  
    126126       elsewhere 
     
    175175          uys1(i,j)=uybar(i,j) 
    176176 
    177           archim = Bsoc(i,j)+H(i,j)*ro/row -sealevel 
     177          archim = Bsoc(i,j)+H(i,j)*ro/row -sealevel_2d(i,j) 
    178178          !      if ((i.eq.132).and.(j.eq.183)) print*,'archim=',archim 
    179179 
     
    189189                if (igrdline.eq.1) then   ! en cas de grounding line prescrite 
    190190                   flot(i,j)=.false. 
    191                    H(i,j)=(10.+sealevel-Bsoc(i,j))*row/ro  ! pour avoir archim=10 
     191                   H(i,j)=(10.+sealevel_2d(i,j)-Bsoc(i,j))*row/ro  ! pour avoir archim=10 
    192192                   new_flot_point(i,j)=.false. 
    193193                endif 
     
    202202                   if (igrdline.eq.1) then   ! en cas de grounding line prescrite 
    203203                      flot(i,j)=.false. 
    204                       H(i,j)=(10.+sealevel-Bsoc(i,j))*row/ro  ! pour avoir archim=10 
     204                      H(i,j)=(10.+sealevel_2d(i,j)-Bsoc(i,j))*row/ro  ! pour avoir archim=10 
    205205                      new_flot_point(i,j)=.false. 
    206206                   endif 
     
    228228             flot(i,j)=.true.  !cdc points ocean sont flot meme sans glace 
    229229             H(i,j)=0. 
    230              S(i,j)=sealevel 
     230             S(i,j)=sealevel_2d(i,j) !afq -- WARNING: est-ce qu'on veut vraiment mettre S a la valeur locale du niveau marin? 
    231231             B(i,j)=S(i,j)-H(i,j) 
    232232          endif arch 
     
    239239 
    240240             surnet=H(i,j)*(1.-ro/row) 
    241              S(i,j)=surnet+sealevel 
     241             S(i,j)=surnet+sealevel_2d(i,j) !afq -- WARNING: est-ce qu'on veut vraiment mettre S a la valeur locale du niveau marin? 
    242242             B(i,j)=S(i,j)-H(i,j) 
    243243          end if 
     
    246246    end do 
    247247    !$OMP END DO 
    248     !~ print*,'flottab 2',S(132,183),H(132,183),BSOC(132,183),B(132,183),sealevel 
    249     !~ print*,'flottab 2',flot(132,183),ice(132,183) 
    250     !~ if (S(132,183).LT.sealevel) print*,'BUGGGGGGGGGGGGG !!!!!!!!!!!!!!!!' 
    251248        
    252249!!$ do i=1,nx 
     
    262259 
    263260 
    264  
    265261    !----------------------------------------------------------------------- 
    266262    !$OMP DO 
     
    298294          !                 sur le demi noeud condition archim < 100 m 
    299295 
    300           archim=(Bsoc(i,j)+Bsoc(i-1,j))*0.5-sealevel+ro/row*Hmx(i,j) 
     296          archim=(Bsoc(i,j)+Bsoc(i-1,j))*0.5-(sealevel_2d(i,j)+sealevel_2d(i-1,j))*0.5+ro/row*Hmx(i,j) 
    301297          gzmx(i,j)=gzmx(i,j).and.(archim.le.100.)  
    302298          cotemx(i,j)=gzmx(i,j) 
     
    337333          !                 sur le demi noeud condition archim > 100 m 
    338334 
    339           archim=(Bsoc(i,j)+Bsoc(i,j-1))*0.5-sealevel+ro/row*Hmy(i,j) 
     335          archim=(Bsoc(i,j)+Bsoc(i,j-1))*0.5-(sealevel_2d(i,j)+sealevel_2d(i,j-1))*0.5+ro/row*Hmy(i,j) 
    340336          gzmy(i,j)=gzmy(i,j).and.(archim.le.100.)  
    341337          cotemy(i,j)=gzmy(i,j) 
     
    605601          ice(:,:)=0 
    606602          H(:,:)=0. 
    607           S(:,:)=H(:,:)*(1.-ro/row) + sealevel 
     603          S(:,:)=H(:,:)*(1.-ro/row) + sealevel_2d(:,:) 
    608604          B(:,:)=S(:,:) - H(:,:) 
    609605       end where 
     
    790786            flgzmy(smax_i,smax_j)=.false. ; flgzmx(smax_i,smax_j+1)=.false. 
    791787      
    792             if (Smax_.le.sealevel) then 
     788            if (Smax_.le.sealevel) then  ! afq -- WARNING: en fait avec le niveau marin local cest plus complique que ca! 
     789            !if (Smax_.le.sealevel_2d(i,j)) then ! afq --WARNING: il faudrait regarder point par point sur la tache... 
    793790              write(num_tracebug,*)'Attention, une ile avec la surface sous l eau' 
    794791              write(num_tracebug,*)'time=',time,'   coord:',smax_i,smax_j 
  • trunk/SOURCES/furst_schoof_mod.f90

    r228 r237  
    8585  gr_line_schoof(:,:)=0 
    8686 
    87   archimtab(:,:) = Bsoc(:,:)+H(:,:)*alpha_flot -sealevel 
     87  archimtab(:,:) = Bsoc(:,:)+H(:,:)*alpha_flot -sealevel_2d(:,:) 
    8888 
    8989  imx_diag_in(:,:) = imx_diag(:,:) 
     
    9393     do i=3,nx-3 
    9494        !archim = Bsoc(i,j)+H(i,j)*ro/row -sealevel 
    95         if ((H(i,j).gt.0.).and.(archimtab(i,j).GE.0.).and.(Bsoc(i,j).LE.sealevel)) then    ! grounded and with ice 
     95        if ((H(i,j).gt.0.).and.(archimtab(i,j).GE.0.).and.(Bsoc(i,j).LE.sealevel_2d(i,j))) then    ! grounded and with ice 
    9696           !           print*,'schoof gr line ij',i,j 
    9797           ! selon x (indice i) 
     
    100100           !if (flot(i-1,j).and.Uxbar(i,j).LT.0.) then    ! grounding line between i-1,i    CAS WEST (ecoulement vers grille West) 
    101101 
    102               call  calc_xgl4schoof(-dx,alpha_flot,H(i,j),Bsoc(i,j),H(i-1,j),Bsoc(i-1,j),Sealevel,xpos,Hgl) 
     102              call  calc_xgl4schoof(-dx,alpha_flot,H(i,j),Bsoc(i,j),sealevel_2d(i,j),H(i-1,j),Bsoc(i-1,j),sealevel_2d(i-1,j),xpos,Hgl) 
    103103              xpos_tab(i,j)=xpos 
    104104              Hglx_tab(i,j)=Hgl 
     
    194194           !else if (flot(i+1,j).and.Uxbar(i+1,j).GT.0.) then    ! grounding line between i,i+1    CAS EST (ecoulement vers grille Est) 
    195195 
    196               call  calc_xgl4schoof(dx,alpha_flot,H(i,j),Bsoc(i,j),H(i+1,j),Bsoc(i+1,j),Sealevel,xpos,Hgl) 
     196              call  calc_xgl4schoof(dx,alpha_flot,H(i,j),Bsoc(i,j),sealevel_2d(i,j),H(i+1,j),Bsoc(i+1,j),sealevel_2d(i+1,j),xpos,Hgl) 
    197197              xpos_tab(i,j)=xpos 
    198198              Hglx_tab(i,j)=Hgl 
     
    285285 
    286286 
    287               call  calc_xgl4schoof(-dy,alpha_flot,H(i,j),Bsoc(i,j),H(i,j-1),Bsoc(i,j-1),Sealevel,ypos,Hgl) 
     287              call  calc_xgl4schoof(-dy,alpha_flot,H(i,j),Bsoc(i,j),sealevel_2d(i,j),H(i,j-1),Bsoc(i,j-1),sealevel_2d(i,j-1),ypos,Hgl) 
    288288              ypos_tab(i,j)=ypos 
    289289              Hgly_tab(i,j)=Hgl 
     
    373373           !else if (flot(i,j+1).and.Uybar(i,j+1).GT.0.) then    ! grounding line between j,j+1    CAS NORD (ecoulement vers grille Nord) 
    374374 
    375               call  calc_xgl4schoof(dy,alpha_flot,H(i,j),Bsoc(i,j),H(i,j+1),Bsoc(i,j+1),Sealevel,ypos,Hgl) 
     375              call  calc_xgl4schoof(dy,alpha_flot,H(i,j),Bsoc(i,j),sealevel_2d(i,j),H(i,j+1),Bsoc(i,j+1),sealevel_2d(i,j+1),ypos,Hgl) 
    376376              ypos_tab(i,j)=ypos 
    377377              Hgly_tab(i,j)=Hgl 
     
    507507!------------------------------------------------------------------------- 
    508508! calcule la position sous maille de la ligne dechouage 
    509  subroutine calc_xgl4schoof(dy,alpha,H_0,B_0,H_1,B_1,SL,xpos,Hgl) 
     509 subroutine calc_xgl4schoof(dy,alpha,H_0,B_0,SL_0,H_1,B_1,SL_1,xpos,Hgl) 
    510510                 
    511511 
     
    516516  real,intent(in)        ::    H_0       !< epaisseur au point pose 
    517517  real,intent(in)        ::    B_0       !< altitude socle au point pose 
     518  real,intent(in)        ::    SL_0      !< sea level au point pose 
    518519  real,intent(in)        ::    H_1       !< epaisseur au point flottant 
    519520  real,intent(in)        ::    B_1       !< altitude socle au point flottant 
    520   real,intent(in)        ::    SL        !< Sea level 
     521  real,intent(in)        ::    SL_1      !< sea level au point flottant 
    521522  real,intent(out)       ::    xpos      !< position de la ligne (en distance depuis le point pose) 
    522   real,intent(out)       ::    Hgl       !< Epaisseur de glace a la grounding line 
     523  real,intent(out)       ::    Hgl       !< epaisseur de glace a la grounding line 
    523524    
    524525   
     
    541542!~      endif 
    542543                 
    543   Cgl   = (B_1-B_0) + alpha * (H_1-H_0) 
     544  Cgl   = (B_1-B_0) + alpha * (H_1-H_0) - (SL_1-SL_0) 
    544545 
    545546  if (abs(Cgl).gt.1.e-5) then    
    546      xpos    = (SL - (B_0 + alpha * H_0)) / Cgl  
     547     xpos    = (SL_0 - (B_0 + alpha * H_0)) / Cgl  
    547548  else 
    548549     xpos    = 0.5  ! verifier 
     
    558559     print*,'calc_xgl4schoof : xpos value error i,j=',i,j 
    559560     print*, 'xpos,Cgl', xpos,Cgl 
    560      print*,'B_0, alpha, H_0, SL', B_0, alpha, H_0, SL 
    561      print*,'archim=',B_1+H_1*alpha - SL 
     561     print*,'B_0, alpha, H_0, SL_0', B_0, alpha, H_0, SL_0 
     562     print*,'archim=',B_1+H_1*alpha - SL_1 
    562563     !stop 
    563564     xpos = min(max(0.,xpos),1.) 
  • trunk/SOURCES/initial-phy-2-job.f90

    r10 r237  
    192192  TJFOR=0.0 
    193193  SEALEVEL=0.0 
     194  sealevel_2d(:,:)=0. 
    194195 
    195196  SECYEAR=365.*24.*3600. 
  • trunk/SOURCES/initial-phy-2.f90

    r233 r237  
    185185  TJFOR=0.0 
    186186  SEALEVEL=0.0 
     187  sealevel_2d(:,:)=0. 
    187188 
    188189  SECYEAR=365.*24.*3600. 
  • trunk/SOURCES/initial2-0.4.f90

    r4 r237  
    7171  do J=1,NY 
    7272     do I=1,NX 
    73         if ((BSOC(I,J)+H(I,J)*RO/ROW -SEALEVEL).LT.0.) then 
     73        if ((BSOC(I,J)+H(I,J)*RO/ROW -SEALEVEL_2D(I,J)).LT.0.) then 
    7474           FLOT(I,J)=.TRUE. 
    7575        else 
  • trunk/SOURCES/isostasie_mod-0.3.f90

    r181 r237  
    103103    do i=1,nx 
    104104       do j=1,ny 
    105           if (ro*h0(i,j).ge.-row*(BSOC0(i,j)-sealevel)) then 
     105          if (ro*h0(i,j).ge.-row*(BSOC0(i,j)-sealevel_2d(i,j))) then 
    106106 
    107107             !        glace ou terre 
     
    110110          else 
    111111             !        ocean 
    112              charge(i,j)=-rowg*(Bsoc0(i,j)-sealevel) 
     112             charge(i,j)=-rowg*(Bsoc0(i,j)-sealevel_2d(i,j)) 
    113113          endif 
    114114       end do 
  • trunk/SOURCES/no_lakes.f90

    r4 r237  
    3131subroutine lake_to_slv 
    3232 
    33 slv(:,:)=sealevel 
     33sealevel_2d(:,:)=sealevel 
    3434 
    3535end subroutine lake_to_slv 
  • trunk/SOURCES/prescribe-H-i2s_mod.f90

    r124 r237  
    198198 
    199199    where (MK_gl0(:,:) .eq. 1) 
    200        hp(:,:) = (sealevel-Bsoc(:,:))*row/ro +1. 
     200       hp(:,:) = (sealevel_2d(:,:)-Bsoc(:,:))*row/ro +1. 
    201201    end where 
    202202 
  • trunk/SOURCES/steps_time_loop.f90

    r182 r237  
    6868     where(.not.flot(:,:)) 
    6969        S(:,:)=Bsoc(:,:)+H(:,:) 
    70         S(:,:)=max(S(:,:),sealevel) 
     70        S(:,:)=max(S(:,:),sealevel_2d(:,:)) !afq -- WARNING: est-ce qu'on veut vraiment mettre S a la valeur locale du niveau marin? 
    7171        B(:,:)=Bsoc(:,:) 
    7272     end where 
  • trunk/SOURCES/steps_time_loop_avec_iterbeta.f90

    r180 r237  
    6969     where(.not.flot(:,:)) 
    7070        S(:,:)=Bsoc(:,:)+H(:,:) 
    71         S(:,:)=max(S(:,:),sealevel) 
     71        S(:,:)=max(S(:,:),sealevel_2d(:,:)) !afq -- WARNING: est-ce qu'on veut vraiment mettre S a la valeur locale du niveau marin? 
    7272        B(:,:)=Bsoc(:,:) 
    7373     end where 
  • trunk/SOURCES/taubed-0.3.f90

    r181 r237  
    6060     do J=1,NY  
    6161        do I=1,NX 
    62            if (RO*H(I,J).ge.-ROW*(BSOC(I,J)-SEALEVEL)) then 
     62           if (RO*H(I,J).ge.-ROW*(BSOC(I,J)-SEALEVEL_2D(i,j))) then 
    6363              !           glace ou terre  
    6464              CHARGE(I,J)=ROG*H(I,J) 
    6565           else 
    6666              !           ocean 
    67               CHARGE(I,J)=-ROWG*(BSOC(I,J)-SEALEVEL) 
     67              CHARGE(I,J)=-ROWG*(BSOC(I,J)-SEALEVEL_2D(i,j)) 
    6868           endif 
    6969        end do 
     
    102102     do J=1,NY 
    103103        do I=1,NX 
    104            if (RO*H(I,J).ge.-ROW*(BSOC(I,J)-SEALEVEL)) then 
     104           if (RO*H(I,J).ge.-ROW*(BSOC(I,J)-SEALEVEL_2D(i,j))) then 
    105105              !           glace ou terre  
    106106              W1(I,J)=RO/ROM*H(I,J) 
    107107           else 
    108108              !           ocean 
    109               W1(I,J)=-ROW/ROM*(BSOC(I,J)-SEALEVEL) 
     109              W1(I,J)=-ROW/ROM*(BSOC(I,J)-SEALEVEL_2D(i,j)) 
    110110           endif 
    111111        end do 
Note: See TracChangeset for help on using the changeset viewer.