Changeset 243

Timestamp:

01/16/19 12:21:03 (5 years ago)

Author:

aquiquet

Message:

Possibility to modulate bmelt following an index in climat_forcage + addition of a forcing file with a glacial index constructed from NGRIP d18O and an oceanic perturbation (sealevel and bmelt) as in the GMD2018 paper

Location:

trunk

Files:

• INPUT/Forcage/forcage-429k-quiquet-et-al_CP2013_glacial-index_Waelbroeck_AMOC.dat (added)
• SOURCES/climat_forcage_mod.f90 (modified) (12 diffs)

trunk/SOURCES/climat_forcage_mod.f90

@@ -6 +6 @@
 ! possibilite d'utiliser autant de snapshots que l'on veut (ntr)
 ! fichiers de forcage, Snapshot et valeurs de lapse rate definis dans fichier param
-  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
+  use module3d_phy,only:nx,ny,sealevel,sealevel_2d,S,S0,Tmois,Tann,Tjuly,acc,num_forc,num_param,num_rep_42,tafor,coefbmshelf,PI,ro,time,dirforcage,dirnameinp
   use netcdf
   use io_netcdf_grisli
@@ -12 +12 @@
 
   implicit none
-  character(len=80) :: filin
+  character(len=100) :: filin
   integer :: nft             !   nft est le nombre de lignes a lire dans le fichier 
                            ! contenant le forcage climatique
   real,dimension(:),allocatable :: tdate          ! time for climate forcing
-  real,dimension(:),allocatable :: alphat
-  real,dimension(:),allocatable :: alphap
+  real,dimension(:),allocatable :: tpert
   real,dimension(:),allocatable :: spert
-
-
-  real :: mincoefbmelt               ! butoirs mini 
-  real :: maxcoefbmelt               ! butoirs maxi de coefbmshelf
+  real,dimension(:),allocatable :: bpert
+
   character(len=100) :: clim_ref_file ! climat de reference
   character(len=100) :: forcage_file1 ! fichier de forcage 1 (climat+topo)
   character(len=100) :: forcage_file2 ! fichier de forcage 2 (climat+topo)
 
-  character(len=80) :: filforc       ! nom du fichier forcage
+  character(len=100) :: filforc       ! nom du fichier forcage
+  real :: coefT                    !< pour modifier l'amplitude de la perturb. T
+  integer :: pertbmb               !< boolean: do we modify the bmelt under ice shelves?
+  real    :: coefbmb               !< if we modify the bmelt, we read a perturbation (based on dD, insolation, etc.) and apply this coefficient
 
 
@@ -153 +153 @@
   end if
 
-  if (.not.allocated(alphat)) THEN
-     allocate(alphat(nft),stat=err)
+  if (.not.allocated(tpert)) THEN
+     allocate(tpert(nft),stat=err)
      if (err/=0) then
         print *,"erreur à l'allocation du tableau tpert",err
@@ -161 +161 @@
   end if
 
-  if (.not.allocated(alphap)) THEN
-     allocate(alphap(nft),stat=err)
+  if (.not.allocated(bpert)) THEN
+     allocate(bpert(nft),stat=err)
      if (err/=0) then
         print *,"erreur à l'allocation du tableau tpert",err
@@ -178 +178 @@
 
   do i=1,nft
-     read(num_forc,*) tdate(i),alphat(i),alphap(i),spert(i)
+     if (pertbmb==1) then
+        read(num_forc,*) tdate(i),spert(i),tpert(i),bpert(i)
+     else
+        read(num_forc,*) tdate(i),spert(i),tpert(i)
+        bpert(i)=0.
+     endif
   end do
 
   close(num_forc)
