Changeset 3605 for LMDZ6/branches/Ocean_skin/libf/phylmd/tracco2i_mod.F90

Timestamp:

Nov 21, 2019, 4:43:45 PM (5 years ago)

Author:

lguez

Message:

Merge revisions 3427:3600 of trunk into branch Ocean_skin

Location:

LMDZ6/branches/Ocean_skin

Files:

• . (modified) (1 prop)
• libf/phylmd/tracco2i_mod.F90 (modified) (14 diffs)

LMDZ6/branches/Ocean_skin/libf/phylmd/tracco2i_mod.F90

@@ -2 +2 @@
 !
 ! This module does the work for the interactive CO2 tracers
+! Authors: Patricia Cadule and Olivier Boucher
+!
+! Purpose and description:
+!  -----------------------
+! Main routine for the interactive carbon cycle
+! Gather all carbon fluxes and emissions from ORCHIDEE, PISCES and fossil fuel
+! Compute the net flux in source field which is used in phytrac
+! Compute global CO2 mixing ratio for radiation scheme if option is activated
+! Redistribute CO2 evenly over the atmosphere if transport is desactivated
 !
 CONTAINS
@@ -10 +19 @@
 
     USE dimphy
-    USE infotrac
-    USE geometry_mod, ONLY : cell_area
-    USE carbon_cycle_mod, ONLY : nbcf_in, fields_in, cfname_in, fco2_ocn_day, fco2_ff, fco2_bb
+    USE infotrac_phy
+    USE geometry_mod, ONLY: cell_area
+    USE carbon_cycle_mod, ONLY: carbon_cycle_init
+    USE carbon_cycle_mod, ONLY: id_CO2, nbcf_in, fields_in, cfname_in
+    USE carbon_cycle_mod, ONLY: fco2_ocn_day, fco2_ff, fco2_bb, fco2_land, fco2_ocean
+    USE carbon_cycle_mod, ONLY: fco2_land_nbp, fco2_land_nep, fco2_land_fLuc
+    USE carbon_cycle_mod, ONLY: fco2_land_fwoodharvest, fco2_land_fHarvest
+    USE carbon_cycle_mod, ONLY: carbon_cycle_cpl, carbon_cycle_tr, carbon_cycle_rad, RCO2_glo, RCO2_tot
     USE mod_grid_phy_lmdz
-    USE mod_phys_lmdz_mpi_data, ONLY :  is_mpi_root
+    USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root
     USE mod_phys_lmdz_para, ONLY: gather, bcast, scatter
     USE phys_cal_mod
+    USE phys_state_var_mod, ONLY: pctsrf
+    USE indice_sol_mod, ONLY: nbsrf, is_ter, is_lic, is_oce, is_sic
 
     IMPLICIT NONE
@@ -45 +61 @@
 !----------------
 
-    INTEGER, PARAMETER :: id_CO2=1              !--temporaire OB -- to be changed 
     INTEGER                        :: it, k, i, nb
     REAL, DIMENSION(klon,klev)     :: m_air     ! mass of air in every grid box [kg]
-    REAL, DIMENSION(klon)          :: co2land   ! surface land CO2 emissions [kg CO2/m2/s]
-    REAL, DIMENSION(klon)          :: co2ocean  ! surface ocean 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)
+    LOGICAL, SAVE :: check_fCO2_nbp_in_cfname
+!$OMP THREADPRIVATE(check_fCO2_nbp_in_cfname)
+    INTEGER, SAVE :: day_pre=-1
+!$OMP THREADPRIVATE(day_pre)
 
     IF (is_mpi_root) THEN
@@ -67 +79 @@
 !--convert 280 ppm into kg CO2 / kg air
     IF (debutphy) THEN
+
+! Initialisation de module carbon_cycle_mod
+      IF (carbon_cycle_cpl) THEN
+        CALL carbon_cycle_init()
+      ENDIF
+
+! Initialisation de tr_seri(id_CO2) si pas initialise
       IF (MAXVAL(tr_seri(:,:,id_CO2)).LT.1.e-15) THEN 
-        tr_seri(:,:,id_CO2)=280.e-6/RMD*RMCO2
+        tr_seri(:,:,id_CO2)=co2_ppm0*1.e-6/RMD*RMCO2 !--initialised from co2_ppm0 in rdem
       ENDIF 
+
+!--check if fCO2_nbp is in 
+      check_fCO2_nbp_in_cfname=.FALSE.
+      DO nb=1, nbcf_in
+        IF (cfname_in(nb)=="fCO2_nbp") check_fCO2_nbp_in_cfname=.TRUE.
+      ENDDO
+
     ENDIF
 
@@ -85 +111 @@
 
 !--retrieving land and ocean CO2 flux
