MODULE routing_native_flow_mod USE ioipsl USE xios_orchidee USE ioipsl_para USE constantes USE constantes_soil USE pft_parameters USE sechiba_io_p USE interpol_help USE grid USE mod_orchidee_para IMPLICIT NONE PRIVATE PUBLIC :: routing_flow_xios_initialize, routing_flow_set, routing_flow_get, routing_flow_main PUBLIC :: routing_flow_initialize, routing_flow_finalize, routing_flow_clear, routing_flow_make_mean PUBLIC :: compute_coastline INTEGER,SAVE :: nbpt !$OMP THREADPRIVATE(nbpt) !! PARAMETERS REAL(r_std), SAVE :: fast_tcst = 3.0 !! Property of the fast reservoir (day/m) !$OMP THREADPRIVATE(fast_tcst) REAL(r_std), SAVE :: slow_tcst = 25.0 !! Property of the slow reservoir (day/m) !$OMP THREADPRIVATE(slow_tcst) REAL(r_std), SAVE :: stream_tcst = 0.24 !! Property of the stream reservoir (day/m) !$OMP THREADPRIVATE(stream_tcst) REAL(r_std),SAVE,ALLOCATABLE :: runoff_mean(:) !$OMP THREADPRIVATE(runoff_mean) REAL(r_std),SAVE,ALLOCATABLE :: drainage_mean(:) !$OMP THREADPRIVATE(drainage_mean) REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED,TARGET :: fast_reservoir_r(:) !! Water amount in the fast reservoir (kg) - (local routing grid) !$OMP THREADPRIVATE(fast_reservoir_r) REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED,TARGET:: slow_reservoir_r(:) !! Water amount in the slow reservoir (kg) - (local routing grid) !$OMP THREADPRIVATE(slow_reservoir_r) REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED,TARGET :: stream_reservoir_r(:) !! Water amount in the stream reservoir (kg) - (local routing grid) !$OMP THREADPRIVATE(stream_reservoir_r) REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED :: riverflow_mean(:) !! Water amount in the stream reservoir (kg) - (local routing grid) !$OMP THREADPRIVATE(riverflow_mean) REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED :: coastalflow_mean(:) !! Water amount in the stream reservoir (kg) - (local routing grid) !$OMP THREADPRIVATE(coastalflow_mean) REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED :: lakeinflow_mean(:) !! Water amount in the stream reservoir (kg) - (local routing grid) !$OMP THREADPRIVATE(lakeinflow_mean) ! ! Specific variables for simple routing ! REAL(r_std),SAVE,ALLOCATABLE :: topoind_r(:) !! Topographic index of the retention time (m) index - (local routing grid) !$OMP THREADPRIVATE(topoind_r) INTEGER,SAVE,ALLOCATABLE :: route_flow_rp1(:) !! flow index from cell to neighboring cell following the trip direction - (local routing grid + halo) !$OMP THREADPRIVATE(route_flow_rp1) LOGICAL,SAVE,ALLOCATABLE :: is_lakeinflow_r(:) !! is lake inflow point - (local routing grid) !$OMP THREADPRIVATE(is_lakeinflow_r) LOGICAL,SAVE,ALLOCATABLE :: is_coastalflow_r(:) !! is coastal flow point - (local routing grid) !$OMP THREADPRIVATE(is_coastalflow_r) LOGICAL,SAVE,ALLOCATABLE :: is_riverflow_r(:) !! is river flow point - (local routing grid) !$OMP THREADPRIVATE(is_riverflow_r) LOGICAL,SAVE,ALLOCATABLE :: is_coastline(:) !! is coastline point - (local native grid) !$OMP THREADPRIVATE(is_coastline) LOGICAL,SAVE,ALLOCATABLE :: is_streamflow_r(:) !! is stream flow point - (local routing grid) !$OMP THREADPRIVATE(is_streamflow_r) LOGICAL,SAVE,ALLOCATABLE,PUBLIC,PROTECTED :: routing_mask_r(:) !! valid routing point - (local routing grid) !$OMP THREADPRIVATE(routing_mask_r) LOGICAL,SAVE,ALLOCATABLE :: coast_mask(:) !! is a coast point - (local native grid) !$OMP THREADPRIVATE(coast_mask) INTEGER,SAVE :: total_coast_points !! global number of coast point - (local native grid) !$OMP THREADPRIVATE(total_coast_points) INTEGER,SAVE :: nbpt_r !! number of point in local routing grid !$OMP THREADPRIVATE(nbpt_r) INTEGER,SAVE :: nbpt_rp1 !! number of point in local routing grid with halo of 1 !$OMP THREADPRIVATE(nbpt_rp1) REAL(r_std),SAVE,ALLOCATABLE :: routing_weight(:) !! Weight to transform runoff and drainage flux to water quantity in a conservative way (local native grid -> local routing grid) !$OMP THREADPRIVATE(routing_weight) REAL(r_std),SAVE,ALLOCATABLE :: routing_weight_in(:) !! Weight to transform runoff and drainage flux to water quantity in a conservative way (local native grid -> local routing grid) !$OMP THREADPRIVATE(routing_weight_in) REAL(r_std),SAVE,ALLOCATABLE :: unrouted_weight(:) !! Weight to transform runoff and drainage flux to water quantity in a conservative way (local native grid -> local routing grid) !$OMP THREADPRIVATE(unrouted_weight) REAL(r_std),SAVE,ALLOCATABLE :: weight_coast_to_coast_r(:) !$OMP THREADPRIVATE(weight_coast_to_coast_r) REAL(r_std),SAVE,ALLOCATABLE :: weight_coast_to_lake_r(:) !$OMP THREADPRIVATE(weight_coast_to_lake_r) REAL(r_std),SAVE,ALLOCATABLE :: weight_lake_to_coast_r(:) !$OMP THREADPRIVATE(weight_lake_to_coast_r) REAL(r_std),SAVE,ALLOCATABLE :: weight_lake_to_lake_r(:) !$OMP THREADPRIVATE(weight_lake_to_lake_r) REAL(r_std),SAVE,ALLOCATABLE :: basins_area_r(:) !$OMP THREADPRIVATE(basins_area_r) !! when doing interpolation from local routing grid to local native grid (river flow+coastal flow) REAL(r_std),SAVE,ALLOCATABLE :: basins_extended_r(:) !! basins riverflow id (local routing grid) !$OMP THREADPRIVATE(basins_extended_r) INTEGER(i_std) :: basins_count !! number of basins (local routing grid) !$OMP THREADPRIVATE(basins_count) INTEGER(i_std) :: basins_out !! number of basins to output for diag !$OMP THREADPRIVATE(basins_out) INTEGER(i_std) :: split_routing !! time spliting for routing !$OMP THREADPRIVATE(split_routing) INTEGER :: nb_station CHARACTER(LEN=60),ALLOCATABLE :: station(:) REAL,ALLOCATABLE :: station_lonlat(:,:) INTEGER,ALLOCATABLE :: station_index(:) INTEGER :: station_ts = 0 ! INTEGER(i_std), PARAMETER :: nb_stations=14 ! REAL(r_std),PARAMETER :: station_lon(nb_stations) = & ! (/ -162.8830, -90.9058, -55.5110, -49.3242, -133.7447, -63.6000, 28.7167, & ! 15.3000, 66.5300, 89.6700, 86.5000, 127.6500, 3.3833, 117.6200 /) ! REAL(r_std),PARAMETER :: station_lat(nb_stations) = & ! (/ 61.9340, 32.3150, -1.9470, -5.1281, 67.4583, 8.1500, 45.2167, & ! -4.3000, 66.5700, 25.1800, 67.4800, 70.7000, 11.8667, 30.7700 /) ! CHARACTER(LEN=17),PARAMETER :: station_name(nb_stations) = & ! (/ "Pilot station ", "Vicksburg ", "Obidos ", & ! "Itupiranga ", "Arctic red river ", "Puente Angostura ", & ! "Ceatal Izmail ", "Kinshasa ", "Salekhard ", & ! "Bahadurabad ", "Igarka ", "Kusur ", & ! "Malanville ", "Datong " /) CONTAINS SUBROUTINE routing_flow_get(slow_reservoir_r, fast_reservoir_r, stream_reservoir_r, riverflow_mean, coastalflow_mean, & lakeinflow_mean) IMPLICIT NONE REAL(r_std),OPTIONAL, INTENT(OUT) :: slow_reservoir_r(:) REAL(r_std),OPTIONAL, INTENT(OUT) :: fast_reservoir_r(:) REAL(r_std),OPTIONAL, INTENT(OUT) :: stream_reservoir_r(:) REAL(r_std),OPTIONAL, INTENT(OUT) :: riverflow_mean(:) REAL(r_std),OPTIONAL, INTENT(OUT) :: coastalflow_mean(:) REAL(r_std),OPTIONAL, INTENT(OUT) :: lakeinflow_mean(:) CALL routing_flow_get_(slow_reservoir_r, fast_reservoir_r, stream_reservoir_r, riverflow_mean, coastalflow_mean, lakeinflow_mean) END SUBROUTINE routing_flow_get SUBROUTINE routing_flow_get_(slow_reservoir_r_, fast_reservoir_r_, stream_reservoir_r_, riverflow_mean_, & coastalflow_mean_, lakeinflow_mean_) IMPLICIT NONE REAL(r_std),OPTIONAL, INTENT(OUT) :: slow_reservoir_r_(:) REAL(r_std),OPTIONAL, INTENT(OUT) :: fast_reservoir_r_(:) REAL(r_std),OPTIONAL, INTENT(OUT) :: stream_reservoir_r_(:) REAL(r_std),OPTIONAL, INTENT(OUT) :: riverflow_mean_(:) REAL(r_std),OPTIONAL, INTENT(OUT) :: coastalflow_mean_(:) REAL(r_std),OPTIONAL, INTENT(OUT) :: lakeinflow_mean_(:) IF (PRESENT(slow_reservoir_r_)) slow_reservoir_r_ = slow_reservoir_r IF (PRESENT(fast_reservoir_r_)) fast_reservoir_r_ = fast_reservoir_r IF (PRESENT(stream_reservoir_r_)) stream_reservoir_r_ = stream_reservoir_r IF (PRESENT(riverflow_mean_)) riverflow_mean_ = riverflow_mean IF (PRESENT(coastalflow_mean_)) coastalflow_mean_ = coastalflow_mean IF (PRESENT(lakeinflow_mean_)) lakeinflow_mean_ = lakeinflow_mean END SUBROUTINE routing_flow_get_ SUBROUTINE routing_flow_set(slow_reservoir_r, fast_reservoir_r, stream_reservoir_r, riverflow_mean, coastalflow_mean, & lakeinflow_mean) IMPLICIT NONE REAL(r_std),OPTIONAL, INTENT(IN) :: slow_reservoir_r(:) REAL(r_std),OPTIONAL, INTENT(IN) :: fast_reservoir_r(:) REAL(r_std),OPTIONAL, INTENT(IN) :: stream_reservoir_r(:) REAL(r_std),OPTIONAL, INTENT(IN) :: riverflow_mean(:) REAL(r_std),OPTIONAL, INTENT(IN) :: coastalflow_mean(:) REAL(r_std),OPTIONAL, INTENT(IN) :: lakeinflow_mean(:) CALL routing_flow_set_(slow_reservoir_r, fast_reservoir_r, stream_reservoir_r, riverflow_mean, coastalflow_mean, lakeinflow_mean) END SUBROUTINE routing_flow_set SUBROUTINE routing_flow_set_(slow_reservoir_r_, fast_reservoir_r_, stream_reservoir_r_, riverflow_mean_, & coastalflow_mean_, lakeinflow_mean_) IMPLICIT NONE REAL(r_std),OPTIONAL, INTENT(IN) :: slow_reservoir_r_(:) REAL(r_std),OPTIONAL, INTENT(IN) :: fast_reservoir_r_(:) REAL(r_std),OPTIONAL, INTENT(IN) :: stream_reservoir_r_(:) REAL(r_std),OPTIONAL, INTENT(IN) :: riverflow_mean_(:) REAL(r_std),OPTIONAL, INTENT(IN) :: coastalflow_mean_(:) REAL(r_std),OPTIONAL, INTENT(IN) :: lakeinflow_mean_(:) IF (PRESENT(slow_reservoir_r_)) slow_reservoir_r = slow_reservoir_r_ IF (PRESENT(fast_reservoir_r_)) fast_reservoir_r = fast_reservoir_r_ IF (PRESENT(stream_reservoir_r_)) stream_reservoir_r = stream_reservoir_r_ IF (PRESENT(riverflow_mean_)) riverflow_mean = riverflow_mean_ IF (PRESENT(coastalflow_mean_)) coastalflow_mean = coastalflow_mean_ IF (PRESENT(lakeinflow_mean_)) lakeinflow_mean = lakeinflow_mean_ END SUBROUTINE routing_flow_set_ !! ============================================================================================================================= !! SUBROUTINE: routing_simple_xios_initialize !! !>\BRIEF Initialize xios dependant defintion before closing context defintion !! !! DESCRIPTION: Initialize xios dependant defintion before closing context defintion. !! This subroutine is called before the xios context is closed. !! !! RECENT CHANGE(S): None !! !! REFERENCE(S): None !! !! FLOWCHART: None !! \n !_ ============================================================================================================================== SUBROUTINE routing_flow_xios_initialize USE xios USE routing, ONLY : routing_names IMPLICIT NONE INTEGER(i_std) ::ib INTEGER :: nbasmax=1 !! 0 Variable and parameter description CHARACTER(LEN=60),ALLOCATABLE :: label(:) LOGICAL :: file_exists IF (is_omp_root) THEN CALL xios_get_axis_attr("basins", n_glo=basins_out) ! get nb basins to output ALLOCATE(label(basins_out)) CALL routing_names(basins_out,label) CALL xios_set_axis_attr("basins", label=label) ! set riverflow basins name INQUIRE(FILE="routing_start.nc", EXIST=file_exists) IF (file_exists) CALL xios_set_file_attr("routing_start", enabled=.TRUE.) ENDIF !! Define XIOS axis size needed for the model output ! Add axis for homogeneity between all routing schemes, these dimensions are currently not used in this scheme CALL xios_orchidee_addaxis("nbhtu", nbasmax, (/(REAL(ib,r_std),ib=1,nbasmax)/)) CALL xios_orchidee_addaxis("nbasmon", 1, (/(REAL(ib,r_std),ib=1,1)/)) END SUBROUTINE routing_flow_xios_initialize SUBROUTINE routing_flow_initialize(kjit, rest_id, nbpt_, dt_routing, contfrac, nbpt_r_, riverflow, coastalflow) IMPLICIT NONE INTEGER ,INTENT(IN) :: kjit INTEGER ,INTENT(IN) :: rest_id INTEGER, INTENT(IN) :: nbpt_ !! nb points orchidee grid REAL(r_std), INTENT(IN) :: dt_routing REAL(r_std), INTENT(IN) :: contfrac(nbpt_) !! fraction of land INTEGER, INTENT(OUT) :: nbpt_r_ !! nb points routing native grid REAL(r_std), INTENT(OUT) :: riverflow(nbpt_) !! Outflow of the major rivers. The flux will be located on the continental grid but this should be a coastal point (kg/dt) REAL(r_std), INTENT(OUT) :: coastalflow(nbpt_) !! Outflow on coastal points by small basins. This is the water which flows in a disperse way into the ocean (kg/dt) INTEGER :: ier CHARACTER(LEN=80) :: var_name !! To store variables names for I/O (unitless) nbpt = nbpt_ CALL routing_flow_init_local(contfrac, nbpt_r_) nbpt_r = nbpt_r_ CALL routing_flow_init_mean(kjit, rest_id) CALL initialize_stations(dt_routing) ! ! Put into the restart file the fluxes so that they can be regenerated at restart. ! ALLOCATE (lakeinflow_mean(nbpt), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for lakeinflow_mean','','') var_name = 'lakeinflow' CALL ioconf_setatt_p('UNITS', 'Kg/dt') CALL ioconf_setatt_p('LONG_NAME','Lake inflow') CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., lakeinflow_mean, "gather", nbp_glo, index_g) CALL setvar_p (lakeinflow_mean, val_exp, 'NO_KEYWORD', zero) ALLOCATE (riverflow_mean(nbpt), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for riverflow_mean','','') var_name = 'riverflow' CALL ioconf_setatt_p('UNITS', 'Kg/dt') CALL ioconf_setatt_p('LONG_NAME','River flux into the sea') CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., riverflow_mean, "gather", nbp_glo, index_g) CALL setvar_p (riverflow_mean, val_exp, 'NO_KEYWORD', zero) ALLOCATE (coastalflow_mean(nbpt), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for coastalflow_mean','','') var_name = 'coastalflow' CALL ioconf_setatt_p('UNITS', 'Kg/dt') CALL ioconf_setatt_p('LONG_NAME','Diffuse flux into the sea') CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., coastalflow_mean, "gather", nbp_glo, index_g) CALL setvar_p (coastalflow_mean, val_exp, 'NO_KEYWORD', zero) riverflow(:) = riverflow_mean(:) coastalflow(:) = coastalflow_mean(:) END SUBROUTINE routing_flow_initialize SUBROUTINE compute_coastline(contfrac, coastline) USE mod_orchidee_para USE grid IMPLICIT NONE REAL(r_std),INTENT(IN) :: contfrac(nbpt) ! INPUT : fraction of continent (unitless) LOGICAL,INTENT(OUT) :: coastline(nbpt) ! OUTPUT : coastline mask (true : on coastaline, else false REAL(r_std) :: contfrac_glo(nbp_glo) REAL(r_std) :: contfrac2D_glo(iim_g*jjm_g) LOGICAL :: coastline_glo(nbp_glo) LOGICAL :: coastline2D_glo(iim_g*jjm_g) INTEGER(i_std) :: i,j,ij,next_i,next_j,next_ij,m,k REAL(r_std), PARAMETER :: epsilon=1e-5 LOGICAL :: is_periodic is_periodic = global CALL gather(contfrac,contfrac_glo) contfrac2D_glo(:)=0 DO i=1,nbp_glo contfrac2D_glo(index_g(i)) = contfrac_glo(i) ENDDO IF (is_mpi_root .AND. is_omp_root) THEN coastline2D_glo(:)=.FALSE. DO j=1,jjm_g DO i=1,iim_g ij=(j-1)*iim_g+i IF (contfrac2D_glo(ij)>epsilon .AND. contfrac2D_glo(ij)<1-epsilon) THEN coastline2D_glo(ij)=.TRUE. ELSE IF (contfrac2D_glo(ij)>=1-epsilon .AND. grid_type/=unstructured ) THEN DO k=-1,1 DO m=-1,1 IF (k==0 .AND. m==0) CYCLE next_i=i+k next_j=j+m IF (next_i==0) THEN ! manage periodicity IF (is_periodic) THEN next_i=iim_g ELSE coastline2D_glo(ij)=.TRUE. CYCLE ENDIF ENDIF IF (next_i==iim_g+1) THEN! manage periodicity IF (is_periodic) THEN next_i=1 !! periodic but not true for limited area ELSE coastline2D_glo(ij)=.TRUE. CYCLE ENDIF ENDIF IF (next_j==0 .OR. next_j==jjm_g+1) THEN IF (.NOT. is_periodic) coastline2D_glo(ij)=.TRUE. CYCLE ENDIF next_ij = (next_j-1)*iim_g+next_i IF (contfrac2D_glo(next_ij)<=epsilon) coastline2D_glo(ij)=.TRUE. ENDDO ENDDO ENDIF ENDDO ENDDO DO i=1,nbp_glo coastline_glo(i) = coastline2D_glo(index_g(i)) ENDDO ENDIF CALL scatter(coastline_glo, coastline) END SUBROUTINE compute_coastline SUBROUTINE routing_flow_correct_riverflow(ni, nj, contfrac, coastline, trip_r, trip_extended_r, topoind_r) USE xios USE grid, ONLY : global IMPLICIT NONE INCLUDE "mpif.h" INTEGER, INTENT(IN) :: ni ! INPUT : size i (longitude) of the local routing domain INTEGER, INTENT(IN) :: nj ! INPUT : size i (latitude) of the local routing domain REAL(r_std), INTENT(IN) :: contfrac(nbp_mpi) ! INPUT : continental fraction (unittless) LOGICAL, INTENT(IN) :: coastline(nbp_mpi) ! INPUT/OUTPUT : coastline mask REAL(r_std), INTENT(INOUT) :: trip_r(ni,nj) ! INPUT/OUTPUT : diection of flow, which will be modified by the routine REAL(r_std), INTENT(INOUT) :: trip_extended_r(ni,nj) ! INPUT : direction of flow extended to ocean points REAL(r_std), INTENT(INOUT) :: topoind_r(ni,nj) ! INPUT/OUTPUT : topographic index which will be modified by the routine if extended to ocean REAL(r_std) :: contfrac_r(ni,nj) REAL(r_std) :: mask_coast(nbp_mpi) REAL(r_std) :: frac_coast_r(ni,nj) REAL(r_std) :: trip_extended_rp1(0:ni+1,0:nj+1) REAL(r_std) :: trip_rp1(0:ni+1,0:nj+1) REAL(r_std) :: state_r(ni,nj) REAL(r_std) :: state_rp1(0:ni+1,0:nj+1) INTEGER, PARAMETER :: is_ter=0 INTEGER, PARAMETER :: is_coast=1 INTEGER, PARAMETER :: is_oce=2 INTEGER :: updated INTEGER :: it, i,j,k,m, nexti,nextj,previ,prevj,ierr INTEGER :: ibegin, jbegin, ni_glo, nj_glo REAL(r_std) :: lonvalue_1d(ni), latvalue_1d(nj) REAL(r_std),PARAMETER :: epsilon=1e-3 TYPE(xios_duration) :: ts TYPE(xios_domain) :: domain_hdl TYPE(xios_domaingroup) :: domain_def_hdl TYPE(xios_field) :: field_hdl TYPE(xios_fieldgroup) :: field_def_hdl TYPE(xios_file) :: file_hdl TYPE(xios_filegroup) :: file_def_hdl TYPE(xios_expand_domain) :: domain_expand_hdl TYPE(xios_date) :: start_date, time_origin CHARACTER(LEN=20) :: calendar_type LOGICAL :: true=.TRUE. CALL xios_get_domain_attr("routing_domain", ibegin=ibegin, jbegin=jbegin, ni_glo=ni_glo,nj_glo=nj_glo) ! get routing domain dimension CALL xios_get_domain_attr("routing_domain", lonvalue_1d=lonvalue_1d, latvalue_1d=latvalue_1d) ! get routing domain dimension ! manage non periodic boundary in case of LAM IF (.NOT. global) THEN IF (ibegin==0) THEN i=1 ; DO j=1,nj IF (trip_r(i,j)>=1 .AND. trip_r(i,j)<=8) trip_r(i,j) = 97 IF (trip_extended_r(i,j)>=1 .AND. trip_extended_r(i,j)<=8) trip_extended_r(i,j) = 97 ENDDO ENDIF IF (ibegin+ni==ni_glo) THEN i=ni ; DO j=1,nj IF (trip_r(i,j)>=1 .AND. trip_r(i,j)<=8) trip_r(i,j) = 97 IF (trip_extended_r(i,j)>=1 .AND. trip_extended_r(i,j)<=8) trip_extended_r(i,j) = 97 ENDDO ENDIF IF (jbegin==0) THEN j=1 ; DO i=1,ni IF (trip_r(i,j)>=1 .AND. trip_r(i,j)<=8) trip_r(i,j) = 97 IF (trip_extended_r(i,j)>=1 .AND. trip_extended_r(i,j)<=8) trip_extended_r(i,j) = 97 ENDDO ENDIF IF (jbegin+nj==nj_glo) THEN j=nj ; DO i=1,ni IF (trip_r(i,j)>=1 .AND. trip_r(i,j)<=8) trip_r(i,j) = 97 IF (trip_extended_r(i,j)>=1 .AND. trip_extended_r(i,j)<=8) trip_extended_r(i,j) = 97 ENDDO ENDIF ENDIF CALL xios_send_field("routing_contfrac",contfrac) CALL xios_recv_field("routing_contfrac_r", contfrac_r) CALL xios_send_field("trip_ext_r",trip_extended_r) CALL xios_recv_field("trip_ext_rp1",trip_extended_rp1) mask_coast=0 WHERE(coastline) mask_coast=1 CALL xios_send_field("mask_coastline",mask_coast) CALL xios_recv_field("frac_coastline_r",frac_coast_r) DO j=1,nj DO i=1,ni IF (contfrac_r(i,j) 2 ) THEN state_r(i,j)=is_oce ELSE IF (contfrac_r(i,j)>1-epsilon) THEN state_r(i,j)=is_ter ELSE state_r(i,j)=is_coast ENDIF IF (frac_coast_r(i,j)>epsilon .AND. frac_coast_r(i,j)<=1+epsilon ) state_r(i,j)=is_coast ENDDO ENDDO DO j=1,nj DO i=1,ni IF (trip_r(i,j)>=100) THEN IF (state_r(i,j) /= is_oce) trip_r(i,j)=trip_extended_r(i,j) ENDIF ENDDO ENDDO ! create a new XIOS context to solve iterative process of correction ! so we can use filter to update halo ! CALL xios_get_calendar_type(calendar_type) ! CALL xios_get_start_date(start_date) ! CALL xios_get_time_origin(time_origin) ! ! CALL xios_context_initialize("orchidee_routing_correction",MPI_COMM_ORCH) ! CALL xios_orchidee_change_context("orchidee_routing_correction") ! ts.second=1 ! calendar_type="gregorian" ! CALL xios_define_calendar(type=calendar_type,start_date=start_date,time_origin=time_origin, timestep=ts ) ! CALL xios_get_handle("domain_definition",domain_def_hdl) ! CALL xios_add_child(domain_def_hdl, domain_hdl, "routing_domain") ! CALL xios_set_attr(domain_hdl, type="rectilinear", ibegin=ibegin, jbegin=jbegin, ni=ni, nj=nj, ni_glo=ni_glo, nj_glo=nj_glo) ! set routing domain dimension ! CALL xios_set_attr(domain_hdl, lonvalue_1d=lonvalue_1d, latvalue_1d=latvalue_1d) ! set routing domain dimension ! CALL xios_add_child(domain_def_hdl, domain_hdl, "routing_domain_expand") ! CALL xios_set_attr(domain_hdl, domain_ref="routing_domain") ! CALL xios_add_child(domain_hdl,domain_expand_hdl) ! CALL xios_set_attr(domain_expand_hdl, type="edge", i_periodic=.TRUE., j_periodic=.TRUE.) ! ! CALL xios_get_handle("field_definition",field_def_hdl) ! CALL xios_add_child(field_def_hdl,field_hdl,"trip_r_init") ! CALL xios_set_attr(field_hdl,domain_ref="routing_domain") ! CALL xios_add_child(field_def_hdl,field_hdl,"trip_r_final") ! CALL xios_set_attr(field_hdl,domain_ref="routing_domain") ! CALL xios_add_child(field_def_hdl,field_hdl,"trip_r") ! CALL xios_set_attr(field_hdl,domain_ref="routing_domain") ! CALL xios_add_child(field_def_hdl,field_hdl,"trip_rp1") ! CALL xios_set_attr(field_hdl, field_ref="trip_r", domain_ref="routing_domain_expand", read_access=.TRUE.) ! CALL xios_add_child(field_def_hdl,field_hdl,"state_r") ! CALL xios_set_attr(field_hdl,domain_ref="routing_domain") ! CALL xios_add_child(field_def_hdl,field_hdl,"state_rp1") ! CALL xios_set_attr(field_hdl, field_ref="state_r", domain_ref="routing_domain_expand", read_access=.TRUE.) ! ! CALL xios_get_handle("file_definition",file_def_hdl) ! CALL xios_add_child(file_def_hdl,file_hdl,"routing_correction") ! CALL xios_set_attr(file_hdl,type="one_file", output_freq=ts, sync_freq=ts, enabled=.TRUE.) ! CALL xios_add_child(file_hdl, field_hdl) ! CALL xios_set_attr(field_hdl, field_ref="trip_r", operation="instant") ! CALL xios_add_child(file_hdl, field_hdl) ! CALL xios_set_attr(field_hdl, field_ref="state_r", operation="instant") ! CALL xios_add_child(file_hdl, field_hdl) ! CALL xios_set_attr(field_hdl, field_ref="trip_r_init", operation="once") ! CALL xios_add_child(file_hdl, field_hdl) ! CALL xios_set_attr(field_hdl, field_ref="state_r_init", operation="once") CALL xios_context_initialize("orchidee_init_routing",MPI_COMM_ORCH) CALL xios_orchidee_change_context("orchidee_init_routing") CALL xios_set_domain_attr("routing_domain", ibegin=ibegin, jbegin=jbegin, ni=ni, nj=nj, ni_glo=ni_glo, nj_glo=nj_glo) ! get routing domain dimension CALL xios_set_domain_attr("routing_domain", lonvalue_1d=lonvalue_1d, latvalue_1d=latvalue_1d) ! get routing domain dimension CALL xios_close_context_definition() CALL xios_update_calendar(1) CALL xios_send_field("trip_r_init",trip_r) CALL xios_send_field("trip_r",trip_r) CALL xios_recv_field("trip_rp1",trip_rp1) CALL xios_send_field("state_r",state_r) CALL xios_recv_field("state_rp1",state_rp1) DO j=1,nj DO i=1,ni IF (trip_rp1(i,j)<=8) THEN CALL next_trip(trip_rp1, ni, nj, i, j, nexti, nextj) IF (trip_rp1(nexti,nextj)>=100) trip_rp1(i,j)=98 ! unrouted point becomes coastal flow ENDIF ENDDO ENDDO state_r(1:ni,1:nj) = state_rp1(1:ni,1:nj) trip_r(1:ni,1:nj) = trip_rp1(1:ni,1:nj) updated=1 it=2 DO WHILE(updated>0 .AND. it=1 .AND. i<=ni .AND. j>=1 .AND. j<=nj) trip_r(i,j)=98 ! => lakinflow will becomes coastalflow ELSE ! ok for riverflow ENDIF ELSE IF (state_rp1(i,j) == is_ter) THEN ! too short IF (trip_rp1(i,j)==97) THEN !ok for lakinflow ELSE ! => for riverflow/coastalflow IF (trip_extended_rp1(i,j)==97) THEN ! definitively a lake that cannot be routed further away IF (i>=1 .AND. i<=ni .AND. j>=1 .AND. j<=nj) THEN trip_r(i,j)=97 updated=updated+1 ENDIF ELSE CALL next_trip(trip_extended_rp1, ni, nj, i, j, nexti, nextj) IF (nexti/=-1 .AND. nextj/=-1) THEN IF (nexti==i .AND. nextj==j) THEN !! routed on itself in midle of land => lake IF (nexti>=1 .AND. nexti<=ni .AND. nextj>=1 .AND. nextj<=nj) trip_r(i,j)=97 ELSE IF (state_rp1(nexti,nextj)==is_coast .OR. state_rp1(nexti,nextj)==is_ter) THEN IF (i>=1 .AND. i<=ni .AND. j>=1 .AND. j<=nj) THEN trip_r(i,j)=trip_extended_rp1(i,j) state_r(i,j)= is_ter topoind_r(i,j) = 1e-10 ! flow instantanesly updated=updated+1 ENDIF IF (nexti>=1 .AND. nexti<=ni .AND. nextj>=1 .AND. nextj<=nj) THEN IF (trip_r(nexti,nextj)/=trip_rp1(i,j)) THEN trip_r(nexti,nextj)=trip_rp1(i,j) updated=updated+1 ENDIF ENDIF ELSE IF (i>=1 .AND. i<=ni .AND. j>=1 .AND. j<=nj) THEN STOP 'not possible' ENDIF ENDIF ENDIF ENDIF ENDIF ENDIF ELSE IF (state_rp1(i,j) == is_oce) THEN ! too long DO k=-1,1 DO m=-1,1 IF (k==0 .AND. m==0) CYCLE prevj = j+k ; previ = i+m CALL next_trip(trip_rp1, ni, nj, previ, prevj, nexti, nextj) ! IF (nexti==i .AND. nextj==j .AND. previ>=1 .AND. previ<=ni .AND. prevj>=1 .AND. prevj<=nj) trip_r(previ,prevj)=trip_rp1(i,j) IF (nexti==i .AND. nextj==j .AND. previ>=1 .AND. previ<=ni .AND. prevj>=1 .AND. prevj<=nj) trip_r(previ,prevj)=trip_rp1(i,j) ENDDO ENDDO IF (i>=1 .AND. i<=ni .AND. j>=1 .AND. j<=nj) THEN trip_r(i,j)=1000 updated=updated+1 ENDIF ENDIF ENDIF ENDDO ENDDO CALL MPI_ALLREDUCE(MPI_IN_PLACE , updated, 1, MPI_INT_ORCH, MPI_SUM, MPI_COMM_ORCH, ierr) it = it +1 ENDDO CALL xios_send_field("trip_r_final",trip_r) CALL xios_context_finalize CALL xios_orchidee_change_context("orchidee") CONTAINS SUBROUTINE next_trip(trip, ni, nj, i, j, nexti, nextj) IMPLICIT NONE REAL, INTENT(IN) :: trip(0:ni+1,0:nj+1) INTEGER, INTENT(IN) :: ni INTEGER, INTENT(IN) :: nj INTEGER, INTENT(IN) :: i INTEGER, INTENT(IN) :: j INTEGER, INTENT(OUT) :: nexti INTEGER, INTENT(OUT) :: nextj IF (i<0 .OR. i>ni+1 .OR. j<0 .OR. j>nj+1) THEN nexti=-1 nextj=-1 RETURN ENDIF SELECT CASE(NINT(trip(i,j))) CASE (1) nextj=j-1 ; nexti=i ! north CASE (2) nextj=j-1 ; nexti=i+1 ! north-east CASE (3) nextj=j ; nexti=i+1 ! east CASE (4) nextj=j+1 ; nexti=i+1 ! south-east CASE (5) nextj=j+1 ; nexti=i ! south CASE(6) nextj=j+1 ; nexti=i-1 ! south-west CASE (7) nextj=j ; nexti=i-1 ! west CASE (8) nextj=j-1 ; nexti=i-1 ! north-west CASE DEFAULT nextj=-1 ; nexti=-1 ! undefined END SELECT IF (nexti<0 .OR. nexti>ni+1) nexti=-1 IF (nextj<0 .OR. nextj>nj+1) nextj=-1 END SUBROUTINE next_trip END SUBROUTINE routing_flow_correct_riverflow SUBROUTINE routing_flow_init_mean(kjit, rest_id) USE ioipsl IMPLICIT NONE INTEGER(i_std), INTENT(in) :: kjit INTEGER(i_std), INTENT(in) :: rest_id !! Restart file identifier (unitless) INTEGER :: ier CHARACTER(LEN=80) :: var_name !! To store variables names for I/O (unitless) ! Get from the restart the fluxes we accumulated. ! ALLOCATE (runoff_mean(nbpt), stat=ier) runoff_mean(:) = 0 IF (ier /= 0) CALL ipslerr_p(3,'routing_flow_init_mean','Pb in allocate for runoff_mean','','') var_name = 'runoff_route' CALL ioconf_setatt_p('UNITS', 'Kg') CALL ioconf_setatt_p('LONG_NAME','Accumulated runoff for routing') CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., runoff_mean, "gather", nbp_glo, index_g) CALL setvar_p (runoff_mean, val_exp, 'NO_KEYWORD', zero) ALLOCATE(drainage_mean(nbpt), stat=ier) drainage_mean(:) = 0 IF (ier /= 0) CALL ipslerr_p(3,'routing_flow_init_mean','Pb in allocate for drainage_mean','','') var_name = 'drainage_route' CALL ioconf_setatt_p('UNITS', 'Kg') CALL ioconf_setatt_p('LONG_NAME','Accumulated drainage for routing') CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., drainage_mean, "gather", nbp_glo, index_g) CALL setvar_p (drainage_mean, val_exp, 'NO_KEYWORD', zero) END SUBROUTINE routing_flow_init_mean SUBROUTINE routing_flow_make_mean(runoff, drainage) IMPLICIT NONE REAL(r_std),INTENT(IN) :: runoff(:) REAL(r_std),INTENT(IN) :: drainage(:) runoff_mean(:) = runoff_mean(:) + runoff drainage_mean(:) = drainage_mean(:) + drainage END SUBROUTINE routing_flow_make_mean SUBROUTINE routing_flow_reset_mean IMPLICIT NONE runoff_mean(:) = 0 drainage_mean(:) = 0 END SUBROUTINE routing_flow_reset_mean SUBROUTINE routing_flow_finalize_mean(kjit, rest_id) USE ioipsl IMPLICIT NONE INTEGER,INTENT(IN) :: kjit INTEGER,INTENT(IN) :: rest_id CALL restput_p (rest_id, 'runoff_route', nbp_glo, 1, 1, kjit, runoff_mean, 'scatter', nbp_glo, index_g) CALL restput_p (rest_id, 'drainage_route', nbp_glo, 1, 1, kjit, drainage_mean, 'scatter', nbp_glo, index_g) DEALLOCATE(runoff_mean) DEALLOCATE(drainage_mean) END SUBROUTINE routing_flow_finalize_mean SUBROUTINE routing_flow_init_local(contfrac, nbpt_r) USE xios USE grid, ONLY : area IMPLICIT NONE INCLUDE "mpif.h" !! 0 Variable and parameter description !! 0.1 Input variables REAL(r_std),INTENT(IN) :: contfrac(nbpt) !! fraction of land INTEGER,INTENT(OUT) :: nbpt_r !! nb points routing native grid !! 0.2 Local variables INTEGER :: ni !! longitude dimension of local routing grid INTEGER :: nj !! latitude dimension of local routing grid REAL(r_std) :: contfrac_mpi(nbp_mpi) REAL(r_std),ALLOCATABLE :: trip_rp1(:,:) !! direction of flow (1-8) or river flow (99) or coastal flow (98) or lake inflow (97) - local routing grid + halo (0:ni+1,0:nj+1) REAL(r_std),ALLOCATABLE :: trip_extended_r(:,:) !! direction of flow (1-8) or river flow (99) or coastal flow (98) or lake inflow (97) - local routing grid (ni,nj) !! routing is artificially computed on sea and endoric basins REAL(r_std),ALLOCATABLE :: diag_r(:) !! fraction of routing cell intesected by coastal cells of native grid REAL(r_std),ALLOCATABLE :: frac_lake_r(:) !! fraction of routing cell intesected by cells of native grid REAL(r_std),ALLOCATABLE :: frac_coast_r(:) REAL(r_std) :: diag(nbp_mpi) !! fraction of routing cell intesected by coastal cells of native grid LOGICAL :: coastline(nbpt) LOGICAL :: coastline_mpi(nbp_mpi) REAL(r_std) :: area_mpi(nbp_mpi) INTEGER :: ij, ij_r, ij_rp1, i,j,jp1,jm1,ip1,im1,jr,ir INTEGER :: basins_count_mpi INTEGER :: nb_coast_points INTEGER :: ierr LOGICAL :: file_exists REAL(r_std) :: epsilon = 1e-5 !_ ================================================================================================================================ ! !> A value for property of each reservoir (in day/m) is given to compute a time constant (in day) !> for each reservoir (product of tcst and topo_resid). !> The value of tcst has been calibrated for the three reservoirs over the Senegal river basin only, !> during the 1 degree NCEP Corrected by Cru (NCC) resolution simulations (Ngo-Duc et al., 2005, Ngo-Duc et al., 2006) and !> generalized for all the basins of the world. The "slow reservoir" and the "fast reservoir" !> have the highest value in order to simulate the groundwater. !> The "stream reservoir", which represents all the water of the stream, has the lowest value. !> Those figures are the same for all the basins of the world. CALL getin_p('SLOW_TCST', slow_tcst) ! !Config Key = FAST_TCST !Config Desc = Time constant for the fast reservoir !Config If = RIVER_ROUTING !Config Def = 3.0 !Config Help = This parameters allows the user to fix the !Config time constant (in days) of the fast reservoir !Config in order to get better river flows for !Config particular regions. !Config Units = [days] CALL getin_p('FAST_TCST', fast_tcst) ! !Config Key = STREAM_TCST !Config Desc = Time constant for the stream reservoir !Config If = RIVER_ROUTING !Config Def = 0.24 !Config Help = This parameters allows the user to fix the !Config time constant (in days) of the stream reservoir !Config in order to get better river flows for !Config particular regions. !Config Units = [days] CALL getin_p('STREAM_TCST', stream_tcst) ! !Config Key = FLOOD_TCST !Config Desc = Time constant for the flood reservoir !Config If = RIVER_ROUTING !Config Def = 4.0 !Config Help = This parameters allows the user to fix the !Config time constant (in days) of the flood reservoir !Config in order to get better river flows for !Config particular regions. !Config Units = [days] split_routing=1 CALL getin_p("SPLIT_ROUTING",split_routing) CALL compute_coastline(contfrac, coastline) CALL gather_omp(contfrac, contfrac_mpi) CALL gather_omp(area, area_mpi) CALL gather_omp(coastline, coastline_mpi) ALLOCATE(is_coastline(nbp_mpi)) is_coastline=coastline_mpi IF (is_omp_root) THEN WHERE(contfrac_mpi<=epsilon) contfrac_mpi=0 WHERE(contfrac_mpi>=epsilon) contfrac_mpi=1 diag=0 WHERE(coastline_mpi) diag=1 CALL xios_send_field("is_coastline", diag) CALL xios_get_domain_attr("routing_domain", ni=ni, nj=nj) ! get routing domain dimension nbpt_r= ni*nj nbpt_rp1= (ni+2)*(nj+2) ! Allocate module variable ALLOCATE(fast_reservoir_r(ni*nj)) ALLOCATE(slow_reservoir_r(ni*nj)) ALLOCATE(stream_reservoir_r(ni*nj)) ALLOCATE(is_lakeinflow_r(ni*nj)) ALLOCATE(is_coastalflow_r(ni*nj)) ALLOCATE(is_riverflow_r(ni*nj)) ALLOCATE(is_streamflow_r(ni*nj)) ALLOCATE(topoind_r(ni*nj)) ALLOCATE(routing_mask_r(ni*nj)) ALLOCATE(basins_extended_r(nbpt_r)) ALLOCATE(route_flow_rp1((ni+2)*(nj+2))) ALLOCATE(routing_weight(nbp_mpi)) ALLOCATE(routing_weight_in(nbp_mpi)) ALLOCATE(unrouted_weight(nbp_mpi)) ALLOCATE(diag_r(nbpt_r)) ! for diags ALLOCATE(frac_coast_r(nbpt_r)) ALLOCATE(frac_lake_r(nbpt_r)) ALLOCATE(weight_coast_to_coast_r(nbpt_r)) ALLOCATE(weight_coast_to_lake_r(nbpt_r)) ALLOCATE(weight_lake_to_coast_r(nbpt_r)) ALLOCATE(weight_lake_to_lake_r(nbpt_r)) ALLOCATE(trip_extended_r(ni,nj)) is_lakeinflow_r(:) = .FALSE. is_coastalflow_r(:) = .FALSE. is_riverflow_r(:) = .FALSE. INQUIRE(FILE="routing_start.nc", EXIST=file_exists) IF (file_exists) THEN CALL xios_recv_field("fast_reservoir_start",fast_reservoir_r) CALL xios_recv_field("slow_reservoir_start",slow_reservoir_r) CALL xios_recv_field("stream_reservoir_start",stream_reservoir_r) ELSE fast_reservoir_r(:)=0 slow_reservoir_r(:)=0 stream_reservoir_r(:)=0 ENDIF ALLOCATE(trip_rp1(0:ni+1,0:nj+1)) trip_rp1(:,:)=1e10 CALL xios_recv_field("trip_r",trip_rp1(1:ni,1:nj)) ! recv trip array with halo of 1 CALL xios_recv_field("trip_extended_r",trip_extended_r) ! recv extended trip array from file CALL xios_recv_field("topoind_r",topoind_r) ! recv topo index array from file CALL xios_recv_field("basins_extended_r",basins_extended_r) ! recv basins index from file CALL routing_flow_correct_riverflow(ni, nj, contfrac_mpi, coastline, trip_rp1(1:ni,1:nj), trip_extended_r, topoind_r) CALL xios_send_field("trip_update_r",trip_rp1(1:ni,1:nj)) ! send to xios trip array to update for halo CALL xios_recv_field("trip_rp1",trip_rp1) ! recv trip array with halo of 1 !! Compute the routing !! Loop on all point point of the local routing grid + halo route_flow_rp1(:)=-1 DO j=0,nj+1 jp1=j+1 jm1=j-1 DO i=0,ni+1 ij_rp1=i+(ni+2)*j+1 ij_r=i+ni*(j-1) ip1=i+1 im1=i-1 IF (trip_rp1(i,j) < 100) THEN IF (i>=1 .AND. i<=ni .AND. j>=1 .AND. j<=nj) routing_mask_r(ij_r)=.TRUE. ir=-1 ; jr=-1 ! -> -1 for 97,98,99 SELECT CASE (NINT(trip_rp1(i,j))) ! get the trip value for each points CASE (1) jr=jm1 ; ir=i ! north CASE (2) jr=jm1 ; ir=ip1 ! north-east CASE (3) jr=j ; ir=ip1 ! east CASE (4) jr=jp1 ; ir=ip1 ! south-east CASE (5) jr=jp1 ; ir=i ! south CASE(6) jr=jp1 ; ir=im1 ! south-west CASE (7) jr=j ; ir=im1 ! west CASE (8) jr=jm1 ; ir=im1 ! north-west CASE (97) IF ( i>0 .AND. i0 .AND. j0 .AND. i0 .AND. j0 .AND. i0 .AND. jni+1 .OR. jr<0 .OR. jr>nj+1) THEN route_flow_rp1(ij_rp1)=-1 ! if route outside my local domain+halo, no routing (will be done by other process) ELSE IF (ir<1 .OR. ir>ni .OR. jr<1 .OR. jr>nj) THEN route_flow_rp1(ij_rp1)=-1 ! if route outside my local domain, no routing (will be done by other process) ELSE route_flow_rp1(ij_rp1)=ir+ni*(jr-1) ! define the cell where to flow IF (trip_rp1(ir,jr)>99) STOP 'Pb point not routed to outflow' ENDIF ENDIF ELSE IF (i>=1 .AND. i<=ni .AND. j>=1 .AND. j<=nj) routing_mask_r(ij_r)=.FALSE. ENDIF ENDDO ENDDO is_streamflow_r(:) = .NOT. (is_lakeinflow_r(:) .OR. is_coastalflow_r(:) .OR. is_riverflow_r(:)) .AND. routing_mask_r(:) diag_r(:)=0 WHERE(routing_mask_r) diag_r=1 CALL xios_send_field("routing_mask_r", diag_r) CALL xios_recv_field("frac_routing", diag) routing_weight_in (:) = 0 unrouted_weight (:) = 0 WHERE (diag > 0) routing_weight_in = contfrac_mpi*area_mpi/diag WHERE (diag == 0 ) unrouted_weight = contfrac_mpi*area_mpi ! compute the number of coast cells to distribute lakeinflow onto diag=0 WHERE (coastline_mpi) diag=1 nb_coast_points=SUM(diag) CALL reduce_sum_mpi(nb_coast_points, total_coast_points) CALL bcast_mpi(total_coast_points) DO ij=1,nbp_mpi routing_weight(ij)=contfrac_mpi(ij)*area_mpi(ij) ENDDO ! looking for basins basins_count_mpi=0 DO ij=1,nbpt_r IF (basins_extended_r(ij) > basins_count_mpi) basins_count_mpi=basins_extended_r(ij) ENDDO CALL MPI_ALLREDUCE(basins_count_mpi,basins_count,1,MPI_INT_ORCH, MPI_MAX,MPI_COMM_ORCH,ierr) diag(:) = 0 WHERE (coastline_mpi) diag=1 CALL xios_send_field("mask_coastal",diag) CALL xios_recv_field("frac_coastal_r",diag_r) WHERE (diag_r > 100.) ! missing_value diag_r=0 ELSE WHERE (diag_r < epsilon) diag_r=0 ELSEWHERE (diag_r > 1-epsilon) diag_r=1 END WHERE frac_coast_r=diag_r ; diag(:) = 0 WHERE (.NOT. coastline_mpi) diag=1 CALL xios_send_field("mask_lake",diag) CALL xios_recv_field("frac_lake_r",diag_r) WHERE (diag_r > 100.) ! missing_value diag_r=0. ELSEWHERE (diag_r < epsilon) diag_r=0. ELSEWHERE (diag_r > 1-epsilon) diag_r=1. END WHERE frac_lake_r = diag_r weight_coast_to_coast_r(:)=0 weight_coast_to_lake_r(:)=0 weight_lake_to_coast_r(:)=0 weight_lake_to_lake_r(:)=0 DO ij=1,nbpt_r IF (is_riverflow_r(ij) .OR. is_coastalflow_r(ij) ) THEN IF (frac_coast_r(ij)==0 .AND. frac_lake_r(ij)==0) THEN PRINT*,"riverflow or costalflow point on routing grid can not be interpolated on orchidee grid, this migh not happen" STOP ELSE IF (frac_coast_r(ij)==0 .AND. frac_lake_r(ij)>0) THEN weight_coast_to_lake_r(ij) = 1./frac_lake_r(ij) weight_coast_to_coast_r(ij) = 0 ELSE IF (frac_coast_r(ij)>0 ) THEN weight_coast_to_coast_r(ij) = 1./(frac_coast_r(ij) + frac_lake_r(ij)) weight_coast_to_lake_r(ij) = 0 ENDIF ELSE IF (is_lakeinflow_r(ij)) THEN IF ( frac_coast_r(ij)==0 .AND. frac_lake_r(ij)==0) THEN PRINT*,"lakeinflow on routing grid can not be interpolated on orchidee grid, this might not happen" STOP ELSE IF (frac_lake_r(ij)==0 .AND. frac_coast_r(ij)>0) THEN weight_lake_to_coast_r(ij) = 1./frac_coast_r(ij) weight_lake_to_lake_r(ij) = 0 ELSE IF (frac_lake_r(ij)>0 ) THEN weight_lake_to_lake_r(ij) = 1./(frac_coast_r(ij) + frac_lake_r(ij)) weight_lake_to_coast_r(ij) = 0 ENDIF ENDIF ENDDO CALL compute_basins_area(nbpt_r,nbpt_rp1) ELSE nbpt_r=0 ENDIF !is_omp_root END SUBROUTINE routing_flow_init_local SUBROUTINE compute_basins_area(nbpt_r,nbpt_rp1 ) USE xios USE routing_native_para IMPLICIT NONE INCLUDE 'mpif.h' INTEGER, INTENT(IN) :: nbpt_r INTEGER, INTENT(IN) :: nbpt_rp1 REAL(r_std),PARAMETER :: PI=atan(1.)*4 REAL(r_std),PARAMETER :: earth_radius=6.371009E6 ! in meter INTEGER :: ni, nj REAL(r_std),ALLOCATABLE :: lon(:), lat(:) REAL(r_std) :: area_rp1(nbpt_rp1) REAL(r_std) :: send_area_rp1(nbpt_rp1) REAL(r_std) :: recv_area_r(nbpt_rp1) REAL :: delta_lon, delta_lat INTEGER :: i,j,ij,it,ig,ierr REAL(r_std) :: sum_area ALLOCATE(basins_area_r(nbpt_r)) CALL xios_get_domain_attr("routing_domain", ni=ni, nj=nj) ALLOCATE(lon(ni),lat(nj)) CALL xios_get_domain_attr("routing_domain",lonvalue_1d=lon, latvalue_1d=lat) ! compute routing cell area delta_lon=ABS(lon(2)-lon(1))*Pi/180. delta_lat=ABS(lat(2)-lat(1))*Pi/180. DO j=1,nj DO i=1,ni ig=(j+1-1)*(ni+2)+i+1 area_rp1(ig) = delta_lon*delta_lat*cos(lat(j)*Pi/180.)*earth_radius*earth_radius ! pretty good approximation ij = (j-1)*ni+i basins_area_r(ij) = delta_lon*delta_lat*cos(lat(j)*Pi/180.)*earth_radius*earth_radius ! pretty good approximation ENDDO ENDDO CALL update_halo(area_rp1) send_area_rp1(:) = area_rp1(:) sum_area=1. DO WHILE (sum_area/=0) recv_area_r(:) = 0 DO ig=1,nbpt_rp1 IF ( route_flow_rp1(ig) > 0 ) THEN basins_area_r(route_flow_rp1(ig)) = basins_area_r(route_flow_rp1(ig)) + send_area_rp1(ig) recv_area_r(route_flow_rp1(ig)) = recv_area_r(route_flow_rp1(ig)) + send_area_rp1(ig) ENDIF ENDDO DO j=1,nj DO i=1,ni ij=((j-1)*ni)+i ig=(j+1-1)*(ni+2)+i+1 send_area_rp1(ig) = recv_area_r(ij) IF (is_riverflow_r(ij) .OR. is_coastalflow_r(ij) .OR. is_lakeinflow_r(ij)) send_area_rp1(ig)=0 ENDDO ENDDO CALL update_halo(send_area_rp1) sum_area = sum(send_area_rp1) CALL MPI_ALLREDUCE(MPI_IN_PLACE, sum_area, 1, MPI_REAL_ORCH, MPI_SUM, MPI_COMM_ORCH, ierr) ENDDO END SUBROUTINE compute_basins_area SUBROUTINE initialize_stations(dt_routing) USE xios IMPLICIT NONE INCLUDE 'mpif.h' REAL(r_std), INTENT (in) :: dt_routing !! Routing time step (s) REAL(r_std),PARAMETER :: PI=atan(1.)*4 REAL(r_std),PARAMETER :: earth_radius=6.371009E6 ! in meter INTEGER :: ni, nj REAL(r_std),ALLOCATABLE :: lon(:), lat(:) INTEGER :: nb_mouth CHARACTER(LEN=60),ALLOCATABLE :: mouth(:) INTEGER,ALLOCATABLE :: mouth_id(:) INTEGER,ALLOCATABLE :: mouth_index(:) REAL,ALLOCATABLE :: station_lonlat(:,:) INTEGER,ALLOCATABLE :: station_prec(:) INTEGER ::i,j,r,k TYPE(xios_duration) :: ts TYPE(xios_scalar) :: scalar_hdl TYPE(xios_scalargroup) :: scalar_def_hdl TYPE(xios_field) :: field_hdl TYPE(xios_fieldgroup) :: field_def_hdl, fieldgroup_hdl TYPE(xios_file) :: file_hdl TYPE(xios_filegroup) :: file_def_hdl TYPE(xios_date) :: start_date, time_origin CHARACTER(LEN=20) :: calendar_type CHARACTER(LEN=256) :: str_station_ind REAL :: lon_a,lat_a,lon_b,lat_b,dist,max_area, max_max_area, min_dist INTEGER :: min_dist_index INTEGER :: ierr CALL xios_get_domain_attr("routing_domain", ni=ni, nj=nj) ALLOCATE(lon(ni),lat(nj)) CALL xios_get_domain_attr("routing_domain",lonvalue_1d=lon, latvalue_1d=lat) nb_station=0 CALL getin("nb_station",nb_station) ALLOCATE(station(nb_station)) ALLOCATE(station_lonlat(nb_station,2)) ALLOCATE(station_index(nb_station)) ALLOCATE(station_prec(nb_station)) DO k=1,nb_station WRITE(str_station_ind,*) k str_station_ind=ADJUSTL(str_station_ind) CALL getin("station"//TRIM(str_station_ind)//"_id",station(k)) CALL getin("station"//TRIM(str_station_ind)//"_coor",station_lonlat(k,:)) station_prec(k)=50000 CALL getin("station"//TRIM(str_station_ind)//"_prec",station_prec(k)) lon_a = station_lonlat(k,1)*Pi/180. lat_a = station_lonlat(k,2)*Pi/180. max_area=0 min_dist=HUGE(min_dist) DO j=1,nj DO i=1,ni r=((j-1)*ni)+i IF (routing_mask_r(r)) THEN lon_b=lon(i)*Pi/180. lat_b=lat(j)*Pi/180. dist=earth_radius*acos(sin(lon_a)*sin(lon_b)+cos(lon_a)*cos(lon_b)*cos(lat_b-lat_a)) IF (dist max_area) THEN max_area = basins_area_r(r) station_index(k) = r ENDIF ENDIF ENDIF ENDDO ENDDO CALL MPI_ALLREDUCE(max_area, max_max_area, 1, MPI_REAL_ORCH, MPI_MAX, MPI_COMM_ORCH, ierr) IF (max_area /= max_max_area) station_index(k)=-1 IF (max_area == 0) station_index(k)=-1 ! PRINT*,"Station ",TRIM(station(k))," old coor :",lon(min_dist_index),lat(min_dist_index)," basin_area (km2)", basins_area_r(min_dist_index)/1000/1000 ! PRINT*,"Station ",TRIM(station(k)),"new coor ",lon(station_index(k)),lat(station_index(k))," basin_area(km2) ", basins_area_r(station_index(k))/1000/1000 ENDDO CALL xios_get_calendar_type(calendar_type) CALL xios_get_start_date(start_date) CALL xios_get_time_origin(time_origin) CALL xios_context_initialize("orchidee_routing_out",MPI_COMM_ORCH) CALL xios_orchidee_change_context("orchidee_routing_out") ts.second=dt_routing calendar_type="gregorian" CALL xios_define_calendar(type=calendar_type,start_date=start_date,time_origin=time_origin, timestep=ts ) CALL xios_get_handle("scalar_definition",scalar_def_hdl) ! CALL xios_get_handle("mouths",field_def_hdl) ! DO k=1,nb_mouth ! CALL xios_add_child(scalar_def_hdl,scalar_hdl,TRIM(mouth(k))//"_scalar") ! CALL xios_set_attr(scalar_hdl,label=TRIM(mouth(k))) ! CALL xios_add_child(field_def_hdl,field_hdl, TRIM(mouth(k))) ! CALL xios_set_attr(field_hdl,scalar_ref=TRIM(mouth(k))//"_scalar") ! CALL xios_set_attr(field_hdl, freq_op = freq_op*xios_second) ! ENDDO CALL xios_get_handle("field_definition",field_def_hdl) CALL xios_add_child(field_def_hdl, fieldgroup_hdl, "stations") DO k=1,nb_station CALL xios_add_child(scalar_def_hdl,scalar_hdl,TRIM(station(k))//"_scalar") CALL xios_set_attr(scalar_hdl,label=TRIM(station(k))) CALL xios_add_child(fieldgroup_hdl,field_hdl, TRIM(station(k))) CALL xios_set_attr(field_hdl,scalar_ref=TRIM(station(k))//"_scalar") ENDDO CALL xios_get_handle("file_definition",file_def_hdl) CALL xios_add_child(file_def_hdl, file_hdl, "stations") CALL xios_set_attr(file_hdl, type="one_file", output_freq=ts, sync_freq=ts, enabled=.TRUE.) CALL xios_add_child(file_hdl, fieldgroup_hdl) CALL xios_set_attr(fieldgroup_hdl, group_ref="stations", operation="average") CALL xios_close_context_definition() CALL xios_orchidee_change_context("orchidee") END SUBROUTINE initialize_stations SUBROUTINE routing_flow_main(dt_routing) ! USE xios USE grid, ONLY : area USE routing_native_para IMPLICIT NONE INCLUDE "mpif.h" !! 0 Variable and parameter description !! 0.1 Input variables REAL(r_std), INTENT (in) :: dt_routing !! Routing time step (s) !! 0.4 Local variables REAL(r_std) :: runoff(nbp_mpi) !! Grid-point runoff (kg/dt) REAL(r_std) :: runoff_in(nbp_mpi) !! Grid-point runoff (kg/dt) REAL(r_std) :: drainage(nbp_mpi) !! Grid-point drainage (kg/dt) REAL(r_std) :: drainage_in(nbp_mpi) !! Grid-point drainage (kg/dt) REAL(r_std) :: riverflow(nbp_mpi) REAL(r_std) :: coastalflow(nbp_mpi) REAL(r_std) :: lakeinflow(nbp_mpi) REAL(r_std) :: fast_diag_mpi(nbp_mpi) REAL(r_std) :: slow_diag_mpi(nbp_mpi) REAL(r_std) :: stream_diag_mpi(nbp_mpi) REAL(r_std) :: area_mpi(nbp_mpi) ! cell area REAL(r_std) :: flow_coast(nbp_mpi) REAL(r_std) :: flow_lake(nbp_mpi) ! from input model -> routing_grid REAL(r_std) :: runoff_r(nbpt_r) !! Grid-point runoff (kg/m^2/dt) REAL(r_std) :: drainage_r(nbpt_r) !! Grid-point drainage (kg/m^2/dt) REAL(r_std), DIMENSION(nbpt_r) :: fast_flow_r !! Outflow from the fast reservoir (kg/dt) REAL(r_std), DIMENSION(nbpt_r) :: slow_flow_r !! Outflow from the slow reservoir (kg/dt) REAL(r_std), DIMENSION(nbpt_r) :: stream_flow_r !! Outflow from the stream reservoir (kg/dt) REAL(r_std), DIMENSION(nbpt_r) :: hydrographs_r !! hydrograph (kg/dt) REAL(r_std), DIMENSION(nbpt_r) :: transport_r !! Water transport between basins (kg/dt) INTEGER(i_std) :: ig !! Indices (unitless) INTEGER(i_std) :: isplit LOGICAL, PARAMETER :: check_reservoir = .TRUE. !! Logical to choose if we write informations when a negative amount of water is occurring in a reservoir (true/false) REAL(r_std), DIMENSION(nbpt_r) :: lakeinflow_r REAL(r_std), DIMENSION(nbpt_r) :: coastalflow_r REAL(r_std), DIMENSION(nbpt_r) :: riverflow_r REAL(r_std), DIMENSION(nbpt_r) :: flow_r REAL(r_std), DIMENSION(nbpt_rp1) :: flow_rp1 REAL(r_std) :: flow !! Outflow computation for the reservoirs (kg/dt) REAL(r_std) :: basins_riverflow_mpi(0:basins_count) REAL(r_std) :: basins_riverflow(0:basins_count) REAL(r_std) :: water_balance_before, water_balance_after REAL(r_std) :: sum_water_before, sum_water_after REAL(r_std) :: value INTEGER(i_std) :: k INTEGER :: ierr !_ ================================================================================================================================ !> The outflow fluxes from the three reservoirs are computed. !> The outflow of volume of water Vi into the reservoir i is assumed to be linearly related to its volume. !> The water travel simulated by the routing scheme is dependent on the water retention index topo_resid !> given by a 0.5 degree resolution map for each pixel performed from a simplification of Manning's formula !> (Dingman, 1994; Ducharne et al., 2003). !> The resulting product of tcst (in day/m) and topo_resid (in m) represents the time constant (day) !> which is an e-folding time, the time necessary for the water amount !> in the stream reservoir to decrease by a factor e. Hence, it gives an order of !> magnitude of the travel time through this reservoir between !> the sub-basin considered and its downstream neighbor. CALL gather_omp(runoff_mean,runoff) CALL gather_omp(drainage_mean, drainage) CALL gather_omp(area, area_mpi) IF (is_omp_root) THEN CALL xios_send_field("routing_basins_area",basins_area_r) hydrographs_r(:)=0 ! water balance before sum_water_before = sum((runoff(:)+drainage(:))*routing_weight(:))+sum(fast_reservoir_r(:)+slow_reservoir_r(:)+stream_reservoir_r(:)) CALL MPI_ALLREDUCE(sum_water_before, water_balance_before,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr) runoff_in=runoff*routing_weight_in CALL xios_send_field("routing_runoff",runoff_in) ! interp conservative model -> routing CALL xios_recv_field("routing_runoff_r",runoff_r) WHERE(.NOT. routing_mask_r) runoff_r = 0 drainage_in=drainage*routing_weight_in CALL xios_send_field("routing_drainage",drainage_in) ! interp conservative model -> routing CALL xios_recv_field("routing_drainage_r",drainage_r) WHERE(.NOT. routing_mask_r) drainage_r = 0 CALL MPI_ALLREDUCE(sum(runoff*routing_weight),sum_water_before,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr) CALL MPI_ALLREDUCE(sum(runoff_r),sum_water_after,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr) IF (sum_water_after/=0) THEN IF (is_mpi_root) PRINT *,"runoff fixer : ", sum_water_before/sum_water_after DO ig=1,nbpt_r runoff_r(ig) = runoff_r(ig) * sum_water_before/sum_water_after ENDDO ENDIF CALL MPI_ALLREDUCE(sum(drainage*routing_weight),sum_water_before,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr) CALL MPI_ALLREDUCE(sum(drainage_r),sum_water_after,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr) IF (sum_water_after/=0) THEN IF (is_mpi_root) PRINT *,"drainage fixer : ", sum_water_before/sum_water_after DO ig=1,nbpt_r drainage_r(ig) = drainage_r(ig) * sum_water_before/sum_water_after ENDDO ENDIF CALL MPI_ALLREDUCE(sum(runoff*routing_weight)+sum(drainage*routing_weight),sum_water_before,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr) CALL MPI_ALLREDUCE(sum(runoff_r) +sum(drainage_r),sum_water_after,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr) IF (is_mpi_root) PRINT *,"Loose water by interpolation ; before : ", sum_water_before," ; after : ",sum_water_after, & " ; delta : ", 100.*(sum_water_after-sum_water_before)/(0.5*(sum_water_after+sum_water_before)),"%" runoff_in(:) = runoff(:)*unrouted_weight(:) drainage_in(:) = drainage(:)*unrouted_weight(:) runoff_r(:) = runoff_r(:) / split_routing drainage_r(:) = drainage_r(:) / split_routing hydrographs_r(:) = 0 transport_r(:)=0 DO isplit=1,split_routing DO ig=1,nbpt_r IF ( routing_mask_r(ig) ) THEN flow = MIN(fast_reservoir_r(ig)/((topoind_r(ig)/1000.)*fast_tcst*one_day/(dt_routing/split_routing)),& & fast_reservoir_r(ig)-min_sechiba) fast_flow_r(ig) = MAX(flow, zero) ! flow = MIN(slow_reservoir_r(ig)/((topoind_r(ig)/1000.)*slow_tcst*one_day/(dt_routing/split_routing)),& & slow_reservoir_r(ig)-min_sechiba) slow_flow_r(ig) = MAX(flow, zero) ! flow = MIN(stream_reservoir_r(ig)/((topoind_r(ig)/1000.)*stream_tcst*one_day/(dt_routing/split_routing)),& & stream_reservoir_r(ig)-min_sechiba) stream_flow_r(ig) = MAX(flow, zero) ! ! ELSE fast_flow_r(ig) = zero slow_flow_r(ig) = zero stream_flow_r(ig) = zero ENDIF ENDDO !- !- Compute the transport !- DO ig=1,nbpt_r flow_rp1(r_to_rp1(ig))=fast_flow_r(ig) + slow_flow_r(ig) + stream_flow_r(ig) ENDDO CALL update_halo(flow_rp1) ! transfert halo DO ig=1,nbpt_rp1 IF ( route_flow_rp1(ig) > 0 ) THEN transport_r(route_flow_rp1(ig))=transport_r(route_flow_rp1(ig))+ flow_rp1(ig) ENDIF ENDDO DO ig=1,nbpt_r IF ( routing_mask_r(ig) ) THEN fast_reservoir_r(ig) = fast_reservoir_r(ig) + runoff_r(ig) - fast_flow_r(ig) slow_reservoir_r(ig) = slow_reservoir_r(ig) + drainage_r(ig) - slow_flow_r(ig) stream_reservoir_r(ig) = stream_reservoir_r(ig) - stream_flow_r(ig) IF (is_streamflow_r(ig)) THEN stream_reservoir_r(ig)= stream_reservoir_r(ig) + transport_r(ig) transport_r(ig) = 0 ENDIF IF ( stream_reservoir_r(ig) .LT. zero ) THEN IF ( check_reservoir ) THEN WRITE(numout,*) "WARNING : negative stream reservoir at :", ig, ". Problem is being corrected." WRITE(numout,*) "stream_reservoir, transport, stream_flow,: ", & stream_reservoir_r(ig), transport_r(ig), stream_flow_r(ig) ENDIF fast_reservoir_r(ig) = fast_reservoir_r(ig) + stream_reservoir_r(ig) stream_reservoir_r(ig) = zero ENDIF ! IF ( fast_reservoir_r(ig) .LT. zero ) THEN IF ( check_reservoir ) THEN WRITE(numout,*) "WARNING : negative fast reservoir at :", ig, ". Problem is being corrected." WRITE(numout,*) "fast_reservoir, runoff, fast_flow : ", fast_reservoir_r(ig), & &runoff_r(ig), fast_flow_r(ig) ENDIF slow_reservoir_r(ig) = slow_reservoir_r(ig) + fast_reservoir_r(ig) fast_reservoir_r(ig) = zero ENDIF IF ( slow_reservoir_r(ig) .LT. - min_sechiba ) THEN WRITE(numout,*) 'WARNING : There is a negative reservoir at :', ig WRITE(numout,*) 'WARNING : slowr, slow_flow, drainage', & & slow_reservoir_r(ig), slow_flow_r(ig), drainage_r(ig) CALL ipslerr_p(2, 'routing_simple_flow', 'WARNING negative slow_reservoir.','','') ENDIF ENDIF ENDDO DO ig=1,nbpt_r IF ( routing_mask_r(ig) ) THEN hydrographs_r(ig)=hydrographs_r(ig)+fast_flow_r(ig)+slow_flow_r(ig)+stream_flow_r(ig) ENDIF ENDDO ENDDO ! isplit lakeinflow_r(:)=0 coastalflow_r(:)=0 riverflow_r(:)=0 basins_riverflow_mpi(:)=0 DO ig=1,nbpt_r IF ( routing_mask_r(ig) ) THEN IF (is_lakeinflow_r(ig)) THEN lakeinflow_r(ig) = transport_r(ig) basins_riverflow_mpi(basins_extended_r(ig)) = basins_riverflow_mpi(basins_extended_r(ig))+lakeinflow_r(ig) ENDIF IF (is_coastalflow_r(ig)) THEN coastalflow_r(ig) = transport_r(ig) basins_riverflow_mpi(basins_extended_r(ig)) = & basins_riverflow_mpi(basins_extended_r(ig))+coastalflow_r(ig) ENDIF IF (is_riverflow_r(ig)) THEN riverflow_r(ig) = transport_r(ig) basins_riverflow_mpi(basins_extended_r(ig)) = & basins_riverflow_mpi(basins_extended_r(ig))+riverflow_r(ig) ENDIF ENDIF ENDDO ! now send riverflow, coastalflow and lakeinflow on orchidee grid coastalflow(:)=0. riverflow(:)=0. lakeinflow(:)=0. CALL xios_send_field("routing_coastalflow_to_coast_r" ,coastalflow_r*weight_coast_to_coast_r) CALL xios_recv_field("routing_coastalflow_to_coast" ,flow_coast) WHERE(.NOT.is_coastline) flow_coast=0 CALL xios_send_field("routing_coastalflow_to_lake_r" ,coastalflow_r*weight_coast_to_lake_r) CALL xios_recv_field("routing_coastalflow_to_lake" ,flow_lake) WHERE(is_coastline) flow_lake=0. CALL MPI_ALLREDUCE(sum(coastalflow_r),sum_water_before,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr) CALL MPI_ALLREDUCE(sum(flow_coast+flow_lake),sum_water_after,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr) IF (sum_water_after/=0) THEN IF (is_mpi_root) PRINT *,"coastalflow fixer : ", sum_water_before/sum_water_after flow_coast(:)=flow_coast(:)*(sum_water_before/sum_water_after) flow_lake(:)=flow_lake(:)*(sum_water_before/sum_water_after) ENDIF coastalflow=coastalflow+flow_coast lakeinflow=lakeinflow+flow_lake CALL xios_send_field("routing_riverflow_to_coast_r" ,riverflow_r*weight_coast_to_coast_r) CALL xios_recv_field("routing_riverflow_to_coast" ,flow_coast) WHERE(.NOT.is_coastline) flow_coast=0 CALL xios_send_field("routing_riverflow_to_lake_r" ,riverflow_r*weight_coast_to_lake_r) CALL xios_recv_field("routing_riverflow_to_lake" ,flow_lake) WHERE(is_coastline) flow_lake=0. CALL MPI_ALLREDUCE(sum(riverflow_r),sum_water_before,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr) CALL MPI_ALLREDUCE(sum(flow_coast+flow_lake),sum_water_after,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr) IF (sum_water_after/=0) THEN IF (is_mpi_root) PRINT *,"riverflow fixer : ", sum_water_before/sum_water_after flow_coast(:)=flow_coast(:)*(sum_water_before/sum_water_after) flow_lake(:)=flow_lake(:)*(sum_water_before/sum_water_after) ENDIF riverflow=riverflow+flow_coast lakeinflow=lakeinflow+flow_lake CALL xios_send_field("routing_lakeinflow_to_coast_r" ,lakeinflow_r*weight_lake_to_coast_r) CALL xios_recv_field("routing_lakeinflow_to_coast" ,flow_coast) WHERE(.NOT.is_coastline) flow_coast=0 CALL xios_send_field("routing_lakeinflow_to_lake_r" ,lakeinflow_r*weight_lake_to_lake_r) CALL xios_recv_field("routing_lakeinflow_to_lake" ,flow_lake) WHERE(is_coastline) flow_lake=0. CALL MPI_ALLREDUCE(sum(lakeinflow_r),sum_water_before,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr) CALL MPI_ALLREDUCE(sum(flow_coast+flow_lake),sum_water_after,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr) IF (sum_water_after/=0) THEN IF (is_mpi_root) PRINT *,"lakeinflow fixer : ", sum_water_before/sum_water_after flow_coast(:)=flow_coast(:)*(sum_water_before/sum_water_after) flow_lake(:)=flow_lake(:)*(sum_water_before/sum_water_after) ENDIF coastalflow=coastalflow+flow_coast lakeinflow=lakeinflow+flow_lake WHERE(is_coastline) coastalflow = coastalflow + runoff_in + drainage_in WHERE(.NOT.is_coastline) lakeinflow = lakeinflow + runoff_in + drainage_in sum_water_after = sum(coastalflow(:) + riverflow(:) + lakeinflow(:))+sum(fast_reservoir_r(:)+slow_reservoir_r(:)+stream_reservoir_r(:)) CALL MPI_ALLREDUCE(sum_water_after, water_balance_after,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr) IF (is_mpi_root) PRINT *,"routing water Balance ; before : ", water_balance_before," ; after : ",water_balance_after, & " ; delta : ", 100*(water_balance_after-water_balance_before)/(0.5*(water_balance_after+water_balance_before)),"%" ! diag CALL xios_send_field("routing_fast_reservoir_r" , fast_reservoir_r) CALL xios_send_field("routing_slow_reservoir_r" , slow_reservoir_r) CALL xios_send_field("routing_stream_reservoir_r" , stream_reservoir_r) CALL xios_send_field("routing_riverflow_r" , riverflow_r) CALL xios_send_field("routing_coastalflow_r" , coastalflow_r) CALL xios_send_field("routing_lakeinflow_r" , lakeinflow_r) CALL xios_send_field("out_flow",lakeinflow+coastalflow+riverflow) CALL xios_send_field("routing_hydrographs_r", (hydrographs_r+lakeinflow_r+coastalflow_r+riverflow_r)/1000./dt_routing) CALL xios_send_field("routing_riverflow" , riverflow) CALL xios_send_field("routing_coastalflow" , coastalflow) CALL xios_send_field("routing_lakeinflow" , lakeinflow) ! ! stations CALL xios_orchidee_change_context("orchidee_routing_out") station_ts = station_ts +1 CALL xios_update_calendar(station_ts) DO k=1,nb_station IF (station_index(k) /=-1) THEN value = hydrographs_r(station_index(k))/1000./dt_routing ELSE value = 0 ENDIF CALL MPI_ALLREDUCE(MPI_IN_PLACE,value,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr) CALL xios_send_field(TRIM(station(k)),value) ENDDO CALL xios_orchidee_change_context("orchidee") ENDIF ! is_omp_root CALL scatter_omp(riverflow,riverflow_mean) CALL scatter_omp(coastalflow,coastalflow_mean) CALL scatter_omp(lakeinflow,lakeinflow_mean) CALL routing_flow_reset_mean END SUBROUTINE routing_flow_main !! ============================================================================================================================= !! SUBROUTINE: routing_simple_finalize !! !>\BRIEF Write to restart file !! !! DESCRIPTION: Write module variables to restart file !! !! RECENT CHANGE(S) !! !! REFERENCE(S) !! !! FLOWCHART !! \n !_ ============================================================================================================================== SUBROUTINE routing_flow_finalize(kjit, rest_id) USE xios USE ioipsl IMPLICIT NONE INTEGER, INTENT(IN) :: kjit INTEGER, INTENT(IN) :: rest_id !_ ================================================================================================================================ IF (is_omp_root) THEN CALL xios_send_field("fast_reservoir_restart",fast_reservoir_r) CALL xios_send_field("slow_reservoir_restart",slow_reservoir_r) CALL xios_send_field("stream_reservoir_restart",stream_reservoir_r) ENDIF CALL restput_p (rest_id, 'riverflow', nbp_glo, 1, 1, kjit, riverflow_mean, 'scatter', nbp_glo, index_g) CALL restput_p (rest_id, 'coastalflow', nbp_glo, 1, 1, kjit, coastalflow_mean, 'scatter', nbp_glo, index_g) CALL restput_p (rest_id, 'lakeinflow', nbp_glo, 1, 1, kjit, lakeinflow_mean, 'scatter', nbp_glo, index_g) CALL routing_flow_finalize_mean(kjit, rest_id) CALL xios_orchidee_change_context("orchidee_routing_out") CALL xios_context_finalize() CALL xios_orchidee_change_context("orchidee") END SUBROUTINE routing_flow_finalize !! ================================================================================================================================ !! SUBROUTINE : routing_simple_clear !! !>\BRIEF This subroutine deallocates the block memory previously allocated. !! !! DESCRIPTION: This subroutine deallocates the block memory previously allocated. !! !! RECENT CHANGE(S): None !! !! MAIN OUTPUT VARIABLE(S): !! !! REFERENCES : None !! !! FLOWCHART :None !! \n !_ ================================================================================================================================ SUBROUTINE routing_flow_clear IMPLICIT NONE IF (is_omp_root) THEN IF (ALLOCATED(topoind_r)) DEALLOCATE(topoind_r) IF (ALLOCATED(route_flow_rp1)) DEALLOCATE(route_flow_rp1) IF (ALLOCATED(routing_mask_r)) DEALLOCATE(routing_mask_r) IF (ALLOCATED(fast_reservoir_r)) DEALLOCATE(fast_reservoir_r) IF (ALLOCATED(slow_reservoir_r)) DEALLOCATE(slow_reservoir_r) IF (ALLOCATED(is_lakeinflow_r)) DEALLOCATE(is_lakeinflow_r) IF (ALLOCATED(is_coastalflow_r)) DEALLOCATE(is_coastalflow_r) IF (ALLOCATED(is_riverflow_r)) DEALLOCATE(is_riverflow_r) ENDIF END SUBROUTINE routing_flow_clear END MODULE routing_native_flow_mod