source: trunk/SOURCES/initial-phy-2.f90 @ 237

!> \file initial-phy-2.f90
!! Initialisation des parametres physique du modele
!<

!> SUBROUTINE: initial_phy
!! Subroutine qui initialise les parametres physique du modele
!! \author ...
!! \date ...
!! @note Used modules:
!! @note    - module3D_phy
!!
!<
subroutine initial_phy

  !     **************************************************
  !     **         initialisation                        *
  !     *                                                *
  !     **************************************************

  use module3D_phy
  use param_phy_mod

  implicit none

  character(len=80) :: filin
  character(len=80) :: comment_run

  namelist/runpar/runname,icompteur,iout,reprcptr,itracebug,num_tracebug,comment_run
  namelist/grdline/igrdline,Schoof,ibmelt_inv
  namelist/timesteps/dtmin,dtmax,dtt,testdiag,tbegin,tend


  ! ouverture du fichier parametre ..._param_list.dat
  !--------------------------------------------------------------
  dirsource='../SOURCES'
  ! dirsource='../../SOURCES'
  ! filin='../'//trim(dirsource)//'/Fichiers-parametres/'//trim(geoplace)//'_param_list.dat'
  filin=trim(geoplace)//'_param_list.dat'

  open(num_param,file=filin) 


  ! lecture des parametres du run                          block runpar
  !--------------------------------------------------------------------

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

!  filin='../Param/param'//runname
  filin='param'//runname
  open(num_rep_42,file=filin)

  ! reecriture des parametres entres dans le fichier 42
  !------------------------------------------------------
  write(num_rep_42,*)'!  Parametres du run :',runname, '    ',geoplace
  write(num_rep_42,*)'!  sources dans  ',trim(dirsource)
  write(num_rep_42,*)
  write(num_rep_42,'(A)')'!___________________________________________________________' 
  write(num_rep_42,'(A)') '&runpar                  ! nom du bloc  parametres du run'
  write(num_rep_42,*)
  write(num_rep_42,*) 'runname      =  "',runname,'"'
  write(num_rep_42,*) 'icompteur    = ', icompteur                              
  write(num_rep_42,*) 'iout         = ', iout 
  write(num_rep_42,*) 'reprcptr     = "',trim(reprcptr),'"'
  write(num_rep_42,*) 'itracebug    = ', itracebug
  write(num_rep_42,*) 'num_tracebug  = ', num_tracebug
  write(num_rep_42,*) 'comment_run   = "',trim(comment_run),'"'
  write(num_rep_42,*)'/'                            
  write(num_rep_42,'(A)') '! runname    : nom de l experience (8 caracteres)  '
  write(num_rep_42,'(A)') '! icompteur  : reprise dans un fichier  0 -> non, 1 -> oui, 2 -> T et Hw'
  write(num_rep_42,'(A)') '! icompteur  : 3-> T seulement'
  write(num_rep_42,'(A)') '! iout   : 1-> sortie cptr pour reprise,2 -> sortie nc pour reprise '
  write(num_rep_42,'(A)') '! reprcptr   : nom du fichier restart'
  write(num_rep_42,'(A)') '! itracebug: 1-> ecriture de traces au debut des routines'
  write(num_rep_42,'(A)') '! num_tracebug   numero d unite ecriture traces'
  write(num_rep_42,'(A)') '! comment_run: commentaire court sur le run'

  write(num_rep_42,*)



  !     switch reprise compteur ou non: 1 => compteur, 0 => pas compteur
  !     icompteur=2 reprise de tout sauf de la topo
  !     icompteur=1 reprise de tout

  !      reprcptr=TRIM('../CPTR-Heino/Heino096+k150.CPTR')




  write(6,*)'_________________________________________________'
  write(6,*)'runname=',runname
  !      write(6,*) dirnameout
  !      write(6,*) TRIM(DIRNAMEOUT)//'short'//runname//'.ritz'


  ! ouverture du fichier short
  !------------------------------

  ! filin='time-series/short'//runname//'.ritz'
  filin='short'//runname//'.ritz'
  filin=TRIM(DIRNAMEOUT)//TRIM(filin)

  open(num_ritz,file=filin)



  ! lecture des parametres du run                          block grdline
  !--------------------------------------------------------------------


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


  write(num_rep_42,'(A)')'!___________________________________________________________'  
  write(num_rep_42,'(A)') '&grdline                 ! bloc grounding line'
  write(num_rep_42,*)
  write(num_rep_42,*) 'igrdline     = ',igrdline
  write(num_rep_42,*) 'Schoof       = ',Schoof
  write(num_rep_42,*) 'ibmelt_inv   = ',ibmelt_inv
  write(num_rep_42,*)'/'                           
  write(num_rep_42,'(A)')'! igrdline   : 1 ligne d echouage fixée, sinon 0'
  write(num_rep_42,'(A)')'! Schoof     : 0 pas de Schoof, 1 flux de Schoof'
  write(num_rep_42,'(A)')'! ibmelt_inv : 0 cas std, 1 inversion du bmelt (avec igrdline=1)'
  write(num_rep_42,*)


  !     grounding line fixee
  !--------------------------
  !     Pour faire des expériences avec l'épaisseur des ice-shelves fixée
  !     igrdline=1 !fixé, sinon (libre) igrdline=0
  !     igrdline=0 



  ! lecture des parametres du run                          block timesteps
  !-----------------------------------------------------------------------
  rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
  read(num_param,timesteps)