-!--fCO2_nep comes in unit of g CO2 m-2 dt_stomate-1 
-!--this needs to be changed in ORCHIDEE
-    co2land(:)=0.0
-    co2ocean(:)=0.0
+    fco2_land(:)=0.0
+    fco2_ocean(:)=0.0
+    fco2_land_nbp(:)=0.
+    fco2_land_nep(:)=0.
+    fco2_land_fLuc(:)=0.
+    fco2_land_fwoodharvest(:)=0.
+    fco2_land_fHarvest(:)=0.
+
     DO nb=1, nbcf_in
-      IF (cfname_in(nb) == "fCO2_nep" )   co2land(:)=fields_in(:,nb)*RMCO2/RMC/86400./1000.
-      !!IF (cfname_in(nb) == "fCO2_fgco2" ) co2ocean(:)=fco2_ocn_day(:) !--for now
+
+      SELECT CASE(cfname_in(nb))
+!--dealing with the different fluxes coming from ORCHIDEE
+!--fluxes come in unit of kg C m-2 s-1 is converted into kg CO2 m-2 s-1
+      CASE("fCO2_nep")
+          fco2_land_nep(:)=fields_in(:,nb)*RMCO2/RMC*pctsrf(:,is_ter)
+      CASE("fCO2_fLuc")
+          fco2_land_fLuc(:)=fields_in(:,nb)*RMCO2/RMC*pctsrf(:,is_ter)
+      CASE("fCO2_fwoodharvest")
+          fco2_land_fwoodharvest(:)=fields_in(:,nb)*RMCO2/RMC*pctsrf(:,is_ter)
+      CASE("fCO2_fHarvest")
+          fco2_land_fHarvest(:)=fields_in(:,nb)*RMCO2/RMC*pctsrf(:,is_ter)
+      CASE("fCO2_nbp")
+          fco2_land_nbp(:)=fields_in(:,nb)*RMCO2/RMC*pctsrf(:,is_ter)
+!--fCO2_fco2_ocn comes in unit of mol C02 m-2 s-1 is converted into kg CO2 m-2 s-1 + change sign
+      CASE("fCO2_fgco2") 
+          fco2_ocean(:)=-1.*fco2_ocn_day(:)*RMCO2/1.e3*(pctsrf(:,is_oce)+pctsrf(:,is_sic))
+      END SELECT
+
     ENDDO
 
-!--preparing the net anthropogenic flux at the surface for mixing layer 
-!--unit kg CO2 / m2 / s
-    source(:,id_CO2)=fco2_ff(:)+fco2_bb(:)+co2land(:)+co2ocean(:)
+!--if fCO2_nbp is transferred we use it, otherwise we use the sum of what has been passed from ORCHIDEE
+    IF (check_fCO2_nbp_in_cfname)  THEN
+       fco2_land(:)=fco2_land_nbp(:)
+    ELSE
+       fco2_land(:)=fco2_land_nep(:)+fco2_land_fLuc(:)+fco2_land_fwoodharvest(:)+fco2_land_fHarvest(:)
+    ENDIF
+
+!!--preparing the net anthropogenic flux at the surface for mixing layer
+!!--unit kg CO2 / m2 / s
+!    PRINT *, 'tracco2i_mod.F90 --- MAXVAL(fco2_ff) ',MAXVAL(fco2_ff)
+!    PRINT *, 'tracco2i_mod.F90 --- MINVAL(fco2_ff) ',MINVAL(fco2_ff)
+!    PRINT *, 'tracco2i_mod.F90 --- MAXVAL(fco2_bb) ',MAXVAL(fco2_bb)
+!    PRINT *, 'tracco2i_mod.F90 --- MINVAL(fco2_bb) ',MINVAL(fco2_bb)
+!    PRINT *, 'tracco2i_mod.F90 --- MAXVAL(fco2_land) ',MAXVAL(fco2_land)
+!    PRINT *, 'tracco2i_mod.F90 --- MINVAL(fco2_land) ',MINVAL(fco2_land)
+!    PRINT *, 'tracco2i_mod.F90 --- MAXVAL(fco2_ocean) ',MAXVAL(fco2_ocean)
+!    PRINT *, 'tracco2i_mod.F90 --- MINVAL(fco2_ocean) ',MINVAL(fco2_ocean)
+!    PRINT *, 'tracco2i_mod.F90 --- MAXVAL(source(:,id_CO2)) ',MAXVAL(source(:,id_CO2))
+!    PRINT *, 'tracco2i_mod.F90 --- MINVAL(source(:,id_CO2)) ',MINVAL(source(:,id_CO2))
+!
+!--build final source term for CO2
+    source(:,id_CO2)=fco2_ff(:)+fco2_bb(:)+fco2_land(:)+fco2_ocean(:)
 
 !--computing global mean CO2 for radiation
-!--every timestep for now but enough every month
-!    IF (debutphy.OR.mth_cur.NE.mth_pre) THEN
+!--for every timestep comment out the IF ENDIF statements 
+!--otherwise this is updated every day
+    IF (debutphy.OR.day_cur.NE.day_pre) THEN
+
       CALL gather(tr_seri(:,:,id_CO2),co2_glo)
       CALL gather(m_air,m_air_glo)
