source: LMDZ6/branches/Amaury_dev/libf/phylmd/lmdz_iotd.f90 @ 5218

!=======================================================================

!   Auteur:  F. Hourdin
!   -------

!   Objet:
!   ------
!   Light interface for netcdf outputs. can be used outside LMDZ

!=======================================================================

MODULE lmdz_iotd
  IMPLICIT NONE; PRIVATE
  PUBLIC iotd_fin, iotd_ecrit, iotd_ini, imax, jmax

  INTEGER imax, jmax, lmax, nid
  INTEGER dim_coord(4)
  REAL iotd_ts, iotd_t0
  INTEGER :: n_names_iotd_def
  CHARACTER*20, DIMENSION(200) :: names_iotd_def
  CHARACTER*20 :: un_nom

  !$OMP THREADPRIVATE(imax, jmax, lmax, nid, dim_coord, iotd_t0, iotd_ts)
  !$OMP THREADPRIVATE(n_names_iotd_def, names_iotd_def)
CONTAINS
  SUBROUTINE iotd_fin
    USE netcdf, ONLY: nf90_close
    IMPLICIT NONE
    INTEGER ierr

    ierr = nf90_close(nid)
  END SUBROUTINE iotd_fin

  SUBROUTINE iotd_ecrit(nom, llm, titre, unite, px)
    !-----------------------------------------------------------------------
    !  ----------
    !      nom  : nom de la variable a sortir (chaine de caracteres)
    !      llm  : nombre de couches
    !      titre: titre de la variable (chaine de caracteres)
    !      unite : unite de la variable (chaine de caracteres)
    !      px : variable a sortir
    !=================================================================

    USE netcdf, ONLY: nf90_put_var, nf90_inq_varid, nf90_enddef, nf90_redef, nf90_sync, nf90_noerr, &
            nf90_float, nf90_def_var
    IMPLICIT NONE

    ! Arguments on input:
    INTEGER llm
    CHARACTER (LEN = *) :: nom, titre, unite
    INTEGER imjmax
    parameter (imjmax = 100000)
    REAL px(imjmax * llm)

    ! Local variables:

    real(kind = 4) date
    real(kind = 4) zx(imjmax * llm)

    INTEGER ierr, ndim, dim_cc(4)
    INTEGER iq
    INTEGER i, j, l

    INTEGER zitau
    CHARACTER firstnom*20
    SAVE firstnom
    SAVE zitau
    SAVE date
    DATA firstnom /'1234567890'/
    DATA zitau /0/

    ! Ajouts
    INTEGER, save :: ntime = 0
    INTEGER :: idim, varid
    CHARACTER (LEN = 50) :: fichnom
    INTEGER, DIMENSION(4) :: id
    INTEGER, DIMENSION(4) :: edges, corner

    IF (n_names_iotd_def>0 .and..not.any(names_iotd_def==nom)) RETURN
    !***************************************************************
    ! Initialisation of 'firstnom' and create/open the "diagfi.nc" NetCDF file
    ! ------------------------------------------------------------------------
    ! (Au tout premier appel de la SUBROUTINE durant le run.)


    !--------------------------------------------------------
    ! Write the variables to output file if it's time to do so
    !--------------------------------------------------------


    ! Compute/write/extend 'time' coordinate (date given in days)
    ! (done every "first call" (at given time level) to writediagfi)
    ! Note: date is incremented as 1 step ahead of physics time
    !--------------------------------------------------------

    zx(1:imax * jmax * llm) = px(1:imax * jmax * llm)
    IF (firstnom =='1234567890') THEN
      firstnom = nom
    endif

    !PRINT*,'nom ',nom,firstnom

    !! Quand on tombe sur la premiere variable on ajoute un pas de temps
    IF (nom==firstnom) THEN
      ! We have identified a "first call" (at given date)

      ntime = ntime + 1 ! increment # of stored time steps

      !!          PRINT*,'ntime ',ntime
      date = iotd_t0 + ntime * iotd_ts
      !PRINT*,'iotd_ecrit ',iotd_ts,ntime, date
      !          date= float (zitau +1)/float (day_step)

      ! compute corresponding date (in days and fractions thereof)
      ! Get NetCDF ID of 'time' variable

      ierr = nf90_sync(nid)

      ierr = nf90_inq_varid(nid, "time", varid)
      ! Write (append) the new date to the 'time' array

      ierr = nf90_put_var(nid, varid, date, [ntime])

      !          PRINT*,'date ',date,ierr,nid
      !        PRINT*,'IOTD Date ,varid,nid,ntime,date',varid,nid,ntime,date

      IF (ierr/=nf90_noerr) THEN
        WRITE(*, *) "***** PUT_VAR matter in writediagfi_nc"
        WRITE(*, *) "***** with time"
        WRITE(*, *) 'ierr=', ierr
      endif

      !          WRITE(6,*)'WRITEDIAGFI: date= ', date
    end if ! of if (nom.EQ.firstnom)


    !Case of a 3D variable
    !---------------------
    IF (llm==lmax) THEN
      ndim = 4
      corner(1) = 1
      corner(2) = 1
      corner(3) = 1
      corner(4) = ntime
      edges(1) = imax
      edges(2) = jmax
      edges(3) = llm
      edges(4) = 1
      dim_cc = dim_coord


      !Case of a 2D variable
      !---------------------

    ELSE IF (llm==1) THEN
      ndim = 3
      corner(1) = 1
      corner(2) = 1
      corner(3) = ntime
      corner(4) = 1
      edges(1) = imax
      edges(2) = jmax
      edges(3) = 1
      edges(4) = 1
      dim_cc(1:2) = dim_coord(1:2)
      dim_cc(3) = dim_coord(4)

    END IF ! of if llm=1 ou llm

    ! AU premier pas de temps, on crée les variables
    !-----------------------------------------------

    IF (ntime==1) THEN
      ierr = nf90_redef (nid)
      ierr = nf90_def_var(nid, nom, nf90_float, dim_cc, varid)
      !PRINT*,'DEF ',nom,nid,varid
      ierr = nf90_enddef(nid)
    ELSE
      ierr = nf90_inq_varid(nid, nom, varid)
      !PRINT*,'INQ ',nom,nid,varid
      ! Commandes pour recuperer automatiquement les coordonnees
      !             ierr= nf90_inq_dimid(nid,"longitude",id(1))
    END IF

    ierr = nf90_put_var(nid, varid, zx, corner, edges)

    IF (ierr/=nf90_noerr) THEN
      WRITE(*, *) "***** PUT_VAR problem in writediagfi"
      WRITE(*, *) "***** with ", nom
      WRITE(*, *) 'ierr=', ierr
    endif

  END

  SUBROUTINE iotd_ini(fichnom, iim, jjm, llm, prlon, prlat, pcoordv, jour0, mois0, an0, t0, timestep, calendrier)
    USE netcdf, ONLY: nf90_enddef, nf90_put_att, nf90_float, nf90_def_var, nf90_redef, &
            nf90_global, nf90_def_dim, nf90_create, nf90_clobber, nf90_unlimited, nf90_put_var
    IMPLICIT NONE

    INTEGER iim, jjm, llm
    REAL prlon(iim), prlat(jjm), pcoordv(llm), timestep, t0
    INTEGER id_FOCE
    INTEGER jour0, mois0, an0
    CHARACTER*(*) calendrier

    INTEGER corner(4), edges(4), ndim
    real  px(1000)
    CHARACTER (LEN = 10) :: nom
    real(kind = 4) rlon(iim), rlat(jjm), coordv(llm)

    !   Local:
    !   ------
    CHARACTER*3, DIMENSION(12) :: cmois = (/'JAN', 'FEB', 'MAR', 'APR', 'MAY', 'JUN', 'JUL', 'AUG', 'SEP', 'OCT', 'NOV', 'DEC'/)
    CHARACTER*10 date0
    CHARACTER*11 date0b

    INTEGER :: ierr

    INTEGER :: nvarid
    INTEGER, DIMENSION(2) :: id

    CHARACTER*(*) fichnom

    REAL pi

    iotd_ts = timestep
    iotd_t0 = t0
    PRINT*, 'iotd_ini, ', timestep, iotd_ts
    imax = iim
    jmax = jjm
    lmax = llm
    ! Utile pour passer en real*4 pour les ecritures
    rlon = prlon
    rlat = prlat
    coordv = pcoordv


    !-----------------------------------------------------------------------
    ! Possibilité de spécifier une liste de variables à sortir
    ! dans iotd.def
    ! Si iotd.def existe et est non vide,
    ! seules les variables faisant à la fois l'objet d'un CALL iotd_ecrit
    ! et étant spécifiées dans iotd.def sont sorties.
    ! Sinon, toutes les variables faisant l'objet d'un CALL iotd_ecrit
    ! sont sorties
    !-----------------------------------------------------------------------
    n_names_iotd_def = 0
    open(99, file = 'iotd.def', form = 'formatted', status = 'old', iostat = ierr)
    IF (ierr==0) THEN
      ierr = 0
      DO while (ierr==0)
        read(99, *, iostat = ierr) un_nom
        IF (ierr==0) THEN
          n_names_iotd_def = n_names_iotd_def + 1
          names_iotd_def(n_names_iotd_def) = un_nom
        endif
      enddo
    endif
    PRINT*, n_names_iotd_def, names_iotd_def(1:n_names_iotd_def)
    close(99)

    pi = 2. * asin(1.)

    ! Define dimensions

    ! Create the NetCDF file
    ierr = nf90_create(fichnom, nf90_clobber, nid)
    ierr = nf90_def_dim(nid, "lon", iim, dim_coord(1))
    ierr = nf90_def_dim(nid, "lat", jjm, dim_coord(2))
    ierr = nf90_def_dim(nid, "lev", llm, dim_coord(3))
    ierr = nf90_def_dim(nid, "time", nf90_unlimited, dim_coord(4))
    ierr = nf90_put_att(nid, nf90_global, 'Conventions', "CF-1.1")
    !ierr = nf90_put_att(nid,nf90_global,'file_name',TRIM(fname))
    ierr = nf90_enddef(nid)

    ! Switch out of NetCDF Define mode

    ierr = nf90_enddef(nid)

    !  Contol parameters for this run
    ! ---- longitude -----------

    ierr = nf90_redef(nid)
    ierr = nf90_def_var(nid, "lon", nf90_float, dim_coord(1), nvarid)
    ierr = nf90_put_att(nid, nvarid, 'axis', 'X')
    ierr = nf90_put_att(nid, nvarid, 'units', "degrees_east")
    ierr = nf90_enddef(nid)
    ierr = nf90_put_var(nid, nvarid, rlon)
    PRINT*, ierr

    ! ---- latitude ------------
    ierr = nf90_redef(nid)
    ierr = nf90_def_var(nid, "lat", nf90_float, dim_coord(2), nvarid)
    ierr = nf90_put_att(nid, nvarid, 'axis', 'Y')
    ierr = nf90_put_att(nid, nvarid, 'units', "degrees_north")
    ierr = nf90_enddef(nid)
    ierr = nf90_put_var(nid, nvarid, rlat)

    ! ---- vertical ------------
    ierr = nf90_redef(nid)
    ierr = nf90_def_var(nid, "lev", nf90_float, dim_coord(3), nvarid)
    ierr = nf90_put_att(nid, nvarid, "long_name", "vert level")
    IF (coordv(2)>coordv(1)) THEN
      ierr = nf90_put_att(nid, nvarid, "long_name", "pseudo-alt")
      ierr = nf90_put_att(nid, nvarid, 'positive', "up")
    else
      ierr = nf90_put_att(nid, nvarid, "long_name", "pressure")
      ierr = nf90_put_att(nid, nvarid, 'positive', "down")
    endif
    ierr = nf90_enddef(nid)
    ierr = nf90_put_var(nid, nvarid, coordv)

    ! ---- time ----------------
    ierr = nf90_redef(nid)
    ! Define the 'time' variable
    ierr = nf90_def_var(nid, "time", nf90_float, dim_coord(4), nvarid)
    !     ! Add attributes
    ierr = nf90_put_att(nid, nvarid, 'axis', 'T')
    ierr = nf90_put_att(nid, nvarid, 'standard_name', 'time')
    WRITE(date0, '(i4.4,"-",i2.2,"-",i2.2)') an0, mois0, jour0
    ierr = nf90_put_att(nid, nvarid, 'units', &
            "seconds since " // date0 // " 00:00:00")
    ierr = nf90_put_att(nid, nvarid, 'calendar', calendrier)
    !ierr = nf90_put_att(nid,nvarid,'calendar','360d')
    ierr = nf90_put_att(nid, nvarid, 'title', 'Time')
    ierr = nf90_put_att(nid, nvarid, 'long_name', 'Time axis')
    WRITE(date0b, '(i4.4,"-",a3,"-",i2.2)') an0, cmois(mois0), jour0
    ierr = nf90_put_att(nid, nvarid, 'time_origin', &
            date0b // ' 00:00:00')
    ierr = nf90_enddef(nid)

  END

END MODULE lmdz_iotd