!!!! Abderrahmane Idelkadi aout 2013 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Module pour definir (au 1er appel) et ecrire les variables dans les fichiers de sortie cosp
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   MODULE cosp_output_write_mod
  
   USE cosp_output_mod
  
   INTEGER, SAVE  :: itau_iocosp
!$OMP THREADPRIVATE(itau_iocosp)

!  INTERFACE histwrite_cosp
!    MODULE PROCEDURE histwrite2d_cosp,histwrite3d_cosp
!  END INTERFACE

   CONTAINS

  SUBROUTINE cosp_output_write(Nlevlmdz, Npoints, Ncolumns, itap, dtime, freq_COSP, cfg, gbx, sglidar, stlidar, isccp)

    USE ioipsl
    USE control_mod

#ifdef CPP_XIOS
    ! ug Pour les sorties XIOS
    USE wxios
#endif

!!! Variables d'entree
  integer               :: itap, Nlevlmdz, Ncolumns, Npoints
  real                  :: freq_COSP, dtime
  type(cosp_config)     :: cfg     ! Control outputs
  type(cosp_gridbox)    :: gbx     ! Gridbox information. Input for COSP
  type(cosp_sglidar)    :: sglidar ! Output from lidar simulator
  type(cosp_isccp)      :: isccp   ! Output from ISCCP simulator
  type(cosp_lidarstats) :: stlidar ! Summary statistics from lidar simulator

!!! Variables locales
  integer               :: icl, iinitend=1
  logical               :: ok_sync
  integer               :: itau_wcosp
  real, dimension(Npoints,PARASOL_NREFL) :: parasolcrefl, Ncref

  include "temps.h"

  IF (MOD(itap,NINT(freq_COSP/dtime)).EQ.0) THEN
       
! A refaire
       itau_wcosp = itau_phy + itap + start_time * day_step / iphysiq

! On le donne a  cosp_output_write_mod pour que les histwrite y aient acces:
       CALL set_itau_iocosp(itau_wcosp)

    ok_sync = .TRUE.
    
    IF(.NOT.cosp_varsdefined) THEN
        iinitend = 2
    ELSE
        iinitend = 1
    ENDIF

! ug la boucle qui suit ne sert qu'une fois, pour l'initialisation, sinon il n'y a toujours qu'un seul passage:
DO iinit=1, iinitend
#ifdef CPP_XIOS
IF (cosp_varsdefined) THEN
     CALL wxios_update_calendar(itau_iocosp)
END IF
#endif

 if (cfg%Llidar_sim) then

! Pb des valeurs indefinies, on les met a 0
! A refaire proprement
  do k = 1,Nlevout
     do ip = 1,Npoints
     if(stlidar%lidarcld(ip,k).eq.R_UNDEF)then
      stlidar%lidarcld(ip,k)=0.
     endif
     enddo

     do ii= 1,SR_BINS
      do ip = 1,Npoints
       if(stlidar%cfad_sr(ip,ii,k).eq.R_UNDEF)then
        stlidar%cfad_sr(ip,ii,k)=0.
       endif
      enddo
     enddo
   enddo

  do ip = 1,Npoints
   do k = 1,Nlevlmdz
     if(sglidar%beta_mol(ip,k).eq.R_UNDEF)then
      sglidar%beta_mol(ip,k)=0.
     endif

     do ii= 1,Ncolumns
       if(sglidar%beta_tot(ip,ii,k).eq.R_UNDEF)then
        sglidar%beta_tot(ip,ii,k)=0.
       endif
     enddo

    enddo    !k = 1,Nlevlmdz
   enddo     !ip = 1,Npoints

   do k = 1,LIDAR_NCAT
    do ip = 1,Npoints
     if(stlidar%cldlayer(ip,k).eq.R_UNDEF)then
      stlidar%cldlayer(ip,k)=0.
     endif
    enddo
   enddo

   CALL histwrite2d_cosp(o_cllcalipso,stlidar%cldlayer(:,1))
   CALL histwrite2d_cosp(o_clhcalipso,stlidar%cldlayer(:,3))
   CALL histwrite2d_cosp(o_clmcalipso,stlidar%cldlayer(:,2)) 
   CALL histwrite2d_cosp(o_cltcalipso,stlidar%cldlayer(:,4))
   CALL histwrite3d_cosp(o_clcalipso,stlidar%lidarcld,nvert)

   do icl=1,SR_BINS
      CALL histwrite3d_cosp(o_cfad_lidarsr532,stlidar%cfad_sr(:,icl,:),nvert,icl)
   enddo

   CALL histwrite3d_cosp(o_parasol_refl,stlidar%parasolrefl,nvertp)

   do k=1,PARASOL_NREFL
    do ip=1, Npoints
     if (stlidar%cldlayer(ip,4).gt.0.01) then
       parasolcrefl(ip,k)=(stlidar%parasolrefl(ip,k)-0.03*(1.-stlidar%cldlayer(ip,4)))/ &
                            stlidar%cldlayer(ip,4)
        Ncref(ip,k) = 1.
     else
        parasolcrefl(ip,k)=0.
        Ncref(ip,k) = 0.
     endif
    enddo
   enddo
   CALL histwrite3d_cosp(o_parasol_crefl,parasolcrefl,nvertp)

   do icl=1,Ncolumns 
      CALL histwrite3d_cosp(o_atb532,sglidar%beta_tot(:,icl,:),nvertmcosp,icl)
   enddo
   CALL histwrite3d_cosp(o_beta_mol532,sglidar%beta_mol,nvertmcosp) 
 endif !Lidar

 if (cfg%Lisccp_sim) then

