source: trunk/SOURCES/Greeneem_files/lect-clim-act-greeneem_mois_mod.f90 @ 4

module lect_clim_act_greeneem_mois                    


use module3D_phy                        
use printtable
use declare_month
use climat_perturb_mois_mod


implicit none



contains


subroutine input_climat_ref()           ! routine qui permet d'initialiser les variables climatiques




real,dimension(nx,ny,mois) :: Tm      ! temperature mensuelle sur la calotte de ref
real,dimension(nx,ny,mois) :: Pm      ! precipitation mensuelle sur la calotte de ref


character(len=200)  :: filin                            ! nom temporaire
character (len=200) :: filtr_topo, filtr_t, filtr_p
real :: T_surf_ref                                      ! variable de travail

real :: ti, tj                     ! variable provisoire de lecture


namelist/clim_init_mois/filtr_topo,filtr_t,filtr_p    

! lecture par namelist 
!---------------------


! formats pour les ecritures dans 42
428 format(A)

rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
read(num_param,clim_init_mois)


write(num_rep_42,428)'!___________________________________________________________' 
write(num_rep_42,428) '&clim_perturb_mois                            ! module climat_perturb_mois_mod'
write(num_rep_42,'(A,A)')   'filtr_topo =', filtr_topo
write(num_rep_42,'(A,A)')   'filtr_t =', filtr_t
write(num_rep_42,'(A,A)')   'filtr_p =', filtr_p
write(num_rep_42,*)'/'                      
write(num_rep_42,428) '! fichiers temperature et precip : 12 mois'
write(num_rep_42,428) '! faire un bloc namelist par snapshot'
write(num_rep_42,*)



call lect_lapserate_months       ! lit les lasperate
                                 ! version actuelle : on utilise un
                                 ! lapserate parametre homogene annuel
                                 ! on le passe en mensuel avec un cosinus



! lecture des fichiers snapshots pour tout geoplace
! -------------------------------------------------
write(6,*) 'fichiers snapshots'
 
   ! lecture de la topo de surface de reference du GCM
   filin=trim(dirnameinp)//trim(filtr_topo)
   write(6,*) filin
   open(20,file=filin)
   do j=1,ny 
      do i=1,nx
         read (20,*) S0CLIM(i,j)
         S0CLIM(i,j)=max(S0CLIM(i,j),slv(i,j))   ! pour etre au niveau des mers
      end do
   end do
   close(20)


   ! temperature : temperature a l'altitude de la calotte de ref (MAR)
   filin=trim(dirnameinp)//trim(filtr_t)
   write(6,*) filin
   open(20,file=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)

! test aurel: on garde la temp parametree ->
!!$   do mo=1,12
!!$      Tm(:,:,mo)=Tm_surf(:,:,mo)
!!$   end do


   ! precipitation a l'altitude de la calotte de ref (MAR)
   filin=trim(dirnameinp)//trim(filtr_p)
   write(6,*) filin
   open(20,file=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)

!!$! test aurel pour restreindre les precipes a la calotte
!!$   do j=1,ny
!!$      do i=1,nx
!!$         if(H(i,j).le.1.) then
!!$            Pm(i,j,:)=0.
!!$         end if
!!$      end do
!!$   end do

   do j=1,ny
      do i=1,nx
! Zs est l'altitude de la surface qu'elle soit mer, glace ou lac
         ZS(i,j)=max(slv(i,j),S0CLIM(I,J)) !slv=sealev pour nolake

         do mo=1,mois
            Tm_0(i,j,mo)=Tm(i,j,mo)      ! temp sur la calotte de ref du GCM a l'instant initial

! on corrige la temperature initiale par l'altitude
            T_surf_ref = Tm(i,j,mo) + lapserate(i,j,mo) * (ZS(i,j)-S0CLIM(i,j))

! on applique le forcage
            T_surf_ref = T_surf_ref+Tafor

            if(retroac.eq.1) then
! on suppose que les precipes sont fournies a la meme altitude que les temperatures
               Pm(i,j,mo)= Pm(i,j,mo)*exp(rappact*(T_surf_ref-Tm(i,j,mo)))
            else if (retroac.eq.0) then
               Pm(i,j,mo)=Pm(i,j,mo)
            end if

            Tm(i,j,mo)=T_surf_ref

         end do
      end do
   end do



do mo =1,12
   Tm_surf(:,:,mo)=Tm(:,:,mo)
   Pm_surf(:,:,mo)=Pm(:,:,mo)
end do
Tm_surf(:,:,mois+1)=Tm_surf(:,:,1)
!Tm_0(:,:,mois+1)=Tm_0(:,:,1)
Pm_surf(:,:,mois+1)=Pm_surf(:,:,1)

