source: tags/ORCHIDEE_1_9_5_2/ORCHIDEE_OL/forcesoil.f90 @ 560

PROGRAM forcesoil
  !---------------------------------------------------------------------
  ! This program reads a forcing file for STOMATE's soil carbon routine
  ! which was created with STOMATE.
  ! It then integrates STOMATE's soil carbon parameterizations
  ! for a given number of years using this forcing.
  !---------------------------------------------------------------------
  !- IPSL (2006)
  !-  This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
  USE netcdf
  !-
  USE defprec
  USE constantes
  USE constantes_veg
  USE constantes_co2
  USE stomate_constants
  USE ioipsl
  USE stomate_soilcarbon
  USE parallel
  !-
  IMPLICIT NONE
  !-
  CHARACTER(LEN=80) :: sto_restname_in,sto_restname_out
  INTEGER(i_std)                             :: iim,jjm

  INTEGER(i_std),PARAMETER                   :: llm = 1
  INTEGER(i_std)                             :: kjpindex

  INTEGER(i_std)                             :: itau_dep,itau_len
  CHARACTER(LEN=30)                         :: time_str
  REAL(r_std)                                :: dt_files
  REAL(r_std)                                :: date0
  INTEGER(i_std)                             :: rest_id_sto
  CHARACTER(LEN=20), SAVE                    :: thecalendar = 'noleap'
  !-
  CHARACTER(LEN=100) :: Cforcing_name
  INTEGER            :: Cforcing_id
  INTEGER            :: v_id
  REAL(r_std)                                :: dt_forcesoil
  INTEGER                                   :: nparan

  INTEGER(i_std),DIMENSION(:),ALLOCATABLE    :: indices
  INTEGER(i_std),DIMENSION(:),ALLOCATABLE    :: indices_g
  REAL(r_std),DIMENSION(:),ALLOCATABLE       :: x_indices_g
  REAL(r_std),DIMENSION(:,:),ALLOCATABLE     :: lon, lat
  REAL(r_std),DIMENSION(llm)                 :: lev


  INTEGER                                   :: i,m,iatt,iv,iyear

  CHARACTER(LEN=80)                         :: var_name
  CHARACTER(LEN=400)                        :: taboo_vars
  REAL(r_std),DIMENSION(1)                   :: xtmp
  INTEGER(i_std),PARAMETER                   :: nbvarmax=300
  INTEGER(i_std)                             :: nbvar
  CHARACTER(LEN=50),DIMENSION(nbvarmax)     :: varnames
  INTEGER(i_std)                             :: varnbdim
  INTEGER(i_std),PARAMETER                   :: varnbdim_max=20
  INTEGER,DIMENSION(varnbdim_max)           :: vardims
  LOGICAL                                   :: l1d
  REAL(r_std),DIMENSION(:,:),ALLOCATABLE     :: var_3d
  REAL(r_std)                                :: x_tmp
  ! string suffix indicating an index
  CHARACTER(LEN=10)  :: part_str
  !
  ! clay fraction
  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)  :: clay_g
  REAL(r_std),DIMENSION(:,:,:,:),ALLOCATABLE :: soilcarbon_input_g
  REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE   :: control_temp_g
  REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE   :: control_moist_g
  REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE   :: carbon_g

  REAL(r_std),ALLOCATABLE :: clay(:)
  REAL(r_std),ALLOCATABLE :: soilcarbon_input(:,:,:,:)
  REAL(r_std),ALLOCATABLE :: control_temp(:,:,:)
  REAL(r_std),ALLOCATABLE :: control_moist(:,:,:)
  REAL(r_std),ALLOCATABLE :: carbon(:,:,:)
  REAL(r_std),DIMENSION(:,:),ALLOCATABLE     :: resp_hetero_soil

  INTEGER(i_std)                             :: ier,iret

  LOGICAL :: debug

  CALL Init_para(.FALSE.) 
  CALL init_timer

!---------------------------------------------------------------------
!-
! set debug to have more information
!-
  !Config  Key  = DEBUG_INFO
  !Config  Desc = Flag for debug information
  !Config  Def  = n
  !Config  Help = This option allows to switch on the output of debug
  !Config         information without recompiling the code.