! Traitement des valeurs indefinies
   do ip = 1,Npoints
    if(isccp%totalcldarea(ip).eq.-1.E+30)then
      isccp%totalcldarea(ip)=0.
    endif
    if(isccp%meanptop(ip).eq.-1.E+30)then
      isccp%meanptop(ip)=0.
    endif
    if(isccp%meantaucld(ip).eq.-1.E+30)then
      isccp%meantaucld(ip)=0.
    endif
    if(isccp%meanalbedocld(ip).eq.-1.E+30)then
      isccp%meanalbedocld(ip)=0.
    endif
    if(isccp%meantb(ip).eq.-1.E+30)then
      isccp%meantb(ip)=0.
    endif
    if(isccp%meantbclr(ip).eq.-1.E+30)then
      isccp%meantbclr(ip)=0.
    endif

    do k=1,7
     do ii=1,7
     if(isccp%fq_isccp(ip,ii,k).eq.-1.E+30)then
      isccp%fq_isccp(ip,ii,k)=0.
     endif
     enddo
    enddo

    do ii=1,Ncolumns
     if(isccp%boxtau(ip,ii).eq.-1.E+30)then
       isccp%boxtau(ip,ii)=0.
     endif
    enddo

    do ii=1,Ncolumns
     if(isccp%boxptop(ip,ii).eq.-1.E+30)then
       isccp%boxptop(ip,ii)=0.
     endif
    enddo
   enddo

   CALL histwrite2d_cosp(o_sunlit,gbx%sunlit)
   do icl=1,7
   CALL histwrite3d_cosp(o_clisccp2,isccp%fq_isccp(:,icl,:),nvertisccp,icl) 
   enddo
   CALL histwrite3d_cosp(o_boxtauisccp,isccp%boxtau,nvertcol)
   CALL histwrite3d_cosp(o_boxptopisccp,isccp%boxptop,nvertcol) 
   CALL histwrite2d_cosp(o_tclisccp,isccp%totalcldarea) 
   CALL histwrite2d_cosp(o_ctpisccp,isccp%meanptop) 
   CALL histwrite2d_cosp(o_tauisccp,isccp%meantaucld) 
   CALL histwrite2d_cosp(o_albisccp,isccp%meanalbedocld) 
   CALL histwrite2d_cosp(o_meantbisccp,isccp%meantb) 
   CALL histwrite2d_cosp(o_meantbclrisccp,isccp%meantbclr)
 endif ! Isccp

 IF(.NOT.cosp_varsdefined) THEN
!$OMP MASTER
            DO iff=1,3
                IF (cosp_outfilekeys(iff)) THEN
                  CALL histend(cosp_nidfiles(iff))
                ENDIF ! cosp_outfilekeys
            ENDDO !  iff