do j=1,ny
   do i=1,nx
      Tann(i,j)=(sum(Tm_surf(i,j,:))-Tm_surf(i,j,1))/12 ! janvier compte en double
      Tjuly(i,j)=Tm_surf(i,j,7)
      !Tjuly(i,j)=(Tm_surf(i,j,6)+Tm_surf(i,j,7)+Tm_surf(i,j,8))/3 ! en fait un Tjja
   end do
end do

call massb_perturb_mois  !modif aurel 16 juil 2010

print*, 'fini input_climat_ref'



end subroutine input_climat_ref




!------------------------------------------------------------------------------------------------------------
subroutine lect_lapserate_months  ! lapserates mensuels mais uniformes spatialement
! cette routine est susceptible de changer en fonction de la facon de lire les lapserates.
! pour le moment : on peut lire un fichier contenant les lapserates non homogènes spatiallement.



character (len=200) :: filtr_lr
integer :: k

! pour la parametrisation de Fausto et al. 2009
real :: d_ann, gamma_ann, c_ann, kappa_ann
real :: d_july, gamma_july, c_july, kappa_july

integer :: gmint  ! choix de la parametrisation temperature
                  ! 0 -> Fausto et al. 2009
                  ! autre -> ancienne version Gmint

real :: ti, tj                     ! variable provisoire de lecture


gmint=0
!gmint=1

! on laisse la possibilite ed lire un fichier avec une carte de lapserate
! mais en realite pas encore operationnel -> 
! on utilise la parametrisation de Fausto et al. 2009
namelist/lapse_month/filtr_lr


! lecture de la namelist
! formats pour les ecritures dans 42
428 format(A)

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_perturb_mois_mod'
write(num_rep_42,'(A,A)')   'filtr_lr =', filtr_lr
write(num_rep_42,*)'/'                      
write(num_rep_42,428) '! lapse rates janvier -> decembre en deg/m'
write(num_rep_42,*)


open(20,file=trim(dirnameinp)//trim(filtr_lr))
do j=1,ny
   do i=1,nx
      read(20,*) ti, tj, (lapserate(i,j,mo),mo=1,12) ! en deg/m
   end do
end do
close(20)

! les runs avec un lapserate non homogène donnent de mauvais résultats
! on choisit des lapserates de Fausto et al. 2009 (actualisation de gmint)
! T = d + gamma x z + c x lat + kappa x long
!annuel et juillet (resp) :
!d (°C)           41.83              14.70
!gamma  (°/km)        -6.309             -5.426
!c (°C/°N)           -0.7189            -0.1585
!kappa (°C/°W)        0.0672             0.0518
! sinon on peut utiliser l'ancienne parametrisation eismint

if(gmint.eq.0) then
d_ann = 41.83
d_july = 14.70
gamma_ann = -6.309/1000
gamma_july = -5.426/1000
c_ann = -0.7189
c_july = -0.1585
kappa_ann = 0.0672 
kappa_july = 0.0518

else
! gmint :
d_ann=49.13
d_july=30.78
gamma_ann=-7.992/1000
gamma_july=-6.277/1000
c_ann=-0.7576
c_july=-0.3262
kappa_ann = 0. 
kappa_july = 0.
end if

! attention aux longitudes, dans Fausto elles sont en °W
do j=1,ny
   do i=1,nx
      ZS(i,j)=max(slv(i,j),S(i,j))
      Tann(i,j)=d_ann+gamma_ann*S(i,j)+c_ann*ylat(i,j)+kappa_ann*(360-xlong(i,j))
      Tjuly(i,j)=d_july+gamma_july*S(i,j)+c_july*ylat(i,j)+kappa_july*(360-xlong(i,j))
   end do
end do


!!$! TEST DE SENSIBILITE SUR LE LAPSERATE :
!!$! JE GARDE LA TEMP PARAMETREE FAUSTO MAIS JE MODIFIE ENSUITE LE GAMMA
!!$gamma_ann = gamma_ann+gamma_ann*30/100 !-6.309/1000
!!$gamma_july = gamma_july+gamma_july*30/100     !-5.426/1000
!!$print*, 'lapse ann et july  ', gamma_ann, gamma_july

! il faut mensualiser maintenant
do mo=1,12
   do j=1,ny
      do i=1,nx
         Tm_surf(:,:,mo)=tann(:,:)+(tjuly(:,:)-tann(:,:))*cos(2*pi*(mo+5)/12) 
         ! on calcule aussi lapserate mensuel, mais homogene spatialement
         lapserate(:,:,mo)=gamma_ann+(gamma_july-gamma_ann)*cos(2*pi*(mo+5)/12)
      end do
   end do
end do

!aurel test : pas de lapse
!lapserate(:,:,:)=0.
! pour les modulos


end subroutine lect_lapserate_months
!--------------------------------------------------------------------------------------------------------



end module lect_clim_act_greeneem_mois