+
+! Warning: when using climat_forcage we should in principle use a tpert which is a glacial index (ranging from about 0 to 1)
+!          here we can still use a coefT to scale the signal if needed
+!          the coefbmb now is copy and paste from what we have done for Antarctica in Quiquet et al. GMD 2018
+
+  tpert(:)=tpert(:)*coefT
+  bpert(:)=max(1.+bpert(:)*coefbmb,0.01) !freezing is not allowed
+
 
 
@@ -191 +204 @@
 
 
-  namelist/clim_forcage/clim_ref_file,ttr,forcage_file1,forcage_file2,typerun,rappact,mincoefbmelt,maxcoefbmelt,filforc,lapserate
+  namelist/clim_forcage/clim_ref_file,ttr,forcage_file1,forcage_file2,typerun,lapserate,rappact,filforc,pertbmb,coeft,coefbmb
 
   rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
@@ -204 +217 @@
   write(num_rep_42,'(A,A,A)') 'forcage_file2 = "',trim(forcage_file2),'"'
   write(num_rep_42,*) 'typerun = ', typerun
+  write(num_rep_42,*) 'lapserate = ', lapserate
   write(num_rep_42,*) 'rappact = ', rappact
-  write(num_rep_42,*) 'mincoefbmelt = ', mincoefbmelt
-  write(num_rep_42,*) 'maxcoefbmelt = ', maxcoefbmelt
   write(num_rep_42,'(A,A,A)') ' filforc      = "',trim(filforc),'"'
-  write(num_rep_42,*) 'lapserate = ', lapserate
+  write(num_rep_42,*) 'pertbmb = ', pertbmb
+  write(num_rep_42,*) 'coefT   = ', coefT
+  write(num_rep_42,*) 'coefbmb = ', coefbmb
   write(num_rep_42,*)'/'                            
   write(num_rep_42,*)
@@ -225 +239 @@
 
   implicit none
-  real :: coeft,coefp,coeftime       !< coeftime is not used
+  real :: coeftime       !< coeftime is not used
+  real :: tafor
   integer :: l                       !< dumm index for loops on snapsots files  l=itr,ntr-1
   integer :: itr                     !< nb of the current snapshot file (change with time)
@@ -247 +262 @@
   if(timeloc.le.tdate(1)) then   ! time avant le debut du fichier forcage
      sealevel=spert(1)
-     coeft=alphat(1)
-     coefp=alphap(1)
+     tafor=tpert(1)
+     coefbmshelf=bpert(1)
      ift=1
 
   else if (timeloc.ge.tdate(nft)) then     ! time apres la fin du fichier forcage
      sealevel=spert(nft)
-     coeft=alphat(nft)
-     coefp=alphap(nft)
+     tafor=tpert(nft)
+     coefbmshelf=bpert(nft)
      ift=nft
 
@@ -266 +281 @@
                 (timeloc-tdate(i))/(tdate(i+1)-tdate(i))
 
-           coeft=alphat(i)+(alphat(i+1)-alphat(i))*      &
+           tafor=tpert(i)+(tpert(i+1)-tpert(i))*      &
                 (timeloc-tdate(i))/(tdate(i+1)-tdate(i))
 
-           coefp=alphap(i)+(alphap(i+1)-alphap(i))*      &
+           coefbmshelf=bpert(i)+(bpert(i+1)-bpert(i))*      &
                 (timeloc-tdate(i))/(tdate(i+1)-tdate(i))
                 
@@ -279 +294 @@
 100 continue
 
+
+!same as for the GMD2018 paper:
+!afq -- as in Pollard and Deconto, melt is more efficient at high temperature... ->
+  if (coefbmshelf.gt.1.) then
+     coefbmshelf = 1. + ( coefbmshelf - 1. ) * 1.5 !afq.... tuning for the degla?...
+  end if
+       
   sealevel_2d(:,:)=sealevel
 
@@ -322 +344 @@
 
 
-! test avec coefbmshelf fonction de coeft (1 a l'actuel et proche de 0 en glaciaire)
-  coefbmshelf=coeft
-!coefbmshelf=1. ! modif pour test de sensibilite (tof 20 avril 01)
-  coefbmshelf=max(coefbmshelf,mincoefbmelt)
-  coefbmshelf=min(coefbmshelf,maxcoefbmelt)
-
 !if (geoplace(1:5).eq.'hemin') then