source: trunk/SOURCES/climat_forcage_mois_mod.f90 @ 23

module climat_forcage_mois_mod                   

! forcage avec champs mensuels
! equilibre climatique avec 1 seul snapshot climatique
! lecture fichiers snapshots climat Temp et Precip mensuels (2 fichiers)
! nouvelle version avec liste des variables utilisées par le module
! C. Dumas 02/2015

USE module3D_phy,only:nx,ny,S,S0,slv,Tann,Tjuly,Tmois,acc,coefbmshelf,ro,num_param,num_rep_42,dirnameinp                        
use netcdf
use io_netcdf_grisli
!USE printtable
       
implicit none

! 1=decalaration variables
!-------------------------

!integer :: nft                 !   NFT est le nombre de lignes a lire dans le fichier contenant le forcage climatique

                        
integer,parameter :: nbmois=12
!integer,parameter :: ntr=1     ! nb of snapshot files : 1 pour stationnaire - ntr is now explicitely specified in the climat module
!integer :: itr        ! index of snapshot                      
!real,dimension(:),allocatable :: TDATE                 ! time for climate forcing
!real,dimension(:),allocatable :: indexco2              ! index en co2 du forcage
!real :: co2                                            ! co2 a l'instant t
!real :: coefindex                                       ! pour le calcul du climat en fn de l'index
!real,dimension(:),allocatable :: alphaT                 ! index temperature
!real,dimension(:),allocatable :: alphaP                 ! index precipitation
!real,dimension(:),allocatable :: SPERT                  ! sealevel forcing

!real,dimension(ntr) :: ttr                             ! le co2 des snapshots  (len=ntr) 
!real :: ttrfile                                         ! pour lecture du co2 des snapshots
!real alphaTTR(ntr) ! Le alphaT de l'index glaciologiq. aux snapshots

real,dimension(nx,ny)      :: Ssnap                 ! altitude de surface dans le snapshot
real,dimension(nx,ny,nbmois) :: Tm                    ! temperature mensuelle
real,dimension(nx,ny,nbmois) :: Pm                    ! precipitation mensuelle

real,dimension(nx,ny,nbmois) :: Tm_surf                 ! surface temperature (after topo. correction) 13 nbmois pour modulo
real,dimension(nx,ny,nbmois) :: Pm_surf                 ! surface precipitation (after topo. correction)

real,dimension(nx,ny) :: ZS                             !< surface topography above sea level

real,dimension(nx,ny,nbmois) :: lapserate                 ! lapse rate
real :: psolid=2.                                      ! temp limit between liquid and solid precip

character(len=200)  :: filin                            ! nom temporaire
!character(len=100) :: file_temporel                    ! forcage temporel
!character(len=100),dimension(ntr) :: filtr_t, filtr_p   ! snapshot file name
character(len=100) :: filtr_t1,filtr_p1                 ! snapshot file name


contains

! 2=lecture des inputs
!--------------------

subroutine input_clim                                   ! routine qui permet d'initialiser les variables climatiques

! variables locales
!-------------------

  implicit none
  integer :: i,j,mo,ti,tj  
  character(len=8) :: control   ! label to check clim. forc. file (filin) is usable
  integer :: l                  ! In snapshot files:the first column is the mask, read but not used
  real :: T_surf_ref            ! variable de travail calcul temp a l'instant t a la surface S
  character(len=100) :: file_ncdf      !< fichier netcdf issue des fichiers .dat
  real*8, dimension(:,:,:), pointer :: data_3D => null() ! donnees lues dans le netcdf

! lecture des fichiers snapshots pour tout geoplace
! -------------------------------------------------
  write(6,*) 'lecture fichiers snapshots'

!temperature
  filin=TRIM(dirnameinp)//TRIM(filtr_t1)
  call Read_ncdf_var('t2m',trim(filin),data_3D)    ! Temperature
  Tm(:,:,:)=data_3D(:,:,:) 
! ancienne lecture fichiers txt
!!$  filin=trim(dirnameinp)//'climat_forcage_mois_file/'//trim(filtr_t1)
!!$  write(6,*) 'lecture fichier ', trim(filin)
!!$
!!$  open(20,file=trim(filin))
!!$  do j=1,ny 
!!$     do i=1,nx
!!$        read(20,*) ti, tj, (Tm(i,j,mo),mo=1,12)
!!$     end do
!!$  end do
!!$  close(20)

!precipitation
  filin=TRIM(dirnameinp)//TRIM(filtr_p1)
  call Read_ncdf_var('pr',trim(filin),data_3D)    ! precipitation
  Pm(:,:,:)=data_3D(:,:,:) 

! conversion en m/an :
! Fichier ERA' en mm/jour
! FIchier NCEP2 en mm/s
  do mo=1,nbmois
     if (mo.eq.1.or.mo.eq.3.or.mo.eq.5.or.mo.eq.7.or.mo.eq.8.or.mo.eq.10.or.mo.eq.12) then
