source: LMDZ6/trunk/libf/phylmd/tracco2i_mod.F90 @ 3361

MODULE tracco2i_mod
!
! This module does the work for the interactive CO2 tracers
!
CONTAINS

  SUBROUTINE tracco2i(pdtphys, debutphy, &
       xlat, xlon, pphis, pphi, &
       t_seri, pplay, paprs, tr_seri, source)

    USE dimphy
    USE infotrac
    USE geometry_mod, ONLY : cell_area
    USE mod_grid_phy_lmdz
    USE mod_phys_lmdz_mpi_data, ONLY :  is_mpi_root
    USE mod_phys_lmdz_para, ONLY: gather, bcast, scatter
    USE phys_cal_mod

    IMPLICIT NONE

    INCLUDE "clesphys.h"
    INCLUDE "YOMCST.h"

! Input argument
!---------------
    REAL,INTENT(IN)    :: pdtphys    ! Pas d'integration pour la physique (seconde)
    LOGICAL,INTENT(IN) :: debutphy   ! le flag de l'initialisation de la physique

    REAL,DIMENSION(klon),INTENT(IN)        :: xlat    ! latitudes pour chaque point 
    REAL,DIMENSION(klon),INTENT(IN)        :: xlon    ! longitudes pour chaque point 
    REAL,DIMENSION(klon),INTENT(IN)        :: pphis   ! geopotentiel du sol
    REAL,DIMENSION(klon,klev),INTENT(IN)   :: pphi    ! geopotentiel de chaque couche

    REAL,DIMENSION(klon,klev),INTENT(IN)   :: t_seri  ! Temperature
    REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay   ! pression pour le mileu de chaque couche (en Pa)
    REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs   ! pression pour chaque inter-couche (en Pa)
    REAL,DIMENSION(klon,nbtr),INTENT(INOUT):: source  ! flux de traceur [U/m2/s]

! Output argument
!----------------
    REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT)  :: tr_seri ! Concentration Traceur [U/kgA]  

! Local variables
!----------------

    INTEGER, PARAMETER :: id_CO2=1                           !--temporaire OB=> PC to be changed 
    REAL, PARAMETER    :: MCO2=44.011   !--g/mol
    INTEGER                                 :: it, k, i
    REAL, DIMENSION(klon,klev)              :: m_air          ! mass of air in every grid box [kg]
    REAL, DIMENSION(klon)                   :: co2ff          ! surface fossil-fuel CO2 emissions [kg CO2/m2/s]
    REAL, DIMENSION(klon)                   :: co2bb          ! surface biomass burning CO2 emissions [kg CO2/m2/s]
    REAL, DIMENSION(klon_glo,klev) :: co2_glo   ! variable temporaire sur la grille global
    REAL, DIMENSION(klon_glo,klev) :: m_air_glo ! variable temporaire sur la grille global


    INTEGER, SAVE :: mth_pre=0
!$OMP THREADPRIVATE(mth_pre)
    REAL, SAVE :: RCO2_glo
!$OMP THREADPRIVATE(RCO2_glo)

    IF (is_mpi_root) THEN
      PRINT *,'in tracco2i: date from phys_cal_mod =',year_cur,'-',mth_cur,'-',day_cur,'-',hour
    ENDIF

!--initialisation of CO2 field if not in restart file
!--dirty way of doing, do it better later
!--convert 280 ppm into kg CO2 / kg air
    IF (debutphy) THEN
      IF (MAXVAL(tr_seri(:,:,id_CO2)).LT.1.e-15) THEN 
        tr_seri(:,:,id_CO2)=280.e-6/RMD*MCO2
      ENDIF 
    ENDIF

!--calculate mass of air in every grid box in kg 
    DO i=1, klon
    DO k=1, klev
      m_air(i,k)=(paprs(i,k)-paprs(i,k+1))/RG*cell_area(i)
    ENDDO
    ENDDO

!--call CO2 emission routine 
!--co2bb is zero for now
    CALL co2_emissions(debutphy,co2ff,co2bb)

