Changeset 367 for branches/ORCHIDEE_EXT/ORCHIDEE_OL/teststomate.f90
- Timestamp:
- 2011-08-01T11:25:14+02:00 (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/ORCHIDEE_EXT/ORCHIDEE_OL/teststomate.f90
r327 r367 21 21 USE slowproc 22 22 USE stomate 23 USE intersurf, ONLY: stom_define_history , intsurf_time23 USE intersurf, ONLY: stom_define_history, stom_ipcc_define_history, intsurf_time, l_first_intersurf, check_time 24 24 USE parallel 25 25 !- … … 28 28 ! Declarations 29 29 !- 30 INTEGER(i_std) :: vegax_id31 30 INTEGER(i_std) :: kjpij,kjpindex 32 31 REAL(r_std) :: dtradia,dt_force 32 33 33 INTEGER(i_std),DIMENSION(:),ALLOCATABLE :: indices 34 34 INTEGER(i_std),DIMENSION(:),ALLOCATABLE :: indexveg … … 41 41 REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: veget_max_force_x 42 42 REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: lai_force_x 43 REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: lon,lat44 43 REAL(r_std),DIMENSION(:),ALLOCATABLE :: t2m,t2m_min,temp_sol 45 44 REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: soiltemp,soilhum … … 53 52 REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: qsintmax_x 54 53 REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: co2_flux 55 !!$ >> DS add for externalised version 54 REAL(r_std),DIMENSION(:),ALLOCATABLE :: fco2_lu 55 56 INTEGER(i_std),DIMENSION(:),ALLOCATABLE :: indices_g 57 REAL(r_std),DIMENSION(:),ALLOCATABLE :: x_indices_g 58 REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: x_neighbours_g 59 60 !!$ >> DS for the externalized version 56 61 REAL(r_std),DIMENSION(:),ALLOCATABLE :: hist_PFTaxis 57 62 !!$ >> DS 58 !!$ >> Add for correct 1.9.5.1 version 01/07/201159 ! TYPE(control_type), SAVE :: control_flags60 REAL(r_std),DIMENSION(:),ALLOCATABLE :: fco2_lu61 !!$>> DS62 63 63 64 64 INTEGER :: ier,iret … … 68 68 & dri_restname_in,dri_restname_out, & 69 69 & sec_restname_in,sec_restname_out, & 70 & sto_restname_in,sto_restname_out, stom_histname 71 INTEGER(i_std) :: iim,jjm 72 INTEGER(i_std),PARAMETER :: llm = 1 73 REAL(r_std),DIMENSION(llm) :: lev 74 REAL(r_std) :: dt_files 70 & sto_restname_in,sto_restname_out, & 71 & stom_histname, stom_ipcc_histname 72 INTEGER(i_std) :: iim,jjm,llm 73 REAL, ALLOCATABLE, DIMENSION(:,:) :: lon, lat 74 REAL, ALLOCATABLE, DIMENSION(:) :: lev 75 LOGICAL :: rectilinear 76 REAL, ALLOCATABLE, DIMENSION(:) :: lon_rect, lat_rect 77 REAL(r_std) :: dt 75 78 INTEGER(i_std) :: itau_dep,itau,itau_len,itau_step 76 79 REAL(r_std) :: date0 77 80 INTEGER(i_std) :: rest_id_sec,rest_id_sto 78 INTEGER(i_std) :: hist_id_sec,hist_id_sec2,hist_id_sto ,hist_id_stom_IPCC81 INTEGER(i_std) :: hist_id_sec,hist_id_sec2,hist_id_stom,hist_id_stom_IPCC 79 82 CHARACTER(LEN=30) :: time_str 80 REAL :: hist_days_stom,hist_dt_stom 83 REAL(r_std) :: dt_slow_ 84 REAL :: hist_days_stom,hist_days_stom_ipcc,hist_dt_stom,hist_dt_stom_ipcc 81 85 REAL(r_std),DIMENSION(10) :: hist_pool_10axis 82 86 REAL(r_std),DIMENSION(100) :: hist_pool_100axis 83 87 REAL(r_std),DIMENSION(11) :: hist_pool_11axis 84 88 REAL(r_std),DIMENSION(101) :: hist_pool_101axis 85 INTEGER :: hist_PFTaxis_id,h ori_id89 INTEGER :: hist_PFTaxis_id,hist_IPCC_PFTaxis_id,hori_id 86 90 INTEGER :: hist_pool_10axis_id 87 91 INTEGER :: hist_pool_100axis_id 88 92 INTEGER :: hist_pool_11axis_id 89 93 INTEGER :: hist_pool_101axis_id 90 INTEGER :: hist_level 91 INTEGER,PARAMETER :: max_hist_level = 10 92 INTEGER(i_std) :: i,j,iv,id 93 CHARACTER*80 :: var_name 94 CHARACTER(LEN=40),DIMENSION(10) :: fluxop 94 INTEGER(i_std) :: i,j,iv 95 95 LOGICAL :: ldrestart_read,ldrestart_write 96 96 LOGICAL :: l1d … … 104 104 REAL(r_std),DIMENSION(1) :: xtmp 105 105 INTEGER :: nsfm,nsft 106 INTEGER :: iisf 107 INTEGER(i_std) :: max_totsize,totsize_1step 108 INTEGER(i_std),DIMENSION(0:2) :: ifirst,ilast 109 INTEGER(i_std) :: iblocks,nblocks 110 INTEGER,PARAMETER :: ndm = 10 111 INTEGER,DIMENSION(ndm) :: start,count 112 INTEGER :: ndim,v_id 113 INTEGER :: force_id 106 INTEGER :: iisf,iiisf 107 INTEGER(i_std) :: max_totsize,totsize_1step,totsize_tmp 108 109 INTEGER :: vid 114 110 CHARACTER(LEN=100) :: forcing_name 115 111 REAL :: x 116 REAL(r_std),DIMENSION(:),ALLOCATABLE :: x_indices 117 REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: x_neighbours 112 118 113 REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: var_3d 114 119 115 !- 120 116 REAL(r_std) :: time_sec,time_step_sec … … 122 118 REAL(r_std),DIMENSION(1) :: r1d 123 119 REAL(r_std) :: julian,djulian 124 ! REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: soiltype 125 !- 126 ! the following variables contain the forcing data 127 !- 128 REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: clay_fm 129 REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: humrel_x_fm 130 REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: litterhum_fm 131 REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: t2m_fm 132 REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: t2m_min_fm 133 REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: temp_sol_fm 134 REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: soiltemp_fm 135 REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: soilhum_fm 136 REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: precip_fm 137 REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: gpp_x_fm 138 REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: veget_force_x_fm 139 REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: veget_max_force_x_fm 140 REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: lai_force_x_fm 141 INTEGER(i_std),ALLOCATABLE,SAVE,DIMENSION(:) :: isf 142 LOGICAL,ALLOCATABLE,SAVE,DIMENSION(:) :: nf_written 143 INTEGER(i_std),ALLOCATABLE,SAVE,DIMENSION(:) :: nf_cumul 144 INTEGER(i_std) :: ji,jv 120 121 INTEGER(i_std) :: ji,jv,l 122 123 LOGICAL :: debug 124 145 125 !--------------------------------------------------------------------- 146 !- No parallelisation yet in teststomate ! 147 #ifdef CPP_PARA 148 CALL ipslerr (3,'teststomate', & 149 & 'You try to run testsomate compiled with parallelisation. (CPP_PARA key)', & 150 & 'But it wasn''t programmed yet and teststomate will stop.','You must compiled it without CPP_PARA key.') 151 #endif 152 126 127 CALL init_para(.FALSE.) 128 CALL init_timer 153 129 !- 154 130 ! calendar … … 156 132 CALL ioconf_calendar ('noleap') 157 133 CALL ioget_calendar (one_year,one_day) 158 !- 159 ! open STOMATE's forcing file to read some basic info 160 !- 161 forcing_name = 'stomate_forcing.nc' 162 CALL getin ('STOMATE_FORCING_NAME',forcing_name) 163 iret = NF90_OPEN (TRIM(forcing_name),NF90_NOWRITE,force_id) 164 IF (iret /= NF90_NOERR) THEN 165 CALL ipslerr (3,'teststomate', & 166 & 'Could not open file : ', & 167 & forcing_name,'(Do you have forget it ?)') 168 ENDIF 169 ier = NF90_GET_ATT (force_id,NF90_GLOBAL,'dtradia',dtradia) 170 ier = NF90_GET_ATT (force_id,NF90_GLOBAL,'dt_slow',dt_force) 171 ier = NF90_GET_ATT (force_id,NF90_GLOBAL,'nsft',x) 172 nsft = NINT(x) 173 ier = NF90_GET_ATT (force_id,NF90_GLOBAL,'kjpij',x) 174 kjpij = NINT(x) 175 ier = NF90_GET_ATT (force_id,NF90_GLOBAL,'kjpindex',x) 176 kjpindex = NINT(x) 177 CALL init_grid( kjpindex ) 178 !- 179 write(*,*) 'ATTENTION',dtradia,dt_force 134 135 IF (is_root_prc) THEN 136 !- 137 ! open STOMATE's forcing file to read some basic info 138 !- 139 forcing_name = 'stomate_forcing.nc' 140 CALL getin ('STOMATE_FORCING_NAME',forcing_name) 141 iret = NF90_OPEN (TRIM(forcing_name),NF90_NOWRITE,forcing_id) 142 IF (iret /= NF90_NOERR) THEN 143 CALL ipslerr (3,'teststomate', & 144 & 'Could not open file : ', & 145 & forcing_name,'(Do you have forget it ?)') 146 ENDIF 147 ier = NF90_GET_ATT (forcing_id,NF90_GLOBAL,'dtradia',dtradia) 148 ier = NF90_GET_ATT (forcing_id,NF90_GLOBAL,'dt_slow',dt_force) 149 ier = NF90_GET_ATT (forcing_id,NF90_GLOBAL,'nsft',x) 150 nsft = NINT(x) 151 ier = NF90_GET_ATT (forcing_id,NF90_GLOBAL,'kjpij',x) 152 kjpij = NINT(x) 153 ier = NF90_GET_ATT (forcing_id,NF90_GLOBAL,'kjpindex',x) 154 nbp_glo = NINT(x) 155 ENDIF 156 CALL bcast(dtradia) 157 CALL bcast(dt_force) 158 CALL bcast(nsft) 159 CALL bcast(nbp_glo) 160 !- 161 write(numout,*) 'ATTENTION',dtradia,dt_force 162 !- 163 ! read info about land points 164 !- 165 IF (is_root_prc) THEN 166 a_er=.FALSE. 167 ALLOCATE (indices_g(nbp_glo),stat=ier) 168 a_er = a_er .OR. (ier.NE.0) 169 IF (a_er) THEN 170 CALL ipslerr (3,'teststomate', & 171 & 'PROBLEM WITH ALLOCATION', & 172 & 'for local variables 1','') 173 ENDIF 174 ! 175 ALLOCATE (x_indices_g(nbp_glo),stat=ier) 176 a_er = a_er .OR. (ier.NE.0) 177 IF (a_er) THEN 178 CALL ipslerr (3,'teststomate', & 179 & 'PROBLEM WITH ALLOCATION', & 180 & 'for global variables 1','') 181 ENDIF 182 ier = NF90_INQ_VARID (forcing_id,'index',vid) 183 IF (ier .NE. 0) THEN 184 CALL ipslerr (3,'teststomate', & 185 & 'PROBLEM WITH ALLOCATION', & 186 & 'for global variables 1','') 187 ENDIF 188 ier = NF90_GET_VAR (forcing_id,vid,x_indices_g) 189 IF (iret /= NF90_NOERR) THEN 190 CALL ipslerr (3,'teststomate', & 191 & 'PROBLEM WITH variable "index" in file ', & 192 & forcing_name,'(check this file)') 193 ENDIF 194 indices_g(:) = NINT(x_indices_g(:)) 195 DEALLOCATE (x_indices_g) 196 ELSE 197 ALLOCATE (indices_g(0)) 198 ENDIF 199 !--------------------------------------------------------------------- 200 !- 201 ! set debug to have more information 202 !- 203 !Config Key = DEBUG_INFO 204 !Config Desc = Flag for debug information 205 !Config Def = n 206 !Config Help = This option allows to switch on the output of debug 207 !Config information without recompiling the code. 208 !- 209 debug = .FALSE. 210 CALL getin_p('DEBUG_INFO',debug) 211 ! 212 !Config Key = LONGPRINT 213 !Config Desc = ORCHIDEE will print more messages 214 !Config Def = n 215 !Config Help = This flag permits to print more debug messages in the run. 216 ! 217 long_print = .FALSE. 218 CALL getin_p('LONGPRINT',long_print) 180 219 181 220 !!$ >> DS added for the externalised version … … 183 222 184 223 ! 1. Read the number of PFTs 185 CALL getin('NVM',nvm) 224 ! 225 !Config Key = NVM 226 !Config Desc = number of PFTs 227 !Config if = ANYTIME 228 !Config Def = 13 229 !Config Help = The number of vegetation types define by the user 230 !Config Units = NONE 231 CALL getin_p('NVM',nvm) 232 186 233 ! 2. Allocate and read the pft parameters stomate specific 187 234 CALL pft_parameters_main 235 !!$ DS 236 237 !- 238 ! activate CO2, STOMATE, but not sechiba 239 !- 240 control%river_routing = .FALSE. 241 control%hydrol_cwrr = .FALSE. 242 control%ok_sechiba = .FALSE. 243 ! 244 control%stomate_watchout = .TRUE. 245 control%ok_co2 = .TRUE. 246 control%ok_stomate = .TRUE. 247 !- 248 ! is DGVM activated? 249 !- 250 control%ok_dgvm = .FALSE. 251 CALL getin_p('STOMATE_OK_DGVM',control%ok_dgvm) 252 WRITE(numout,*) 'LPJ is activated: ',control%ok_dgvm 253 254 !!$ >> DS now we search the values for the scalar parameters in some .def files 255 256 ! 07/2011 257 CALL activate_sub_models(control%ok_sechiba,control%river_routing,control%ok_stomate) 258 CALL veget_config 259 ! 188 260 CALL getin_stomate_pft_parameters 189 !!$ DS 190 191 !- 192 ! allocate arrays 193 !- 261 CALL getin_co2_parameters 262 CALL getin_stomate_parameters 263 IF (control%ok_dgvm) THEN 264 CALL getin_dgvm_parameters 265 ENDIF 266 !!$ >> DS 267 268 !- 269 ! restart files 270 !- 271 IF (is_root_prc) THEN 272 ! Sechiba's restart files 273 sec_restname_in = 'sechiba_start.nc' 274 CALL getin('SECHIBA_restart_in',sec_restname_in) 275 WRITE(numout,*) 'SECHIBA INPUT RESTART_FILE: ',TRIM(sec_restname_in) 276 IF ( TRIM(sec_restname_in) .EQ. 'NONE' ) THEN 277 STOP 'Need a restart file for Sechiba' 278 ENDIF 279 sec_restname_out = 'sechiba_rest_out.nc' 280 CALL getin('SECHIBA_rest_out',sec_restname_out) 281 WRITE(numout,*) 'SECHIBA OUTPUT RESTART_FILE: ',TRIM(sec_restname_out) 282 ! Stomate's restart files 283 sto_restname_in = 'stomate_start.nc' 284 CALL getin('STOMATE_RESTART_FILEIN',sto_restname_in) 285 WRITE(numout,*) 'STOMATE INPUT RESTART_FILE: ',TRIM(sto_restname_in) 286 sto_restname_out = 'stomate_rest_out.nc' 287 CALL getin('STOMATE_RESTART_FILEOUT',sto_restname_out) 288 WRITE(numout,*) 'STOMATE OUTPUT RESTART_FILE: ',TRIM(sto_restname_out) 289 290 !- 291 ! We need to know iim, jjm. 292 ! Get them from the restart files themselves. 293 !- 294 iret = NF90_OPEN (sec_restname_in,NF90_NOWRITE,ncfid) 295 IF (iret /= NF90_NOERR) THEN 296 CALL ipslerr (3,'teststomate', & 297 & 'Could not open file : ', & 298 & sec_restname_in,'(Do you have forget it ?)') 299 ENDIF 300 iret = NF90_INQUIRE_DIMENSION (ncfid,1,len=iim_g) 301 iret = NF90_INQUIRE_DIMENSION (ncfid,2,len=jjm_g) 302 iret = NF90_CLOSE (ncfid) 303 304 ENDIF 305 CALL bcast(iim_g) 306 CALL bcast(jjm_g) 307 ! 308 ! Parallelization : 309 ! 310 CALL init_data_para(iim_g,jjm_g,nbp_glo,indices_g) 311 kjpindex=nbp_loc 312 jjm=jj_nb 313 iim=iim_g 314 kjpij=iim*jjm 315 IF (debug) WRITE(numout,*) "Local grid : ",kjpindex,iim,jjm 316 !- 317 !- 318 ! read info about grids 319 !- 320 !- 321 llm=1 322 ALLOCATE(lev(llm)) 323 IF (is_root_prc) THEN 324 !- 325 ier = NF90_INQ_VARID (forcing_id,'lalo',vid) 326 ier = NF90_GET_VAR (forcing_id,vid,lalo_g) 327 !- 328 ALLOCATE (x_neighbours_g(nbp_glo,8),stat=ier) 329 ier = NF90_INQ_VARID (forcing_id,'neighbours',vid) 330 ier = NF90_GET_VAR (forcing_id,vid,x_neighbours_g) 331 neighbours_g(:,:) = NINT(x_neighbours_g(:,:)) 332 DEALLOCATE (x_neighbours_g) 333 !- 334 ier = NF90_INQ_VARID (forcing_id,'resolution',vid) 335 ier = NF90_GET_VAR (forcing_id,vid,resolution_g) 336 !- 337 ier = NF90_INQ_VARID (forcing_id,'contfrac',vid) 338 ier = NF90_GET_VAR (forcing_id,vid,contfrac_g) 339 340 lon_g(:,:) = 0.0 341 lat_g(:,:) = 0.0 342 lev(1) = 0.0 343 !- 344 CALL restini & 345 & (sec_restname_in, iim_g, jjm_g, lon_g, lat_g, llm, lev, & 346 & sec_restname_out, itau_dep, date0, dt, rest_id_sec) 347 !- 348 IF ( dt .NE. dtradia ) THEN 349 WRITE(numout,*) 'dt',dt 350 WRITE(numout,*) 'dtradia',dtradia 351 CALL ipslerr (3,'teststomate', & 352 & 'PROBLEM with time steps.', & 353 & sec_restname_in,'(dt .NE. dtradia)') 354 ENDIF 355 !- 356 CALL restini & 357 & (sto_restname_in, iim_g, jjm_g, lon_g, lat_g, llm, lev, & 358 & sto_restname_out, itau_dep, date0, dt, rest_id_sto) 359 !- 360 IF ( dt .NE. dtradia ) THEN 361 WRITE(numout,*) 'dt',dt 362 WRITE(numout,*) 'dtradia',dtradia 363 CALL ipslerr (3,'teststomate', & 364 & 'PROBLEM with time steps.', & 365 & sto_restname_in,'(dt .NE. dtradia)') 366 ENDIF 367 ENDIF 368 CALL bcast(rest_id_sec) 369 CALL bcast(rest_id_sto) 370 CALL bcast(date0) 371 CALL bcast(dt) 372 CALL bcast(lev) 373 !--- 374 !--- Create the index table 375 !--- 376 !--- This job return a LOCAL kindex 377 !--- 378 ALLOCATE (indices(kjpindex),stat=ier) 379 IF (debug .AND. is_root_prc) WRITE(numout,*) "indices_g = ",indices_g(1:nbp_glo) 380 CALL scatter(indices_g,indices) 381 indices(1:kjpindex)=indices(1:kjpindex)-(jj_begin-1)*iim_g 382 IF (debug) WRITE(numout,*) "indices = ",indices(1:kjpindex) 383 384 !--- 385 !--- initialize global grid 386 !--- 387 CALL init_grid( kjpindex ) 388 CALL grid_stuff (nbp_glo, iim_g, jjm_g, lon_g, lat_g, indices_g) 389 390 !--- 391 !--- initialize local grid 392 !--- 393 jlandindex = (((indices(1:kjpindex)-1)/iim) + 1) 394 if (debug) WRITE(numout,*) "jlandindex = ",jlandindex(1:kjpindex) 395 ilandindex = (indices(1:kjpindex) - (jlandindex(1:kjpindex)-1)*iim) 396 IF (debug) WRITE(numout,*) "ilandindex = ",ilandindex(1:kjpindex) 397 ALLOCATE(lon(iim,jjm)) 398 ALLOCATE(lat(iim,jjm)) 399 CALL scatter2D(lon_g,lon) 400 CALL scatter2D(lat_g,lat) 401 IF (debug) THEN 402 IF (is_root_prc) THEN 403 WRITE(numout,*) mpi_rank, lat_g 404 WRITE(numout,*) mpi_rank, lon_g 405 ENDIF 406 WRITE(numout,*) mpi_rank, lat 407 WRITE(numout,*) mpi_rank, lon 408 ENDIF 409 410 DO ji=1,kjpindex 411 412 j = jlandindex(ji) 413 i = ilandindex(ji) 414 415 !- Create the internal coordinate table 416 !- 417 lalo(ji,1) = lat(i,j) 418 lalo(ji,2) = lon(i,j) 419 ENDDO 420 CALL scatter(neighbours_g,neighbours) 421 CALL scatter(resolution_g,resolution) 422 CALL scatter(contfrac_g,contfrac) 423 CALL scatter(area_g,area) 424 !- 425 !- Check if we have by any change a rectilinear grid. This would allow us to 426 !- simplify the output files. 427 ! 428 rectilinear = .FALSE. 429 IF (is_root_prc) THEN 430 IF ( ALL(lon_g(:,:) == SPREAD(lon_g(:,1), 2, SIZE(lon_g,2))) .AND. & 431 & ALL(lat_g(:,:) == SPREAD(lat_g(1,:), 1, SIZE(lat_g,1))) ) THEN 432 rectilinear = .TRUE. 433 ENDIF 434 ENDIF 435 CALL bcast(rectilinear) 436 IF (rectilinear) THEN 437 ALLOCATE(lon_rect(iim),stat=ier) 438 IF (ier .NE. 0) THEN 439 WRITE (numout,*) ' error in lon_rect allocation. We stop. We need iim words = ',iim 440 STOP 'intersurf_history' 441 ENDIF 442 ALLOCATE(lat_rect(jjm),stat=ier) 443 IF (ier .NE. 0) THEN 444 WRITE (numout,*) ' error in lat_rect allocation. We stop. We need jjm words = ',jjm 445 STOP 'intersurf_history' 446 ENDIF 447 lon_rect(:) = lon(:,1) 448 lat_rect(:) = lat(1,:) 449 ENDIF 450 !- 451 ! allocate arrays 452 !- 453 ! 194 454 a_er = .FALSE. 195 455 !!$ >> DS added for the externalised version … … 200 460 a_er = a_er .OR. (ier.NE.0) 201 461 !!$ end add DS 202 ALLOCATE (indices(kjpindex),stat=ier)203 a_er = a_er .OR. (ier.NE.0)204 462 ALLOCATE (indexveg(kjpindex*nvm), stat=ier) 205 463 a_er = a_er .OR. (ier.NE.0) … … 252 510 ALLOCATE (co2_flux(kjpindex,nvm),stat=ier) 253 511 a_er = a_er .OR. (ier.NE.0) 254 IF ( a_er ) STOP 'PROBLEM WITH ALLOCATION' 512 ALLOCATE (fco2_lu(kjpindex),stat=ier) 513 a_er = a_er .OR. (ier.NE.0) 514 IF (a_er) THEN 515 CALL ipslerr (3,'teststomate', & 516 & 'PROBLEM WITH ALLOCATION', & 517 & 'for local variables 1','') 518 ENDIF 255 519 ! 256 520 ! prepare forcing 257 521 ! 258 522 max_totsize = 50 259 CALL getin ('STOMATE_FORCING_MEMSIZE',max_totsize)523 CALL getin_p ('STOMATE_FORCING_MEMSIZE',max_totsize) 260 524 max_totsize = max_totsize * 1000000 525 261 526 totsize_1step = SIZE(soiltype(:,3))*KIND(soiltype(:,3)) + & 262 SIZE(humrel_x)*KIND(humrel_x) + & 263 SIZE(litterhum)*KIND(litterhum) + & 264 SIZE(t2m)*KIND(t2m) + & 265 SIZE(t2m_min)*KIND(t2m_min) + & 266 SIZE(temp_sol)*KIND(temp_sol) + & 267 SIZE(soiltemp)*KIND(soiltemp) + & 268 SIZE(soilhum)*KIND(soilhum) + & 269 SIZE(precip_rain)*KIND(precip_rain) + & 270 SIZE(precip_snow)*KIND(precip_snow) + & 271 SIZE(gpp_x)*KIND(gpp_x) + & 272 SIZE(veget_force_x)*KIND(veget_force_x) + & 273 SIZE(veget_max_force_x)*KIND(veget_max_force_x) + & 274 SIZE(lai_force_x)*KIND(lai_force_x) 527 SIZE(humrel_x)*KIND(humrel_x) + & 528 SIZE(litterhum)*KIND(litterhum) + & 529 SIZE(t2m)*KIND(t2m) + & 530 SIZE(t2m_min)*KIND(t2m_min) + & 531 SIZE(temp_sol)*KIND(temp_sol) + & 532 SIZE(soiltemp)*KIND(soiltemp) + & 533 SIZE(soilhum)*KIND(soilhum) + & 534 SIZE(precip_rain)*KIND(precip_rain) + & 535 SIZE(precip_snow)*KIND(precip_snow) + & 536 SIZE(gpp_x)*KIND(gpp_x) + & 537 SIZE(veget_force_x)*KIND(veget_force_x) + & 538 SIZE(veget_max_force_x)*KIND(veget_max_force_x) + & 539 SIZE(lai_force_x)*KIND(lai_force_x) 540 CALL reduce_sum(totsize_1step,totsize_tmp) 541 CALL bcast(totsize_tmp) 542 totsize_1step=totsize_tmp 543 275 544 ! total number of forcing steps 276 nsft = INT(one_year/(dt_force/one_day)) 545 IF ( nsft .NE. INT(one_year/(dt_force/one_day)) ) THEN 546 CALL ipslerr (3,'teststomate', & 547 & 'stomate: error with total number of forcing steps', & 548 & 'nsft','teststomate computation different with forcing file value.') 549 ENDIF 277 550 ! number of forcing steps in memory 278 nsfm = MIN(nsft,MAX(1,NINT(FLOAT(max_totsize)/FLOAT(totsize_1step)))) 551 nsfm = MIN(nsft, & 552 & MAX(1,NINT( REAL(max_totsize,r_std) & 553 & /REAL(totsize_1step,r_std)))) 279 554 !- 280 555 WRITE(numout,*) 'Offline forcing of Stomate:' … … 282 557 WRITE(numout,*) ' Number of forcing steps in memory:',nsfm 283 558 !- 284 a_er = .FALSE. 285 ALLOCATE (clay_fm(kjpindex,nsfm),stat=ier) 286 a_er = a_er.OR.(ier.NE.0) 287 ALLOCATE (humrel_x_fm(kjpindex,nvm,nsfm),stat=ier) 288 a_er = a_er.OR.(ier.NE.0) 289 ALLOCATE (litterhum_fm(kjpindex,nsfm),stat=ier) 290 a_er = a_er.OR.(ier.NE.0) 291 ALLOCATE (t2m_fm(kjpindex,nsfm),stat=ier) 292 a_er = a_er.OR.(ier.NE.0) 293 ALLOCATE (t2m_min_fm(kjpindex,nsfm),stat=ier) 294 a_er = a_er.OR.(ier.NE.0) 295 ALLOCATE (temp_sol_fm(kjpindex,nsfm),stat=ier) 296 a_er = a_er.OR.(ier.NE.0) 297 ALLOCATE (soiltemp_fm(kjpindex,nbdl,nsfm),stat=ier) 298 a_er = a_er.OR.(ier.NE.0) 299 ALLOCATE (soilhum_fm(kjpindex,nbdl,nsfm),stat=ier) 300 a_er = a_er.OR.(ier.NE.0) 301 ALLOCATE (precip_fm(kjpindex,nsfm),stat=ier) 302 a_er = a_er.OR.(ier.NE.0) 303 ALLOCATE (gpp_x_fm(kjpindex,nvm,nsfm),stat=ier) 304 a_er = a_er.OR.(ier.NE.0) 305 ALLOCATE (veget_force_x_fm(kjpindex,nvm,nsfm),stat=ier) 306 a_er = a_er.OR.(ier.NE.0) 307 ALLOCATE (veget_max_force_x_fm(kjpindex,nvm,nsfm),stat=ier) 308 a_er = a_er.OR.(ier.NE.0) 309 ALLOCATE (lai_force_x_fm(kjpindex,nvm,nsfm),stat=ier) 310 a_er = a_er.OR.(ier.NE.0) 311 ALLOCATE (isf(nsfm),stat=ier) 312 a_er = a_er.OR.(ier.NE.0) 313 ALLOCATE (nf_written(nsft),stat=ier) 314 a_er = a_er.OR.(ier.NE.0) 315 ALLOCATE (nf_cumul(nsft),stat=ier) 316 a_er = a_er.OR.(ier.NE.0) 317 IF (a_er) THEN 318 STOP 'stomate: error in memory allocation for forcing data' 319 ENDIF 559 CALL init_forcing(kjpindex,nsfm,nsft) 560 !- 320 561 ! ensure that we read all new forcing states 321 562 iisf = nsfm … … 323 564 ! of the forcing states that will be in memory 324 565 isf(:) = (/ (i,i=1,nsfm) /) 325 !- 326 ! read info about grids 327 !- 328 contfrac(:) = 1.0 329 !- 330 ALLOCATE (x_indices(kjpindex),stat=ier) 331 ier = NF90_INQ_VARID (force_id,'index',v_id) 332 ier = NF90_GET_VAR (force_id,v_id,x_indices) 333 indices(:) = NINT(x_indices(:)) 334 DEALLOCATE (x_indices) 335 !- 336 DO ji=1,kjpindex 337 DO jv=1,nvm 338 indexveg((jv-1)*kjpindex+ji) = indices(ji)+(jv-1)*kjpij 339 ENDDO 340 ENDDO 341 !- 342 ier = NF90_INQ_VARID (force_id,'lalo',v_id) 343 ier = NF90_GET_VAR (force_id,v_id,lalo) 344 !- 345 ALLOCATE (x_neighbours(kjpindex,8),stat=ier) 346 ier = NF90_INQ_VARID (force_id,'neighbours',v_id) 347 ier = NF90_GET_VAR (force_id,v_id,x_neighbours) 348 neighbours(:,:) = NINT(x_neighbours(:,:)) 349 DEALLOCATE (x_neighbours) 350 !- 351 ier = NF90_INQ_VARID (force_id,'resolution',v_id) 352 ier = NF90_GET_VAR (force_id,v_id,resolution) 353 !- 354 ier = NF90_INQ_VARID (force_id,'contfrac',v_id) 355 ier = NF90_GET_VAR (force_id,v_id,contfrac) 356 !- 357 ! activate CO2, STOMATE, but not sechiba 358 !- 359 control%river_routing = .FALSE. 360 control%hydrol_cwrr = .FALSE. 361 control%ok_sechiba = .FALSE. 362 ! 363 control%stomate_watchout = .TRUE. 364 control%ok_co2 = .TRUE. 365 control%ok_stomate = .TRUE. 366 !!$ >> DS now we search the values for the scalar parameters in some .def files 367 ! 07/2011 368 CALL activate_sub_models(control%ok_sechiba,control%river_routing,control%ok_stomate) 369 CALL veget_config 370 ! 371 CALL getin_co2_parameters 372 CALL getin_stomate_parameters 373 !!$ >> DS 374 !- 375 ! is DGVM activated? 376 !- 377 control%ok_dgvm = .FALSE. 378 CALL getin('STOMATE_OK_DGVM',control%ok_dgvm) 379 WRITE(*,*) 'LPJ is activated: ',control%ok_dgvm 380 !!$ >> DS for the externalisation reading the scalar parameters when ok_dgvm is set to true 381 IF (control%ok_dgvm) THEN 382 CALL getin_dgvm_parameters 383 ENDIF 384 !!$ >> DS 385 !- 386 ! restart files 387 !- 388 ! Sechiba's restart files 389 sec_restname_in = 'sechiba_start.nc' 390 CALL getin('SECHIBA_restart_in',sec_restname_in) 391 WRITE(*,*) 'SECHIBA INPUT RESTART_FILE: ',TRIM(sec_restname_in) 392 IF ( TRIM(sec_restname_in) .EQ. 'NONE' ) THEN 393 STOP 'Need a restart file for Sechiba' 394 ENDIF 395 sec_restname_out = 'sechiba_restart.nc' 396 CALL getin('SECHIBA_rest_out',sec_restname_out) 397 WRITE(*,*) 'SECHIBA OUTPUT RESTART_FILE: ',TRIM(sec_restname_out) 398 ! Stomate's restart files 399 sto_restname_in = 'stomate_start.nc' 400 CALL getin ('STOMATE_RESTART_FILEIN',sto_restname_in) 401 WRITE(*,*) 'STOMATE INPUT RESTART_FILE: ',TRIM(sto_restname_in) 402 sto_restname_out = 'stomate_restart.nc' 403 CALL getin ('STOMATE_RESTART_FILEOUT',sto_restname_out) 404 WRITE(*,*) 'STOMATE OUTPUT RESTART_FILE: ',TRIM(sto_restname_out) 405 !- 406 ! We need to know iim, jjm. 407 ! Get them from the restart files themselves. 408 !- 409 iret = NF90_OPEN (sec_restname_in,NF90_NOWRITE,ncfid) 410 IF (iret /= NF90_NOERR) THEN 411 CALL ipslerr (3,'teststomate', & 412 & 'Could not open file : ', & 413 & sec_restname_in,'(Do you have forget it ?)') 414 ENDIF 415 iret = NF90_INQUIRE_DIMENSION (ncfid,1,len=iim) 416 iret = NF90_INQUIRE_DIMENSION (ncfid,2,len=jjm) 417 iret = NF90_CLOSE (ncfid) 418 ! Allocate longitudes and latitudes 419 ALLOCATE (lon(iim,jjm),stat=ier) 420 a_er = a_er.OR.(ier.NE.0) 421 ALLOCATE (lat(iim,jjm),stat=ier) 422 a_er = a_er.OR.(ier.NE.0) 423 lon(:,:) = 0.0 424 lat(:,:) = 0.0 425 lev(1) = 0.0 426 !- 427 CALL restini & 428 & (sec_restname_in, iim, jjm, lon, lat, llm, lev, & 429 & sec_restname_out, itau_dep, date0, dt_files, rest_id_sec) 430 !- 431 CALL restini & 432 & (sto_restname_in, iim, jjm, lon, lat, llm, lev, & 433 & sto_restname_out, itau_dep, date0, dt_files, rest_id_sto) 434 !- 435 IF ( dt_files .NE. dtradia ) THEN 436 WRITE(*,*) 'dt_files',dt_files 437 WRITE(*,*) 'dtradia',dtradia 438 STOP 'PROBLEM with time steps.' 439 ENDIF 566 567 nf_written(:) = .FALSE. 568 nf_written(isf(:)) = .TRUE. 569 440 570 !- 441 571 ! a time step for STOMATE corresponds to itau_step timesteps in SECHIBA 442 572 !- 443 573 itau_step = NINT(dt_force/dtradia) 574 IF (debug) WRITE(numout,*) "dtradia, dt_rest, dt_force, itau_step",dtradia, dt, dt_force, itau_step 444 575 ! 445 576 CALL ioconf_startdate(date0) … … 451 582 CALL getin ('TIME_LENGTH', time_str) 452 583 ! transform into itau 453 CALL tlen2itau(time_str, dt _files, date0, itau_len)584 CALL tlen2itau(time_str, dt, date0, itau_len) 454 585 ! itau_len*dtradia must be a multiple of dt_force 455 586 itau_len = NINT( MAX(1.,FLOAT(NINT(itau_len*dtradia/dt_force))) & 456 & *dt_force/dtradia)457 !-458 ! set up STOMATE history file459 460 461 462 463 464 465 466 467 468 469 !-587 & *dt_force/dtradia) 588 !- 589 ! set up STOMATE history file 590 !- 591 !Config Key = STOMATE_OUTPUT_FILE 592 !Config Desc = Name of file in which STOMATE's output is going 593 !Config to be written 594 !Config Def = stomate_history.nc 595 !Config Help = This file is going to be created by the model 596 !Config and will contain the output from the model. 597 !Config This file is a truly COADS compliant netCDF file. 598 !Config It will be generated by the hist software from 599 !Config the IOIPSL package. 600 !- 470 601 stom_histname='stomate_history.nc' 471 CALL getin ('STOMATE_OUTPUT_FILE', stom_histname)472 WRITE( *,*) 'STOMATE_OUTPUT_FILE', TRIM(stom_histname)473 474 475 476 477 478 !-602 CALL getin_p ('STOMATE_OUTPUT_FILE', stom_histname) 603 WRITE(numout,*) 'STOMATE_OUTPUT_FILE', TRIM(stom_histname) 604 !- 605 !Config Key = STOMATE_HIST_DT 606 !Config Desc = STOMATE history time step (d) 607 !Config Def = 10. 608 !Config Help = Time step of the STOMATE history file 609 !- 479 610 hist_days_stom = 10. 480 CALL getin ('STOMATE_HIST_DT', hist_days_stom)611 CALL getin_p ('STOMATE_HIST_DT', hist_days_stom) 481 612 IF ( hist_days_stom == -1. ) THEN 482 613 hist_dt_stom = -1. … … 488 619 ENDIF 489 620 !- 490 ! initialize 491 CALL histbeg(stom_histname, iim, lon, jjm, lat, 1, iim, 1, jjm, & 492 & itau_dep, date0, dt_files, hori_id, hist_id_sto) 493 ! define PFT axis 621 !- 622 !- initialize 623 WRITE(numout,*) "before histbeg : ",date0,dt 624 IF ( rectilinear ) THEN 625 #ifdef CPP_PARA 626 CALL histbeg(stom_histname, iim, lon_rect, jjm, lat_rect, 1, iim, 1, jjm, & 627 & itau_dep, date0, dt, hori_id, hist_id_stom, domain_id=orch_domain_id) 628 #else 629 CALL histbeg(stom_histname, iim, lon_rect, jjm, lat_rect, 1, iim, 1, jjm, & 630 & itau_dep, date0, dt, hori_id, hist_id_stom) 631 #endif 632 ELSE 633 #ifdef CPP_PARA 634 CALL histbeg(stom_histname, iim, lon, jjm, lat, 1, iim, 1, jjm, & 635 & itau_dep, date0, dt, hori_id, hist_id_stom, domain_id=orch_domain_id) 636 #else 637 CALL histbeg(stom_histname, iim, lon, jjm, lat, 1, iim, 1, jjm, & 638 & itau_dep, date0, dt, hori_id, hist_id_stom) 639 #endif 640 ENDIF 641 !- define PFT axis 494 642 hist_PFTaxis = (/ ( REAL(i,r_std), i=1,nvm ) /) 495 !declare this axis496 CALL histvert (hist_id_sto , 'PFT', 'Plant functional type', &497 & '-', nvm, hist_PFTaxis, hist_PFTaxis_id)643 !- declare this axis 644 CALL histvert (hist_id_stom, 'PFT', 'Plant functional type', & 645 & '1', nvm, hist_PFTaxis, hist_PFTaxis_id) 498 646 !- define Pool_10 axis 499 647 hist_pool_10axis = (/ ( REAL(i,r_std), i=1,10 ) /) 500 648 !- declare this axis 501 CALL histvert (hist_id_sto, 'P10', 'Pool 10 years', & 502 & '-', 10, hist_pool_10axis, hist_pool_10axis_id) 649 CALL histvert (hist_id_stom, 'P10', 'Pool 10 years', & 650 & '1', 10, hist_pool_10axis, hist_pool_10axis_id) 651 503 652 !- define Pool_100 axis 504 653 hist_pool_100axis = (/ ( REAL(i,r_std), i=1,100 ) /) 505 654 !- declare this axis 506 CALL histvert (hist_id_sto, 'P100', 'Pool 100 years', & 507 & '-', 100, hist_pool_100axis, hist_pool_100axis_id) 655 CALL histvert (hist_id_stom, 'P100', 'Pool 100 years', & 656 & '1', 100, hist_pool_100axis, hist_pool_100axis_id) 657 508 658 !- define Pool_11 axis 509 659 hist_pool_11axis = (/ ( REAL(i,r_std), i=1,11 ) /) 510 660 !- declare this axis 511 CALL histvert (hist_id_sto, 'P11', 'Pool 10 years + 1', & 512 & '-', 11, hist_pool_11axis, hist_pool_11axis_id) 661 CALL histvert (hist_id_stom, 'P11', 'Pool 10 years + 1', & 662 & '1', 11, hist_pool_11axis, hist_pool_11axis_id) 663 513 664 !- define Pool_101 axis 514 665 hist_pool_101axis = (/ ( REAL(i,r_std), i=1,101 ) /) 515 666 !- declare this axis 516 CALL histvert (hist_id_sto, 'P101', 'Pool 100 years + 1', & 517 & '-', 101, hist_pool_101axis, hist_pool_101axis_id) 518 ! define STOMATE history file 519 CALL stom_define_history (hist_id_sto, nvm, iim, jjm, & 520 & dt_files, hist_dt_stom, hori_id, hist_PFTaxis_id, & 667 CALL histvert (hist_id_stom, 'P101', 'Pool 100 years + 1', & 668 & '1', 101, hist_pool_101axis, hist_pool_101axis_id) 669 670 !- define STOMATE history file 671 CALL stom_define_history (hist_id_stom, nvm, iim, jjm, & 672 & dt, hist_dt_stom, hori_id, hist_PFTaxis_id, & 521 673 & hist_pool_10axis_id, hist_pool_100axis_id, & 522 674 & hist_pool_11axis_id, hist_pool_101axis_id) 523 ! end definition 524 CALL histend(hist_id_sto) 675 676 !- end definition 677 CALL histend(hist_id_stom) 678 !- 679 !- 680 ! STOMATE IPCC OUTPUTS IS ACTIVATED 681 !- 682 !Config Key = STOMATE_IPCC_OUTPUT_FILE 683 !Config Desc = Name of file in which STOMATE's output is going 684 !Config to be written 685 !Config Def = stomate_ipcc_history.nc 686 !Config Help = This file is going to be created by the model 687 !Config and will contain the output from the model. 688 !Config This file is a truly COADS compliant netCDF file. 689 !Config It will be generated by the hist software from 690 !Config the IOIPSL package. 691 !- 692 stom_ipcc_histname='stomate_ipcc_history.nc' 693 CALL getin_p('STOMATE_IPCC_OUTPUT_FILE', stom_ipcc_histname) 694 WRITE(numout,*) 'STOMATE_IPCC_OUTPUT_FILE', TRIM(stom_ipcc_histname) 695 !- 696 !Config Key = STOMATE_IPCC_HIST_DT 697 !Config Desc = STOMATE IPCC history time step (d) 698 !Config Def = 0. 699 !Config Help = Time step of the STOMATE IPCC history file 700 !- 701 hist_days_stom_ipcc = zero 702 CALL getin_p('STOMATE_IPCC_HIST_DT', hist_days_stom_ipcc) 703 IF ( hist_days_stom_ipcc == moins_un ) THEN 704 hist_dt_stom_ipcc = moins_un 705 WRITE(numout,*) 'output frequency for STOMATE IPCC history file (d): one month.' 706 ELSE 707 hist_dt_stom_ipcc = NINT( hist_days_stom_ipcc ) * one_day 708 WRITE(numout,*) 'output frequency for STOMATE IPCC history file (d): ', & 709 hist_dt_stom_ipcc/one_day 710 ENDIF 711 712 ! test consistency between STOMATE_IPCC_HIST_DT and DT_SLOW parameters 713 dt_slow_ = one_day 714 CALL getin_p('DT_SLOW', dt_slow_) 715 IF ( hist_days_stom_ipcc > zero ) THEN 716 IF (dt_slow_ > hist_dt_stom_ipcc) THEN 717 WRITE(numout,*) "DT_SLOW = ",dt_slow_," , STOMATE_IPCC_HIST_DT = ",hist_dt_stom_ipcc 718 CALL ipslerr (3,'intsurf_history', & 719 & 'Problem with DT_SLOW > STOMATE_IPCC_HIST_DT','', & 720 & '(must be less or equal)') 721 ENDIF 722 ENDIF 723 724 IF ( hist_dt_stom_ipcc == 0 ) THEN 725 hist_id_stom_ipcc = -1 726 ELSE 727 !- 728 !- initialize 729 IF ( rectilinear ) THEN 730 #ifdef CPP_PARA 731 CALL histbeg(stom_ipcc_histname, iim, lon_rect, jjm, lat_rect, 1, iim, 1, jjm, & 732 & itau_dep, date0, dt, hori_id, hist_id_stom_ipcc,domain_id=orch_domain_id) 733 #else 734 CALL histbeg(stom_ipcc_histname, iim, lon_rect, jjm, lat_rect, 1, iim, 1, jjm, & 735 & itau_dep, date0, dt, hori_id, hist_id_stom_ipcc) 736 #endif 737 ELSE 738 #ifdef CPP_PARA 739 CALL histbeg(stom_ipcc_histname, iim, lon, jjm, lat, 1, iim, 1, jjm, & 740 & itau_dep, date0, dt, hori_id, hist_id_stom_ipcc,domain_id=orch_domain_id) 741 #else 742 CALL histbeg(stom_ipcc_histname, iim, lon, jjm, lat, 1, iim, 1, jjm, & 743 & itau_dep, date0, dt, hori_id, hist_id_stom_ipcc) 744 #endif 745 ENDIF 746 !- declare this axis 747 CALL histvert (hist_id_stom_IPCC, 'PFT', 'Plant functional type', & 748 & '1', nvm, hist_PFTaxis, hist_IPCC_PFTaxis_id) 749 750 !- define STOMATE history file 751 CALL stom_IPCC_define_history (hist_id_stom_IPCC, nvm, iim, jjm, & 752 & dt, hist_dt_stom_ipcc, hori_id, hist_IPCC_PFTaxis_id) 753 754 !- end definition 755 CALL histend(hist_id_stom_IPCC) 756 757 ENDIF 758 ! 759 CALL histwrite(hist_id_stom, 'Areas', itau_dep+itau_step, area, kjpindex, indices) 760 IF ( hist_id_stom_IPCC > 0 ) THEN 761 CALL histwrite(hist_id_stom_IPCC, 'Areas', itau_dep+itau_step, area, kjpindex, indices) 762 ENDIF 763 ! 525 764 hist_id_sec = -1 526 765 hist_id_sec2 = -1 527 hist_id_stom_IPCC = -1 528 !- 529 ! read some variables we need from SECHIBA's restart file 530 !- 766 !- 767 ! first call of slowproc to initialize variables 768 !- 769 itau = itau_dep 770 ! 771 DO ji=1,kjpindex 772 DO jv=1,nvm 773 indexveg((jv-1)*kjpindex + ji) = indices(ji) + (jv-1)*kjpij 774 ENDDO 775 ENDDO 776 !- 777 !MM Problem here with dpu which depends on soil type 778 DO l = 1, nbdl-1 779 ! first 2.0 is dpu 780 ! second 2.0 is average 781 diaglev(l) = dpu_cste/(2**(nbdl-1) -1) * ( ( 2**(l-1) -1) + ( 2**(l) -1) ) / 2.0 782 ENDDO 783 diaglev(nbdl) = dpu_cste 784 ! 785 786 ! >> DS : getin the value of slowproc parameters 787 ! 788 !Config Key = CLAYFRACTION_DEFAULT 789 !Config Desc = 790 !Config If = OK_SECHIBA 791 !Config Def = 0.2 792 !Config Help = 793 !Config Units = NONE 794 CALL getin_p('CLAYFRACTION_DEFAULT',clayfraction_default) 795 ! 796 !Config Key = MIN_VEGFRAC 797 !Config Desc = Minimal fraction of mesh a vegetation type can occupy 798 !Config If = OK_SECHIBA 799 !Config Def = 0.001 800 !Config Help = 801 !Config Units = NONE 802 CALL getin_p('MIN_VEGFRAC',min_vegfrac) 803 ! 804 !Config Key = STEMPDIAG_BID 805 !Config Desc = only needed for an initial LAI if there is no restart file 806 !Config If = OK_SECHIBA 807 !Config Def = 280. 808 !Config Help = 809 !Config Units = 810 CALL getin_p('STEMPDIAG_BID',stempdiag_bid) 811 ! 812 !Config Key = SOILTYPE_DEFAULT 813 !Config Desc = Default soil texture distribution in the following order : sand, loam and clay 814 !Config If = OK_SECHIBA 815 !Config Def = 0.0, 1.0, 0.0 816 !Config Help = 817 !Config Units = NONE 818 CALL getin('SOILTYPE_DEFAULT',soiltype_default) 819 ! >> DS 820 531 821 CALL ioget_expval(val_exp) 532 !-533 ! first call of slowproc to initialize variables534 !-535 itau = itau_dep536 822 ldrestart_read = .TRUE. 537 823 ldrestart_write = .FALSE. 538 !- 539 !MM Problem here with dpu which depends on soil type 540 DO jv = 1, nbdl-1 541 ! first 2.0 is dpu 542 ! second 2.0 is average 543 diaglev(jv) = 2.0/(2**(nbdl-1) -1) * ( ( 2**(jv-1) -1) + ( 2**(jv) -1) ) / 2.0 544 ENDDO 545 diaglev(nbdl) = 2.0 546 !- 547 ! For sequential use only, we must initialize data_para : 548 ! 549 CALL init_para(.FALSE.) 550 ! 551 CALL init_data_para(iim,jjm,kjpindex,indices) 552 ! 553 !- global index index_g is the index_l of root proc 554 IF (is_root_prc) index_g(:)=indices(1:kjpindex) 555 !- 556 557 ! >> DS : getin the value of slowproc parameters 558 CALL getin('CLAYFRACTION_DEFAULT',clayfraction_default) 559 CALL getin('SOILTYPE_DEFAULT',soiltype_default) 560 ! >> DS 561 824 !- 825 ! read some variables we need from SECHIBA's restart file 826 !- 562 827 CALL slowproc_main & 563 828 & (itau, kjpij, kjpindex, dt_force, date0, & … … 568 833 & deadleaf_cover, assim_param_x, lai_x, height_x, veget_x, & 569 834 & frac_nobio, veget_max_x, totfrac_nobio, qsintmax_x, & 570 & rest_id_sec, hist_id_sec, hist_id_sec2, rest_id_sto, hist_id_sto, hist_id_stom_IPCC, co2_flux, & 571 !!>> DS add 01/07/2011 correction due to some modifications in the trunk 572 & fco2_lu) 835 & rest_id_sec, hist_id_sec, hist_id_sec2, rest_id_sto, hist_id_stom, hist_id_stom_IPCC, co2_flux, fco2_lu) 573 836 ! correct date 574 837 day_counter = one_day - dt_force … … 577 840 ! time loop 578 841 !- 842 IF (debug) check_time=.TRUE. 843 CALL intsurf_time( itau_dep+itau_step, date0, dtradia ) 844 l_first_intersurf=.FALSE. 845 ! 579 846 DO itau = itau_dep+itau_step,itau_dep+itau_len,itau_step 847 ! 580 848 !-- next forcing state 581 849 iisf = iisf+1 … … 587 855 !---- determine blocks of forcing states that are contiguous in memory 588 856 !----- 589 nblocks = 0 590 ifirst(:) = 1 591 ilast(:) = 1 592 !----- 593 DO iisf=1,nsfm 594 IF ( (nblocks .NE. 0) ) THEN 595 IF ( (isf(iisf) .EQ. isf(ilast(nblocks))+1) ) THEN 596 !-------- element is contiguous with last element found 597 ilast(nblocks) = iisf 598 ELSE 599 !-------- found first element of new block 600 nblocks = nblocks+1 601 IF (nblocks .GT. nsfm) THEN 602 ! IF (nblocks .GT. 2) THEN 603 STOP 'Problem in teststomate' 604 ENDIF 605 ifirst(nblocks) = iisf 606 ilast(nblocks) = iisf 607 ENDIF 608 ELSE 609 !-------- found first element of new block 610 nblocks = nblocks+1 611 IF (nblocks .GT. nsfm) THEN 612 ! IF (nblocks .GT. 2) THEN 613 STOP 'Problem in teststomate' 614 ENDIF 615 ifirst(nblocks) = iisf 616 ilast(nblocks) = iisf 617 ENDIF 618 ENDDO 619 !----- 620 DO iblocks=1,nblocks 621 IF (ifirst(iblocks) .NE. ilast(iblocks)) THEN 622 ndim = 2; 623 start(:) = 1; start(ndim) = isf(ifirst(iblocks)); 624 count(1:ndim) = SHAPE(clay_fm) 625 count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1 626 ier = NF90_INQ_VARID (force_id,'clay',v_id) 627 a_er = a_er.OR.(ier.NE.0) 628 ier = NF90_GET_VAR (force_id,v_id, & 629 & clay_fm(:,ifirst(iblocks):ilast(iblocks)), & 630 & start=start(1:ndim), count=count(1:ndim)) 631 a_er = a_er.OR.(ier.NE.0) 632 !--------- 633 ndim = 3; 634 start(:) = 1; start(ndim) = isf(ifirst(iblocks)); 635 count(1:ndim) = SHAPE(humrel_x_fm) 636 count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1 637 ier = NF90_INQ_VARID (force_id,'humrel',v_id) 638 a_er = a_er.OR.(ier.NE.0) 639 ier = NF90_GET_VAR (force_id,v_id, & 640 & humrel_x_fm(:,:,ifirst(iblocks):ilast(iblocks)), & 641 & start=start(1:ndim), count=count(1:ndim)) 642 a_er = a_er.OR.(ier.NE.0) 643 !--------- 644 ndim = 2; 645 start(:) = 1; start(ndim) = isf(ifirst(iblocks)); 646 count(1:ndim) = SHAPE(litterhum_fm) 647 count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1 648 ier = NF90_INQ_VARID (force_id,'litterhum',v_id) 649 a_er = a_er.OR.(ier.NE.0) 650 ier = NF90_GET_VAR (force_id,v_id, & 651 & litterhum_fm(:,ifirst(iblocks):ilast(iblocks)), & 652 & start=start(1:ndim), count=count(1:ndim)) 653 a_er = a_er.OR.(ier.NE.0) 654 !--------- 655 ndim = 2; 656 start(:) = 1; start(ndim) = isf(ifirst(iblocks)); 657 count(1:ndim) = SHAPE(t2m_fm) 658 count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1 659 ier = NF90_INQ_VARID (force_id,'t2m',v_id) 660 a_er = a_er.OR.(ier.NE.0) 661 ier = NF90_GET_VAR (force_id,v_id, & 662 & t2m_fm(:,ifirst(iblocks):ilast(iblocks)), & 663 & start=start(1:ndim), count=count(1:ndim)) 664 a_er = a_er.OR.(ier.NE.0) 665 !--------- 666 ndim = 2; 667 start(:) = 1; start(ndim) = isf(ifirst(iblocks)); 668 count(1:ndim) = SHAPE(t2m_min_fm) 669 count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1 670 ier = NF90_INQ_VARID (force_id,'t2m_min',v_id) 671 a_er = a_er.OR.(ier.NE.0) 672 ier = NF90_GET_VAR (force_id,v_id, & 673 & t2m_min_fm(:,ifirst(iblocks):ilast(iblocks)), & 674 & start=start(1:ndim), count=count(1:ndim)) 675 a_er = a_er.OR.(ier.NE.0) 676 !--------- 677 ndim = 2; 678 start(:) = 1; start(ndim) = isf(ifirst(iblocks)); 679 count(1:ndim) = SHAPE(temp_sol_fm) 680 count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1 681 ier = NF90_INQ_VARID (force_id,'tsurf',v_id) 682 a_er = a_er.OR.(ier.NE.0) 683 ier = NF90_GET_VAR (force_id,v_id, & 684 & temp_sol_fm(:,ifirst(iblocks):ilast(iblocks)), & 685 & start=start(1:ndim), count=count(1:ndim)) 686 a_er = a_er.OR.(ier.NE.0) 687 !--------- 688 ndim = 3; 689 start(:) = 1; start(ndim) = isf(ifirst(iblocks)); 690 count(1:ndim) = SHAPE(soiltemp_fm) 691 count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1 692 ier = NF90_INQ_VARID (force_id,'tsoil',v_id) 693 a_er = a_er.OR.(ier.NE.0) 694 ier = NF90_GET_VAR (force_id,v_id, & 695 & soiltemp_fm(:,:,ifirst(iblocks):ilast(iblocks)), & 696 & start=start(1:ndim), count=count(1:ndim)) 697 a_er = a_er.OR.(ier.NE.0) 698 !--------- 699 ndim = 3; 700 start(:) = 1; start(ndim) = isf(ifirst(iblocks)); 701 count(1:ndim) = SHAPE(soilhum_fm) 702 count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1 703 ier = NF90_INQ_VARID (force_id,'soilhum',v_id) 704 a_er = a_er.OR.(ier.NE.0) 705 ier = NF90_GET_VAR (force_id,v_id, & 706 & soilhum_fm(:,:,ifirst(iblocks):ilast(iblocks)), & 707 & start=start(1:ndim), count=count(1:ndim)) 708 a_er = a_er.OR.(ier.NE.0) 709 !--------- 710 ndim = 2; 711 start(:) = 1; start(ndim) = isf(ifirst(iblocks)); 712 count(1:ndim) = SHAPE(precip_fm) 713 count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1 714 ier = NF90_INQ_VARID (force_id,'precip',v_id) 715 a_er = a_er.OR.(ier.NE.0) 716 ier = NF90_GET_VAR (force_id,v_id, & 717 & precip_fm(:,ifirst(iblocks):ilast(iblocks)), & 718 & start=start(1:ndim), count=count(1:ndim)) 719 a_er = a_er.OR.(ier.NE.0) 720 !--------- 721 ndim = 3; 722 start(:) = 1; start(ndim) = isf(ifirst(iblocks)); 723 count(1:ndim) = SHAPE(gpp_x_fm) 724 count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1 725 ier = NF90_INQ_VARID (force_id,'gpp',v_id) 726 a_er = a_er.OR.(ier.NE.0) 727 ier = NF90_GET_VAR (force_id,v_id, & 728 & gpp_x_fm(:,:,ifirst(iblocks):ilast(iblocks)), & 729 & start=start(1:ndim), count=count(1:ndim)) 730 a_er = a_er.OR.(ier.NE.0) 731 !--------- 732 ndim = 3; 733 start(:) = 1; start(ndim) = isf(ifirst(iblocks)); 734 count(1:ndim) = SHAPE(veget_force_x_fm) 735 count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1 736 ier = NF90_INQ_VARID (force_id,'veget',v_id) 737 a_er = a_er.OR.(ier.NE.0) 738 ier = NF90_GET_VAR (force_id,v_id, & 739 & veget_force_x_fm(:,:,ifirst(iblocks):ilast(iblocks)), & 740 & start=start(1:ndim), count=count(1:ndim)) 741 a_er = a_er.OR.(ier.NE.0) 742 !--------- 743 ndim = 3; 744 start(:) = 1; start(ndim) = isf(ifirst(iblocks)); 745 count(1:ndim) = SHAPE(veget_max_force_x_fm) 746 count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1 747 ier = NF90_INQ_VARID (force_id,'veget_max',v_id) 748 a_er = a_er.OR.(ier.NE.0) 749 ier = NF90_GET_VAR (force_id,v_id, & 750 & veget_max_force_x_fm(:,:,ifirst(iblocks):ilast(iblocks)), & 751 & start=start(1:ndim), count=count(1:ndim)) 752 a_er = a_er.OR.(ier.NE.0) 753 !--------- 754 ndim = 3; 755 start(:) = 1; start(ndim) = isf(ifirst(iblocks)); 756 count(1:ndim) = SHAPE(lai_force_x_fm) 757 count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1 758 ier = NF90_INQ_VARID (force_id,'lai',v_id) 759 a_er = a_er.OR.(ier.NE.0) 760 ier = NF90_GET_VAR (force_id,v_id, & 761 & lai_force_x_fm(:,:,ifirst(iblocks):ilast(iblocks)), & 762 & start=start(1:ndim), count=count(1:ndim)) 763 a_er = a_er.OR.(ier.NE.0) 764 ENDIF 765 ENDDO 857 CALL forcing_read(forcing_id,nsfm) 858 859 !-------------------------- 860 766 861 !----- 767 862 !---- determine which forcing states must be read next time … … 769 864 isf(1) = isf(nsfm)+1 770 865 IF ( isf(1) .GT. nsft ) isf(1) = 1 771 DO iisf = 2, nsfm 772 isf(iisf) = isf(iisf-1)+1 773 IF ( isf(iisf) .GT. nsft ) isf(iisf) = 1 774 ENDDO 866 DO iiisf = 2, nsfm 867 isf(iiisf) = isf(iiisf-1)+1 868 IF ( isf(iiisf) .GT. nsft ) isf(iiisf) = 1 869 ENDDO 870 nf_written(isf(:)) = .TRUE. 775 871 !---- start again at first forcing state 776 iisf = 1 777 ENDIF 778 soiltype(:,3) = clay_fm(:,iisf) 779 humrel_x(:,:) = humrel_x_fm(:,:,iisf) 780 litterhum(:) = litterhum_fm(:,iisf) 781 t2m(:) = t2m_fm(:,iisf) 782 t2m_min(:) = t2m_min_fm(:,iisf) 783 temp_sol(:) = temp_sol_fm(:,iisf) 784 soiltemp(:,:) = soiltemp_fm(:,:,iisf) 785 soilhum(:,:) = soilhum_fm(:,:,iisf) 786 precip_rain(:) = precip_fm(:,iisf) 787 gpp_x(:,:) = gpp_x_fm(:,:,iisf) 788 veget_force_x(:,:) = veget_force_x_fm(:,:,iisf) 789 veget_max_force_x(:,:) = veget_max_force_x_fm(:,:,iisf) 790 lai_force_x(:,:) = lai_force_x_fm(:,:,iisf) 791 WHERE ( t2m(:) .LT. ZeroCelsius ) 792 precip_snow(:) = precip_rain(:) 793 precip_rain(:) = 0. 794 ELSEWHERE 795 precip_snow(:) = 0. 796 ENDWHERE 872 iisf = 1 873 ENDIF 874 ! Bug here ! soiltype(:,3) != clay 875 ! soiltype(:,3) = clay_fm(:,iisf) 876 humrel_x(:,:) = humrel_daily_fm(:,:,iisf) 877 litterhum(:) = litterhum_daily_fm(:,iisf) 878 t2m(:) = t2m_daily_fm(:,iisf) 879 t2m_min(:) = t2m_min_daily_fm(:,iisf) 880 temp_sol(:) = tsurf_daily_fm(:,iisf) 881 soiltemp(:,:) = tsoil_daily_fm(:,:,iisf) 882 soilhum(:,:) = soilhum_daily_fm(:,:,iisf) 883 precip_rain(:) = precip_fm(:,iisf) 884 gpp_x(:,:) = gpp_daily_fm(:,:,iisf) 885 veget_force_x(:,:) = veget_fm(:,:,iisf) 886 veget_max_force_x(:,:) = veget_max_fm(:,:,iisf) 887 lai_force_x(:,:) = lai_fm(:,:,iisf) 888 WHERE ( t2m(:) .LT. ZeroCelsius ) 889 precip_snow(:) = precip_rain(:) 890 precip_rain(:) = 0. 891 ELSEWHERE 892 precip_snow(:) = 0. 893 ENDWHERE 797 894 !--- 798 895 !-- scale GPP to new lai and veget_max 799 896 !--- 800 WHERE ( lai_x(:,:) .EQ. 0.0 ) gpp_x(:,:) = 0.0897 WHERE ( lai_x(:,:) .EQ. 0.0 ) gpp_x(:,:) = 0.0 801 898 !-- scale GPP to new LAI 802 WHERE (lai_force_x(:,:) .GT. 0.0 )803 gpp_x(:,:) = gpp_x(:,:)*atan(2.*lai_x(:,:)) &804 & / atan( 2.*MAX(lai_force_x(:,:),0.01))805 ENDWHERE899 WHERE (lai_force_x(:,:) .GT. 0.0 ) 900 gpp_x(:,:) = gpp_x(:,:)*ATAN(2.*lai_x(:,:)) & 901 & /ATAN( 2.*MAX(lai_force_x(:,:),0.01)) 902 ENDWHERE 806 903 !-- scale GPP to new veget_max 807 WHERE (veget_max_force_x(:,:) .GT. 0.0 )808 gpp_x(:,:) = gpp_x(:,:)*veget_max_x(:,:)/veget_max_force_x(:,:)809 ENDWHERE904 WHERE (veget_max_force_x(:,:) .GT. 0.0 ) 905 gpp_x(:,:) = gpp_x(:,:)*veget_max_x(:,:)/veget_max_force_x(:,:) 906 ENDWHERE 810 907 !--- 811 908 !-- number crunching 812 909 !--- 813 CALL intsurf_time( itau, date0, dtradia ) 814 ldrestart_read = .FALSE. 815 ldrestart_write = .FALSE. 816 CALL slowproc_main & 910 ldrestart_read = .FALSE. 911 ldrestart_write = .FALSE. 912 CALL slowproc_main & 817 913 & (itau, kjpij, kjpindex, dt_force, date0, & 818 914 & ldrestart_read, ldrestart_write, .FALSE., .TRUE., & … … 822 918 & deadleaf_cover, assim_param_x, lai_x, height_x, veget_x, & 823 919 & frac_nobio, veget_max_x, totfrac_nobio, qsintmax_x, & 824 & rest_id_sec, hist_id_sec, hist_id_sec2, rest_id_sto, hist_id_sto , hist_id_stom_IPCC, co2_flux, &825 !!>> DS add 01/07/2011 correction due to some modifications in the trunk 826 & fco2_lu)827 day_counter = one_day - dt_force920 & rest_id_sec, hist_id_sec, hist_id_sec2, rest_id_sto, hist_id_stom, hist_id_stom_IPCC, co2_flux, fco2_lu) 921 day_counter = one_day - dt_force 922 ! 923 CALL intsurf_time( itau+itau_step, date0, dtradia ) 828 924 ENDDO ! end of the time loop 829 925 !- 830 itau = itau -itau_step926 itau = itau -itau_step 831 927 !- 832 928 ! write restart files 833 929 !- 930 IF (is_root_prc) THEN 834 931 ! first, read and write variables that are not managed otherwise 835 taboo_vars = & 836 & '$lat$ $lon$ $lev$ $veget_year$ '// & 837 & '$height$ $day_counter$ $veget$ $veget_max$ $frac_nobio$ '// & 838 & '$lai$ $soiltype_frac$ $clay_frac$ '// & 839 & '$nav_lon$ $nav_lat$ $nav_lev$ $time$ $time_steps$' 840 !- 841 CALL ioget_vname(rest_id_sec, nbvar, varnames) 842 !!$!- 843 !!$! read and write some special variables (1D or variables that we need) 844 !!$!- 845 !!$ var_name = 'day_counter' 846 !!$ CALL restget (rest_id_sto, var_name, 1, 1, 1, itau_dep, .TRUE., xtmp) 847 !!$ CALL restput (rest_id_sto, var_name, 1, 1, 1, itau_dep, xtmp) 848 !!$!- 849 !!$ var_name = 'dt_days' 850 !!$ CALL restget (rest_id_sto, var_name, 1, 1, 1, itau_dep, .TRUE., xtmp) 851 !!$ CALL restput (rest_id_sto, var_name, 1, 1, 1, itau_dep, xtmp) 852 !!$!- 853 !!$ var_name = 'date' 854 !!$ CALL restget (rest_id_sto, var_name, 1, 1, 1, itau_dep, .TRUE., xtmp) 855 !!$ CALL restput (rest_id_sto, var_name, 1, 1, 1, itau_dep, xtmp) 856 !- 857 DO iv = 1, nbvar 932 taboo_vars = & 933 & '$lat$ $lon$ $lev$ $veget_year$ '// & 934 & '$height$ $day_counter$ $veget$ $veget_max$ $frac_nobio$ '// & 935 & '$lai$ $soiltype_frac$ $clay_frac$ '// & 936 & '$nav_lon$ $nav_lat$ $nav_lev$ $time$ $time_steps$' 937 !- 938 CALL ioget_vname(rest_id_sec, nbvar, varnames) 939 !- 940 DO iv = 1, nbvar 858 941 !-- check if the variable is to be written here 859 IF (INDEX( taboo_vars,'$'//TRIM(varnames(iv))//'$') .EQ. 0 ) THEN942 IF (INDEX( taboo_vars,'$'//TRIM(varnames(iv))//'$') .EQ. 0 ) THEN 860 943 !---- get variable dimensions, especially 3rd dimension 861 CALL ioget_vdim &862 & (rest_id_sec, varnames(iv), varnbdim_max, varnbdim, vardims)863 l1d = ALL(vardims(1:varnbdim) .EQ. 1)864 ALLOCATE(var_3d(kjpindex,vardims(3)),stat=ier)865 IF (ier .NE. 0) STOP 'ALLOCATION PROBLEM'944 CALL ioget_vdim & 945 & (rest_id_sec, varnames(iv), varnbdim_max, varnbdim, vardims) 946 l1d = ALL(vardims(1:varnbdim) .EQ. 1) 947 ALLOCATE(var_3d(nbp_glo,vardims(3)),stat=ier) 948 IF (ier .NE. 0) STOP 'ALLOCATION PROBLEM' 866 949 !---- read it 867 IF (l1d) THEN868 CALL restget (rest_id_sec,TRIM(varnames(iv)), &869 870 ELSE871 CALL restget (rest_id_sec,TRIM(varnames(iv)), &872 kjpindex,vardims(3),1,itau_dep,.TRUE.,var_3d, &873 "gather",kjpindex,indices)874 ENDIF950 IF (l1d) THEN 951 CALL restget (rest_id_sec,TRIM(varnames(iv)), & 952 1,vardims(3),1,itau_dep,.TRUE.,var_3d) 953 ELSE 954 CALL restget (rest_id_sec,TRIM(varnames(iv)), & 955 nbp_glo,vardims(3),1,itau_dep,.TRUE.,var_3d, & 956 "gather",nbp_glo,indices_g) 957 ENDIF 875 958 !---- write it 876 IF (l1d) THEN 877 CALL restput (rest_id_sec,TRIM(varnames(iv)), & 878 1,vardims(3),1,itau,var_3d) 879 ELSE 880 CALL restput (rest_id_sec,TRIM(varnames(iv)), & 881 kjpindex,vardims(3),1,itau,var_3d, & 882 'scatter',kjpindex,indices) 883 ENDIF 884 DEALLOCATE (var_3d) 885 ENDIF 886 ENDDO 959 IF (l1d) THEN 960 CALL restput (rest_id_sec,TRIM(varnames(iv)), & 961 1,vardims(3),1,itau,var_3d) 962 ELSE 963 CALL restput (rest_id_sec,TRIM(varnames(iv)), & 964 nbp_glo,vardims(3),1,itau,var_3d, & 965 'scatter',nbp_glo,indices_g) 966 ENDIF 967 DEALLOCATE (var_3d) 968 ENDIF 969 ENDDO 970 ENDIF 971 CALL barrier_para() 972 887 973 ! call slowproc to write restart files 888 974 ldrestart_read = .FALSE. … … 897 983 & deadleaf_cover, assim_param_x, lai_x, height_x, veget_x, & 898 984 & frac_nobio, veget_max_x, totfrac_nobio, qsintmax_x, & 899 & rest_id_sec, hist_id_sec, hist_id_sec2, rest_id_sto, hist_id_sto, hist_id_stom_IPCC, co2_flux, & 900 !!>> DS add 01/07/2011 correction due to some modifications in the trunk 901 & fco2_lu) 902 985 & rest_id_sec, hist_id_sec, hist_id_sec2, rest_id_sto, hist_id_stom, hist_id_stom_IPCC, co2_flux, fco2_lu) 903 986 !- 904 987 ! close files 905 988 !- 906 CALL restclo 989 IF (is_root_prc) & 990 CALL restclo 907 991 CALL histclo 908 ier = NF90_CLOSE (force_id) 909 !- 910 ! write a new driver restart file with correct time step 911 !- 912 write(*,*) 'teststomate: writing driver restart file with correct time step.' 913 dri_restname_in = 'driver_start.nc' 914 CALL getin ('RESTART_FILEIN',dri_restname_in) 915 dri_restname_out = 'driver_restart.nc' 916 CALL getin ('RESTART_FILEOUT',dri_restname_out) 917 CALL SYSTEM & 918 & ('cp '//TRIM(dri_restname_in)//' '//TRIM(dri_restname_out)) 919 !- 920 iret = NF90_OPEN (TRIM(sec_restname_out),NF90_NOWRITE,ncfid) 921 iret = NF90_INQ_VARID (ncfid,'time',v_id) 922 iret = NF90_GET_VAR (ncfid,v_id,r1d) 923 time_sec = r1d(1) 924 iret = NF90_INQ_VARID (ncfid,'time_steps',v_id) 925 iret = NF90_GET_VAR (ncfid,v_id,time_step_sec) 926 iret = NF90_CLOSE (ncfid) 927 !- 928 iret = NF90_OPEN (TRIM(dri_restname_out),NF90_WRITE,ncfid) 929 iret = NF90_INQ_VARID (ncfid,'time',v_id) 930 iret = NF90_GET_VAR (ncfid,v_id,r1d) 931 time_dri = r1d(1) 932 r1d(1) = time_sec 933 iret = NF90_PUT_VAR (ncfid,v_id,r1d) 934 iret = NF90_INQ_VARID (ncfid,'time_steps',v_id) 935 iret = NF90_GET_VAR (ncfid,v_id,time_step_dri) 936 iret = NF90_PUT_VAR (ncfid,v_id,time_step_sec) 937 iret = NF90_INQ_VARID (ncfid,'julian',v_id) 938 iret = NF90_GET_VAR (ncfid,v_id,r1d) 939 julian = r1d(1) 940 djulian = (time_step_sec-time_step_dri)*dtradia/one_day 941 julian = julian & 942 & +djulian-FLOAT(INT((julian+djulian)/one_year))*one_year 943 r1d(1) = julian 944 iret = NF90_PUT_VAR (ncfid,v_id,r1d) 945 iret = NF90_CLOSE (ncfid) 946 947 CALL getin_dump 992 993 IF (is_root_prc) & 994 ier = NF90_CLOSE (forcing_id) 995 996 IF (is_root_prc) THEN 997 CALL getin_dump() 998 IF ( debug ) WRITE(numout,*) 'REST CLOSED' 999 ENDIF 1000 #ifdef CPP_PARA 1001 CALL MPI_FINALIZE(ier) 1002 #endif 1003 WRITE(numout,*) "End of teststomate." 1004 948 1005 !--------------- 949 1006 END PROGRAM teststomate
Note: See TracChangeset
for help on using the changeset viewer.
