source: trunk/SOURCES/GrIce2sea_files/massb-GrIce2sea_RCM.f90

!> \file massb-GrIce2sea_RCM.f90
!!  Calcule le mass balance d'apres un fichier de snapshots
!!  avec la temperature parametree
!<

!> SUBROUTINE: massb_Ice2sea_RCM   
!! \author Cat
!! \date 1 er juillet 2012
!! @note  Routine qui utilise un fichier de snapshots
!! et la temperature en surface avec un lapse rate
!! version groenland voir plus tard si convient Antarctique
!! 
!! @note Used modules:
!! @note    - use module3D_phy
!! @note    - use climat_Grice2sea_mod
!<



subroutine massb_Ice2sea_RCM                ! calcule le mass balance 

use module3D_phy
use climat_Grice2sea_years_mod

implicit none
integer             :: k_snap               ! pour calculer les indices de temps
real                :: time_bis             ! pour repliquer les dernieres annees
!     surface temperature et accumulation

Tann (:,:) = Ta0 (:,:) + T_lapse_rate * (S(:,:)-S0(:,:)) 
Ts(:,:)    = Tann(:,:)



! calcule bm_time par interpolation entre deux snapshots
! avant prend la valeur de reference
! apres prend la derniere valeur

! en general les snapshots vont de 2000 a 2200


  if(time.lt.time_snap(1)) then              ! time avant le forcage
     bm_time(:,:) = bm_0(:,:)
     k_snap       = 1
     S_ref(:,:)   = S(:,:)                   ! du coup sera la surface de reference avant le forcage
     icum         = 0
     i_moy        = 0
     bm_ref_t(:,:)=  bm_0(:,:)               ! bilan de masse de reference
     time_prec    = time


  else if (time.ge.time_snap(nb_snap)) then  ! time apres le forcage
     k_snap       = nb_snap
     bm_time(:,:) =  smb_snap (:,:,k_snap)

     if (abs(time-time_prec-1.).lt.dt) then   ! 
        time_prec = time_prec + 1
        icum      = 1
     else
        icum      = 0
     endif


  else                                       ! cas general
 
     do k = 1 , nb_snap-1

        if((time.ge.time_snap(k)).and.(time.lt.time_snap(k+1))) then ! entre k et k+1 
           bm_time(:,:) = smb_snap(:,:,k)+(smb_snap(:,:,k+1)-smb_snap(:,:,k)) *   &
                            (time-time_snap(k))/(time_snap(k+1)-time_snap(k))


! exactement sur le snapshot et avec un ecart 1 an par rapport au precedent stockage
!           write(6,*) 'time,tests',k,time,time-time_snap(k),time-time_prec-1.
             if ((abs(time-time_snap(k)).le.dt).and.(abs(time-time_prec-1.).lt.dt)) then    
                k_snap    = k
                icum      = 1
                time_prec = time_snap(k)                                ! time_prec est le temps du precedent
                                                                        ! snapshot garde
             else
                icum    = 0
             endif

           exit
        endif

     end do

  endif


call  grad_smb              !-----------------------------> A faire

! ajuste bm en fonction du temps et du gradient

bm(:,:) = bm_time(:,:) + grad_bm(:,:) *(S(:,:) - S_ref(:,:))

write(6,897) time, time_prec, icum, i_moy
897   format('test temps smb   ',2(f0.3,1x),2(i0,1x))

! garde les 10 dernieres annees et calcule la moyenne

if (icum.eq.1) then                           ! stockage dans le tableau bm_ref_10
   do k = 9,1,-1
      bm_ref_10(:,:,k+1) = bm_ref_10(:,:,k)   ! on decale tous les elements
   end do
   bm_ref_10(:,:,k+1) = bm(:,:)               ! le plus recent est en position 1
   i_moy              = i_moy +1              ! compte combien il y en a pour la moyenne
   i_moy              = min(i_moy,10)
   
    bm_ref_t(:,:)   = 0.
   do k = 1,i_moy
      bm_ref_t(:,:) = bm_ref_t(:,:) + bm_ref_10(:,:,k)
   end do
   bm_ref_t(:,:)    = bm_ref_t(:,:)/i_moy

   write(6,898) time, time_prec, icum, i_moy
898   format('cumul pour gradient  ',2(f0.3,1x),2(i0,1x))
end if


!     ablation (fonction de T et acc) est maintenant appelee dans le main
!     et dans steps_time_loop.
!     On fait donc une routine fake ici

