!### macros ##################################################### #define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call goErr #define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if #define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if #include "tm5.inc" !################################################################# module emission_common ! While the averagers could be put into each flux module, managing ! these routines is easier when all the repetitive code is gathered ! here in one place. use GO, only: gol, goErr, goPr, goLabel use Dims, only: nlon360,nlat180 public :: co2_ocean, co2_fires, co2_ff, co2_bio, co2_all, co2_gpp, field2d_eps public :: co2_ocean_m, co2_fires_m, co2_ff_m, co2_bio_m, co2_all_m, co2_gpp_m, field2d_eps_m public :: co2_ocean_n, co2_fires_n, co2_ff_n, co2_bio_n, co2_all_n, co2_gpp_n, field2d_eps_n public :: c13_bio, c13_sig,c13_gpp, c13_res, c13_ocean, c13_ff, c13_fires, c13_biodis, c13_oceandis, c13_all public :: c13_bio_m, c13_gpp_m, c13_res_m, c13_ocean_m, c13_ff_m, c13_fires_m, c13_biodis_m, c13_oceandis_m, c13_all_m public :: c13_bio_n, c13_gpp_n, c13_res_n, c13_ocean_n, c13_ff_n, c13_fires_n, c13_biodis_n, c13_oceandis_n, c13_all_n public :: declare_emission_averagers,free_emission_averagers,reset_emission_averagers, calc_emission_averagers, save_emission_averagers public :: enkf_is_inverse, enkf_is_forward, enkf_set_inverse, enkf_set_forward public :: flux_means_are_priors, EmisInputDir,C13InputDir public :: ReadFromFile,ReadFromFile3D public :: fluxmultiplier, fluxmultiplier_biodis,fluxmultiplier_epsbio,GetParameterMaps, ParamInputFile character(len=*), parameter :: mname = 'emission_common' private:: mname real,dimension(nlon360,nlat180),target :: co2_ocean real,dimension(nlon360,nlat180),target :: co2_fires real,dimension(nlon360,nlat180),target :: co2_bio real,dimension(nlon360,nlat180),target :: co2_ff real,dimension(nlon360,nlat180),target :: co2_all real,dimension(nlon360,nlat180),target :: co2_gpp real,dimension(nlon360,nlat180),target :: c13_all real,dimension(nlon360,nlat180),target :: c13_bio real,dimension(nlon360,nlat180),target :: c13_gpp real,dimension(nlon360,nlat180),target :: c13_res real,dimension(nlon360,nlat180),target :: c13_ocean real,dimension(nlon360,nlat180),target :: c13_ff real,dimension(nlon360,nlat180),target :: c13_fires real,dimension(nlon360,nlat180),target :: c13_biodis real,dimension(nlon360,nlat180),target :: c13_oceandis real,dimension(nlon360,nlat180),target :: field2d_eps real,dimension(nlon360,nlat180),target :: c13_sig ! weekly average flux accumulators. "_m" is the flux field itself, ! "_n" is the number of accumulations. real,allocatable,dimension(:,:) :: co2_ocean_m real,allocatable,dimension(:,:) :: co2_fires_m real,allocatable,dimension(:,:) :: co2_bio_m real,allocatable,dimension(:,:) :: co2_ff_m real,allocatable,dimension(:,:) :: co2_all_m real,allocatable,dimension(:,:) :: co2_gpp_m real,allocatable,dimension(:,:) :: c13_all_m real,allocatable,dimension(:,:) :: c13_bio_m real,allocatable,dimension(:,:) :: c13_gpp_m real,allocatable,dimension(:,:) :: c13_res_m real,allocatable,dimension(:,:) :: c13_ocean_m real,allocatable,dimension(:,:) :: c13_ff_m real,allocatable,dimension(:,:) :: c13_fires_m real,allocatable,dimension(:,:) :: c13_biodis_m real,allocatable,dimension(:,:) :: c13_oceandis_m real,allocatable,dimension(:,:) :: field2d_eps_m !WP! This will hold multiplier values for each member of the inverse run, or !WP! for the mean of the forward run real,allocatable,dimension(:,:,:) :: fluxmultiplier real,allocatable,dimension(:,:,:) :: fluxmultiplier_biodis real,allocatable,dimension(:,:,:) :: fluxmultiplier_epsbio integer :: co2_ocean_n integer :: co2_fires_n integer :: co2_bio_n integer :: co2_ff_n integer :: co2_all_n integer :: co2_gpp_n integer :: c13_all_n integer :: c13_bio_n integer :: c13_gpp_n integer :: c13_res_n integer :: c13_ocean_n integer :: c13_ff_n integer :: c13_fires_n integer :: c13_biodis_n integer :: c13_oceandis_n integer :: field2d_eps_n ! Are we running in inverse or forward mode? Set enkf mode in Emis_Init(). integer,parameter :: enkf_mode_unknown = 0 integer,parameter :: enkf_mode_inverse = 1 integer,parameter :: enkf_mode_forward = 2 integer :: enkf_mode = enkf_mode_unknown character(len=256) :: EmisInputDir character(len=256) :: ParamInputFile character(len=256) :: C13InputDir character(len=1024) :: FluxOutFile = '' ! Weekly mean fluxes can either be means of the raw priors, or means ! of the optimized fluxes. Default is priors. This was the ! behavior in previous versions, since one needs to verify what the ! flux modules are doing, and it is possible to retrieve the scaled ! fluxes in postprocessing. This switch removes the need for that ! postprocessing, and may be useful in a future with a more complex ! parameter model. ! ! flux_means_are_priors = .TRUE. -> report weekly means of prior fluxes ! ! flux_means_are_priors = .FALSE. -> report weekly means of scaled/parameterized/optimized fluxes ! logical :: flux_means_are_priors = .TRUE. contains logical function enkf_is_inverse() enkf_is_inverse = .FALSE. if(enkf_mode .EQ. enkf_mode_inverse) then enkf_is_inverse = .TRUE. endif end function enkf_is_inverse logical function enkf_is_forward() enkf_is_forward = .FALSE. if(enkf_mode .EQ. enkf_mode_forward) then enkf_is_forward = .TRUE. endif end function enkf_is_forward ! Comment on flux_means_are_priors (Andy Jacobson, 29 Sep 2007) ! ! Wouter and I decided that the inverse run will output weekly-average ! prior fluxes by default, and that a follow-on forward run (which is ! already pretty much a requirement for various output) will write out ! the posterior fluxes. Thus the default values for flux_means_as_priors ! below. Note that in the inv and fwd versions of emission.F90, there ! is a commented-out ReadRC call that could be used to change this ! hard-coded default behavior. Use of old RC files with flux_means_are_priors ! lines may cause some confusion, but the variables written to the ! output file actually are given explicit names so it will be tough ! to make a real mistake. ! ! It is convenient to look at posterior fluxes while an inverse run is ! underway. For this, we have two suggestions. First, the current model ! for fluxes is that the posterior is a scaled version of the prior. The ! scale factors are the EnKF parameters, and are written to the savestate ! file along with the prior fluxes in the inverse run. Thus, the weekly ! mean posterior fluxes can be constructed from those outputs. This will ! not work if the parameter model changes such that we are optimizing ! parameters which affect the flux in a nonlinear fashion. Second, ! the forward run is generally much, much cheaper to run and should be ! able to generate "real" posterior fluxes quite quickly. subroutine enkf_set_inverse enkf_mode=enkf_mode_inverse flux_means_are_priors=.TRUE. ! print *,'EnKF mode set to inverse' end subroutine enkf_set_inverse subroutine enkf_set_forward enkf_mode=enkf_mode_forward ! print *,'EnKF mode set to forward' flux_means_are_priors=.TRUE. end subroutine enkf_set_forward subroutine declare_emission_averagers ! print *,"DEBUG: in declare_emission_averagers" allocate(co2_bio_m(nlon360,nlat180)) allocate(co2_fires_m(nlon360,nlat180)) allocate(co2_ocean_m(nlon360,nlat180)) allocate(co2_ff_m(nlon360,nlat180)) allocate(co2_all_m(nlon360,nlat180)) allocate(co2_gpp_m(nlon360,nlat180)) allocate(c13_all_m(nlon360,nlat180)) allocate(c13_bio_m(nlon360,nlat180)) allocate(c13_gpp_m(nlon360,nlat180)) allocate(c13_res_m(nlon360,nlat180)) allocate(c13_ocean_m(nlon360,nlat180)) allocate(c13_ff_m(nlon360,nlat180)) allocate(c13_fires_m(nlon360,nlat180)) allocate(c13_biodis_m(nlon360,nlat180)) allocate(c13_oceandis_m(nlon360,nlat180)) allocate(field2d_eps_m(nlon360,nlat180)) call reset_emission_averagers end subroutine declare_emission_averagers subroutine free_emission_averagers if(allocated(co2_all_m)) then deallocate(co2_all_m) endif if(allocated(co2_ocean_m)) then deallocate(co2_ocean_m) endif if(allocated(co2_fires_m)) then deallocate(co2_fires_m) endif if(allocated(co2_bio_m)) then deallocate(co2_bio_m) endif if(allocated(field2d_eps_m)) then deallocate(field2d_eps_m) endif if(allocated(co2_ff_m)) then deallocate(co2_ff_m) endif if(allocated(co2_gpp_m)) then deallocate(co2_gpp_m) endif if(allocated(c13_all_m)) then deallocate(c13_all_m) endif if(allocated(c13_bio_m)) then deallocate(c13_bio_m) endif if(allocated(c13_gpp_m)) then deallocate(c13_gpp_m) endif if(allocated(c13_res_m)) then deallocate(c13_res_m) endif if(allocated(c13_ocean_m)) then deallocate(c13_ocean_m) endif if(allocated(c13_ff_m)) then deallocate(c13_ff_m) endif if(allocated(c13_fires_m)) then deallocate(c13_fires_m) endif if(allocated(c13_biodis_m)) then deallocate(c13_biodis_m) endif if(allocated(c13_oceandis_m)) then deallocate(c13_oceandis_m) endif end subroutine free_emission_averagers subroutine calc_emission_averagers if(allocated(co2_all_m)) then if(co2_all_n .GT. 0) then co2_all_m=co2_all_m/co2_all_n endif endif if(allocated(co2_ff_m)) then if(co2_ff_n .GT. 0) then co2_ff_m=co2_ff_m/co2_ff_n endif endif if(allocated(co2_bio_m)) then if(co2_bio_n .GT. 0) then co2_bio_m=co2_bio_m/co2_bio_n endif endif if(allocated(field2d_eps_m)) then if(field2d_eps_n .GT. 0) then field2d_eps_m=field2d_eps_m/field2d_eps_n endif endif if(allocated(co2_fires_m)) then if(co2_fires_n .GT. 0) then co2_fires_m=co2_fires_m/co2_fires_n endif endif if(allocated(co2_ocean_m)) then if(co2_ocean_n .GT. 0) then co2_ocean_m=co2_ocean_m/co2_ocean_n endif endif if(allocated(co2_gpp_m)) then if(co2_gpp_n .GT. 0) then co2_gpp_m=co2_gpp_m/co2_gpp_n endif endif if(allocated(c13_all_m)) then if(c13_all_n .GT. 0) then c13_all_m=c13_all_m/c13_all_n endif endif if(allocated(c13_bio_m)) then if(c13_bio_n .GT. 0) then c13_bio_m=c13_bio_m/c13_bio_n endif endif if(allocated(c13_gpp_m)) then if(c13_gpp_n .GT. 0) then c13_gpp_m=c13_gpp_m/c13_gpp_n endif endif if(allocated(c13_res_m)) then if(c13_res_n .GT. 0) then c13_res_m=c13_res_m/c13_res_n endif endif if(allocated(c13_ocean_m)) then if(c13_ocean_n .GT. 0) then c13_ocean_m=c13_ocean_m/c13_ocean_n endif endif if(allocated(c13_ff_m)) then if(c13_ff_n .GT. 0) then c13_ff_m=c13_ff_m/c13_ff_n endif endif if(allocated(c13_fires_m)) then if(c13_fires_n .GT. 0) then c13_fires_m=c13_fires_m/c13_fires_n endif endif if(allocated(c13_biodis_m)) then if(c13_biodis_n .GT. 0) then c13_biodis_m=c13_biodis_m/c13_biodis_n endif endif if(allocated(c13_oceandis_m)) then if(c13_oceandis_n .GT. 0) then c13_oceandis_m=c13_oceandis_m/c13_oceandis_n endif endif end subroutine calc_emission_averagers subroutine save_emission_averagers(status) use ParTools, only: myid, root use MDF , only : MDF_Open use MDF , only : MDF_Create, MDF_Close, MDF_EndDef use MDF , only : MDF_NETCDF, MDF_REPLACE, MDF_GLOBAL, MDF_CHAR, MDF_INT, MDF_FLOAT, MDF_UNLIMITED, MDF_READ use MDF , only : MDF_Put_Att use MDF , only : MDF_Def_Dim use MDF , only : MDF_Inq_Dimid, MDF_Inquire_Dimension use MDF , only : MDF_Inq_VarID use MDF , only : MDF_Def_Var, MDF_Put_Var, MDF_Get_Var use global_data, only : outdir use dims, only : itaui, itaue use datetime, only : tau2date integer,intent(inout) :: status character(len=200) :: dsname,attname,comp,pp character(len=*), parameter :: rname = mname//'/save_emission_averagers' ! local integer :: i integer :: hid integer :: dim_lon_id, dim_lat_id integer :: var_flux_c13_all, var_flux_c13_bio, var_flux_c13_gpp, var_flux_co2_gpp, var_flux_c13_res, var_flux_c13_ocean, var_flux_c13_ff, var_flux_c13_fires, var_flux_c13_biodis, var_flux_c13_oceandis, var_flux_all,var_flux_ff, var_flux_bio, var_eps,var_flux_ocean, var_flux_fires, var_lon, var_lat, var_multiplier,var_multiplier_epsbio integer, dimension(6) :: idatei,idatee,idatef character(len=1024) :: attstring if(myid .ne. root) return pp='post' if(flux_means_are_priors .eqv. .True.) then pp='prior' endif call tau2date(itaui,idatei) call tau2date(itaue,idatee) write (FluxOutFile,'(a,"/flux1x1_",i4.4,3i2.2,"_",i4.4,3i2.2,".nc")') & trim(outdir), idatei(1:4), idatee(1:4) call MDF_Create( trim(FluxOutFile), MDF_NETCDF, MDF_REPLACE, hid, status ) IF_NOTOK_RETURN(status=1) ! Global Attributes write(attstring,'(i4.4,"/",i2.2,"/",i2.2," ",i2.2,":",i2.2,":",i2.2, " UTC")') idatei call MDF_Put_Att(hid, MDF_GLOBAL,"model_start_date",values=trim(attstring),status=status) IF_NOTOK_RETURN(status=1) write(attstring,'(i4.4,"/",i2.2,"/",i2.2," ",i2.2,":",i2.2,":",i2.2, " UTC")') idatee call MDF_Put_Att(hid, MDF_GLOBAL,"model_end_date",values=trim(attstring),status=status) IF_NOTOK_RETURN(status=1) ! define dimensions: call MDF_Def_Dim( hid, 'lon', 360 , dim_lon_id, status ) IF_NOTOK_RETURN(status=1) call MDF_Def_Dim( hid, 'lat', 180 , dim_lat_id, status ) IF_NOTOK_RETURN(status=1) ! dimension variables: call MDF_Def_Var( hid, 'lon', MDF_FLOAT, (/dim_lon_id/), var_lon, status ) IF_NOTOK_RETURN(status=1) call MDF_Put_Att(hid, var_lon,"units",values="degrees East",status=status) IF_NOTOK_RETURN(status=1) call MDF_Put_Att(hid, var_lon,"comment",values="center gridbox",status=status) IF_NOTOK_RETURN(status=1) call MDF_Def_Var( hid, 'lat', MDF_FLOAT, (/dim_lat_id/), var_lat, status ) IF_NOTOK_RETURN(status=1) call MDF_Put_Att(hid, var_lat,"units",values="degrees North",status=status) IF_NOTOK_RETURN(status=1) call MDF_Put_Att(hid, var_lat,"comment",values="center gridbox",status=status) IF_NOTOK_RETURN(status=1) ! define variables: dsname='flux_multiplier_m' call MDF_Def_Var( hid, trim(dsname), MDF_FLOAT, (/dim_lon_id,dim_lat_id/), var_multiplier, status ) IF_NOTOK_RETURN(status=1) call MDF_Put_Att(hid, var_multiplier,"units",values="dimensionless",status=status) IF_NOTOK_RETURN(status=1) dsname='flux_multiplier_epsbio_m' call MDF_Def_Var( hid, trim(dsname), MDF_FLOAT, (/dim_lon_id,dim_lat_id/), var_multiplier_epsbio, status ) IF_NOTOK_RETURN(status=1) call MDF_Put_Att(hid, var_multiplier_epsbio,"units",values="dimensionless",status=status) IF_NOTOK_RETURN(status=1) if(allocated(co2_all_m)) then comp='all' dsname='flux_'//trim(comp)//'_'//trim(pp)//'_mean' attname='flux_'//trim(comp)//'_'//trim(pp)//'_n' call MDF_Def_Var( hid, trim(dsname), MDF_FLOAT, (/dim_lon_id,dim_lat_id/), var_flux_all, status ) IF_NOTOK_RETURN(status=1) call MDF_Put_Att(hid, var_flux_all,"units",values="mol m-2 s-1",status=status) IF_NOTOK_RETURN(status=1) call MDF_Put_Att(hid, var_flux_all,attname,values=co2_all_n,status=status) IF_NOTOK_RETURN(status=1) endif if(allocated(co2_ff_m)) then comp='ff' dsname='flux_'//trim(comp)//'_'//trim(pp)//'_mean' attname='flux_'//trim(comp)//'_'//trim(pp)//'_n' call MDF_Def_Var( hid, trim(dsname), MDF_FLOAT, (/dim_lon_id,dim_lat_id/), var_flux_ff, status ) IF_NOTOK_RETURN(status=1) call MDF_Put_Att(hid, var_flux_ff,"units",values="mol m-2 s-1",status=status) IF_NOTOK_RETURN(status=1) call MDF_Put_Att(hid, var_flux_ff,attname,values=co2_ff_n,status=status) IF_NOTOK_RETURN(status=1) endif if(allocated(co2_bio_m)) then comp='bio' dsname='flux_'//trim(comp)//'_'//trim(pp)//'_mean' attname='flux_'//trim(comp)//'_'//trim(pp)//'_n' call MDF_Def_Var( hid, trim(dsname), MDF_FLOAT, (/dim_lon_id,dim_lat_id/), var_flux_bio, status ) IF_NOTOK_RETURN(status=1) call MDF_Put_Att(hid, var_flux_bio,"units",values="mol m-2 s-1",status=status) IF_NOTOK_RETURN(status=1) call MDF_Put_Att(hid, var_flux_bio,attname,values=co2_bio_n,status=status) IF_NOTOK_RETURN(status=1) endif if(allocated(field2d_eps_m)) then comp='eps' dsname='flux_'//trim(comp)//'_'//trim(pp)//'_mean' attname='flux_'//trim(comp)//'_'//trim(pp)//'_n' call MDF_Def_Var( hid, trim(dsname), MDF_FLOAT, (/dim_lon_id,dim_lat_id/), var_eps, status ) IF_NOTOK_RETURN(status=1) call MDF_Put_Att(hid, var_eps,"units",values="permil",status=status) IF_NOTOK_RETURN(status=1) call MDF_Put_Att(hid, var_eps,attname,values=field2d_eps_n,status=status) IF_NOTOK_RETURN(status=1) endif if(allocated(co2_fires_m)) then comp='fires' dsname='flux_'//trim(comp)//'_'//trim(pp)//'_mean' attname='flux_'//trim(comp)//'_'//trim(pp)//'_n' call MDF_Def_Var( hid, trim(dsname), MDF_FLOAT, (/dim_lon_id,dim_lat_id/), var_flux_fires, status ) IF_NOTOK_RETURN(status=1) call MDF_Put_Att(hid, var_flux_fires,"units",values="mol m-2 s-1",status=status) IF_NOTOK_RETURN(status=1) call MDF_Put_Att(hid, var_flux_fires,attname,values=co2_fires_n,status=status) IF_NOTOK_RETURN(status=1) endif if(allocated(co2_ocean_m)) then comp='ocean' dsname='flux_'//trim(comp)//'_'//trim(pp)//'_mean' attname='flux_'//trim(comp)//'_'//trim(pp)//'_n' call MDF_Def_Var( hid, trim(dsname), MDF_FLOAT, (/dim_lon_id,dim_lat_id/), var_flux_ocean, status ) IF_NOTOK_RETURN(status=1) call MDF_Put_Att(hid, var_flux_ocean,"units",values="mol m-2 s-1",status=status) IF_NOTOK_RETURN(status=1) call MDF_Put_Att(hid, var_flux_ocean,attname,values=co2_ocean_n,status=status) IF_NOTOK_RETURN(status=1) endif if(allocated(c13_all_m)) then comp='c13_all' dsname='flux_'//trim(comp)//'_'//trim(pp)//'_mean' attname='flux_'//trim(comp)//'_'//trim(pp)//'_n' call MDF_Def_Var( hid, trim(dsname), MDF_FLOAT, (/dim_lon_id,dim_lat_id/), var_flux_c13_all, status ) IF_NOTOK_RETURN(status=1) call MDF_Put_Att(hid, var_flux_c13_all,"units",values="mol m-2 s-1",status=status) IF_NOTOK_RETURN(status=1) call MDF_Put_Att(hid, var_flux_c13_all,attname,values=c13_all_n,status=status) IF_NOTOK_RETURN(status=1) endif if(allocated(c13_bio_m)) then comp='c13_bio' dsname='flux_'//trim(comp)//'_'//trim(pp)//'_mean' attname='flux_'//trim(comp)//'_'//trim(pp)//'_n' call MDF_Def_Var( hid, trim(dsname), MDF_FLOAT, (/dim_lon_id,dim_lat_id/), var_flux_c13_bio, status ) IF_NOTOK_RETURN(status=1) call MDF_Put_Att(hid, var_flux_c13_bio,"units",values="mol m-2 s-1",status=status) IF_NOTOK_RETURN(status=1) call MDF_Put_Att(hid, var_flux_c13_bio,attname,values=c13_bio_n,status=status) IF_NOTOK_RETURN(status=1) endif if(allocated(c13_gpp_m)) then comp='c13_gpp' dsname='flux_'//trim(comp)//'_'//trim(pp)//'_mean' attname='flux_'//trim(comp)//'_'//trim(pp)//'_n' call MDF_Def_Var( hid, trim(dsname), MDF_FLOAT, (/dim_lon_id,dim_lat_id/), var_flux_c13_gpp, status ) IF_NOTOK_RETURN(status=1) call MDF_Put_Att(hid, var_flux_c13_gpp,"units",values="mol m-2 s-1",status=status) IF_NOTOK_RETURN(status=1) call MDF_Put_Att(hid, var_flux_c13_gpp,attname,values=c13_gpp_n,status=status) IF_NOTOK_RETURN(status=1) endif if(allocated(co2_gpp_m)) then comp='gpp' dsname='flux_'//trim(comp)//'_'//trim(pp)//'_mean' attname='flux_'//trim(comp)//'_'//trim(pp)//'_n' call MDF_Def_Var( hid, trim(dsname), MDF_FLOAT, (/dim_lon_id,dim_lat_id/), var_flux_co2_gpp, status ) IF_NOTOK_RETURN(status=1) call MDF_Put_Att(hid, var_flux_co2_gpp,"units",values="mol m-2 s-1",status=status) IF_NOTOK_RETURN(status=1) call MDF_Put_Att(hid, var_flux_co2_gpp,attname,values=co2_gpp_n,status=status) IF_NOTOK_RETURN(status=1) endif if(allocated(c13_res_m)) then comp='c13_res' dsname='flux_'//trim(comp)//'_'//trim(pp)//'_mean' attname='flux_'//trim(comp)//'_'//trim(pp)//'_n' call MDF_Def_Var( hid, trim(dsname), MDF_FLOAT, (/dim_lon_id,dim_lat_id/), var_flux_c13_res, status ) IF_NOTOK_RETURN(status=1) call MDF_Put_Att(hid, var_flux_c13_res,"units",values="mol m-2 s-1",status=status) IF_NOTOK_RETURN(status=1) call MDF_Put_Att(hid, var_flux_c13_res,attname,values=c13_res_n,status=status) IF_NOTOK_RETURN(status=1) endif if(allocated(c13_ocean_m)) then comp='c13_ocean' dsname='flux_'//trim(comp)//'_'//trim(pp)//'_mean' attname='flux_'//trim(comp)//'_'//trim(pp)//'_n' call MDF_Def_Var( hid, trim(dsname), MDF_FLOAT, (/dim_lon_id,dim_lat_id/), var_flux_c13_ocean, status ) IF_NOTOK_RETURN(status=1) call MDF_Put_Att(hid, var_flux_c13_ocean,"units",values="mol m-2 s-1",status=status) IF_NOTOK_RETURN(status=1) call MDF_Put_Att(hid, var_flux_c13_ocean,attname,values=c13_ocean_n,status=status) IF_NOTOK_RETURN(status=1) endif if(allocated(c13_ff_m)) then comp='c13_ff' dsname='flux_'//trim(comp)//'_'//trim(pp)//'_mean' attname='flux_'//trim(comp)//'_'//trim(pp)//'_n' call MDF_Def_Var( hid, trim(dsname), MDF_FLOAT, (/dim_lon_id,dim_lat_id/), var_flux_c13_ff, status ) IF_NOTOK_RETURN(status=1) call MDF_Put_Att(hid, var_flux_c13_ff,"units",values="mol m-2 s-1",status=status) IF_NOTOK_RETURN(status=1) call MDF_Put_Att(hid, var_flux_c13_ff,attname,values=c13_ff_n,status=status) IF_NOTOK_RETURN(status=1) endif if(allocated(c13_fires_m)) then comp='c13_fires' dsname='flux_'//trim(comp)//'_'//trim(pp)//'_mean' attname='flux_'//trim(comp)//'_'//trim(pp)//'_n' call MDF_Def_Var( hid, trim(dsname), MDF_FLOAT, (/dim_lon_id,dim_lat_id/), var_flux_c13_fires, status ) IF_NOTOK_RETURN(status=1) call MDF_Put_Att(hid, var_flux_c13_fires,"units",values="mol m-2 s-1",status=status) IF_NOTOK_RETURN(status=1) call MDF_Put_Att(hid, var_flux_c13_fires,attname,values=c13_fires_n,status=status) IF_NOTOK_RETURN(status=1) endif if(allocated(c13_biodis_m)) then comp='c13_biodis' dsname='flux_'//trim(comp)//'_'//trim(pp)//'_mean' attname='flux_'//trim(comp)//'_'//trim(pp)//'_n' call MDF_Def_Var( hid, trim(dsname), MDF_FLOAT, (/dim_lon_id,dim_lat_id/), var_flux_c13_biodis, status ) IF_NOTOK_RETURN(status=1) call MDF_Put_Att(hid, var_flux_c13_biodis,"units",values="mol m-2 s-1",status=status) IF_NOTOK_RETURN(status=1) call MDF_Put_Att(hid, var_flux_c13_biodis,attname,values=c13_biodis_n,status=status) IF_NOTOK_RETURN(status=1) endif if(allocated(c13_oceandis_m)) then comp='c13_oceandis' dsname='flux_'//trim(comp)//'_'//trim(pp)//'_mean' attname='flux_'//trim(comp)//'_'//trim(pp)//'_n' call MDF_Def_Var( hid, trim(dsname), MDF_FLOAT, (/dim_lon_id,dim_lat_id/), var_flux_c13_oceandis, status ) IF_NOTOK_RETURN(status=1) call MDF_Put_Att(hid, var_flux_c13_oceandis,"units",values="mol m-2 s-1",status=status) IF_NOTOK_RETURN(status=1) call MDF_Put_Att(hid, var_flux_c13_oceandis,attname,values=c13_oceandis_n,status=status) IF_NOTOK_RETURN(status=1) endif ! finished definition: call MDF_EndDef( hid, status ) IF_NOTOK_RETURN(status=1) call MDF_Put_Var( hid, var_lon, (/(i+0.5,i=-180,179)/), status) IF_NOTOK_RETURN(status=1) call MDF_Put_Var( hid, var_lat, (/(i+0.5,i=-90,89)/), status) IF_NOTOK_RETURN(status=1) call MDF_Put_Var( hid, var_multiplier, fluxmultiplier(:,:,1), status) IF_NOTOK_RETURN(status=1) call MDF_Put_Var( hid, var_multiplier_epsbio, fluxmultiplier_epsbio(:,:,1), status) IF_NOTOK_RETURN(status=1) if(allocated(co2_all_m)) then call MDF_Put_Var( hid, var_flux_all, co2_all_m, status) IF_NOTOK_RETURN(status=1) endif if(allocated(co2_ff_m)) then call MDF_Put_Var( hid, var_flux_ff, co2_ff_m, status) IF_NOTOK_RETURN(status=1) endif if(allocated(co2_bio_m)) then call MDF_Put_Var( hid, var_flux_bio, co2_bio_m, status) IF_NOTOK_RETURN(status=1) endif if(allocated(field2d_eps_m)) then call MDF_Put_Var( hid, var_eps, field2d_eps_m, status) IF_NOTOK_RETURN(status=1) endif if(allocated(co2_ocean_m)) then call MDF_Put_Var( hid, var_flux_ocean, co2_ocean_m, status) IF_NOTOK_RETURN(status=1) endif if(allocated(co2_fires_m)) then call MDF_Put_Var( hid, var_flux_fires, co2_fires_m, status) IF_NOTOK_RETURN(status=1) endif if(allocated(c13_all_m)) then call MDF_Put_Var( hid, var_flux_c13_all, c13_all_m, status) IF_NOTOK_RETURN(status=1) endif if(allocated(c13_bio_m)) then call MDF_Put_Var( hid, var_flux_c13_bio, c13_bio_m, status) IF_NOTOK_RETURN(status=1) endif if(allocated(c13_gpp_m)) then call MDF_Put_Var( hid, var_flux_c13_gpp, c13_gpp_m, status) IF_NOTOK_RETURN(status=1) endif if(allocated(co2_gpp_m)) then call MDF_Put_Var( hid, var_flux_co2_gpp, co2_gpp_m, status) IF_NOTOK_RETURN(status=1) endif if(allocated(c13_res_m)) then call MDF_Put_Var( hid, var_flux_c13_res, c13_res_m, status) IF_NOTOK_RETURN(status=1) endif if(allocated(c13_ocean_m)) then call MDF_Put_Var( hid, var_flux_c13_ocean, c13_ocean_m, status) IF_NOTOK_RETURN(status=1) print*,'PUT VAR',sum(c13_ocean_m) endif if(allocated(c13_ff_m)) then call MDF_Put_Var( hid, var_flux_c13_ff, c13_ff_m, status) IF_NOTOK_RETURN(status=1) print*,'PUT VAR',sum(c13_ff_m) endif if(allocated(c13_fires_m)) then call MDF_Put_Var( hid, var_flux_c13_fires, c13_fires_m, status) IF_NOTOK_RETURN(status=1) print*,'PUT VAR',sum(c13_fires_m) endif if(allocated(c13_biodis_m)) then call MDF_Put_Var( hid, var_flux_c13_biodis, c13_biodis_m, status) IF_NOTOK_RETURN(status=1) print*,'PUT VAR',sum(c13_biodis_m) endif if(allocated(c13_oceandis_m)) then call MDF_Put_Var( hid, var_flux_c13_oceandis, c13_oceandis_m, status) IF_NOTOK_RETURN(status=1) print*,'PUT VAR',sum(c13_oceandis_m) endif ! close file: call MDF_Close( hid, status ) IF_NOTOK_RETURN(status=1) write (gol,'("[save_emission_averagers] Done writing arrays and closing output file.")'); call goPr end subroutine save_emission_averagers subroutine reset_emission_averagers if(allocated(co2_all_m)) then co2_all_m=0. co2_all_n=0 endif if(allocated(co2_ocean_m)) then co2_ocean_m=0. co2_ocean_n=0 endif if(allocated(co2_fires_m)) then co2_fires_m=0. co2_fires_n=0 endif if(allocated(co2_bio_m)) then co2_bio_m=0. co2_bio_n=0 endif if(allocated(field2d_eps_m)) then field2d_eps_m=0. field2d_eps_n=0 endif if(allocated(co2_ff_m)) then co2_ff_m=0. co2_ff_n=0 endif if(allocated(c13_all_m)) then c13_all_m=0. c13_all_n=0 endif if(allocated(c13_bio_m)) then c13_bio_m=0. c13_bio_n=0 endif if(allocated(c13_gpp_m)) then c13_gpp_m=0. c13_gpp_n=0 endif if(allocated(co2_gpp_m)) then co2_gpp_m=0. co2_gpp_n=0 endif if(allocated(c13_res_m)) then c13_res_m=0. c13_res_n=0 endif if(allocated(c13_ocean_m)) then c13_ocean_m=0. c13_ocean_n=0 endif if(allocated(c13_ff_m)) then c13_ff_m=0. c13_ff_n=0 endif if(allocated(c13_fires_m)) then c13_fires_m=0. c13_fires_n=0 endif if(allocated(c13_biodis_m)) then c13_biodis_m=0. c13_biodis_n=0 endif if(allocated(c13_oceandis_m)) then c13_oceandis_m=0. c13_oceandis_n=0 endif end subroutine reset_emission_averagers subroutine ReadFromFile3D( filename, sds_name , data_array, status, netcdf4) use GO , only : gol, goErr, goPr, goLabel, goUpCase !use file_hdf use MDF !I/O character(len=*),intent(in) :: filename character(len=*),intent(in) :: sds_name real,dimension(:,:,:),intent(out):: data_array logical,optional:: netcdf4 integer, intent (inout) :: status character(len=64) :: varname integer :: hid integer :: n, k, l character(len=*), parameter :: rname = mname//'/ReadFromFile' if (present(netcdf4)) then call MDF_Open(filename=trim(filename), ftype=MDF_NETCDF4, mode=MDF_READ, hid=hid, status=status) IF_NOTOK_RETURN(status=1) else call MDF_Open(filename=trim(filename), ftype=MDF_HDF4, mode=MDF_READ, hid=hid, status=status) IF_NOTOK_RETURN(status=1) endif call MDF_Inquire(hid, status, nVariables = n ) IF_NOTOK_RETURN(status=1) k = 0 do k = k + 1 if ( k > n ) then write(gol,*) 'SDS not found' ; call goPr write(gol,*) 'filename : ', trim(filename) ; call goPr IF_NOTOK_RETURN(status=1) endif call MDF_Inquire_Variable( hid, k, status, name=varname ) IF_NOTOK_RETURN(status=1) if (trim(sds_name) .ne. trim(varname) ) cycle ! name not correct ! all lights green, go read call MDF_Get_Var( hid, k , data_array, status ) IF_ERROR_RETURN(status=1) exit ! break from loop end do call MDF_Close( hid, status ) IF_ERROR_RETURN(status=1) status=0 end subroutine ReadFromFile3D subroutine ReadFromFile( filename, sds_name , data_array, status, sds_date, netcdf4) use GO , only : gol, goErr, goPr, goLabel, goUpCase !use file_hdf use MDF !I/O character(len=*),intent(in) :: filename character(len=*),intent(in) :: sds_name real,dimension(:,:),intent(out):: data_array integer,dimension(6),optional:: sds_date logical,optional:: netcdf4 integer, intent (inout) :: status character(len=64) :: varname integer :: hid integer :: n, k, l character(len=*), parameter :: rname = mname//'/ReadFromFile' integer,dimension(6) :: idate_f if (present(netcdf4)) then call MDF_Open(filename=trim(filename), ftype=MDF_NETCDF4, mode=MDF_READ, hid=hid, status=status) IF_NOTOK_RETURN(status=1) else call MDF_Open(filename=trim(filename), ftype=MDF_HDF4, mode=MDF_READ, hid=hid, status=status) IF_NOTOK_RETURN(status=1) endif call MDF_Inquire(hid, status, nVariables = n ) IF_NOTOK_RETURN(status=1) k = 0 do k = k + 1 if ( k > n ) then write(gol,*) 'SDS not found' ; call goPr write(gol,*) 'filename : ', trim(filename) ; call goPr write(gol,*) 'sds_name : ', trim(sds_name) ; call goPr if (present(sds_date)) write(gol,*) 'sds_date :',sds_date ; call goPr IF_NOTOK_RETURN(status=1) endif call MDF_Inquire_Variable( hid, k, status, name=varname ) IF_NOTOK_RETURN(status=1) if (trim(sds_name) .ne. trim(varname) ) cycle ! name not correct if (present(sds_date)) then call MDF_Get_Att( hid, k, 'idate', idate_f, status) IF_NOTOK_RETURN(status=1) do l = 1,6 if (sds_date(l) .ne. idate_f(l) ) status = 1 enddo if (status == 1) cycle endif ! all lights green, go read call MDF_Get_Var( hid, k , data_array, status ) IF_ERROR_RETURN(status=1) exit ! break from loop end do call MDF_Close( hid, status ) IF_ERROR_RETURN(status=1) status=0 end subroutine ReadFromFile subroutine GetParameterMaps(status) use GO , only : goLoCase,gol, goErr, goPr, goBug, goNum2str use MDF , only : MDF_Open use MDF , only : MDF_Create, MDF_Close, MDF_EndDef use MDF , only : MDF_NETCDF4, MDF_REPLACE, MDF_GLOBAL, MDF_CHAR, MDF_INT, MDF_FLOAT, MDF_UNLIMITED, MDF_READ use MDF , only : MDF_Put_Att use MDF , only : MDF_Def_Dim use MDF , only : MDF_Inq_Dimid, MDF_Inquire_Dimension use MDF , only : MDF_Inq_VarID use MDF , only : MDF_Def_Var, MDF_Put_Var, MDF_Get_Var use Partools, only: ntracetloc, ntracet_ar, myid use chem_param, only : nmembersloc,names implicit none ! ------------- I/O ------------------ integer, intent(out) :: status ! ------------- local ------------------ character(len=*), parameter :: rname = mname//'/GetParameterMaps' character(len=1024) :: inFile = '' logical :: file_exists integer :: n, offsetn, nloop integer :: hid,dimid, varid !WP! This routine opens a file with 1x1 degree parameters to be used in !WP! CarbonTracker. For now, these are assumed to be scaling factors that are !WP! applied to the ocean and biosphere at each time step. We read one file !WP! for every ensemble member (inv) or one just for the mean multipliers !WP! (fwd). !WP! !WP! The file names are constructed from a root file name that comes from !WP! the rc-file and is read into the ParamInputFile variable during Emis_Init. !WP! We add a qualifier with the tracer number to the filename to read one !WP! file per member for the inv runs. if (enkf_is_inverse()) then write (gol,'("[GetParameterMaps] Parameter files will be read for each member ")') ; call goPr nloop=nmembersloc !ntracetloc !nloop=ntracetloc else write (gol,'("[GetParameterMaps] Parameter file will be read for mean ")') ; call goPr nloop=1 endif do n=1,nloop ! loop over tracers (members) on this processor !IV! This version of carbontracker uses two kinds of tracers: CO2 and 13CO2. !IV! It is designed to use an even number of processors. The first half of it are used !IV! for CO2 ensembles, 10 ensembles per processors, and the other half are used for 13CO2 ensembles. !IV! That requires that same parameters files must be used for both kinds of tracers. As such, for first 13CO2 !IV! ensemble the offsetn is set back to zero. offsetn = sum(ntracet_ar(0:myid-1)) if (goLoCase(trim(names(offsetn+1)))=='co2c13'.or.goLoCase(trim(names(offsetn+1)))=='co2c13_bg') then offsetn=offsetn-sum(ntracet_ar)/2 endif if (trim(ParamInputFile)=='None') then write (gol,'("[GetParameterMaps] ..... File specified is: ",a,".")') trim(ParamInputFile); call goPr write (gol,'("[GetParameterMaps] ..... Using parametervalues of 1.0 as default")'); call goPr fluxmultiplier = 1.0 fluxmultiplier_biodis = 1.0 fluxmultiplier_epsbio = 1.0 status=0 return endif inFile = trim(ParamInputFile)//'.'//trim(goNum2str(offsetn+n-1,fmt='(i3.3)'))//'.nc' INQUIRE(FILE=trim(inFile), EXIST=file_exists) if (.not.file_exists) then write (gol,'("[GetParameterMaps] ..... File to open does not exist ",a,".")') trim(inFile); call goErr status=1 IF_ERROR_RETURN(status=1) endif ! open existing file: call MDF_Open( trim(inFile), MDF_NETCDF4, MDF_READ, hid=hid, status=status ) IF_NOTOK_RETURN(status=1) call MDF_Inq_VarID(hid, "parametermap", varid, status) if (status /= 0) then write (gol,'("[GetParameterMaps] ..... required variable parametermap does not exist")') ; call goErr status=1 IF_NOTOK_RETURN(status=1) endif call MDF_Get_Var(hid,varid,fluxmultiplier(:,:,n),status) if (status /= 0) then write (gol,'("[GetParameterMaps] ..... failed to read variable parametermap")') ; call goErr status=1 IF_NOTOK_RETURN(status=1) endif ! !IV 2nd set parameter files call MDF_Inq_VarID(hid, "parametermap_epsbio", varid, status) if (status /= 0) then write (gol,'("[GetParameterMaps] ..... required variable parametermap does not exist")') ; call goErr status=1 IF_NOTOK_RETURN(status=1) endif call MDF_Get_Var(hid,varid,fluxmultiplier_epsbio(:,:,n),status) if (status /= 0) then write (gol,'("[GetParameterMaps] ..... failed to read variable parametermap")') ; call goErr status=1 IF_NOTOK_RETURN(status=1) endif call MDF_Close(hid,status) IF_NOTOK_RETURN(status=1) end do write (gol,'("[GetParameterMaps] ..... All parametermaps read successfully")') ; call goPr !ok status = 0 end subroutine GetParameterMaps subroutine get_ntimesteps(ntimesteps) use dims, only: idate implicit none integer, intent(out) :: ntimesteps if(idate(2) .eq. 1 .or. idate(2) .eq. 3 .or. idate(2) .eq. 5 .or. idate(2) .eq. 7 .or. idate(2) .eq. 8 .or. idate(2) .eq. 10 .or. idate(2) .eq. 12) then ntimesteps = 248 else if (idate(2) .eq. 4 .or. idate(2) .eq. 6 .or. idate(2) .eq. 9 .or. idate(2) .eq. 11) then ntimesteps = 240 else if (idate(2) .eq. 2) then ntimesteps = 224 endif end subroutine get_ntimesteps end module emission_common