Changeset 347 for trunk/SOURCES/climat_forcage_mod.f90

Timestamp:

05/12/21 12:07:40 (3 years ago)

Author:

dumas

Message:

climat_forcage_mod is now compatible with more than 2 GCM snapshots and ntr (number of snapshots) is now defined in param_list.dat

File:

• trunk/SOURCES/climat_forcage_mod.f90 (modified) (8 diffs)

trunk/SOURCES/climat_forcage_mod.f90

@@ -19 +19 @@
   real,dimension(:),allocatable :: spert
   real,dimension(:),allocatable :: bpert
+  
+  integer :: ntr         ! nb of snapshot files now explicitely specified 
 
   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=70), dimension(5) :: forcage_file ! fichier de forcage (climat+topo) => dimension >= ntr
+!  character(len=100) :: forcage_file2 ! fichier de forcage 2 (climat+topo)
 
   character(len=100) :: filforc       ! nom du fichier forcage
@@ -31 +33 @@
 
 
-  integer,parameter :: ntr=2         ! nb of snapshot files now explicitely specified 
-  character(len=200) :: filtr(ntr)   ! fichiers de forcage
-
-  real,dimension(ntr) :: ttr         !< date des tranches (snapshots) 
-  real,dimension(nx,ny,ntr) :: delta, deltj, rapact
+
+!  character(len=200) :: filtr(ntr)   ! fichiers de forcage
+
+  real,dimension(:), allocatable :: ttr         !< date des tranches (snapshots) 
   real,dimension(nx,ny) :: delTatime, delTjtime, rapactime
   real,dimension(nx,ny) :: Zs        !< surface topography above sea level
@@ -49 +50 @@
   
 ! climat des snapshots :
-  real,dimension(nx,ny,12,ntr) :: t2m_ss  !< t2m monthly for every GCM snapshot 
-  real,dimension(nx,ny,12,ntr) :: pr_ss   !< precip monthly for every GCM snapshot
-  real,dimension(nx,ny,ntr) :: orog_ss    !< topo for every GCM snapshot 
-  
-  real,dimension(nx,ny,12,ntr) :: t2m_sstopo0 !< t2m GCM snapshot on orog0
-  real,dimension(nx,ny,12,ntr) :: pr_sstopo0  !< pr GCM snapshot on orog0
+  real,dimension(:,:,:,:), allocatable :: t2m_ss  !< t2m 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 :: pr_sstopo0  !< pr GCM snapshot on orog0
+
+
+!~   real,dimension(nx,ny,12,ntr) :: t2m_ss  !< t2m monthly for every GCM snapshot 
+!~   real,dimension(nx,ny,12,ntr) :: pr_ss   !< precip monthly for every GCM snapshot
+!~   real,dimension(nx,ny,ntr) :: orog_ss    !< topo for every GCM snapshot 
+  
+!~   real,dimension(nx,ny,12,ntr) :: t2m_sstopo0 !< t2m GCM snapshot on orog0
+!~   real,dimension(nx,ny,12,ntr) :: pr_sstopo0  !< pr GCM snapshot on orog0
   
   real :: PYY                        !< constante pour calcul de la fraction solide des precips
@@ -90 +99 @@
   orog0(:,:)=tab2d(:,:)
   
-  filtr(1)=TRIM(DIRNAMEINP)//'Snapshots-GCM/'//TRIM(forcage_file1)
-  filtr(2)=TRIM(DIRNAMEINP)//'Snapshots-GCM/'//TRIM(forcage_file2)
+!  filtr(1)=TRIM(DIRNAMEINP)//'Snapshots-GCM/'//TRIM(forcage_file1)
+!  filtr(2)=TRIM(DIRNAMEINP)//'Snapshots-GCM/'//TRIM(forcage_file2)
 
 ! fichiers donnant l'evolution temporelle
@@ -105 +114 @@
 ! -------------------------------------------------------------
   do k=1,ntr
