[65] | 1 | PROGRAM forcesoil |
---|
| 2 | !--------------------------------------------------------------------- |
---|
| 3 | ! This program reads a forcing file for STOMATE's soil carbon routine |
---|
| 4 | ! which was created with STOMATE. |
---|
| 5 | ! It then integrates STOMATE's soil carbon parameterizations |
---|
| 6 | ! for a given number of years using this forcing. |
---|
| 7 | !--------------------------------------------------------------------- |
---|
| 8 | !- IPSL (2006) |
---|
| 9 | !- This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC |
---|
| 10 | USE netcdf |
---|
| 11 | !- |
---|
| 12 | USE defprec |
---|
| 13 | USE constantes |
---|
| 14 | USE pft_parameters |
---|
| 15 | USE stomate_data |
---|
| 16 | USE ioipsl |
---|
| 17 | USE stomate_soilcarbon |
---|
| 18 | USE parallel |
---|
| 19 | !- |
---|
| 20 | IMPLICIT NONE |
---|
| 21 | !- |
---|
| 22 | CHARACTER(LEN=80) :: sto_restname_in,sto_restname_out,var_name |
---|
| 23 | INTEGER(i_std) :: iim,jjm |
---|
| 24 | INTEGER(i_std),PARAMETER :: llm = 1 |
---|
| 25 | INTEGER(i_std) :: kjpindex |
---|
| 26 | INTEGER(i_std) :: itau_dep,itau_len |
---|
| 27 | CHARACTER(LEN=30) :: time_str |
---|
| 28 | INTEGER(i_std) :: ier,iret |
---|
| 29 | REAL(r_std) :: dt_files |
---|
| 30 | REAL(r_std) :: date0 |
---|
| 31 | INTEGER(i_std) :: rest_id_sto |
---|
| 32 | INTEGER(i_std) :: ncfid |
---|
| 33 | REAL(r_std) :: dt_force,dt_forcesoil |
---|
| 34 | INTEGER :: nparan |
---|
| 35 | INTEGER,PARAMETER :: nparanmax=36 |
---|
| 36 | REAL(r_std) :: xbid1,xbid2 |
---|
| 37 | INTEGER(i_std) :: ibid |
---|
| 38 | INTEGER(i_std),DIMENSION(:),ALLOCATABLE :: indices |
---|
| 39 | REAL(r_std),DIMENSION(llm) :: lev |
---|
| 40 | REAL(r_std),DIMENSION(:,:,:,:),ALLOCATABLE :: soilcarbon_input |
---|
| 41 | REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: & |
---|
| 42 | & carbon,control_moist,control_temp |
---|
| 43 | REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: & |
---|
| 44 | & lon,lat,resp_hetero_soil,var_3d |
---|
| 45 | REAL(r_std),DIMENSION(:),ALLOCATABLE :: & |
---|
| 46 | & x_indices |
---|
| 47 | REAL(r_std) :: time |
---|
| 48 | INTEGER :: i,j,m,iatt,iv |
---|
| 49 | CHARACTER(LEN=400) :: taboo_vars |
---|
| 50 | REAL(r_std),DIMENSION(1) :: xtmp |
---|
| 51 | INTEGER(i_std),PARAMETER :: nbvarmax=300 |
---|
| 52 | INTEGER(i_std) :: nbvar |
---|
| 53 | CHARACTER(LEN=50),DIMENSION(nbvarmax) :: varnames |
---|
| 54 | INTEGER(i_std) :: varnbdim |
---|
| 55 | INTEGER(i_std),PARAMETER :: varnbdim_max=20 |
---|
| 56 | INTEGER,DIMENSION(varnbdim_max) :: vardims |
---|
| 57 | LOGICAL :: l1d |
---|
| 58 | REAL(r_std) :: x_tmp |
---|
| 59 | ! clay fraction |
---|
| 60 | REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:) :: clay |
---|
| 61 | !- |
---|
| 62 | ! string suffix indicating an index |
---|
| 63 | CHARACTER(LEN=10) :: part_str |
---|
| 64 | ! |
---|
| 65 | CHARACTER(LEN=100) :: Cforcing_name |
---|
| 66 | INTEGER :: Cforcing_id |
---|
| 67 | INTEGER :: v_id |
---|
| 68 | |
---|
| 69 | REAL(r_std),ALLOCATABLE :: clay_loc(:) |
---|
| 70 | REAL(r_std),ALLOCATABLE :: soilcarbon_input_loc(:,:,:,:) |
---|
| 71 | REAL(r_std),ALLOCATABLE :: control_temp_loc(:,:,:) |
---|
| 72 | REAL(r_std),ALLOCATABLE :: control_moist_loc(:,:,:) |
---|
| 73 | REAL(r_std),ALLOCATABLE :: carbon_loc(:,:,:) |
---|
| 74 | INTEGER :: ierr |
---|
| 75 | |
---|
| 76 | CALL Init_para(.FALSE.) |
---|
| 77 | |
---|
| 78 | CALL getin('NVM',nvm) |
---|
| 79 | |
---|
| 80 | !- |
---|
| 81 | ! Stomate's restart files |
---|
| 82 | !- |
---|
| 83 | IF (is_root_prc) THEN |
---|
| 84 | sto_restname_in = 'stomate_start.nc' |
---|
| 85 | CALL getin ('STOMATE_RESTART_FILEIN',sto_restname_in) |
---|
| 86 | WRITE(*,*) 'STOMATE INPUT RESTART_FILE: ',TRIM(sto_restname_in) |
---|
| 87 | sto_restname_out = 'stomate_restart.nc' |
---|
| 88 | CALL getin ('STOMATE_RESTART_FILEOUT',sto_restname_out) |
---|
| 89 | WRITE(*,*) 'STOMATE OUTPUT RESTART_FILE: ',TRIM(sto_restname_out) |
---|
| 90 | !- |
---|
| 91 | ! We need to know iim, jjm. |
---|
| 92 | ! Get them from the restart files themselves. |
---|
| 93 | !- |
---|
| 94 | iret = NF90_OPEN (sto_restname_in, NF90_NOWRITE, ncfid) |
---|
| 95 | iret = NF90_INQUIRE_DIMENSION (ncfid,1,len=iim) |
---|
| 96 | iret = NF90_INQUIRE_DIMENSION (ncfid,2,len=jjm) |
---|
| 97 | iret = NF90_CLOSE (ncfid) |
---|
| 98 | !- |
---|
| 99 | ! Allocate longitudes and latitudes |
---|
| 100 | !- |
---|
| 101 | ALLOCATE (lon(iim,jjm)) |
---|
| 102 | ALLOCATE (lat(iim,jjm)) |
---|
| 103 | lon(:,:) = 0.0 |
---|
| 104 | lat(:,:) = 0.0 |
---|
| 105 | lev(1) = 0.0 |
---|
| 106 | !- |
---|
| 107 | CALL restini & |
---|
| 108 | & (sto_restname_in, iim, jjm, lon, lat, llm, lev, & |
---|
| 109 | & sto_restname_out, itau_dep, date0, dt_files, rest_id_sto) |
---|
| 110 | ENDIF |
---|
| 111 | !- |
---|
| 112 | ! calendar |
---|
| 113 | !- |
---|
| 114 | CALL bcast(date0) |
---|
| 115 | !!! MM : Ã revoir : choix du calendrier dans forcesoil ?? Il est dans le restart de stomate ! |
---|
| 116 | ! CALL ioconf_calendar ('noleap') |
---|
| 117 | CALL ioget_calendar (one_year,one_day) |
---|
| 118 | |
---|
| 119 | CALL ioconf_startdate(date0) |
---|
| 120 | |
---|
| 121 | IF (is_root_prc) THEN |
---|
| 122 | !- |
---|
| 123 | ! open FORCESOIL's forcing file to read some basic info |
---|
| 124 | !- |
---|
| 125 | Cforcing_name = 'stomate_Cforcing.nc' |
---|
| 126 | CALL getin ('STOMATE_CFORCING_NAME',Cforcing_name) |
---|
| 127 | !- |
---|
| 128 | ier = NF90_OPEN (TRIM(Cforcing_name),NF90_NOWRITE,Cforcing_id) |
---|
| 129 | !- |
---|
| 130 | ier = NF90_GET_ATT (Cforcing_id,NF90_GLOBAL,'kjpindex',x_tmp) |
---|
| 131 | kjpindex = NINT(x_tmp) |
---|
| 132 | ier = NF90_GET_ATT (Cforcing_id,NF90_GLOBAL,'nparan',x_tmp) |
---|
| 133 | nparan = NINT(x_tmp) |
---|
| 134 | !- |
---|
| 135 | ALLOCATE (indices(kjpindex)) |
---|
| 136 | ALLOCATE (clay(kjpindex)) |
---|
| 137 | !- |
---|
| 138 | ALLOCATE (x_indices(kjpindex),stat=ier) |
---|
| 139 | ier = NF90_INQ_VARID (Cforcing_id,'index',v_id) |
---|
| 140 | ier = NF90_GET_VAR (Cforcing_id,v_id,x_indices) |
---|
| 141 | indices(:) = NINT(x_indices(:)) |
---|
| 142 | DEALLOCATE (x_indices) |
---|
| 143 | !- |
---|
| 144 | ier = NF90_INQ_VARID (Cforcing_id,'clay',v_id) |
---|
| 145 | ier = NF90_GET_VAR (Cforcing_id,v_id,clay) |
---|
| 146 | !- |
---|
| 147 | ! time step of forcesoil |
---|
| 148 | !- |
---|
| 149 | dt_forcesoil = one_year / FLOAT(nparan) |
---|
| 150 | WRITE(*,*) 'time step (d): ',dt_forcesoil |
---|
| 151 | !- |
---|
| 152 | ! read (and partially write) the restart file |
---|
| 153 | !- |
---|
| 154 | ! read and write the variables we do not need |
---|
| 155 | !- |
---|
| 156 | taboo_vars = '$lon$ $lat$ $lev$ $nav_lon$ $nav_lat$ $nav_lev$ $time$ $time_steps$ '// & |
---|
| 157 | & '$day_counter$ $dt_days$ $date$ '// & |
---|
| 158 | & '$carbon_01$ $carbon_02$ $carbon_03$ $carbon_04$ $carbon_05$'// & |
---|
| 159 | & '$carbon_06$ $carbon_07$ $carbon_08$ $carbon_09$ $carbon_10$'// & |
---|
| 160 | & '$carbon_11$ $carbon_12$ $carbon_13$' |
---|
| 161 | !- |
---|
| 162 | CALL ioget_vname(rest_id_sto, nbvar, varnames) |
---|
| 163 | !- |
---|
| 164 | ! read and write some special variables (1D or variables that we need) |
---|
| 165 | !- |
---|
| 166 | var_name = 'day_counter' |
---|
| 167 | CALL restget (rest_id_sto, var_name, 1, 1, 1, itau_dep, .TRUE., xtmp) |
---|
| 168 | CALL restput (rest_id_sto, var_name, 1, 1, 1, itau_dep, xtmp) |
---|
| 169 | !- |
---|
| 170 | var_name = 'dt_days' |
---|
| 171 | CALL restget (rest_id_sto, var_name, 1, 1, 1, itau_dep, .TRUE., xtmp) |
---|
| 172 | CALL restput (rest_id_sto, var_name, 1, 1, 1, itau_dep, xtmp) |
---|
| 173 | !- |
---|
| 174 | var_name = 'date' |
---|
| 175 | CALL restget (rest_id_sto, var_name, 1, 1, 1, itau_dep, .TRUE., xtmp) |
---|
| 176 | CALL restput (rest_id_sto, var_name, 1, 1, 1, itau_dep, xtmp) |
---|
| 177 | !- |
---|
| 178 | DO iv=1,nbvar |
---|
| 179 | !-- check if the variable is to be written here |
---|
| 180 | IF (INDEX(taboo_vars,'$'//TRIM(varnames(iv))//'$') == 0 ) THEN |
---|
| 181 | !---- get variable dimensions, especially 3rd dimension |
---|
| 182 | CALL ioget_vdim & |
---|
| 183 | & (rest_id_sto, varnames(iv), varnbdim_max, varnbdim, vardims) |
---|
| 184 | l1d = ALL(vardims(1:varnbdim) == 1) |
---|
| 185 | !---- |
---|
| 186 | ALLOCATE( var_3d(kjpindex,vardims(3)), stat=ier) |
---|
| 187 | IF (ier /= 0) STOP 'ALLOCATION PROBLEM' |
---|
| 188 | !---- read it |
---|
| 189 | IF (l1d) THEN |
---|
| 190 | CALL restget & |
---|
| 191 | & (rest_id_sto, TRIM(varnames(iv)), 1, vardims(3), & |
---|
| 192 | & 1, itau_dep, .TRUE., xtmp) |
---|
| 193 | ELSE |
---|
| 194 | CALL restget & |
---|
| 195 | & (rest_id_sto, TRIM(varnames(iv)), kjpindex, vardims(3), & |
---|
| 196 | & 1, itau_dep, .TRUE., var_3d, "gather", kjpindex, indices) |
---|
| 197 | ENDIF |
---|
| 198 | !---- write it |
---|
| 199 | IF (l1d) THEN |
---|
| 200 | CALL restput & |
---|
| 201 | & (rest_id_sto, TRIM(varnames(iv)), 1, vardims(3), & |
---|
| 202 | & 1, itau_dep, xtmp) |
---|
| 203 | ELSE |
---|
| 204 | CALL restput & |
---|
| 205 | & (rest_id_sto, TRIM(varnames(iv)), kjpindex, vardims(3), & |
---|
| 206 | & 1, itau_dep, var_3d, 'scatter', kjpindex, indices) |
---|
| 207 | ENDIF |
---|
| 208 | !---- |
---|
| 209 | DEALLOCATE(var_3d) |
---|
| 210 | ENDIF |
---|
| 211 | ENDDO |
---|
| 212 | !- |
---|
| 213 | ! read soil carbon |
---|
| 214 | !- |
---|
| 215 | ALLOCATE(carbon(kjpindex,ncarb,nvm)) |
---|
| 216 | carbon(:,:,:) = val_exp |
---|
| 217 | DO m = 1, nvm |
---|
| 218 | WRITE (part_str, '(I2)') m |
---|
| 219 | IF (m<10) part_str(1:1)='0' |
---|
| 220 | var_name = 'carbon_'//part_str(1:LEN_TRIM(part_str)) |
---|
| 221 | CALL restget & |
---|
| 222 | & (rest_id_sto, var_name, kjpindex, ncarb , 1, itau_dep, & |
---|
| 223 | & .TRUE., carbon(:,:,m), 'gather', kjpindex, indices) |
---|
| 224 | IF (ALL(carbon(:,:,m) == val_exp)) carbon(:,:,m) = zero |
---|
| 225 | !-- do not write this variable: it will be modified. |
---|
| 226 | ENDDO |
---|
| 227 | !- |
---|
| 228 | ! Length of run |
---|
| 229 | !- |
---|
| 230 | WRITE(time_str,'(a)') '10000Y' |
---|
| 231 | CALL getin('TIME_LENGTH', time_str) |
---|
| 232 | ! transform into itau |
---|
| 233 | CALL tlen2itau(time_str, dt_forcesoil*one_year, date0, itau_len) |
---|
| 234 | write(*,*) 'Number of time steps to do: ',itau_len |
---|
| 235 | !- |
---|
| 236 | ! read the rest of the forcing file and store forcing in an array. |
---|
| 237 | ! We read an average year. |
---|
| 238 | !- |
---|
| 239 | ALLOCATE(soilcarbon_input(kjpindex,ncarb,nvm,nparan)) |
---|
| 240 | ALLOCATE(control_temp(kjpindex,nlevs,nparan)) |
---|
| 241 | ALLOCATE(control_moist(kjpindex,nlevs,nparan)) |
---|
| 242 | !- |
---|
| 243 | ier = NF90_INQ_VARID (Cforcing_id,'soilcarbon_input',v_id) |
---|
| 244 | ier = NF90_GET_VAR (Cforcing_id,v_id,soilcarbon_input) |
---|
| 245 | ier = NF90_INQ_VARID (Cforcing_id, 'control_moist',v_id) |
---|
| 246 | ier = NF90_GET_VAR (Cforcing_id,v_id,control_moist) |
---|
| 247 | ier = NF90_INQ_VARID (Cforcing_id, 'control_temp',v_id) |
---|
| 248 | ier = NF90_GET_VAR (Cforcing_id,v_id,control_temp) |
---|
| 249 | !- |
---|
| 250 | ier = NF90_CLOSE (Cforcing_id) |
---|
| 251 | !- |
---|
| 252 | !MM Problem here with dpu which depends on soil type |
---|
| 253 | DO iv = 1, nbdl-1 |
---|
| 254 | ! first 2.0 is dpu |
---|
| 255 | ! second 2.0 is average |
---|
| 256 | diaglev(iv) = 2.0/(2**(nbdl-1) -1) * ( ( 2**(iv-1) -1) + ( 2**(iv) -1) ) / 2.0 |
---|
| 257 | ENDDO |
---|
| 258 | diaglev(nbdl) = 2.0 |
---|
| 259 | !- |
---|
| 260 | ! For sequential use only, we must initialize data_para : |
---|
| 261 | ! |
---|
| 262 | ! |
---|
| 263 | ENDIF |
---|
| 264 | |
---|
| 265 | CALL bcast(iim) |
---|
| 266 | CALL bcast(jjm) |
---|
| 267 | call bcast(kjpindex) |
---|
| 268 | CALL init_data_para(iim,jjm,kjpindex,indices) |
---|
| 269 | |
---|
| 270 | !- |
---|
| 271 | !- |
---|
| 272 | ! there we go: time loop |
---|
| 273 | !- |
---|
| 274 | CALL bcast(nparan) |
---|
| 275 | ALLOCATE(clay_loc(nbp_loc)) |
---|
| 276 | ALLOCATE(soilcarbon_input_loc(nbp_loc,ncarb,nvm,nparan)) |
---|
| 277 | ALLOCATE(control_temp_loc(nbp_loc,nlevs,nparan)) |
---|
| 278 | ALLOCATE(control_moist_loc(nbp_loc,nlevs,nparan)) |
---|
| 279 | ALLOCATE(carbon_loc(nbp_loc,ncarb,nvm)) |
---|
| 280 | ALLOCATE(resp_hetero_soil(nbp_loc,nvm)) |
---|
| 281 | iatt = 0 |
---|
| 282 | |
---|
| 283 | CALL bcast(itau_len) |
---|
| 284 | CALL bcast(nparan) |
---|
| 285 | CALL bcast(dt_forcesoil) |
---|
| 286 | CALL Scatter(clay,clay_loc) |
---|
| 287 | CALL Scatter(soilcarbon_input,soilcarbon_input_loc) |
---|
| 288 | CALL Scatter(control_temp,control_temp_loc) |
---|
| 289 | CALL Scatter(control_moist,control_moist_loc) |
---|
| 290 | CALL Scatter(carbon,carbon_loc) |
---|
| 291 | |
---|
| 292 | DO i=1,itau_len |
---|
| 293 | iatt = iatt+1 |
---|
| 294 | IF (iatt > nparan) iatt = 1 |
---|
| 295 | |
---|
| 296 | !!$ DS : calling the new_values of soilcarbon parameters |
---|
| 297 | ! |
---|
| 298 | CALL getin('FRAC_CARB_AA',frac_carb_aa) |
---|
| 299 | CALL getin('FRAC_CARB_AP',frac_carb_ap) |
---|
| 300 | CALL getin('FRAC_CARB_SS',frac_carb_ss) |
---|
| 301 | CALL getin('FRAC_CARB_SA',frac_carb_sa) |
---|
| 302 | CALL getin('FRAC_CARB_SP',frac_carb_sp) |
---|
| 303 | CALL getin('FRAC_CARB_PP',frac_carb_pp) |
---|
| 304 | CALL getin('FRAC_CARB_PA',frac_carb_pa) |
---|
| 305 | CALL getin('FRAC_CARB_PS',frac_carb_ps) |
---|
| 306 | ! |
---|
| 307 | CALL getin('ACTIVE_TO_PASS_CLAY_FRAC',active_to_pass_clay_frac) |
---|
| 308 | CALL getin('CARBON_TAU_IACTIVE',carbon_tau_iactive) |
---|
| 309 | CALL getin('CARBON_TAU_ISLOW',carbon_tau_islow) |
---|
| 310 | CALL getin('CARBON_TAU_IPASSIVE',carbon_tau_ipassive) |
---|
| 311 | CALL getin('FLUX_TOT_COEFF',flux_tot_coeff) |
---|
| 312 | |
---|
| 313 | CALL soilcarbon & |
---|
| 314 | & (nbp_loc, dt_forcesoil, clay_loc, & |
---|
| 315 | & soilcarbon_input_loc(:,:,:,iatt), & |
---|
| 316 | & control_temp_loc(:,:,iatt), control_moist_loc(:,:,iatt), & |
---|
| 317 | & carbon_loc, resp_hetero_soil) |
---|
| 318 | ENDDO |
---|
| 319 | |
---|
| 320 | CALL Gather(carbon_loc,carbon) |
---|
| 321 | !- |
---|
| 322 | ! write new carbon into restart file |
---|
| 323 | !- |
---|
| 324 | IF (is_root_prc) THEN |
---|
| 325 | DO m=1,nvm |
---|
| 326 | WRITE (part_str, '(I2)') m |
---|
| 327 | IF (m<10) part_str(1:1)='0' |
---|
| 328 | var_name = 'carbon_'//part_str(1:LEN_TRIM(part_str)) |
---|
| 329 | CALL restput & |
---|
| 330 | & (rest_id_sto, var_name, kjpindex, ncarb , 1, itau_dep, & |
---|
| 331 | & carbon(:,:,m), 'scatter', kjpindex, indices) |
---|
| 332 | ENDDO |
---|
| 333 | !- |
---|
| 334 | CALL getin_dump |
---|
| 335 | CALL restclo |
---|
| 336 | ENDIF |
---|
| 337 | #ifdef CPP_PARA |
---|
| 338 | CALL MPI_FINALIZE(ierr) |
---|
| 339 | #endif |
---|
| 340 | !-------------------- |
---|
| 341 | END PROGRAM forcesoil |
---|