#ifdef CPP_XIOS
            !On finalise l'initialisation:
            CALL wxios_closedef()
#endif
!$OMP END MASTER
!$OMP BARRIER
            cosp_varsdefined = .TRUE.
 END IF
END DO  !! iinit

!    IF(cosp_varsdefined) THEN
! On synchronise les fichiers pour IOIPSL
!$OMP MASTER
   DO iff=1,3
         IF (ok_sync .AND. cosp_outfilekeys(iff)) THEN
             CALL histsync(cosp_nidfiles(iff))
         ENDIF
   END DO
!$OMP END MASTER

   ENDIF ! if freq_COSP

    END SUBROUTINE cosp_output_write

! ug Routine pour definir itau_iocosp depuis cosp_output_write_mod:
  SUBROUTINE set_itau_iocosp(ito)
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: ito
      itau_iocosp = ito
  END SUBROUTINE

  SUBROUTINE histdef2d_cosp (iff,var)

    USE ioipsl
    USE dimphy
    use iophy
    USE mod_phys_lmdz_para
#ifdef CPP_XIOS
  USE wxios
#endif

    IMPLICIT NONE

    INCLUDE "dimensions.h"
    INCLUDE "temps.h"

    INTEGER                          :: iff
    TYPE(ctrl_outcosp)               :: var

    REAL zstophym
    CHARACTER(LEN=20) :: typeecrit

    ! ug On récupère le type écrit de la structure:
    !       Assez moche, Ã|  refaire si meilleure méthode...
    IF (INDEX(var%cosp_typeecrit(iff), "once") > 0) THEN
       typeecrit = 'once'
    ELSE IF(INDEX(var%cosp_typeecrit(iff), "t_min") > 0) THEN
       typeecrit = 't_min(X)'
    ELSE IF(INDEX(var%cosp_typeecrit(iff), "t_max") > 0) THEN
       typeecrit = 't_max(X)'
    ELSE IF(INDEX(var%cosp_typeecrit(iff), "inst") > 0) THEN
       typeecrit = 'inst(X)'
    ELSE
       typeecrit = cosp_outfiletypes(iff)
    ENDIF

    IF (typeecrit=='inst(X)'.OR.typeecrit=='once') THEN
       zstophym=zoutm_cosp(iff)
    ELSE
       zstophym=zdtimemoy_cosp
    ENDIF

#ifdef CPP_XIOS
        CALL wxios_add_field_to_file(var%name, 2, cosp_nidfiles(iff), cosp_outfilenames(iff), &
        var%description, var%unit, 1, typeecrit)
#endif
       IF ( var%cles(iff) ) THEN
          CALL histdef (cosp_nidfiles(iff), var%name, var%description, var%unit, &
               iim,jj_nb,nhoricosp(iff), 1,1,1, -99, 32, &
               typeecrit, zstophym,zoutm_cosp(iff))
       ENDIF

  END SUBROUTINE histdef2d_cosp

 SUBROUTINE histdef3d_cosp (iff,var,nvertsave,ncols)
    USE ioipsl
    USE dimphy
    use iophy
    USE mod_phys_lmdz_para

#ifdef CPP_XIOS
  USE wxios
#endif


    IMPLICIT NONE

    INCLUDE "dimensions.h"
    INCLUDE "temps.h"

    INTEGER                        :: iff, klevs
    INTEGER, INTENT(IN), OPTIONAL  :: ncols ! ug RUSTINE POUR LES variables 4D
    INTEGER, INTENT(IN)           :: nvertsave
    TYPE(ctrl_outcosp)             :: var

    REAL zstophym
    CHARACTER(LEN=20) :: typeecrit, nomi
    CHARACTER(LEN=20) :: nom
    character(len=2) :: str2

! Axe vertical
      IF (nvertsave.eq.nvertp(iff)) THEN
          klevs=PARASOL_NREFL
      ELSE IF (nvertsave.eq.nvertisccp(iff)) THEN
          klevs=7
      ELSE IF (nvertsave.eq.nvertcol(iff)) THEN
          klevs=Ncolout
      ELSE
          klevs=Nlevout
      ENDIF
      
         