!        Pm(:,:,mo)= Pm(:,:,mo)*31.  ! ERA int
        Pm(:,:,mo)= Pm(:,:,mo)*86400.*31./1000. ! NCEP2
     elseif (mo.eq.4.or.mo.eq.6.or.mo.eq.9.or.mo.eq.11) then
!        Pm(:,:,mo)= Pm(:,:,mo)*30.  ! ERA int
        Pm(:,:,mo)= Pm(:,:,mo)*86400.*30./1000. ! NCEP2
     elseif (mo.eq.2) then
!        Pm(:,:,mo)= Pm(:,:,mo)*28.25  ! ERA int
        Pm(:,:,mo)= Pm(:,:,mo)*86400.*28.25/1000. ! NCEP2
     else
        print*,'ERREUR CONVERSION Precip input_clim climat_forcage_mois_mod'
     endif
  enddo


! ancienne lecture fichiers txt
!!$  filin=trim(dirnameinp)//'climat_forcage_mois_file/'//trim(filtr_p1)
!!$  write(6,*) 'lecture fichier ', trim(filin)
!!$
!!$  open(20,file=trim(filin))
!!$  do j=1,ny 
!!$     do i=1,nx
!!$        read(20,*) ti, tj, (Pm(i,j,mo),mo=1,12)
!!$     end do
!!$  end do
!!$  close(20)



!print*,'dans climat-forcage-mois'
!print*,'Tm et Pm', Tm(100,100,3,1), Pm(100,100,3,1)

!  Correction d'altitude pour que la temperature des snapshots soit celle
!  sur l'altitude de reference S0

!  ATTENTION. Pour l'instant Ssnap(:,:,itr)= S topo initiale
!  ensuite il faudra lire Ssnap dans le  fichier forcage
  Ssnap(:,:)=S0(:,:)

end subroutine input_clim



!-------------------------------------------------------------------------------
! lecture de la namelist
subroutine init_forclim

  real,dimension(12) :: lect_lapse       ! lecture lapse rate
  integer :: mo

! fichiers snapshots
  namelist/snap_forcage_mois/filtr_t1,filtr_p1

! dmr moved up since namelist declaration should be in the header
  namelist/lapse_month/lect_lapse

! formats pour les ecritures dans 42
428 format(A)
  rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
  read(num_param,snap_forcage_mois)
  write(num_rep_42,428)'!___________________________________________________________' 
  write(num_rep_42,428) '&snap_forcage_mois                            ! module climat_forcage_mois_mod'
  write(num_rep_42,'(A,A)')   'filtr_t1 = ', filtr_t1
  write(num_rep_42,'(A,A)')   'filtr_p1 = ', filtr_p1
  write(num_rep_42,*)'/'                      
  write(num_rep_42,428) '! fichiers temperature et precip : 12 mois'
  write(num_rep_42,*)

! lapserate
  rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
  read(num_param,lapse_month)
  write(num_rep_42,428)'!___________________________________________________________' 
  write(num_rep_42,428) '&lapse_month                                   !  module climat_forcage_mois_mod'
  write(num_rep_42,'(A,12(f0.2,","))')  'lapse_month  = ', lect_lapse
  write(num_rep_42,*)'/'                      
  write(num_rep_42,428) '! laspe rates janvier -> decembre en deg/km'
  write(num_rep_42,*)

! pour repasser en deg/m et copier dans lapserate
  do mo=1,nbmois
     lapserate(:,:,mo)=lect_lapse(mo)/1000.
  enddo

! coefbmshelf coefficient pour la fusion basale sous les ice shelves
! ici pas d'evolution temporelle donc coefbmshelf est fixe
  coefbmshelf=1.

end subroutine init_forclim


!---------------------------------------------------------------------
!forcage climatique au cours du temps
subroutine forclim
  
  integer i,j,mo

!  Correction d'altitude pour la temperature y compris sur les lacs
! Zs est l'altitude de la surface qu'elle soit mer, glace ou lac

  do j=1,ny
     do i=1,nx
        ZS(I,J)=max(slv(i,j),S(I,J))
        do mo=1,nbmois
           Tm_surf(i,j,mo)= - lapserate(i,j,mo) * (Zs(i,j)-S0(i,j)) + Tm(i,j,mo) 
           Pm_surf(i,j,mo)= Pm(i,j,mo)*exp(0.05*(Tm_surf(i,j,mo)-Tm(i,j,mo)))
        end do
     end do
  end do

!print*,'dans climat-forcage-mois aprec correction altitude surface'
!print*,'Tm_surf et Pm_surf :', Tm_surf(100,100,3), Pm_surf(100,100,3)


! Tm_surf dans Tmois pour visibilite dans tous le code
  Tmois(:,:,:)=tm_surf(:,:,:)

! calcul de Tann et Tjuly pour les sorties :
  Tann(:,:)=sum(Tmois,dim=3)/12.     ! moy annuelle
  Tjuly(:,:)=Tmois(:,:,7)            ! temp juillet

  acc(:,:)=sum(Pm_surf,dim=3,mask=Tmois < psolid) ! /12.
  acc(:,:)=acc(:,:)*1000./ro

END subroutine forclim


end module climat_forcage_mois_mod