PROGRAM teststomate !--------------------------------------------------------------------- ! This program reads a forcing file for STOMATE ! which was created with STOMATE coupled to SECHIBA. ! It then integrates STOMATE for a ! given number of years using this forcing. ! The GPP is scaled to the new calculated LAI... ! Nothing is done for soil moisture which should ! actually evolve with the vegetation. !--------------------------------------------------------------------- !- IPSL (2006) !- This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC USE netcdf !- USE defprec USE constantes USE pft_parameters USE stomate_data USE ioipsl USE grid USE slowproc USE stomate USE intersurf, ONLY: stom_define_history , intsurf_time USE parallel !- IMPLICIT NONE !- ! Declarations !- INTEGER(i_std) :: vegax_id INTEGER(i_std) :: kjpij,kjpindex REAL(r_std) :: dtradia,dt_force INTEGER(i_std),DIMENSION(:),ALLOCATABLE :: indices INTEGER(i_std),DIMENSION(:),ALLOCATABLE :: indexveg REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: soiltype, veget_x REAL(r_std),DIMENSION(:),ALLOCATABLE :: totfrac_nobio REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: frac_nobio REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: veget_max_x REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: lai_x REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: veget_force_x REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: veget_max_force_x REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: lai_force_x REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: lon,lat REAL(r_std),DIMENSION(:),ALLOCATABLE :: t2m,t2m_min,temp_sol REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: soiltemp,soilhum REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: humrel_x REAL(r_std),DIMENSION(:),ALLOCATABLE :: litterhum REAL(r_std),DIMENSION(:),ALLOCATABLE :: precip_rain,precip_snow REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: gpp_x REAL(r_std),DIMENSION(:),ALLOCATABLE :: deadleaf_cover REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: assim_param_x REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: height_x REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: qsintmax_x REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: co2_flux !!$ >> DS add for externalised version REAL(r_std),DIMENSION(:),ALLOCATABLE :: hist_PFTaxis !!$ >> DS !!$ >> Add for correct 1.9.5.1 version 01/07/2011 ! TYPE(control_type), SAVE :: control_flags REAL(r_std),DIMENSION(:),ALLOCATABLE :: fco2_lu !!$>> DS INTEGER :: ier,iret INTEGER :: ncfid LOGICAL :: a_er CHARACTER(LEN=80) :: & & dri_restname_in,dri_restname_out, & & sec_restname_in,sec_restname_out, & & sto_restname_in,sto_restname_out, stom_histname INTEGER(i_std) :: iim,jjm INTEGER(i_std),PARAMETER :: llm = 1 REAL(r_std),DIMENSION(llm) :: lev REAL(r_std) :: dt_files INTEGER(i_std) :: itau_dep,itau,itau_len,itau_step REAL(r_std) :: date0 INTEGER(i_std) :: rest_id_sec,rest_id_sto INTEGER(i_std) :: hist_id_sec,hist_id_sec2,hist_id_sto,hist_id_stom_IPCC CHARACTER(LEN=30) :: time_str REAL :: hist_days_stom,hist_dt_stom REAL(r_std),DIMENSION(10) :: hist_pool_10axis REAL(r_std),DIMENSION(100) :: hist_pool_100axis REAL(r_std),DIMENSION(11) :: hist_pool_11axis REAL(r_std),DIMENSION(101) :: hist_pool_101axis INTEGER :: hist_PFTaxis_id,hori_id INTEGER :: hist_pool_10axis_id INTEGER :: hist_pool_100axis_id INTEGER :: hist_pool_11axis_id INTEGER :: hist_pool_101axis_id INTEGER :: hist_level INTEGER,PARAMETER :: max_hist_level = 10 INTEGER(i_std) :: i,j,iv,id CHARACTER*80 :: var_name CHARACTER(LEN=40),DIMENSION(10) :: fluxop LOGICAL :: ldrestart_read,ldrestart_write LOGICAL :: l1d INTEGER(i_std),PARAMETER :: nbvarmax=200 INTEGER(i_std) :: nbvar CHARACTER*50,DIMENSION(nbvarmax) :: varnames INTEGER(i_std) :: varnbdim INTEGER(i_std),PARAMETER :: varnbdim_max=20 INTEGER,DIMENSION(varnbdim_max) :: vardims CHARACTER*200 :: taboo_vars REAL(r_std),DIMENSION(1) :: xtmp INTEGER :: nsfm,nsft INTEGER :: iisf INTEGER(i_std) :: max_totsize,totsize_1step INTEGER(i_std),DIMENSION(0:2) :: ifirst,ilast INTEGER(i_std) :: iblocks,nblocks INTEGER,PARAMETER :: ndm = 10 INTEGER,DIMENSION(ndm) :: start,count INTEGER :: ndim,v_id INTEGER :: force_id CHARACTER(LEN=100) :: forcing_name REAL :: x REAL(r_std),DIMENSION(:),ALLOCATABLE :: x_indices REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: x_neighbours REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: var_3d !- REAL(r_std) :: time_sec,time_step_sec REAL(r_std) :: time_dri,time_step_dri REAL(r_std),DIMENSION(1) :: r1d REAL(r_std) :: julian,djulian ! REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: soiltype !- ! the following variables contain the forcing data !- REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: clay_fm REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: humrel_x_fm REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: litterhum_fm REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: t2m_fm REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: t2m_min_fm REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: temp_sol_fm REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: soiltemp_fm REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: soilhum_fm REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: precip_fm REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: gpp_x_fm REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: veget_force_x_fm REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: veget_max_force_x_fm REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: lai_force_x_fm INTEGER(i_std),ALLOCATABLE,SAVE,DIMENSION(:) :: isf LOGICAL,ALLOCATABLE,SAVE,DIMENSION(:) :: nf_written INTEGER(i_std),ALLOCATABLE,SAVE,DIMENSION(:) :: nf_cumul INTEGER(i_std) :: ji,jv !--------------------------------------------------------------------- !- No parallelisation yet in teststomate ! #ifdef CPP_PARA CALL ipslerr (3,'teststomate', & & 'You try to run testsomate compiled with parallelisation. (CPP_PARA key)', & & 'But it wasn''t programmed yet and teststomate will stop.','You must compiled it without CPP_PARA key.') #endif !- ! calendar !- CALL ioconf_calendar ('noleap') CALL ioget_calendar (one_year,one_day) !- ! open STOMATE's forcing file to read some basic info !- forcing_name = 'stomate_forcing.nc' CALL getin ('STOMATE_FORCING_NAME',forcing_name) iret = NF90_OPEN (TRIM(forcing_name),NF90_NOWRITE,force_id) IF (iret /= NF90_NOERR) THEN CALL ipslerr (3,'teststomate', & & 'Could not open file : ', & & forcing_name,'(Do you have forget it ?)') ENDIF ier = NF90_GET_ATT (force_id,NF90_GLOBAL,'dtradia',dtradia) ier = NF90_GET_ATT (force_id,NF90_GLOBAL,'dt_slow',dt_force) ier = NF90_GET_ATT (force_id,NF90_GLOBAL,'nsft',x) nsft = NINT(x) ier = NF90_GET_ATT (force_id,NF90_GLOBAL,'kjpij',x) kjpij = NINT(x) ier = NF90_GET_ATT (force_id,NF90_GLOBAL,'kjpindex',x) kjpindex = NINT(x) CALL init_grid( kjpindex ) !- write(*,*) 'ATTENTION',dtradia,dt_force !!$ >> DS added for the externalised version !!$ we need to know the pft parameters values used for stomate ! 1. Read the number of PFTs CALL getin('NVM',nvm) ! 2. Allocate and read the pft parameters stomate specific CALL pft_parameters_main CALL getin_stomate_pft_parameters !!$ DS !- ! allocate arrays !- a_er = .FALSE. !!$ >> DS added for the externalised version ALLOCATE (hist_PFTaxis(nvm),stat=ier) a_er = a_er .OR. (ier.NE.0) ! 01/07/2011 allocate fco2_lu ALLOCATE (fco2_lu(kjpindex),stat=ier) a_er = a_er .OR. (ier.NE.0) !!$ end add DS ALLOCATE (indices(kjpindex),stat=ier) a_er = a_er .OR. (ier.NE.0) ALLOCATE (indexveg(kjpindex*nvm), stat=ier) a_er = a_er .OR. (ier.NE.0) ALLOCATE (soiltype(kjpindex,nstm),stat=ier) a_er = a_er .OR. (ier.NE.0) ALLOCATE (veget_x(kjpindex,nvm),stat=ier) a_er = a_er .OR. (ier.NE.0) ALLOCATE (totfrac_nobio(kjpindex),stat=ier) a_er = a_er .OR. (ier.NE.0) ALLOCATE (frac_nobio(kjpindex,nnobio),stat=ier) a_er = a_er .OR. (ier.NE.0) ALLOCATE (veget_max_x(kjpindex,nvm),stat=ier) a_er = a_er .OR. (ier.NE.0) ALLOCATE (lai_x(kjpindex,nvm),stat=ier) a_er = a_er .OR. (ier.NE.0) ALLOCATE (veget_force_x(kjpindex,nvm),stat=ier) a_er = a_er .OR. (ier.NE.0) ALLOCATE (veget_max_force_x(kjpindex,nvm),stat=ier) a_er = a_er .OR. (ier.NE.0) ALLOCATE (lai_force_x(kjpindex,nvm),stat=ier) a_er = a_er .OR. (ier.NE.0) ALLOCATE (t2m(kjpindex),stat=ier) a_er = a_er .OR. (ier.NE.0) ALLOCATE (t2m_min(kjpindex),stat=ier) a_er = a_er .OR. (ier.NE.0) ALLOCATE (temp_sol(kjpindex),stat=ier) a_er = a_er .OR. (ier.NE.0) ALLOCATE (soiltemp(kjpindex,nbdl),stat=ier) a_er = a_er .OR. (ier.NE.0) ALLOCATE (soilhum(kjpindex,nbdl),stat=ier) a_er = a_er .OR. (ier.NE.0) ALLOCATE (humrel_x(kjpindex,nvm),stat=ier) a_er = a_er .OR. (ier.NE.0) ALLOCATE (litterhum(kjpindex),stat=ier) a_er = a_er .OR. (ier.NE.0) ALLOCATE (precip_rain(kjpindex),stat=ier) a_er = a_er .OR. (ier.NE.0) ALLOCATE (precip_snow(kjpindex),stat=ier) a_er = a_er .OR. (ier.NE.0) ALLOCATE (gpp_x(kjpindex,nvm),stat=ier) a_er = a_er .OR. (ier.NE.0) ALLOCATE (deadleaf_cover(kjpindex),stat=ier) a_er = a_er .OR. (ier.NE.0) ALLOCATE (assim_param_x(kjpindex,nvm,npco2),stat=ier) a_er = a_er .OR. (ier.NE.0) ALLOCATE (height_x(kjpindex,nvm),stat=ier) a_er = a_er .OR. (ier.NE.0) ALLOCATE (qsintmax_x(kjpindex,nvm),stat=ier) a_er = a_er .OR. (ier.NE.0) ALLOCATE (co2_flux(kjpindex,nvm),stat=ier) a_er = a_er .OR. (ier.NE.0) IF ( a_er ) STOP 'PROBLEM WITH ALLOCATION' ! ! prepare forcing ! max_totsize = 50 CALL getin ('STOMATE_FORCING_MEMSIZE',max_totsize) max_totsize = max_totsize * 1000000 totsize_1step = SIZE(soiltype(:,3))*KIND(soiltype(:,3)) + & SIZE(humrel_x)*KIND(humrel_x) + & SIZE(litterhum)*KIND(litterhum) + & SIZE(t2m)*KIND(t2m) + & SIZE(t2m_min)*KIND(t2m_min) + & SIZE(temp_sol)*KIND(temp_sol) + & SIZE(soiltemp)*KIND(soiltemp) + & SIZE(soilhum)*KIND(soilhum) + & SIZE(precip_rain)*KIND(precip_rain) + & SIZE(precip_snow)*KIND(precip_snow) + & SIZE(gpp_x)*KIND(gpp_x) + & SIZE(veget_force_x)*KIND(veget_force_x) + & SIZE(veget_max_force_x)*KIND(veget_max_force_x) + & SIZE(lai_force_x)*KIND(lai_force_x) ! total number of forcing steps nsft = INT(one_year/(dt_force/one_day)) ! number of forcing steps in memory nsfm = MIN(nsft,MAX(1,NINT(FLOAT(max_totsize)/FLOAT(totsize_1step)))) !- WRITE(numout,*) 'Offline forcing of Stomate:' WRITE(numout,*) ' Total number of forcing steps:',nsft WRITE(numout,*) ' Number of forcing steps in memory:',nsfm !- a_er = .FALSE. ALLOCATE (clay_fm(kjpindex,nsfm),stat=ier) a_er = a_er.OR.(ier.NE.0) ALLOCATE (humrel_x_fm(kjpindex,nvm,nsfm),stat=ier) a_er = a_er.OR.(ier.NE.0) ALLOCATE (litterhum_fm(kjpindex,nsfm),stat=ier) a_er = a_er.OR.(ier.NE.0) ALLOCATE (t2m_fm(kjpindex,nsfm),stat=ier) a_er = a_er.OR.(ier.NE.0) ALLOCATE (t2m_min_fm(kjpindex,nsfm),stat=ier) a_er = a_er.OR.(ier.NE.0) ALLOCATE (temp_sol_fm(kjpindex,nsfm),stat=ier) a_er = a_er.OR.(ier.NE.0) ALLOCATE (soiltemp_fm(kjpindex,nbdl,nsfm),stat=ier) a_er = a_er.OR.(ier.NE.0) ALLOCATE (soilhum_fm(kjpindex,nbdl,nsfm),stat=ier) a_er = a_er.OR.(ier.NE.0) ALLOCATE (precip_fm(kjpindex,nsfm),stat=ier) a_er = a_er.OR.(ier.NE.0) ALLOCATE (gpp_x_fm(kjpindex,nvm,nsfm),stat=ier) a_er = a_er.OR.(ier.NE.0) ALLOCATE (veget_force_x_fm(kjpindex,nvm,nsfm),stat=ier) a_er = a_er.OR.(ier.NE.0) ALLOCATE (veget_max_force_x_fm(kjpindex,nvm,nsfm),stat=ier) a_er = a_er.OR.(ier.NE.0) ALLOCATE (lai_force_x_fm(kjpindex,nvm,nsfm),stat=ier) a_er = a_er.OR.(ier.NE.0) ALLOCATE (isf(nsfm),stat=ier) a_er = a_er.OR.(ier.NE.0) ALLOCATE (nf_written(nsft),stat=ier) a_er = a_er.OR.(ier.NE.0) ALLOCATE (nf_cumul(nsft),stat=ier) a_er = a_er.OR.(ier.NE.0) IF (a_er) THEN STOP 'stomate: error in memory allocation for forcing data' ENDIF ! ensure that we read all new forcing states iisf = nsfm ! initialize that table that contains the indices ! of the forcing states that will be in memory isf(:) = (/ (i,i=1,nsfm) /) !- ! read info about grids !- contfrac(:) = 1.0 !- ALLOCATE (x_indices(kjpindex),stat=ier) ier = NF90_INQ_VARID (force_id,'index',v_id) ier = NF90_GET_VAR (force_id,v_id,x_indices) indices(:) = NINT(x_indices(:)) DEALLOCATE (x_indices) !- DO ji=1,kjpindex DO jv=1,nvm indexveg((jv-1)*kjpindex+ji) = indices(ji)+(jv-1)*kjpij ENDDO ENDDO !- ier = NF90_INQ_VARID (force_id,'lalo',v_id) ier = NF90_GET_VAR (force_id,v_id,lalo) !- ALLOCATE (x_neighbours(kjpindex,8),stat=ier) ier = NF90_INQ_VARID (force_id,'neighbours',v_id) ier = NF90_GET_VAR (force_id,v_id,x_neighbours) neighbours(:,:) = NINT(x_neighbours(:,:)) DEALLOCATE (x_neighbours) !- ier = NF90_INQ_VARID (force_id,'resolution',v_id) ier = NF90_GET_VAR (force_id,v_id,resolution) !- ier = NF90_INQ_VARID (force_id,'contfrac',v_id) ier = NF90_GET_VAR (force_id,v_id,contfrac) !- ! activate CO2, STOMATE, but not sechiba !- control%river_routing = .FALSE. control%hydrol_cwrr = .FALSE. control%ok_sechiba = .FALSE. ! control%stomate_watchout = .TRUE. control%ok_co2 = .TRUE. control%ok_stomate = .TRUE. !!$ >> DS now we search the values for the scalar parameters in some .def files ! 07/2011 CALL activate_sub_models(control%ok_sechiba,control%river_routing,control%ok_stomate) CALL veget_config ! CALL getin_co2_parameters CALL getin_stomate_parameters !!$ >> DS !- ! is DGVM activated? !- control%ok_dgvm = .FALSE. CALL getin('STOMATE_OK_DGVM',control%ok_dgvm) WRITE(*,*) 'LPJ is activated: ',control%ok_dgvm !!$ >> DS for the externalisation reading the scalar parameters when ok_dgvm is set to true IF (control%ok_dgvm) THEN CALL getin_dgvm_parameters ENDIF !!$ >> DS !- ! restart files !- ! Sechiba's restart files sec_restname_in = 'sechiba_start.nc' CALL getin('SECHIBA_restart_in',sec_restname_in) WRITE(*,*) 'SECHIBA INPUT RESTART_FILE: ',TRIM(sec_restname_in) IF ( TRIM(sec_restname_in) .EQ. 'NONE' ) THEN STOP 'Need a restart file for Sechiba' ENDIF sec_restname_out = 'sechiba_restart.nc' CALL getin('SECHIBA_rest_out',sec_restname_out) WRITE(*,*) 'SECHIBA OUTPUT RESTART_FILE: ',TRIM(sec_restname_out) ! Stomate's restart files sto_restname_in = 'stomate_start.nc' CALL getin ('STOMATE_RESTART_FILEIN',sto_restname_in) WRITE(*,*) 'STOMATE INPUT RESTART_FILE: ',TRIM(sto_restname_in) sto_restname_out = 'stomate_restart.nc' CALL getin ('STOMATE_RESTART_FILEOUT',sto_restname_out) WRITE(*,*) 'STOMATE OUTPUT RESTART_FILE: ',TRIM(sto_restname_out) !- ! We need to know iim, jjm. ! Get them from the restart files themselves. !- iret = NF90_OPEN (sec_restname_in,NF90_NOWRITE,ncfid) IF (iret /= NF90_NOERR) THEN CALL ipslerr (3,'teststomate', & & 'Could not open file : ', & & sec_restname_in,'(Do you have forget it ?)') ENDIF iret = NF90_INQUIRE_DIMENSION (ncfid,1,len=iim) iret = NF90_INQUIRE_DIMENSION (ncfid,2,len=jjm) iret = NF90_CLOSE (ncfid) ! Allocate longitudes and latitudes ALLOCATE (lon(iim,jjm),stat=ier) a_er = a_er.OR.(ier.NE.0) ALLOCATE (lat(iim,jjm),stat=ier) a_er = a_er.OR.(ier.NE.0) lon(:,:) = 0.0 lat(:,:) = 0.0 lev(1) = 0.0 !- CALL restini & & (sec_restname_in, iim, jjm, lon, lat, llm, lev, & & sec_restname_out, itau_dep, date0, dt_files, rest_id_sec) !- CALL restini & & (sto_restname_in, iim, jjm, lon, lat, llm, lev, & & sto_restname_out, itau_dep, date0, dt_files, rest_id_sto) !- IF ( dt_files .NE. dtradia ) THEN WRITE(*,*) 'dt_files',dt_files WRITE(*,*) 'dtradia',dtradia STOP 'PROBLEM with time steps.' ENDIF !- ! a time step for STOMATE corresponds to itau_step timesteps in SECHIBA !- itau_step = NINT(dt_force/dtradia) ! CALL ioconf_startdate(date0) CALL intsurf_time( itau_dep, date0, dtradia ) !- ! Length of integration !- WRITE(time_str,'(a)') '1Y' CALL getin ('TIME_LENGTH', time_str) ! transform into itau CALL tlen2itau(time_str, dt_files, date0, itau_len) ! itau_len*dtradia must be a multiple of dt_force itau_len = NINT( MAX(1.,FLOAT(NINT(itau_len*dtradia/dt_force))) & & *dt_force/dtradia) !- ! set up STOMATE history file !- !Config Key = STOMATE_OUTPUT_FILE !Config Desc = Name of file in which STOMATE's output is going !Config to be written !Config Def = stomate_history.nc !Config Help = This file is going to be created by the model !Config and will contain the output from the model. !Config This file is a truly COADS compliant netCDF file. !Config It will be generated by the hist software from !Config the IOIPSL package. !- stom_histname='stomate_history.nc' CALL getin ('STOMATE_OUTPUT_FILE', stom_histname) WRITE(*,*) 'STOMATE_OUTPUT_FILE', TRIM(stom_histname) !- !Config Key = STOMATE_HIST_DT !Config Desc = STOMATE history time step (d) !Config Def = 10. !Config Help = Time step of the STOMATE history file !- hist_days_stom = 10. CALL getin ('STOMATE_HIST_DT', hist_days_stom) IF ( hist_days_stom == -1. ) THEN hist_dt_stom = -1. WRITE(numout,*) 'output frequency for STOMATE history file (d): one month.' ELSE hist_dt_stom = NINT( hist_days_stom ) * one_day WRITE(numout,*) 'output frequency for STOMATE history file (d): ', & hist_dt_stom/one_day ENDIF !- ! initialize CALL histbeg(stom_histname, iim, lon, jjm, lat, 1, iim, 1, jjm, & & itau_dep, date0, dt_files, hori_id, hist_id_sto) ! define PFT axis hist_PFTaxis = (/ ( REAL(i,r_std), i=1,nvm ) /) ! declare this axis CALL histvert (hist_id_sto, 'PFT', 'Plant functional type', & & '-', nvm, hist_PFTaxis, hist_PFTaxis_id) !- define Pool_10 axis hist_pool_10axis = (/ ( REAL(i,r_std), i=1,10 ) /) !- declare this axis CALL histvert (hist_id_sto, 'P10', 'Pool 10 years', & & '-', 10, hist_pool_10axis, hist_pool_10axis_id) !- define Pool_100 axis hist_pool_100axis = (/ ( REAL(i,r_std), i=1,100 ) /) !- declare this axis CALL histvert (hist_id_sto, 'P100', 'Pool 100 years', & & '-', 100, hist_pool_100axis, hist_pool_100axis_id) !- define Pool_11 axis hist_pool_11axis = (/ ( REAL(i,r_std), i=1,11 ) /) !- declare this axis CALL histvert (hist_id_sto, 'P11', 'Pool 10 years + 1', & & '-', 11, hist_pool_11axis, hist_pool_11axis_id) !- define Pool_101 axis hist_pool_101axis = (/ ( REAL(i,r_std), i=1,101 ) /) !- declare this axis CALL histvert (hist_id_sto, 'P101', 'Pool 100 years + 1', & & '-', 101, hist_pool_101axis, hist_pool_101axis_id) ! define STOMATE history file CALL stom_define_history (hist_id_sto, nvm, iim, jjm, & & dt_files, hist_dt_stom, hori_id, hist_PFTaxis_id, & & hist_pool_10axis_id, hist_pool_100axis_id, & & hist_pool_11axis_id, hist_pool_101axis_id) ! end definition CALL histend(hist_id_sto) hist_id_sec = -1 hist_id_sec2 = -1 hist_id_stom_IPCC = -1 !- ! read some variables we need from SECHIBA's restart file !- CALL ioget_expval(val_exp) !- ! first call of slowproc to initialize variables !- itau = itau_dep ldrestart_read = .TRUE. ldrestart_write = .FALSE. !- !MM Problem here with dpu which depends on soil type DO jv = 1, nbdl-1 ! first 2.0 is dpu ! second 2.0 is average diaglev(jv) = 2.0/(2**(nbdl-1) -1) * ( ( 2**(jv-1) -1) + ( 2**(jv) -1) ) / 2.0 ENDDO diaglev(nbdl) = 2.0 !- ! For sequential use only, we must initialize data_para : ! CALL init_para(.FALSE.) ! CALL init_data_para(iim,jjm,kjpindex,indices) ! !- global index index_g is the index_l of root proc IF (is_root_prc) index_g(:)=indices(1:kjpindex) !- ! >> DS : getin the value of slowproc parameters CALL getin('CLAYFRACTION_DEFAULT',clayfraction_default) CALL getin('SOILTYPE_DEFAULT',soiltype_default) ! >> DS CALL slowproc_main & & (itau, kjpij, kjpindex, dt_force, date0, & & ldrestart_read, ldrestart_write, .FALSE., .TRUE., & & indices, indexveg, lalo, neighbours, resolution, contfrac, soiltype, & & t2m, t2m_min, temp_sol, soiltemp, & & humrel_x, soilhum, litterhum, precip_rain, precip_snow, gpp_x, & & deadleaf_cover, assim_param_x, lai_x, height_x, veget_x, & & frac_nobio, veget_max_x, totfrac_nobio, qsintmax_x, & & rest_id_sec, hist_id_sec, hist_id_sec2, rest_id_sto, hist_id_sto, hist_id_stom_IPCC, co2_flux, & !!>> DS add 01/07/2011 correction due to some modifications in the trunk & fco2_lu) ! correct date day_counter = one_day - dt_force date=1 !- ! time loop !- DO itau = itau_dep+itau_step,itau_dep+itau_len,itau_step !-- next forcing state iisf = iisf+1 !--- IF (iisf .GT. nsfm) THEN !----- !---- we have to read new forcing states from the file !----- !---- determine blocks of forcing states that are contiguous in memory !----- nblocks = 0 ifirst(:) = 1 ilast(:) = 1 !----- DO iisf=1,nsfm IF ( (nblocks .NE. 0) ) THEN IF ( (isf(iisf) .EQ. isf(ilast(nblocks))+1) ) THEN !-------- element is contiguous with last element found ilast(nblocks) = iisf ELSE !-------- found first element of new block nblocks = nblocks+1 IF (nblocks .GT. nsfm) THEN ! IF (nblocks .GT. 2) THEN STOP 'Problem in teststomate' ENDIF ifirst(nblocks) = iisf ilast(nblocks) = iisf ENDIF ELSE !-------- found first element of new block nblocks = nblocks+1 IF (nblocks .GT. nsfm) THEN ! IF (nblocks .GT. 2) THEN STOP 'Problem in teststomate' ENDIF ifirst(nblocks) = iisf ilast(nblocks) = iisf ENDIF ENDDO !----- DO iblocks=1,nblocks IF (ifirst(iblocks) .NE. ilast(iblocks)) THEN ndim = 2; start(:) = 1; start(ndim) = isf(ifirst(iblocks)); count(1:ndim) = SHAPE(clay_fm) count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1 ier = NF90_INQ_VARID (force_id,'clay',v_id) a_er = a_er.OR.(ier.NE.0) ier = NF90_GET_VAR (force_id,v_id, & & clay_fm(:,ifirst(iblocks):ilast(iblocks)), & & start=start(1:ndim), count=count(1:ndim)) a_er = a_er.OR.(ier.NE.0) !--------- ndim = 3; start(:) = 1; start(ndim) = isf(ifirst(iblocks)); count(1:ndim) = SHAPE(humrel_x_fm) count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1 ier = NF90_INQ_VARID (force_id,'humrel',v_id) a_er = a_er.OR.(ier.NE.0) ier = NF90_GET_VAR (force_id,v_id, & & humrel_x_fm(:,:,ifirst(iblocks):ilast(iblocks)), & & start=start(1:ndim), count=count(1:ndim)) a_er = a_er.OR.(ier.NE.0) !--------- ndim = 2; start(:) = 1; start(ndim) = isf(ifirst(iblocks)); count(1:ndim) = SHAPE(litterhum_fm) count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1 ier = NF90_INQ_VARID (force_id,'litterhum',v_id) a_er = a_er.OR.(ier.NE.0) ier = NF90_GET_VAR (force_id,v_id, & & litterhum_fm(:,ifirst(iblocks):ilast(iblocks)), & & start=start(1:ndim), count=count(1:ndim)) a_er = a_er.OR.(ier.NE.0) !--------- ndim = 2; start(:) = 1; start(ndim) = isf(ifirst(iblocks)); count(1:ndim) = SHAPE(t2m_fm) count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1 ier = NF90_INQ_VARID (force_id,'t2m',v_id) a_er = a_er.OR.(ier.NE.0) ier = NF90_GET_VAR (force_id,v_id, & & t2m_fm(:,ifirst(iblocks):ilast(iblocks)), & & start=start(1:ndim), count=count(1:ndim)) a_er = a_er.OR.(ier.NE.0) !--------- ndim = 2; start(:) = 1; start(ndim) = isf(ifirst(iblocks)); count(1:ndim) = SHAPE(t2m_min_fm) count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1 ier = NF90_INQ_VARID (force_id,'t2m_min',v_id) a_er = a_er.OR.(ier.NE.0) ier = NF90_GET_VAR (force_id,v_id, & & t2m_min_fm(:,ifirst(iblocks):ilast(iblocks)), & & start=start(1:ndim), count=count(1:ndim)) a_er = a_er.OR.(ier.NE.0) !--------- ndim = 2; start(:) = 1; start(ndim) = isf(ifirst(iblocks)); count(1:ndim) = SHAPE(temp_sol_fm) count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1 ier = NF90_INQ_VARID (force_id,'tsurf',v_id) a_er = a_er.OR.(ier.NE.0) ier = NF90_GET_VAR (force_id,v_id, & & temp_sol_fm(:,ifirst(iblocks):ilast(iblocks)), & & start=start(1:ndim), count=count(1:ndim)) a_er = a_er.OR.(ier.NE.0) !--------- ndim = 3; start(:) = 1; start(ndim) = isf(ifirst(iblocks)); count(1:ndim) = SHAPE(soiltemp_fm) count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1 ier = NF90_INQ_VARID (force_id,'tsoil',v_id) a_er = a_er.OR.(ier.NE.0) ier = NF90_GET_VAR (force_id,v_id, & & soiltemp_fm(:,:,ifirst(iblocks):ilast(iblocks)), & & start=start(1:ndim), count=count(1:ndim)) a_er = a_er.OR.(ier.NE.0) !--------- ndim = 3; start(:) = 1; start(ndim) = isf(ifirst(iblocks)); count(1:ndim) = SHAPE(soilhum_fm) count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1 ier = NF90_INQ_VARID (force_id,'soilhum',v_id) a_er = a_er.OR.(ier.NE.0) ier = NF90_GET_VAR (force_id,v_id, & & soilhum_fm(:,:,ifirst(iblocks):ilast(iblocks)), & & start=start(1:ndim), count=count(1:ndim)) a_er = a_er.OR.(ier.NE.0) !--------- ndim = 2; start(:) = 1; start(ndim) = isf(ifirst(iblocks)); count(1:ndim) = SHAPE(precip_fm) count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1 ier = NF90_INQ_VARID (force_id,'precip',v_id) a_er = a_er.OR.(ier.NE.0) ier = NF90_GET_VAR (force_id,v_id, & & precip_fm(:,ifirst(iblocks):ilast(iblocks)), & & start=start(1:ndim), count=count(1:ndim)) a_er = a_er.OR.(ier.NE.0) !--------- ndim = 3; start(:) = 1; start(ndim) = isf(ifirst(iblocks)); count(1:ndim) = SHAPE(gpp_x_fm) count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1 ier = NF90_INQ_VARID (force_id,'gpp',v_id) a_er = a_er.OR.(ier.NE.0) ier = NF90_GET_VAR (force_id,v_id, & & gpp_x_fm(:,:,ifirst(iblocks):ilast(iblocks)), & & start=start(1:ndim), count=count(1:ndim)) a_er = a_er.OR.(ier.NE.0) !--------- ndim = 3; start(:) = 1; start(ndim) = isf(ifirst(iblocks)); count(1:ndim) = SHAPE(veget_force_x_fm) count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1 ier = NF90_INQ_VARID (force_id,'veget',v_id) a_er = a_er.OR.(ier.NE.0) ier = NF90_GET_VAR (force_id,v_id, & & veget_force_x_fm(:,:,ifirst(iblocks):ilast(iblocks)), & & start=start(1:ndim), count=count(1:ndim)) a_er = a_er.OR.(ier.NE.0) !--------- ndim = 3; start(:) = 1; start(ndim) = isf(ifirst(iblocks)); count(1:ndim) = SHAPE(veget_max_force_x_fm) count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1 ier = NF90_INQ_VARID (force_id,'veget_max',v_id) a_er = a_er.OR.(ier.NE.0) ier = NF90_GET_VAR (force_id,v_id, & & veget_max_force_x_fm(:,:,ifirst(iblocks):ilast(iblocks)), & & start=start(1:ndim), count=count(1:ndim)) a_er = a_er.OR.(ier.NE.0) !--------- ndim = 3; start(:) = 1; start(ndim) = isf(ifirst(iblocks)); count(1:ndim) = SHAPE(lai_force_x_fm) count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1 ier = NF90_INQ_VARID (force_id,'lai',v_id) a_er = a_er.OR.(ier.NE.0) ier = NF90_GET_VAR (force_id,v_id, & & lai_force_x_fm(:,:,ifirst(iblocks):ilast(iblocks)), & & start=start(1:ndim), count=count(1:ndim)) a_er = a_er.OR.(ier.NE.0) ENDIF ENDDO !----- !---- determine which forcing states must be read next time !----- isf(1) = isf(nsfm)+1 IF ( isf(1) .GT. nsft ) isf(1) = 1 DO iisf = 2, nsfm isf(iisf) = isf(iisf-1)+1 IF ( isf(iisf) .GT. nsft ) isf(iisf) = 1 ENDDO !---- start again at first forcing state iisf = 1 ENDIF soiltype(:,3) = clay_fm(:,iisf) humrel_x(:,:) = humrel_x_fm(:,:,iisf) litterhum(:) = litterhum_fm(:,iisf) t2m(:) = t2m_fm(:,iisf) t2m_min(:) = t2m_min_fm(:,iisf) temp_sol(:) = temp_sol_fm(:,iisf) soiltemp(:,:) = soiltemp_fm(:,:,iisf) soilhum(:,:) = soilhum_fm(:,:,iisf) precip_rain(:) = precip_fm(:,iisf) gpp_x(:,:) = gpp_x_fm(:,:,iisf) veget_force_x(:,:) = veget_force_x_fm(:,:,iisf) veget_max_force_x(:,:) = veget_max_force_x_fm(:,:,iisf) lai_force_x(:,:) = lai_force_x_fm(:,:,iisf) WHERE ( t2m(:) .LT. ZeroCelsius ) precip_snow(:) = precip_rain(:) precip_rain(:) = 0. ELSEWHERE precip_snow(:) = 0. ENDWHERE !--- !-- scale GPP to new lai and veget_max !--- WHERE ( lai_x(:,:) .EQ. 0.0 ) gpp_x(:,:) = 0.0 !-- scale GPP to new LAI WHERE (lai_force_x(:,:) .GT. 0.0 ) gpp_x(:,:) = gpp_x(:,:)*atan(2.*lai_x(:,:)) & & /atan( 2.*MAX(lai_force_x(:,:),0.01)) ENDWHERE !-- scale GPP to new veget_max WHERE (veget_max_force_x(:,:) .GT. 0.0 ) gpp_x(:,:) = gpp_x(:,:)*veget_max_x(:,:)/veget_max_force_x(:,:) ENDWHERE !--- !-- number crunching !--- CALL intsurf_time( itau, date0, dtradia ) ldrestart_read = .FALSE. ldrestart_write = .FALSE. CALL slowproc_main & & (itau, kjpij, kjpindex, dt_force, date0, & & ldrestart_read, ldrestart_write, .FALSE., .TRUE., & & indices, indexveg, lalo, neighbours, resolution, contfrac, soiltype, & & t2m, t2m_min, temp_sol, soiltemp, & & humrel_x, soilhum, litterhum, precip_rain, precip_snow, gpp_x, & & deadleaf_cover, assim_param_x, lai_x, height_x, veget_x, & & frac_nobio, veget_max_x, totfrac_nobio, qsintmax_x, & & rest_id_sec, hist_id_sec, hist_id_sec2, rest_id_sto, hist_id_sto, hist_id_stom_IPCC, co2_flux, & !!>> DS add 01/07/2011 correction due to some modifications in the trunk & fco2_lu) day_counter = one_day - dt_force ENDDO ! end of the time loop !- itau = itau -itau_step !- ! write restart files !- ! first, read and write variables that are not managed otherwise taboo_vars = & & '$lat$ $lon$ $lev$ $veget_year$ '// & & '$height$ $day_counter$ $veget$ $veget_max$ $frac_nobio$ '// & & '$lai$ $soiltype_frac$ $clay_frac$ '// & & '$nav_lon$ $nav_lat$ $nav_lev$ $time$ $time_steps$' !- CALL ioget_vname(rest_id_sec, nbvar, varnames) !!$!- !!$! read and write some special variables (1D or variables that we need) !!$!- !!$ var_name = 'day_counter' !!$ CALL restget (rest_id_sto, var_name, 1, 1, 1, itau_dep, .TRUE., xtmp) !!$ CALL restput (rest_id_sto, var_name, 1, 1, 1, itau_dep, xtmp) !!$!- !!$ var_name = 'dt_days' !!$ CALL restget (rest_id_sto, var_name, 1, 1, 1, itau_dep, .TRUE., xtmp) !!$ CALL restput (rest_id_sto, var_name, 1, 1, 1, itau_dep, xtmp) !!$!- !!$ var_name = 'date' !!$ CALL restget (rest_id_sto, var_name, 1, 1, 1, itau_dep, .TRUE., xtmp) !!$ CALL restput (rest_id_sto, var_name, 1, 1, 1, itau_dep, xtmp) !- DO iv = 1, nbvar !-- check if the variable is to be written here IF (INDEX( taboo_vars,'$'//TRIM(varnames(iv))//'$') .EQ. 0 ) THEN !---- get variable dimensions, especially 3rd dimension CALL ioget_vdim & & (rest_id_sec, varnames(iv), varnbdim_max, varnbdim, vardims) l1d = ALL(vardims(1:varnbdim) .EQ. 1) ALLOCATE(var_3d(kjpindex,vardims(3)),stat=ier) IF (ier .NE. 0) STOP 'ALLOCATION PROBLEM' !---- read it IF (l1d) THEN CALL restget (rest_id_sec,TRIM(varnames(iv)), & 1,vardims(3),1,itau_dep,.TRUE.,var_3d) ELSE CALL restget (rest_id_sec,TRIM(varnames(iv)), & kjpindex,vardims(3),1,itau_dep,.TRUE.,var_3d, & "gather",kjpindex,indices) ENDIF !---- write it IF (l1d) THEN CALL restput (rest_id_sec,TRIM(varnames(iv)), & 1,vardims(3),1,itau,var_3d) ELSE CALL restput (rest_id_sec,TRIM(varnames(iv)), & kjpindex,vardims(3),1,itau,var_3d, & 'scatter',kjpindex,indices) ENDIF DEALLOCATE (var_3d) ENDIF ENDDO ! call slowproc to write restart files ldrestart_read = .FALSE. ldrestart_write = .TRUE. !- CALL slowproc_main & & (itau, kjpij, kjpindex, dt_force, date0, & & ldrestart_read, ldrestart_write, .FALSE., .TRUE., & & indices, indexveg, lalo, neighbours, resolution, contfrac, soiltype, & & t2m, t2m_min, temp_sol, soiltemp, & & humrel_x, soilhum, litterhum, precip_rain, precip_snow, gpp_x, & & deadleaf_cover, assim_param_x, lai_x, height_x, veget_x, & & frac_nobio, veget_max_x, totfrac_nobio, qsintmax_x, & & rest_id_sec, hist_id_sec, hist_id_sec2, rest_id_sto, hist_id_sto, hist_id_stom_IPCC, co2_flux, & !!>> DS add 01/07/2011 correction due to some modifications in the trunk & fco2_lu) !- ! close files !- CALL restclo CALL histclo ier = NF90_CLOSE (force_id) !- ! write a new driver restart file with correct time step !- write(*,*) 'teststomate: writing driver restart file with correct time step.' dri_restname_in = 'driver_start.nc' CALL getin ('RESTART_FILEIN',dri_restname_in) dri_restname_out = 'driver_restart.nc' CALL getin ('RESTART_FILEOUT',dri_restname_out) CALL SYSTEM & & ('cp '//TRIM(dri_restname_in)//' '//TRIM(dri_restname_out)) !- iret = NF90_OPEN (TRIM(sec_restname_out),NF90_NOWRITE,ncfid) iret = NF90_INQ_VARID (ncfid,'time',v_id) iret = NF90_GET_VAR (ncfid,v_id,r1d) time_sec = r1d(1) iret = NF90_INQ_VARID (ncfid,'time_steps',v_id) iret = NF90_GET_VAR (ncfid,v_id,time_step_sec) iret = NF90_CLOSE (ncfid) !- iret = NF90_OPEN (TRIM(dri_restname_out),NF90_WRITE,ncfid) iret = NF90_INQ_VARID (ncfid,'time',v_id) iret = NF90_GET_VAR (ncfid,v_id,r1d) time_dri = r1d(1) r1d(1) = time_sec iret = NF90_PUT_VAR (ncfid,v_id,r1d) iret = NF90_INQ_VARID (ncfid,'time_steps',v_id) iret = NF90_GET_VAR (ncfid,v_id,time_step_dri) iret = NF90_PUT_VAR (ncfid,v_id,time_step_sec) iret = NF90_INQ_VARID (ncfid,'julian',v_id) iret = NF90_GET_VAR (ncfid,v_id,r1d) julian = r1d(1) djulian = (time_step_sec-time_step_dri)*dtradia/one_day julian = julian & & +djulian-FLOAT(INT((julian+djulian)/one_year))*one_year r1d(1) = julian iret = NF90_PUT_VAR (ncfid,v_id,r1d) iret = NF90_CLOSE (ncfid) CALL getin_dump !--------------- END PROGRAM teststomate ! !=== !