! ug RUSTINE POUR LES Champs 4D
      IF (PRESENT(ncols)) THEN
               write(str2,'(i2.2)')ncols
               nomi=var%name
               nom="c"//str2//"_"//nomi
      ELSE
               nom=var%name
      END IF

    ! ug On récupère le type écrit de la structure:
    !       Assez moche, Ã|  refaire si meilleure méthode...
    IF (INDEX(var%cosp_typeecrit(iff), "once") > 0) THEN
       typeecrit = 'once'
    ELSE IF(INDEX(var%cosp_typeecrit(iff), "t_min") > 0) THEN
       typeecrit = 't_min(X)'
    ELSE IF(INDEX(var%cosp_typeecrit(iff), "t_max") > 0) THEN
       typeecrit = 't_max(X)'
    ELSE IF(INDEX(var%cosp_typeecrit(iff), "inst") > 0) THEN
       typeecrit = 'inst(X)'
    ELSE
       typeecrit = cosp_outfiletypes(iff)
    ENDIF

    IF (typeecrit=='inst(X)'.OR.typeecrit=='once') THEN
       zstophym=zoutm_cosp(iff)
    ELSE
       zstophym=zdtimemoy_cosp
    ENDIF

#ifdef CPP_XIOS
        CALL wxios_add_field_to_file(nom, 3, cosp_nidfiles(iff), cosp_outfilenames(iff), &
        var%description, var%unit, 1, typeecrit)
#endif
       IF ( var%cles(iff) ) THEN
          CALL histdef (cosp_nidfiles(iff), nom, var%description, var%unit, &
               iim, jj_nb, nhoricosp(iff), klevs, 1, &
               klevs, nvertsave, 32, typeecrit, &
               zstophym, zoutm_cosp(iff))
       ENDIF
  END SUBROUTINE histdef3d_cosp

 SUBROUTINE histwrite2d_cosp(var,field)
  USE dimphy
  USE mod_phys_lmdz_para
  USE ioipsl
  use iophy

#ifdef CPP_XIOS
  USE wxios
#endif

  IMPLICIT NONE
  INCLUDE 'dimensions.h'
  INCLUDE 'iniprint.h'

    TYPE(ctrl_outcosp), INTENT(IN) :: var
    REAL, DIMENSION(:), INTENT(IN) :: field

    INTEGER :: iff

    REAL,DIMENSION(klon_mpi) :: buffer_omp
    INTEGER, allocatable, DIMENSION(:) :: index2d
    REAL :: Field2d(iim,jj_nb)
    CHARACTER(LEN=20) ::  nomi, nom
    character(len=2) :: str2

    IF (prt_level >= 9) WRITE(lunout,*)'Begin histrwrite2d ',var%name

  ! On regarde si on est dans la phase de définition ou d'écriture:
  IF(.NOT.cosp_varsdefined) THEN
!$OMP MASTER
      !Si phase de définition.... on définit
      CALL conf_cospoutputs(var%name,var%cles)
      DO iff=1, 3
         IF (cosp_outfilekeys(iff)) THEN
            CALL histdef2d_cosp(iff, var)
         ENDIF
      ENDDO
!$OMP END MASTER
  ELSE
    !Et sinon on.... écrit
    IF (SIZE(field)/=klon) &
  CALL abort_gcm('iophy::histwrite2d_cosp','Field first DIMENSION not equal to klon',1) 

    CALL Gather_omp(field,buffer_omp)
!$OMP MASTER
    CALL grid1Dto2D_mpi(buffer_omp,Field2d)

! La boucle sur les fichiers:
      DO iff=1, 3
           IF (var%cles(iff) .AND. cosp_outfilekeys(iff)) THEN
                ALLOCATE(index2d(iim*jj_nb))
        CALL histwrite(cosp_nidfiles(iff),var%name,itau_iocosp,Field2d,iim*jj_nb,index2d) 
#ifdef CPP_XIOS
                IF (iff == 1) THEN
                   CALL wxios_write_2D(var%name, Field2d)
                ENDIF
#endif

                deallocate(index2d)
           ENDIF !levfiles
      ENDDO
!$OMP END MASTER   
  ENDIF ! vars_defined
  IF (prt_level >= 9) WRITE(lunout,*)'End histrwrite2d_cosp ',nom
  END SUBROUTINE histwrite2d_cosp

