[7541] | 1 | ! ================================================================================================================================= |
---|
| 2 | ! PROGRAM : forcesoil |
---|
| 3 | ! |
---|
| 4 | ! CONTACT : orchidee-help _at_ listes.ipsl.fr |
---|
| 5 | ! |
---|
| 6 | ! LICENCE : IPSL (2006). This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC |
---|
| 7 | ! |
---|
| 8 | !>\BRIEF This subroutine runs the soilcarbon submodel using specific initial conditions |
---|
| 9 | !! and driving variables in order to obtain soil carbon stocks closed to the steady-state values |
---|
| 10 | !! quicker than when using the ''full'' ORCHIDEE. |
---|
| 11 | !! |
---|
| 12 | !!\n DESCRIPTION: None |
---|
| 13 | !! This subroutine computes the soil carbon stocks by calling the soilcarbon routine at each time step. \n |
---|
| 14 | !! The aim is to obtain soil carbon stocks closed to the steady-state values and ultimately to create |
---|
| 15 | !! an updated stomate restart file for the stomate component. The state variables of the subsystem are the clay content |
---|
| 16 | !! (fixed value) and the soil carbon stocks. Initial conditions for the state variables are read in an |
---|
| 17 | !! input stomate restart file. Driving variables are Soil carbon input, Water and Temperature stresses on |
---|
| 18 | !! Organic Matter decomposition. Driving variables are read from a specific forcing file produced by a former run of ORCHIDEE |
---|
| 19 | !! (SECHIBA+STOMATE). \n |
---|
| 20 | !! The FORCESOIL program first consists in reading a set of input files, allocating variables and |
---|
| 21 | !! preparing output stomate restart file. \n |
---|
| 22 | !! Then, a loop over time is performed in which the soilcarbon routine is called at each time step. \n |
---|
| 23 | !! Last, final values of the soil carbon stocks are written into the output stomate restart file. \n |
---|
| 24 | !! No flag is associated with the use of the FORCESOIL program. \n |
---|
| 25 | !! |
---|
| 26 | !! RECENT CHANGE(S): None |
---|
| 27 | !! |
---|
| 28 | !! REFERENCE(S) : None |
---|
| 29 | !! |
---|
| 30 | !! FLOWCHART : None |
---|
| 31 | !! |
---|
| 32 | !! SVN : |
---|
| 33 | !! $HeadURL: $ |
---|
| 34 | !! $Date: $ |
---|
| 35 | !! $Revision: $ |
---|
| 36 | !! \n |
---|
| 37 | !_ ================================================================================================================================= |
---|
| 38 | |
---|
| 39 | PROGRAM forcesoil |
---|
| 40 | |
---|
| 41 | USE netcdf |
---|
| 42 | !- |
---|
| 43 | USE defprec |
---|
| 44 | USE constantes |
---|
| 45 | USE constantes_mtc |
---|
| 46 | USE pft_parameters |
---|
| 47 | USE stomate_data |
---|
| 48 | USE ioipsl_para |
---|
| 49 | USE mod_orchidee_para |
---|
| 50 | USE stomate_soilcarbon |
---|
| 51 | |
---|
| 52 | !- |
---|
| 53 | IMPLICIT NONE |
---|
| 54 | !- |
---|
| 55 | CHARACTER(LEN=80) :: sto_restname_in,sto_restname_out |
---|
| 56 | INTEGER(i_std) :: iim,jjm !! Indices (unitless) |
---|
| 57 | |
---|
| 58 | INTEGER(i_std),PARAMETER :: llm = 1 !! Vertical Layers (requested by restini routine) (unitless) |
---|
| 59 | INTEGER(i_std) :: kjpindex !! Domain size (unitless) |
---|
| 60 | |
---|
| 61 | INTEGER(i_std) :: itau_dep,itau_len !! Time step read in the restart file (?) |
---|
| 62 | !! and number of time steps of the simulation (unitless) |
---|
| 63 | CHARACTER(LEN=30) :: time_str !! Length of the simulation (year) |
---|
| 64 | REAL(r_std) :: dt_files !! time step between two successive itaus (?) |
---|
| 65 | !! (requested by restini routine) (seconds) |
---|
| 66 | REAL(r_std) :: date0 !! Time at which itau = 0 (requested by restini routine) (?) |
---|
| 67 | INTEGER(i_std) :: rest_id_sto !! ID of the input restart file (unitless) |
---|
| 68 | CHARACTER(LEN=20), SAVE :: thecalendar = 'noleap' !! Type of calendar defined in the input restart file |
---|
| 69 | !! (unitless) |
---|
| 70 | !- |
---|
| 71 | CHARACTER(LEN=100) :: Cforcing_name !! Name of the forcing file (unitless) |
---|
| 72 | INTEGER :: Cforcing_id !! ID of the forcing file (unitless) |
---|
| 73 | INTEGER :: v_id !! ID of the variable 'Index' stored in the forcing file |
---|
| 74 | !! (unitless) |
---|
| 75 | REAL(r_std) :: dt_forcesoil !! Time step at which soilcarbon routine is called (days) |
---|
| 76 | INTEGER :: nparan !! Number of values stored per year in the forcing file |
---|
| 77 | !! (unitless) |
---|
| 78 | INTEGER :: nbyear |
---|
| 79 | INTEGER(i_std),DIMENSION(:),ALLOCATABLE :: indices !! Grid Point Index used per processor (unitless) |
---|
| 80 | INTEGER(i_std),DIMENSION(:),ALLOCATABLE :: indices_g !! Grid Point Index for all processor (unitless) |
---|
| 81 | REAL(r_std),DIMENSION(:),ALLOCATABLE :: x_indices_g !! Grid Point Index for all processor (unitless) |
---|
| 82 | REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: lon, lat !! Longitude and Latitude of each grid point defined |
---|
| 83 | !! in lat/lon (2D) (degrees) |
---|
| 84 | REAL(r_std),DIMENSION(llm) :: lev !! Number of level (requested by restini routine) (unitless) |
---|
| 85 | |
---|
| 86 | |
---|
| 87 | INTEGER :: i,m,iatt,iv,iyear !! counters (unitless) |
---|
| 88 | |
---|
| 89 | CHARACTER(LEN=80) :: var_name |
---|
| 90 | CHARACTER(LEN=800) :: taboo_vars !! string used for storing the name of the variables |
---|
| 91 | !! of the stomate restart file that are not automatically |
---|
| 92 | !! duplicated from input to output restart file (unitless) |
---|
| 93 | REAL(r_std),DIMENSION(1) :: xtmp !! scalar read/written in restget/restput routines (unitless) |
---|
| 94 | INTEGER(i_std),PARAMETER :: nbvarmax=300 !! maximum # of variables assumed in the stomate restart file |
---|
| 95 | !! (unitless) |
---|
| 96 | INTEGER(i_std) :: nbvar !! # of variables effectively present |
---|
| 97 | !! in the stomate restart file (unitless) |
---|
| 98 | CHARACTER(LEN=50),DIMENSION(nbvarmax) :: varnames !! list of the names of the variables stored |
---|
| 99 | !! in the stomate restart file (unitless) |
---|
| 100 | INTEGER(i_std) :: varnbdim !! # of dimensions of a given variable |
---|
| 101 | !! of the stomate restart file |
---|
| 102 | INTEGER(i_std),PARAMETER :: varnbdim_max=20 !! maximal # of dimensions assumed for any variable |
---|
| 103 | !! of the stomate restart file |
---|
| 104 | INTEGER,DIMENSION(varnbdim_max) :: vardims !! length of each dimension of a given variable |
---|
| 105 | !! of the stomate restart file |
---|
| 106 | LOGICAL :: l1d !! boolean : TRUE if all dimensions of a given variable |
---|
| 107 | !! of the stomate restart file are of length 1 (ie scalar) |
---|
| 108 | !! (unitless) |
---|
| 109 | REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: var_3d !! matrix read/written in restget/restput routines (unitless) |
---|
| 110 | REAL(r_std) :: x_tmp !! temporary variable used to store return value |
---|
| 111 | !! from nf90_get_att (unitless) |
---|
| 112 | CHARACTER(LEN=10) :: part_str !! string suffix indicating the index of a PFT |
---|
| 113 | REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:) :: clay_g !! clay fraction (nbpglo) (unitless) |
---|
| 114 | REAL(r_std),DIMENSION(:,:,:,:),ALLOCATABLE :: soilcarbon_input_g !! soil carbon input (nbpglob,ncarb,nvm,time) |
---|
| 115 | !! (\f$gC m^{-2} dt_forcesoil^{-1}\f$) |
---|
| 116 | REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: control_temp_g !! Temperature control (nbp_glo,above/below,time) on OM decomposition |
---|
| 117 | !! (unitless) |
---|
| 118 | REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: control_moist_g !! Moisture control (nbp_glo,abo/below,time) on OM decomposition |
---|
| 119 | !! ?? Should be defined per PFT as well (unitless) |
---|
| 120 | REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: carbon_g !! Soil carbon stocks (nbp_glo,ncarb,nvm) (\f$gC m^{-2}\f$) |
---|
| 121 | |
---|
| 122 | REAL(r_std),ALLOCATABLE :: clay(:) !! clay fraction (nbp_loc) (unitless) |
---|
| 123 | REAL(r_std),ALLOCATABLE :: soilcarbon_input(:,:,:,:) !! soil carbon input (nbp_loc,ncarb,nvm,time) |
---|
| 124 | !! (\f$gC m^{-2} dt_forcesoil^{-1}\f$) |
---|
| 125 | REAL(r_std),ALLOCATABLE :: control_temp(:,:,:) !! Temperature control (nbp_loc,above/below,time) on OM decomposition |
---|
| 126 | !! (unitless) |
---|
| 127 | REAL(r_std),ALLOCATABLE :: control_moist(:,:,:) !! Moisture control (nbp_loc,abo/below,time) on OM decomposition |
---|
| 128 | !! ?? Should be defined per PFT as well (unitless) |
---|
| 129 | REAL(r_std),ALLOCATABLE :: carbon(:,:,:) !! Soil carbon stocks (nbp_loc,ncarb,nvm) (\f$gC m^{-2}\f$) |
---|
| 130 | REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: resp_hetero_soil !! Heterotrophic respiration (\f$gC m^{-2} dt_forcesoil^{-1}\f$) |
---|
| 131 | !! (requested by soilcarbon routine but not used here) |
---|
| 132 | |
---|
| 133 | INTEGER(i_std) :: printlev_loc !! Local write level |
---|
| 134 | INTEGER(i_std) :: ier,iret !! Used for hangling errors |
---|
| 135 | |
---|
| 136 | CHARACTER(LEN=30) :: temp_name |
---|
| 137 | LOGICAL :: l_error !! boolean for memory allocation |
---|
| 138 | REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:,:) :: matrixA !! Carbon fluxes matrix |
---|
| 139 | |
---|
| 140 | !_ ================================================================================================================================= |
---|
| 141 | |
---|
| 142 | CALL Init_orchidee_para |
---|
| 143 | CALL init_timer |
---|
| 144 | |
---|
| 145 | ! Set specific write level to forcesoil using PRINTLEV_forcesoil=[0-4] in run.def. |
---|
| 146 | ! The global printlev is used as default value. |
---|
| 147 | printlev_loc=get_printlev('forcesoil') |
---|
| 148 | |
---|
| 149 | !- |
---|
| 150 | ! Configure the number of PFTS |
---|
| 151 | !- |
---|
| 152 | |
---|
| 153 | ! 1. Read the number of PFTs |
---|
| 154 | ! |
---|
| 155 | !Config Key = NVM |
---|
| 156 | !Config Desc = number of PFTs |
---|
| 157 | !Config If = OK_SECHIBA or OK_STOMATE |
---|
| 158 | !Config Def = 13 |
---|
| 159 | !Config Help = The number of vegetation types define by the user |
---|
| 160 | !Config Units = [-] |
---|
| 161 | CALL getin_p('NVM',nvm) |
---|
| 162 | |
---|
| 163 | ! 2. Allocation |
---|
| 164 | l_error = .FALSE. |
---|
| 165 | ALLOCATE(pft_to_mtc(nvm),stat=ier) |
---|
| 166 | l_error = l_error .OR. (ier .NE. 0) |
---|
| 167 | IF (l_error) THEN |
---|
| 168 | STOP 'pft_to_mtc (forcesoil only) : error in memory allocation' |
---|
| 169 | ENDIF |
---|
| 170 | |
---|
| 171 | ! 3. Initialisation of the correspondance table |
---|
| 172 | pft_to_mtc(:) = undef_int |
---|
| 173 | |
---|
| 174 | ! 4.Reading of the conrrespondance table in the .def file |
---|
| 175 | ! |
---|
| 176 | !Config Key = PFT_TO_MTC |
---|
| 177 | !Config Desc = correspondance array linking a PFT to MTC |
---|
| 178 | !Config if = OK_SECHIBA or OK_STOMATE |
---|
| 179 | !Config Def = 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13 |
---|
| 180 | !Config Help = |
---|
| 181 | !Config Units = [-] |
---|
| 182 | CALL getin_p('PFT_TO_MTC',pft_to_mtc) |
---|
| 183 | |
---|
| 184 | ! 4.1 if nothing is found, we use the standard configuration |
---|
| 185 | IF(nvm == nvmc ) THEN |
---|
| 186 | IF(pft_to_mtc(1) == undef_int) THEN |
---|
| 187 | WRITE(numout,*) 'Note to the user : we will use ORCHIDEE to its standard configuration' |
---|
| 188 | pft_to_mtc(:) = (/ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13 /) |
---|
| 189 | ENDIF |
---|
| 190 | ELSE |
---|
| 191 | IF(pft_to_mtc(1) == undef_int) THEN |
---|
| 192 | WRITE(numout,*)' The array PFT_TO_MTC is empty : we stop' |
---|
| 193 | ENDIF |
---|
| 194 | ENDIF |
---|
| 195 | |
---|
| 196 | ! 4.2 What happened if pft_to_mtc(j) > nvmc (if the mtc doesn't exist)? |
---|
| 197 | DO i = 1, nvm |
---|
| 198 | IF(pft_to_mtc(i) > nvmc) THEN |
---|
| 199 | WRITE(numout,*) "the MTC you chose doesn't exist" |
---|
| 200 | STOP 'we stop reading pft_to_mtc' |
---|
| 201 | ENDIF |
---|
| 202 | ENDDO |
---|
| 203 | |
---|
| 204 | ! 4.3 Check if pft_to_mtc(1) = 1 |
---|
| 205 | IF(pft_to_mtc(1) /= 1) THEN |
---|
| 206 | WRITE(numout,*) 'the first pft has to be the bare soil' |
---|
| 207 | STOP 'we stop reading next values of pft_to_mtc' |
---|
| 208 | ELSE |
---|
| 209 | DO i = 2,nvm |
---|
| 210 | IF(pft_to_mtc(i) == 1) THEN |
---|
| 211 | WRITE(numout,*) 'only pft_to_mtc(1) has to be the bare soil' |
---|
| 212 | STOP 'we stop reading pft_to_mtc' |
---|
| 213 | ENDIF |
---|
| 214 | ENDDO |
---|
| 215 | ENDIF |
---|
| 216 | |
---|
| 217 | ! 5. Allocate and initialize natural ans is_c4 |
---|
| 218 | |
---|
| 219 | ! 5.1 Memory allocation |
---|
| 220 | l_error = .FALSE. |
---|
| 221 | ALLOCATE(natural(nvm),stat=ier) |
---|
| 222 | l_error = l_error .OR. (ier .NE. 0) |
---|
| 223 | ALLOCATE(is_c4(nvm),stat=ier) |
---|
| 224 | |
---|
| 225 | IF (l_error) THEN |
---|
| 226 | STOP 'natural or is_c4 (forcesoil only) : error in memory allocation' |
---|
| 227 | ENDIF |
---|
| 228 | |
---|
| 229 | ! 5.2 Initialisation |
---|
| 230 | DO i = 1, nvm |
---|
| 231 | natural(i) = natural_mtc(pft_to_mtc(i)) |
---|
| 232 | is_c4(i) = is_c4_mtc(pft_to_mtc(i)) |
---|
| 233 | ENDDO |
---|
| 234 | |
---|
| 235 | !!- |
---|
| 236 | !! 1. Initialisation stage |
---|
| 237 | !! Reading a set of input files, allocating variables and preparing output restart file. |
---|
| 238 | !!- |
---|
| 239 | ! Define restart file name |
---|
| 240 | ! for reading initial conditions (sto_restname_in var) and for writting final conditions (sto_restname_out var). |
---|
| 241 | ! User values are used if present in the .def file. |
---|
| 242 | ! If not present, default values (stomate_start.nc and stomate_rest_out.c) are used. |
---|
| 243 | !- |
---|
| 244 | IF (is_root_prc) THEN |
---|
| 245 | sto_restname_in = 'stomate_start.nc' |
---|
| 246 | CALL getin ('STOMATE_RESTART_FILEIN',sto_restname_in) |
---|
| 247 | WRITE(numout,*) 'STOMATE INPUT RESTART_FILE: ',TRIM(sto_restname_in) |
---|
| 248 | sto_restname_out = 'stomate_rest_out.nc' |
---|
| 249 | CALL getin ('STOMATE_RESTART_FILEOUT',sto_restname_out) |
---|
| 250 | WRITE(numout,*) 'STOMATE OUTPUT RESTART_FILE: ',TRIM(sto_restname_out) |
---|
| 251 | !- |
---|
| 252 | ! Open the input file and Get some Dimension and Attributes ID's |
---|
| 253 | !- |
---|
| 254 | iret = NF90_OPEN (sto_restname_in, NF90_NOWRITE, rest_id_sto) |
---|
| 255 | iret = NF90_INQUIRE_DIMENSION (rest_id_sto,1,len=iim_g) |
---|
| 256 | iret = NF90_INQUIRE_DIMENSION (rest_id_sto,2,len=jjm_g) |
---|
| 257 | iret = NF90_INQ_VARID (rest_id_sto, "time", iv) |
---|
| 258 | iret = NF90_GET_ATT (rest_id_sto, iv, 'calendar',thecalendar) |
---|
| 259 | iret = NF90_CLOSE (rest_id_sto) |
---|
| 260 | i=INDEX(thecalendar,ACHAR(0)) |
---|
| 261 | IF ( i > 0 ) THEN |
---|
| 262 | thecalendar(i:20)=' ' |
---|
| 263 | ENDIF |
---|
| 264 | !- |
---|
| 265 | ! Allocate longitudes and latitudes |
---|
| 266 | !- |
---|
| 267 | ALLOCATE (lon(iim_g,jjm_g)) |
---|
| 268 | ALLOCATE (lat(iim_g,jjm_g)) |
---|
| 269 | lon(:,:) = zero |
---|
| 270 | lat(:,:) = zero |
---|
| 271 | lev(1) = zero |
---|
| 272 | !- |
---|
| 273 | CALL restini & |
---|
| 274 | & (sto_restname_in, iim_g, jjm_g, lon, lat, llm, lev, & |
---|
| 275 | & sto_restname_out, itau_dep, date0, dt_files, rest_id_sto) |
---|
| 276 | ENDIF |
---|
| 277 | |
---|
| 278 | CALL bcast(date0) |
---|
| 279 | CALL bcast(thecalendar) |
---|
| 280 | WRITE(numout,*) "calendar = ",thecalendar |
---|
| 281 | !- |
---|
| 282 | ! calendar |
---|
| 283 | !- |
---|
| 284 | CALL ioconf_calendar (thecalendar) |
---|
| 285 | CALL ioget_calendar (one_year,one_day) |
---|
| 286 | CALL ioconf_startdate(date0) |
---|
| 287 | ! |
---|
| 288 | !! For master process only |
---|
| 289 | ! |
---|
| 290 | IF (is_root_prc) THEN |
---|
| 291 | !- |
---|
| 292 | ! define forcing file's name (Cforcing_name var) |
---|
| 293 | ! User value is used if present in the .def file |
---|
| 294 | ! If not, default (NONE) is used |
---|
| 295 | !- |
---|
| 296 | Cforcing_name = 'NONE' |
---|
| 297 | CALL getin ('STOMATE_CFORCING_NAME',Cforcing_name) |
---|
| 298 | !- |
---|
| 299 | ! Open FORCESOIL's forcing file to read some basic info (dimensions, variable ID's) |
---|
| 300 | ! and allocate variables. |
---|
| 301 | !- |
---|
| 302 | iret = NF90_OPEN (TRIM(Cforcing_name),NF90_NOWRITE,Cforcing_id) |
---|
| 303 | IF (iret /= NF90_NOERR) THEN |
---|
| 304 | CALL ipslerr (3,'forcesoil', & |
---|
| 305 | & 'Could not open file : ', & |
---|
| 306 | & Cforcing_name,'(Do you have forget it ?)') |
---|
| 307 | ENDIF |
---|
| 308 | !- |
---|
| 309 | ! Total Domain size is stored in nbp_glo variable |
---|
| 310 | !- |
---|
| 311 | ier = NF90_GET_ATT (Cforcing_id,NF90_GLOBAL,'kjpindex',x_tmp) |
---|
| 312 | nbp_glo = NINT(x_tmp) |
---|
| 313 | !- |
---|
| 314 | ! Number of values stored per year in the forcing file is stored in nparan var. |
---|
| 315 | !- |
---|
| 316 | ier = NF90_GET_ATT (Cforcing_id,NF90_GLOBAL,'nparan',x_tmp) |
---|
| 317 | nparan = NINT(x_tmp) |
---|
| 318 | ier = NF90_GET_ATT (Cforcing_id,NF90_GLOBAL,'nbyear',x_tmp) |
---|
| 319 | nbyear = NINT(x_tmp) |
---|
| 320 | !- |
---|
| 321 | ALLOCATE (indices_g(nbp_glo)) |
---|
| 322 | ALLOCATE (clay_g(nbp_glo)) |
---|
| 323 | !- |
---|
| 324 | ALLOCATE (x_indices_g(nbp_glo),stat=ier) |
---|
| 325 | ier = NF90_INQ_VARID (Cforcing_id,'index',v_id) |
---|
| 326 | ier = NF90_GET_VAR (Cforcing_id,v_id,x_indices_g) |
---|
| 327 | indices_g(:) = NINT(x_indices_g(:)) |
---|
| 328 | WRITE(numout,*) mpi_rank,"indices globaux : ",indices_g |
---|
| 329 | DEALLOCATE (x_indices_g) |
---|
| 330 | !- |
---|
| 331 | ier = NF90_INQ_VARID (Cforcing_id,'clay',v_id) |
---|
| 332 | ier = NF90_GET_VAR (Cforcing_id,v_id,clay_g) |
---|
| 333 | !- |
---|
| 334 | ! time step of forcesoil program (in days) |
---|
| 335 | !- |
---|
| 336 | dt_forcesoil = one_year / FLOAT(nparan) |
---|
| 337 | ! Initialize dt_sechiba in module constantes_var. This is needed in carbonsoil. |
---|
| 338 | dt_sechiba=dt_forcesoil*one_day |
---|
| 339 | WRITE(numout,*) 'time step (d): ',dt_forcesoil |
---|
| 340 | WRITE(numout,*) 'nparan: ',nparan |
---|
| 341 | WRITE(numout,*) 'nbyear: ',nbyear |
---|
| 342 | !- |
---|
| 343 | ! read and write the variables in the output restart file we do not modify within the Forcesoil program |
---|
| 344 | ! ie all variables stored in the input restart file except those stored in taboo_vars |
---|
| 345 | !- |
---|
| 346 | taboo_vars ='$lon$ $lat$ $lev$ $nav_lon$ $nav_lat$ $nav_lev$ $time$ $time_steps$ '// & |
---|
| 347 | & '$day_counter$ $dt_days$ $date$ ' |
---|
| 348 | !- |
---|
| 349 | DO m = 1,nvm |
---|
| 350 | WRITE(part_str,'(I2)') m |
---|
| 351 | IF (m < 10) part_str(1:1) = '0' |
---|
| 352 | temp_name = '$carbon_'//part_str(1:LEN_TRIM(part_str))//'$' |
---|
| 353 | taboo_vars = TRIM(taboo_vars)//' '//TRIM(temp_name) |
---|
| 354 | ENDDO |
---|
| 355 | !- |
---|
| 356 | CALL ioget_vname(rest_id_sto, nbvar, varnames) |
---|
| 357 | !- |
---|
| 358 | ! read and write some special variables (1D or variables that we need) |
---|
| 359 | !- |
---|
| 360 | var_name = 'day_counter' |
---|
| 361 | CALL restget (rest_id_sto, var_name, 1, 1, 1, itau_dep, .TRUE., xtmp) |
---|
| 362 | CALL restput (rest_id_sto, var_name, 1, 1, 1, itau_dep, xtmp) |
---|
| 363 | !- |
---|
| 364 | var_name = 'dt_days' |
---|
| 365 | CALL restget (rest_id_sto, var_name, 1, 1, 1, itau_dep, .TRUE., xtmp) |
---|
| 366 | CALL restput (rest_id_sto, var_name, 1, 1, 1, itau_dep, xtmp) |
---|
| 367 | !- |
---|
| 368 | var_name = 'date' |
---|
| 369 | CALL restget (rest_id_sto, var_name, 1, 1, 1, itau_dep, .TRUE., xtmp) |
---|
| 370 | CALL restput (rest_id_sto, var_name, 1, 1, 1, itau_dep, xtmp) |
---|
| 371 | !- |
---|
| 372 | DO iv=1,nbvar |
---|
| 373 | !-- check if the variable is to be written here |
---|
| 374 | IF (INDEX(taboo_vars,'$'//TRIM(varnames(iv))//'$') == 0 ) THEN |
---|
| 375 | !---- get variable dimensions, especially 3rd dimension |
---|
| 376 | CALL ioget_vdim & |
---|
| 377 | & (rest_id_sto, varnames(iv), varnbdim_max, varnbdim, vardims) |
---|
| 378 | l1d = ALL(vardims(1:varnbdim) == 1) |
---|
| 379 | !---- read it |
---|
| 380 | IF (l1d) THEN |
---|
| 381 | CALL restget & |
---|
| 382 | & (rest_id_sto, TRIM(varnames(iv)), 1, vardims(3), & |
---|
| 383 | & 1, itau_dep, .TRUE., xtmp) |
---|
| 384 | ELSE |
---|
| 385 | ALLOCATE( var_3d(nbp_glo,vardims(3)), stat=ier) |
---|
| 386 | IF (ier /= 0) STOP 'ALLOCATION PROBLEM' |
---|
| 387 | !---- |
---|
| 388 | CALL restget & |
---|
| 389 | & (rest_id_sto, TRIM(varnames(iv)), nbp_glo, vardims(3), & |
---|
| 390 | & 1, itau_dep, .TRUE., var_3d, "gather", nbp_glo, indices_g) |
---|
| 391 | ENDIF |
---|
| 392 | !---- write it |
---|
| 393 | IF (l1d) THEN |
---|
| 394 | CALL restput & |
---|
| 395 | & (rest_id_sto, TRIM(varnames(iv)), 1, vardims(3), & |
---|
| 396 | & 1, itau_dep, xtmp) |
---|
| 397 | ELSE |
---|
| 398 | CALL restput & |
---|
| 399 | & (rest_id_sto, TRIM(varnames(iv)), nbp_glo, vardims(3), & |
---|
| 400 | & 1, itau_dep, var_3d, 'scatter', nbp_glo, indices_g) |
---|
| 401 | !---- |
---|
| 402 | DEALLOCATE(var_3d) |
---|
| 403 | ENDIF |
---|
| 404 | ENDIF |
---|
| 405 | ENDDO |
---|
| 406 | !- |
---|
| 407 | ! read soil carbon stocks values stored in the input restart file |
---|
| 408 | !- |
---|
| 409 | ALLOCATE(carbon_g(nbp_glo,ncarb,nvm)) |
---|
| 410 | carbon_g(:,:,:) = val_exp |
---|
| 411 | DO m = 1, nvm |
---|
| 412 | WRITE (part_str, '(I2)') m |
---|
| 413 | IF (m<10) part_str(1:1)='0' |
---|
| 414 | var_name = 'carbon_'//part_str(1:LEN_TRIM(part_str)) |
---|
| 415 | CALL restget & |
---|
| 416 | & (rest_id_sto, var_name, nbp_glo, ncarb , 1, itau_dep, & |
---|
| 417 | & .TRUE., carbon_g(:,:,m), 'gather', nbp_glo, indices_g) |
---|
| 418 | IF (ALL(carbon_g(:,:,m) == val_exp)) carbon_g(:,:,m) = zero |
---|
| 419 | !-- do not write this variable: it will be modified. |
---|
| 420 | ENDDO |
---|
| 421 | WRITE(numout,*) "date0 : ",date0, itau_dep |
---|
| 422 | !- |
---|
| 423 | ! Analytical spinup is set to false |
---|
| 424 | ! |
---|
| 425 | spinup_analytic = .FALSE. |
---|
| 426 | |
---|
| 427 | ! Length of the run (in Years) |
---|
| 428 | ! User value is used if present in the .def file |
---|
| 429 | ! If not, default value (10000 Years) is used |
---|
| 430 | !- |
---|
| 431 | WRITE(time_str,'(a)') '10000Y' |
---|
| 432 | CALL getin('TIME_LENGTH', time_str) |
---|
| 433 | write(numout,*) 'Number of years for carbon spinup : ',time_str |
---|
| 434 | ! transform into itau |
---|
| 435 | CALL tlen2itau(time_str, dt_forcesoil*one_day, date0, itau_len) |
---|
| 436 | write(numout,*) 'Number of time steps to do: ',itau_len |
---|
| 437 | !- |
---|
| 438 | ! read soil carbon inputs, water and temperature stresses on OM decomposition |
---|
| 439 | ! into the forcing file - We read an average year. |
---|
| 440 | !- |
---|
| 441 | ALLOCATE(soilcarbon_input_g(nbp_glo,ncarb,nvm,nparan*nbyear)) |
---|
| 442 | ALLOCATE(control_temp_g(nbp_glo,nlevs,nparan*nbyear)) |
---|
| 443 | ALLOCATE(control_moist_g(nbp_glo,nlevs,nparan*nbyear)) |
---|
| 444 | !- |
---|
| 445 | ier = NF90_INQ_VARID (Cforcing_id,'soilcarbon_input',v_id) |
---|
| 446 | ier = NF90_GET_VAR (Cforcing_id,v_id,soilcarbon_input_g) |
---|
| 447 | ier = NF90_INQ_VARID (Cforcing_id, 'control_moist',v_id) |
---|
| 448 | ier = NF90_GET_VAR (Cforcing_id,v_id,control_moist_g) |
---|
| 449 | ier = NF90_INQ_VARID (Cforcing_id, 'control_temp',v_id) |
---|
| 450 | ier = NF90_GET_VAR (Cforcing_id,v_id,control_temp_g) |
---|
| 451 | !- |
---|
| 452 | ier = NF90_CLOSE (Cforcing_id) |
---|
| 453 | !- |
---|
| 454 | ENDIF |
---|
| 455 | CALL bcast(nparan) |
---|
| 456 | CALL bcast(nbyear) |
---|
| 457 | CALL bcast(dt_forcesoil) |
---|
| 458 | CALL bcast(iim_g) |
---|
| 459 | CALL bcast(jjm_g) |
---|
| 460 | CALL bcast(nbp_glo) |
---|
| 461 | CALL bcast(itau_dep) |
---|
| 462 | CALL bcast(itau_len) |
---|
| 463 | IF (.NOT. ALLOCATED(indices_g)) ALLOCATE (indices_g(nbp_glo)) |
---|
| 464 | CALL bcast(indices_g) |
---|
| 465 | |
---|
| 466 | ! |
---|
| 467 | ! We must initialize data_para : |
---|
| 468 | CALL init_orchidee_data_para_driver(nbp_glo,indices_g) |
---|
| 469 | |
---|
| 470 | kjpindex=nbp_loc |
---|
| 471 | jjm=jj_nb |
---|
| 472 | iim=iim_g |
---|
| 473 | IF (printlev_loc>=3) WRITE(numout,*) "Local grid : ",kjpindex,iim,jjm |
---|
| 474 | |
---|
| 475 | !--- |
---|
| 476 | !--- Create the index table |
---|
| 477 | !--- |
---|
| 478 | !--- This job returns a LOCAL kindex. |
---|
| 479 | !--- |
---|
| 480 | ALLOCATE (indices(kjpindex),stat=ier) |
---|
| 481 | ! |
---|
| 482 | !! scattering to all processes in parallel mode |
---|
| 483 | ! |
---|
| 484 | CALL scatter(indices_g,indices) |
---|
| 485 | indices(1:kjpindex)=indices(1:kjpindex)-(jj_begin-1)*iim_g |
---|
| 486 | IF (printlev_loc>=3) WRITE(numout,*) mpi_rank,"indices locaux = ",indices(1:kjpindex) |
---|
| 487 | !- |
---|
| 488 | ! Allocation of the variables for a processor |
---|
| 489 | !- |
---|
| 490 | ALLOCATE(clay(kjpindex)) |
---|
| 491 | ALLOCATE(soilcarbon_input(kjpindex,ncarb,nvm,nparan*nbyear)) |
---|
| 492 | ALLOCATE(control_temp(kjpindex,nlevs,nparan*nbyear)) |
---|
| 493 | ALLOCATE(control_moist(kjpindex,nlevs,nparan*nbyear)) |
---|
| 494 | ALLOCATE(carbon(kjpindex,ncarb,nvm)) |
---|
| 495 | ALLOCATE(resp_hetero_soil(kjpindex,nvm)) |
---|
| 496 | ALLOCATE(matrixA(kjpindex,nvm,nbpools,nbpools)) |
---|
| 497 | DO i = 1,nbpools |
---|
| 498 | matrixA(:,:,i,i) = un |
---|
| 499 | ENDDO |
---|
| 500 | iatt = 0 |
---|
| 501 | |
---|
| 502 | !- |
---|
| 503 | ! Initialization of the variables for a processor |
---|
| 504 | !- |
---|
| 505 | CALL Scatter(clay_g,clay) |
---|
| 506 | CALL Scatter(soilcarbon_input_g,soilcarbon_input) |
---|
| 507 | CALL Scatter(control_temp_g,control_temp) |
---|
| 508 | CALL Scatter(control_moist_g,control_moist) |
---|
| 509 | CALL Scatter(carbon_g,carbon) |
---|
| 510 | |
---|
| 511 | !- |
---|
| 512 | ! Configuration of the parameters |
---|
| 513 | !- |
---|
| 514 | |
---|
| 515 | !Config Key = FRAC_CARB_AP |
---|
| 516 | !Config Desc = frac carb coefficients from active pool: depends on clay content |
---|
| 517 | !Config if = OK_STOMATE |
---|
| 518 | !Config Def = 0.004 |
---|
| 519 | !Config Help = fraction of the active pool going to the passive pool |
---|
| 520 | !Config Units = [-] |
---|
| 521 | CALL getin_p('FRAC_CARB_AP',frac_carb_ap) |
---|
| 522 | ! |
---|
| 523 | !Config Key = FRAC_CARB_SA |
---|
| 524 | !Config Desc = frac_carb_coefficients from slow pool |
---|
| 525 | !Config if = OK_STOMATE |
---|
| 526 | !Config Def = 0.42 |
---|
| 527 | !Config Help = fraction of the slow pool going to the active pool |
---|
| 528 | !Config Units = [-] |
---|
| 529 | CALL getin_p('FRAC_CARB_SA',frac_carb_sa) |
---|
| 530 | ! |
---|
| 531 | !Config Key = FRAC_CARB_SP |
---|
| 532 | !Config Desc = frac_carb_coefficients from slow pool |
---|
| 533 | !Config if = OK_STOMATE |
---|
| 534 | !Config Def = 0.03 |
---|
| 535 | !Config Help = fraction of the slow pool going to the passive pool |
---|
| 536 | !Config Units = [-] |
---|
| 537 | CALL getin_p('FRAC_CARB_SP',frac_carb_sp) |
---|
| 538 | ! |
---|
| 539 | !Config Key = FRAC_CARB_PA |
---|
| 540 | !Config Desc = frac_carb_coefficients from passive pool |
---|
| 541 | !Config if = OK_STOMATE |
---|
| 542 | !Config Def = 0.45 |
---|
| 543 | !Config Help = fraction of the passive pool going to the passive pool |
---|
| 544 | !Config Units = [-] |
---|
| 545 | CALL getin_p('FRAC_CARB_PA',frac_carb_pa) |
---|
| 546 | ! |
---|
| 547 | !Config Key = FRAC_CARB_PS |
---|
| 548 | !Config Desc = frac_carb_coefficients from passive pool |
---|
| 549 | !Config if = OK_STOMATE |
---|
| 550 | !Config Def = 0.0 |
---|
| 551 | !Config Help = fraction of the passive pool going to the passive pool |
---|
| 552 | !Config Units = [-] |
---|
| 553 | CALL getin_p('FRAC_CARB_PS',frac_carb_ps) |
---|
| 554 | ! |
---|
| 555 | !Config Key = ACTIVE_TO_PASS_CLAY_FRAC |
---|
| 556 | !Config Desc = |
---|
| 557 | !Config if = OK_STOMATE |
---|
| 558 | !Config Def = .68 |
---|
| 559 | !Config Help = |
---|
| 560 | !Config Units = [-] |
---|
| 561 | CALL getin_p('ACTIVE_TO_PASS_CLAY_FRAC',active_to_pass_clay_frac) |
---|
| 562 | ! |
---|
| 563 | !Config Key = CARBON_TAU_IACTIVE |
---|
| 564 | !Config Desc = residence times in carbon pools |
---|
| 565 | !Config if = OK_STOMATE |
---|
| 566 | !Config Def = 0.149 |
---|
| 567 | !Config Help = |
---|
| 568 | !Config Units = [days] |
---|
| 569 | CALL getin_p('CARBON_TAU_IACTIVE',carbon_tau_iactive) |
---|
| 570 | ! |
---|
| 571 | !Config Key = CARBON_TAU_ISLOW |
---|
| 572 | !Config Desc = residence times in carbon pools |
---|
| 573 | !Config if = OK_STOMATE |
---|
| 574 | !Config Def = 5.48 |
---|
| 575 | !Config Help = |
---|
| 576 | !Config Units = [days] |
---|
| 577 | CALL getin_p('CARBON_TAU_ISLOW',carbon_tau_islow) |
---|
| 578 | ! |
---|
| 579 | !Config Key = CARBON_TAU_IPASSIVE |
---|
| 580 | !Config Desc = residence times in carbon pools |
---|
| 581 | !Config if = OK_STOMATE |
---|
| 582 | !Config Def = 241. |
---|
| 583 | !Config Help = |
---|
| 584 | !Config Units = [days] |
---|
| 585 | CALL getin_p('CARBON_TAU_IPASSIVE',carbon_tau_ipassive) |
---|
| 586 | ! |
---|
| 587 | !Config Key = FLUX_TOT_COEFF |
---|
| 588 | !Config Desc = |
---|
| 589 | !Config if = OK_STOMATE |
---|
| 590 | !Config Def = 1.2, 1.4,.75 |
---|
| 591 | !Config Help = |
---|
| 592 | !Config Units = [days] |
---|
| 593 | CALL getin_p('FLUX_TOT_COEFF',flux_tot_coeff) |
---|
| 594 | |
---|
| 595 | !!- |
---|
| 596 | !! 2. Computational step |
---|
| 597 | !! Loop over time - Call of soilcarbon routine at each time step |
---|
| 598 | !! Updated soil carbon stocks are stored into carbon variable |
---|
| 599 | !! We only keep the last value of carbon variable (no time dimension). |
---|
| 600 | !!- |
---|
| 601 | iyear=1 |
---|
| 602 | DO i=1,itau_len |
---|
| 603 | iatt = iatt+1 |
---|
| 604 | IF (iatt > nparan*nbyear) THEN |
---|
| 605 | IF (printlev_loc>=3) WRITE(numout,*) iyear |
---|
| 606 | iatt = 1 |
---|
| 607 | iyear=iyear+1 |
---|
| 608 | ENDIF |
---|
| 609 | CALL soilcarbon & |
---|
| 610 | & (kjpindex, clay, & |
---|
| 611 | & soilcarbon_input(:,:,:,iatt), & |
---|
| 612 | & control_temp(:,:,iatt), control_moist(:,:,iatt), & |
---|
| 613 | & carbon, resp_hetero_soil, & |
---|
| 614 | & matrixA) |
---|
| 615 | ENDDO |
---|
| 616 | WRITE(numout,*) "End of soilcarbon LOOP." |
---|
| 617 | |
---|
| 618 | ! |
---|
| 619 | !! Gathering of variables towards main processor in parallel mode |
---|
| 620 | ! |
---|
| 621 | CALL Gather(carbon,carbon_g) |
---|
| 622 | !!- |
---|
| 623 | !! 3. write new carbon stocks into the ouput restart file |
---|
| 624 | !!- |
---|
| 625 | IF (is_root_prc) THEN |
---|
| 626 | DO m=1,nvm |
---|
| 627 | WRITE (part_str, '(I2)') m |
---|
| 628 | IF (m<10) part_str(1:1)='0' |
---|
| 629 | var_name = 'carbon_'//part_str(1:LEN_TRIM(part_str)) |
---|
| 630 | CALL restput & |
---|
| 631 | & (rest_id_sto, var_name, nbp_glo, ncarb , 1, itau_dep, & |
---|
| 632 | & carbon_g(:,:,m), 'scatter', nbp_glo, indices_g) |
---|
| 633 | ENDDO |
---|
| 634 | !- |
---|
| 635 | CALL getin_dump |
---|
| 636 | CALL restclo |
---|
| 637 | ENDIF |
---|
| 638 | #ifdef CPP_PARA |
---|
| 639 | CALL MPI_FINALIZE(ier) |
---|
| 640 | #endif |
---|
| 641 | WRITE(numout,*) "End of forcesoil." |
---|
| 642 | !-------------------- |
---|
| 643 | END PROGRAM forcesoil |
---|