! ! Daily update of leaf area index etc. ! !< $HeadURL$ !< $Date$ !< $Author$ !< $Revision$ !! IPSL (2006) !! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC ! ! MODULE slowproc ! modules used: USE defprec USE constantes USE constantes_veg USE constantes_co2 USE ioipsl USE sechiba_io USE interpol_help USE stomate USE stomate_constants USE grid USE parallel ! USE Write_Field_p IMPLICIT NONE PRIVATE PUBLIC slowproc_main,slowproc_clear ! To use OLD or NEW iterpollation schemes : INTERFACE slowproc_interlai MODULE PROCEDURE slowproc_interlai_OLD, slowproc_interlai_NEW END INTERFACE INTERFACE slowproc_interpol MODULE PROCEDURE slowproc_interpol_OLD, slowproc_interpol_NEW END INTERFACE INTERFACE slowproc_interpol_g MODULE PROCEDURE slowproc_interpol_OLD_g, slowproc_interpol_NEW_g END INTERFACE ! ! variables used inside slowproc module : declaration and initialisation ! LOGICAL, SAVE :: l_first_slowproc = .TRUE.!! Initialisation has to be done one time REAL(r_std), SAVE :: dt_slow !! time step of slow processes and STOMATE ! REAL(r_std), SAVE :: clayfraction_default = 0.2 INTEGER(i_std) , SAVE :: veget_update=0 !! update frequency in years for landuse INTEGER(i_std) , SAVE :: veget_year_orig=0 !! first year for landuse LOGICAL, SAVE :: land_use = .FALSE. ! Land Use LOGICAL, SAVE :: veget_reinit=.FALSE. !! To change LAND USE file in a run. ! LOGICAL, SAVE :: read_lai = .FALSE. ! Lai Map LOGICAL, SAVE :: old_lai = .FALSE. ! Old Lai Map interpolation LOGICAL, SAVE :: impveg = .FALSE. LOGICAL, SAVE :: impsoilt = .FALSE. LOGICAL, SAVE :: old_veget = .FALSE. ! Old veget Map interpolation ! REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: clayfraction REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: laimap ! REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: veget_nextyear !! next year fraction of vegetation type REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: frac_nobio_nextyear !! next year fraction of ice+lakes+cities+... REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: totfrac_nobio_nextyear !! next year total fraction of ice+lakes+cities+... CONTAINS SUBROUTINE slowproc_main (kjit, kjpij, kjpindex, dtradia, date0, & ldrestart_read, ldrestart_write, ldforcing_write, ldcarbon_write, & IndexLand, indexveg, lalo, neighbours, resolution, contfrac, soiltype, & t2m, t2m_min, temp_sol, stempdiag, & humrel, shumdiag, litterhumdiag, precip_rain, precip_snow, gpp, & deadleaf_cover, & assim_param, & lai, height, veget, frac_nobio, veget_max, totfrac_nobio, qsintmax, & rest_id, hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC, & co2_flux) ! interface description ! input scalar INTEGER(i_std), INTENT(in) :: kjit !! Time step number INTEGER(i_std), INTENT(in) :: kjpij !! Total size of the un-compressed grid INTEGER(i_std),INTENT (in) :: kjpindex !! Domain size REAL(r_std),INTENT (in) :: dtradia !! Time step REAL(r_std),INTENT (in) :: date0 !! Initial date LOGICAL, INTENT(in) :: ldrestart_read !! Logical for _restart_ file to read LOGICAL, INTENT(in) :: ldrestart_write !! Logical for _restart_ file to write LOGICAL, INTENT(in) :: ldforcing_write !! Logical for _forcing_ file to write LOGICAL, INTENT(in) :: ldcarbon_write !! Logical for _carbon_forcing_ file to write INTEGER(i_std),INTENT (in) :: rest_id,hist_id !! _Restart_ file and _history_ file identifier INTEGER(i_std),INTENT (in) :: hist2_id !! _history_ file 2 identifier INTEGER(i_std),INTENT (in) :: rest_id_stom !! STOMATE's _Restart_ file file identifier INTEGER(i_std),INTENT (in) :: hist_id_stom !! STOMATE's _history_ file file identifier INTEGER(i_std),INTENT(in) :: hist_id_stom_IPCC !! STOMATE's IPCC _history_ file file identifier ! input fields INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: IndexLand !! Indeces of the points on the map INTEGER(i_std),DIMENSION (kjpindex*nvm), INTENT (in) :: indexveg !! Indeces of the points on the 3D map REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geogr. coordinates (latitude,longitude) (degrees) INTEGER(i_std), DIMENSION (kjpindex,8), INTENT(in) :: neighbours !! neighoring grid points if land REAL(r_std), DIMENSION (kjpindex,2), INTENT(in) :: resolution !! size in x an y of the grid (m) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: contfrac !! Fraction of continent in the grid REAL(r_std), DIMENSION (kjpindex,nvm), INTENT (in) :: humrel !! Relative humidity ("moisture stress") REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: t2m !! 2 m air temperature (K) REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: t2m_min !! min. 2 m air temp. during forcing time step (K) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol !! Surface temperature REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (in) :: stempdiag !! Soil temperature REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (in) :: shumdiag !! Relative soil moisture REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: litterhumdiag !! Litter humidity REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: precip_rain !! Rain precipitation REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: precip_snow !! Snow precipitation REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: gpp !! GPP (gC/(m**2 of total ground)/time step) ! output scalar ! output fields REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(out) :: co2_flux !! CO2 flux in gC/m**2 of average ground/second ! modified scalar ! modified fields REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout):: lai !! Surface foliaire REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout):: height !! height (m) REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout):: veget !! Fraction of vegetation type REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (inout):: frac_nobio !! Fraction of ice, lakes, cities etc. in the mesh REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout):: veget_max !! Max fraction of vegetation type REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: totfrac_nobio !! Total fraction of ice+lakes+cities etc. in the mesh REAL(r_std), DIMENSION (kjpindex,nstm), INTENT(inout):: soiltype !! fraction of soil type REAL(r_std),DIMENSION (kjpindex,nvm,npco2), INTENT (inout):: assim_param!! min+max+opt temps, vcmax, vjmax for photosynthesis REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: deadleaf_cover !! fraction of soil covered by dead leaves REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout):: qsintmax !! Maximum water on vegetation for interception ! local declaration INTEGER(i_std) :: j, jv INTEGER(i_std), SAVE :: lcanop !! soil level used for LAI REAL(r_std) :: tmp_day(1) CHARACTER(LEN=80) :: var_name !! To store variables names for I/O ! REAL(r_std), DIMENSION(kjpindex,nvm) :: resp_maint !! autotrophic resp. (gC/(m**2 of total ground)/time step) REAL(r_std), DIMENSION(kjpindex,nvm) :: resp_hetero !! heterotrophic resp. (gC/(m**2 of total ground)/time step) REAL(r_std), DIMENSION(kjpindex,nvm) :: resp_growth !! growth resp. (gC/(m**2 of total ground)/time step) REAL(r_std), DIMENSION(kjpindex,nvm) :: npp !! Net Ecosystem Exchange (gC/(m**2 of total ground)/time step) ! INTEGER(i_std) , SAVE :: veget_year !! year for landuse LOGICAL :: land_use_updated ! LOGICAL, PARAMETER :: check = .FALSE. REAL(r_std), SAVE :: sec_old = zero ! ! do initialisation ! IF (l_first_slowproc) THEN ! ! 1.1 allocation, file definitions. Restart file read for Sechiba. Set flags. ! IF (long_print) WRITE (numout,*) ' l_first_slowproc : call slowproc_init ' CALL slowproc_init (kjit, ldrestart_read, dtradia, date0, kjpindex, IndexLand, lalo, neighbours, resolution, contfrac, & & rest_id, read_lai, lai, veget, frac_nobio, totfrac_nobio, soiltype, veget_max, height, lcanop,& & veget_update, veget_year) ! ! Time step in days for stomate dt_days = dt_slow / one_day ! resp_maint(:,:) = zero resp_hetero(:,:) = zero resp_growth(:,:) = zero ! ! 1.2 check time step ! IF ( dt_slow .LT. dtradia ) THEN WRITE(numout,*) 'slow_processes: time step smaller than forcing time step.' STOP 'slowproc_main' ENDIF ! IF ( control%stomate_watchout .OR. control%ok_stomate ) THEN ! ! 1.3 call STOMATE for initialization ! CALL stomate_main (kjit, kjpij, kjpindex, dtradia, dt_slow, & ldrestart_read, ldrestart_write, ldforcing_write, ldcarbon_write, & IndexLand, lalo, neighbours, resolution, contfrac, totfrac_nobio, clayfraction, & t2m, t2m_min, temp_sol, stempdiag, & humrel, shumdiag, litterhumdiag, precip_rain, precip_snow, gpp, & deadleaf_cover, & assim_param, & lai, height, veget, veget_max, qsintmax, & veget_nextyear, totfrac_nobio_nextyear, & hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC, & co2_flux,resp_maint,resp_hetero,resp_growth) ! ENDIF ! IF ( .NOT. control%ok_stomate ) THEN ! ! 1.4 initialize some variables ! STOMATE diagnoses some variables for SECHIBA: height, deadleaf_cover, etc. ! IF SECHIBA is not coupled to STOMATE, then we must set these values otherwise. ! CALL slowproc_derivvar (kjpindex, veget, lai, & qsintmax, deadleaf_cover, assim_param, height) ENDIF RETURN ENDIF ! !!$ ! Land USE for next year !!$ land_use_updated=.FALSE. !!$ IF ( (land_use) .AND. (veget_update .GT. 0) ) THEN !!$ ! if next iteration is divisibled by veget_update !!$ IF ( mod(kjit+1, veget_update) .le. min_sechiba) THEN !!$ ! !!$ veget_year=veget_year+veget_year_add !!$ WRITE(numout,*) 'We are updating veget for year =' , veget_year !!$ ! !!$ ! Save veget !!$ ! !!$ veget_lastyear(:,:) = veget_max(:,:) !!$ ! !!$ CALL slowproc_update(kjpindex, lalo, neighbours, resolution, contfrac, & !!$ & veget_max, frac_nobio, veget_year) !!$ !!$ land_use_updated=.TRUE. !!$ ENDIF !!$ ENDIF ! ! prepares restart file for the next simulation ! IF (ldrestart_write) THEN IF (long_print) WRITE (numout,*) ' we have to complete restart file with SLOWPROC variables ' tmp_day(1) = day_counter var_name= 'day_counter' IF (is_root_prc) CALL restput (rest_id, var_name, 1 , 1 , 1, kjit, tmp_day) var_name= 'veget' CALL restput_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, veget, 'scatter', nbp_glo, index_g) ! var_name= 'veget_max' CALL restput_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, veget_max, 'scatter', nbp_glo, index_g) ! var_name= 'lai' CALL restput_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, lai, 'scatter', nbp_glo, index_g) ! var_name= 'frac_nobio' CALL restput_p (rest_id, var_name, nbp_glo, nnobio, 1, kjit, frac_nobio, 'scatter', nbp_glo, index_g) ! var_name= 'soiltype_frac' CALL restput_p (rest_id, var_name, nbp_glo, nstm, 1, kjit, soiltype, 'scatter', nbp_glo, index_g) ! var_name= 'clay_frac' CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, clayfraction, 'scatter', nbp_glo, index_g) ! ! The height of the vegetation could in principle be recalculated at the beginning of the run. ! However, this is very tedious, as many special cases have to be taken into account. This variable ! is therefore saved in the restart file. var_name= 'height' CALL restput_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, height, 'scatter', nbp_glo, index_g) ! IF (read_lai) THEN var_name= 'laimap' CALL restput_p (rest_id, var_name, nbp_glo, nvm, 12, kjit, laimap) ENDIF ! IF (land_use) THEN tmp_day(1) = REAL(veget_year,r_std) var_name='veget_year' IF (is_root_prc) CALL restput (rest_id, var_name, 1 , 1 , 1, kjit, tmp_day) ENDIF ! ! call STOMATE to write RESTART files ! IF ( control%stomate_watchout .OR. control%ok_stomate ) THEN CALL stomate_main (kjit, kjpij, kjpindex, dtradia, dt_slow, & ldrestart_read, ldrestart_write, ldforcing_write, ldcarbon_write, & IndexLand, lalo, neighbours, resolution, contfrac, totfrac_nobio, clayfraction, & t2m, t2m_min, temp_sol, stempdiag, & humrel, shumdiag, litterhumdiag, precip_rain, precip_snow, gpp, & deadleaf_cover, & assim_param, & lai, height, veget, veget_max, qsintmax, & veget_nextyear, totfrac_nobio_nextyear, & hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC, & co2_flux,resp_maint,resp_hetero,resp_growth) ENDIF RETURN END IF ! ! update day counter day_counter = day_counter + dtradia IF (check) WRITE(*,*) "slowproc: day_counter 3",day_counter ! ! ! 1. Compute date ! ! Test each day and assert all slow processes (days and years) ! IF ( sec_old >= one_day - dtradia .AND. sec >= zero ) THEN ! ! reset counter day_counter = zero IF (check) WRITE(*,*) "slowproc: day_counter 2",day_counter ! ! Active slow processes do_slow = .TRUE. ! ! count days date = date + nint(dt_days) IF (check) WRITE(numout,*) "New date : ",date, 'year_length ',year_length,kjit ! ! is one year over? ! EndOfYear must be true once per year ! during a call of stomate_season. ! IF ( month == 1 .AND. day == 1 .AND. sec .LT. dtradia ) THEN EndOfYear = .TRUE. ELSE EndOfYear = .FALSE. ENDIF ELSE do_slow = .FALSE. EndOfYear = .FALSE. ENDIF sec_old = sec IF ( EndOfYear ) THEN WRITE(numout,*) 'slowproc: EndOfYear is activated.' ENDIF ! Land USE for next year land_use_updated=.FALSE. IF ( (land_use) .AND. (veget_update .GT. 0) ) THEN ! if next iteration is divisibled by veget_update IF ( EndOfYear ) THEN ! veget_year = veget_year + 1 ! IF ( MOD(veget_year - veget_year_orig, veget_update) == 0 ) THEN WRITE(numout,*) 'We are updating land use veget for year =' , veget_year ! CALL slowproc_update(kjpindex, lalo, neighbours, resolution, contfrac, & & veget_max, frac_nobio, veget_nextyear, frac_nobio_nextyear, veget_year) ! ! If some PFT has disapeared in new map WHERE(veget_nextyear(:,:).LT.veget(:,:)) veget(:,:) = veget_nextyear(:,:) ENDWHERE ! DO j = 1, nnobio totfrac_nobio_nextyear(:) = totfrac_nobio_nextyear(:) + frac_nobio_nextyear(:,j) ENDDO land_use_updated=.TRUE. ! ENDIF ! ENDIF ENDIF ! Land Use part IF (EndOfYear .AND. .NOT. land_use_updated) THEN lcchange=.FALSE. ENDIF IF ( control%stomate_watchout .OR. control%ok_stomate ) THEN ! ! 2. call STOMATE, either because we want to keep track of long-term variables or ! because STOMATE is activated ! CALL stomate_main (kjit, kjpij, kjpindex, dtradia, dt_slow, & ldrestart_read, ldrestart_write, ldforcing_write, ldcarbon_write, & IndexLand, lalo, neighbours, resolution, contfrac, totfrac_nobio, clayfraction, & t2m, t2m_min, temp_sol, stempdiag, & humrel, shumdiag, litterhumdiag, precip_rain, precip_snow, gpp, & deadleaf_cover, & assim_param, & lai, height, veget, veget_max, qsintmax, & veget_nextyear, totfrac_nobio_nextyear, & hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC, & co2_flux,resp_maint,resp_hetero,resp_growth) IF ( control%ok_stomate .AND. control%ok_sechiba ) THEN CALL histwrite(hist_id, 'maint_resp', kjit, resp_maint, kjpindex*nvm, indexveg) CALL histwrite(hist_id, 'hetero_resp', kjit, resp_hetero, kjpindex*nvm, indexveg) CALL histwrite(hist_id, 'growth_resp', kjit, resp_growth, kjpindex*nvm, indexveg) npp(:,1)=zero DO j = 2,nvm npp(:,j) = gpp(:,j) - resp_growth(:,j) - resp_maint(:,j) ENDDO CALL histwrite(hist_id, 'npp', kjit, npp, kjpindex*nvm, indexveg) ENDIF IF ( hist2_id > 0 ) THEN IF ( control%ok_stomate ) THEN CALL histwrite(hist2_id, 'maint_resp', kjit, resp_maint, kjpindex*nvm, indexveg) CALL histwrite(hist2_id, 'hetero_resp', kjit, resp_hetero, kjpindex*nvm, indexveg) CALL histwrite(hist2_id, 'growth_resp', kjit, resp_growth, kjpindex*nvm, indexveg) CALL histwrite(hist2_id, 'npp', kjit, npp, kjpindex*nvm, indexveg) ENDIF ENDIF ENDIF ! IF ( .NOT. control%ok_stomate ) THEN ! ! 2 STOMATE is not activated: we have to guess lai etc. ! ! ! 2.2 test if we have work to do ! IF ( do_slow .OR. land_use_updated ) THEN ! ! 2.2.1 do daily processes if necessary ! IF (long_print) WRITE (numout,*) 'slowproc_main : We update the daily variables' ! 2.2.2 updates lai CALL slowproc_lai (kjpindex, lcanop,stempdiag, & lalo,resolution,lai,month,day,read_lai,laimap) ! IF (land_use_updated) THEN veget_max(:,:)=veget_nextyear(:,:) frac_nobio(:,:)=frac_nobio_nextyear(:,:) ENDIF ! ! 2.2.3 updates veget CALL slowproc_veget (kjpindex, lai, frac_nobio, veget_max, veget) totfrac_nobio(:) = zero DO jv = 1, nnobio totfrac_nobio(:) = totfrac_nobio(:) + frac_nobio(:,jv) ENDDO ! 2.2.4 updates qsintmax and other derived variables CALL slowproc_derivvar (kjpindex, veget, lai, & qsintmax, deadleaf_cover, assim_param, height) ELSE ! IF (land_use_updated) THEN frac_nobio(:,:)=frac_nobio_nextyear(:,:) ENDIF ! END IF ! ! 2.3 some output fields ! co2_flux(:,:) = zero ENDIF IF (long_print) WRITE (numout,*) ' slowproc_main done ' END SUBROUTINE slowproc_main !! !! !! SUBROUTINE slowproc_init (kjit, ldrestart_read, dtradia, date0, kjpindex, IndexLand, lalo, neighbours, resolution, contfrac, & & rest_id, read_lai, lai, veget, frac_nobio, totfrac_nobio, soiltype, veget_max, height, lcanop,& & veget_update, veget_year) ! interface description ! input scalar INTEGER(i_std), INTENT (in) :: kjit !! Time step number LOGICAL, INTENT (in) :: ldrestart_read !! Logical for _restart_ file to read REAL(r_std),INTENT (in) :: dtradia !! Time step REAL(r_std), INTENT (in) :: date0 !! intial date INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size INTEGER(i_std), INTENT (in) :: rest_id !! _Restart_ file identifier ! input fields INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: IndexLand !! Indeces of the points on the map REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geogr. coordinates (latitude,longitude) (degrees) INTEGER(i_std), DIMENSION (kjpindex,8), INTENT(in) :: neighbours !! neighoring grid points if land REAL(r_std), DIMENSION (kjpindex,2), INTENT(in) :: resolution !! size in x an y of the grid (m) REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: contfrac !! Fraction of continent in the grid ! output scalar INTEGER(i_std), INTENT(out) :: lcanop !! soil level used for LAI INTEGER(i_std), INTENT(out) :: veget_update !! update frequency in timesteps for landuse INTEGER(i_std), INTENT(out) :: veget_year !! first year for landuse ! output fields REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: lai !! Surface foliere REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: veget !! Fraction of vegetation type REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (out) :: frac_nobio !! Fraction of ice,lakes,cities, ... REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: totfrac_nobio !! Total fraction of ice+lakes+cities+... REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: veget_max !! Max fraction of vegetation type REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: height !! Height of vegetation REAL(r_std), DIMENSION (kjpindex,nstm), INTENT(out) :: soiltype !! fraction of soil type ! local declaration REAL(r_std) :: tmp_day(1) REAL(r_std) :: tmp_veget_year(1) REAL(r_std) :: zcanop !! soil depth taken for canopy INTEGER(i_std) :: vtmp(1) REAL(r_std), DIMENSION(nbdl) :: zsoil !! soil depths at diagnostic levels INTEGER(i_std) :: j,l !! Index CHARACTER(LEN=80) :: var_name !! To store variables names for I/O INTEGER(i_std) :: ji, jv, ier LOGICAL, INTENT(out) :: read_lai REAL(r_std) :: frac_nobio1 !! temporary REAL(r_std) :: stempdiag_bid !! only needed for an initial LAI !! if there is no restart file REAL(r_std), DIMENSION(kjpindex,nbdl) :: stempdiag2_bid !! matrix to store stempdiag_bid CHARACTER(LEN=4) :: vegsoil_dist !! Flag to choose the soil/vegetation distribution ! CHARACTER(LEN=30), SAVE :: veget_str !! update frequency for landuse REAL(r_std),DIMENSION (nbp_glo,nnobio) :: frac_nobio_g !! Fraction of ice, lakes, cities etc. in the mesh (global) REAL(r_std),DIMENSION (nbp_glo,nvm) :: veget_max_g !! Fraction of vegetation type (globa) LOGICAL, PARAMETER :: check = .FALSE. ! ! 1 allocation ! ALLOCATE (clayfraction(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in clayfraction allocation. We stop. We need kjpindex words = ',kjpindex STOP 'slowproc_init' END IF clayfraction(:)=undef_sechiba ! veget_max(:,1) = un veget_max(:,2:nvm) = zero frac_nobio(:,:) = zero totfrac_nobio(:) = zero ! ier=-1 ALLOCATE(veget_nextyear(kjpindex, nvm), STAT=ier) IF (ier/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of veget_nextyear : ",ier STOP ENDIF veget_nextyear(:,1) = un veget_nextyear(:,2:nvm) = zero ! ier=-1 ALLOCATE(frac_nobio_nextyear(kjpindex, nnobio), STAT=ier) IF (ier/=0) THEN PRINT *,"ERROR IN ALLOCATION of frac_nobio_nextyear : ",ier STOP ENDIF frac_nobio_nextyear(:,:) = zero ! ier=-1 ALLOCATE(totfrac_nobio_nextyear(kjpindex), STAT=ier) IF (ier/=0) THEN PRINT *,"ERROR IN ALLOCATION of totfrac_nobio_nextyear : ",ier STOP ENDIF !MM must be corrected when nnobio > 1 !! totfrac_nobio_nextyear(:) = nnobio*un ! ! 2 read restart file ! var_name= 'day_counter' CALL ioconf_setatt('UNITS', 'd') CALL ioconf_setatt('LONG_NAME','Fraction of computed day') IF (is_root_prc) THEN CALL restget (rest_id, var_name, 1 , 1 , 1, kjit, .TRUE., tmp_day) IF (tmp_day(1) == val_exp) THEN day_counter = zero ELSE day_counter = tmp_day(1) ENDIF ENDIF CALL bcast(day_counter) ! get restart value if none were found in the restart file ! !Config Key = SECHIBA_DAY !Config Desc = Time within the day simulated !Config Def = 0.0 !Config Help = This is the time spent simulating the current day. This variable is !Config prognostic as it will trigger all the computations which are !Config only done once a day. ! CALL setvar_p (day_counter, val_exp, 'SECHIBA_DAY', zero) ! !Config Key = LAI_MAP !Config Desc = Read the LAI map !Config Def = n !Config Help = It is possible to read a 12 month LAI map which will !Config then be interpolated to daily values as needed. ! read_lai = .FALSE. CALL getin_p('LAI_MAP',read_lai) ! var_name= 'veget' CALL ioconf_setatt('UNITS', '-') CALL ioconf_setatt('LONG_NAME','Vegetation fraction') CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., veget, "gather", nbp_glo, index_g) ! var_name= 'veget_max' CALL ioconf_setatt('UNITS', '-') CALL ioconf_setatt('LONG_NAME','Maximum vegetation fraction') CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., veget_max, "gather", nbp_glo, index_g) ! frac_nobio(:,:) = val_exp var_name= 'frac_nobio' CALL ioconf_setatt('UNITS', '-') CALL ioconf_setatt('LONG_NAME','Special soil type fraction') CALL restget_p (rest_id, var_name, nbp_glo, nnobio, 1, kjit, .TRUE., frac_nobio, "gather", nbp_glo, index_g) ! !Config Key = LAND_USE !Config Desc = Read a land_use vegetation map !Config Def = n !Config Help = pft values are needed, max time axis is 293 ! land_use = .FALSE. veget_update=0 CALL getin_p('LAND_USE',land_use) IF (land_use) THEN ! !Config Key = VEGET_YEAR !Config Desc = Year of the land_use vegetation map to be read (0 == NO TIME AXIS) !Config If = LAND_USE !Config Def = 282 !Config Help = First year for landuse vegetation (2D map by pft). !Config Help If VEGET_YEAR == 0, this means there is no time axis. ! veget_year_orig=282 CALL getin_p('VEGET_YEAR', veget_year_orig) ! !Config Key = VEGET_REINIT !Config Desc = booleen to indicate that a new LAND USE file will be used. !Config If = LAND_USE !Config Def = n !Config Help = The parameter is used to bypass veget_year count !Config Help and reinitialize it with VEGET_YEAR parameter. !Config Help Then it is possible to change LAND USE file. ! veget_reinit = .FALSE. CALL getin_p('VEGET_REINIT', veget_reinit) ! ! var_name= 'veget_year' CALL ioconf_setatt('UNITS', '-') CALL ioconf_setatt('LONG_NAME','Last year get in Land Use file.') IF (is_root_prc) THEN CALL restget (rest_id, var_name, 1 , 1 , 1, kjit, .TRUE., tmp_veget_year) ! IF (tmp_veget_year(1) == val_exp) THEN veget_year=veget_year_orig ELSE IF (veget_reinit) THEN veget_year=veget_year_orig ELSE veget_year=INT(tmp_veget_year(1)) ENDIF ENDIF ENDIF CALL bcast(veget_year) ! !Config Key = VEGET_UPDATE !Config Desc = Update vegetation frequency !Config If = LAND_USE !Config Def = 0Y !Config Help = The veget datas will be update each this time step. ! veget_update=0 WRITE(veget_str,'(a)') '0Y' CALL getin_p('VEGET_UPDATE', veget_str) ! ! l=INDEX(TRIM(veget_str),'Y') READ(veget_str(1:(l-1)),"(I2.2)") veget_update WRITE(numout,*) "Update frequency for land use in years :",veget_update ! !Config Key = LAND_COVER_CHANGE !Config Desc = treat land use modifications !Config If = LAND_USE !Config Def = y !Config Help = With this variable, you can use a Land Use map !Config to simulate anthropic modifications such as !Config deforestation. ! lcchange = .TRUE. CALL getin_p('LAND_COVER_CHANGE', lcchange) IF ( veget_update == 0 .AND. lcchange ) THEN CALL ipslerr (2,'slowproc_init', & & 'You have asked for LAND_COVER_CHANGE activated with VEGET_UPDATE = 0Y.',& & 'We can''t use this land cover change model if veget is not updated.', & & 'We have disabled it.') lcchange=.FALSE. ENDIF ENDIF ! var_name= 'soiltype_frac' CALL ioconf_setatt('UNITS', '-') CALL ioconf_setatt('LONG_NAME','Fraction of each soil type') CALL restget_p (rest_id, var_name, nbp_glo, nstm, 1, kjit, .TRUE., soiltype, "gather", nbp_glo, index_g) ! var_name= 'clay_frac' CALL ioconf_setatt('UNITS', '-') CALL ioconf_setatt('LONG_NAME','Fraction of clay in each mesh') CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., clayfraction, "gather", nbp_glo, index_g) ! var_name= 'lai' CALL ioconf_setatt('UNITS', '-') CALL ioconf_setatt('LONG_NAME','Leaf area index') CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., lai, "gather", nbp_glo, index_g) ! ! The height of the vegetation could in principle be recalculated at the beginning of the run. ! However, this is very tedious, as many special cases have to be taken into account. This variable ! is therefore saved in the restart file. var_name= 'height' CALL ioconf_setatt('UNITS', 'm') CALL ioconf_setatt('LONG_NAME','Height of vegetation') CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., height, "gather", nbp_glo, index_g) ! IF (read_lai)THEN ! ALLOCATE (laimap(kjpindex,nvm,12),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in laimap allocation. We stop. We need kjpindex*nvm*12 words = ',kjpindex*nvm*12 STOP 'slowproc_init' END IF laimap(:,:,:) = val_exp ! var_name= 'laimap' CALL ioconf_setatt('UNITS', '-') CALL ioconf_setatt('LONG_NAME','Leaf area index read') CALL restget_p (rest_id, var_name, nbp_glo, nvm, 12, kjit, .TRUE., laimap) ! ELSE ! ALLOCATE (laimap(1,1,1)) ENDIF ! !Config Key = SECHIBA_ZCANOP !Config Desc = Soil level (m) used for canopy development (if STOMATE disactivated) !Config Def = 0.5 !Config Help = The temperature at this soil depth is used to determine the LAI when !Config STOMATE is not activated. ! zcanop = 0.5_r_std CALL setvar_p (zcanop, val_exp, 'SECHIBA_ZCANOP', 0.5_r_std) ! !Config Key = HYDROL_SOIL_DEPTH !Config Desc = Total depth of soil reservoir !Config Def = 2. ! dpu_cste=2. CALL getin_p ("HYDROL_SOIL_DEPTH", dpu_cste) dpu(:)=dpu_cste ! !Config Key = HYDROL_HUMCSTE !Config Desc = Root profile !Config Def = 5., .8, .8, 1., .8, .8, 1., 1., .8, 4., 4., 4., 4. !Config Help = Default values were defined for 2 meters soil depth. !Config For 4 meters soil depth, you may use those ones : !Config 5., .4, .4, 1., .8, .8, 1., 1., .8, 4., 1., 4., 1. ! humcste(:)= & & (/5., .8, .8, 1., .8, .8, 1., 1., .8, 4., 4., 4., 4./) CALL getin_p ("HYDROL_HUMCSTE", humcste) !MM, T. d'O. : before in constantes_soil : ! diaglev = & ! & (/ 0.001, 0.004, 0.01, 0.018, 0.045, & ! & 0.092, 0.187, 0.375, 0.750, 1.5, & ! & 2.0 /) ! Place here because it is used in sechiba and stomate (then in teststomate) ! to be in sechiba when teststomate will have disapeared. !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) = dpu_cste/(2**(nbdl-1) -1) * ( ( 2**(jv-1) -1) + ( 2**(jv) -1) ) / 2.0 ENDDO diaglev(nbdl) = dpu_cste ! depth at center of the levels zsoil(1) = diaglev(1) / 2. DO l = 2, nbdl zsoil(l) = ( diaglev(l) + diaglev(l-1) ) / 2. ENDDO ! index of this level vtmp = MINLOC ( ABS ( zcanop - zsoil(:) ) ) lcanop = vtmp(1) ! ! Interception reservoir coefficient ! !Config Key = 'SECHIBA_QSINT' !Config Desc = Interception reservoir coefficient !Config Def = 0.1 !Config Help = Transforms leaf area index into size of interception reservoir !Config for slowproc_derivvar or stomate qsintcst = 0.1 CALL getin_p('SECHIBA_QSINT', qsintcst) WRITE(numout, *)' SECHIBA_QSINT, qsintcst = ', qsintcst ! ! Time step of STOMATE and LAI update ! !Config Key = DT_SLOW !Config Desc = Time step of STOMATE and other slow processes !Config Def = one_day !Config Help = Time step (s) of regular update of vegetation !Config cover, LAI etc. This is also the time step !Config of STOMATE. dt_slow = one_day CALL getin_p('DT_SLOW', dt_slow) ! !Config Key = SLOWPROC_LAI_TEMPDIAG !Config Desc = Temperature used for the initial guess of LAI !Config Def = 280. !Config Help = If there is no LAI in the restart file, we may need !Config a temperature that is used to guess the initial LAI. ! stempdiag_bid = 280. CALL getin_p('SLOWPROC_LAI_TEMPDIAG',stempdiag_bid) ! ! ! get restart value if none were found in the restart file ! !Config Key = AGRICULTURE !Config Desc = agriculture allowed? !Config Def = y !Config Help = With this variable, you can determine !Config whether agriculture is allowed ! agriculture = .TRUE. CALL getin_p('AGRICULTURE', agriculture) IF ( .NOT. agriculture .AND. land_use ) THEN CALL ipslerr (2,'slowproc_init', & & 'Problem with agriculture desactivated and Land Use activated.',& & 'Are you sure ?', & & '(check your parameters).') ENDIF ! !Config Key = IMPOSE_VEG !Config Desc = Should the vegetation be prescribed !Config Def = n !Config Help = This flag allows the user to impose a vegetation distribution !Config and its characterisitcs. It is espacially interesting for 0D !Config simulations. On the globe it does not make too much sense as !Config it imposes the same vegetation everywhere ! impveg = .FALSE. CALL getin_p('IMPOSE_VEG', impveg) ! IF ( impveg ) THEN ! ! We are on a point and thus we can read the information from the run.def ! !Config Key = SECHIBA_VEG !Config Desc = Vegetation distribution within the mesh (0-dim mode) !Config If = IMPOSE_VEG !Config Def = 0.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.8, 0.0, 0.0, 0.0 !Config Help = The fraction of vegetation is read from the restart file. If !Config it is not found there we will use the values provided here. ! CALL setvar_p (veget, val_exp, 'SECHIBA_VEG', veget_ori_fixed_test_1) ! !Config Key = SECHIBA_VEGMAX !Config Desc = Maximum vegetation distribution within the mesh (0-dim mode) !Config If = IMPOSE_VEG !Config Def = 0.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.8, 0.0, 0.0, 0.0 !Config Help = The fraction of vegetation is read from the restart file. If !Config it is not found there we will use the values provided here. ! CALL setvar_p (veget_max, val_exp, 'SECHIBA_VEGMAX', veget_ori_fixed_test_1) ! !Config Key = SECHIBA_FRAC_NOBIO !Config Desc = Fraction of other surface types within the mesh (0-dim mode) !Config If = IMPOSE_VEG !Config Def = 0.0 !Config Help = The fraction of ice, lakes, etc. is read from the restart file. If !Config it is not found there we will use the values provided here. !Config For the moment, there is only ice. ! ! laisser ca tant qu'il n'y a que la glace. Pb avec setvar? frac_nobio1 = frac_nobio(1,1) CALL setvar_p (frac_nobio1, val_exp, 'SECHIBA_FRAC_NOBIO', frac_nobio_fixed_test_1) frac_nobio(:,:) = frac_nobio1 ! CALL setvar (frac_nobio, val_exp, 'SECHIBA_FRAC_NOBIO', frac_nobio_fixed_test_1) ! !Config Key = SECHIBA_LAI !Config Desc = LAI for all vegetation types (0-dim mode) !Config Def = 0., 8., 8., 4., 4.5, 4.5, 4., 4.5, 4., 2., 2., 2., 2. !Config If = IMPOSE_VEG !Config Help = The maximum LAI used in the 0dim mode. The values should be found !Config in the restart file. The new values of LAI will be computed anyway !Config at the end of the current day. The need for this variable is caused !Config by the fact that the model may stop during a day and thus we have not !Config yet been through the routines which compute the new surface conditions. ! CALL setvar_p (lai, val_exp, 'SECHIBA_LAI', llaimax) ! !Config Key = IMPOSE_SOILT !Config Desc = Should the soil typ be prescribed !Config Def = n !Config If = IMPOSE_VEG !Config Help = This flag allows the user to impose a soil type distribution. !Config It is espacially interesting for 0D !Config simulations. On the globe it does not make too much sense as !Config it imposes the same soil everywhere ! impsoilt = .FALSE. CALL getin_p('IMPOSE_SOILT', impsoilt) IF (impsoilt) THEN !Config Key = SOIL_FRACTIONS !Config Desc = Fraction of the 3 soil types (0-dim mode) !Config Def = 0.28, 0.52, 0.20 !Config If = IMPOSE_VEG !Config If = IMPOSE_SOILT !Config Help = Determines the fraction for the 3 soil types !Config in the mesh in the following order : sand loam and clay. ! CALL setvar_p (soiltype, val_exp, 'SOIL_FRACTIONS', soiltype_default) !Config Key = CLAY_FRACTION !Config Desc = Fraction of the clay fraction (0-dim mode) !Config Def = 0.2 !Config If = IMPOSE_VEG !Config If = IMPOSE_SOILT !Config Help = Determines the fraction of clay in the grid box. ! CALL setvar_p (clayfraction, val_exp, 'CLAY_FRACTION', clayfraction_default) ELSE IF ( MINVAL(soiltype) .EQ. MAXVAL(soiltype) .AND. MAXVAL(soiltype) .EQ. val_exp .OR.& & MINVAL(clayfraction) .EQ. MAXVAL(clayfraction) .AND. MAXVAL(clayfraction) .EQ. val_exp) THEN CALL slowproc_soilt(kjpindex, lalo, neighbours, resolution, contfrac, soiltype, clayfraction) ENDIF ENDIF ! !Config Key = SLOWPROC_HEIGHT !Config Desc = Height for all vegetation types (m) !Config Def = 0., 30., 30., 20., 20., 20., 15., 15., 15., .5, .6, 1.0, 1.0 !Config If = IMPOSE_VEG !Config Help = The height used in the 0dim mode. The values should be found !Config in the restart file. The new values of height will be computed anyway !Config at the end of the current day. The need for this variable is caused !Config by the fact that the model may stop during a day and thus we have not !Config yet been through the routines which compute the new surface conditions. ! CALL setvar_p (height, val_exp, 'SLOWPROC_HEIGHT', height_presc) ELSE ! ! We are in the full 2-D case thus we need to do the interpolation if the potential vegetation ! is not available ! IF ( ALL( veget_max(:,:) .EQ. val_exp ) .OR. ALL( frac_nobio(:,:) .EQ. val_exp ) ) THEN IF ( .NOT. land_use ) THEN !Config Key = SLOWPROC_VEGET_OLD_INTERPOL !Config Desc = Flag to use old "interpolation" of vegetation map. !Config If = NOT IMPOSE_VEG and NOT LAND_USE !Config Def = FALSE !Config Help = If you want to recover the old (ie orchidee_1_2 branch) !Config "interpolation" of vegetation map. ! old_veget = .FALSE. CALL getin_p('SLOWPROC_VEGET_OLD_INTERPOL',old_veget) ! The interpolation of vegetation has changed. IF (is_root_prc) THEN IF ( .NOT. old_veget ) THEN ! NEW slowproc interpol : CALL slowproc_interpol_g(nbp_glo, lalo_g, neighbours_g, resolution_g, contfrac_g, veget_max_g, frac_nobio_g) ELSE ! OLD slowproc interpol : CALL slowproc_interpol_g(nbp_glo, lalo_g, neighbours_g, resolution_g, veget_max_g, frac_nobio_g) ENDIF ENDIF CALL scatter(veget_max_g,veget_max) CALL scatter(frac_nobio_g, frac_nobio) ! IF ( control%ok_dgvm ) THEN ! ! If we are dealing with dynamic vegetation then all ! natural PFTs should be set to veget_max = 0 ! And in case no agriculture is desired, agriculture PFTS should be ! set to 0 as well IF (agriculture) THEN veget_max(:,2:nvm-2) = zero DO ji = 1, kjpindex veget_max(ji,1) = un - SUM(veget_max(ji,nvm-1:nvm)) & - SUM(frac_nobio(ji,:)) ENDDO ELSE veget_max(:,:) = zero DO ji = 1, kjpindex veget_max(ji,1) = un - SUM(frac_nobio(ji,:)) ENDDO ENDIF ! ENDIF ELSE WRITE(numout,*) 'We initialize land use veget for year =' , veget_year ! If restart doesn''t contain veget, then it is the first computation CALL slowproc_update(kjpindex, lalo, neighbours, resolution, contfrac, & & veget_nextyear, frac_nobio_nextyear, veget_max, frac_nobio, veget_year, init=.TRUE.) ! IF ( control%ok_dgvm ) THEN ! ! If we are dealing with dynamic vegetation then all ! natural PFTs should be set to veget_max = 0 ! And in case no agriculture is desired, agriculture PFTS should be ! set to 0 as well IF (agriculture) THEN veget_max(:,2:nvm-2) = zero DO ji = 1, kjpindex veget_max(ji,1) = un - SUM(veget_max(ji,nvm-1:nvm)) & - SUM(frac_nobio(ji,:)) ENDDO ELSE veget_max(:,:) = zero DO ji = 1, kjpindex veget_max(ji,1) = un - SUM(frac_nobio(ji,:)) ENDDO ENDIF ! ENDIF ! ENDIF ! ELSE ! WITH restarts for vegetation and DGVM and NO AGRICULTURE IF ( control%ok_dgvm ) THEN ! ! If we are dealing with dynamic vegetation then all ! natural PFTs should be set to veget_max = 0 ! And in case no agriculture is desired, agriculture PFTS should be ! set to 0 as well ! IF (.NOT. agriculture) THEN DO ji = 1, kjpindex veget_max(ji,1) = veget_max(ji,1) + SUM(veget_max(ji,nvm-1:nvm)) ENDDO veget_max(ji,nvm-1:nvm) = zero ENDIF ! ENDIF ! ENDIF ! totfrac_nobio(:) = zero DO j = 1, nnobio totfrac_nobio(:) = totfrac_nobio(:) + frac_nobio(:,j) ENDDO ! ! IF (read_lai) THEN !Config Key = SLOWPROC_LAI_OLD_INTERPOL !Config Desc = Flag to use old "interpolation" of LAI !Config If = LAI_MAP !Config Def = FALSE !Config Help = If you want to recover the old (ie orchidee_1_2 branch) !Config "interpolation" of LAI map. ! old_lai = .FALSE. CALL getin_p('SLOWPROC_LAI_OLD_INTERPOL',old_lai) ! ! In case the LAI map was not found in the restart then we need to recompute it ! IF ( ALL( laimap(:,:,:) .EQ. val_exp) ) THEN ! The interpolation of LAI has changed. IF ( .NOT. old_lai ) THEN ! NEW slowproc interlai : CALL slowproc_interlai (kjpindex, lalo, resolution, neighbours, contfrac, laimap) ELSE ! OLD slowproc interlai : CALL slowproc_interlai(kjpindex, lalo, resolution, laimap) ENDIF ! ! Compute the current LAI just to be sure ! stempdiag2_bid(1:kjpindex,1:nbdl) = stempdiag_bid CALL slowproc_lai (kjpindex, lcanop, stempdiag2_bid, & lalo,resolution,lai,month,day,read_lai,laimap) ! ! Make sure that we redo the computation when we will be back in slowproc_main day_counter = dt_slow - dtradia ! ENDIF ! ENDIF ! IF ( MINVAL(lai) .EQ. MAXVAL(lai) .AND. MAXVAL(lai) .EQ. val_exp) THEN ! ! Get a first guess at the LAI ! IF ( read_lai ) THEN IF ( ALL( laimap(:,:,:) .EQ. val_exp) ) THEN ! The interpolation of LAI has changed. IF ( .NOT. old_lai ) THEN ! NEW slowproc interlai : CALL slowproc_interlai (kjpindex, lalo, resolution, neighbours, contfrac, laimap) ELSE ! OLD slowproc interlai : CALL slowproc_interlai(kjpindex, lalo, resolution, laimap) ENDIF ENDIF ENDIF ! stempdiag2_bid(1:kjpindex,1:nbdl) = stempdiag_bid CALL slowproc_lai (kjpindex, lcanop, stempdiag2_bid, & lalo,resolution,lai,month,day,read_lai,laimap) ! Make sure that we redo the computation when we will be back in slowproc_main day_counter = dt_slow - dtradia ENDIF IF ( MINVAL(veget) .EQ. MAXVAL(veget) .AND. MAXVAL(veget) .EQ. val_exp) THEN ! Impose height DO jv = 1, nvm height(:,jv) = height_presc(jv) ENDDO ! Have a first guess at the vegetation fraction CALL slowproc_veget (kjpindex, lai, frac_nobio, veget_max, veget) ENDIF IF ( MINVAL(soiltype) .EQ. MAXVAL(soiltype) .AND. MAXVAL(soiltype) .EQ. val_exp .OR.& & MINVAL(clayfraction) .EQ. MAXVAL(clayfraction) .AND. MAXVAL(clayfraction) .EQ. val_exp) THEN CALL slowproc_soilt(kjpindex, lalo, neighbours, resolution, contfrac, soiltype, clayfraction) ENDIF ENDIF ! ! Select the preferences for the distribution of the 13 PFTs on the 3 soil types. ! vegsoil_dist='EQUI' ! SELECT CASE(vegsoil_dist) ! ! A reasonable choice ! CASE('MAXR') pref_soil_veg(:,1) = (/ 1, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2 /) pref_soil_veg(:,2) = (/ 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3 /) pref_soil_veg(:,3) = (/ 3, 1, 1, 1, 1, 1 ,1 ,1 ,1 ,1 ,1 ,1, 1 /) ! ! Current default : equidistribution. ! CASE('EQUI') ! pref_soil_veg(:,:) = zero ! CASE DEFAULT ! WRITE(numout,*) 'The vegetation soil type preference you have chosen does not exist' WRITE(numout,*) 'You chose :', vegsoil_dist STOP 'slowproc_init' ! END SELECT ! ! ! calculate total fraction of the mesh that is covered by particular surface types: ice, lakes, ... ! totfrac_nobio(:) = zero DO jv = 1, nnobio totfrac_nobio(:) = totfrac_nobio(:) + frac_nobio(:,jv) ENDDO l_first_slowproc = .FALSE. IF (long_print) WRITE (numout,*) ' slowproc_init done ' END SUBROUTINE slowproc_init !! !! Clear Memory !! SUBROUTINE slowproc_clear l_first_slowproc = .TRUE. IF (ALLOCATED (clayfraction)) DEALLOCATE (clayfraction) IF (ALLOCATED (laimap)) DEALLOCATE (laimap) ! ! IF (ALLOCATED (veget_nextyear)) DEALLOCATE (veget_nextyear) IF (ALLOCATED (frac_nobio_nextyear)) DEALLOCATE (frac_nobio_nextyear) IF (ALLOCATED (totfrac_nobio_nextyear)) DEALLOCATE (totfrac_nobio_nextyear) ! CALL stomate_clear ! END SUBROUTINE slowproc_clear !! !! Derive some variables !! SUBROUTINE slowproc_derivvar (kjpindex, veget, lai, & qsintmax, deadleaf_cover, assim_param, height) ! interface description ! input scalar INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size ! input fields REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Fraction of vegetation type REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: lai !! Surface foliere ! output scalar ! output fields REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: qsintmax !! Maximum water on vegetation for interception REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: deadleaf_cover!! fraction of soil covered by dead leaves REAL(r_std), DIMENSION (kjpindex,nvm,npco2), INTENT (out) :: assim_param !! min+max+opt temps & vmax for photosynthesis REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: height !! height (m) ! ! local declaration ! INTEGER(i_std) :: ji, jv ! ! 1 Assimilation parameters ! DO jv = 1, nvm assim_param(:,jv,itmin) = co2_tmin_fix(jv) + tp_00 assim_param(:,jv,itopt) = co2_topt_fix(jv) + tp_00 assim_param(:,jv,itmax) = co2_tmax_fix(jv) + tp_00 assim_param(:,jv,ivcmax) = vcmax_fix(jv) assim_param(:,jv,ivjmax) = vjmax_fix(jv) ENDDO ! ! 2 fraction of soil covered by dead leaves ! deadleaf_cover(:) = zero ! ! 3 height ! DO jv = 1, nvm height(:,jv) = height_presc(jv) ENDDO ! ! 4 qsintmax ! qsintmax(:,:) = qsintcst * veget(:,:) * lai(:,:) ! Ajout Nathalie - Juillet 2006 qsintmax(:,1) = zero END SUBROUTINE slowproc_derivvar !! !! !! SUBROUTINE slowproc_mean (npts, n_dim2, dt_tot, dt, ldmean, field_in, field_mean) ! Accumulates field_in over a period of dt_tot. ! Has to be called at every time step (dt). ! Mean value is calculated if ldmean=.TRUE. ! field_mean must be initialized outside of this routine! ! ! 0 declarations ! ! 0.1 input ! Domain size INTEGER(i_std), INTENT(in) :: npts ! 2nd dimension (1 or npft) INTEGER(i_std), INTENT(in) :: n_dim2 ! Time step of STOMATE (days) REAL(r_std), INTENT(in) :: dt_tot ! Time step in days REAL(r_std), INTENT(in) :: dt ! Calculate mean ? LOGICAL, INTENT(in) :: ldmean ! Daily field REAL(r_std), DIMENSION(npts,n_dim2), INTENT(in) :: field_in ! 0.2 modified field ! Annual field REAL(r_std), DIMENSION(npts,n_dim2), INTENT(inout) :: field_mean ! ========================================================================= ! ! 1 accumulation ! field_mean(:,:) = field_mean(:,:) + field_in(:,:) * dt ! ! 2 mean fields ! IF (ldmean) THEN field_mean(:,:) = field_mean(:,:) / dt_tot ENDIF END SUBROUTINE slowproc_mean !! !! !! SUBROUTINE slowproc_long (npts, n_dim2, dt, tau, field_inst, field_long) ! Calculates a temporally smoothed field (field_long) from instantaneous ! input fields. ! Time constant tau determines the strength of the smoothing. ! For tau -> infty, field_long becomes the true mean value of field_inst (but ! the spinup becomes infinietly long, too). ! field_long must be initialized outside of this routine! ! ! 0 declarations ! ! 0.1 input ! Domain size INTEGER(i_std), INTENT(in) :: npts ! 2nd dimension (1 or npft) INTEGER(i_std), INTENT(in) :: n_dim2 ! Time step REAL(r_std), INTENT(in) :: dt ! Integration time constant (has to have same unit as dt!) REAL(r_std), INTENT(in) :: tau ! Instantaneous field REAL(r_std), DIMENSION(npts,n_dim2), INTENT(in) :: field_inst ! 0.2 modified field ! Long-term field REAL(r_std), DIMENSION(npts,n_dim2), INTENT(inout) :: field_long ! ========================================================================= ! ! 1 test coherence ! IF ( ( tau .LT. dt ) .OR. ( dt .LE. zero ) .OR. ( tau .LE. zero ) ) THEN WRITE(numout,*) 'slowproc_long: Problem with time steps' WRITE(numout,*) 'dt=',dt WRITE(numout,*) 'tau=',tau ENDIF ! ! 2 integration ! field_long(:,:) = ( field_inst(:,:)*dt + field_long(:,:)*(tau-dt) ) / tau END SUBROUTINE slowproc_long !! !! !! SUBROUTINE slowproc_veget (kjpindex, lai, frac_nobio, veget_max, veget) ! ! 0. Declarations ! ! 0.1 Input INTEGER(i_std), INTENT(in) :: kjpindex REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: lai ! 0.2 Modified REAL(r_std), DIMENSION(kjpindex,nnobio), INTENT(inout) :: frac_nobio REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(inout) :: veget_max ! 0.3 Output REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(out) :: veget ! 0.4 Local REAL(r_std), DIMENSION(kjpindex) :: fracsum INTEGER(i_std) :: nbad INTEGER(i_std) :: ji, jv ! Test Nathalie REAL(r_std) :: SUMveg !!$ REAL(r_std) :: trans_veg ! ! 1. Sum up veget_max and frac_nobio and test if sum is equal to 1 ! ! ! 1.1 Sum up ! fracsum(:) = zero DO jv = 1, nnobio DO ji = 1, kjpindex fracsum(ji) = fracsum(ji) + frac_nobio(ji,jv) ENDDO ENDDO DO jv = 1, nvm DO ji = 1, kjpindex fracsum(ji) = fracsum(ji) + veget_max(ji,jv) ENDDO ENDDO ! ! 1.2 stop if there is a severe problem, that is an error above 0.01% ! The rest will be scaled ! nbad = COUNT( ABS(fracsum(:)-un) .GT. 0.0001 ) IF ( nbad .GT. 0 ) THEN WRITE(numout,*) 'Problem with total surface areas.' DO ji = 1, kjpindex IF ( ABS(fracsum(ji)-un) .GT. 0.0001 ) THEN WRITE(numout,*) 'Point :', ji WRITE(numout,*) ' frac_nobio :', frac_nobio(ji,:) WRITE(numout,*) ' Veget_max :', veget_max(ji,:) WRITE(numout,*) ' Fracsum :', fracsum(ji), EPSILON(un) WRITE(numout,*) ' The error is :', un - ( SUM(frac_nobio(ji,:)) + SUM(veget_max(ji,:)) ) STOP 'slowproc_veget' ENDIF ENDDO ENDIF ! ! 1.3 correction at places where the problem is surely precision-related ! nbad = COUNT( ABS(fracsum(:)-un) .GT. EPSILON(un) ) ! IF ( nbad .GT. 0 ) THEN ! IF ( long_print ) THEN WRITE(numout,*) 'WARNING! scaling frac_nobio and veget_max at', nbad, ' points' ENDIF ! DO jv = 1, nnobio DO ji = 1, kjpindex IF ( ABS(fracsum(ji)-un) .GT. EPSILON(un) ) THEN frac_nobio(ji,jv) = frac_nobio(ji,jv)/fracsum(ji) ENDIF ENDDO ENDDO ! DO jv = 1, nvm DO ji = 1, kjpindex IF ( ABS(fracsum(ji)-un) .GT. EPSILON(un) ) THEN veget_max(ji,jv) = veget_max(ji,jv)/fracsum(ji) ENDIF ENDDO ENDDO ! ENDIF ! ! 2. Set veget equal to veget_max ! DO jv = 1, nvm DO ji = 1, kjpindex veget(ji,jv) = veget_max(ji,jv) ENDDO ENDDO ! ! 3. if lai of a vegetation type (jv > 1) is small, increase soil part ! ! Nathalie - nouveau calcul de veget !!$ DO jv = 1,nvm !!$ DO ji = 1, kjpindex !!$ !!$ IF ((jv .GT. 1) .AND. (lai(ji,jv) .LT. 0.5)) THEN !!$ trans_veg = 2.*(0.5 - lai(ji,jv)) * veget(ji,jv) !!$ ! We limit the smallest fraction to 0.5% !!$ IF ( veget(ji,jv) - trans_veg .GT. min_vegfrac ) THEN !!$ veget(ji,1) = veget(ji,1) + trans_veg !!$ veget(ji,jv) = veget(ji,jv) - trans_veg !!$ ELSE !!$ veget(ji,1) = veget(ji,1) + veget(ji,jv) !!$ veget(ji,jv) = zero !!$ ENDIF !!$ ENDIF !!$ !!$ ENDDO !!$ ENDDO ! Ajout Nouveau calcul (stomate-like) DO ji = 1, kjpindex SUMveg = zero veget(ji,1) = veget_max(ji,1) DO jv = 2, nvm veget(ji,jv) = veget_max(ji,jv) * ( un - exp( - lai(ji,jv) * ext_coef(jv) ) ) veget(ji,1) = veget(ji,1) + (veget_max(ji,jv) - veget(ji,jv)) SUMveg = SUMveg + veget(ji,jv) ENDDO SUMveg = SUMveg + veget(ji,1) + SUM(frac_nobio(ji,:)) IF (SUMveg .LT. 0.99999) THEN WRITE(numout,*)' ATTENTION, en ji, SUMveg LT 1: ', ji, SUMveg WRITE(numout,*)' frac_nobio = ',SUM(frac_nobio(ji,:)) WRITE(numout,*) ' ',veget(ji,:) ENDIF ENDDO ! ! 4. Sum up surface fractions and test if the sum is equal to 1 ! ! ! 4.1 Sum up ! fracsum(:) = zero DO jv = 1, nnobio DO ji = 1, kjpindex fracsum(ji) = fracsum(ji) + frac_nobio(ji,jv) ENDDO ENDDO DO jv = 1, nvm DO ji = 1, kjpindex fracsum(ji) = fracsum(ji) + veget_max(ji,jv) ENDDO ENDDO ! ! 4.2 stop if there is a severe problem ! nbad = COUNT( ABS(fracsum(:)-un) .GT. (REAL(nvm+nnobio,r_std)*EPSILON(un)) ) ! IF ( nbad .GT. 0 ) THEN WRITE(numout,*) 'Problem with veget or frac_nobio.' DO ji = 1, kjpindex IF ( ABS(fracsum(ji)-un) .GT. (10.*EPSILON(un)) ) THEN WRITE(numout,*) 'Point :', ji WRITE(numout,*) ' frac_nobio :', frac_nobio(ji,:) WRITE(numout,*) ' Veget :', veget(ji,:) WRITE(numout,*) ' The error is :', un - (SUM(frac_nobio(ji,:)) + SUM(veget(ji,:))) STOP 'slowproc_veget' ENDIF ENDDO ENDIF ! ! 4.3 correction at places where the problem is surely precision-related ! nbad = COUNT( ABS(fracsum(:)-un) .GT. EPSILON(un) ) ! IF ( nbad .GT. 0 ) THEN ! IF ( long_print ) THEN WRITE(numout,*) 'slowproc_veget: scaling veget at', nbad, ' points' ENDIF ! DO jv = 1, nvm DO ji = 1, kjpindex IF ( ABS(fracsum(ji)-un) .GT. EPSILON(un) ) THEN veget(ji,jv) = veget(ji,jv) / fracsum(ji) ENDIF ENDDO ENDDO ! DO jv = 1, nnobio DO ji = 1, kjpindex IF ( ABS(fracsum(ji)-un) .GT. EPSILON(un) ) THEN frac_nobio(ji,jv) = frac_nobio(ji,jv) / fracsum(ji) ENDIF ENDDO ENDDO ! ENDIF END SUBROUTINE slowproc_veget !! !! !! SUBROUTINE slowproc_lai (kjpindex,lcanop,stempdiag,lalo,resolution,lai,mm,dd,read_lai,laimap) ! ! 0. Declarations ! ! 0.1 Input INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size INTEGER(i_std), INTENT(in) :: lcanop !! soil level used for LAI INTEGER(i_std), INTENT(in) :: mm, dd REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (in) :: stempdiag !! Soil temperature REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geogr. coordinates (latitude,longitude) (degrees) REAL(r_std), DIMENSION (kjpindex,2), INTENT(in) :: resolution !! size in x an y of the grid (m) REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: laimap !! LAI lue LOGICAL, INTENT(in) :: read_lai ! 0.2 Output REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(out) :: lai !! LAI ! 0.3 Local INTEGER(i_std) :: ji,jv ! Test Nathalie. On impose LAI PFT 1 a 0 ! On boucle sur 2,nvm au lieu de 1,nvm lai(: ,1) = zero DO jv = 2,nvm !!$ DO jv = 1,nvm SELECT CASE (type_of_lai(jv)) CASE ("mean ") ! ! 1. do the interpolation between laimax and laimin ! IF ( .NOT. read_lai ) THEN lai(:,jv) = undemi * (llaimax(jv) + llaimin(jv)) ELSE DO ji = 1,kjpindex lai(ji,jv) = MAXVAL(laimap(ji,jv,:)) ENDDO ENDIF ! CASE ("inter") ! ! 2. do the interpolation between laimax and laimin ! IF ( .NOT. read_lai ) THEN DO ji = 1,kjpindex lai(ji,jv) = llaimin(jv) + tempfunc(stempdiag(ji,lcanop)) * (llaimax(jv) - llaimin(jv)) ENDDO ELSE IF (mm .EQ. 1 ) THEN IF (dd .LE. 15) THEN lai(:,jv) = laimap(:,jv,12)*(1-(dd+15)/30.) + laimap(:,jv,1)*((dd+15)/30.) ELSE lai(:,jv) = laimap(:,jv,1)*(1-(dd-15)/30.) + laimap(:,jv,2)*((dd-15)/30.) ENDIF ELSE IF (mm .EQ. 12) THEN IF (dd .LE. 15) THEN lai(:,jv) = laimap(:,jv,11)*(1-(dd+15)/30.) + laimap(:,jv,12)*((dd+15)/30.) ELSE lai(:,jv) = laimap(:,jv,12)*(1-(dd-15)/30.) + laimap(:,jv,1)*((dd-15)/30.) ENDIF ELSE IF (dd .LE. 15) THEN lai(:,jv) = laimap(:,jv,mm-1)*(1-(dd+15)/30.) + laimap(:,jv,mm)*((dd+15)/30.) ELSE lai(:,jv) = laimap(:,jv,mm)*(1-(dd-15)/30.) + laimap(:,jv,mm+1)*((dd-15)/30.) ENDIF ENDIF ENDIF ! CASE default ! ! 3. Problem ! WRITE (numout,*) 'This kind of lai choice is not possible. '// & ' We stop with type_of_lai ',jv,' = ', type_of_lai(jv) STOP 'slowproc_lai' END SELECT ENDDO END SUBROUTINE slowproc_lai !! !! Interpolate the LAI map to the grid of the model !MM TAG 1.6 model ! !! SUBROUTINE slowproc_interlai_OLD(nbpt, lalo, resolution, laimap) ! ! ! ! 0.1 INPUT ! INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longitudes (beware of the order !) REAL(r_std), INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid-box in X and Y ! ! 0.2 OUTPUT ! REAL(r_std), INTENT(out) :: laimap(nbpt,nvm,12) ! lai read variable and re-dimensioned ! ! 0.3 LOCAL ! REAL(r_std), PARAMETER :: min_sechiba=1.E-8 ! ! CHARACTER(LEN=80) :: filename INTEGER(i_std) :: iml, jml, ijml, i, j, ik, lml, tml, fid, ib, jb,ip, jp, vid, ai,iki,jkj REAL(r_std) :: lev(1), date, dt, coslat, pi INTEGER(i_std) :: itau(1) REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: mask_lu REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_lu, lon_lu, mask REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_ful, lon_ful REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:) :: laimaporig REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:,:) :: laimap_lu REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lon_up, lon_low, lat_up, lat_low INTEGER, DIMENSION(nbpt) :: n_origlai INTEGER, DIMENSION(nbpt) :: n_found REAL(r_std), DIMENSION(nbpt,nvm) :: frac_origlai CHARACTER(LEN=80) :: meter REAL(r_std) :: prog, sumf LOGICAL :: found INTEGER :: idi,jdi, ilast, jlast, jj, ii, jv, inear, iprog REAL(r_std) :: domaine_lon_min, domaine_lon_max, domaine_lat_min, domaine_lat_max ! pi = 4. * ATAN(1.) ! !Config Key = LAI_FILE !Config Desc = Name of file from which the vegetation map is to be read !Config If = !LAI_MAP !Config Def = ../surfmap/lai2D.nc !Config Help = The name of the file to be opened to read the LAI !Config map is to be given here. Usualy SECHIBA runs with a 5kmx5km !Config map which is derived from a Nicolas VIOVY one. ! filename = 'lai2D.nc' CALL getin_p('LAI_FILE',filename) ! IF (is_root_prc) CALL flininfo(filename, iml, jml, lml, tml, fid) CALL bcast(iml) CALL bcast(jml) CALL bcast(lml) CALL bcast(tml) ! ! ALLOCATE(lon_lu(iml)) ALLOCATE(lat_lu(jml)) ALLOCATE(laimap_lu(iml,jml,nvm,tml)) ALLOCATE(mask_lu(iml,jml)) ! WRITE(numout,*) 'slowproc_interlai : Reading the LAI file' ! IF (is_root_prc) THEN CALL flinget(fid, 'longitude', iml, 0, 0, 0, 1, 1, lon_lu) CALL flinget(fid, 'latitude', jml, 0, 0, 0, 1, 1, lat_lu) CALL flinget(fid, 'LAI', iml, jml, nvm, tml, 1, 12, laimap_lu) CALL flinget(fid, 'mask', iml, jml, 0, 0, 0, 1, mask_lu) ! CALL flinclo(fid) ENDIF CALL bcast(lon_lu) CALL bcast(lat_lu) CALL bcast(laimap_lu) CALL bcast(mask_lu) ! WRITE(numout,*) 'slowproc_interlai : ', lon_lu(1), lon_lu(iml),lat_lu(1), lat_lu(jml) ! ! ijml=iml*jml ALLOCATE(lon_ful(ijml)) ALLOCATE(lat_ful(ijml)) ALLOCATE(laimaporig(ijml,nvm,tml)) ALLOCATE(mask(ijml)) DO i=1,iml DO j=1,jml iki=(j-1)*iml+i lon_ful(iki)=lon_lu(i) lat_ful(iki)=lat_lu(j) laimaporig(iki,:,:)=laimap_lu(i,j,:,:) mask(iki)=mask_lu(i,j) ENDDO ENDDO ! WHERE ( laimaporig(:,:,:) .LT. 0 ) laimaporig(:,:,:) = zero ENDWHERE ! ! ALLOCATE(lon_up(nbpt)) ALLOCATE(lon_low(nbpt)) ALLOCATE(lat_up(nbpt)) ALLOCATE(lat_low(nbpt)) ! DO ib =1, nbpt ! ! We find the 4 limits of the grid-box. As we transform the resolution of the model ! into longitudes and latitudes we do not have the problem of periodicity. ! coslat is a help variable here ! ! coslat = MAX(COS(lalo(ib,1) * pi/180. ), 0.001 )*pi/180. * R_Earth ! lon_up(ib) = lalo(ib,2) + resolution(ib,1)/(2.0*coslat) lon_low(ib) = lalo(ib,2) - resolution(ib,1)/(2.0*coslat) ! coslat = pi/180. * R_Earth ! lat_up(ib) = lalo(ib,1) + resolution(ib,2)/(2.0*coslat) lat_low(ib) = lalo(ib,1) - resolution(ib,2)/(2.0*coslat) ! ! ! ENDDO lon_up = NINT( lon_up * 1E6 ) * 1E-6 lon_low = NINT( lon_low * 1E6 ) * 1E-6 lat_up = NINT( lat_up * 1E6 ) * 1E-6 lat_low = NINT( lat_low * 1E6 ) * 1E-6 ! ! Get the limits of the integration domaine so that we can speed up the calculations ! domaine_lon_min = MINVAL(lon_low) domaine_lon_max = MAXVAL(lon_up) domaine_lat_min = MINVAL(lat_low) domaine_lat_max = MAXVAL(lat_up) ! !!$ WRITE(*,*) 'DOMAINE lon :', domaine_lon_min, domaine_lon_max !!$ WRITE(*,*) 'DOMAINE lat :', domaine_lat_min, domaine_lat_max ! ! Ensure that the fine grid covers the whole domain WHERE ( lon_ful(:) .LT. domaine_lon_min ) lon_ful(:) = lon_ful(:) + 360. ENDWHERE ! WHERE ( lon_ful(:) .GT. domaine_lon_max ) lon_ful(:) = lon_ful(:) - 360. ENDWHERE lon_ful = NINT( lon_ful * 1E6 ) * 1E-6 lat_ful = NINT( lat_ful * 1E6 ) * 1E-6 ! WRITE(numout,*) 'Interpolating the lai map :' WRITE(numout,'(2a40)')'0%--------------------------------------', & & '------------------------------------100%' ! ilast = 1 n_origlai(:) = 0 laimap(:,:,:) = zero ! DO ip=1,ijml ! ! Give a progress meter ! ! prog = ip/float(ijml)*79. ! IF ( ABS(prog - NINT(prog)) .LT. 1/float(ijml)*79. ) THEN ! meter(NINT(prog)+1:NINT(prog)+1) = 'x' ! WRITE(numout, advance="no", FMT='(a)') ACHAR(13) ! WRITE(numout, advance="no", FMT='(a80)') meter ! ENDIF iprog = NINT(float(ip)/float(ijml)*79.) - NINT(float(ip-1)/float(ijml)*79.) IF ( iprog .NE. 0 ) THEN WRITE(numout,'(a1,$)') 'y' ENDIF ! ! Only start looking for its place in the smaler grid if we are within the domaine ! That should speed up things ! ! IF ( ( lon_ful(ip) .GE. domaine_lon_min ) .AND. & ( lon_ful(ip) .LE. domaine_lon_max ) .AND. & ( lat_ful(ip) .GE. domaine_lat_min ) .AND. & ( lat_ful(ip) .LE. domaine_lat_max ) ) THEN ! ! look for point on GCM grid which this point on fine grid belongs to. ! First look at the point on the model grid where we arrived just before. There is ! a good chace that neighbouring points on the fine grid fall into the same model ! grid box. ! IF ( ( lon_ful(ip) .GE. lon_low(ilast) ) .AND. & ( lon_ful(ip) .LT. lon_up(ilast) ) .AND. & ( lat_ful(ip) .GE. lat_low(ilast) ) .AND. & ( lat_ful(ip) .LT. lat_up(ilast) ) ) THEN ! ! We were lucky ! IF (mask(ip) .GT. 0) THEN n_origlai(ilast) = n_origlai(ilast) + 1 DO i=1,nvm DO j=1,12 laimap(ilast,i,j) = laimap(ilast,i,j) + laimaporig(ip,i,j) ENDDO ENDDO ENDIF ! ELSE ! ! Otherwise, look everywhere. ! Begin close to last grid point. ! found = .FALSE. idi = 1 ! DO WHILE ( (idi .LT. nbpt) .AND. ( .NOT. found ) ) ! ! forward and backward ! DO ii = 1,2 ! IF ( ii .EQ. 1 ) THEN ib = ilast - idi ELSE ib = ilast + idi ENDIF ! IF ( ( ib .GE. 1 ) .AND. ( ib .LE. nbpt ) ) THEN IF ( ( lon_ful(ip) .GE. lon_low(ib) ) .AND. & ( lon_ful(ip) .LT. lon_up(ib) ) .AND. & ( lat_ful(ip) .GE. lat_low(ib) ) .AND. & ( lat_ful(ip) .LT. lat_up(ib) ) ) THEN ! IF (mask(ip) .gt. 0) THEN DO i=1,nvm DO j=1,12 laimap(ib,i,j) = laimap(ib,i,j) + laimaporig(ip,i,j) ENDDO ENDDO n_origlai(ib) = n_origlai(ib) + 1 ENDIF ilast = ib found = .TRUE. ! ENDIF ENDIF ! ENDDO ! idi = idi + 1 ! ENDDO ! ENDIF ! lucky/not lucky ! ENDIF ! in the domain ENDDO ! determine fraction of LAI points in each box of the coarse grid DO ip=1,nbpt IF ( n_origlai(ip) .GT. 0 ) THEN DO jv =1,nvm laimap(ip,jv,:) = laimap(ip,jv,:)/REAL(n_origlai(ip),r_std) ENDDO ELSE ! ! Now we need to handle some exceptions ! IF ( lalo(ip,1) .LT. -56.0) THEN ! Antartica DO jv =1,nvm laimap(ip,jv,:) = zero ENDDO ! ELSE IF ( lalo(ip,1) .GT. 70.0) THEN ! Artica DO jv =1,nvm laimap(ip,jv,:) = zero ENDDO ! ELSE IF ( lalo(ip,1) .GT. 55.0 .AND. lalo(ip,2) .GT. -65.0 .AND. lalo(ip,2) .LT. -20.0) THEN ! Greenland DO jv =1,nvm laimap(ip,jv,:) = zero ENDDO ! ELSE ! WRITE(numout,*) 'PROBLEM, no point in the lai map found for this grid box' WRITE(numout,*) 'Longitude range : ', lon_low(ip), lon_up(ip) WRITE(numout,*) 'Latitude range : ', lat_low(ip), lat_up(ip) ! WRITE(numout,*) 'Looking for nearest point on the lai map file' CALL slowproc_nearest (ijml, lon_ful, lat_ful, & lalo(ip,2), lalo(ip,1), inear) WRITE(numout,*) 'Coordinates of the nearest point, ',inear,' :', & lon_ful(inear),lat_ful(inear) ! DO jv =1,nvm laimap(ip,jv,:) = laimaporig(inear,jv,:) ENDDO ENDIF ENDIF ENDDO ! WRITE(numout,*) 'slowproc_interlai : Interpolation Done' ! ! ! DEALLOCATE(lon_up) DEALLOCATE(lon_low) DEALLOCATE(lat_up) DEALLOCATE(lat_low) DEALLOCATE(lat_ful) DEALLOCATE(lon_ful) DEALLOCATE(lat_lu) DEALLOCATE(lon_lu) DEALLOCATE(laimap_lu) DEALLOCATE(laimaporig) DEALLOCATE(mask_lu) DEALLOCATE(mask) ! RETURN ! END SUBROUTINE slowproc_interlai_OLD !! !! Interpolate the LAI map to the grid of the model !! SUBROUTINE slowproc_interlai_NEW(nbpt, lalo, resolution, neighbours, contfrac, laimap) ! ! ! ! 0.1 INPUT ! INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longitudes (beware of the order !) REAL(r_std), INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid-box in X and Y ! INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) ! Vector of neighbours for each grid point (1=N, 2=E, 3=S, 4=W) REAL(r_std), INTENT(in) :: contfrac(nbpt) ! Fraction of land in each grid box. ! ! 0.2 OUTPUT ! REAL(r_std), INTENT(out) :: laimap(nbpt,nvm,12) ! lai read variable and re-dimensioned ! ! 0.3 LOCAL ! ! CHARACTER(LEN=80) :: filename INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, ip, jp, it, jj, jv REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_lu, lon_lu REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lat, lon REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: sub_area INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:) :: sub_index REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:,:) :: laimap_lu INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: mask ! REAL(r_std) :: totarea INTEGER(i_std) :: idi, nbvmax CHARACTER(LEN=30) :: callsign ! LOGICAL :: ok_interpol ! optionnal return of aggregate_2d ! INTEGER :: ALLOC_ERR ! !Config Key = LAI_FILE !Config Desc = Name of file from which the vegetation map is to be read !Config If = LAI_MAP !Config Def = ../surfmap/lai2D.nc !Config Help = The name of the file to be opened to read the LAI !Config map is to be given here. Usualy SECHIBA runs with a 5kmx5km !Config map which is derived from a Nicolas VIOVY one. ! filename = 'lai2D.nc' CALL getin_p('LAI_FILE',filename) ! IF (is_root_prc) CALL flininfo(filename, iml, jml, lml, tml, fid) CALL bcast(iml) CALL bcast(jml) CALL bcast(lml) CALL bcast(tml) ! ! ALLOC_ERR=-1 ALLOCATE(lon_lu(iml), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of lon_lu : ",ALLOC_ERR STOP ENDIF ALLOC_ERR=-1 ALLOCATE(lat_lu(jml), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of lat_lu : ",ALLOC_ERR STOP ENDIF ALLOC_ERR=-1 ALLOCATE(laimap_lu(iml,jml,nvm,tml), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of laimap_lu : ",ALLOC_ERR STOP ENDIF ! ! IF (is_root_prc) THEN CALL flinget(fid, 'longitude', iml, 0, 0, 0, 1, 1, lon_lu) CALL flinget(fid, 'latitude', jml, 0, 0, 0, 1, 1, lat_lu) CALL flinget(fid, 'LAI', iml, jml, nvm, tml, 1, 12, laimap_lu) ! WHERE (laimap_lu(:,:,:,:) < zero ) laimap_lu(:,:,:,:) = zero ENDWHERE ! CALL flinclo(fid) ENDIF CALL bcast(lon_lu) CALL bcast(lat_lu) CALL bcast(laimap_lu) ! ALLOC_ERR=-1 ALLOCATE(lon(iml,jml), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of lon : ",ALLOC_ERR STOP ENDIF ALLOC_ERR=-1 ALLOCATE(lat(iml,jml), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of lat : ",ALLOC_ERR STOP ENDIF ! DO ip=1,iml lat(ip,:) = lat_lu(:) ENDDO DO jp=1,jml lon(:,jp) = lon_lu(:) ENDDO ! ALLOC_ERR=-1 ALLOCATE(mask(iml,jml), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of mask : ",ALLOC_ERR STOP ENDIF ! ! Consider all points a priori ! !!$ mask(:,:) = 1 mask(:,:) = 0 ! DO ip=1,iml DO jp=1,jml ! ! Exclude the points where there is never a LAI value. It is probably ! an ocean point. ! !!$ IF (ABS(SUM(laimap_lu(ip,jp,:,:))) < EPSILON(laimap_lu) ) THEN !!$ mask(ip,jp) = 0 !!$ ENDIF ! IF ( ANY(laimap_lu(ip,jp,:,:) < 20.) ) THEN mask(ip,jp) = 1 ENDIF ! ENDDO ENDDO ! nbvmax = 20 ! callsign = 'LAI map' ! ok_interpol = .FALSE. DO WHILE ( .NOT. ok_interpol ) WRITE(numout,*) "Projection arrays for ",callsign," : " WRITE(numout,*) "nbvmax = ",nbvmax ! ALLOC_ERR=-1 ALLOCATE(sub_index(nbpt, nbvmax, 2), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of sub_index : ",ALLOC_ERR STOP ENDIF sub_index(:,:,:)=0 ALLOC_ERR=-1 ALLOCATE(sub_area(nbpt, nbvmax), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of sub_area : ",ALLOC_ERR STOP ENDIF sub_area(:,:)=zero ! CALL aggregate_p(nbpt, lalo, neighbours, resolution, contfrac, & & iml, jml, lon, lat, mask, callsign, & & nbvmax, sub_index, sub_area, ok_interpol) ! IF ( .NOT. ok_interpol ) THEN DEALLOCATE(sub_area) DEALLOCATE(sub_index) ENDIF ! nbvmax = nbvmax * 2 ENDDO ! laimap(:,:,:) = zero ! DO ib=1,nbpt idi = COUNT(sub_area(ib,:) > zero) IF ( idi > 0 ) THEN totarea = zero DO jj=1,idi ip = sub_index(ib,jj,1) jp = sub_index(ib,jj,2) DO jv=1,nvm DO it=1,12 laimap(ib,jv,it) = laimap(ib,jv,it) + laimap_lu(ip,jp,jv,it)*sub_area(ib,jj) ENDDO ENDDO totarea = totarea + sub_area(ib,jj) ENDDO ! ! Normalize ! laimap(ib,:,:) = laimap(ib,:,:)/totarea !!$ IF ( ANY( laimap(ib,:,:) > 10 ) ) THEN !!$ WRITE(numout,*) "For point ",ib !!$ WRITE(numout,*) lalo(ib,1), lalo(ib,2) !!$ WRITE(numout,*) "For ib=",ib," we have LAI ",laimap(ib,:,1:idi) !!$ WRITE(numout,*) "Detail of projection :" !!$ WRITE(numout,*) sub_area(ib,1:idi) !!$ WRITE(numout,*) "Total for projection :" ,totarea !!$ ENDIF ! ELSE WRITE(numout,*) 'On point ', ib, ' no points where found for interpolating the LAI map.' WRITE(numout,*) 'Location : ', lalo(ib,2), lalo(ib,1) DO jv=1,nvm laimap(ib,jv,:) = (llaimax(jv)+llaimin(jv))/deux ENDDO WRITE(numout,*) 'Solved by putting the average LAI for the PFT all year long' ENDIF ENDDO ! WRITE(numout,*) 'slowproc_interlai : Interpolation Done' ! ! ! DEALLOCATE(lat_lu) DEALLOCATE(lon_lu) DEALLOCATE(lon) DEALLOCATE(lat) DEALLOCATE(laimap_lu) DEALLOCATE(mask) DEALLOCATE(sub_area) DEALLOCATE(sub_index) ! RETURN ! END SUBROUTINE slowproc_interlai_NEW !! !! Interpolate a vegetation map (by pft) !! !MM modif TAG 1.4 : ! SUBROUTINE slowproc_update(nbpt, lalo, resolution, vegetmap, frac_nobiomap) !MM modif TAG 1.9.2 : ! SUBROUTINE slowproc_update(nbpt, lalo, neighbours, resolution, contfrac, vegetmax, frac_nobio, veget_year, init) SUBROUTINE slowproc_update(nbpt, lalo, neighbours, resolution, contfrac, & & veget_last, frac_nobio_last, & & veget_next, frac_nobio_next, veget_year, init) ! ! ! ! 0.1 INPUT ! INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs ! to be interpolated REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: lalo ! Vector of latitude and longitudes (beware of the order !) !MM modif TAG 1.4 : add grid variables for aggregate INTEGER(i_std), DIMENSION(nbpt,8), INTENT(in) :: neighbours ! Vector of neighbours for each grid point ! (1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW) REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: resolution ! The size in km of each grid-box in X and Y REAL(r_std), DIMENSION(nbpt), INTENT(in) :: contfrac ! Fraction of continent in the grid ! REAL(r_std), DIMENSION(nbpt,nvm), INTENT(in) :: veget_last ! old max vegetfrac REAL(r_std), DIMENSION(nbpt,nnobio), INTENT(in) :: frac_nobio_last ! old fraction of the mesh which is ! covered by ice, lakes, ... ! INTEGER(i_std), INTENT(in) :: veget_year ! first year for landuse (0 == NO TIME AXIS) LOGICAL, OPTIONAL, INTENT(in) :: init ! initialisation : in case of dgvm, it forces update of all PFTs ! ! 0.2 OUTPUT ! REAL(r_std), DIMENSION(nbpt,nvm), INTENT(out) :: veget_next ! new max vegetfrac REAL(r_std), DIMENSION(nbpt,nnobio), INTENT(out) :: frac_nobio_next ! new fraction of the mesh which is ! covered by ice, lakes, ... ! ! 0.3 LOCAL ! ! CHARACTER(LEN=80) :: filename INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, ip, jp, inobio, jv INTEGER(i_std) :: nb_coord,nb_var, nb_gat,nb_dim REAL(r_std) :: date, dt INTEGER(i_std), ALLOCATABLE, DIMENSION(:) :: itau REAL(r_std), DIMENSION(1) :: time_counter REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_lu, lon_lu INTEGER,DIMENSION(flio_max_var_dims) :: l_d_w, i_d_w LOGICAL :: exv, l_ex ! !MM modif TAG 1.4 : suppression of all time axis reading and interpolation, replaced by each year 2D reading. ! REAL(r_std), INTENT(inout) :: vegetmap(nbpt,nvm,293) ! vegetfrac read variable and re-dimensioned ! REAL(r_std), INTENT(inout) :: frac_nobiomap(nbpt,nnobio,293) ! Fraction of the mesh which is covered by ice, lakes, ... REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:,:) :: vegmap ! last coord is time with only one value = 1 REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:) :: vegmap_1 ! last coord is time with only one value = 1 (IF VEGET_YEAR == 0 , NO TIME AXIS) REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lat_ful, lon_ful REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: sub_area INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:) :: sub_index INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: mask ! REAL(r_std) :: sumf, err, norm REAL(r_std) :: totarea INTEGER(i_std) :: idi, nbvmax CHARACTER(LEN=30) :: callsign ! LOGICAL :: ok_interpol ! optionnal return of aggregate_2d ! ! for DGVM case : REAL(r_std) :: sum_veg ! sum of vegets REAL(r_std) :: sum_nobio ! sum of nobios REAL(r_std) :: sumvAnthro_old, sumvAnthro ! last an new sum of antrhopic vegets REAL(r_std) :: rapport ! (S-B) / (S-A) LOGICAL :: partial_update ! if TRUE, partialy update PFT (only anthropic ones) ! e.g. in case of DGVM and not init (optional parameter) ! LOGICAL, PARAMETER :: debug = .FALSE. ! INTEGER :: ALLOC_ERR ! !Config Key = VEGETATION_FILE !Config Desc = Name of file from which the vegetation map is to be read !Config If = LAND_USE !Config Def = pft_new.nc !Config Help = The name of the file to be opened to read a vegetation !Config map (in pft) is to be given here. ! filename = 'pft_new.nc' CALL getin_p('VEGETATION_FILE',filename) ! IF (is_root_prc) THEN IF (debug) THEN WRITE(numout,*) "Entering slowproc_update. Debug mode." WRITE (*,'(/," --> fliodmpf")') CALL fliodmpf (TRIM(filename)) WRITE (*,'(/," --> flioopfd")') ENDIF CALL flioopfd (TRIM(filename),fid,nb_dim=nb_coord,nb_var=nb_var,nb_gat=nb_gat) IF (debug) THEN WRITE (*,'(" Number of coordinate in the file : ",I2)') nb_coord WRITE (*,'(" Number of variables in the file : ",I2)') nb_var WRITE (*,'(" Number of global attributes in the file : ",I2)') nb_gat ENDIF ENDIF CALL bcast(nb_coord) CALL bcast(nb_var) CALL bcast(nb_gat) IF ( veget_year > 0 ) THEN IF (is_root_prc) & CALL flioinqv (fid,v_n="time_counter",l_ex=l_ex,nb_dims=nb_dim,len_dims=l_d_w) CALL bcast(l_d_w) tml=l_d_w(1) WRITE(numout,*) "tml =",tml ENDIF IF (is_root_prc) & CALL flioinqv (fid,v_n="lon",l_ex=l_ex,nb_dims=nb_dim,len_dims=l_d_w) CALL bcast(l_d_w) iml=l_d_w(1) WRITE(numout,*) "iml =",iml IF (is_root_prc) & CALL flioinqv (fid,v_n="lat",l_ex=l_ex,nb_dims=nb_dim,len_dims=l_d_w) CALL bcast(l_d_w) jml=l_d_w(1) WRITE(numout,*) "jml =",jml IF (is_root_prc) & CALL flioinqv (fid,v_n="veget",l_ex=l_ex,nb_dims=nb_dim,len_dims=l_d_w) CALL bcast(l_d_w) lml=l_d_w(1) IF (lml /= nvm) & CALL ipslerr (3,'slowproc_update', & & 'Problem with vegetation file for Land Use','lml /= nvm', & & '(number of pft must be equal)') ! ALLOC_ERR=-1 ALLOCATE(lat_lu(jml), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of lat_lu : ",ALLOC_ERR STOP ENDIF ALLOC_ERR=-1 ALLOCATE(lon_lu(iml), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of lon_lu : ",ALLOC_ERR STOP ENDIF IF ( veget_year > 0 ) THEN IF (tml > 0) THEN ALLOC_ERR=-1 ALLOCATE(itau(tml), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of itau : ",ALLOC_ERR STOP ENDIF ENDIF ! IF (is_root_prc) THEN IF (tml > 0) THEN CALL fliogstc (fid, t_axis=itau,x_axis=lon_lu,y_axis=lat_lu) IF (veget_year <= tml) THEN CALL fliogetv (fid,"time_counter",time_counter, start=(/ veget_year /), count=(/ 1 /)) WRITE(numout,*) "slowproc_update LAND_USE : the date read for vegetmax is ",time_counter ELSE CALL fliogetv (fid,"time_counter",time_counter, start=(/ tml /), count=(/ 1 /)) WRITE(numout,*) "slowproc_update LAND_USE : You try to update vegetmax with a the date greater than in the file." WRITE(numout,*) " We will keep the last one :",time_counter ENDIF ELSE CALL fliogstc (fid, x_axis=lon_lu,y_axis=lat_lu) ENDIF ENDIF ENDIF IF (tml > 0) THEN CALL bcast(itau) ENDIF CALL bcast(lon_lu) CALL bcast(lat_lu) ! ALLOC_ERR=-1 ALLOCATE(lat_ful(iml,jml), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of lat_ful : ",ALLOC_ERR STOP ENDIF ALLOC_ERR=-1 ALLOCATE(lon_ful(iml,jml), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of lon_ful : ",ALLOC_ERR STOP ENDIF ! DO ip=1,iml lon_ful(ip,:)=lon_lu(ip) ENDDO DO jp=1,jml lat_ful(:,jp)=lat_lu(jp) ENDDO ! ! WRITE(numout,*) 'Reading the LAND USE vegetation file' ! ALLOC_ERR=-1 ALLOCATE(vegmap(iml,jml,nvm,1), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of vegmap : ",ALLOC_ERR STOP ENDIF IF ( veget_year == 0 ) THEN IF (is_root_prc) THEN ALLOC_ERR=-1 ALLOCATE(vegmap_1(iml,jml,nvm), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of vegmap_1 : ",ALLOC_ERR STOP ENDIF ENDIF ENDIF ! !!$ CALL flinopen & !!$ & (filename, .FALSE., iml, jml, lml, lon_ful, lat_ful, & !!$ & lev_ful, tml, itau, date, dt, fid) !=> FATAL ERROR FROM ROUTINE flinopen ! --> No time axis found !MM modif TAG 1.4 : ! CALL flinget(fid, 'lon', iml, 0, 0, 0, 1, 1, lon_lu) ! CALL flinget(fid, 'lat', jml, 0, 0, 0, 1, 1, lat_lu) ! CALL flinget(fid, 'maxvegetfrac', iml, jml, nvm, tml, 1, 293, vegmap_lu) !FATAL ERROR FROM ROUTINE flinopen ! --> No variable lon ! We get only the right year !!$ CALL flinget(fid, 'maxvegetfrac', iml, jml, nvm, tml, veget_year, veget_year, vegmap) !!$ ! !!$ CALL flinclo(fid) IF (is_root_prc) & CALL flioinqv (fid,"maxvegetfrac",exv,nb_dims=nb_dim,len_dims=l_d_w,id_dims=i_d_w) CALL bcast(exv) CALL bcast(nb_dim) CALL bcast(l_d_w) CALL bcast(i_d_w) IF (exv) THEN IF (debug) THEN WRITE (numout,'(" Number of dimensions : ",I2)') nb_dim WRITE (numout,'(" Dimensions :",/,5(1X,I7,:))') l_d_w(1:nb_dim) WRITE (numout,'(" Identifiers :",/,5(1X,I7,:))') i_d_w(1:nb_dim) ENDIF ! IF ( veget_year > 0 ) THEN IF (is_root_prc) THEN IF (veget_year <= tml) THEN CALL fliogetv (fid,"maxvegetfrac", vegmap, start=(/ 1, 1, 1, veget_year /), count=(/ iml, jml, nvm, 1 /)) ELSE CALL fliogetv (fid,"maxvegetfrac", vegmap, start=(/ 1, 1, 1, tml /), count=(/ iml, jml, nvm, 1 /)) ENDIF ENDIF ELSE IF (is_root_prc) THEN CALL fliogetv (fid,"maxvegetfrac", vegmap_1, start=(/ 1, 1, 1 /), count=(/ iml, jml, nvm /)) vegmap(:,:,:,1)=vegmap_1(:,:,:) DEALLOCATE(vegmap_1) ENDIF ENDIF CALL bcast(vegmap) IF (is_root_prc) CALL flioclo (fid) ELSE CALL ipslerr (3,'slowproc_update', & & 'Problem with vegetation file for Land Use.', & & "Variable maxvegetfrac doesn't exist.", & & '(verify your land use file.)') ENDIF ! ! Mask of permitted variables. ! ALLOC_ERR=-1 ALLOCATE(mask(iml,jml), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of mask : ",ALLOC_ERR STOP ENDIF ! mask(:,:) = 0 DO ip=1,iml DO jp=1,jml sum_veg=SUM(vegmap(ip,jp,:,1)) IF ( sum_veg .GE. min_sechiba .AND. sum_veg .LE. 1.-1.e-7) THEN mask(ip,jp) = 1 IF (debug) THEN WRITE(numout,*) "update : SUM(vegmap(",ip,jp,")) = ",sum_veg ENDIF ELSEIF ( sum_veg .GT. 1.-1.e-7 .AND. sum_veg .LE. 2.) THEN ! normalization vegmap(ip,jp,:,1) = vegmap(ip,jp,:,1) / sum_veg mask(ip,jp) = 1 IF (debug) THEN WRITE(numout,*) "update : SUM(vegmap(",ip,jp,"))_c = ",SUM(vegmap(ip,jp,:,1)) ENDIF ENDIF ENDDO ENDDO ! ! ! The number of maximum vegetation map points in the GCM grid should ! also be computed and not imposed here. ! nbvmax = 200 ! callsign="Land Use Vegetation map" ! ok_interpol = .FALSE. DO WHILE ( .NOT. ok_interpol ) WRITE(numout,*) "Projection arrays for ",callsign," : " WRITE(numout,*) "nbvmax = ",nbvmax ! ALLOC_ERR=-1 ALLOCATE(sub_index(nbpt, nbvmax,2), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of sub_index : ",ALLOC_ERR STOP ENDIF sub_index(:,:,:)=0 ALLOC_ERR=-1 ALLOCATE(sub_area(nbpt, nbvmax), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of sub_area : ",ALLOC_ERR STOP ENDIF sub_area(:,:)=zero ! CALL aggregate_p(nbpt, lalo, neighbours, resolution, contfrac, & & iml, jml, lon_ful, lat_ful, mask, callsign, & & nbvmax, sub_index, sub_area, ok_interpol) ! IF ( .NOT. ok_interpol ) THEN DEALLOCATE(sub_area) DEALLOCATE(sub_index) ENDIF ! nbvmax = nbvmax * 2 ENDDO ! ! Compute the logical for partial (only anthropic) PTFs update IF (PRESENT(init)) THEN partial_update = control%ok_dgvm .AND. .NOT. init ELSE partial_update = control%ok_dgvm ENDIF ! IF ( .NOT. partial_update ) THEN ! veget_next(:,:)=zero ! DO ib = 1, nbpt idi=1 sumf=zero DO WHILE ( sub_area(ib,idi) > zero ) ip = sub_index(ib,idi,1) jp = sub_index(ib,idi,2) veget_next(ib,:) = veget_next(ib,:) + vegmap(ip,jp,:,1)*sub_area(ib,idi) sumf=sumf + sub_area(ib,idi) idi = idi +1 ENDDO !!$ ! !!$ ! Limit the smalest vegetation fraction to 0.5% !!$ ! !!$ DO jv = 1, nvm !!$ IF ( veget_next(ib,jv) .LT. min_vegfrac ) THEN !!$ veget_next(ib,jv) = zero !!$ ENDIF !!$ ENDDO ! ! Normalize ! IF (sumf > min_sechiba) THEN veget_next(ib,:) = veget_next(ib,:) / sumf ELSE WRITE(numout,*) "slowproc_update : No land point in the map for point ",& ib,",(",lalo(ib,1),",",lalo(ib,2),")" CALL ipslerr (2,'slowproc_update', & & 'Problem with vegetation file for Land Use.', & & "No land point in the map for point", & & 'Keep old values. (verify your land use file.)') !!$ CALL slowproc_nearest (iml, lon_ful, lat_ful, & !!$ lalo(ib,2), lalo(ib,1), inear) IF (PRESENT(init)) THEN IF (init) THEN veget_next(ib,:) = (/ 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0. /) ELSE veget_next(ib,:) = veget_last(ib,:) ENDIF ELSE veget_next(ib,:) = veget_last(ib,:) ENDIF ENDIF ! IF (debug) THEN WRITE(numout,*) "SUM(veget_next(",ib,")) = ",SUM(veget_next(ib,:)) ENDIF ENDDO ELSE DO ib = 1, nbpt ! last veget for this point sum_veg=SUM(veget_last(ib,:)) IF (debug) THEN WRITE(*,*) "SUM(veget_last(",ib,")) = ",sum_veg ENDIF ! ! If the DGVM is activated, only anthropiques PFT are utpdated, ! other are copied veget_next(ib,:) = veget_last(ib,:) ! ! natural ones are initialized to zero. DO jv = 2, nvm ! If the DGVM is activated, only anthropiques PFT are utpdated IF ( .NOT. natural(jv) ) THEN veget_next(ib,jv) = zero ENDIF ENDDO ! idi=1 sumf=zero DO WHILE ( sub_area(ib,idi) > zero ) ip = sub_index(ib,idi,1) jp = sub_index(ib,idi,2) ! If the DGVM is activated, only anthropic PFTs are utpdated DO jv = 2, nvm IF ( .NOT. natural(jv) ) THEN veget_next(ib,jv) = veget_next(ib,jv) + vegmap(ip,jp,jv,1)*sub_area(ib,idi) ENDIF ENDDO sumf=sumf + sub_area(ib,idi) idi = idi +1 ENDDO !!$ ! !!$ ! Limit the smalest vegetation fraction to 0.5% !!$ ! !!$ DO jv = 2, nvm !!$ ! On anthropic and natural PFTs ? !!$ IF ( veget_next(ib,jv) .LT. min_vegfrac ) THEN !!$ veget_next(ib,jv) = zero !!$ ENDIF !!$ ENDDO ! ! Normalize ! ! Proposition de Pierre : ! apres modification de la surface des PFTs anthropiques, ! on doit conserver la proportion des PFTs naturels. ! ie la somme des vegets est conservee ! et PFT naturel / (somme des vegets - somme des vegets anthropiques) ! est conservee. ! Sum veget_next = old (sum veget_next Naturel) + (sum veget_next Anthropic) ! = new (sum veget_next Naturel) + (sum veget_next Anthropic) ! a / (S-A) = e / (S-B) ; b/(S-A) = f/(S-B) IF (sumf > min_sechiba) THEN sumvAnthro_old = zero sumvAnthro = zero DO jv = 2, nvm IF ( .NOT. natural(jv) ) THEN veget_next(ib,jv) = veget_next(ib,jv) / sumf sumvAnthro = sumvAnthro + veget_last(ib,jv) sumvAnthro_old = sumvAnthro_old + veget_last(ib,jv) ENDIF ENDDO ! conservation : rapport = ( sum_veg - sumvAnthro ) / ( sum_veg - sumvAnthro_old ) veget_next(ib,1) = veget_last(ib,1) * rapport DO jv = 2, nvm IF ( .NOT. natural(jv) ) THEN veget_next(ib,jv) = veget_last(ib,jv) * rapport ENDIF ENDDO ! test IF ( ABS( SUM(veget_next(ib,:)) - sum_veg ) > EPSILON(un) ) THEN WRITE(numout,*) "No conservation of sum of veget for point ",ib,",(",lalo(ib,1),",",lalo(ib,2),")" WRITE(numout,*) "last sum of veget ",sum_veg," new sum of veget ",SUM(veget_next(ib,:))," error : ",& & SUM(veget_next(ib,:)) - sum_veg WRITE(numout,*) "Anthropic modificaztions : last ",sumvAnthro_old," new ",sumvAnthro CALL ipslerr (3,'slowproc_update', & & 'No conservation of sum of veget_next', & & "The sum of veget_next is different after reading Land Use map.", & & '(verify the dgvm case model.)') ENDIF ELSE WRITE(numout,*) "No land point in the map for point ",ib,",(",lalo(ib,1),",",lalo(ib,2),")" ! CALL ipslerr (3,'slowproc_update', & CALL ipslerr (2,'slowproc_update', & & 'Problem with vegetation file for Land Use.', & & "No land point in the map for point", & & '(verify your land use file.)') veget_next(ib,:) = veget_last(ib,:) ENDIF ENDDO ENDIF ! frac_nobio_next (:,:) = un ! !MM ! Work only for one nnobio !! (ie ice) DO inobio=1,nnobio DO jv=1,nvm ! DO ib = 1, nbpt frac_nobio_next(ib,inobio) = frac_nobio_next(ib,inobio) - veget_next(ib,jv) ENDDO ENDDO ! ENDDO ! DO ib = 1, nbpt sum_veg = SUM(veget_next(ib,:)) sum_nobio = SUM(frac_nobio_next(ib,:)) IF (sum_nobio < 0.) THEN frac_nobio_next(ib,:) = zero veget_next(ib,1) = veget_next(ib,1) - sum_nobio sum_veg = SUM(veget_next(ib,:)) ENDIF sumf = sum_veg + sum_nobio IF (sumf > min_sechiba) THEN veget_next(ib,:) = veget_next(ib,:) / sumf frac_nobio_next(ib,:) = frac_nobio_next(ib,:) / sumf norm=SUM(veget_next(ib,:))+SUM(frac_nobio_next(ib,:)) err=norm-un IF (debug) & WRITE(numout,*) "ib ",ib," SUM(veget_next(ib,:)+frac_nobio_next(ib,:))-un, sumf",err,sumf IF (abs(err) > -EPSILON(un)) THEN !MM 1.9.3 ! IF (abs(err) > 0.) THEN IF ( SUM(frac_nobio_next(ib,:)) > min_sechiba ) THEN frac_nobio_next(ib,1) = frac_nobio_next(ib,1) - err ELSE veget_next(ib,1) = veget_next(ib,1) - err ENDIF norm=SUM(veget_next(ib,:))+SUM(frac_nobio_next(ib,:)) err=norm-un IF (debug) & WRITE(numout,*) "ib ",ib," SUM(veget_next(ib,:)+frac_nobio_next(ib,:))-un",err IF (abs(err) > EPSILON(un)) THEN !MM 1.9.3 ! IF (abs(err) > 0.) THEN WRITE(numout,*) "update : Problem with point ",ib,",(",lalo(ib,1),",",lalo(ib,2),")" WRITE(numout,*) " err(sum-1.) = ",abs(err) CALL ipslerr (2,'slowproc_update', & & 'Problem with sum vegetation + sum fracnobio for Land Use.', & & "sum not equal to 1.", & & '(verify your land use file.)') ENDIF ENDIF ELSE WRITE(numout,*) "No vegetation nor frac_nobio for point ",ib,",(",lalo(ib,1),",",lalo(ib,2),")" WRITE(numout,*)"Replaced by bare_soil !! " veget_next(ib,1) = un veget_next(ib,2:nvm) = zero frac_nobio_next(ib,:) = zero !!$ CALL ipslerr (3,'slowproc_update', & !!$ & 'Problem with vegetation file for Land Use.', & !!$ & "No vegetation nor frac_nobio for point ", & !!$ & '(verify your land use file.)') ENDIF ENDDO ! WRITE(numout,*) 'slowproc_update : Interpolation Done' ! DEALLOCATE(vegmap) DEALLOCATE(lat_lu,lon_lu) DEALLOCATE(lat_ful,lon_ful) DEALLOCATE(mask) DEALLOCATE(sub_index,sub_area) ! RETURN ! END SUBROUTINE slowproc_update !! !! Interpolate the IGBP vegetation map to the grid of the model !MM TAG 1.6 model ! !! SUBROUTINE slowproc_interpol_OLD(nbpt, lalo, neighbours, resolution, veget, frac_nobio ) ! ! ! ! 0.1 INPUT ! INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longitudes (beware of the order !) INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) ! Vector of neighbours for each grid point ! (1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW) REAL(r_std), INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid-box in X and Y ! ! 0.2 OUTPUT ! REAL(r_std), INTENT(out) :: veget(nbpt,nvm) ! Vegetation fractions REAL(r_std), INTENT(out) :: frac_nobio(nbpt,nnobio) ! Fraction of the mesh which is covered by ice, lakes, ... ! ! 0.3 LOCAL ! REAL(r_std), PARAMETER :: min_sechiba=1.E-8 INTEGER(i_std), PARAMETER :: nolson = 94 ! Number of Olson classes ! ! CHARACTER(LEN=80) :: filename INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, ip, jp, vid REAL(r_std) :: lev(1), date, dt, coslat, pi INTEGER(i_std) :: itau(1) REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_ful, lon_ful, vegmap REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lon_up, lon_low, lat_up, lat_low INTEGER, DIMENSION(nbpt,nolson) :: n_origveg INTEGER, DIMENSION(nbpt) :: n_found REAL(r_std), DIMENSION(nbpt,nolson) :: frac_origveg REAL(r_std) :: vegcorr(nolson,nvm) REAL(r_std) :: nobiocorr(nolson,nnobio) CHARACTER(LEN=80) :: meter REAL(r_std) :: prog, sumf LOGICAL :: found INTEGER :: idi, ilast, ii, jv, inear, iprog REAL(r_std) :: domaine_lon_min, domaine_lon_max, domaine_lat_min, domaine_lat_max ! pi = 4. * ATAN(1.) ! CALL get_vegcorr (nolson,vegcorr,nobiocorr) ! !Config Key = VEGETATION_FILE !Config Desc = Name of file from which the vegetation map is to be read !Config If = !IMPOSE_VEG !Config Def = ../surfmap/carteveg5km.nc !Config Help = The name of the file to be opened to read the vegetation !Config map is to be given here. Usualy SECHIBA runs with a 5kmx5km !Config map which is derived from the IGBP one. We assume that we have !Config a classification in 87 types. This is Olson modified by Viovy. ! filename = '../surfmap/carteveg5km.nc' CALL getin_p('VEGETATION_FILE',filename) ! if (is_root_prc) CALL flininfo(filename, iml, jml, lml, tml, fid) CALL bcast(iml) CALL bcast(jml) CALL bcast(lml) CALL bcast(tml) ! ! ALLOCATE(lat_ful(iml)) ALLOCATE(lon_ful(iml)) ALLOCATE(vegmap(iml)) ! WRITE(numout,*) 'Reading the vegetation file' ! IF (is_root_prc) THEN CALL flinget(fid, 'longitude', iml, jml, lml, tml, 1, 1, lon_ful) CALL flinget(fid, 'latitude', iml, jml, lml, tml, 1, 1, lat_ful) CALL flinget(fid, 'vegetation_map', iml, jml, lml, tml, 1, 1, vegmap) ! CALL flinclo(fid) ENDIF CALL bcast(lon_ful) CALL bcast(lat_ful) CALL bcast(vegmap) ! IF (MAXVAL(vegmap) .LT. nolson) THEN WRITE(*,*) 'WARNING -- WARNING' WRITE(*,*) 'The vegetation map has to few vegetation types.' WRITE(*,*) 'If you are lucky it will work but please check' ELSE IF ( MAXVAL(vegmap) .GT. nolson) THEN WRITE(*,*) 'More vegetation types in file than the code can' WRITE(*,*) 'deal with.: ', MAXVAL(vegmap), nolson STOP 'slowproc_interpol' ENDIF ! ALLOCATE(lon_up(nbpt)) ALLOCATE(lon_low(nbpt)) ALLOCATE(lat_up(nbpt)) ALLOCATE(lat_low(nbpt)) ! DO ib =1, nbpt ! ! We find the 4 limits of the grid-box. As we transform the resolution of the model ! into longitudes and latitudes we do not have the problem of periodicity. ! coslat is a help variable here ! ! coslat = MAX(COS(lalo(ib,1) * pi/180. ), 0.001 )*pi/180. * R_Earth ! lon_up(ib) = lalo(ib,2) + resolution(ib,1)/(2.0*coslat) lon_low(ib) = lalo(ib,2) - resolution(ib,1)/(2.0*coslat) ! coslat = pi/180. * R_Earth ! lat_up(ib) = lalo(ib,1) + resolution(ib,2)/(2.0*coslat) lat_low(ib) = lalo(ib,1) - resolution(ib,2)/(2.0*coslat) ! ! veget(ib,:) = zero frac_nobio (ib,:) = zero ! ENDDO ! ! Get the limits of the integration domaine so that we can speed up the calculations ! domaine_lon_min = MINVAL(lon_low) domaine_lon_max = MAXVAL(lon_up) domaine_lat_min = MINVAL(lat_low) domaine_lat_max = MAXVAL(lat_up) ! !!$ WRITE(*,*) 'DOMAINE lon :', domaine_lon_min, domaine_lon_max !!$ WRITE(*,*) 'DOMAINE lat :', domaine_lat_min, domaine_lat_max ! ! Ensure that the fine grid covers the whole domain WHERE ( lon_ful(:) .LT. domaine_lon_min ) lon_ful(:) = lon_ful(:) + 360. ENDWHERE ! WHERE ( lon_ful(:) .GT. domaine_lon_max ) lon_ful(:) = lon_ful(:) - 360. ENDWHERE ! WRITE(numout,*) 'Interpolating the vegetation map :' WRITE(numout,'(2a40)')'0%--------------------------------------', & & '------------------------------------100%' ! ilast = 1 n_origveg(:,:) = 0 ! DO ip=1,iml ! ! Give a progress meter ! ! prog = ip/float(iml)*79. ! IF ( ABS(prog - NINT(prog)) .LT. 1/float(iml)*79. ) THEN ! meter(NINT(prog)+1:NINT(prog)+1) = 'x' ! WRITE(numout, advance="no", FMT='(a)') ACHAR(13) ! WRITE(numout, advance="no", FMT='(a80)') meter ! ENDIF iprog = NINT(float(ip)/float(iml)*79.) - NINT(float(ip-1)/float(iml)*79.) IF ( iprog .NE. 0 ) THEN WRITE(numout,'(a1,$)') 'x' ENDIF ! ! Only start looking for its place in the smaler grid if we are within the domaine ! That should speed up things ! ! IF ( ( lon_ful(ip) .GE. domaine_lon_min ) .AND. & ( lon_ful(ip) .LE. domaine_lon_max ) .AND. & ( lat_ful(ip) .GE. domaine_lat_min ) .AND. & ( lat_ful(ip) .LE. domaine_lat_max ) ) THEN ! ! look for point on GCM grid which this point on fine grid belongs to. ! First look at the point on the model grid where we arrived just before. There is ! a good chace that neighbouring points on the fine grid fall into the same model ! grid box. ! ! ! THERE IS A BUG HERE !!! IF THE GCM GRID SITS ON THE DATE LINE WE WILL HAVE FOR INSTANCE ! LON_LOW = -182 AND LON_UP = -178. THUS WE WILL ONLY PICK UP HALF THE POINTS NEEDED. ! IF ( ( lon_ful(ip) .GT. lon_low(ilast) ) .AND. & ( lon_ful(ip) .LT. lon_up(ilast) ) .AND. & ( lat_ful(ip) .GT. lat_low(ilast) ) .AND. & ( lat_ful(ip) .LT. lat_up(ilast) ) ) THEN ! ! We were lucky ! n_origveg(ilast,NINT(vegmap(ip))) = n_origveg(ilast,NINT(vegmap(ip))) + 1 ! ELSE ! ! Otherwise, look everywhere. ! Begin close to last grid point. ! found = .FALSE. idi = 1 ! DO WHILE ( (idi .LT. nbpt) .AND. ( .NOT. found ) ) ! ! forward and backward ! DO ii = 1,2 ! IF ( ii .EQ. 1 ) THEN ib = ilast - idi ELSE ib = ilast + idi ENDIF ! IF ( ( ib .GE. 1 ) .AND. ( ib .LE. nbpt ) ) THEN IF ( ( lon_ful(ip) .GT. lon_low(ib) ) .AND. & ( lon_ful(ip) .LT. lon_up(ib) ) .AND. & ( lat_ful(ip) .GT. lat_low(ib) ) .AND. & ( lat_ful(ip) .LT. lat_up(ib) ) ) THEN ! n_origveg(ib,NINT(vegmap(ip))) = n_origveg(ib,NINT(vegmap(ip))) + 1 ilast = ib found = .TRUE. ! ENDIF ENDIF ! ENDDO ! idi = idi + 1 ! ENDDO ! ENDIF ! lucky/not lucky ! ENDIF ! in the domain ENDDO ! ! Now we know how many points of which Olson type from the fine grid fall ! into each box of the (coarse) model grid: n_origveg(nbpt,nolson) ! ! ! determine number of points of the fine grid which fall into each box of the ! coarse grid ! DO ib = 1, nbpt n_found(ib) = SUM( n_origveg(ib,:) ) ENDDO ! ! determine fraction of Olson vegetation type in each box of the coarse grid ! DO vid = 1, nolson WHERE ( n_found(:) .GT. 0 ) frac_origveg(:,vid) = REAL(n_origveg(:,vid),r_std) / REAL(n_found(:),r_std) ELSEWHERE frac_origveg(:,vid) = zero ENDWHERE ENDDO ! ! now finally calculate coarse vegetation map ! Find which model vegetation corresponds to each Olson type ! DO vid = 1, nolson ! DO jv = 1, nvm veget(:,jv) = veget(:,jv) + frac_origveg(:,vid) * vegcorr(vid,jv) ENDDO ! DO jv = 1, nnobio frac_nobio(:,jv) = frac_nobio(:,jv) + frac_origveg(:,vid) * nobiocorr(vid,jv) ENDDO ! ENDDO ! ! WRITE(numout,*) WRITE(numout,*) 'Interpolation Done' ! ! Clean up the point of the map ! DO ib = 1, nbpt ! ! Let us see if all points found something in the 5km map ! ! IF ( n_found(ib) .EQ. 0 ) THEN ! ! Now we need to handle some exceptions ! IF ( lalo(ib,1) .LT. -56.0) THEN ! Antartica frac_nobio(ib,:) = zero frac_nobio(ib,iice) = un veget(ib,:) = zero ! ELSE IF ( lalo(ib,1) .GT. 70.0) THEN ! Artica frac_nobio(ib,:) = zero frac_nobio(ib,iice) = un veget(ib,:) = zero ! ELSE IF ( lalo(ib,1) .GT. 55.0 .AND. lalo(ib,2) .GT. -65.0 .AND. lalo(ib,2) .LT. -20.0) THEN ! Greenland frac_nobio(ib,:) = zero frac_nobio(ib,iice) = un veget(ib,:) = zero ! ELSE ! WRITE(numout,*) 'PROBLEM, no point in the 5km map found for this grid box' WRITE(numout,*) 'Longitude range : ', lon_low(ib), lon_up(ib) WRITE(numout,*) 'Latitude range : ', lat_low(ib), lat_up(ib) ! WRITE(numout,*) 'Looking for nearest point on the 5 km map' CALL slowproc_nearest (iml, lon_ful, lat_ful, & lalo(ib,2), lalo(ib,1), inear) WRITE(numout,*) 'Coordinates of the nearest point:', & lon_ful(inear),lat_ful(inear) ! DO jv = 1, nvm veget(ib,jv) = vegcorr(NINT(vegmap(inear)),jv) ENDDO ! DO jv = 1, nnobio frac_nobio(ib,jv) = nobiocorr(NINT(vegmap(inear)),jv) ENDDO ! ENDIF ! ENDIF ! ! ! Limit the smalest vegetation fraction to 0.5% ! DO vid = 1, nvm IF ( veget(ib,vid) .LT. min_vegfrac ) THEN veget(ib,vid) = zero ENDIF ENDDO ! sumf = SUM(frac_nobio(ib,:))+SUM(veget(ib,:)) frac_nobio(ib,:) = frac_nobio(ib,:)/sumf veget(ib,:) = veget(ib,:)/sumf ! ! ENDDO ! DEALLOCATE(lon_up) DEALLOCATE(lon_low) DEALLOCATE(lat_up) DEALLOCATE(lat_low) DEALLOCATE(lat_ful) DEALLOCATE(lon_ful) DEALLOCATE(vegmap) ! RETURN ! END SUBROUTINE slowproc_interpol_OLD !! !! Interpolate the IGBP vegetation map to the grid of the model !! SUBROUTINE slowproc_interpol_NEW(nbpt, lalo, neighbours, resolution, contfrac, veget, frac_nobio ) ! ! ! ! 0.1 INPUT ! INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longitudes (beware of the order !) INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) ! Vector of neighbours for each grid point ! (1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW) REAL(r_std), INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid-box in X and Y REAL(r_std),DIMENSION (nbpt), INTENT (in) :: contfrac !! Fraction of continent in the grid ! ! 0.2 OUTPUT ! REAL(r_std), INTENT(out) :: veget(nbpt,nvm) ! Vegetation fractions REAL(r_std), INTENT(out) :: frac_nobio(nbpt,nnobio) ! Fraction of the mesh which is covered by ice, lakes, ... ! LOGICAL :: ok_interpol ! optionnal return of aggregate_vec ! ! 0.3 LOCAL ! INTEGER(i_std), PARAMETER :: nolson = 94 ! Number of Olson classes ! ! CHARACTER(LEN=80) :: filename INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, ip, vid REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_ful, lon_ful, vegmap REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: sub_area INTEGER(i_std),ALLOCATABLE, DIMENSION(:,:) :: sub_index REAL(r_std), DIMENSION(nbpt,nolson) :: n_origveg REAL(r_std), DIMENSION(nbpt) :: n_found REAL(r_std), DIMENSION(nbpt,nolson) :: frac_origveg REAL(r_std) :: vegcorr(nolson,nvm) REAL(r_std) :: nobiocorr(nolson,nnobio) CHARACTER(LEN=40) :: callsign REAL(r_std) :: sumf, resol_lon, resol_lat INTEGER(i_std) :: idi, jv, inear, nbvmax ! INTEGER :: ALLOC_ERR ! n_origveg(:,:) = zero n_found(:) = zero ! CALL get_vegcorr (nolson,vegcorr,nobiocorr) ! !Config Key = VEGETATION_FILE !Config Desc = Name of file from which the vegetation map is to be read !Config If = !IMPOSE_VEG !Config If = !LAND_USE !Config Def = ../surfmap/carteveg5km.nc !Config Help = The name of the file to be opened to read the vegetation !Config map is to be given here. Usualy SECHIBA runs with a 5kmx5km !Config map which is derived from the IGBP one. We assume that we have !Config a classification in 87 types. This is Olson modified by Viovy. ! filename = '../surfmap/carteveg5km.nc' CALL getin_p('VEGETATION_FILE',filename) ! if (is_root_prc) CALL flininfo(filename, iml, jml, lml, tml, fid) CALL bcast(iml) CALL bcast(jml) CALL bcast(lml) CALL bcast(tml) ! ! ALLOC_ERR=-1 ALLOCATE(lat_ful(iml), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of lat_ful : ",ALLOC_ERR STOP ENDIF ALLOC_ERR=-1 ALLOCATE(lon_ful(iml), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of lon_ful : ",ALLOC_ERR STOP ENDIF ALLOC_ERR=-1 ALLOCATE(vegmap(iml), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of vegmap : ",ALLOC_ERR STOP ENDIF ! WRITE(numout,*) 'Reading the OLSON type vegetation file' ! IF (is_root_prc) THEN CALL flinget(fid, 'longitude', iml, jml, lml, tml, 1, 1, lon_ful) CALL flinget(fid, 'latitude', iml, jml, lml, tml, 1, 1, lat_ful) CALL flinget(fid, 'vegetation_map', iml, jml, lml, tml, 1, 1, vegmap) ! CALL flinclo(fid) ENDIF CALL bcast(lon_ful) CALL bcast(lat_ful) CALL bcast(vegmap) ! IF (MAXVAL(vegmap) .LT. nolson) THEN WRITE(numout,*) 'WARNING -- WARNING' WRITE(numout,*) 'The vegetation map has to few vegetation types.' WRITE(numout,*) 'If you are lucky it will work but please check' ELSE IF ( MAXVAL(vegmap) .GT. nolson) THEN WRITE(numout,*) 'More vegetation types in file than the code can' WRITE(numout,*) 'deal with.: ', MAXVAL(vegmap), nolson STOP 'slowproc_interpol' ENDIF ! ! Some assumptions on the vegetation file. This information should be ! be computed or read from the file. ! It is the reolution in meters of the grid of the vegetation file. ! resol_lon = 5000. resol_lat = 5000. ! ! The number of maximum vegetation map points in the GCM grid should ! also be computed and not imposed here. nbvmax = iml/nbpt ! callsign="Vegetation map" ! ok_interpol = .FALSE. DO WHILE ( .NOT. ok_interpol ) WRITE(numout,*) "Projection arrays for ",callsign," : " WRITE(numout,*) "nbvmax = ",nbvmax ! ALLOC_ERR=-1 ALLOCATE(sub_index(nbpt, nbvmax), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of sub_index : ",ALLOC_ERR STOP ENDIF sub_index(:,:)=0 ALLOC_ERR=-1 ALLOCATE(sub_area(nbpt, nbvmax), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of sub_area : ",ALLOC_ERR STOP ENDIF sub_area(:,:)=zero ! CALL aggregate_p (nbpt, lalo, neighbours, resolution, contfrac, & & iml, lon_ful, lat_ful, resol_lon, resol_lat, callsign, & & nbvmax, sub_index, sub_area, ok_interpol) ! IF ( .NOT. ok_interpol ) THEN DEALLOCATE(sub_area) DEALLOCATE(sub_index) ! nbvmax = nbvmax * 2 ELSE ! DO ib = 1, nbpt idi=1 DO WHILE ( sub_area(ib,idi) > zero ) ip = sub_index(ib,idi) n_origveg(ib,NINT(vegmap(ip))) = n_origveg(ib,NINT(vegmap(ip))) + sub_area(ib,idi) n_found(ib) = n_found(ib) + sub_area(ib,idi) idi = idi +1 ENDDO ENDDO ! ENDIF ENDDO ! ! Now we know how many points of which Olson type from the fine grid fall ! into each box of the (coarse) model grid: n_origveg(nbpt,nolson) ! ! ! determine fraction of Olson vegetation type in each box of the coarse grid ! DO vid = 1, nolson WHERE ( n_found(:) .GT. 0 ) frac_origveg(:,vid) = n_origveg(:,vid) / n_found(:) ELSEWHERE frac_origveg(:,vid) = zero ENDWHERE ENDDO ! ! now finally calculate coarse vegetation map ! Find which model vegetation corresponds to each Olson type ! veget(:,:) = zero frac_nobio(:,:) = zero ! DO vid = 1, nolson ! DO jv = 1, nvm veget(:,jv) = veget(:,jv) + frac_origveg(:,vid) * vegcorr(vid,jv) ENDDO ! DO jv = 1, nnobio frac_nobio(:,jv) = frac_nobio(:,jv) + frac_origveg(:,vid) * nobiocorr(vid,jv) ENDDO ! ENDDO ! WRITE(numout,*) 'slowproc_interpol : Interpolation Done' ! ! Clean up the point of the map ! DO ib = 1, nbpt ! ! Let us see if all points found something in the 5km map ! ! IF ( n_found(ib) .EQ. 0 ) THEN ! ! Now we need to handle some exceptions ! IF ( lalo(ib,1) .LT. -56.0) THEN ! Antartica frac_nobio(ib,:) = zero frac_nobio(ib,iice) = un veget(ib,:) = zero ! ELSE IF ( lalo(ib,1) .GT. 70.0) THEN ! Artica frac_nobio(ib,:) = zero frac_nobio(ib,iice) = un veget(ib,:) = zero ! ELSE IF ( lalo(ib,1) .GT. 55.0 .AND. lalo(ib,2) .GT. -65.0 .AND. lalo(ib,2) .LT. -20.0) THEN ! Greenland frac_nobio(ib,:) = zero frac_nobio(ib,iice) = un veget(ib,:) = zero ! ELSE ! WRITE(numout,*) 'PROBLEM, no point in the 5km map found for this grid box',ib WRITE(numout,*) 'Longitude range : ', lalo(ib,2) WRITE(numout,*) 'Latitude range : ', lalo(ib,1) ! WRITE(numout,*) 'Looking for nearest point on the 5 km map' CALL slowproc_nearest (iml, lon_ful, lat_ful, & lalo(ib,2), lalo(ib,1), inear) WRITE(numout,*) 'Coordinates of the nearest point:', & lon_ful(inear),lat_ful(inear) ! DO jv = 1, nvm veget(ib,jv) = vegcorr(NINT(vegmap(inear)),jv) ENDDO ! DO jv = 1, nnobio frac_nobio(ib,jv) = nobiocorr(NINT(vegmap(inear)),jv) ENDDO ! ENDIF ! ENDIF ! ! ! Limit the smalest vegetation fraction to 0.5% ! DO vid = 1, nvm IF ( veget(ib,vid) .LT. min_vegfrac ) THEN veget(ib,vid) = zero ENDIF ENDDO ! sumf = SUM(frac_nobio(ib,:))+SUM(veget(ib,:)) frac_nobio(ib,:) = frac_nobio(ib,:)/sumf veget(ib,:) = veget(ib,:)/sumf ! ! ENDDO ! DEALLOCATE(vegmap) DEALLOCATE(lat_ful, lon_ful) DEALLOCATE(sub_index) DEALLOCATE(sub_area) ! RETURN ! END SUBROUTINE slowproc_interpol_NEW !! !! Interpolate the IGBP vegetation map to the grid of the model !MM TAG 1.6 model ! !! SUBROUTINE slowproc_interpol_OLD_g(nbpt, lalo, neighbours, resolution, veget, frac_nobio ) ! ! ! ! 0.1 INPUT ! INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longitudes (beware of the order !) INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) ! Vector of neighbours for each grid point ! (1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW) REAL(r_std), INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid-box in X and Y ! ! 0.2 OUTPUT ! REAL(r_std), INTENT(out) :: veget(nbpt,nvm) ! Vegetation fractions REAL(r_std), INTENT(out) :: frac_nobio(nbpt,nnobio) ! Fraction of the mesh which is covered by ice, lakes, ... ! ! 0.3 LOCAL ! REAL(r_std), PARAMETER :: min_sechiba=1.E-8 INTEGER(i_std), PARAMETER :: nolson = 94 ! Number of Olson classes ! ! CHARACTER(LEN=80) :: filename INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, ip, jp, vid REAL(r_std) :: lev(1), date, dt, coslat, pi INTEGER(i_std) :: itau(1) REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_ful, lon_ful, vegmap REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lon_up, lon_low, lat_up, lat_low INTEGER, DIMENSION(nbpt,nolson) :: n_origveg INTEGER, DIMENSION(nbpt) :: n_found REAL(r_std), DIMENSION(nbpt,nolson) :: frac_origveg REAL(r_std) :: vegcorr(nolson,nvm) REAL(r_std) :: nobiocorr(nolson,nnobio) CHARACTER(LEN=80) :: meter REAL(r_std) :: prog, sumf LOGICAL :: found INTEGER :: idi, ilast, ii, jv, inear, iprog REAL(r_std) :: domaine_lon_min, domaine_lon_max, domaine_lat_min, domaine_lat_max ! pi = 4. * ATAN(1.) ! CALL get_vegcorr (nolson,vegcorr,nobiocorr) ! !Config Key = VEGETATION_FILE !Config Desc = Name of file from which the vegetation map is to be read !Config If = !IMPOSE_VEG !Config Def = ../surfmap/carteveg5km.nc !Config Help = The name of the file to be opened to read the vegetation !Config map is to be given here. Usualy SECHIBA runs with a 5kmx5km !Config map which is derived from the IGBP one. We assume that we have !Config a classification in 87 types. This is Olson modified by Viovy. ! filename = '../surfmap/carteveg5km.nc' CALL getin('VEGETATION_FILE',filename) ! CALL flininfo(filename, iml, jml, lml, tml, fid) ! ! ALLOCATE(lat_ful(iml)) ALLOCATE(lon_ful(iml)) ALLOCATE(vegmap(iml)) ! WRITE(numout,*) 'Reading the vegetation file' ! CALL flinget(fid, 'longitude', iml, jml, lml, tml, 1, 1, lon_ful) CALL flinget(fid, 'latitude', iml, jml, lml, tml, 1, 1, lat_ful) CALL flinget(fid, 'vegetation_map', iml, jml, lml, tml, 1, 1, vegmap) ! CALL flinclo(fid) ! IF (MAXVAL(vegmap) .LT. nolson) THEN WRITE(*,*) 'WARNING -- WARNING' WRITE(*,*) 'The vegetation map has to few vegetation types.' WRITE(*,*) 'If you are lucky it will work but please check' ELSE IF ( MAXVAL(vegmap) .GT. nolson) THEN WRITE(*,*) 'More vegetation types in file than the code can' WRITE(*,*) 'deal with.: ', MAXVAL(vegmap), nolson STOP 'slowproc_interpol' ENDIF ! ALLOCATE(lon_up(nbpt)) ALLOCATE(lon_low(nbpt)) ALLOCATE(lat_up(nbpt)) ALLOCATE(lat_low(nbpt)) ! DO ib =1, nbpt ! ! We find the 4 limits of the grid-box. As we transform the resolution of the model ! into longitudes and latitudes we do not have the problem of periodicity. ! coslat is a help variable here ! ! coslat = MAX(COS(lalo(ib,1) * pi/180. ), 0.001 )*pi/180. * R_Earth ! lon_up(ib) = lalo(ib,2) + resolution(ib,1)/(2.0*coslat) lon_low(ib) = lalo(ib,2) - resolution(ib,1)/(2.0*coslat) ! coslat = pi/180. * R_Earth ! lat_up(ib) = lalo(ib,1) + resolution(ib,2)/(2.0*coslat) lat_low(ib) = lalo(ib,1) - resolution(ib,2)/(2.0*coslat) ! ! veget(ib,:) = zero frac_nobio (ib,:) = zero ! ENDDO ! ! Get the limits of the integration domaine so that we can speed up the calculations ! domaine_lon_min = MINVAL(lon_low) domaine_lon_max = MAXVAL(lon_up) domaine_lat_min = MINVAL(lat_low) domaine_lat_max = MAXVAL(lat_up) ! !!$ WRITE(*,*) 'DOMAINE lon :', domaine_lon_min, domaine_lon_max !!$ WRITE(*,*) 'DOMAINE lat :', domaine_lat_min, domaine_lat_max ! ! Ensure that the fine grid covers the whole domain WHERE ( lon_ful(:) .LT. domaine_lon_min ) lon_ful(:) = lon_ful(:) + 360. ENDWHERE ! WHERE ( lon_ful(:) .GT. domaine_lon_max ) lon_ful(:) = lon_ful(:) - 360. ENDWHERE ! WRITE(numout,*) 'Interpolating the vegetation map :' WRITE(numout,'(2a40)')'0%--------------------------------------', & & '------------------------------------100%' ! ilast = 1 n_origveg(:,:) = 0 ! DO ip=1,iml ! ! Give a progress meter ! ! prog = ip/float(iml)*79. ! IF ( ABS(prog - NINT(prog)) .LT. 1/float(iml)*79. ) THEN ! meter(NINT(prog)+1:NINT(prog)+1) = 'x' ! WRITE(numout, advance="no", FMT='(a)') ACHAR(13) ! WRITE(numout, advance="no", FMT='(a80)') meter ! ENDIF iprog = NINT(float(ip)/float(iml)*79.) - NINT(float(ip-1)/float(iml)*79.) IF ( iprog .NE. 0 ) THEN WRITE(numout,'(a1,$)') 'x' ENDIF ! ! Only start looking for its place in the smaler grid if we are within the domaine ! That should speed up things ! ! IF ( ( lon_ful(ip) .GE. domaine_lon_min ) .AND. & ( lon_ful(ip) .LE. domaine_lon_max ) .AND. & ( lat_ful(ip) .GE. domaine_lat_min ) .AND. & ( lat_ful(ip) .LE. domaine_lat_max ) ) THEN ! ! look for point on GCM grid which this point on fine grid belongs to. ! First look at the point on the model grid where we arrived just before. There is ! a good chace that neighbouring points on the fine grid fall into the same model ! grid box. ! ! ! THERE IS A BUG HERE !!! IF THE GCM GRID SITS ON THE DATE LINE WE WILL HAVE FOR INSTANCE ! LON_LOW = -182 AND LON_UP = -178. THUS WE WILL ONLY PICK UP HALF THE POINTS NEEDED. ! IF ( ( lon_ful(ip) .GT. lon_low(ilast) ) .AND. & ( lon_ful(ip) .LT. lon_up(ilast) ) .AND. & ( lat_ful(ip) .GT. lat_low(ilast) ) .AND. & ( lat_ful(ip) .LT. lat_up(ilast) ) ) THEN ! ! We were lucky ! n_origveg(ilast,NINT(vegmap(ip))) = n_origveg(ilast,NINT(vegmap(ip))) + 1 ! ELSE ! ! Otherwise, look everywhere. ! Begin close to last grid point. ! found = .FALSE. idi = 1 ! DO WHILE ( (idi .LT. nbpt) .AND. ( .NOT. found ) ) ! ! forward and backward ! DO ii = 1,2 ! IF ( ii .EQ. 1 ) THEN ib = ilast - idi ELSE ib = ilast + idi ENDIF ! IF ( ( ib .GE. 1 ) .AND. ( ib .LE. nbpt ) ) THEN IF ( ( lon_ful(ip) .GT. lon_low(ib) ) .AND. & ( lon_ful(ip) .LT. lon_up(ib) ) .AND. & ( lat_ful(ip) .GT. lat_low(ib) ) .AND. & ( lat_ful(ip) .LT. lat_up(ib) ) ) THEN ! n_origveg(ib,NINT(vegmap(ip))) = n_origveg(ib,NINT(vegmap(ip))) + 1 ilast = ib found = .TRUE. ! ENDIF ENDIF ! ENDDO ! idi = idi + 1 ! ENDDO ! ENDIF ! lucky/not lucky ! ENDIF ! in the domain ENDDO ! ! Now we know how many points of which Olson type from the fine grid fall ! into each box of the (coarse) model grid: n_origveg(nbpt,nolson) ! ! ! determine number of points of the fine grid which fall into each box of the ! coarse grid ! DO ib = 1, nbpt n_found(ib) = SUM( n_origveg(ib,:) ) ENDDO ! ! determine fraction of Olson vegetation type in each box of the coarse grid ! DO vid = 1, nolson WHERE ( n_found(:) .GT. 0 ) frac_origveg(:,vid) = REAL(n_origveg(:,vid),r_std) / REAL(n_found(:),r_std) ELSEWHERE frac_origveg(:,vid) = zero ENDWHERE ENDDO ! ! now finally calculate coarse vegetation map ! Find which model vegetation corresponds to each Olson type ! DO vid = 1, nolson ! DO jv = 1, nvm veget(:,jv) = veget(:,jv) + frac_origveg(:,vid) * vegcorr(vid,jv) ENDDO ! DO jv = 1, nnobio frac_nobio(:,jv) = frac_nobio(:,jv) + frac_origveg(:,vid) * nobiocorr(vid,jv) ENDDO ! ENDDO ! ! WRITE(numout,*) WRITE(numout,*) 'Interpolation Done' ! ! Clean up the point of the map ! DO ib = 1, nbpt ! ! Let us see if all points found something in the 5km map ! ! IF ( n_found(ib) .EQ. 0 ) THEN ! ! Now we need to handle some exceptions ! IF ( lalo(ib,1) .LT. -56.0) THEN ! Antartica frac_nobio(ib,:) = zero frac_nobio(ib,iice) = un veget(ib,:) = zero ! ELSE IF ( lalo(ib,1) .GT. 70.0) THEN ! Artica frac_nobio(ib,:) = zero frac_nobio(ib,iice) = un veget(ib,:) = zero ! ELSE IF ( lalo(ib,1) .GT. 55.0 .AND. lalo(ib,2) .GT. -65.0 .AND. lalo(ib,2) .LT. -20.0) THEN ! Greenland frac_nobio(ib,:) = zero frac_nobio(ib,iice) = un veget(ib,:) = zero ! ELSE ! WRITE(numout,*) 'PROBLEM, no point in the 5km map found for this grid box' WRITE(numout,*) 'Longitude range : ', lon_low(ib), lon_up(ib) WRITE(numout,*) 'Latitude range : ', lat_low(ib), lat_up(ib) ! WRITE(numout,*) 'Looking for nearest point on the 5 km map' CALL slowproc_nearest (iml, lon_ful, lat_ful, & lalo(ib,2), lalo(ib,1), inear) WRITE(numout,*) 'Coordinates of the nearest point:', & lon_ful(inear),lat_ful(inear) ! DO jv = 1, nvm veget(ib,jv) = vegcorr(NINT(vegmap(inear)),jv) ENDDO ! DO jv = 1, nnobio frac_nobio(ib,jv) = nobiocorr(NINT(vegmap(inear)),jv) ENDDO ! ENDIF ! ENDIF ! ! ! Limit the smalest vegetation fraction to 0.5% ! DO vid = 1, nvm IF ( veget(ib,vid) .LT. min_vegfrac ) THEN veget(ib,vid) = zero ENDIF ENDDO ! sumf = SUM(frac_nobio(ib,:))+SUM(veget(ib,:)) frac_nobio(ib,:) = frac_nobio(ib,:)/sumf veget(ib,:) = veget(ib,:)/sumf ! ! ENDDO ! DEALLOCATE(lon_up) DEALLOCATE(lon_low) DEALLOCATE(lat_up) DEALLOCATE(lat_low) DEALLOCATE(lat_ful) DEALLOCATE(lon_ful) DEALLOCATE(vegmap) ! RETURN ! END SUBROUTINE slowproc_interpol_OLD_g !! !! Interpolate the IGBP vegetation map to the grid of the model !! SUBROUTINE slowproc_interpol_NEW_g(nbpt, lalo, neighbours, resolution, contfrac, veget, frac_nobio ) ! ! ! ! 0.1 INPUT ! INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longitudes (beware of the order !) INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) ! Vector of neighbours for each grid point ! (1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW) REAL(r_std), INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid-box in X and Y REAL(r_std),DIMENSION (nbpt), INTENT (in) :: contfrac !! Fraction of continent in the grid ! ! 0.2 OUTPUT ! REAL(r_std), INTENT(out) :: veget(nbpt,nvm) ! Vegetation fractions REAL(r_std), INTENT(out) :: frac_nobio(nbpt,nnobio) ! Fraction of the mesh which is covered by ice, lakes, ... ! LOGICAL :: ok_interpol ! optionnal return of aggregate_vec ! ! 0.3 LOCAL ! INTEGER(i_std), PARAMETER :: nolson = 94 ! Number of Olson classes ! ! CHARACTER(LEN=80) :: filename INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, ip, vid REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_ful, lon_ful, vegmap REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: sub_area INTEGER(i_std),ALLOCATABLE, DIMENSION(:,:) :: sub_index REAL(r_std), DIMENSION(nbpt,nolson) :: n_origveg REAL(r_std), DIMENSION(nbpt) :: n_found REAL(r_std), DIMENSION(nbpt,nolson) :: frac_origveg REAL(r_std) :: vegcorr(nolson,nvm) REAL(r_std) :: nobiocorr(nolson,nnobio) CHARACTER(LEN=40) :: callsign REAL(r_std) :: sumf, resol_lon, resol_lat INTEGER(i_std) :: idi, jv, inear, nbvmax ! INTEGER :: ALLOC_ERR ! n_origveg(:,:) = zero n_found(:) = zero ! CALL get_vegcorr (nolson,vegcorr,nobiocorr) ! !Config Key = VEGETATION_FILE !Config Desc = Name of file from which the vegetation map is to be read !Config If = !IMPOSE_VEG !Config If = !LAND_USE !Config Def = ../surfmap/carteveg5km.nc !Config Help = The name of the file to be opened to read the vegetation !Config map is to be given here. Usualy SECHIBA runs with a 5kmx5km !Config map which is derived from the IGBP one. We assume that we have !Config a classification in 87 types. This is Olson modified by Viovy. ! filename = '../surfmap/carteveg5km.nc' CALL getin('VEGETATION_FILE',filename) ! CALL flininfo(filename, iml, jml, lml, tml, fid) ! ! ALLOC_ERR=-1 ALLOCATE(lat_ful(iml), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of lat_ful : ",ALLOC_ERR STOP ENDIF ALLOC_ERR=-1 ALLOCATE(lon_ful(iml), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of lon_ful : ",ALLOC_ERR STOP ENDIF ALLOC_ERR=-1 ALLOCATE(vegmap(iml), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of vegmap : ",ALLOC_ERR STOP ENDIF ! WRITE(numout,*) 'Reading the OLSON type vegetation file' ! CALL flinget(fid, 'longitude', iml, jml, lml, tml, 1, 1, lon_ful) CALL flinget(fid, 'latitude', iml, jml, lml, tml, 1, 1, lat_ful) CALL flinget(fid, 'vegetation_map', iml, jml, lml, tml, 1, 1, vegmap) ! CALL flinclo(fid) ! IF (MAXVAL(vegmap) .LT. nolson) THEN WRITE(numout,*) 'WARNING -- WARNING' WRITE(numout,*) 'The vegetation map has to few vegetation types.' WRITE(numout,*) 'If you are lucky it will work but please check' ELSE IF ( MAXVAL(vegmap) .GT. nolson) THEN WRITE(numout,*) 'More vegetation types in file than the code can' WRITE(numout,*) 'deal with.: ', MAXVAL(vegmap), nolson STOP 'slowproc_interpol' ENDIF ! ! Some assumptions on the vegetation file. This information should be ! be computed or read from the file. ! It is the reolution in meters of the grid of the vegetation file. ! resol_lon = 5000. resol_lat = 5000. ! ! The number of maximum vegetation map points in the GCM grid should ! also be computed and not imposed here. nbvmax = iml/nbpt ! callsign="Vegetation map" ! ok_interpol = .FALSE. DO WHILE ( .NOT. ok_interpol ) WRITE(numout,*) "Projection arrays for ",callsign," : " WRITE(numout,*) "nbvmax = ",nbvmax ! ALLOC_ERR=-1 ALLOCATE(sub_index(nbpt, nbvmax), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of sub_index : ",ALLOC_ERR STOP ENDIF sub_index(:,:)=0 ALLOC_ERR=-1 ALLOCATE(sub_area(nbpt, nbvmax), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of sub_area : ",ALLOC_ERR STOP ENDIF sub_area(:,:)=zero ! CALL aggregate (nbpt, lalo, neighbours, resolution, contfrac, & & iml, lon_ful, lat_ful, resol_lon, resol_lat, callsign, & & nbvmax, sub_index, sub_area, ok_interpol) ! IF ( .NOT. ok_interpol ) THEN DEALLOCATE(sub_area) DEALLOCATE(sub_index) ! nbvmax = nbvmax * 2 ELSE ! DO ib = 1, nbpt idi=1 DO WHILE ( sub_area(ib,idi) > zero ) ip = sub_index(ib,idi) n_origveg(ib,NINT(vegmap(ip))) = n_origveg(ib,NINT(vegmap(ip))) + sub_area(ib,idi) n_found(ib) = n_found(ib) + sub_area(ib,idi) idi = idi +1 ENDDO ENDDO ! ENDIF ENDDO ! ! Now we know how many points of which Olson type from the fine grid fall ! into each box of the (coarse) model grid: n_origveg(nbpt,nolson) ! ! ! determine fraction of Olson vegetation type in each box of the coarse grid ! DO vid = 1, nolson WHERE ( n_found(:) .GT. 0 ) frac_origveg(:,vid) = n_origveg(:,vid) / n_found(:) ELSEWHERE frac_origveg(:,vid) = zero ENDWHERE ENDDO ! ! now finally calculate coarse vegetation map ! Find which model vegetation corresponds to each Olson type ! veget(:,:) = zero frac_nobio(:,:) = zero ! DO vid = 1, nolson ! DO jv = 1, nvm veget(:,jv) = veget(:,jv) + frac_origveg(:,vid) * vegcorr(vid,jv) ENDDO ! DO jv = 1, nnobio frac_nobio(:,jv) = frac_nobio(:,jv) + frac_origveg(:,vid) * nobiocorr(vid,jv) ENDDO ! ENDDO ! WRITE(numout,*) 'slowproc_interpol : Interpolation Done' ! ! Clean up the point of the map ! DO ib = 1, nbpt ! ! Let us see if all points found something in the 5km map ! ! IF ( n_found(ib) .EQ. 0 ) THEN ! ! Now we need to handle some exceptions ! IF ( lalo(ib,1) .LT. -56.0) THEN ! Antartica frac_nobio(ib,:) = zero frac_nobio(ib,iice) = un veget(ib,:) = zero ! ELSE IF ( lalo(ib,1) .GT. 70.0) THEN ! Artica frac_nobio(ib,:) = zero frac_nobio(ib,iice) = un veget(ib,:) = zero ! ELSE IF ( lalo(ib,1) .GT. 55.0 .AND. lalo(ib,2) .GT. -65.0 .AND. lalo(ib,2) .LT. -20.0) THEN ! Greenland frac_nobio(ib,:) = zero frac_nobio(ib,iice) = un veget(ib,:) = zero ! ELSE ! WRITE(numout,*) 'PROBLEM, no point in the 5km map found for this grid box',ib WRITE(numout,*) 'Longitude range : ', lalo(ib,2) WRITE(numout,*) 'Latitude range : ', lalo(ib,1) ! WRITE(numout,*) 'Looking for nearest point on the 5 km map' CALL slowproc_nearest (iml, lon_ful, lat_ful, & lalo(ib,2), lalo(ib,1), inear) WRITE(numout,*) 'Coordinates of the nearest point:', & lon_ful(inear),lat_ful(inear) ! DO jv = 1, nvm veget(ib,jv) = vegcorr(NINT(vegmap(inear)),jv) ENDDO ! DO jv = 1, nnobio frac_nobio(ib,jv) = nobiocorr(NINT(vegmap(inear)),jv) ENDDO ! ENDIF ! ENDIF ! ! ! Limit the smalest vegetation fraction to 0.5% ! DO vid = 1, nvm IF ( veget(ib,vid) .LT. min_vegfrac ) THEN veget(ib,vid) = zero ENDIF ENDDO ! sumf = SUM(frac_nobio(ib,:))+SUM(veget(ib,:)) frac_nobio(ib,:) = frac_nobio(ib,:)/sumf veget(ib,:) = veget(ib,:)/sumf ! ! ENDDO ! DEALLOCATE(vegmap) DEALLOCATE(lat_ful, lon_ful) DEALLOCATE(sub_index) DEALLOCATE(sub_area) ! RETURN ! END SUBROUTINE slowproc_interpol_NEW_g !! !! looks for nearest grid point on the fine map !! SUBROUTINE slowproc_nearest(iml, lon5, lat5, lonmod, latmod, inear) INTEGER(i_std), INTENT(in) :: iml REAL(r_std), DIMENSION(iml), INTENT(in) :: lon5, lat5 REAL(r_std), INTENT(in) :: lonmod, latmod INTEGER(i_std), INTENT(out) :: inear REAL(r_std) :: pi REAL(r_std) :: pa, p REAL(r_std) :: coscolat, sincolat REAL(r_std) :: cospa, sinpa REAL(r_std), ALLOCATABLE, DIMENSION(:) :: cosang INTEGER(i_std) :: i INTEGER(i_std), DIMENSION(1) :: ineartab INTEGER :: ALLOC_ERR ALLOC_ERR=-1 ALLOCATE(cosang(iml), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of cosang : ",ALLOC_ERR STOP ENDIF pi = 4.0 * ATAN(1.0) pa = pi/2.0 - latmod*pi/180.0 ! dist. entre pole n et point a cospa = COS(pa) sinpa = SIN(pa) DO i = 1, iml sincolat = SIN( pi/2.0 - lat5(i)*pi/180.0 ) coscolat = COS( pi/2.0 - lat5(i)*pi/180.0 ) p = (lonmod-lon5(i))*pi/180.0 ! angle entre a et b (leurs meridiens) ! dist(i) = ACOS( cospa*coscolat + sinpa*sincolat*COS(p)) cosang(i) = cospa*coscolat + sinpa*sincolat*COS(p) ENDDO ineartab = MAXLOC( cosang(:) ) inear = ineartab(1) DEALLOCATE(cosang) END SUBROUTINE slowproc_nearest !! !! Interpolate the Zobler soil type map !! SUBROUTINE slowproc_soilt(nbpt, lalo, neighbours, resolution, contfrac, soiltype, clayfraction) ! ! ! This subroutine should read the Zobler map and interpolate to the model grid. The method ! is to get fraction of the three main soiltypes for each grid box. ! The soil fraction are going to be put into the array soiltype in the following order : ! coarse, medium and fine. ! ! ! 0.1 INPUT ! INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longitudes (beware of the order !) INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) ! Vector of neighbours for each grid point ! (1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW) REAL(r_std), INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid-box in X and Y REAL(r_std), INTENT(in) :: contfrac(nbpt) ! Fraction of land in each grid box. ! ! 0.2 OUTPUT ! REAL(r_std), INTENT(out) :: soiltype(nbpt, nstm) ! Soil type map to be created from the Zobler map REAL(r_std), INTENT(out) :: clayfraction(nbpt) ! The fraction of clay as used by STOMATE ! ! ! 0.3 LOCAL ! INTEGER(i_std) :: nbvmax ! CHARACTER(LEN=80) :: filename INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, ip, jp, fopt, ilf, nbexp REAL(r_std) :: lev(1), date, dt INTEGER(i_std) :: itau(1) REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lat_rel, lon_rel, soiltext INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: mask REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: sub_area INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:) :: sub_index INTEGER(i_std), ALLOCATABLE, DIMENSION(:) :: solt REAL(r_std) :: sgn CHARACTER(LEN=30) :: callsign ! ! Number of texture classes in Zobler ! INTEGER(i_std), PARAMETER :: classnb = 7 REAL(r_std) :: textfrac_table(classnb, nstm) ! LOGICAL :: ok_interpol ! optionnal return of aggregate_2d ! INTEGER :: ALLOC_ERR ! ! CALL get_soilcorr (classnb, textfrac_table) ! ! Needs to be a configurable variable ! ! !Config Key = SOILTYPE_FILE !Config Desc = Name of file from which soil types are read !Config Def = ../surfmap/soils_param.nc !Config If = !IMPOSE_VEG !Config Help = The name of the file to be opened to read the soil types. !Config The data from this file is then interpolated to the grid of !Config of the model. The aim is to get fractions for sand loam and !Config clay in each grid box. This information is used for soil hydrology !Config and respiration. ! filename = '../surfmap/soils_param.nc' CALL getin_p('SOILTYPE_FILE',filename) ! IF (is_root_prc) THEN CALL flininfo(filename,iml, jml, lml, tml, fid) CALL flinclo(fid) ENDIF CALL bcast(iml) CALL bcast(jml) CALL bcast(lml) CALL bcast(tml) ! ! soils_param.nc file is 1° soit texture file. ! ALLOC_ERR=-1 ALLOCATE(lat_rel(iml,jml), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of lat_rel : ",ALLOC_ERR STOP ENDIF ALLOC_ERR=-1 ALLOCATE(lon_rel(iml,jml), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of lon_rel : ",ALLOC_ERR STOP ENDIF ALLOC_ERR=-1 ALLOCATE(mask(iml,jml), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of mask : ",ALLOC_ERR STOP ENDIF ALLOC_ERR=-1 ALLOCATE(soiltext(iml,jml), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of soiltext : ",ALLOC_ERR STOP ENDIF ! IF (is_root_prc) CALL flinopen(filename, .FALSE., iml, jml, lml, lon_rel, lat_rel, lev, tml, itau, date, dt, fid) CALL bcast(lon_rel) CALL bcast(lat_rel) CALL bcast(itau) CALL bcast(date) CALL bcast(dt) ! IF (is_root_prc) CALL flinget(fid, 'soiltext', iml, jml, lml, tml, 1, 1, soiltext) CALL bcast(soiltext) ! IF (is_root_prc) CALL flinclo(fid) ! nbexp = 0 ! ! ! Mask of permitted variables. ! mask(:,:) = zero DO ip=1,iml DO jp=1,jml IF (soiltext(ip,jp) .GT. min_sechiba) THEN mask(ip,jp) = un ENDIF ENDDO ENDDO ! nbvmax = 200 ! callsign = "Soil types" ! ok_interpol = .FALSE. DO WHILE ( .NOT. ok_interpol ) WRITE(numout,*) "Projection arrays for ",callsign," : " WRITE(numout,*) "nbvmax = ",nbvmax ! ALLOC_ERR=-1 ALLOCATE(solt(nbvmax), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of solt : ",ALLOC_ERR STOP ENDIF ALLOC_ERR=-1 ALLOCATE(sub_index(nbpt,nbvmax,2), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of sub_index : ",ALLOC_ERR STOP ENDIF sub_index(:,:,:)=0 ALLOC_ERR=-1 ALLOCATE(sub_area(nbpt,nbvmax), STAT=ALLOC_ERR) IF (ALLOC_ERR/=0) THEN WRITE(numout,*) "ERROR IN ALLOCATION of sub_area : ",ALLOC_ERR STOP ENDIF sub_area(:,:)=zero ! CALL aggregate_p(nbpt, lalo, neighbours, resolution, contfrac, & & iml, jml, lon_rel, lat_rel, mask, callsign, & & nbvmax, sub_index, sub_area, ok_interpol) ! IF ( .NOT. ok_interpol ) THEN DEALLOCATE(sub_area) DEALLOCATE(sub_index) DEALLOCATE(solt) ! nbvmax = nbvmax * 2 ENDIF ENDDO ! DO ib =1, nbpt ! soiltype(ib,:) = zero clayfraction(ib) = zero ! ! GO through the point we have found ! ! fopt = COUNT(sub_area(ib,:) > zero) ! ! Check that we found some points ! IF ( fopt .EQ. 0) THEN nbexp = nbexp + 1 soiltype(ib,:) = soiltype_default(:) clayfraction(ib) = clayfraction_default ELSE ! DO ilf = 1,fopt solt(ilf) = soiltext(sub_index(ib,ilf,1),sub_index(ib,ilf,2)) ENDDO ! sgn = zero ! ! Compute the average bare soil albedo parameters ! DO ilf = 1,fopt ! ! We have to take care of two exceptions here : type 6 = glacier and type 0 = ocean ! IF ( (solt(ilf) .LE. classnb) .AND. (solt(ilf) .GT. 0) .AND.& & (solt(ilf) .NE. 6)) THEN SELECTCASE(solt(ilf)) CASE(1) soiltype(ib,1) = soiltype(ib,1) + sub_area(ib,ilf) CASE(2) soiltype(ib,2) = soiltype(ib,2) + sub_area(ib,ilf) CASE(3) soiltype(ib,2) = soiltype(ib,2) + sub_area(ib,ilf) CASE(4) soiltype(ib,2) = soiltype(ib,2) + sub_area(ib,ilf) CASE(5) soiltype(ib,3) = soiltype(ib,3) + sub_area(ib,ilf) CASE(7) soiltype(ib,2) = soiltype(ib,2) + sub_area(ib,ilf) CASE DEFAULT WRITE(numout,*) 'We should not be here, an impossible case appeared' STOP 'slowproc_soilt' END SELECT clayfraction(ib) = clayfraction(ib) + & & textfrac_table(solt(ilf),3) * sub_area(ib,ilf) sgn = sgn + sub_area(ib,ilf) ELSE IF (solt(ilf) .GT. classnb) THEN WRITE(numout,*) 'The file contains a soil color class which is incompatible with this program' STOP 'slowproc_soilt' ENDIF ENDIF ! ENDDO ! ! Normalize the surface ! IF ( sgn .LT. min_sechiba) THEN nbexp = nbexp + 1 soiltype(ib,:) = soiltype_default(:) clayfraction(ib) = clayfraction_default ELSE soiltype(ib,:) = soiltype(ib,:)/sgn clayfraction(ib) = clayfraction(ib)/sgn ENDIF ! ENDIF ! ENDDO ! IF ( nbexp .GT. 0 ) THEN WRITE(numout,*) 'slowproc_soilt : The interpolation of the bare soil albedo had ', nbexp WRITE(numout,*) 'slowproc_soilt : points without data. This are either coastal points or' WRITE(numout,*) 'slowproc_soilt : ice covered land.' WRITE(numout,*) 'slowproc_soilt : The problem was solved by using the default soil types.' ENDIF ! DEALLOCATE (lat_rel) DEALLOCATE (lon_rel) DEALLOCATE (mask) DEALLOCATE (sub_area) DEALLOCATE (sub_index) DEALLOCATE (soiltext) DEALLOCATE (solt) ! ! RETURN ! END SUBROUTINE slowproc_soilt ! END MODULE slowproc