! ================================================================================================================================= ! PROGRAM : forcesoil ! ! CONTACT : orchidee-help _at_ ipsl.jussieu.fr ! ! LICENCE : IPSL (2006). This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC ! !>\BRIEF This subroutine runs the soilcarbon submodel using specific initial conditions !! and driving variables in order to obtain soil carbon stocks closed to the steady-state values !! quicker than when using the ''full'' ORCHIDEE. !! !!\n DESCRIPTION: None !! This subroutine computes the soil carbon stocks by calling the soilcarbon routine at each time step. \n !! The aim is to obtain soil carbon stocks closed to the steady-state values and ultimately to create !! an updated stomate restart file for the stomate component. The state variables of the subsystem are the clay content !! (fixed value) and the soil carbon stocks. Initial conditions for the state variables are read in an !! input stomate restart file. Driving variables are Soil carbon input, Water and Temperature stresses on !! Organic Matter decomposition. Driving variables are read from a specific forcing file produced by a former run of ORCHIDEE !! (SECHIBA+STOMATE). \n !! The FORCESOIL program first consists in reading a set of input files, allocating variables and !! preparing output stomate restart file. \n !! Then, a loop over time is performed in which the soilcarbon routine is called at each time step. \n !! Last, final values of the soil carbon stocks are written into the output stomate restart file. \n !! No flag is associated with the use of the FORCESOIL program. \n !! !! RECENT CHANGE(S): None !! !! REFERENCE(S) : None !! !! FLOWCHART : None !! !! SVN : !! $HeadURL: $ !! $Date: $ !! $Revision: $ !! \n !_ ================================================================================================================================= PROGRAM forcesoil USE netcdf !- USE defprec USE constantes USE constantes_veg USE constantes_co2 USE stomate_constants USE ioipsl USE stomate_soilcarbon USE parallel !- IMPLICIT NONE !- CHARACTER(LEN=80) :: sto_restname_in,sto_restname_out INTEGER(i_std) :: iim,jjm !! Indices (unitless) INTEGER(i_std),PARAMETER :: llm = 1 !! Vertical Layers (requested by restini routine) (unitless) INTEGER(i_std) :: kjpindex !! Domain size (unitless) INTEGER(i_std) :: itau_dep,itau_len !! Time step read in the restart file (?) !! and number of time steps of the simulation (unitless) CHARACTER(LEN=30) :: time_str !! Length of the simulation (year) REAL(r_std) :: dt_files !! time step between two successive itaus (?) !! (requested by restini routine) (seconds) REAL(r_std) :: date0 !! Time at which itau = 0 (requested by restini routine) (?) INTEGER(i_std) :: rest_id_sto !! ID of the input restart file (unitless) CHARACTER(LEN=20), SAVE :: thecalendar = 'noleap' !! Type of calendar defined in the input restart file !! (unitless) !- CHARACTER(LEN=100) :: Cforcing_name !! Name of the forcing file (unitless) INTEGER :: Cforcing_id !! ID of the forcing file (unitless) INTEGER :: v_id !! ID of the variable 'Index' stored in the forcing file !! (unitless) REAL(r_std) :: dt_forcesoil !! Time step at which soilcarbon routine is called (days) INTEGER :: nparan !! Number of values stored per year in the forcing file !! (unitless) INTEGER(i_std),DIMENSION(:),ALLOCATABLE :: indices !! Grid Point Index used per processor (unitless) INTEGER(i_std),DIMENSION(:),ALLOCATABLE :: indices_g !! Grid Point Index for all processor (unitless) REAL(r_std),DIMENSION(:),ALLOCATABLE :: x_indices_g !! Grid Point Index for all processor (unitless) REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: lon, lat !! Longitude and Latitude of each grid point defined !! in lat/lon (2D) (degrees) REAL(r_std),DIMENSION(llm) :: lev !! Number of level (requested by restini routine) (unitless) INTEGER :: i,m,iatt,iv,iyear !! counters (unitless) CHARACTER(LEN=80) :: var_name !! temporary variable to store string (unitless) CHARACTER(LEN=400) :: taboo_vars !! string used for storing the name of the variables !! of the stomate restart file that are not automatically !! duplicated from input to output restart file (unitless) REAL(r_std),DIMENSION(1) :: xtmp !! scalar read/written in restget/restput routines (unitless) INTEGER(i_std),PARAMETER :: nbvarmax=300 !! maximum # of variables assumed in the stomate restart file !! (unitless) INTEGER(i_std) :: nbvar !! # of variables effectively present !! in the stomate restart file (unitless) CHARACTER(LEN=50),DIMENSION(nbvarmax) :: varnames !! list of the names of the variables stored !! in the stomate restart file (unitless) INTEGER(i_std) :: varnbdim !! # of dimensions of a given variable !! of the stomate restart file INTEGER(i_std),PARAMETER :: varnbdim_max=20 !! maximal # of dimensions assumed for any variable !! of the stomate restart file INTEGER,DIMENSION(varnbdim_max) :: vardims !! length of each dimension of a given variable !! of the stomate restart file LOGICAL :: l1d !! boolean : TRUE if all dimensions of a given variable !! of the stomate restart file are of length 1 (ie scalar) !! (unitless) REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: var_3d !! matrix read/written in restget/restput routines (unitless) REAL(r_std) :: x_tmp !! temporary variable used to store return value !! from nf90_get_att (unitless) CHARACTER(LEN=10) :: part_str !! string suffix indicating the index of a PFT REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:) :: clay_g !! clay fraction (nbpglo) (unitless) REAL(r_std),DIMENSION(:,:,:,:),ALLOCATABLE :: soilcarbon_input_g !! soil carbon input (nbpglob,ncarb,nvm,time) !! (\f$gC m^{-2} dt_forcesoil^{-1}\f$) REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: control_temp_g !! Temperature control (nbp_glo,above/below,time) on OM decomposition !! (unitless) REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: control_moist_g !! Moisture control (nbp_glo,abo/below,time) on OM decomposition !! ?? Should be defined per PFT as well (unitless) REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: carbon_g !! Soil carbon stocks (nbp_glo,ncarb,nvm) (\f$gC m^{-2}\f$) REAL(r_std),ALLOCATABLE :: clay(:) !! clay fraction (nbp_loc) (unitless) REAL(r_std),ALLOCATABLE :: soilcarbon_input(:,:,:,:) !! soil carbon input (nbp_loc,ncarb,nvm,time) !! (\f$gC m^{-2} dt_forcesoil^{-1}\f$) REAL(r_std),ALLOCATABLE :: control_temp(:,:,:) !! Temperature control (nbp_loc,above/below,time) on OM decomposition !! (unitless) REAL(r_std),ALLOCATABLE :: control_moist(:,:,:) !! Moisture control (nbp_loc,abo/below,time) on OM decomposition !! ?? Should be defined per PFT as well (unitless) REAL(r_std),ALLOCATABLE :: carbon(:,:,:) !! Soil carbon stocks (nbp_loc,ncarb,nvm) (\f$gC m^{-2}\f$) REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: resp_hetero_soil !! Heterotrophic respiration (\f$gC m^{-2} dt_forcesoil^{-1}\f$) !! (requested by soilcarbon routine but not used here) INTEGER(i_std) :: ier,iret !! Used for hangling errors LOGICAL :: debug !! boolean used for printing messages CALL Init_para(.FALSE.) CALL init_timer !--------------------------------------------------------------------- !- ! set debug to have more information !- !Config Key = DEBUG_INFO !Config Desc = Flag for debug information !Config Def = n !Config Help = This option allows to switch on the output of debug !Config information without recompiling the code. !- debug = .FALSE. CALL getin_p('DEBUG_INFO',debug) !!- !! 1. Initialisation stage !! Reading a set of input files, allocating variables and preparing output restart file. !!- ! Define restart file name ! for reading initial conditions (sto_restname_in var) and for writting final conditions (sto_restname_out var). ! User values are used if present in the .def file. ! If not present, default values (stomate_start.nc and stomate_rest_out.c) are used. !- IF (is_root_prc) THEN sto_restname_in = 'stomate_start.nc' CALL getin ('STOMATE_RESTART_FILEIN',sto_restname_in) WRITE(numout,*) 'STOMATE INPUT RESTART_FILE: ',TRIM(sto_restname_in) sto_restname_out = 'stomate_rest_out.nc' CALL getin ('STOMATE_RESTART_FILEOUT',sto_restname_out) WRITE(numout,*) 'STOMATE OUTPUT RESTART_FILE: ',TRIM(sto_restname_out) !- ! Open the input file and Get some Dimension and Attributes ID's !- iret = NF90_OPEN (sto_restname_in, NF90_NOWRITE, rest_id_sto) iret = NF90_INQUIRE_DIMENSION (rest_id_sto,1,len=iim_g) iret = NF90_INQUIRE_DIMENSION (rest_id_sto,2,len=jjm_g) iret = NF90_INQ_VARID (rest_id_sto, "time", iv) iret = NF90_GET_ATT (rest_id_sto, iv, 'calendar',thecalendar) iret = NF90_CLOSE (rest_id_sto) i=INDEX(thecalendar,ACHAR(0)) IF ( i > 0 ) THEN thecalendar(i:20)=' ' ENDIF !- ! Allocate longitudes and latitudes !- ALLOCATE (lon(iim_g,jjm_g)) ALLOCATE (lat(iim_g,jjm_g)) lon(:,:) = 0.0 lat(:,:) = 0.0 lev(1) = 0.0 !- CALL restini & & (sto_restname_in, iim_g, jjm_g, lon, lat, llm, lev, & & sto_restname_out, itau_dep, date0, dt_files, rest_id_sto) ENDIF CALL bcast(date0) CALL bcast(thecalendar) WRITE(numout,*) "calendar = ",thecalendar !- ! calendar !- CALL ioconf_calendar (thecalendar) CALL ioget_calendar (one_year,one_day) CALL ioconf_startdate(date0) ! !! For master process only ! IF (is_root_prc) THEN !- ! define forcing file's name (Cforcing_name var) ! User value is used if present in the .def file ! If not, default (NONE) is used !- Cforcing_name = 'NONE' CALL getin ('STOMATE_CFORCING_NAME',Cforcing_name) !- ! Open FORCESOIL's forcing file to read some basic info (dimensions, variable ID's) ! and allocate variables. !- iret = NF90_OPEN (TRIM(Cforcing_name),NF90_NOWRITE,Cforcing_id) IF (iret /= NF90_NOERR) THEN CALL ipslerr (3,'forcesoil', & & 'Could not open file : ', & & Cforcing_name,'(Do you have forget it ?)') ENDIF !- ! Total Domain size is stored in nbp_glo variable !- ier = NF90_GET_ATT (Cforcing_id,NF90_GLOBAL,'kjpindex',x_tmp) nbp_glo = NINT(x_tmp) !- ! Number of values stored per year in the forcing file is stored in nparan var. !- ier = NF90_GET_ATT (Cforcing_id,NF90_GLOBAL,'nparan',x_tmp) nparan = NINT(x_tmp) !- ALLOCATE (indices_g(nbp_glo)) ALLOCATE (clay_g(nbp_glo)) !- ALLOCATE (x_indices_g(nbp_glo),stat=ier) ier = NF90_INQ_VARID (Cforcing_id,'index',v_id) ier = NF90_GET_VAR (Cforcing_id,v_id,x_indices_g) indices_g(:) = NINT(x_indices_g(:)) WRITE(numout,*) mpi_rank,"indices globaux : ",indices_g DEALLOCATE (x_indices_g) !- ier = NF90_INQ_VARID (Cforcing_id,'clay',v_id) ier = NF90_GET_VAR (Cforcing_id,v_id,clay_g) !- ! time step of forcesoil program (in days) !- dt_forcesoil = one_year / FLOAT(nparan) WRITE(numout,*) 'time step (d): ',dt_forcesoil !- ! read and write the variables in the output restart file we do not modify within the Forcesoil program ! ie all variables stored in the input restart file except those stored in taboo_vars !- taboo_vars = '$lon$ $lat$ $lev$ $nav_lon$ $nav_lat$ $nav_lev$ $time$ $time_steps$ '// & & '$day_counter$ $dt_days$ $date$ '// & & '$carbon_01$ $carbon_02$ $carbon_03$ $carbon_04$ $carbon_05$'// & & '$carbon_06$ $carbon_07$ $carbon_08$ $carbon_09$ $carbon_10$'// & & '$carbon_11$ $carbon_12$ $carbon_13$' !- CALL ioget_vname(rest_id_sto, nbvar, varnames) !- ! read and write some special variables (1D or variables that we need) !- var_name = 'day_counter' CALL restget (rest_id_sto, var_name, 1, 1, 1, itau_dep, .TRUE., xtmp) CALL restput (rest_id_sto, var_name, 1, 1, 1, itau_dep, xtmp) !- var_name = 'dt_days' CALL restget (rest_id_sto, var_name, 1, 1, 1, itau_dep, .TRUE., xtmp) CALL restput (rest_id_sto, var_name, 1, 1, 1, itau_dep, xtmp) !- var_name = 'date' CALL restget (rest_id_sto, var_name, 1, 1, 1, itau_dep, .TRUE., xtmp) CALL restput (rest_id_sto, var_name, 1, 1, 1, itau_dep, xtmp) !- DO iv=1,nbvar !-- check if the variable is to be written here IF (INDEX(taboo_vars,'$'//TRIM(varnames(iv))//'$') == 0 ) THEN !---- get variable dimensions, especially 3rd dimension CALL ioget_vdim & & (rest_id_sto, varnames(iv), varnbdim_max, varnbdim, vardims) l1d = ALL(vardims(1:varnbdim) == 1) !---- read it IF (l1d) THEN CALL restget & & (rest_id_sto, TRIM(varnames(iv)), 1, vardims(3), & & 1, itau_dep, .TRUE., xtmp) ELSE ALLOCATE( var_3d(nbp_glo,vardims(3)), stat=ier) IF (ier /= 0) STOP 'ALLOCATION PROBLEM' !---- CALL restget & & (rest_id_sto, TRIM(varnames(iv)), nbp_glo, vardims(3), & & 1, itau_dep, .TRUE., var_3d, "gather", nbp_glo, indices_g) ENDIF !---- write it IF (l1d) THEN CALL restput & & (rest_id_sto, TRIM(varnames(iv)), 1, vardims(3), & & 1, itau_dep, xtmp) ELSE CALL restput & & (rest_id_sto, TRIM(varnames(iv)), nbp_glo, vardims(3), & & 1, itau_dep, var_3d, 'scatter', nbp_glo, indices_g) !---- DEALLOCATE(var_3d) ENDIF ENDIF ENDDO !- ! read soil carbon stocks values stored in the input restart file !- ALLOCATE(carbon_g(nbp_glo,ncarb,nvm)) carbon_g(:,:,:) = val_exp DO m = 1, nvm WRITE (part_str, '(I2)') m IF (m<10) part_str(1:1)='0' var_name = 'carbon_'//part_str(1:LEN_TRIM(part_str)) CALL restget & & (rest_id_sto, var_name, nbp_glo, ncarb , 1, itau_dep, & & .TRUE., carbon_g(:,:,m), 'gather', nbp_glo, indices_g) IF (ALL(carbon_g(:,:,m) == val_exp)) carbon_g(:,:,m) = zero !-- do not write this variable: it will be modified. ENDDO WRITE(numout,*) "date0 : ",date0, itau_dep !- ! Length of the run (in Years) ! User value is used if present in the .def file ! If not, default value (10000 Years) is used !- WRITE(time_str,'(a)') '10000Y' CALL getin('TIME_LENGTH', time_str) write(numout,*) 'Number of years for carbon spinup : ',time_str ! transform into itau CALL tlen2itau(time_str, dt_forcesoil*one_day, date0, itau_len) write(numout,*) 'Number of time steps to do: ',itau_len !- ! read soil carbon inputs, water and temperature stresses on OM decomposition ! into the forcing file - We read an average year. !?! It should not be an average year !?! ORCHIDEE is capable to produce multi-annual Cforcing files !?! time dimension for soilcarbon_input_g, control_moist_g, control_temp_g !?! should not be nparan but nparan*nb_year !!?? I completely agree. !- ALLOCATE(soilcarbon_input_g(nbp_glo,ncarb,nvm,nparan)) ALLOCATE(control_temp_g(nbp_glo,nlevs,nparan)) ALLOCATE(control_moist_g(nbp_glo,nlevs,nparan)) !- ier = NF90_INQ_VARID (Cforcing_id,'soilcarbon_input',v_id) ier = NF90_GET_VAR (Cforcing_id,v_id,soilcarbon_input_g) ier = NF90_INQ_VARID (Cforcing_id, 'control_moist',v_id) ier = NF90_GET_VAR (Cforcing_id,v_id,control_moist_g) ier = NF90_INQ_VARID (Cforcing_id, 'control_temp',v_id) ier = NF90_GET_VAR (Cforcing_id,v_id,control_temp_g) !- ier = NF90_CLOSE (Cforcing_id) !- ENDIF CALL bcast(nparan) CALL bcast(dt_forcesoil) CALL bcast(iim_g) CALL bcast(jjm_g) call bcast(nbp_glo) CALL bcast(itau_dep) CALL bcast(itau_len) ! ! We must initialize data_para : ! ! CALL init_data_para(iim_g,jjm_g,nbp_glo,indices_g) kjpindex=nbp_loc jjm=jj_nb iim=iim_g IF (debug) WRITE(numout,*) "Local grid : ",kjpindex,iim,jjm !--- !--- Create the index table !--- !--- This job returns a LOCAL kindex. !--- ALLOCATE (indices(kjpindex),stat=ier) ! !! scattering to all processes in parallel mode ! CALL scatter(indices_g,indices) indices(1:kjpindex)=indices(1:kjpindex)-(jj_begin-1)*iim_g IF (debug) WRITE(numout,*) mpi_rank,"indices locaux = ",indices(1:kjpindex) !- ! Allocation of the variables for a processor !- ALLOCATE(clay(kjpindex)) !?! It should not be an average year !?! ORCHIDEE is capable to produce multi-annual Cforcing files !?! time dimension for soilcarbon_input, control_moist, control_temp !?! should not be nparan but nparan*nb_year ALLOCATE(soilcarbon_input(kjpindex,ncarb,nvm,nparan)) ALLOCATE(control_temp(kjpindex,nlevs,nparan)) ALLOCATE(control_moist(kjpindex,nlevs,nparan)) ALLOCATE(carbon(kjpindex,ncarb,nvm)) ALLOCATE(resp_hetero_soil(kjpindex,nvm)) iatt = 0 !- ! Initialization of the variables for a processor !- CALL Scatter(clay_g,clay) CALL Scatter(soilcarbon_input_g,soilcarbon_input) CALL Scatter(control_temp_g,control_temp) CALL Scatter(control_moist_g,control_moist) CALL Scatter(carbon_g,carbon) !!- !! 2. Computational step !! Loop over time - Call of soilcarbon routine at each time step !! Updated soil carbon stocks are stored into carbon variable !! We only keep the last value of carbon variable (no time dimension). !!- iyear=1 DO i=1,itau_len iatt = iatt+1 IF (iatt > nparan) THEN IF (debug) WRITE(numout,*) iyear iatt = 1 iyear=iyear+1 ENDIF !?! Time index should not be iatt !?! but iatt+(iyear-1)*nparan or index should be i CALL soilcarbon & & (kjpindex, dt_forcesoil, clay, & & soilcarbon_input(:,:,:,iatt), & & control_temp(:,:,iatt), control_moist(:,:,iatt), & & carbon, resp_hetero_soil) ENDDO WRITE(numout,*) "End of soilcarbon LOOP." ! !! Gathering of variables towards main processor in parallel mode ! CALL Gather(carbon,carbon_g) !!- !! 3. write new carbon stocks into the ouput restart file !!- IF (is_root_prc) THEN DO m=1,nvm WRITE (part_str, '(I2)') m IF (m<10) part_str(1:1)='0' var_name = 'carbon_'//part_str(1:LEN_TRIM(part_str)) CALL restput & & (rest_id_sto, var_name, nbp_glo, ncarb , 1, itau_dep, & & carbon_g(:,:,m), 'scatter', nbp_glo, indices_g) ENDDO !- CALL getin_dump CALL restclo ENDIF #ifdef CPP_PARA CALL MPI_FINALIZE(ier) #endif WRITE(numout,*) "End of forcesoil." !-------------------- END PROGRAM forcesoil