+
 !$OMP MASTER
-!--conversion from kg CO2/kg air into ppm
+
+!--compute a global mean CO2 value and print its value in ppm
        IF (is_mpi_root) THEN
-         RCO2_glo=SUM(co2_glo*m_air_glo)/SUM(m_air_glo)*1.e6*RMD/RMCO2
+         RCO2_tot=SUM(co2_glo*m_air_glo)  !--unit kg CO2
+         RCO2_glo=RCO2_tot/SUM(m_air_glo) !--unit kg CO2 / kg air
+         PRINT *,'tracco2i: global CO2 in ppm =', RCO2_glo*1.e6*RMD/RMCO2
+         PRINT *,'tracco2i: total CO2 in kg =', RCO2_tot
        ENDIF
-       PRINT *,'toto in tracco2i: global CO2 in ppm =', RCO2_glo
 !$OMP END MASTER
        CALL bcast(RCO2_glo)
-       mth_pre=mth_cur
-!    ENDIF
+       day_pre=day_cur
+!--if not carbon_cycle_tr, then we reinitialize the CO2 each day to its global mean value
+       IF (.NOT.carbon_cycle_tr) THEN 
+         tr_seri(:,:,id_CO2)=RCO2_glo
+       ENDIF
+    ENDIF
 
   END SUBROUTINE tracco2i
@@ -119 +195 @@
 
     USE dimphy
-    USE infotrac
+    USE infotrac_phy
     USE geometry_mod, ONLY : cell_area
     USE mod_grid_phy_lmdz
@@ -129 +205 @@
     USE netcdf, ONLY: nf90_get_var, nf90_noerr, nf90_nowrite
 
-    USE carbon_cycle_mod, ONLY : fco2_ff, fco2_bb
+    USE carbon_cycle_mod, ONLY : fco2_ff, fco2_bb, fco2_land, fco2_ocean
 
     IMPLICIT NONE
@@ -150 +226 @@
 !! may be controlled via the .def later on
 !! also co2bb for now comes from ORCHIDEE
-    LOGICAL, PARAMETER :: readco2ff=.TRUE., readco2bb=.FALSE. 
+    LOGICAL, PARAMETER :: readco2ff=.TRUE.
+!! this should be left to FALSE for now
+    LOGICAL, PARAMETER :: readco2bb=.FALSE. 
+
+    CHARACTER (len = 20) :: modname = 'tracco2i.co2_emissions'
+    CHARACTER (len = 80) :: abort_message
 
     IF (debutphy) THEN 
@@ -173 +254 @@
         n_glo = size(vector)
         IF (n_glo.NE.klon_glo) THEN
-           PRINT *,'sflx_lmdz_co2_ff: le nombre de points n est pas egal a klon_glo'
-           STOP
+           abort_message='sflx_lmdz_co2_ff: le nombre de points n est pas egal a klon_glo'
+           CALL abort_physic(modname,abort_message,1)
         ENDIF
 
@@ -181 +262 @@
         n_month = size(time)
         IF (n_month.NE.12) THEN
-           PRINT *,'sflx_lmdz_co2_ff: le nombre de month n est pas egal a 12'
-           STOP
+           abort_message='sflx_lmdz_co2_ff: le nombre de month n est pas egal a 12'
+           CALL abort_physic(modname,abort_message,1)
         ENDIF
 
@@ -196 +277 @@
 
 !--reading CO2 biomass burning emissions
+!--using it will be inconsistent with treatment in ORCHIDEE
       IF (readco2bb) THEN 
 
@@ -205 +287 @@
       n_glo = size(vector)
       IF (n_glo.NE.klon_glo) THEN
-         PRINT *,'sflx_lmdz_co2_bb: le nombre de points n est pas egal a klon_glo'
-         STOP
+         abort_message='sflx_lmdz_co2_bb: le nombre de points n est pas egal a klon_glo'
+         CALL abort_physic(modname,abort_message,1)
       ENDIF
 
@@ -213 +295 @@
       n_month = size(time)
       IF (n_month.NE.12) THEN
-         PRINT *,'sflx_lmdz_co2_bb: le nombre de month n est pas egal a 12'
-         STOP
+         abort_message='sflx_lmdz_co2_bb: le nombre de month n est pas egal a 12'
+         CALL abort_physic(modname,abort_message,1)
       ENDIF
 
@@ -247 +329 @@
     PRINT *,'probleme avec le mois dans co2_ini =', mth_cur
   ENDIF
-  IF (.NOT.ALLOCATED(fco2_ff)) ALLOCATE(fco2_ff(klon))
-  IF (.NOT.ALLOCATED(fco2_bb)) ALLOCATE(fco2_bb(klon))
+
   fco2_ff(:) = flx_co2ff(:,mth_cur)
   fco2_bb(:) = flx_co2bb(:,mth_cur)