! $Id: stratH2O_methox.F90 3677 2020-05-06 15:18:32Z oboucher $

SUBROUTINE stratH2O_methox(debutphy,paprs,dq_ch4mmr)

! output: CH4VMR in MMR/s (mass mixing ratio/s or kg H2O/kg air/s)

  USE netcdf95, ONLY: nf95_close, nf95_gw_var, nf95_inq_dimid, & 
                      nf95_inq_varid, nf95_inquire_dimension, nf95_open
  USE netcdf, ONLY: nf90_get_var, nf90_noerr, nf90_nowrite

  USE lmdz_grid_phy
  USE lmdz_phys_mpi_data, ONLY:  is_mpi_root
  USE lmdz_phys_omp_data, ONLY:  is_omp_root

  USE lmdz_phys_para
  USE dimphy
  USE phys_cal_mod, ONLY: mth_cur
  USE infotrac_phy
  USE aerophys
  USE lmdz_yomcst
  USE strataer_local_var_mod, ONLY: flag_newclim_file
  
  IMPLICIT NONE

  include "dimensions.h"
  
! Input variables
  REAL paprs(klon,klev+1)
  LOGICAL, INTENT(IN) :: debutphy   ! flag for first physiq step
! Output variables
! tendency buffer used in add_phys_tend SUBROUTINE (in physiq_mod)
  REAL, INTENT(INOUT), DIMENSION(klon,klev)  :: dq_ch4mmr
  
! Local variables
  INTEGER n_lat   ! number of latitudes in the input data
  INTEGER n_lon   ! number of longitudes in the input data
  INTEGER, SAVE :: n_lev   ! number of levels in the input data
!$OMP THREADPRIVATE(n_lev)
  INTEGER n_mth   ! number of months in the input data
  INTEGER, SAVE :: mth_pre
!$OMP THREADPRIVATE(mth_pre)
  
! Reconstitutes fields
  REAL paprs_glo(klon_glo,klev+1)
  
  REAL, ALLOCATABLE:: latitude(:)
! (of input data sorted in strictly ascending order)
  REAL, ALLOCATABLE:: longitude(:)
! (of input data sorted in strictly ascending order)
  REAL, ALLOCATABLE:: time(:)
! (of input data sorted in strictly ascending order)
  REAL, ALLOCATABLE:: lev(:)
! levels of input data
  
! Stratospheric H2O source from CH4 oxidation from fixed climatos
! (H2O production in VMR/sec)
  REAL, ALLOCATABLE :: CH4RVMR_in(:, :, :, :)
  REAL, ALLOCATABLE :: CH4RVMR_mth(:, :, :)
  REAL, ALLOCATABLE :: CH4RVMR_tmp(:, :)
  REAL CH4RVMR_glo(klon_glo,klev)

  INTEGER i, k, kk, j
  
! For NetCDF:
  INTEGER ncid_in  ! IDs for input files
  INTEGER varid, ncerr
    
  INTEGER, PARAMETER :: lev_input=17
!--pressure at interfaces of input data (in Pa)
  REAL, DIMENSION(lev_input+1), PARAMETER ::          & 
                    paprs_input=(/                    &
  1.00000002e+05,   6.06530673e+04,   3.67879449e+04, &
  2.23130165e+04,   1.35335286e+04,   8.20850004e+03, &
  4.97870695e+03,   3.01973841e+03,   1.83156393e+03, &
  1.11089968e+03,   6.73794715e+02,   4.08677153e+02, &
  2.47875223e+02,   1.50343923e+02,   9.11881985e+01, &
  5.53084382e+01,   3.35462635e+01,   0.0           /)

  IF (debutphy .OR. mth_cur /= mth_pre) THEN
     
!--preparation of global fields
     CALL gather(paprs, paprs_glo)
     
     IF (is_mpi_root.AND.is_omp_root) THEN
        
        print *,'In stratH2O_methox read file: ch4r_annual_lmdz.nc'
        
        !--reading strat. H2O source (CH4 oxidation) files
        CALL nf95_open("ch4r_annual_lmdz.nc", nf90_nowrite, ncid_in)
        
        CALL  nf95_inq_varid(ncid_in, "LEV", varid)
        CALL nf95_gw_var(ncid_in, varid, lev)
        n_lev = size(lev)
        
        CALL nf95_inq_varid(ncid_in, "LAT", varid)
        CALL nf95_gw_var(ncid_in, varid, latitude)
        n_lat = size(latitude)
        
        CALL nf95_inq_varid(ncid_in, "LON", varid)
        CALL nf95_gw_var(ncid_in, varid, longitude)
        n_lon = size(longitude)
        
        CALL nf95_inq_varid(ncid_in, "TIME", varid)
        CALL nf95_gw_var(ncid_in, varid, time)
        n_mth = size(time)
        
        IF (.NOT.ALLOCATED(CH4RVMR_in)) ALLOCATE(CH4RVMR_in(n_lon, n_lat, n_lev, n_mth))
        
        CALL nf95_inq_varid(ncid_in, "CH4RVMR", varid)
        ncerr = nf90_get_var(ncid_in, varid, CH4RVMR_in)
        print *,'code erreur CH4RVMR=', ncerr, varid
        
        CALL nf95_close(ncid_in)
        
        IF (.NOT.ALLOCATED(CH4RVMR_mth)) ALLOCATE(CH4RVMR_mth(n_lon, n_lat, n_lev))
        IF (.NOT.ALLOCATED(CH4RVMR_tmp)) ALLOCATE(CH4RVMR_tmp(klon_glo, n_lev))
        
        !---select the correct month,
        !---correct latitudinal order,convert input from volume mixing ratio to mass mixing ratio
        DO j=1,n_lat
           ! convert VMR/s in MMR/s (mass mixing ratio/s or kg H2O/kg air/s)
           ! x2 because CH4->2*H2O
           CH4RVMR_mth(:,j,:) = 2*CH4RVMR_in(:,n_lat+1-j,:,mth_cur)*mH2Omol/mAIRmol
        ENDDO
        
        !---reduce to a klon_glo grid but keep the levels
        CALL grid2dTo1d_glo(CH4RVMR_mth,CH4RVMR_tmp)
        
        !---regrid weighted climatologies
        DO i=1, klon_glo
           DO k=1, klev
              
              CH4RVMR_glo(i,k)=0.0
              
              DO kk=1, n_lev

                 CH4RVMR_glo(i,k)=CH4RVMR_glo(i,k)+ &
                      MAX(0.0,MIN(paprs_glo(i,k),paprs_input(kk)) &
                      -MAX(paprs_glo(i,k+1),paprs_input(kk+1))) &
                      *CH4RVMR_tmp(i,kk)/(paprs_glo(i,k)-paprs_glo(i,k+1))
              ENDDO ! kk loop
           ENDDO ! k loop
        ENDDO ! i loop
     ENDIF !--is_mpi_root and is_omp_root
     
     !--keep memory of previous month
     mth_pre=mth_cur
     
     !--scatter global fields around
     CALL scatter(CH4RVMR_glo, dq_ch4mmr)
     
     IF (is_mpi_root.AND.is_omp_root) THEN
        DEALLOCATE(CH4RVMR_in,CH4RVMR_mth,CH4RVMR_tmp)
     ENDIF !--is_mpi_root and is_omp_root
  ENDIF ! debutphy.OR.new month
  

  
END SUBROUTINE stratH2O_methox