!      massb_Ice2sea_fixe est appeles par forclim du module climat_Grice2sea_mod

      
debug_3D(:,:,29)  = Tann(:,:)-Ta0(:,:)
debug_3D(:,:,30)  = Acc(:,:)-precip(:,:)
debug_3D(:,:,97) = bm_ref_t(:,:)
debug_3D(:,:,98) = grad_bm(:,:) 
debug_3D(:,:,99) = bm_time(:,:)    
return
end subroutine massb_Ice2sea_RCM
!------------------------------------------------------------------------------
!> initialise le calcul du gradient vertical de smb
subroutine init_grad_smb 


use module3D_phy
use climat_Grice2sea_years_mod
implicit none

character(len=120)       ::  file_grad_smb ! nom du fichier gradients de Tamsin
character(len=40)        ::  fin_ligne     ! fin de la ligne
real                     ::  grad
        
namelist/grad_smb/file_grad_smb

! Dans lequel sont :
! grad_N_smb_pos,grad_N_smb_neg,grad_S_smb_pos,grad_S_smb_neg,lim_lat

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


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

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

    write(num_rep_42,428)'!___________________________________________________________' 
    write(num_rep_42,428)'!  gradient smb GrIce2sea gradients_18_07_2012                                  '
    write(num_rep_42,grad_smb)
    write(num_rep_42,428)'!___________________________________________________________' 


    file_grad_smb=trim(dirnameinp)//'SMB-H-Feedback/gradients_18_07_2012/'//trim(file_grad_smb)
    open(622,file=file_grad_smb)
    do i=1,4
       read(622,'(f9.4,A)') grad,fin_ligne
       write(6,*) grad,fin_ligne

       if (index(fin_ligne, "North").ne.0) then       ! North est dans la ligne
          if (index(fin_ligne, "<").ne.0)  then       ! smb negatif
             grad_N_smb_neg = grad
          else if (index(fin_ligne, ">=").ne.0) then  ! smb positif
             grad_N_smb_pos = grad
          else
             write(6,*) 'probleme lecture North fichier smb', file_grad_smb
             STOP
          end if

       else if (index(fin_ligne, "South").ne.0) then   !South est dans la ligne
          if (index(fin_ligne, "<").ne.0)  then        ! smb negatif
             grad_S_smb_neg = grad
          else if (index(fin_ligne, ">=").ne.0) then   ! smb positif
             grad_S_smb_pos = grad
          else
             write(6,*) 'probleme lecture South fichier smb  ', file_grad_smb
             STOP
          end if

       else 
          write(6,*) 'probleme lecture ni North ni South fichier smb  ', file_grad_smb
          write(6,*) 'fin_ligne',fin_ligne,'   index North', index(fin_ligne, "North"),'   index South', index(fin_ligne, "South")
          STOP
       end if
    end do

    write(6,*) 'coefficients lus '
    write(6,*) 'grad_N_smb_pos', grad_N_smb_pos
    write(6,*) 'grad_N_smb_neg', grad_N_smb_neg
    write(6,*) 'grad_S_smb_pos', grad_S_smb_pos
    write(6,*) 'grad_S_smb_neg', grad_S_smb_neg

    grad_N_smb_pos = coef_smb_unit * grad_N_smb_pos
    grad_N_smb_neg = coef_smb_unit * grad_N_smb_neg
    grad_S_smb_pos = coef_smb_unit * grad_S_smb_pos
    grad_S_smb_neg = coef_smb_unit * grad_S_smb_neg

    write(6,*) 'coefficients *coef_smb_unit '
    write(6,*) 'grad_N_smb_pos', grad_N_smb_pos
    write(6,*) 'grad_N_smb_neg', grad_N_smb_neg
    write(6,*) 'grad_S_smb_pos', grad_S_smb_pos
    write(6,*) 'grad_S_smb_neg', grad_S_smb_neg


end subroutine init_grad_smb


!------------------------------------------------------------------------------
!> Calcule le gradient vertical de smb

subroutine grad_smb 

use module3D_phy
use climat_Grice2sea_years_mod
implicit none

do j = 1,ny
   do i = 1,nx
      if (Ylat(i,j).gt.77.) then                     ! region nord
         if (bm_ref_t(i,j).lt. 0.) then              ! smb negatif
            grad_bm (i,j) = grad_N_smb_neg

         else if (bm_ref_t(i,j).ge. 0.) then         ! smb positif
            grad_bm (i,j) = grad_N_smb_pos 

         end if

      else if (Ylat(i,j).le.77.) then                ! region sud
         if (bm_ref_t(i,j).lt. 0.) then              ! smb negatif
            grad_bm (i,j) = grad_S_smb_neg

         else if (bm_ref_t(i,j).ge. 0.) then         ! smb positif
            grad_bm (i,j) = grad_S_smb_pos 

         end if
      end if
   end do
end do

end subroutine grad_smb
!------------------------------------------------------------------------------