source: branches/iLoveclim/SOURCES/climat-perturb_mod-0.4.f90 @ 77

!> \file climat-perturb_mod-0.4.f90
!! Module pour les variations temporelles des variables climatiques
!<

!> \namespace  climat_perturb_mod
!!Module pour les variations temporelles des variables climatiques
!! \author C. Dumas
!! \date 12/2015
!! @note Used modules:
!! @note   - use module3D_phy
!<
module  climat_perturb_mod


use module3d_phy,only:nx,ny,S,S0,Tann,Tjuly,precip,acc,Ylat,num_forc,num_param,num_rep_42,dirforcage,dirnameinp,tafor,time,sealevel,coefbmshelf

implicit none

integer nft             !< nombre de lignes a lire dans le fichier forcage climatique
real,dimension(:),allocatable :: tdate          !< time for climate forcing
real,dimension(:),allocatable :: tpert          !< temperature for climate forcing
real,dimension(:),allocatable :: spert          !< sea surface perturbation

real,dimension(nx,ny) :: ta0          !< initial air temperature at sea level annual

real :: coefT                   !< pour modifier l'amplitude de la perturb. T
real :: rappact                 !< pour le calcul du rapport 'accumulation
integer :: retroac              !< 1-> full retroactions accum
real :: rapbmshelf              !< pour calcul coefbmshelf
real :: mincoefbmelt            !< butoirs mini 
real :: maxcoefbmelt            !< butoirs maxi de coefbmelt
character(len=80) :: filforc    !< nom du fichier forcage

! Pour l'instant tafor est global meme si inutile si on utilise
! un forcage  (variation spatiales)

contains

!--------------------------------------------------------------------------------
!> SUBROUTINE: input_clim  
!! Routine qui permet d'initialiser les variations temporelles des variables climatiques
!>
  subroutine input_clim !routine qui permet d'initialiser les variations temporelles des variables climatiques

    implicit none
    character(len=8) :: control      !label to check clim. forc. file (filin) is usable
    character(len=80):: filin
    integer ::  err                       !< pour l'allocation des tableaux
    integer :: i

    ! Lecture du forcage
    !-----------------------
    ! Le fichier de forcage est lu dans le fichier entree parametres

    filin=trim(dirforcage)//trim(filforc)

    open(num_forc,file=filin,status='old')

    read(num_forc,*) control,nft

    ! Determination of file size (line nb), allocation of perturbation array

    if (control.ne.'nb_lines') then
       write(6,*) filin,'indiquer le nb de ligne en debut de fichier:'
       write(6,*) 'le nb de lignes et le label de control nb_lines'
       stop
    endif

    ! Dimensionnement des tableaux tdate, ....

    if (.not.allocated(tdate)) then
       allocate(tdate(nft),stat=err)
       if (err/=0) then
          print *,"erreur a l'allocation du tableau Tdate",err
          stop 4
       end if
    end if

    if (.not.allocated(spert)) then
       allocate(spert(nft),stat=err)
       if (err/=0) then
          print *,"erreur a l'allocation du tableau Spert",err
          stop 4
       end if
    end if


    if (.not.allocated(tpert)) then
       allocate(tpert(nft),stat=err)
       if (err/=0) then
          print *,"erreur a l'allocation du tableau Tpert",err
          stop 4
       end if
    end if

    do i=1,nft
       read(num_forc,*) tdate(i),spert(i),tpert(i)
    end do
    close(num_forc)

    tpert(:)=tpert(:)*coefT


  end subroutine input_clim

!--------------------------------------------------------------------------------
!> SUBROUTINE: init_forclim
!! Routine qui permet d'initialiser les variables climatiques au cours du temps
!>
subroutine init_forclim


  character(len=80):: filin
  integer :: i,j
  integer :: ivo,jvo

namelist/clim_pert/coefT,rappact,retroac,rapbmshelf,mincoefbmelt,maxcoefbmelt,filforc

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

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

write(num_rep_42,428)'!___________________________________________________________' 
write(num_rep_42,428) '&clim_pert                      ! module climat_perturb_mod '
write(num_rep_42,*)
write(num_rep_42,*) 'coefT        = ', coefT 
write(num_rep_42,*) 'rappact      = ', rappact 
write(num_rep_42,*) 'rretroac     = ', retroac 
write(num_rep_42,*) 'rapbmshelf   = ', rapbmshelf  
write(num_rep_42,*) 'mincoefbmelt = ', mincoefbmelt
write(num_rep_42,*) 'maxcoefbmelt = ', maxcoefbmelt
write(num_rep_42,'(A,A)') ' filforc      = ', filforc
write(num_rep_42,*)'/'                            
write(num_rep_42,*)



!!!!!!!! ATTENTION AJOUTE POUR TEST MAIS A REMETTRE AU PROPRE PLUS TARD C. DUMAS !!!!!!!!!!
!!!!!!!! ancien input_climat_ref de lect_clim_act_anteis
!     accumulation de Philippe
      filin='accumHUY40km.dat'
      call lect_eis(nx,ny,precip,filin,DIRNAMEINP)
!====================================== La reponse est 42 ===========
      write(num_rep_42,*) 'fichier accum : ', filin 

