!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Interface pour ecrire en netcdf avec les routines d'enseignement
! iotd de Frederic Hourdin
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      subroutine iophys_ecrit(nom,lllm,titre,unite,px)

      USE mod_phys_lmdz_para, ONLY: klon_omp, is_mpi_root
      USE mod_phys_lmdz_transfert_para, ONLY: gather
      USE mod_grid_phy_lmdz, ONLY: klon_glo, nbp_lon, nbp_lat, grid1dto2d_glo
      IMPLICIT NONE



!  Ecriture de variables diagnostiques au choix dans la physique 
!  dans un fichier NetCDF nomme  'diagfi'. Ces variables peuvent etre
!  3d (ex : temperature), 2d (ex : temperature de surface), ou
!  0d (pour un scalaire qui ne depend que du temps : ex : la longitude
!  solaire)
!  (ou encore 1d, dans le cas de testphys1d, pour sortir une colonne)
!  La periode d'ecriture est donnee par 
!  "ecritphy " regle dans le fichier de controle de run :  run.def
!
!    writediagfi peut etre appele de n'importe quelle subroutine
!    de la physique, plusieurs fois. L'initialisation et la creation du
!    fichier se fait au tout premier appel.
!
! WARNING : les variables dynamique (u,v,t,q,ps)
!  sauvees par writediagfi avec une
! date donnee sont legerement differentes que dans le fichier histoire car 
! on ne leur a pas encore ajoute de la dissipation et de la physique !!!
! IL est  RECOMMANDE d'ajouter les tendance physique a ces variables
! avant l'ecriture dans diagfi (cf. physiq.F)
!  
! Modifs: Aug.2010 Ehouarn: enforce outputs to be real*4
!
!  parametres (input) :
!  ----------
!      unit : unite logique du fichier de sortie (toujours la meme)
!      nom  : nom de la variable a sortir (chaine de caracteres)
!      titre: titre de la variable (chaine de caracteres)
!      unite : unite de la variable (chaine de caracteres)
!      px : variable a sortir (real 0, 1, 2, ou 3d)
!
!=================================================================


! Arguments on input:
      integer lllm
      character (len=*) :: nom,titre,unite
      integer imjmax
      parameter (imjmax=100000)
      real px(klon_omp,lllm)
      real xglo(klon_glo,lllm)
      real zx(nbp_lon,nbp_lat,lllm)



      CALL Gather(px,xglo)
!$OMP MASTER
      IF (is_mpi_root) THEN       
        CALL Grid1Dto2D_glo(xglo,zx)
        call iotd_ecrit(nom,lllm,titre,unite,zx)
      ENDIF
!$OMP END MASTER

      return
      end subroutine iophys_ecrit

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Version avec reindexation pour appeler depuis les routines internes
! à la sous surface
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


    subroutine iophys_ecrit_index(nom,lllm,titre,unite,knon,knindex,px)

    USE mod_phys_lmdz_para, ONLY: klon_omp
    USE dimphy, ONLY : klon
    USE mod_grid_phy_lmdz, ONLY: klon_glo
    IMPLICIT NONE

! This subroutine returns the sea surface temperature already read from limit.nc
!

! Arguments on input:
    INTEGER lllm
    CHARACTER (len=*) :: nom,titre,unite
    REAL px(klon_omp,lllm)
    INTEGER, INTENT(IN)                  :: knon     ! nomber of points on compressed grid
    INTEGER, DIMENSION(klon), INTENT(IN) :: knindex  ! grid point number for compressed grid
    REAL, DIMENSION(klon,lllm) :: varout

    INTEGER :: i,l

    IF (klon/=klon_omp) THEN
      print*,'klon, klon_omp',klon,klon_omp
      CALL abort_physic('iophys_ecrit','probleme de dimension parallele',1)
    ENDIF

    varout(1:klon,1:lllm)=0.
    DO l = 1, lllm
    DO i = 1, knon
       varout(knindex(i),l) = px(i,l)
    END DO
    END DO
    CALL iophys_ecrit(nom,lllm,titre,unite,varout)

  END SUBROUTINE iophys_ecrit_index

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      SUBROUTINE iophys_ini(timestep,nlev)
      USE dimphy, ONLY: klev
      USE mod_phys_lmdz_para, ONLY: is_mpi_root
      USE vertical_layers_mod, ONLY: presnivs
      USE regular_lonlat_mod, ONLY: lon_reg, lat_reg
      USE mod_grid_phy_lmdz, ONLY: klon_glo
      USE time_phylmdz_mod, ONLY : annee_ref, day_ref, day_ini
      USE phys_cal_mod, ONLY : calend

      USE yomcst_mod_h
IMPLICIT NONE