!--preparing the net anthropogenic flux at the surface for mixing layer 
!--unit kg CO2 / m2 / s
    DO i=1, klon
      source(i,id_CO2)=co2ff(i)+co2bb(i)
    ENDDO

!--computing global mean CO2 for radiation
!    IF (debutphy.OR.mth_cur.NE.mth_pre) THEN
      CALL gather(tr_seri(:,:,id_CO2),co2_glo)
      CALL gather(m_air,m_air_glo)
!$OMP MASTER
       IF (is_mpi_root) THEN
         RCO2_glo=SUM(co2_glo*m_air_glo)/SUM(m_air_glo)*1.e6*RMD/MCO2
       ENDIF
       PRINT *,'in tracco2i: global CO2 in ppm =', RCO2_glo
!$OMP END MASTER
       CALL bcast(RCO2_glo)
       mth_pre=mth_cur
!    ENDIF

  END SUBROUTINE tracco2i

  SUBROUTINE co2_emissions(debutphy,co2ff,co2bb)

    USE dimphy
    USE infotrac
    USE geometry_mod, ONLY : cell_area
    USE mod_grid_phy_lmdz
    USE mod_phys_lmdz_mpi_data, ONLY :  is_mpi_root
    USE mod_phys_lmdz_para, ONLY: gather, scatter
    USE phys_cal_mod
    USE YOMCST

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

    IMPLICIT NONE

    LOGICAL,INTENT(IN) :: debutphy
    REAL,DIMENSION(klon),INTENT(out) :: co2ff            !  fossil-fuel CO2 emissions
    REAL,DIMENSION(klon),INTENT(out) :: co2bb            !  biomass burning CO2 emissions

! For NetCDF:
    integer ncid_in  ! IDs for input files
    integer varid, ncerr

    INTEGER :: n_glo, n_month
    REAL, POINTER:: vector(:), time(:)
    REAL,ALLOCATABLE       :: flx_co2_glo(:,:) !  fossil-fuel CO2 
    REAL,ALLOCATABLE, SAVE :: flx_co2(:,:)     !  fossil-fuel CO2 
!$OMP THREADPRIVATE(flx_co2)

  IF (debutphy) THEN 

  ALLOCATE(flx_co2(klon,12))

!$OMP MASTER
  IF (is_mpi_root) THEN
    ! ... Open the file
    CALL nf95_open("sflx_lmdz_co2.nc", nf90_nowrite, ncid_in)

    CALL nf95_inq_varid(ncid_in, "vector", varid)
    CALL nf95_gw_var(ncid_in, varid, vector)
    n_glo = size(vector)
    IF (n_glo.NE.klon_glo) THEN
       print *,'Le nombre de points n est pas egal a klon_glo'
       STOP
    ENDIF

    CALL nf95_inq_varid(ncid_in, "time", varid)
    CALL nf95_gw_var(ncid_in, varid, time)
    n_month = size(time)
    IF (n_month.NE.12) THEN
       print *,'Le nombre de month n est pas egal a 12'
       STOP
    ENDIF

    IF (.NOT.ALLOCATED(flx_co2_glo)) ALLOCATE(flx_co2_glo(n_glo,n_month))

!--reading flx_co2
    CALL nf95_inq_varid(ncid_in, "flx_co2", varid)
    ncerr = nf90_get_var(ncid_in, varid, flx_co2_glo)

    CALL nf95_close(ncid_in)
  ENDIF
!$OMP END MASTER

!--scatter on all proc
  CALL scatter(flx_co2_glo,flx_co2)

!$OMP MASTER
  IF (is_mpi_root) THEN
     DEALLOCATE(flx_co2_glo)
  ENDIF
!$OMP END MASTER

  ENDIF !--end debuthy

!---select the correct month
  IF (mth_cur.LT.1.OR.mth_cur.GT.12) THEN
    print *,'probleme avec le mois dans co2_ini =', mth_cur
  ENDIF
  co2ff(:) = flx_co2(:,mth_cur)
  co2bb(:) = 0.0

  END SUBROUTINE co2_emissions

END MODULE tracco2i_mod