-     write(6,*) 'Read climate file :',trim(filtr(k))
-     call Read_Ncdf_var('pr',trim(filtr(k)),tab3d)
+     write(6,*) 'Read climate file :',trim(DIRNAMEINP)//'Snapshots-GCM/'//trim(forcage_file(k))
+     call Read_Ncdf_var('pr',trim(DIRNAMEINP)//'Snapshots-GCM/'//trim(forcage_file(k)),tab3d)
      pr_ss(:,:,:,k)=tab3d(:,:,:)
-     call Read_Ncdf_var('tas',trim(filtr(k)),tab3d)
+     call Read_Ncdf_var('tas',trim(DIRNAMEINP)//'Snapshots-GCM/'//trim(forcage_file(k)),tab3d)
      t2m_ss(:,:,:,k)=tab3d(:,:,:)
-     call Read_Ncdf_var('orog',trim(filtr(k)),tab2d)
+     call Read_Ncdf_var('orog',trim(DIRNAMEINP)//'Snapshots-GCM/'//trim(forcage_file(k)),tab2d)
      orog_ss(:,:,k)=tab2d(:,:)
 !~      read(num_forc,*) ttr(k)
@@ -138 +147 @@
 !        print*,nft
   read(num_forc,*) control,nft
-  print*,'control',control,nft
+!  print*,'control :',control,nft
 ! determination of file size (line nb), allocation of perturbation array
 
@@ -205 +214 @@
 subroutine init_forclim
 
-
-  namelist/clim_forcage/clim_ref_file,ttr,forcage_file1,forcage_file2,typerun,lapserate,rappact,filforc,pertbmb,coeft,coefbmb,r_atmvar
+  real,dimension(5) :: 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
 
   rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
   read(num_param,clim_forcage)
-
+  
+  if (ntr.le.size(ttr_temp,dim=1)) then
+   print*,'climat_forcage : we will read ',ntr,' climate snapshots'
+   do k=ntr+1,size(ttr_temp,dim=1)
+      ttr_temp(k)= ttr_temp(ntr)
+      forcage_file(k)=' '
+   enddo
+  else
+   print*,'climat_forcage ERROR : ntr > ttr_temp size'
+   print*,'ttr_temp & forcage_file size must be >= ntr'
+   stop
+  end if
+  
   write(num_rep_42,'(A)')'!___________________________________________________________' 
-  write(num_rep_42,'(A)') '&clim_forcage                                ! nom du bloc '
-  write(num_rep_42,*)
-  write(num_rep_42,'(A,A,A)') 'clim_ref_file = "',trim(clim_ref_file),'"'
-  write(num_rep_42,*) 'ttr = ', ttr
-  write(num_rep_42,'(A,A,A)') 'forcage_file1 = "',trim(forcage_file1),'"'
-  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,'(A,A,A)') ' filforc      = "',trim(filforc),'"'
-  write(num_rep_42,*) 'pertbmb = ', pertbmb
-  write(num_rep_42,*) 'coefT   = ', coefT
-  write(num_rep_42,*) 'coefbmb = ', coefbmb
-  write(num_rep_42,*) 'r_atmvar = ', r_atmvar
-  write(num_rep_42,*)'/'                            
-  write(num_rep_42,*)
+  write(num_rep_42,clim_forcage)
+
+
+  
+  if (.not.allocated(ttr)) THEN
+   allocate(ttr(ntr),stat=err)
+   if (err/=0) then
+      print *,"Erreur à l'allocation du tableau ttr",err
+      stop 4
+   end if
+  end if
+
+  if (.not.allocated(t2m_ss)) THEN
+   allocate(t2m_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
+   allocate(pr_ss(nx,ny,12,ntr),stat=err)
+   if (err/=0) then
+      print *,"Erreur à l'allocation du tableau pr_ss",err
+      stop 4
+   end if
+  end if
+  
+  if (.not.allocated(orog_ss)) THEN
+   allocate(orog_ss(nx,ny,ntr),stat=err)
+   if (err/=0) then
+      print *,"Erreur à l'allocation du tableau orog_ss",err
+      stop 4
+   end if
+  end if 
+  
+  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(pr_sstopo0)) THEN
+   allocate(pr_sstopo0(nx,ny,12,ntr),stat=err)
+   if (err/=0) then
+      print *,"Erreur à l'allocation du tableau pr_sstopo0",err
+      stop 4
+   end if
+  end if
+  
+  do k=1,ntr
+   ttr(k)=ttr_temp(k)
+  enddo
+   
+
+!~   write(num_rep_42,'(A)')'!___________________________________________________________' 
+!~   write(num_rep_42,'(A)') '&clim_forcage                                ! nom du bloc '
+!~   write(num_rep_42,*)
+!~   write(num_rep_42,'(A,A,A)') 'clim_ref_file = "',trim(clim_ref_file),'"'
+!~   write(num_rep_42,*) 'ttr = ', ttr
+!~ !  write(num_rep_42,'(A,A,A)') 'forcage_file= "',trim(forcage_file),'"'
+!~   write(num_rep_42,'(A,A,A)') 'forcage_file1 = "',trim(forcage_file1),'"'
+!~   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,'(A,A,A)') ' filforc      = "',trim(filforc),'"'
+!~   write(num_rep_42,*) 'pertbmb = ', pertbmb
+!~   write(num_rep_42,*) 'coefT   = ', coefT
+!~   write(num_rep_42,*) 'coefbmb = ', coefbmb
+!~   write(num_rep_42,*) 'r_atmvar = ', r_atmvar
+!~   write(num_rep_42,*)'/'                            
+!~   write(num_rep_42,*)
 
   PYY=2.*PI/50.
@@ -362 +445 @@
 ! 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
+    Tmois(:,:,m)=  (t2m_sstopo0(:,:,m,itr)*(1-coeftime) + t2m_sstopo0(:,:,m,itr+1)*coeftime) - t2m_sstopo0(:,:,m,ntr) + t2m0(:,:,m) + deltaZs(:,:)
+    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 :
-  if (itr.LT.ntr) then
+  else if (itr.LT.ntr) then
     !$OMP WORKSHARE
     Tmois(:,:,m)=  (t2m_sstopo0(:,:,m,itr) - t2m_sstopo0(:,:,m,itr+1))*(1-coeftime) + t2m0(:,:,m) + deltaZs(:,:)