!=======================================================================
!
!   Auteur:  L. Fairhead  ,  P. Le Van, Y. Wanherdrick, F. Forget
!   -------
!
!   Objet:
!   ------
!
!   'Initialize' the diagfi.nc file: write down dimensions as well
!   as time-independent fields (e.g: geopotential, mesh area, ...)
!
!=======================================================================
!-----------------------------------------------------------------------
!   Declarations:
!   -------------

integer, intent(in) :: nlev
real, intent(in) :: timestep

real pi
INTEGER nlat_eff
INTEGER jour0,mois0,an0
REAL t0
CHARACTER(len=20) :: calendrier
integer ilev
real coord_vert(nlev)

!   Arguments:
!   ----------


!$OMP MASTER
    IF (is_mpi_root) THEN       

! Bidouille pour gerer le fait que lat_reg contient deux latitudes
! en version uni-dimensionnelle (chose qui pourrait être résolue
! par ailleurs)
IF (klon_glo==1) THEN
   nlat_eff=1
ELSE
   nlat_eff=size(lat_reg)
ENDIF
pi=2.*asin(1.)

! print*,'day_ini,annee_ref,day_ref',day_ini,annee_ref,day_ref
! print*,'jD_ref,jH_ref,start_time, calend',jD_ref,jH_ref,start_time, calend

! Attention : les lignes ci dessous supposent un calendrier en 360 jours
! Pourrait être retravaillé

jour0=day_ref-30*(day_ref/30)
mois0=day_ref/30+1
an0=annee_ref
!FH BIZARE QUAND 1D ...  t0=(day_ini-1)*RDAY
t0=0.
calendrier=calend
if ( calendrier == "earth_360d" ) calendrier="360_day"

print*,'iophys_ini annee_ref day_ref',annee_ref,day_ref,day_ini,calend,t0

if ( nlev == klev ) then
     coord_vert=presnivs
print*,'ON EST LA '
else
     do ilev=1,nlev
        coord_vert(ilev)=ilev
     enddo
endif
print*,'nlev=',nlev
print*,'coord_vert',coord_vert
call iotd_ini('phys.nc', &
size(lon_reg),nlat_eff,nlev,lon_reg(:)*180./pi,lat_reg*180./pi,coord_vert,jour0,mois0,an0,t0,timestep,calendrier)
 !  SUBROUTINE iotd_ini(fichnom,iim,jjm,llm,prlon,prlat,pcoordv,jour0,mois0,an0,t0,timestep,calendrier)
!   -------
    ENDIF
!$OMP END MASTER

      END SUBROUTINE iophys_ini

#ifdef und
      SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
      IMPLICIT none

!=======================================================================
      INTEGER nfield,nlon,iim,jjmp1, jjm
      REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)

      INTEGER i, n, ig

      jjm = jjmp1 - 1
      DO n = 1, nfield
         DO i=1,iim
            ecrit(i,n) = fi(1,n)
            ecrit(i+jjm*iim,n) = fi(nlon,n)
         ENDDO
         DO ig = 1, nlon - 2
           ecrit(iim+ig,n) = fi(1+ig,n)
         ENDDO
      ENDDO
      RETURN
      END SUBROUTINE gr_fi_ecrit

#endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Interface pour ecrire en netcdf avec les routines d'enseignement
! iotd de Frederic Hourdin
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      SUBROUTINE iotd_ecrit_seq(nom,lllm,titre,unite,px)
!call iotd_ecrit_seq('f0',1,'f0 in thermcell_plume_6A',' ',f0(1:ngrid))

        USE iotd_mod_h

      IMPLICIT NONE


! Arguments on input:
      integer lllm
      character (len=*) :: nom,titre,unite
      integer imjmax
      parameter (imjmax=100000)
      real px(imjmax*lllm)
      real, allocatable :: zx(:,:,:)
      integer i,j,l,ijl

      !print*,'iotd_ecrit_seq ,nom,lllm,titre,unite,px',nom,lllm,titre,unite,px
      allocate(zx(imax,jmax,lllm))

      ijl=0
      do l=1,lllm
         ! Pole nord
         ijl=ijl+1
         do i=1,imax
            zx(i,1,l)=px(ijl)
         enddo
         ! Grille normale
         do j=2,jmax-1
            do i=1,imax
               ijl=ijl+1
               zx(i,j,l)=px(ijl)
            enddo
         enddo
         ! Pole sud
         if ( jmax > 1 ) then
            ijl=ijl+1
            do i=1,imax
               zx(i,jmax,l)=px(ijl)
            enddo
         endif
      enddo

      call iotd_ecrit(nom,lllm,titre,unite,zx)
      deallocate(zx)

      return
      end subroutine iotd_ecrit_seq