!     cas particulier de Vostok
      ivo=101
      jvo=62
      do j=jvo-1,jvo+1
        do i=ivo-1,ivo+1
            precip(i,j)=0.02         ! valeur plus faible a Vostok.
        end do
      end do
      acc(:,:)=precip(:,:)

!    temperature en surface :
!    parametrisation de Fortuin pour la temperature annuelle.
      do j=1,ny
        do i=1,nx

            if (s0(i,j).le.200.) then   ! shelfs
               tann(i,j)=49.642-0.943*abs(ylat(i,j))
            else if ((s0(i,j).gt.200.).and.(s0(i,j).lt.1500.)) then ! pente
               tann(i,j)=36.689-0.005102*s0(i,j)-0.725*abs(ylat(i,j))
            else if (s0(i,j).ge.1500.) then        ! plateau
               tann(i,j)=7.405-0.014285*s0(i,j)-0.180*abs(ylat(i,j))
            endif

            ta0(i,j)=tann(i,j)
!           pour la temperature d'ete, idem parametrisation huybrechts
            tjuly(i,j)=tann(i,j)-17.65+0.00222*s0(i,j)&
                         +0.40802*abs(ylat(i,j))
        end do
      end do
!!!!!!!! FIN MODIF  TEMPORAIRE !!!!!!!!!!


return
end subroutine init_forclim

!--------------------------------------------------------------------------------
!> SUBROUTINE: forclim
!!  
!!  Routine qui permet le calcule climatiques au cours du temps
!!  @note Au temps considere (time) attribue les scalaires
!!  @note   - tafor : forcage en temperature
!!  @note   - sealevel : forcage niveau des mers
!!  @note   - coefbmelt : forcage fusion basale ice shelves 
!>
subroutine forclim               !  au temps considere (time) attribue les scalaires
  !  tafor : forcage en temperature
  !  sealevel : forcage niveau des mers
  !  coefbmelt : forcage fusion basale ice shelves 

!  use module3d_phy
  implicit none
  integer :: i,ift

  !       time en dehors des limites du fichier forcage
  if(time.lt.tdate(1)) then
     tafor=tpert(1)
     sealevel=spert(1)
     ift=1

  else if (time.ge.tdate(nft)) then
     tafor=tpert(nft)
     sealevel=spert(nft)
     ift=nft

  else
     do i=1,nft-1 
        if((time.ge.tdate(i)).and.(time.lt.tdate(i+1))) then ! entre i et i+1 : cas general
           tafor=tpert(i)+(tpert(i+1)-tpert(i))*       &
                (time-tdate(i))/(tdate(i+1)-tdate(i))
           sealevel=spert(i)+(spert(i+1)-spert(i))*    &
                (time-tdate(i))/(tdate(i+1)-tdate(i))
           ift=i
           goto 100
        endif
     end do
  endif
100 continue


  !  coefmshelf est un coefficient qui fait varier bmgrz et bmshelf
  !  en fonction de tafor 

  !   coefbmshelf=(1.+tafor/10.)           ! coefbmshelf=0 pour tafor=-10deg

  if ((tafor.le.0).and.(tafor.gt.-5.)) then
     coefbmshelf=(1.+tafor/rapbmshelf)     ! coefbmshelf=0 pour tafor=-7 standard precedent

  else if (tafor.le.-5) then               ! lineaire en 0 a -10 et la valeur a -5
     coefbmshelf=(1.-5./rapbmshelf)/5.*(tafor+10)

  else
     coefbmshelf=(1.+5.*tafor/rapbmshelf)  ! 5 fois plus efficace vers le chaud
  endif

  coefbmshelf=max(coefbmshelf,mincoefbmelt)  
  coefbmshelf=min(coefbmshelf,maxcoefbmelt)

        call massb_perturb_Tparam  
  !  massb_perturb_Tparam remplace desormais massb_anteis_perturb 
  !  le code est le meme pour l'Antarctique (simplement passe en format libre)
  !  dans le fichier massb-ant_perturb_Tparam.f90


end subroutine forclim

subroutine massb_perturb_Tparam               ! calcule le mass balance en mode perturbation
                                              ! avec la temperature parametree
                                              ! version pour l'antarctique
                                              ! simple copie nettoyee de massb_anteis_perturb
  implicit none

  integer :: i,j

!     surface temperature et accumulation

  do j=1,ny
     do i=1,nx

        if(retroac.eq.1) then

           tann(i,j)=ta0(i,j)-0.00914*(s(i,j)-s0(i,j))+tafor
           tjuly(i,j)=tann(i,j)-17.65+0.00222*s(i,j)  &
                +0.40802*abs(ylat(i,j))
         
           acc(i,j)=precip(i,j)*exp(rappact*(tann(i,j)-ta0(i,j)))

        else if(retroac.eq.0) then

           tann(i,j)=ta0(i,j)
           tjuly(i,j)=tann(i,j)
           acc(i,j)=precip(i,j)*exp(rappact*(tann(i,j)-ta0(i,j)))
        endif
     end do
  end do

end subroutine massb_perturb_Tparam


end module  climat_perturb_mod