! dmr&aurel ##  namelist/timesteps/dtmin,dtmax,dtt,testdiag,tbegin,tend
  write(num_rep_42,'(A)')'!___________________________________________________________'  
  write(num_rep_42,'(A)') '&timesteps                 ! bloc time steps'
  write(num_rep_42,*)
  write(num_rep_42,*) 'tend      = ',tend
  write(num_rep_42,*) 'tbegin    = ',tbegin
  write(num_rep_42,*) 'dtmin     = ',dtmin
  write(num_rep_42,*) 'dtmax     = ',dtmax
  write(num_rep_42,*) 'dtt       = ',dtt
  write(num_rep_42,*) 'testdiag  = ',testdiag


  write(num_rep_42,*)'/'                           
  write(num_rep_42,'(A)')'! tous les temps en annees. tbegin et tend : debut et fin du run '
  write(num_rep_42,'(A)')'! pour equation masse, pas de temps mini -> dtmin, maxi -> dtmax'      
  write(num_rep_42,'(A)')'! dtt : pas de temps long'
  write(num_rep_42,'(A)')'! testdiag, pour gérer le pas de temps dynamique dt'
  write(num_rep_42,'(A)')'! ordres de grandeur (a moduler selon dx) : '
  write(num_rep_42,'(A)')'! 40 km dtmin=2.e-2, dtmax=1., dtt=5., tesdiag=0.02'

  write(num_rep_42,*)


  !      si shelf rapide
  !        testdiag=0.016

  !---------------------------------------------------------------------
  dt=dtmin           ! sera réattribue a chaque pas de temps
  ntmax=90000000     ! nombre de tours maxi dans la boucle temps
  dttest=dtmin       ! sert dans plein de tests 
  time=tbegin 
  nt=-1



  ! ----------------parametres dont le traitement en namelist est remis a plus tard


  ! Tous les parametres liés à l'ablation devraient être dans un module ablation
  ! A changer plus tard

  !     temperature forcing : annual, july. 
  TAFOR=0.
  TJFOR=0.0
  SEALEVEL=0.0
  sealevel_2d(:,:)=0.

  SECYEAR=365.*24.*3600.
  secyear= 31556926 ! s /an   pour Heino

!!$!     *** ABLATION
!!$!     integrating step for positive degree days (degrees)
!!$      DTP=2.0
!!$!     number of months in 1 year, st. dev. for temp *)
!!$      NYEAR=12
!!$      SIGMA=5.0
!!$!     proportion of melted water that can refreeze *)
!!$      CSI=0.6
!!$!     melting factors for snow and ice 
!!$      Csnow=0.003
!!$      Cice=0.008
!!$!     ct for PDD calculation
!!$      S22=0.5/SIGMA/SIGMA
!!$      PY=2*PI/NYEAR
!!$      PDDCT=DTP/SIGMA/sqrt(2.*PI)/NYEAR*365.


! calcul de certain parametres de param_phy car ro et g ne sont plus parameter

rog      = ro*g                !<  ro*g   (glace)
rowg     = row*g               !<  row*g  (ocean) 
romg     = rom*g               !<  rom*g  (asthenosphere)
rofreshg = rofresh*g           !<  rofresh*g 
d00      = ro**3*g**3          !<  ro**3*g**3
dice     = ro/row              !<  ro/row

coef_Sflot = (Row-Ro)/Row      !<   S = coef_Sflot * H + sealevel pour les shelves
coef_Bflot = -Ro/Row           !<   B = coef_Bflot * H + sealevel pour les shelves

  !--------------------- 


  !     *** TEMPERATURE IN ICE AND MANTLE
  !     Total number of vertical grid points, grid step in mantle
  NZZ=NZ+NZM

  !     specific latent heat of fusion of ice J/Kg
  CL=3.35E5

  !     DA mantle diffusion
  DA=4.E7
  !     switch 
  NICE=1
  isynchro=0


end subroutine initial_phy