!-
  debug = .FALSE.
  CALL getin_p('DEBUG_INFO',debug)
  !-
  ! Stomate's restart files
  !-
  IF (is_root_prc) THEN
     sto_restname_in = 'stomate_start.nc'
     CALL getin ('STOMATE_RESTART_FILEIN',sto_restname_in)
     WRITE(numout,*) 'STOMATE INPUT RESTART_FILE: ',TRIM(sto_restname_in)
     sto_restname_out = 'stomate_rest_out.nc'
     CALL getin ('STOMATE_RESTART_FILEOUT',sto_restname_out)
     WRITE(numout,*) 'STOMATE OUTPUT RESTART_FILE: ',TRIM(sto_restname_out)
     !-
     ! We need to know iim_g, jjm.
     ! Get them from the restart files themselves.
     !-
     iret = NF90_OPEN (sto_restname_in, NF90_NOWRITE, rest_id_sto)
     iret = NF90_INQUIRE_DIMENSION (rest_id_sto,1,len=iim_g)
     iret = NF90_INQUIRE_DIMENSION (rest_id_sto,2,len=jjm_g)
     iret = NF90_INQ_VARID (rest_id_sto, "time", iv)
     iret = NF90_GET_ATT (rest_id_sto, iv, 'calendar',thecalendar)
     iret = NF90_CLOSE (rest_id_sto)
     i=INDEX(thecalendar,ACHAR(0))
     IF ( i > 0 ) THEN
        thecalendar(i:20)=' '
     ENDIF
     !-
     ! Allocate longitudes and latitudes
     !-
     ALLOCATE (lon(iim_g,jjm_g))
     ALLOCATE (lat(iim_g,jjm_g))
     lon(:,:) = 0.0
     lat(:,:) = 0.0
     lev(1)   = 0.0
     !-
     CALL restini &
          & (sto_restname_in, iim_g, jjm_g, lon, lat, llm, lev, &
          &  sto_restname_out, itau_dep, date0, dt_files, rest_id_sto)
  ENDIF
  CALL bcast(date0)
  CALL bcast(thecalendar)
  WRITE(numout,*) "calendar = ",thecalendar
  !-
  ! calendar
  !-
  CALL ioconf_calendar (thecalendar)
  CALL ioget_calendar  (one_year,one_day)
  CALL ioconf_startdate(date0)
  !
  IF (is_root_prc) THEN
     !-
     ! open FORCESOIL's forcing file to read some basic info
     !-
     Cforcing_name = 'NONE'
     CALL getin ('STOMATE_CFORCING_NAME',Cforcing_name)
     !-
     iret = NF90_OPEN (TRIM(Cforcing_name),NF90_NOWRITE,Cforcing_id)
     IF (iret /= NF90_NOERR) THEN
        CALL ipslerr (3,'forcesoil', &
             &        'Could not open file : ', &
             &          Cforcing_name,'(Do you have forget it ?)')
     ENDIF
     !-
     ier = NF90_GET_ATT (Cforcing_id,NF90_GLOBAL,'kjpindex',x_tmp)
     nbp_glo = NINT(x_tmp)
     ier = NF90_GET_ATT (Cforcing_id,NF90_GLOBAL,'nparan',x_tmp)
     nparan = NINT(x_tmp)
     !-
     ALLOCATE (indices_g(nbp_glo))
     ALLOCATE (clay_g(nbp_glo))
     !-
     ALLOCATE (x_indices_g(nbp_glo),stat=ier)
     ier = NF90_INQ_VARID (Cforcing_id,'index',v_id)
     ier = NF90_GET_VAR   (Cforcing_id,v_id,x_indices_g)
     indices_g(:) = NINT(x_indices_g(:))
     WRITE(numout,*) mpi_rank,"indices globaux : ",indices_g
     DEALLOCATE (x_indices_g)
     !-
     ier = NF90_INQ_VARID (Cforcing_id,'clay',v_id)
     ier = NF90_GET_VAR   (Cforcing_id,v_id,clay_g)
     !-
     ! time step of forcesoil
     !-
     dt_forcesoil = one_year / FLOAT(nparan)
     WRITE(numout,*) 'time step (d): ',dt_forcesoil
     !-
     ! read (and partially write) the restart file
     !-
     ! read and write the variables we do not need
     !-
     taboo_vars = '$lon$ $lat$ $lev$ $nav_lon$ $nav_lat$ $nav_lev$ $time$ $time_steps$ '// &
          &             '$day_counter$ $dt_days$ $date$ '// &
          &             '$carbon_01$ $carbon_02$ $carbon_03$ $carbon_04$ $carbon_05$'// &
          &             '$carbon_06$ $carbon_07$ $carbon_08$ $carbon_09$ $carbon_10$'// &
          &             '$carbon_11$ $carbon_12$ $carbon_13$'
     !-
     CALL ioget_vname(rest_id_sto, nbvar, varnames)
     !-
     ! read and write some special variables (1D or variables that we need)
     !-
     var_name = 'day_counter'
     CALL restget (rest_id_sto, var_name, 1, 1, 1, itau_dep, .TRUE., xtmp)
     CALL restput (rest_id_sto, var_name, 1, 1, 1, itau_dep, xtmp)
     !-
     var_name = 'dt_days'
     CALL restget (rest_id_sto, var_name, 1, 1, 1, itau_dep, .TRUE., xtmp)
     CALL restput (rest_id_sto, var_name, 1, 1, 1, itau_dep, xtmp)
     !-
     var_name = 'date'
     CALL restget (rest_id_sto, var_name, 1, 1, 1, itau_dep, .TRUE., xtmp)
     CALL restput (rest_id_sto, var_name, 1, 1, 1, itau_dep, xtmp)
     !-
     DO iv=1,nbvar
        !-- check if the variable is to be written here
        IF (INDEX(taboo_vars,'$'//TRIM(varnames(iv))//'$') == 0 ) THEN
           !---- get variable dimensions, especially 3rd dimension
           CALL ioget_vdim &
                &      (rest_id_sto, varnames(iv), varnbdim_max, varnbdim, vardims)
           l1d = ALL(vardims(1:varnbdim) == 1)
           !---- read it
           IF (l1d) THEN
              CALL restget &
                   &        (rest_id_sto, TRIM(varnames(iv)), 1, vardims(3), &
                   &         1, itau_dep, .TRUE., xtmp)
           ELSE
              ALLOCATE( var_3d(nbp_glo,vardims(3)), stat=ier)
              IF (ier /= 0) STOP 'ALLOCATION PROBLEM'
              !----
              CALL restget &
                   &        (rest_id_sto, TRIM(varnames(iv)), nbp_glo, vardims(3), &
                   &         1, itau_dep, .TRUE., var_3d, "gather", nbp_glo, indices_g)
           ENDIF
           !---- write it
           IF (l1d) THEN
              CALL restput &
                   &        (rest_id_sto, TRIM(varnames(iv)), 1, vardims(3), &
                   &         1, itau_dep, xtmp)
           ELSE
              CALL restput &
                   &        (rest_id_sto, TRIM(varnames(iv)), nbp_glo, vardims(3), &
                   &         1, itau_dep, var_3d, 'scatter',  nbp_glo, indices_g)
              !----
              DEALLOCATE(var_3d)
           ENDIF
        ENDIF
     ENDDO
     !-
     ! read soil carbon
     !-
     ALLOCATE(carbon_g(nbp_glo,ncarb,nvm))
     carbon_g(:,:,:) = val_exp
     DO m = 1, nvm
        WRITE (part_str, '(I2)') m
        IF (m<10) part_str(1:1)='0'
        var_name = 'carbon_'//part_str(1:LEN_TRIM(part_str))
        CALL restget &
             &    (rest_id_sto, var_name, nbp_glo, ncarb , 1, itau_dep, &
             &     .TRUE., carbon_g(:,:,m), 'gather', nbp_glo, indices_g)
        IF (ALL(carbon_g(:,:,m) == val_exp)) carbon_g(:,:,m) = zero
        !-- do not write this variable: it will be modified.
     ENDDO
     WRITE(numout,*) "date0 : ",date0, itau_dep
     !-
     ! Length of run
     !-
     WRITE(time_str,'(a)') '10000Y'
     CALL getin('TIME_LENGTH', time_str)
     write(numout,*) 'Number of years for carbon spinup : ',time_str
     ! transform into itau
     CALL tlen2itau(time_str, dt_forcesoil*one_day, date0, itau_len)
     write(numout,*) 'Number of time steps to do: ',itau_len
     !-
     ! read the rest of the forcing file and store forcing in an array.
     ! We read an average year.
     !-
     ALLOCATE(soilcarbon_input_g(nbp_glo,ncarb,nvm,nparan))
     ALLOCATE(control_temp_g(nbp_glo,nlevs,nparan))
     ALLOCATE(control_moist_g(nbp_glo,nlevs,nparan))
     !-
     ier = NF90_INQ_VARID (Cforcing_id,'soilcarbon_input',v_id)
     ier = NF90_GET_VAR   (Cforcing_id,v_id,soilcarbon_input_g)
     ier = NF90_INQ_VARID (Cforcing_id,   'control_moist',v_id)
     ier = NF90_GET_VAR   (Cforcing_id,v_id,control_moist_g)
     ier = NF90_INQ_VARID (Cforcing_id,    'control_temp',v_id)
     ier = NF90_GET_VAR   (Cforcing_id,v_id,control_temp_g)
     !-
     ier = NF90_CLOSE (Cforcing_id)
     !-
  ENDIF
  CALL bcast(nparan)
  CALL bcast(dt_forcesoil)
  CALL bcast(iim_g)
  CALL bcast(jjm_g)
  call bcast(nbp_glo)
  CALL bcast(itau_dep)
  CALL bcast(itau_len)
  !
  ! We must initialize data_para :
     !
     !
  CALL init_data_para(iim_g,jjm_g,nbp_glo,indices_g)

  kjpindex=nbp_loc
  jjm=jj_nb
  iim=iim_g
  IF (debug) WRITE(numout,*) "Local grid : ",kjpindex,iim,jjm

  !---
  !--- Create the index table
  !---
  !--- This job return a LOCAL kindex
  !---
  ALLOCATE (indices(kjpindex),stat=ier)
  CALL scatter(indices_g,indices)
  indices(1:kjpindex)=indices(1:kjpindex)-(jj_begin-1)*iim_g
  IF (debug) WRITE(numout,*) mpi_rank,"indices locaux = ",indices(1:kjpindex)
  !-
  !-
  ! there we go: time loop
  !-
  ALLOCATE(clay(kjpindex))
  ALLOCATE(soilcarbon_input(kjpindex,ncarb,nvm,nparan))
  ALLOCATE(control_temp(kjpindex,nlevs,nparan))
  ALLOCATE(control_moist(kjpindex,nlevs,nparan))
  ALLOCATE(carbon(kjpindex,ncarb,nvm))
  ALLOCATE(resp_hetero_soil(kjpindex,nvm))
  iatt = 0

  CALL Scatter(clay_g,clay)
  CALL Scatter(soilcarbon_input_g,soilcarbon_input)
  CALL Scatter(control_temp_g,control_temp)
  CALL Scatter(control_moist_g,control_moist)
  CALL Scatter(carbon_g,carbon)

  iyear=1
  DO i=1,itau_len
     iatt = iatt+1
     IF (iatt > nparan) THEN
        IF (debug) WRITE(numout,*) iyear
        iatt = 1
        iyear=iyear+1
     ENDIF
     CALL soilcarbon &
          &    (kjpindex, dt_forcesoil, clay, &
          &     soilcarbon_input(:,:,:,iatt), &
          &     control_temp(:,:,iatt), control_moist(:,:,iatt), &
          &     carbon, resp_hetero_soil)
  ENDDO
  WRITE(numout,*) "End of soilcarbon LOOP."
  CALL Gather(carbon,carbon_g)
  !-
  ! write new carbon into restart file
  !-
  IF (is_root_prc) THEN
     DO m=1,nvm
        WRITE (part_str, '(I2)') m
        IF (m<10) part_str(1:1)='0'
        var_name = 'carbon_'//part_str(1:LEN_TRIM(part_str))
        CALL restput &
             &    (rest_id_sto, var_name, nbp_glo, ncarb , 1, itau_dep, &
             &     carbon_g(:,:,m), 'scatter', nbp_glo, indices_g)
     ENDDO
     !-
     CALL getin_dump
     CALL restclo
  ENDIF
#ifdef CPP_PARA
  CALL MPI_FINALIZE(ier)
#endif
  WRITE(numout,*) "End of forcesoil."
  !--------------------
END PROGRAM forcesoil