Changeset 366 for trunk/SOURCES/climat_forcage_mod.f90

Timestamp:

02/22/23 16:28:12 (15 months ago)

Author:

dumas

Message:

ablation_mod.f90 : SMB can be computed using ITM with pdd_type=3

File:

• trunk/SOURCES/climat_forcage_mod.f90 (modified) (13 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,Tmois,Tann,Tjuly,acc,num_forc,num_param,num_rep_42,tafor,coefbmshelf,PI,ro,time,dirforcage,dirnameinp
+  use module3d_phy,only:nx,ny,sealevel,sealevel_2d,S,Tmois,Tann,Tjuly,acc,swmois,prmois,num_forc,num_param,num_rep_42,tafor,coefbmshelf,PI,ro,time,dirforcage,dirnameinp
   use netcdf
   use io_netcdf_grisli
@@ -23 +23 @@
 
   character(len=100) :: clim_ref_file ! climat de reference
-  character(len=70), dimension(5) :: forcage_file ! fichier de forcage (climat+topo) => dimension >= ntr
+  character(len=100), dimension(12) :: forcage_file ! fichier de forcage (climat+topo) => dimension >= ntr
 !  character(len=100) :: forcage_file2 ! fichier de forcage 2 (climat+topo)
 
@@ -32 +32 @@
   real    :: r_atmvar              !< a ratio to account fast atmos variability between the 2 snapshots (in %: 0=>no fast var,1=>only fast var)
 
-
-
 !  character(len=200) :: filtr(ntr)   ! fichiers de forcage
 
@@ -43 +41 @@
                      ! 2 for steady state using snap 2 ...
                      ! 0 for transient
+
+  integer :: smbcomp ! 0 for PDD
+                     ! 1 for ITM
   
 ! variables pour climato mensuelle :  
   real,dimension(nx,ny,12) :: t2m0   !< initial monthly air temperature (climato)
+  real,dimension(nx,ny,12) :: sw0   !< initial monthly downward SW radiations at the surface (climato)
   real,dimension(nx,ny,12) :: pr0    !< initial monthly precip (climato)
   real,dimension(nx,ny) :: orog0     !< topo de la climato
@@ -51 +53 @@
 ! climat des snapshots :
   real,dimension(:,:,:,:), allocatable :: t2m_ss  !< t2m monthly for every GCM snapshot 
+  real,dimension(:,:,:,:), allocatable :: sw_ss  !< sw radiations monthly for every GCM snapshot
   real,dimension(:,:,:,:), allocatable :: pr_ss   !< precip monthly for every GCM snapshot
   real,dimension(:,:,:), allocatable :: orog_ss    !< topo for every GCM snapshot 
   
   real,dimension(:,:,:,:), allocatable :: t2m_sstopo0 !< t2m GCM snapshot on orog0
+  real,dimension(:,:,:,:), allocatable :: sw_sstopo0  !< sw GCM snapshot on orog0
   real,dimension(:,:,:,:), allocatable :: pr_sstopo0  !< pr GCM snapshot on orog0
-
 
 !~   real,dimension(nx,ny,12,ntr) :: t2m_ss  !< t2m monthly for every GCM snapshot 
@@ -141 +144 @@
   enddo
 
+  if (smbcomp .EQ. 1) then
+     ! lecture du climat de reference : climato mensuelle
+     call Read_Ncdf_var('rsds',trim(DIRNAMEINP)//TRIM(clim_ref_file),tab3d)
+     sw0(:,:,:)=tab3d(:,:,:)
+
+     ! lecture des fichiers snapshots pour n'importe quel geoplace
+     do k=1,ntr
+        call Read_Ncdf_var('rsds',trim(DIRNAMEINP)//'Snapshots-GCM/'//trim(forcage_file(k)),tab3d)
+        sw_ss(:,:,:,k)=tab3d(:,:,:)
+     end do
+
+     ! calcul de sw sur la topo des obs :
+     do m=1,12
+        do k=1,ntr
+           !sw_sstopo0(:,:,m,k)=sw_ss(:,:,m,k)
+           sw_sstopo0(:,:,m,k)=sw_ss(:,:,m,k)+sw_ss(:,:,m,k)*0.00006*(orog0(:,:)-orog_ss(:,:,k)) ! Adapted from Robinson et al., 2011
+        enddo
+     enddo
+
+  end if
 
 ! lecture des fichiers d'evolution pour tout geoplace
@@ -214 +237 @@
 subroutine init_forclim
 
-  real,dimension(5) :: ttr_temp         !< date des tranches (snapshots) 
+  real,dimension(12) :: ttr_temp         !< date des tranches (snapshots) 
   integer :: err                  !< recuperation erreur
   integer :: k
   
-  namelist/clim_forcage/clim_ref_file,ntr,ttr_temp,forcage_file,typerun,lapserate,rappact,filforc,pertbmb,coeft,coefbmb,r_atmvar
+  namelist/clim_forcage/clim_ref_file,ntr,ttr_temp,forcage_file,typerun,lapserate,rappact,filforc,pertbmb,coeft,coefbmb,r_atmvar,smbcomp
 
   rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
@@ -255 +278 @@
    end if
   end if
+
+  if (.not.allocated(sw_ss)) THEN
+   allocate(sw_ss(nx,ny,12,ntr),stat=err)
+   if (err/=0) then
+      print *,"Erreur à l'allocation du tableau t2m_ss",err
+      stop 4
+   end if
+  end if
   
   if (.not.allocated(pr_ss)) THEN
@@ -274 +305 @@
   if (.not.allocated(t2m_sstopo0)) THEN
    allocate(t2m_sstopo0(nx,ny,12,ntr),stat=err)
+   if (err/=0) then
+      print *,"Erreur à l'allocation du tableau t2m_sstopo0",err
+      stop 4
+   end if
+  end if
+
+  if (.not.allocated(sw_sstopo0)) THEN
+   allocate(sw_sstopo0(nx,ny,12,ntr),stat=err)
    if (err/=0) then
       print *,"Erreur à l'allocation du tableau t2m_sstopo0",err
@@ -333 +372 @@
   real :: TEMP                       !< calcul du nbr de jour < psolid
   real,dimension(nx,ny) :: deltaZs   ! diff temp par rapport a topo orog0 
-  real,dimension(nx,ny,12) :: prmois ! precip sur topo Zs
-  
   real (kind=kind(0.d0)) ::  timeloc
   
@@ -451 +488 @@
     prmois(:,:,m)= pr0(:,:,m)*((pr_sstopo0(:,:,m,itr)*(1-coeftime) + pr_sstopo0(:,:,m,itr+1)*coeftime)/ pr_sstopo0(:,:,m,ntr)) * exp(rappact*(deltaZs(:,:)))
     !$OMP END WORKSHARE
-
 ! Temp et precip fonction de coeftime :
   else if (itr.LT.ntr) then
@@ -466 +502 @@
 enddo 
 
-
+if (smbcomp .EQ. 1) then
+   ! correction d'altitude
+   do m=1,12 
+      if ((itr.LT.ntr-1).and.(ntr.ge.3)) then
+         ! if the is more than 2 snapshots : time climate is mean between 2 snapshots and anomaly with snapshot ntr (ctrl)  
+         !$OMP WORKSHARE
+         swmois(:,:,m)= (sw_sstopo0(:,:,m,itr)*(1-coeftime) + sw_sstopo0(:,:,m,itr+1)*coeftime) - sw_sstopo0(:,:,m,ntr) + sw0(:,:,m)+sw0(:,:,m)*0.00006*(Zs(:,:)-orog0(:,:))
+         !$OMP END WORKSHARE
+      ! Temp et precip fonction de coeftime :
+      elseif (itr.LT.ntr) then
+         !$OMP WORKSHARE
+         swmois(:,:,m)= (sw_sstopo0(:,:,m,itr) - sw_sstopo0(:,:,m,itr+1))*(1-coeftime) + sw0(:,:,m)+sw0(:,:,m)*0.00006*(Zs(:,:)-orog0(:,:))
+         !$OMP END WORKSHARE
+      else ! itr=ntr (time>tsnapmax)
+         !$OMP WORKSHARE
+         !swmois(:,:,m) = sw0(:,:,m)
+         swmois(:,:,m) = sw0(:,:,m)+sw0(:,:,m)*0.00006*(Zs(:,:)-orog0(:,:))
+         !$OMP END WORKSHARE
+      endif
+   enddo
+end if
 
   !$OMP WORKSHARE
@@ -472 +528 @@
   Tjuly=Tmois(:,:,7)
   
-
-
 
 ! calcul de l'accumulation de neige :