! ug NOUVELLE VERSION DES WRITE AVEC LA BOUCLE DO RENTREE
! AI sept 2013
  SUBROUTINE histwrite3d_cosp(var, field, nverts, ncols)
  USE dimphy
  USE mod_phys_lmdz_para
  USE ioipsl
  use iophy

#ifdef CPP_XIOS
 USE WXIOS
#endif


  IMPLICIT NONE
  INCLUDE 'dimensions.h'
  INCLUDE 'iniprint.h'

    TYPE(ctrl_outcosp), INTENT(IN)    :: var
    REAL, DIMENSION(:,:), INTENT(IN)  :: field ! --> field(klon,:)
    INTEGER, INTENT(IN), OPTIONAL     :: ncols ! ug RUSTINE POUR LES Champs 4D.....
    INTEGER, DIMENSION(3), INTENT(IN) :: nverts

    INTEGER :: iff, k

    REAL,DIMENSION(klon_mpi,SIZE(field,2)) :: buffer_omp
    REAL :: Field3d(iim,jj_nb,SIZE(field,2))
    INTEGER :: ip, n, nlev
    INTEGER, ALLOCATABLE, DIMENSION(:) :: index3d
    CHARACTER(LEN=20) ::  nomi, nom
    character(len=2) :: str2

  IF (prt_level >= 9) write(lunout,*)'Begin histrwrite3d ',var%name

! ug RUSTINE POUR LES STD LEVS.....
      IF (PRESENT(ncols)) THEN
              write(str2,'(i2.2)')ncols
              nomi=var%name
              nom="c"//str2//"_"//nomi
      ELSE
               nom=var%name
      END IF
  ! On regarde si on est dans la phase de définition ou d'écriture:
  IF(.NOT.cosp_varsdefined) THEN
      !Si phase de définition.... on définit
!$OMP MASTER
      CALL conf_cospoutputs(var%name,var%cles)
      DO iff=1, 3
        IF (cosp_outfilekeys(iff)) THEN
          CALL histdef3d_cosp(iff, var, nverts(iff), ncols)
        ENDIF
      ENDDO
!$OMP END MASTER
  ELSE
    !Et sinon on.... écrit
    IF (SIZE(field,1)/=klon) &
   CALL abort_gcm('iophy::histwrite3d','Field first DIMENSION not equal to klon',1)                                  
    nlev=SIZE(field,2)


    CALL Gather_omp(field,buffer_omp)
!$OMP MASTER
    CALL grid1Dto2D_mpi(buffer_omp,field3d)


! BOUCLE SUR LES FICHIERS
     DO iff=1, 3
        IF (var%cles(iff) .AND. cosp_outfilekeys(iff)) THEN
           ALLOCATE(index3d(iim*jj_nb*nlev))
    CALL histwrite(cosp_nidfiles(iff),nom,itau_iocosp,Field3d,iim*jj_nb*nlev,index3d) 

#ifdef CPP_XIOS
           IF (iff == 1) THEN
               CALL wxios_write_3D(nom, Field3d(:,:,1:klev))
           ENDIF
#endif
         deallocate(index3d)
        ENDIF
      ENDDO
!$OMP END MASTER   
  ENDIF ! vars_defined
  IF (prt_level >= 9) write(lunout,*)'End histrwrite3d_cosp ',nom
  END SUBROUTINE histwrite3d_cosp

  SUBROUTINE conf_cospoutputs(nam_var,cles_var)
!!! Lecture des noms et cles de sortie des variables dans config.def
    !   en utilisant les routines getin de IOIPSL  
    use ioipsl

    IMPLICIT NONE
    include 'iniprint.h'

   CHARACTER(LEN=20)               :: nam_var, nnam_var
   LOGICAL, DIMENSION(3)           :: cles_var

! Lecture dans config.def ou output.def de cles_var et name_var
    CALL getin('cles_'//nam_var,cles_var)
    CALL getin('name_'//nam_var,nam_var)
    IF(prt_level>10) WRITE(lunout,*)'nam_var cles_var ',nam_var,cles_var(:)

  END SUBROUTINE conf_cospoutputs

 END MODULE cosp_output_write_mod