#ifdef USE_LPJ
module lpj_driver_mod
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! !!
!! GNU General Public License !!
!! !!
!! This file is part of the Flexible Modeling System (FMS). !!
!! !!
!! FMS is free software; you can redistribute it and/or modify !!
!! it and are expected to follow the terms of the GNU General Public !!
!! License as published by the Free Software Foundation. !!
!! !!
!! FMS is distributed in the hope that it will be useful, !!
!! but WITHOUT ANY WARRANTY; without even the implied warranty of !!
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !!
!! GNU General Public License for more details. !!
!! !!
!! You should have received a copy of the GNU General Public License !!
!! along with FMS; if not, write to: !!
!! Free Software Foundation, Inc. !!
!! 59 Temple Place, Suite 330 !!
!! Boston, MA 02111-1307 USA !!
!! or see: !!
!! http://www.gnu.org/licenses/gpl.txt !!
!! !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Stefan Petri
!
!
!
!
!
! Wrapper for the vegetation model LPJ
!
!
! This module is the Fortran90 side of the wrapper for the
! vegetation model LPJ. It wraps the replacement for the main program of LPJ.
!
!#include
#define __ERROR__(message) \
call my_error(module_name,message,FATAL,thisfile,__LINE__)
#define __NOTE__(message) \
call my_error(module_name,message,NOTE,thisfile,__LINE__)
! CHECK is like ASSERT, but independent of NDEBUG
#define __CHECK__(x, message) \
if(.NOT.(x))call my_error(module_name,message,FATAL,thisfile,__LINE__)
#ifndef NDEBUG
#define __ASSERT__(x, message) \
if(.NOT.(x))call my_error(module_name,message,FATAL,thisfile,__LINE__)
#else
#define __ASSERT(x, message) ! empty
#endif
use time_manager_mod, only: time_type, operator(+), operator(==), assignment(=), get_time, get_date, increment_time
use mpp_domains_mod, only: domain1D, domain2d, mpp_get_layout, mpp_define_layout, &
mpp_define_domains, mpp_get_compute_domain, &
CYCLIC_GLOBAL_DOMAIN, mpp_get_compute_domains, &
mpp_get_domain_components, mpp_get_pelist, &
mpp_define_mosaic, mpp_global_field
use fms_mod, only: stdlog, open_namelist_file, check_nml_error, &
write_version_number, error_mesg, FATAL, NOTE, &
mpp_pe, mpp_npes, mpp_root_pe, close_file, &
mpp_clock_id, mpp_clock_begin, mpp_clock_end, &
clock_flag_default, CLOCK_COMPONENT, CLOCK_SUBCOMPONENT, &
field_size, read_data, field_exist, get_mosaic_tile_grid
use mpp_mod, only: mpp_error, FATAL, mpp_get_current_pelist, stdout, mpp_chksum, mpp_sum
use diag_manager_mod, only: diag_axis_init, register_static_field, register_diag_field, send_data, diag_manager_end
use field_manager_mod, only: MODEL_LAND
use tracer_manager_mod, only: register_tracers, get_tracer_index, NO_TRACER
use mosaic_mod, only: get_mosaic_ntiles, calc_mosaic_grid_area, &
get_mosaic_xgrid_size, get_mosaic_xgrid
use constants_mod, only: PI, SECONDS_PER_DAY, KELVIN, RAD_TO_DEG, DEG_TO_RAD, &
WTMAIR, & ! 28.96440 [AMU] molecular weight of air
WTMC, & ! 12.00000 [AMU] molecular weight of carbon
WTMCO2 ! 44.00995 [AMU] molecular weight of carbon dioxide
use land_types_mod, only: land_data_type
use numerics_mod, only: is_latlon
use topography_mod, only: get_topog_mean
!! To enable strict floating point exceptions inside this module.
!! Taken from the example in the Intel Fortran Compiler Reference Manual.
!! Most probably will not work with any other compiler
use ifcore
use ifport
implicit none
private
! ==== public interfaces =====================================================
public lpj_driver_init
public lpj_driver_update_fast ! , lpj_driver_update_slow
public update_lpj_bnd_fast, update_lpj_bnd_slow
public lpj_driver_end
!public lpj_driver_restart
! ==== some names, for information only ======================================
character(len=*), parameter :: &
module_name = 'lpj_driver_mod', &
version = '$Id$', &
tagname = '$Name$', &
thisfile = __FILE__
logical :: module_is_initialized = .FALSE.
! ==== public data type =====================================================
!
!
! The type describing the state of the LPJ model wrapper. There is only one
! variable of this type, called "lpj_state", which is used in the model to
! hold the data.
! Strictly speaking, this type is not necessary, but there is no harm in
! having it and it is possible to image a situation where a single entity
! describing the lpj model wrapper state would be useful.
!
type lpj_wrapper_type
private ! it's an opaque type: all data members are private
!
! Domain of calculations
!
type(domain2d) :: domain ! domain of calculations
! The update subroutines of the land model do not get the current time as arguments.
! Instead they must save the start time and time steps and advance the current time
! on each invocation.
type(time_type) :: current_time, dt_fast, dt_slow
integer :: dt_fast_seconds, dt_slow_seconds ! unpacked into seconds
! LPJ requires daily values. land_lad typically uses sub-daily time steps.
! Thus we aggregate here over one day's time steps
integer :: yesterday = -1
! The following data arrays are on the local compute domain, unless noted otherwise.
! In these, the input data from land_lad is accumulated over the time steps of one day
real, dimension(:,:), pointer :: &
!lprec => NULL(), & ! [kg/m2] == [mm] daily sum liquid precipitation
!fprec => NULL(), & ! [kg/m2] == [mm] daily sum frozen precipitation
! Currently, LPJ requires total precipitation and computes liquid and snow fractions
! on its own, depending on temperature.
! It is planned that LPJ can also make use of separate Xprec input directly.
! But that still needs to be implemented by th LPJ team.
prec => NULL(), & ! [kg/m2] == [mm] daily sum of lprec and fprec
temp_mean => NULL(), & ! [degC] daily mean
temp_max => NULL(), & ! [degC] daily max
temp_min => NULL(), & ! [degC] daily min
swdown => NULL(), & ! [W/m2] daily mean
lwnet => NULL(), & ! [W/m2] daily mean
t_flux => NULL(), & ! [W/m2] Sensible heat flux
qflux => NULL(), & ! [kg/m2/s] water vapor flux near surface
dedq => NULL(), & ! [kg/m2/s] Moisture flux humidity sensitivity
dhdt => NULL(), & ! derivative of sens over T
drdt => NULL(), & ! derivative of LW radiation over T
drag_q => NULL(), & ! [m/s] Product of drag coefficient for water vapor by wind
p_surf => NULL(), & ! [N/m2] == [Pa] == [kg/m/s2] Surface pressure
wind => NULL(), & ! [m/s] daily mean
co2 => NULL() ! CO2 concentration [kg CO2/kg wet air]; Note: LPJ wants [ppm] daily mean dry volume mixing ratio
! TODO: lightning
real, dimension(:,:), pointer :: lad_land_mask => NULL() ! land fraction, 0 = ocean
real, dimension(:,:), pointer :: lad_warea => NULL() ! Area covered by water per grid cell
real, dimension(:,:), pointer :: lpj_land_mask => NULL()
! LPJ output data fields
real, dimension(:,:), pointer :: &
co2_flux_out => NULL(), & ! daily co2 flux from LPJ to atmosphere [kg CO2/m2/sec]
q_ca_out => NULL(), & !
runoff_out => NULL(), & !
roughness_length_out => NULL(), & !
surface_temp_out => NULL(), & !
albedo_out => NULL(), & !
evap1_out => NULL(), & !
dedt_out => NULL(), & !
gc_out => NULL(), & !
exlpj_out => NULL()
! TODO: folia protective cover [%]
! TODO: LAI
! the local runoff_out maps are assembled into a global map
! TODO: possible optimisations:
! 1. each task stores only the coast cell information for those cells
! which pour out into ocean cells within its own compute domain.
! 2. instead of global runoff map, it might be sufficient
! to store in each task a data domain with a halo width corresponding
! to the maximum redirection distance used in subroutine
! prepare_coastal_runoff_routing(), currently 4.
! The FMS domains code has some flexibility in the halowidth,
! but it is yet unclear wether this really increases performance.
real, dimension(:,:), pointer :: discharge => NULL() ! local map of runoff _into_ the ocean
! list of coast cells, global grid indices for runoff source and destination
integer :: ncoastcells = 0
integer, dimension(:), pointer :: &
coast_src_i => NULL(), &
coast_src_j => NULL(), &
coast_dst_i => NULL(), &
coast_dst_j => NULL()
real, dimension(:,:), pointer :: inland_sinks => NULL()
end type lpj_wrapper_type
!
!! To enable strict floating point exceptions inside this module.
!! Taken from the example in the Intel Fortran Compiler Reference Manual.
!! Most probably will not work with any other compiler
integer(4) :: original_fpe_flags
integer(4) :: strict_fpe_flags = FPE_M_TRAP_UND + FPE_M_TRAP_OVF + FPE_M_TRAP_DIV0 &
+ FPE_M_TRAP_INV + FPE_M_ABRUPT_UND + FPE_M_ABRUPT_DMZ
!! Some Fortran-operatiosn apparently abuse floating point registers for fast (?) transfer
!! of non-floating-point data. This can trigger unwanted exceptions.
!! Lets try suppress some of that.
integer(4) :: loose_fpe_flags = FPE_M_TRAP_DIV0 + FPE_M_ABRUPT_UND + FPE_M_ABRUPT_DMZ
! ==== local module variables ================================================
type(lpj_wrapper_type),save :: lpj_state ! internal state of the model ! do we really need save here?
integer :: lpj_init_clock_id, lpj_update_fast_clock_id, lpj_update_slow_clock_id, lpj_update_cpl_clock_id
! wrapper function in fmsclocks.F90, to pull in the object file, so that it is available
! in the LPJ C interface
INTERFACE
FUNCTION FMSCLOCK_ID(NAME)
CHARACTER(*), INTENT(IN) :: NAME
INTEGER(KIND=4) :: FMSCLOCK_ID
END FUNCTION FMSCLOCK_ID
END INTERFACE
integer :: isphum ! number of specific humidity in the tracer array, or NO_TRACER
integer :: ico2 ! number of CO2 in the tracer array, or NO_TRACER
! cold start CO2 concentration [kg CO2 / kg wet air]
! 434.5e-06 [kg/kg] corresponds to 285.9 [ppmv] in dry air, i.e. q_ca = 0
real :: co2_init = 434.5e-06 ! [kg/kg wet air] also used as constant CO2 concentration if no CO2 tracer is defined in the field table
! Taken from land_lad2/canopy/cana_tile.F90
! The value 10.0 is set in the namelist of the ESM2M_pi-control_C2 example setup
real :: canopy_air_mass_for_tracers = 10.0 ! mass of wet air in the canopy air
! space for tracers other than water vapor, [kg wet air/m^2]
! Number of local and global grid cells in each dimensions.
! We require that FMS/LAD and LPJ run on the same grid resolution.
integer :: nlons, nlats, nlons_global, nlats_global
integer :: is, ie, js, je ! local domain boundary indices
logical :: &
do_lpj = .true., & ! invoke LPJ at all?
do_check_land_masks = .true., & ! set to false for production runs
use_co2_from_lpj = .false., &
use_q_ca_from_lpj = .false., &
use_runoff_from_lpj = .false., &
! use_runoff_antarctica_from_lpj is obsolete
use_runoff_antarctica_from_lpj = .false., &
use_roughness_length_from_lpj = .false., &
use_surface_temp_from_lpj = .false., &
use_albedo_from_lpj = .false., &
enable_fpe = .false.
real :: min_water_frac = 1e-4 ! min water area fraction for the point to be valid discharge destination
real :: min_land_frac = 1e-4 ! min land fraction for the cell to be considered
! "land" for discharge destination calculations
integer :: & ! diag manager field ids for input data on land_lad resolution
id_lad_mask, &
id_fprec, id_lprec, id_prec, &
id_temp, id_temp_mean, id_temp_min, id_temp_max, &
id_swdown, id_lwnet, id_t_flux, id_wind, id_qflux, id_dedq, id_dhdt, id_drdt,&
id_co2, id_fco2, id_Dfco2Dq, id_co2_ppmv
integer :: & ! diag manager field ids for output data, on LPJ units and daily time step
id_lpj_domainmask, id_lpj_mask, &
id_carbon_flux_out, id_co2_flux_out, id_q_ca_out, id_roughness_length_out, id_runoff_out, id_surface_temp_out, id_albedo_out, id_evap1_out, id_dedt_out, id_gc_out, id_exlpj_out
! ... and on land_lad units and time step
integer :: id_co2_out_lad, id_q_ca_out_lad, id_roughness_length_out_lad, id_runoff_out_lad, id_surface_temp_out_lad, id_albedo_out_lad, id_discharge_out_lad, id_discharge_out_lad_aw
!
!
! Override land_lad-produced data with LPJ-produced data
!
!
! Override land_lad-produced data with LPJ-produced data
!
!
! Override land_lad-produced data with LPJ-produced data
!
!
! obsolete, just to keep old namelists valid.
!
!
! Override land_lad-produced data with LPJ-produced data
!
!
! Override land_lad-produced data with LPJ-produced data
!
!
! Override land_lad-produced data with LPJ-produced data
!
!
! Mass of wet air in the canopy air space for tracers other than water vapor.
! The default value 10.0 is taken from the ESM2M_pi-control_C2 example setup input.nml
!
!
! Initial CO2 concentration.
! Unfortunately, this is not passed from the atmosphere.
! 434.5e-06 [kg/kg] corresponds to 285.9 [ppmv] in dry air.
!
!
! Minimum water area fraction for the point to be a
! valid discharge destination
!
!
! "Land" for discharge destination calculations
!
namelist /lpj_driver_nml/ &
do_lpj, &
use_co2_from_lpj, &
use_q_ca_from_lpj, &
use_runoff_from_lpj, &
use_runoff_antarctica_from_lpj, &
use_roughness_length_from_lpj, &
use_surface_temp_from_lpj, &
use_albedo_from_lpj, &
canopy_air_mass_for_tracers, &
co2_init, &
enable_fpe, &
min_water_frac, &
min_land_frac
!
contains
! ============================================================================
!
!
! Initializes the LPJ driver
!
!
! ...
! Sigh. land_lad uses two very similar named variables,
! n_tiles and ntiles,
! but to my understanding they are indeed very independant from
! each other. One seems to reflect the grid type (1 = spherical lat/lon,
! 6 = so-called cubed sphere, else unknown).
! The other has yet unknown meaning.
! Here we require that both must be 1.
!
!
subroutine lpj_driver_init &
( time_init, time, dt_fast, dt_slow, &
lad_domain, lad_layout, &
glon, glat, &
n_tiles, ntiles, &
lad_id_lon, lad_id_lat, &
lad_blon, lad_blat, lad_garea, lad_gfrac, isphum_in, ico2_in, &
bnd)
!USE, INTRINSIC :: IEEE_EXCEPTIONS
!type(atmos_land_boundary_type), intent(inout) :: atmos2land ! land boundary data
type(time_type), intent(in) :: time_init ! initial time of simulation (?)
type(time_type), intent(in) :: time ! current time
type(time_type), intent(in) :: dt_fast ! fast time step
type(time_type), intent(in) :: dt_slow ! slow time step
type(domain2d), intent(in) :: lad_domain ! domain of computations
integer, dimension(1:2), intent(in) :: lad_layout ! PE layout of lad domain
real, intent(in) :: glon(:,:), glat(:,:) ! global cell centers, used for mapping between LPJ and FMS domains
integer, intent(in) :: n_tiles, ntiles ! sigh. am I confused, or was the land_lad programmer? S.P.
integer, intent(in) :: lad_id_lon, lad_id_lat ! diag_manager axis ids for land_lad coordinates
real, intent(in) :: lad_blon(:,:) ! lon corners of the local grid cells in radians
real, intent(in) :: lad_blat(:,:) ! lat corners of the local grid cells in radians
real, intent(in) :: lad_garea(:,:) ! global cell area
real, intent(in) :: lad_gfrac(:,:) ! global map of fraction of land in the grid cell, 0 for ocean - or should we use land area instead?
integer, intent(in) :: isphum_in, ico2_in ! index in tracer array for specific humidity and CO2, or NO_TRACER
type(land_data_type), intent(in) :: bnd ! land boundary data, to share the initial conditions from land_lad
!
integer :: i,j,k, unit, ierr, io, used
! some magical values for testing the interface between Fortran90 and C++ code.
! we pass them from Fortran90 to C++ functions and test in there if they were
! received correctly.
integer, parameter :: deadbeef = X'eadbeef' ! X'deadbeef' yields overflow error on 32bit machines
integer, parameter :: fse = 4711
real, parameter :: onetwothree = 1234567809.0987654321D0, &
pitest = PI ! from constants_mod ! 3.1415926535897932384626433832795029D0 ! from linux
logical :: true = .true., false = .false.
integer :: kind_i = KIND(deadbeef), kind_r = KIND(pitest), kind_l = KIND(true)
real, dimension(2,3,4) :: testarray3d
integer :: atmos_comm ! the MPI communicator for the atmosphere
integer,allocatable,dimension(:) :: lad_pelist ! needed to obtain the peset
integer :: Time_init_year, Time_init_month, Time_init_day, Time_init_hour, Time_init_minute, Time_init_second
integer :: Time_year, Time_month, Time_day, Time_hour, Time_minute, Time_second
!integer :: Time_step_days, Time_step_seconds
integer :: log_unit
integer :: pe, npes
integer :: lad_nlons, lad_nlats ! number of cells in lad resolution
real, dimension(:), allocatable :: lad_lons, lad_lats ! lad cell center coordinates
! info for assembling the global runoff routing map
integer :: global_number_lpj_cells = 0 ! number of LPJ cells (land cells)
real, allocatable, dimension(:) :: lpj_cell_lons, lpj_cell_lats
integer, allocatable, dimension(:) :: lpj_cell_next
integer :: id_runoff_map = -1, id_runoff_coastcells = -1
integer :: id_runoff_dest_lon = -1, id_runoff_dest_lat = -1
integer :: id_runoff_dcount = -1, id_runoff_basin = -1
integer :: id_runin_coastcells = -1
integer :: id_inland_sinks = -1 ! map of inland sink cells
!! Water content of all rivers is written to lpj_runoff.
!! Water content of inland sinks can be computed by
!! applying inland_sinks as mask to lpj_runoff.
integer :: id_warea = -1
! just for diagnostic output, as helper to prepare runoff routing map
integer :: id_horo = -1
real, allocatable, dimension(:,:) :: HORO ! Orography height, averaged over grid cell
integer :: id_lad_domainmask
real, allocatable, dimension(:,:) :: lad_domainmask
integer(4) :: previous_fpe_flags
!CALL IEEE_SET_HALTING_MODE(IEEE_DIVIDE_BY_ZERO, .FALSE.)
if (module_is_initialized) return
!----- write version and namelist to log file -----
call write_version_number ( version, tagname )
write(*,*) 'lpj_driver_init', version
unit = open_namelist_file()
ierr=1;
do while (ierr /= 0)
read (unit, nml=lpj_driver_nml, iostat=io, end=10)
ierr = check_nml_error(io,'lpj_driver_nml') ! also initializes nml error codes
enddo
10 call close_file (unit)
if ( mpp_pe() == mpp_root_pe() ) then
log_unit = stdlog()
write (log_unit, nml=lpj_driver_nml)
write (*, nml=lpj_driver_nml)
endif
! for preparation of LPJ configuration we need the topography map on
! LPJ resolution and land-sea-mask even before we can actually start and use LPJ.
! Thus write it out here, before checking for do_lpj .
call mpp_get_compute_domain(lad_domain, is, ie, js, je)
id_horo = register_static_field('lpj', "horo", (/lad_id_lon,lad_id_lat/), "Orographic Height from FMS", "m")
if (id_horo > 0) then
allocate(HORO(ie-is+1, je-js+1))
write(*,*) 'for get_topog_mean(): shape(HORO)', shape(HORO)
write(*,*) ' vs. (/size(lad_blon,1)-1,size(lad_blon,2)-1/): ', size(lad_blon,1)-1, size(lad_blon,2)-1
used = get_topog_mean(lad_blon, lad_blat, HORO)
if (.not. used) &
call error_mesg (module_name, 'could not get topography mean HORO', FATAL)
used = send_data(id_horo, HORO, Time)
deallocate(HORO)
end if
if (.not. do_lpj) then
write (*,*) 'do_lpj = .false. -- Not invoking LPJ!'
return
endif
if (enable_fpe) then
!! To enable strict floating point exceptions inside this module.
!! Taken from the example in the Intel Fortran Compiler Reference Manual.
!! Most probably will not work with any other compiler
call clearstatusfpqq() ! do not inherit signalled FP exceptions from FMS
original_fpe_flags = FOR_SET_FPE(strict_fpe_flags)
endif
lpj_init_clock_id = mpp_clock_id( 'LPJ_init', flags=clock_flag_default, grain=CLOCK_SUBCOMPONENT )
lpj_update_slow_clock_id = mpp_clock_id( 'LPJ_update_slow', flags=clock_flag_default, grain=CLOCK_SUBCOMPONENT )
lpj_update_fast_clock_id = mpp_clock_id( 'LPJ_update_fast', flags=clock_flag_default, grain=CLOCK_SUBCOMPONENT )
lpj_update_cpl_clock_id = fmsclock_id( 'LPJ_update_cpl' ) ! to force the linker to load the fmsclock wrapper object file
call mpp_clock_begin(lpj_init_clock_id)
if (n_tiles /= 1) then
call error_mesg('lpj_driver_init', 'n_tiles /= 1', FATAL)
endif
if (ntiles /= 1) then
call error_mesg('lpj_driver_init', 'ntiles /= 1', FATAL)
endif
if (.not. is_latlon(lad_blon, lad_blat)) then
call error_mesg('lpj_driver_init', 'FMS grid is not lat-lon type', FATAL)
endif
! TODO: collect LPJ parameters from namelist and pass them to C-side of LPJ driver
pe = mpp_pe()
npes = mpp_npes()
allocate(lad_pelist(npes))
call mpp_get_pelist(lad_domain, lad_pelist)
call mpp_get_current_pelist(lad_pelist, commID=atmos_comm)
deallocate(lad_pelist)
lpj_state%domain = lad_domain
lad_nlons = size(lad_blon, 1)-1 ! local domain
lad_nlats = size(lad_blat, 2)-1
!write(*,*) 'lad_blon lbound', lbound(lad_blon), 'ubound', ubound(lad_blon)
!write(*,*) 'lad_blat lbound', lbound(lad_blat), 'ubound', ubound(lad_blat)
write(*,'(a,i4,a,i4,a,i4,a,i4,a,i5,a,i5)') 'is ', is, ' ie ', ie, ' js ', js, ' je ', je, ' lad_nlons ', lad_nlons, ' lad_nlats ', lad_nlats
if (ie-is+1 /= lad_nlons) &
call error_mesg('lpj_driver_init', 'dimension of land_lad lon boundaries does not match boundary indices', FATAL)
if (je-js+1 /= lad_nlats) &
call error_mesg('lpj_driver_init', 'dimension of land_lad lat boundaries does not match boundary indices', FATAL)
write(*,*) 'lad_gfrac bound ',lbound(lad_gfrac),ubound(lad_gfrac)
allocate(lad_lons(is:ie), lad_lats(js:je))
! sigh. The array boundary values are lost during the subroutine call.
! Thus, lad_blon and lad_blat always start at 1 here, even though with domain decomposition
! they were allocated with different boundaries by our caller.
lad_lons = ((lad_blon(1:lad_nlons,1)+lad_blon(2:lad_nlons+1,1))/2)*RAD_TO_DEG ! lad cell center longitudes in degrees
lad_lats = ((lad_blat(1,1:lad_nlats)+lad_blat(1,2:lad_nlats+1))/2)*RAD_TO_DEG ! lad cell center latitudes in degrees
lpj_state%current_time = Time
lpj_state%dt_fast = dt_fast
lpj_state%dt_slow = dt_slow
call get_date(Time_init, Time_init_year, Time_init_month, Time_init_day, &
Time_init_hour, Time_init_minute, Time_init_second)
call get_date(Time, Time_year, Time_month, Time_day, &
Time_hour, Time_minute, Time_second)
call get_time(dt_fast, lpj_state%dt_fast_seconds)
!lpj_state%dt_fast_seconds = Time_step_seconds
call get_time(dt_slow, lpj_state%dt_slow_seconds)
!lpj_state%dt_slow_seconds = Time_step_seconds
allocate( & ! data fields
!lpj_state%lprec(is:ie,js:je), &
!lpj_state%fprec(is:ie,js:je), &
lpj_state%prec(is:ie,js:je), &
lpj_state%temp_mean(is:ie,js:je), &
lpj_state%temp_max(is:ie,js:je), &
lpj_state%temp_min(is:ie,js:je), &
lpj_state%swdown(is:ie,js:je), &
lpj_state%lwnet(is:ie,js:je), &
lpj_state%t_flux(is:ie,js:je), &
lpj_state%qflux(is:ie,js:je), &
lpj_state%dedq(is:ie,js:je), &
lpj_state%dhdt(is:ie,js:je), &
lpj_state%drdt(is:ie,js:je), &
lpj_state%drag_q(is:ie,js:je), &
lpj_state%p_surf(is:ie,js:je), &
lpj_state%wind(is:ie,js:je), &
lpj_state%co2(is:ie,js:je), &
!lpj_state%fco2(is:ie,js:je), &
!lpj_state%Dfco2Dq(is:ie,js:je) &
STAT=ierr)
if (0 /= ierr) &
call error_mesg('lpj_driver_init', 'allocate input fields failed', FATAL)
allocate( &
lpj_state%co2_flux_out(is:ie,js:je), &
lpj_state%q_ca_out(is:ie,js:je), &
lpj_state%runoff_out(is:ie,js:je), &
lpj_state%roughness_length_out(is:ie,js:je), &
lpj_state%surface_temp_out(is:ie,js:je), &
lpj_state%albedo_out(is:ie,js:je), &
lpj_state%evap1_out(is:ie,js:je), &
lpj_state%dedt_out(is:ie,js:je), &
lpj_state%gc_out(is:ie,js:je), &
lpj_state%exlpj_out(is:ie,js:je), &
STAT=ierr)
if (0 /= ierr) &
call error_mesg('lpj_driver_init', 'allocate output fields failed', FATAL)
! TODO: check which of the lpj_state variables need to be saved in RESTART files
! Make sure that we have nice values also over ocean and antarctica!
! else netcdf output will complain bitterly.
lpj_state%co2_flux_out = 0
lpj_state%q_ca_out = 0
lpj_state%runoff_out = 0
lpj_state%roughness_length_out = bnd%rough_mom(:,:,1)
!lpj_state%surface_temp_out = 288.0 - KELVIN ! default initial condition from soil.F90 ???
lpj_state%surface_temp_out = bnd%t_surf(:,:,1) ! copy the initial conditions that were set in soil.F90
lpj_state%albedo_out = 0
lpj_state%evap1_out = 0
lpj_state%dedt_out = 0
lpj_state%gc_out = 0
lpj_state%exlpj_out = 0
allocate(lpj_state%lad_land_mask(is:ie, js:je))
lpj_state%lad_land_mask = lad_gfrac(is:ie,js:je)
allocate(lpj_state%lad_warea(is:ie, js:je))
lpj_state%lad_warea = lad_garea(is:ie,js:je) * (1-lad_gfrac(is:ie,js:je))
! now prepare to call out into the C - implemented heart of LPJ
! first assemble some parameters with well-known values to verify
! that the F90<->C++ interface is OK.
do i=1,2
do j=1,3
do k=1,4
testarray3d(i,j,k) = i*100.0+j+k/10.0
end do
end do
end do
! The arrays in lpj_state%XXX are allocate()d, they might behave different than
! testarray3d. We put remarkable values at the corners of such array to test
! that we can find those values on the C side in the expected locations.
lpj_state%co2(is,js) = 42.42
lpj_state%co2(ie,js) = 43.43
lpj_state%co2(is,je) = 44.44
lpj_state%co2(ie,je) = 45.45
!write(*,*) 'co2 test array', lpj_state%co2
write(*,*) 'MAXEXPONENT ', MAXEXPONENT(pitest)
! Get back from LPJ information about its grid resolution.
! TODO: Check if it matches the land_lad grid resolution.
! We pass the boundary box (cell borders) of the FMS compute domain to LPJ.
! On the LPJ side we then can compute a grid in LPJ-resolution inside those boxes,
! and can compute the mappings between LPJ grid cell indices and FMS lat-lon-array indices.
! This slightly complicated approach is necessary because in LPJ only the land cells
! are represented, while the boundary box of the FMS domain might also contain ocean-only
! cells along its borders.
call lpj_init( &
! first some arguments for checking the F90->C++ interface
deadbeef, pitest, &
true, false, &
testarray3d, &
lpj_state%co2, &
kind_i, kind_r, kind_l, &
! now the ``real'' arguments
! the FMS time stamps are broken out into days and seconds
Time_init_year, Time_init_month, Time_init_day, &
Time_init_hour, Time_init_minute, Time_init_second, &
Time_year, Time_month, Time_day, &
Time_hour, Time_minute, Time_second, &
lpj_state%dt_fast_seconds, lpj_state%dt_slow_seconds, &
! information about domain partitioning
atmos_comm, pe, mpp_root_pe(), npes, &
lad_layout, & ! of domain decomposition, number of PEs in lon and lat directions
! cell centers of corners of global grid
glon(1,1)*RAD_TO_DEG, glon(size(glon,1),size(glon,2))*RAD_TO_DEG, &
glat(1,1)*RAD_TO_DEG, glat(size(glat,1),size(glat,2))*RAD_TO_DEG, &
! the boundary indices of the local compute domain
is, ie, js, je, &
! grid axis info of local domain
lad_nlons, lad_nlats, &
!!blon, blat, & pass only 1-D arrays here because we require rectilinear lat-lon grid
lad_blon(1,1)*RAD_TO_DEG, lad_blon(lad_nlons+1,lad_nlats+1)*RAD_TO_DEG, &
lad_blat(1,1)*RAD_TO_DEG, lad_blat(lad_nlons+1,lad_nlats+1)*RAD_TO_DEG, & ! boundary box of FMS compute domain
!lons, lats, &
lpj_update_cpl_clock_id, &
! LPJ-specific options TODO
1, & !argc,
!argv, &
! output parameters
global_number_lpj_cells, &
! finally interface safety parameters.
onetwothree, fse &
)
! now, after the co2 array has been (ab)used for the F90->C interface test,
! we can set the initial value for the CO2 concentration
lpj_state%co2 = co2_init ! [kg/kg wet air]
isphum = isphum_in
ico2 = ico2_in
if (ico2 /= NO_TRACER) bnd%tr(:,:,1,ico2) = co2_init ! remains [kg/kg wet air]
! get information to set up the horizontal interpolation
! (not really needed until we allow different resolution between land_lad and LPJ)
nlons = lad_nlons
nlats = lad_nlats
nlons_global = size(glon,1)
nlats_global = size(glon,2)
write(*,'(a,i4,a,i4,a,i4,a,i4)') 'nlons ',nlons,' nlats ',nlats, &
' nlons_global ',nlons_global,' nlats_global ',nlats_global
allocate(lpj_state%lpj_land_mask(is:ie,js:je))
allocate(lpj_state%discharge(is:ie,js:je), STAT=ierr)
if (0 /= ierr) &
call error_mesg('lpj_driver_init', 'allocate discharge failed', FATAL)
lpj_state%discharge = 0.0
! The arrays in lpj_state%XXX are allocate()d, they might behave different than
! testarray3d. We put remarkable values at the corners of such array to test
! that we can find those values on the C side in the expected locations.
lpj_state%lpj_land_mask(is,js) = 4242.42
lpj_state%lpj_land_mask(ie,js) = 4343.43
lpj_state%lpj_land_mask(is,je) = 4444.44
lpj_state%lpj_land_mask(ie,je) = 4545.45
allocate(lpj_cell_lons(global_number_lpj_cells))
allocate(lpj_cell_lats(global_number_lpj_cells))
allocate(lpj_cell_next(global_number_lpj_cells))
lpj_cell_lons = 0;
lpj_cell_lats = 0;
lpj_cell_next = 0;
!allocate(lpj_state%global_gfrac(nlons_global, nlats_global))
!lpj_state%global_gfrac = lad_gfrac
call lpj_init_grid_interpolation( &
! output parameters
lpj_state%lpj_land_mask, &
lpj_cell_lons, &
lpj_cell_lats, &
lpj_cell_next, &
! interface safety parameters
fse &
)
!if (do_check_land_masks) & ! here do it unconditionally
call check_land_masks(lpj_state%lad_land_mask, lpj_state%lpj_land_mask, 'LPJ', 'land_lad', 'LPJ')
id_lad_mask = register_static_field('lpj', 'lad_land_mask', &
(/lad_id_lon,lad_id_lat/), &
long_name='land fraction as obtained from land_lad (gfrac)', &
missing_value = 0.0 &
)
id_lad_domainmask = register_static_field('lpj', 'lad_domain_mask', &
(/lad_id_lon,lad_id_lat/), long_name='MPI task numbers of lad cells', &
units="lad task number", missing_value = -1.0 &
)
id_lpj_domainmask = register_static_field('lpj', 'lpj_domain_mask', &
(/lad_id_lon,lad_id_lat/), long_name='MPI task numbers of LPJ cells', &
units="lpj task number", missing_value = -1.0 &
)
id_lpj_mask = register_static_field('lpj', 'lpj_land_mask', &
(/lad_id_lon,lad_id_lat/), long_name='land mask as obtained from LPJ', &
missing_value = -1.0 &
)
id_warea = register_static_field('lpj', 'lad_warea', &
(/lad_id_lon,lad_id_lat/), long_name='ocean area per grid cell', &
units = 'm^2', missing_value = -1.0 &
)
! Note that the static fields are actually written out into the file only in
! diag_manager_end()
if (id_lad_mask > 0) &
used = send_data(id_lad_mask, lpj_state%lad_land_mask)
if (id_lad_domainmask > 0) then
allocate(lad_domainmask(is:ie,js:je))
lad_domainmask(:,:) = pe
used = send_data(id_lad_domainmask, lad_domainmask)
deallocate(lad_domainmask)
end if
if (id_lpj_domainmask > 0) &
used = send_data(id_lpj_domainmask, lpj_state%lpj_land_mask)
if (id_lpj_mask > 0) then
where (lpj_state%lpj_land_mask > 0)
lpj_state%lpj_land_mask = 1
elsewhere
lpj_state%lpj_land_mask = 0
end where
used = send_data(id_lpj_mask, lpj_state%lpj_land_mask)
end if
if (id_warea > 0) &
used = send_data(id_warea, lpj_state%lad_warea)
id_runoff_map = register_static_field('lpj', 'lpj_runoffmap', &
(/lad_id_lon,lad_id_lat/), &
long_name='runoff routing map (destination cell index)', &
missing_value = -10.0 & ! FMS ocean
)
id_runoff_dest_lon = register_static_field('lpj', 'lpj_runoff_dest_lon', &
(/lad_id_lon,lad_id_lat/), &
long_name='runoff routing map (destination cell longitude)', &
units = 'degrees_E', &
missing_value = -1000.0 & ! FMS ocean
)
id_runoff_dest_lat = register_static_field('lpj', 'lpj_runoff_dest_lat', &
(/lad_id_lon,lad_id_lat/), &
long_name='runoff routing map (destination cell latitude)', &
units = 'degrees_N', &
missing_value = -1000.0 & ! FMS ocean
)
id_runoff_coastcells = register_static_field('lpj', 'lpj_coastcells', &
(/lad_id_lon,lad_id_lat/), &
long_name='land coastcells which route into ocean cells', &
missing_value = 0.0 &
)
id_runoff_dcount = register_static_field('lpj', 'lpj_dcount', &
(/lad_id_lon,lad_id_lat/), &
long_name='count of cells which route into this destination point', &
missing_value = 0.0 &
)
id_runoff_basin = register_static_field('lpj', 'lpj_basin', &
(/lad_id_lon,lad_id_lat/), &
long_name='river basins', &
missing_value = 0.0 &
)
id_runin_coastcells = register_static_field('lpj', 'lpj_beachcells', &
(/lad_id_lon,lad_id_lat/), &
long_name='ocean coastcells showing count of cells which route into this destination point', &
missing_value = -1.0 &
)
id_inland_sinks = register_static_field('lpj', 'lpj_inland_sinks', &
(/lad_id_lon,lad_id_lat/), &
long_name='inland river sinks', &
missing_value = 0.0 &
)
call prepare_coastal_runoff_routing(global_number_lpj_cells, &
lpj_cell_lons, lpj_cell_lats, lpj_cell_next, &
lad_gfrac, &
id_runoff_map, id_runoff_coastcells, &
id_runoff_dest_lon, id_runoff_dest_lat, &
id_runoff_dcount, id_runoff_basin, &
id_runin_coastcells, id_inland_sinks)
deallocate(lpj_cell_lons)
deallocate(lpj_cell_lats)
deallocate(lpj_cell_next)
!!! to get the above-written static fields out into the file now,
!!! and not only after the entire program ended successfully.
!call diag_manager_end(Time) !!!
!flush(stdout())
!! init diagnostic output
! input data on land_lad grid
id_prec = register_diag_field('lpj', 'prec', (/lad_id_lon,lad_id_lat/), Time, &
long_name='daily sum of liquid and frozen precipitation', &
!standard_name='', &
units='mm', missing_value=-999.0 &
)
id_temp_mean = register_diag_field('lpj', 'temp_mean', (/lad_id_lon,lad_id_lat/), Time, &
long_name='daily mean temperature', &
!standard_name='', &
units='degC', missing_value=-999.0 &
)
id_temp_max = register_diag_field('lpj', 'temp_max', (/lad_id_lon,lad_id_lat/), Time, &
long_name='daily maximum temperature', &
!standard_name='', &
units='degC', missing_value=-999.0 &
)
id_temp_min = register_diag_field('lpj', 'temp_min', (/lad_id_lon,lad_id_lat/), Time, &
long_name='daily minimum temperature', &
units='degC', missing_value=-999.0 &
)
id_swdown = register_diag_field('lpj', 'swdown', (/lad_id_lon,lad_id_lat/), Time, &
long_name='daily mean shortwave downwelling radiation', &
!standard_name='', &
units='W/m2', missing_value=-999.0 &
)
id_lwnet = register_diag_field('lpj', 'lwnet', (/lad_id_lon,lad_id_lat/), Time, &
long_name='daily mean longwave net radiation', &
!standard_name='', &
units='W/m2', missing_value=-999.0 &
)
id_t_flux = register_diag_field('lpj', 't_flux', (/lad_id_lon,lad_id_lat/), Time, &
long_name='daily mean sensible heat flux', &
!standard_name='', &
units='W/m2', missing_value=-999.0 &
)
id_qflux = register_diag_field('lpj', 'qflux', (/lad_id_lon,lad_id_lat/), Time, &
long_name='water vapor flux near surface', &
!standard_name='', &
units='kg/m2/s', missing_value=-999.0 &
)
id_dedq = register_diag_field('lpj', 'dedq', (/lad_id_lon,lad_id_lat/), Time, &
long_name='Moisture flux humidity sensitivity', &
!standard_name='', &
units='kg/m2/s', missing_value=-999.0 &
)
id_dhdt = register_diag_field('lpj', 'dhdt', (/lad_id_lon,lad_id_lat/), Time, &
long_name='derivative of sens over T', &
!standard_name='', &
units='kg/m2/s', missing_value=-999.0 &
)
id_drdt = register_diag_field('lpj', 'drdt', (/lad_id_lon,lad_id_lat/), Time, &
long_name='derivative of LW radiation over T', &
!standard_name='', &
units='kg/m2/s', missing_value=-999.0 &
)
id_wind = register_diag_field('lpj', 'wind', (/lad_id_lon,lad_id_lat/), Time, &
long_name='daily mean absolute wind speed', &
!standard_name='', &
units='m/s', missing_value=-999.0 &
)
id_fco2 = register_diag_field('lpj', 'fco2', (/lad_id_lon,lad_id_lat/), Time, &
long_name='daily mean atmospheric CO2 flux', &
!standard_name='', &
units='kg CO2/(m2 s)', missing_value=-999.0 &
)
id_Dfco2Dq = register_diag_field('lpj', 'Dfco2Dq', (/lad_id_lon,lad_id_lat/), Time, &
long_name='daily mean atmospheric CO2 d(tracer flux)/d(surf tracer)', &
!standard_name='', &
units='kg wet air/(m2 s)', missing_value=-999.0 &
)
id_co2 = register_diag_field('lpj', 'co2', (/lad_id_lon,lad_id_lat/), Time, &
long_name='daily integrated atmospheric CO2 concentration', &
!standard_name='', &
units='kg CO2/kg wet air', missing_value=-999.0 &
)
id_co2_ppmv = register_diag_field('lpj', 'co2_ppmv', (/lad_id_lon,lad_id_lat/), Time, &
long_name='daily inegrated atmospheric CO2 concentration', &
!standard_name='', &
units='ppmv', missing_value=-999.0 &
)
! daily output data in LPJ units
id_carbon_flux_out = register_diag_field('lpj', 'carbon_flux_out', (/lad_id_lon,lad_id_lat/), Time, &
long_name='LPJ daily carbon flux to atmosphere output', &
!standard_name='', &
units='g C/(m^2 day)', missing_value=-999.0 &
)
id_co2_flux_out = register_diag_field('lpj', 'co2_flux_out', (/lad_id_lon,lad_id_lat/), Time, &
long_name='LPJ daily CO2 flux to atmosphere output', &
!standard_name='', &
units='kg CO2/(m^2 sec)', missing_value=-999.0 &
)
id_q_ca_out = register_diag_field('lpj', 'q_ca_out', (/lad_id_lon,lad_id_lat/), Time, &
long_name='surface/canopy humidity', &
!standard_name='', &
units='kg/kg', missing_value=-999.0 &
)
id_roughness_length_out = register_diag_field('lpj', 'roughness_length_out', (/lad_id_lon,lad_id_lat/), Time, &
long_name='LPJ roughness length output', &
!standard_name='', &
units='??', missing_value=-999.0 &
)
id_runoff_out = register_diag_field('lpj', 'runoff_out', (/lad_id_lon,lad_id_lat/), Time, &
long_name='LPJ runoff', &
!standard_name='', &
units='kg/s', missing_value=-999.0 &
)
id_surface_temp_out = register_diag_field('lpj', 'surface_temp_out', (/lad_id_lon,lad_id_lat/), Time, &
long_name='LPJ surface temperature output', &
!standard_name='', &
units='degC', missing_value=-999.0 &
)
id_albedo_out = register_diag_field('lpj', 'albedo_out', (/lad_id_lon,lad_id_lat/), Time, &
long_name='LPJ albedo output', &
!standard_name='', &
units='1', missing_value=-999.0 &
)
id_evap1_out = register_diag_field('lpj', 'evap1_out', (/lad_id_lon,lad_id_lat/), Time, &
long_name='LPJ evapotranspiration output', &
!standard_name='', &
units='kg/m2/s', missing_value=-999.0 &
)
id_dedt_out = register_diag_field('lpj', 'dedt_out', (/lad_id_lon,lad_id_lat/), Time, &
long_name='LPJ evapotranspiration output', &
!standard_name='', &
units='kg/m2/s', missing_value=-999.0 &
)
id_gc_out = register_diag_field('lpj', 'gc_out', (/lad_id_lon,lad_id_lat/), Time, &
long_name='LPJ evapotranspiration output', &
!standard_name='', &
units='kg/m2/s', missing_value=-999.0 &
)
id_exlpj_out = register_diag_field('lpj', 'exlpj_out', (/lad_id_lon,lad_id_lat/), Time, &
long_name='LPJ evapotranspiration output', &
!standard_name='', &
units='kg/m2/s', missing_value=-999.0 &
)
! output data interpolated onto sub-daily time steps and converted to land_lad units
! should be output into a different file than above
id_co2_out_lad = register_diag_field('lpj_lad', 'co2_out', (/lad_id_lon,lad_id_lat/), Time, &
long_name='LPJ CO2 output on land_lad resolution', &
!standard_name='', &
units='kg/kg wet air', missing_value=-999.0 &
)
id_q_ca_out_lad = register_diag_field('lpj_lad', 'q_ca_out', (/lad_id_lon,lad_id_lat/), Time, &
long_name='surface/canopy humidity', &
!standard_name='', &
units='kg/kg', missing_value=-999.0 &
)
id_roughness_length_out_lad = register_diag_field('lpj_lad', 'roughness_length_out', (/lad_id_lon,lad_id_lat/), Time, &
long_name='LPJ roughness length output on land_lad resolution', &
!standard_name='', &
units='??', missing_value=-999.0 &
)
id_runoff_out_lad = register_diag_field('lpj_lad', 'runoff_out', (/lad_id_lon,lad_id_lat/), Time, &
long_name='LPJ runoff on land_lad resolution', &
!standard_name='', &
units='kg/s', missing_value=-999.0 &
)
id_surface_temp_out_lad = register_diag_field('lpj_lad', 'surface_temp_out', (/lad_id_lon,lad_id_lat/), Time, &
long_name='LPJ surface temperature output on land_lad resolution', &
!standard_name='', &
units='K', missing_value=-999.0 &
)
id_albedo_out_lad = register_diag_field('lpj_lad', 'albedo_out', (/lad_id_lon,lad_id_lat/), Time, &
long_name='LPJ albedo output on land_lad resolution', &
!standard_name='', &
units='1', missing_value=-999.0 &
)
id_discharge_out_lad = register_diag_field('lpj_lad', 'discharge_out', (/lad_id_lon,lad_id_lat/), Time, &
long_name='LPJ runoff discharged into the ocean on land_lad resolution', &
!standard_name='', &
units='kg/s', missing_value=-999.0 &
)
id_discharge_out_lad_aw = register_diag_field('lpj_lad', 'discharge_out_aw', (/lad_id_lon,lad_id_lat/), Time, &
long_name='LPJ area-weighted runoff discharged into the ocean on land_lad resolution', &
!standard_name='', &
units='kg/m2/s', missing_value=-999.0 &
)
deallocate(lad_lons, lad_lats)
!{ begin init_output
!
! Init daily output from values that are saved in the LPJ restart file.
! We must pass back meaningful values to land_lad/FMS already in the first time step,
! even before the invocation of lpj_update.
! This can be done only after all else has been initialised above.
! Note that the same code for writing out diagnostics and units conversion
! appears in lpj_update below.
call lpj_init_output( &
! get back output from LPJ
lpj_state%co2_flux_out, & ! returned as carbon flux in [g C/ (m^2/day)], must be converted to [kg CO2/(m^2 sec)]
lpj_state%q_ca_out, &
lpj_state%runoff_out, & ! returned in [m3/s]
lpj_state%roughness_length_out, &
lpj_state%surface_temp_out, & ! returned in [degC]
lpj_state%albedo_out, &
lpj_state%evap1_out, &
lpj_state%dedt_out, &
lpj_state%gc_out, &
lpj_state%exlpj_out, &
! input from FMS to LPJ
lpj_state%inland_sinks, &
! size of the arrays
size(lpj_state%q_ca_out, 1), &
size(lpj_state%q_ca_out, 2), &
! interface safety parameters
fse &
)
deallocate(lpj_state%inland_sinks)
lpj_state%inland_sinks => NULL()
! experimental, to prevent instabilities
where (lpj_state%q_ca_out>0.4)
lpj_state%q_ca_out=0.4
endwhere
where (lpj_state%q_ca_out<0)
lpj_state%q_ca_out=0
endwhere
write(*,'(a,i4,i4,a,i3,i3,i4,i4)') 'lpj_init_output returned size q_ca_out ', size(lpj_state%q_ca_out,1), size(lpj_state%q_ca_out,2),' bounds q_ca_out ', lbound(lpj_state%q_ca_out,1), lbound(lpj_state%q_ca_out,2), ubound(lpj_state%q_ca_out,1), ubound(lpj_state%q_ca_out,2)
write(*,'(a,f8.3,a,f8.3,a,f8.3,a,Z20)') ' q_ca_out min ',minval(lpj_state%q_ca_out(:,:)),' max ',maxval(lpj_state%q_ca_out(:,:)),' mean ',sum(lpj_state%q_ca_out(:,:))/size(lpj_state%q_ca_out(:,:)), ' checksum', mpp_chksum(lpj_state%q_ca_out(:,:))
call pour_into_ocean ! for runoff -> discharge
!
! LPJ output
!
!if (id_carbon_flux_out > 0) & ! [g Carbon/(m^2*day)]
! used = send_data(id_carbon_flux_out, lpj_state%co2_flux_out, lpj_state%current_time, &
! rmask=lpj_state%lpj_land_mask)
write(*,*) 'sending to id_q_ca_out ', id_q_ca_out, ' at ', Time_year, Time_month, Time_day, Time_hour, Time_minute, Time_second
if (id_q_ca_out > 0) &
used = send_data(id_q_ca_out, lpj_state%q_ca_out, lpj_state%current_time &
!, rmask=lpj_state%lpj_land_mask &
)
!if (maxval(lpj_state%runoff_out) == 0.0) __NOTE__('lpj_init_output: runoff_out is 0')
if (id_runoff_out > 0) &
used = send_data(id_runoff_out, lpj_state%runoff_out, lpj_state%current_time)
!if (id_roughness_length_out > 0) &
! used = send_data(id_roughness_length_out, lpj_state%roughness_length_out, lpj_state%current_time, &
! rmask=lpj_state%lpj_land_mask)
!if (id_surface_temp_out > 0) &
! used = send_data(id_surface_temp_out, lpj_state%surface_temp_out, lpj_state%current_time, &
! rmask=lpj_state%lpj_land_mask)
!if (id_albedo_out > 0) &
! used = send_data(id_albedo_out, lpj_state%albedo_out, lpj_state%current_time, &
! rmask=lpj_state%lpj_land_mask)
! now convert to correct units for land_lad
!lpj_state%co2_flux_out = lpj_state%co2_flux_out / SECONDS_PER_DAY * WTMCO2/WTMC / 1000 ! from [g Carbon/(m^2*day)] to [kg CO2/(m^2*sec)]
!if (id_co2_flux_out > 0) &
! used = send_data(id_co2_flux_out, lpj_state%co2_flux_out, lpj_state%current_time, &
! rmask=lpj_state%lpj_land_mask)
lpj_state%q_ca_out = lpj_state%q_ca_out ! from daily sum [kg/m^2/day]==[liter/day] to rate [kg/m^2/sec]
write(*,'(a)') 'lpj_init_output converted to FMS units '
write(*,'(a,g11.3,a,g11.3,a,g11.3,a,Z20)') ' q_ca_out min ',minval(lpj_state%q_ca_out(:,:)),' max ',maxval(lpj_state%q_ca_out(:,:)),' mean ',sum(lpj_state%q_ca_out(:,:))/size(lpj_state%q_ca_out(:,:)), ' checksum', mpp_chksum(lpj_state%q_ca_out(:,:))
!lpj_state%runoff_out = lpj_state%runoff_out / SECONDS_PER_DAY ! from daily sum [kg/m^2/day] to rate [kg/m^2/sec]
! roughness_length_out
!lpj_state%surface_temp_out = lpj_state%surface_temp_out + KELVIN ! from [C] to [K]
! albedo_out
!} end init_output
!!! to get the above-written static fields out into the file now,
!!! and not only after the entire program ended successfully.
!call diag_manager_end(Time) !!!
module_is_initialized = .true.
if (enable_fpe) then
!! To enable strict floating point exceptions inside this module.
!! Taken from the example in the Intel Fortran Compiler Reference Manual.
!! Most probably will not work with any other compiler
original_fpe_flags = FOR_SET_FPE(original_fpe_flags)
endif
call mpp_clock_end(lpj_init_clock_id)
end subroutine lpj_driver_init
!
! ============================================================================
!
!
! clean-up for the LPJ driver
!
!
! ...
!
!
subroutine lpj_driver_end( &
)
if (.not. do_lpj) return
call lpj_end()
!deallocate(lpj_state%lprec)
!deallocate(lpj_state%fprec)
deallocate(lpj_state%prec)
deallocate(lpj_state%temp_mean)
deallocate(lpj_state%temp_max)
deallocate(lpj_state%temp_min)
deallocate(lpj_state%swdown)
deallocate(lpj_state%lwnet)
deallocate(lpj_state%t_flux)
deallocate(lpj_state%qflux)
deallocate(lpj_state%dedq)
deallocate(lpj_state%dhdt)
deallocate(lpj_state%drdt)
deallocate(lpj_state%drag_q)
deallocate(lpj_state%p_surf)
deallocate(lpj_state%wind)
deallocate(lpj_state%co2)
!deallocate(lpj_state%fco2)
!deallocate(lpj_state%Dfco2Dq)
deallocate(lpj_state%co2_flux_out)
deallocate(lpj_state%q_ca_out)
deallocate(lpj_state%runoff_out)
deallocate(lpj_state%roughness_length_out)
deallocate(lpj_state%surface_temp_out)
deallocate(lpj_state%albedo_out)
deallocate(lpj_state%evap1_out)
deallocate(lpj_state%dedt_out)
deallocate(lpj_state%gc_out)
deallocate(lpj_state%exlpj_out)
deallocate(lpj_state%discharge)
deallocate(lpj_state%lad_land_mask)
deallocate(lpj_state%lad_warea)
deallocate(lpj_state%lpj_land_mask)
deallocate( &
lpj_state%coast_src_i, &
lpj_state%coast_src_j, &
lpj_state%coast_dst_i, &
lpj_state%coast_dst_j)
module_is_initialized = .false.
end subroutine lpj_driver_end
!
subroutine check_land_masks(regridded_mask, native_mask, resname, srcname, dstname)
real, dimension(:,:), intent(in) :: regridded_mask, native_mask
character(*), intent(in) :: resname, srcname, dstname
real, dimension(size(regridded_mask, 1),size(regridded_mask,2)) :: mask_check
integer :: n
! Check matching of regridded interpolation output mask with ``real'' native destination land mask
mask_check = regridded_mask
where(native_mask > 0) mask_check = 0.0 ! remaining are cells where regridded src cells fall into dst native ocean cells
n = count(mask_check > 0)
if (n > 0) &
write(stdout(),*) 'WARNING: for ', n, ' cells on ', resname, &
' resolution the input from ', srcname, ' vanishes into ocean/antarctica'
! check if there are native dst cells which do not get input because
! the interpolated-from-src grid assumes they blong to the ocean
mask_check = native_mask
where (regridded_mask > 0) mask_check = 0.0 ! remaining are dst cells that do not get input
n = count(mask_check > 0)
if (n > 0) &
write(stdout(),'(a,i5,a,a,a,a,a)') 'WARNING: ', n, ' ', dstname, ' cells do not get input from ', srcname, ' because of horizontal interpolation errors'
!! TODO: make this an one-time static field
!if (id_lpj_mask_out > 0) &
! used = send_data(id_lpj_mask_out, regridded_mask, lpj_state%current_time)
end subroutine check_land_masks
! ============================================================================
!
!
! Accumulates the FMS data of one fast time step as input data for LPJ
!
!
! CO2 input data is delivered by the coupler on every fast time step.
! Thus we must accumulated it here, and save it for lpj_driver_update_slow()
!
! FMS land model does not pass the current time to the update functions,
! thus we must use start time and time step information to keep track
! of the current time here.
!
!
subroutine lpj_driver_update_fast( &
! input variables
time, & ! from soil
! on lad resolution
q_ca, & ! mandatory
qflux, dedq, & ! mandatory
evap1, &
dedt, &
gc, &
exlpj, &
dhdt, drdt, & ! mandatory
drag_q, & ! mandatory
p_surf, & ! mandatory
lprec, fprec, temp, swdown, lwnet, t_flux, wind, &
tr_flux_co2, dfdtr_co2 & ! optional
! no output variables here
)
type(time_type), intent(in) :: time
! instant values in FMS units. must be aggregated over each day for LPJ,
! and converted into the units expected by LPJ
real, dimension(:,:), intent(out) :: & ! on lad resolution
evap1, &
dedt, &
gc, &
exlpj
real, dimension(:,:), intent(in) :: & ! on lad resolution
q_ca, & ! [kg/kg] Canopy air specific humidity, needed to convert between CO2 dry ppm and kg CO2/kg wet air
qflux, & ! [kg/m2/s] water vapor flux near surface
dedq, & ! [kg/m2/s] Moisture flux humidity sensitivity
dhdt, & ! derivative of sens over T
drdt, & ! derivative of LW radiation over T
drag_q, & ! [m/s] Product of drag coefficient for water vapor by wind
p_surf ! [N/m2] == [Pa] == [kg/m/s2] Surface pressure
real, dimension(:,:), intent(in) :: &
lprec, & ! [kg/m2/s] liquid precipitation; LPJ needs daily sum [mm] == [kg/m2]
fprec, & ! [kg/m2/s] frozen precipitation; LPJ needs daily sum [mm] == [kg/m2]
temp, & ! [K] temperature; LPJ needs [C]
swdown, & ! [W/m2] shortwave downwelling; LPJ needs same
lwnet, & ! [W/m2] longwave net; LPJ needs same
t_flux, & ! [W/m2] Sensible heat flux; LPJ needs the same
wind ! [m/s] wind speed; LPJ needs same
real, dimension(:,:), intent(in), optional :: & ! on lad resolution
tr_flux_co2, & ! [kg CO2/m2/s]
dfdtr_co2 ! [kg wet air/m2/s] d(tracer_flux)/d(surface tracer)
integer :: Time_year, Time_month, Time_day, Time_hour, Time_minute, Time_second, fortytwo = 42
logical :: goodmorning, goodnight
integer :: used
! dont refer to optional argument. strange crashes occur if tr_flux_co2 isnot present
!real, dimension(size(tr_flux_co2,1),size(tr_flux_co2,2)) :: tfco2, dco2
real, dimension(size(wind,1),size(wind,2)) :: co2_ppmv ! CO2 concentration in ppmv
if (.not. do_lpj) return
call mpp_clock_begin(lpj_update_fast_clock_id)
if (enable_fpe) then
!! To enable strict floating point exceptions inside this module.
!! Taken from the example in the Intel Fortran Compiler Reference Manual.
!! Most probably will not work with any other compiler
call clearstatusfpqq() ! do not inherit signalled FP exceptions from FMS
original_fpe_flags = FOR_SET_FPE(strict_fpe_flags)
endif
!lpj_state%current_time = lpj_state%current_time + lpj_state%dt_fast
lpj_state%current_time = time
call get_date(lpj_state%current_time, Time_year, Time_month, Time_day, Time_hour, Time_minute, Time_second)
goodmorning = lpj_state%yesterday /= Time_day;
goodnight = Time_hour*3600+Time_minute*60+Time_second + lpj_state%dt_fast_seconds >= SECONDS_PER_DAY
! write(*,*)
! write(*,*)
! write(*,*)
!write(*,'(a,i4,a,i2.2,a,i2.2,a,i2.2,a,i2.2,a,i2.2,a,l2,a,l2)') &
! 'lpj_driver_update_fast ',Time_year,'.',Time_month,'.',Time_day, &
! ' ',Time_hour,':',Time_minute,':',Time_Second, &
! ' goodmorning',goodmorning,' goodnight', goodnight
if (goodmorning) then
! write(*,'(a,i4,a,i2.2,a,i2.2,a,i2.2,a,i2.2,a,i2.2)') &
! 'Model time:',Time_year,'.',Time_month,'.',Time_day, &
! ' ',Time_hour,':',Time_minute,':',Time_Second
!write(*,*) ' - goodmorning: zero aggregations'
!lpj_state%nfaststepstoday = 0
!lpj_state%lprec = 0.0
!lpj_state%fprec = 0.0
lpj_state%prec = 0.0
!lpj_state%lpj_temp_mean = 288.0-KELVIN ! default initial condition from soil.F90
!lpj_state%lpj_temp_max = 288.0-KELVIN ! default initial condition from soil.F90
!lpj_state%lpj_temp_min = 288.0-KELVIN ! default initial condition from soil.F90
lpj_state%temp_mean = 0.0
lpj_state%temp_max = -1000.0
lpj_state%temp_min = 1000.0
lpj_state%swdown = 0.0
lpj_state%lwnet = 0.0
lpj_state%t_flux = 0.0
lpj_state%qflux = 0.0
lpj_state%dedq = 0.0
lpj_state%dhdt = 0.0
lpj_state%drdt = 0.0
lpj_state%drag_q = 0.0
lpj_state%p_surf = 0.0
lpj_state%wind = 0.0
! just continue the integration of flux values over the coming day
!lpj_state%co2 = 0.0
end if ! goodmorning
! aggregate in FMS-provided units
!lpj_state%nfaststepstoday = lpj_state%nfaststepstoday + 1
if (present(tr_flux_co2)) then
!write(*,*) ' - collecting aggregation for CO2'
! This yields kg CO2/kg wet air. LPJ needs [ppm]
!tfco2 = tr_flux_co2(:,:,1)
!dco2 = dfdtr_co2(:,:,1)
!lpj_state%co2 = lpj_state%co2 + tfco2 / (canopy_air_mass_for_tracers/lpj_state%dt_fast_seconds + dco2)
lpj_state%co2 = lpj_state%co2 &
+ ( lpj_state%co2_flux_out - tr_flux_co2) &
/ (canopy_air_mass_for_tracers/lpj_state%dt_fast_seconds + dfdtr_co2)
if (id_fco2 > 0) &
used = send_data(id_fco2, tr_flux_co2, lpj_state%current_time, &
rmask=lpj_state%lad_land_mask)
if (id_Dfco2Dq > 0) &
used = send_data(id_Dfco2Dq, dfdtr_co2, lpj_state%current_time, &
rmask=lpj_state%lad_land_mask)
end if ! present(tr_flux_co2)
! init output variables with default / missing values
! TODO: missing values and masks
lpj_state%co2_flux_out = 0
lpj_state%q_ca_out = 0
lpj_state%albedo_out = 0
if (goodnight) lpj_state%runoff_out = 0 ! recomputed only on goodnight steps
lpj_state%roughness_length_out = 0
lpj_state%surface_temp_out = 288.0-KELVIN ! default initial condition from soil.F90 ???
!write(*,*) 'before lpj update'
! tracer fluxes (humidity and CO2) are aggregated on fast time step
call lpj_update(Time_year, Time_month, Time_day, Time_hour, Time_minute, Time_second, &
! pass climate status information
!lpj_state%lpj_fprec, lpj_state%fpj_fprec, &
(fprec + lprec)* lpj_state%dt_fast_seconds, &
temp - KELVIN, &
swdown, lwnet, &
t_flux, &
qflux, dedq, &
dhdt, drdt, &
drag_q, p_surf, &
wind, &
co2_ppmv, &
! get back output from LPJ
lpj_state%co2_flux_out, & ! returned as carbon flux in [g C/ (m^2/day)], must be converted to [kg CO2/(m^2 sec)]
lpj_state%q_ca_out, &
lpj_state%runoff_out, & ! recomputed only on goodnight steps, returned in [l/s]
lpj_state%roughness_length_out, &
lpj_state%surface_temp_out, & ! returned in [degC]
lpj_state%albedo_out, &
lpj_state%evap1_out, &
lpj_state%dedt_out, &
lpj_state%gc_out, &
lpj_state%exlpj_out, &
!lpj_state%dt_fast_seconds, &
! last arg, minimal safety against interface changes
fortytwo &
)
!write(*,*) 'update lpj done'
!if (goodnight) &
!write(*,*) 'SURFTEMP_LPJOUT=', maxval(lpj_state%surface_temp_out), 'min=', minval(lpj_state%surface_temp_out), 'mean=', sum(lpj_state%surface_temp_out)/size(lpj_state%surface_temp_out)
evap1 = lpj_state%evap1_out
dedt = lpj_state%dedt_out
gc = lpj_state%gc_out
exlpj = lpj_state%exlpj_out
! diagnostic output
!
! input to LPJ
!
if (id_prec > 0) &
used = send_data(id_prec, lpj_state%prec, lpj_state%current_time, &
rmask=lpj_state%lad_land_mask)
if (id_temp_mean > 0) &
used = send_data(id_temp_mean, lpj_state%temp_mean, lpj_state%current_time, &
rmask=lpj_state%lad_land_mask)
if (id_temp_min > 0) &
used = send_data(id_temp_min, lpj_state%temp_min, lpj_state%current_time, &
rmask=lpj_state%lad_land_mask)
if (id_temp_max > 0) &
used = send_data(id_temp_max, lpj_state%temp_max, lpj_state%current_time, &
rmask=lpj_state%lad_land_mask)
if (id_swdown > 0) &
used = send_data(id_swdown, lpj_state%swdown, lpj_state%current_time, &
rmask=lpj_state%lad_land_mask)
if (id_lwnet > 0) &
used = send_data(id_lwnet, lpj_state%lwnet, lpj_state%current_time, &
rmask=lpj_state%lad_land_mask)
if (id_t_flux > 0) &
used = send_data(id_t_flux, lpj_state%t_flux, lpj_state%current_time, &
rmask=lpj_state%lad_land_mask)
if (id_qflux > 0) &
used = send_data(id_qflux, lpj_state%qflux, lpj_state%current_time, &
rmask=lpj_state%lad_land_mask)
if (id_dedq > 0) &
used = send_data(id_dedq, lpj_state%dedq, lpj_state%current_time, &
rmask=lpj_state%lad_land_mask)
if (id_dhdt > 0) &
used = send_data(id_dhdt, lpj_state%dhdt, lpj_state%current_time, &
rmask=lpj_state%lad_land_mask)
if (id_drdt > 0) &
used = send_data(id_drdt, lpj_state%drdt, lpj_state%current_time, &
rmask=lpj_state%lad_land_mask)
if (id_wind > 0) &
used = send_data(id_wind, lpj_state%wind, lpj_state%current_time, &
rmask=lpj_state%lad_land_mask)
! CO2 is updated on every time step.
!if (id_co2 > 0) & ! write out CO2 input from coupler in [kg CO2/kg wet air]
! used = send_data(id_co2, lpj_state%co2, lpj_state%current_time, &
! rmask=lpj_state%lad_land_mask)
if (id_co2_ppmv > 0) & ! write out CO2 input in ppmv dry air
used = send_data(id_co2_ppmv, co2_ppmv, lpj_state%current_time, &
rmask=lpj_state%lad_land_mask)
!
! LPJ output
!
if (id_carbon_flux_out > 0) & ! [g Carbon/(m^2*day)]
used = send_data(id_carbon_flux_out, lpj_state%co2_flux_out, lpj_state%current_time, &
rmask=lpj_state%lpj_land_mask)
!write(*,*) 'sending to id_q_ca_out ', id_q_ca_out, ' at ', Time_year, Time_month, Time_day, Time_hour, Time_minute, Time_second
if (id_q_ca_out > 0) &
used = send_data(id_q_ca_out, lpj_state%q_ca_out, lpj_state%current_time &
!, rmask=lpj_state%lpj_land_mask &
)
!if (maxval(lpj_state%runoff_out) == 0.0) __NOTE__('lpj_update (fast): runoff_out is 0')
if (id_runoff_out > 0) &
used = send_data(id_runoff_out, lpj_state%runoff_out, lpj_state%current_time)
if (id_roughness_length_out > 0) &
used = send_data(id_roughness_length_out, lpj_state%roughness_length_out, lpj_state%current_time, &
rmask=lpj_state%lpj_land_mask)
if (id_surface_temp_out > 0) &
used = send_data(id_surface_temp_out, lpj_state%surface_temp_out, lpj_state%current_time, &
rmask=lpj_state%lpj_land_mask)
if (id_albedo_out > 0) &
used = send_data(id_albedo_out, lpj_state%albedo_out, lpj_state%current_time, &
rmask=lpj_state%lpj_land_mask)
if (id_evap1_out > 0) &
used = send_data(id_evap1_out, lpj_state%evap1_out, lpj_state%current_time, &
rmask=lpj_state%lpj_land_mask)
if (id_dedt_out > 0) &
used = send_data(id_dedt_out, lpj_state%dedt_out, lpj_state%current_time, &
rmask=lpj_state%lpj_land_mask)
if (id_gc_out > 0) &
used = send_data(id_gc_out, lpj_state%gc_out, lpj_state%current_time, &
rmask=lpj_state%lpj_land_mask)
if (id_exlpj_out > 0) &
used = send_data(id_exlpj_out, lpj_state%exlpj_out, lpj_state%current_time, &
rmask=lpj_state%lpj_land_mask)
! Convert to land_lad / FMS conforming units
!lpj_state%co2_flux_out = lpj_state%co2_flux_out / SECONDS_PER_DAY * WTMCO2/WTMC / 1000 ! from [g Carbon/(m^2*day)] to [kg CO2/(m^2*sec)]
!if (id_co2_flux_out > 0) &
! used = send_data(id_co2_flux_out, lpj_state%co2_flux_out, lpj_state%current_time, &
! rmask=lpj_state%lpj_land_mask)
! runoff is recomputed only at goodnight timesteps.
if (goodnight) then
call pour_into_ocean
end if ! goodnight
! roughness_length_out
lpj_state%surface_temp_out = lpj_state%surface_temp_out + KELVIN ! from [C] to [K]
! albedo_out
! if (.false.) then !from old version, can be deleted
!
! ! output CO2 concentration on every fast timestep,
! ! because it is (potentially) updated on every fast timestep.
! if (id_co2 > 0) & ! write out CO2 input to LPJ from coupler in [kg CO2/kg wet air]
! used = send_data(id_co2, lpj_state%co2, lpj_state%current_time, &
! rmask=lpj_state%lad_land_mask)
! write(*,*) 'dtslowsteps', lpj_state%dt_slow_seconds
! write(*,*) ' - goodnight: finish aggregations'
! write(*,*) 'nslowtimestep=',lpj_state%nslowstepstoday
!
! write(*,*) 'dtslowsteps', lpj_state%dt_slow_seconds
! write(*,*) ' - goodnight: finish aggregations'
! write(*,*) 'nslowtimestep=',lpj_state%nslowstepstoday
!
!! LPJ wants ppm based on volume, not on mass.
! ! modelled after cold start in land_lad2/canopy_air/canopy_air.F90
! co2_ppmv = lpj_state%co2 * WTMAIR / (WTMCO2 * (1-q_ca)) * 1e6 ! [kg CO2/kg wet air] -> ppmv dry air
!
! if (do_check_land_masks) &
! call check_land_masks(lpj_state%lad_land_mask, lpj_state%lpj_land_mask, 'LPJ', 'land_lad', 'LPJ')
!
! call lpj_update(...)
!
! if (do_check_land_masks) &
! call check_land_masks(lpj_state%lpj_land_mask, lpj_state%lad_land_mask, 'land_lad', 'LPJ', 'land_lad')
!
!
!
!!!! XXX for debugging
! !call diag_manager_end(time)
! !stop
!
!
! else ! not goodnight
! ! Invoke C++ side to advance time and things, but avoid the overhead for aggregation and interpolation.
! ! The actual climate info is not accessed for those maintenance works,
! ! thus just pass NULL.
! ! Sigh deeply. Passing NULL() is rejected by ifort version 15.
! ! Workaround: pass the actual arrays; in this invocation they will not be used inside lpj_update() anyway.
!
! if (.false.) &
! call lpj_update(...)
!
! endif ! goodnight
lpj_state%yesterday = Time_day
if (enable_fpe) then
!! To enable strict floating point exceptions inside this module.
!! Taken from the example in the Intel Fortran Compiler Reference Manual.
!! Most probably will not work with any other compiler
original_fpe_flags = FOR_SET_FPE(original_fpe_flags)
endif
where (lpj_state%q_ca_out>0.2)
lpj_state%q_ca_out=0.2
endwhere
where (lpj_state%q_ca_out<0)
lpj_state%q_ca_out=0
endwhere
call mpp_clock_end(lpj_update_fast_clock_id)
end subroutine lpj_driver_update_fast
!
! ============================================================================
!
!
! Advances the LPJ model by one (slow!) time step
!
!
!
!
! subroutine lpj_driver_update_slow( &
! ! input variables
! !time, &
! ! climate status information on lad resolution
! !!lprec, fprec, temp, swdown, lwnet, t_flux, wind &
! ! qflux and dedq are accumulated on fast time step
! ! , co2 & ! CO2 is accumulated on fast time step
! ! output variables on lad resolution
! ! no output variables here, output is written back in the _update_bnd_ subroutines.
! !co2_flux_out, q_ca_out, runoff_out, roughness_length_out, surface_temp_out &
! ! TODO: folia protective cover [%] &
! ! TODO: LAI &
! )
!
! !type(time_type), intent(in) :: time
! ! instant values in FMS units. must be aggregated over each day for LPJ,
! ! and converted into the units expected by LPJ
! !!real, dimension(:,:), intent(in) :: &
! !! lprec, & ! [kg/m2/s] liquid precipitation; LPJ needs daily sum [mm] == [kg/m2]
! !! fprec, & ! [kg/m2/s] frozen precipitation; LPJ needs daily sum [mm] == [kg/m2]
! !! temp, & ! [K] temperature; LPJ needs [C]
! !! swdown, & ! [W/m2] shortwave downwelling; LPJ needs same
! !! lwnet, & ! [W/m2] longwave net; LPJ needs same
! !! t_flux, & ! [W/m2] Sensible heat flux; LPJ needs the same
! !! wind ! [m/s] wind speed; LPJ needs same
!
! if (.not. do_lpj) return
! call mpp_clock_begin(lpj_update_slow_clock_id)
!
! if (enable_fpe) then
! !! To enable strict floating point exceptions inside this module.
! !! Taken from the example in the Intel Fortran Compiler Reference Manual.
! !! Most probably will not work with any other compiler
! call clearstatusfpqq() ! do not inherit signalled FP exceptions from FMS
! original_fpe_flags = FOR_SET_FPE(strict_fpe_flags)
! endif
! !write(*,*) 'swdown=', sum(swdown)/size(swdown)
! write(*,*) 'lpj_driver_update_slow - collecting aggregation'
! !write(*,'(a,f7.2,a,f7.2,a,f7.2)') ' temp min ',minval(temp),' max ',maxval(temp),' mean ',sum(temp)/size(temp)
!
! ! Accumulate the climate status information for LPJ over all slow time steps in this day.
! ! Should be done here and not on the C-side to minimize interpolation overhead.
! ! aggregate in FMS-provided units
! ! re-set to 0 is done at goodmorning time steps in lpj_driver_update_fast()
! lpj_state%nslowstepstoday = lpj_state%nslowstepstoday + 1
! !lpj_state%lprec = lpj_state%lprec + lprec * lpj_state%dt_slow_seconds
! !lpj_state%fprec = lpj_state%fprec + fprec * lpj_state%dt_slow_seconds
!
!
!! lpj_state%prec = lpj_state%prec + (fprec + lprec) ! * lpj_state%dt_slow_seconds
!! lpj_state%temp_mean = lpj_state%temp_mean + temp
!! lpj_state%temp_max = max(lpj_state%temp_max, temp)
!! lpj_state%temp_min = min(lpj_state%temp_min, temp)
!! lpj_state%swdown = lpj_state%swdown + swdown
!! lpj_state%lwnet = lpj_state%lwnet + lwnet
!! lpj_state%t_flux = lpj_state%t_flux + t_flux
!! lpj_state%wind = lpj_state%wind + wind
!
!
!!lpj_state%co2 = lpj_state%co2 + co2 ! on fast time step
!
! if (enable_fpe) then
! !! To enable strict floating point exceptions inside this module.
! !! Taken from the example in the Intel Fortran Compiler Reference Manual.
! !! Most probably will not work with any other compiler
! original_fpe_flags = FOR_SET_FPE(original_fpe_flags)
! endif
!
! call mpp_clock_end(lpj_update_slow_clock_id)
!
! end subroutine lpj_driver_update_slow
!
! ============================================================================
!
!
! Override land_lad land boundary data data with LPJ data
!
!
! Updates land boundary data (the ones that atmosphere sees) on the
! fast time scale.
! The results of the most recent lpj_update() invocation are time-interpolated
! onto the land_lad fast time steps and written into the land_data_type structure.
!
!
!
subroutine update_lpj_bnd_fast(bnd)
type(land_data_type), intent(inout) :: bnd ! state to update
integer :: used
!integer :: Time_year, Time_month, Time_day, Time_hour, Time_minute, Time_second
!call get_date(lpj_state%current_time, Time_year, Time_month, Time_day, Time_hour, Time_minute, Time_second)
!write(*,'(a,i4,a,i2.2,a,i2.2,a,i2.2,a,i2.2,a,i2.2)') &
! 'update_lpj_bnd_fast at', Time_year,'.',Time_month,'.',Time_day, &
! ' ',Time_hour,':',Time_minute,':',Time_second
if (.not. do_lpj) return
if (enable_fpe) then
!! To enable strict floating point exceptions inside this module.
!! Taken from the example in the Intel Fortran Compiler Reference Manual.
!! Most probably will not work with any other compiler
call clearstatusfpqq() ! do not inherit signalled FP exceptions from FMS
original_fpe_flags = FOR_SET_FPE(strict_fpe_flags)
endif
! now time-interpolate the most recent daily output onto the individual time steps
! of the following (!!) day
if (use_co2_from_lpj .and. ico2 /= NO_TRACER) then
! put CO2 into tracer
!!!!! TODO: according to the comments for land_lad2/canopy/cana_tile.F90::cana_prog_type%co2,
!!!!! we must return co2 concentration in [kg CO2/kg of wet air] .
!!!!! That conforms to the document unit of evaporation as humidity tracer in [kg/kg] .
bnd%tr(:,:,1,ico2) = lpj_state%co2(:,:)
end if
if (id_co2_out_lad > 0 .and. ico2 /= NO_TRACER) then
used = send_data(id_co2_out_lad, lpj_state%co2(:,:), lpj_state%current_time, &
rmask=lpj_state%lad_land_mask)
end if
!!
!! override values which are set in update_vegetation_bnd_fast()
!!
!write(*,*) 'computed from land_lad'
!write(*,'(a,i4,i4,a,i3,i3,i4,i4)') ' size tr(:,:,1,isphum) ', size(bnd%tr(:,:,1,isphum),1), size(bnd%tr(:,:,1,isphum),2), &
! ' bounds tr(:,:,1,isphum) ', lbound(bnd%tr(:,:,1,isphum),1), lbound(bnd%tr(:,:,1,isphum),2), &
! ubound(bnd%tr(:,:,1,isphum),1), ubound(bnd%tr(:,:,1,isphum),2)
!write(*,'(a,g11.3,a,g11.3,a,g11.3,a,Z20)') ' tr(:,:,1,isphum) min ',minval(bnd%tr(:,:,1,isphum)),' max ',maxval(bnd%tr(:,:,1,isphum)),' mean ',sum(bnd%tr(:,:,1,isphum))/size(bnd%tr(:,:,1,isphum)), ' checksum', mpp_chksum(bnd%tr(:,:,1,isphum))
if (use_q_ca_from_lpj) then
! put evaporation into humidity tracer
! 3rd dimension n_tiles must be fixed to 1
!write(*,'(a,i4,i4,a,i3,i3,i4,i4)') 'overriding from lpj with size q_ca_out ', size(lpj_state%q_ca_out(:,:),1), size(lpj_state%q_ca_out(:,:),2),' bounds q_ca_out ', lbound(lpj_state%q_ca_out(:,:),1), lbound(lpj_state%q_ca_out(:,:),2), ubound(lpj_state%q_ca_out(:,:),1), ubound(lpj_state%q_ca_out(:,:),2)
!write(*,'(a,g11.3,a,g11.3,a,g11.3,a,Z20)') ' q_ca_out min ',minval(lpj_state%q_ca_out(:,:)),' max ',maxval(lpj_state%q_ca_out(:,:)),' mean ',sum(lpj_state%q_ca_out(:,:))/size(lpj_state%q_ca_out(:,:)), ' checksum', mpp_chksum(lpj_state%q_ca_out(:,:))
! where(lpj_state%q_ca_out(:,:)<0)
! lpj_state%q_ca_out(:,:)=0;
! endwhere
! where(lpj_state%q_ca_out(:,:)>0.1)
! lpj_state%q_ca_out(:,:)=0.1
! endwhere
bnd%tr(:,:,1,isphum) = lpj_state%q_ca_out(:,:)
!write(*,'(a,g11.3,a,g11.3,a,g11.3,a,Z20)') ' tr(:,:,1,isphum) min ',minval(bnd%tr(:,:,1,isphum)),' max ',maxval(bnd%tr(:,:,1,isphum)),' mean ',sum(bnd%tr(:,:,1,isphum))/size(bnd%tr(:,:,1,isphum)), ' checksum', mpp_chksum(bnd%tr(:,:,1,isphum))
else
!write(*,*) 'NOT overriding from lpj with'
!write(*,'(a,g11.3,a,g11.3,a,g11.3,a,Z20)') ' q_ca_out min ',minval(lpj_state%q_ca_out(:,:)),' max ',maxval(lpj_state%q_ca_out(:,:)),' mean ',sum(lpj_state%q_ca_out(:,:))/size(lpj_state%q_ca_out(:,:)), ' checksum', mpp_chksum(lpj_state%q_ca_out(:,:))
end if
if (id_q_ca_out_lad > 0) then
used = send_data(id_q_ca_out_lad, lpj_state%q_ca_out, lpj_state%current_time, &
rmask=lpj_state%lad_land_mask)
end if
!!
!! end of override values which are set in update_vegetation_bnd_fast()
!!
!!
!! override values which are set in update_soil_bnd_fast()
!!
!write(*,*) 'computed from land_lad'
!write(*,'(a,i3,i3,a,i3,i3,i3,i3)') ' size t_surf ', size(bnd%t_surf(:,:,1),1), size(bnd%t_surf(:,:,1),2), &
! ' bounds t_surf ', lbound(bnd%t_surf(:,:,1),1), lbound(bnd%t_surf(:,:,1),2), &
! ubound(bnd%t_surf(:,:,1),1), ubound(bnd%t_surf(:,:,1),2)
!write(*,'(a,i3,i3,a,i3,i3,i3,i3)')' size t_ca ', size(bnd%t_ca(:,:,1),1), size(bnd%t_ca(:,:,1),2), &
! ' bounds t_ca ', lbound(bnd%t_ca(:,:,1),1), lbound(bnd%t_ca(:,:,1),2), &
! ubound(bnd%t_ca(:,:,1),1), ubound(bnd%t_ca(:,:,1),2)
!write(*,'(a,f7.2,a,f7.2,a,f7.2,a,Z20)') ' t_surf min ',minval(bnd%t_surf(:,:,1)),' max ',maxval(bnd%t_surf(:,:,1)),' mean ',sum(bnd%t_surf(:,:,1))/size(bnd%t_surf(:,:,1)), ' checksum', mpp_chksum(bnd%t_surf(:,:,1))
!write(*,'(a,f7.2,a,f7.2,a,f7.2,a,Z20)') ' t_ca min ',minval(bnd%t_ca(:,:,1)),' max ',maxval(bnd%t_ca(:,:,1)),' mean ',sum(bnd%t_ca(:,:,1))/size(bnd%t_ca(:,:,1)), ' checksum', mpp_chksum(bnd%t_ca(:,:,1))
if (use_surface_temp_from_lpj) then
!lpj_state%surface_temp_out(:,:) = bnd%t_surf(:,:,1) ! for testing
where (lpj_state%surface_temp_out > 400)
lpj_state%surface_temp_out = 400
endwhere
where (lpj_state%surface_temp_out < 190)
lpj_state%surface_temp_out = 190
endwhere
bnd%t_surf(:,:,1) = lpj_state%surface_temp_out(:,:)
bnd%t_ca(:,:,1) = lpj_state%surface_temp_out(:,:)
! write(*,'(a,i3,i3,a,i3,i3,i3,i3)') 'overriding from lpj with size s_t_o ', size(lpj_state%surface_temp_out(:,:)), ' bounds s_t_o ', lbound(lpj_state%surface_temp_out(:,:)), ubound(lpj_state%surface_temp_out(:,:))
! write(*,'(a,f7.2,a,f7.2,a,f7.2,a,Z20)') ' temp_out min ',minval(lpj_state%surface_temp_out(:,:)),' max ',maxval(lpj_state%surface_temp_out(:,:)),' mean ',sum(lpj_state%surface_temp_out(:,:))/size(lpj_state%surface_temp_out(:,:)), ' checksum', mpp_chksum(lpj_state%surface_temp_out(:,:))
! write(*,'(a,f7.2,a,f7.2,a,f7.2,a,Z20)') ' t_surf min ',minval(bnd%t_surf(:,:,1)),' max ',maxval(bnd%t_surf(:,:,1)),' mean ',sum(bnd%t_surf(:,:,1))/size(bnd%t_surf(:,:,1)), ' checksum', mpp_chksum(bnd%t_surf(:,:,1))
! write(*,'(a,f7.2,a,f7.2,a,f7.2,a,Z20)') ' t_ca min ',minval(bnd%t_ca(:,:,1)),' max ',maxval(bnd%t_ca(:,:,1)),' mean ',sum(bnd%t_ca(:,:,1))/size(bnd%t_ca(:,:,1)), ' checksum', mpp_chksum(bnd%t_ca(:,:,1))
else
! write(*,*) 'NOT overriding from lpj with'
! write(*,'(a,f7.2,a,f7.2,a,f7.2,a,Z20)') ' temp_out min ',minval(lpj_state%surface_temp_out(:,:)),' max ',maxval(lpj_state%surface_temp_out(:,:)),' mean ',sum(lpj_state%surface_temp_out(:,:))/size(lpj_state%surface_temp_out(:,:)), ' checksum', mpp_chksum(lpj_state%surface_temp_out(:,:))
end if
if (id_surface_temp_out_lad > 0) then
used = send_data(id_surface_temp_out_lad, lpj_state%surface_temp_out, lpj_state%current_time &
! , rmask=lpj_state%lad_land_mask &
)
end if
!bnd % albedo = soil % albedo
!bnd % albedo_vis_dir = soil%albedo_vis_dir
!bnd % albedo_nir_dir = soil%albedo_nir_dir
!bnd % albedo_vis_dif = soil%albedo_vis_dif
!bnd % albedo_nir_dif = soil%albedo_nir_dif
!TODO: albedo values from LPJ must be merged/blended with glacier / ice / snow albedos from land_lad
if (use_albedo_from_lpj) then
bnd%albedo(:,:,1) = lpj_state%albedo_out(:,:)
!for now set all to the same value, TODO: caculate seperate values within lpjml
!it looks like land lad also sets one value for all components
bnd % albedo_vis_dir(:,:,1) = lpj_state%albedo_out(:,:)
bnd % albedo_nir_dir(:,:,1) = lpj_state%albedo_out(:,:)
bnd % albedo_vis_dif(:,:,1) = lpj_state%albedo_out(:,:)
bnd % albedo_nir_dif(:,:,1) = lpj_state%albedo_out(:,:)
end if
if (id_albedo_out_lad > 0) then
used = send_data(id_albedo_out_lad, lpj_state%albedo_out, lpj_state%current_time, &
rmask=lpj_state%lad_land_mask)
!!
! ============================================================================
!
!
! Pour the water from LPJ runoff into the ocean cells.
!
!
! Parts of this code was taken from / inspired by land_lad/soil/river.F90
!
! Actually pour the water from LPJ cells into ocean cells,
! using the prepared lpj_runoff_coast_map_i/j .
! lpj_state%runoff_out was updated at goodnight timestep in lpj_update_fast()
! Here, it is assembled into a global map, and then forwarded into discharge,
! which then is copied into the corresponding FMS data structure in
! subroutine update_lpj_bnd_slow(). There also the division by
! per-cell ocean area is done, as in land_lad update_rivers_bnd_slow()
!
! This subroutine is invoked during init to pour out the runoff that
! is read from the LPJ restart file,
! and at goodnight timesteps.
!
subroutine pour_into_ocean
integer :: i, cellspoured
real, dimension(1:nlons_global,1:nlats_global) :: global_runoff_out
!lpj_state%runoff_out(:,:) = lpj_state%runoff_out(:,:) / SECONDS_PER_DAY ! from daily sum [kg/m^2/day] to rate [kg/m^2/sec]
global_runoff_out(:,:) = 0.0
!write(*,'(a,4(i4),2(2(i4)))') 'is,ie,js,je,bounds(runoff_out)', &
! is, ie, js, je, &
! lbound(lpj_state%runoff_out), ubound(lpj_state%runoff_out)
!write(*,'(a,4(i4),2(2(i4)))') 'is,ie,js,je,bounds(discharge)', &
! is, ie, js, je, &
! lbound(lpj_state%discharge), ubound(lpj_state%discharge)
!write(*,'(a,2(2(i4)))') 'bounds(global_runoff_out)', lbound(global_runoff_out), ubound(global_runoff_out)
!if (maxval(lpj_state%runoff_out) == 0.0) __NOTE__('pour_into_ocean: runoff_out is 0')
global_runoff_out(is:ie,js:je) = lpj_state%runoff_out(is:ie,js:je)
!call mpp_sum(global_runoff_out, size(global_runoff_out))
call mpp_global_field(lpj_state%domain, lpj_state%runoff_out, global_runoff_out)
! now pour the water into the oceans
lpj_state%discharge(:,:) = 0
cellspoured = 0
do i=1,lpj_state%ncoastcells
if (lpj_state%coast_dst_i(i) >= is &
.and. lpj_state%coast_dst_i(i) <= ie &
.and. lpj_state%coast_dst_j(i) >= js &
.and. lpj_state%coast_dst_j(i) <= je) then
cellspoured = cellspoured + 1
lpj_state%discharge(lpj_state%coast_dst_i(i),lpj_state%coast_dst_j(i)) &
= lpj_state%discharge(lpj_state%coast_dst_i(i),lpj_state%coast_dst_j(i)) &
+ global_runoff_out(lpj_state%coast_src_i(i),lpj_state%coast_src_j(i))
!write(*,'(a,i5,2i4,a,2i4,2g12.5)') &
! 'pour ',i,lpj_state%coast_src_i(i),lpj_state%coast_src_j(i), &
! ' -> ',lpj_state%coast_dst_i(i),lpj_state%coast_dst_j(i), &
! global_runoff_out(lpj_state%coast_src_i(i),lpj_state%coast_src_j(i)), &
! lpj_state%discharge(lpj_state%coast_dst_i(i),lpj_state%coast_dst_j(i))
end if
end do
!write(*,*) 'poured from ',cellspoured, ' cells'
!if (maxval(lpj_state%discharge) == 0.0) __NOTE__('pour_into_ocean: discharge is 0')
end subroutine pour_into_ocean
!
! ============================================================================
!
!
! Override land_lad land boundary data data with LPJ data
!
!
! Updates land boundary data (the ones that atmosphere sees) on the
! slow time scale.
! The results of the most recent lpj_update() invocation are time-interpolated
! onto the land_lad slow time steps and written into the land_data_type structure.
!
!
!
subroutine update_lpj_bnd_slow(bnd)
type(land_data_type), intent(inout) :: bnd ! state to update
integer :: used
if (.not. do_lpj) return
! now time-interpolate the most recent daily output onto the individual time steps
! of the following (!!) day
if (enable_fpe) then
!! To enable strict floating point exceptions inside this module.
!! Taken from the example in the Intel Fortran Compiler Reference Manual.
!! Most probably will not work with any other compiler
call clearstatusfpqq() ! do not inherit signalled FP exceptions from FMS
original_fpe_flags = FOR_SET_FPE(strict_fpe_flags)
endif
!!
!! override values which are set in update_vegetation_bnd_slow()
!!
! no operation
!!
!! end of override values which are set in update_vegetation_bnd_slow()
!!
!!
!! override values which are set in update_soil_bnd_slow()
!!
if (use_roughness_length_from_lpj) then
!! TODO: is this the correct variable pairing ?
!write(*,*) 'roughness_lad=', maxval(bnd%rough_mom(:,:,1)), 'min=', minval(bnd%rough_mom(:,:,1)), 'mean=', sum(bnd%rough_mom(:,:,1))/size(bnd%rough_mom(:,:,1))
!write(*,*) 'roughness_lad_heat=', maxval(bnd%rough_heat(:,:,1)), 'min=', minval(bnd%rough_heat(:,:,1)), 'mean=', sum(bnd%rough_heat(:,:,1))/size(bnd%rough_heat(:,:,1))
!write(*,*) 'roughness_lad_scale=', maxval(bnd%rough_scale(:,:,1)), 'min=', minval(bnd%rough_scale(:,:,1)), 'mean=', sum(bnd%rough_scale(:,:,1))/size(bnd%rough_scale(:,:,1))
bnd%rough_mom(:,:,1) = lpj_state%roughness_length_out(:,:)
!set the same value for now. TODO. calculate all values within LPJ
bnd % rough_heat(:,:,1) = lpj_state%roughness_length_out(:,:)
bnd % rough_scale(:,:,1) = lpj_state%roughness_length_out(:,:)
!write(*,*) 'roughness=', maxval(bnd%rough_mom(:,:,1)), 'min=', minval(bnd%rough_mom(:,:,1)), 'mean=', sum(bnd%rough_mom(:,:,1))/size(bnd%rough_mom(:,:,1))
end if
if (id_roughness_length_out_lad > 0) then
used = send_data(id_roughness_length_out_lad, lpj_state%roughness_length_out, lpj_state%current_time, &
rmask=lpj_state%lad_land_mask)
end if
!! 0.0 )
! adjust partial-ocean-cell discharge values for area fraction
! and recalculate to kg/(m2 s)
bnd % discharge = bnd % discharge / lpj_state%lad_warea ! [l/s] -> [l/m^2/s] = [kg/(m2 s)]
!bnd % discharge_snow = bnd % discharge_snow / rivers%
elsewhere
bnd % discharge = 0
!bnd % discharge_snow = 0
endwhere
if (id_discharge_out_lad_aw > 0) &
used = send_data(id_discharge_out_lad_aw, bnd%discharge, lpj_state%current_time)
end if
#if 0
! for debugging, pass on in parallel to land_lad runoff
bnd%lpj_discharge(:,:) = lpj_state%discharge(is:ie,js:je)
where ( lpj_state%lad_warea > 0.0 )
! adjust partial-ocean-cell discharge values for area fraction
! and recalculate to kg/(m2 s)
bnd % lpj_discharge = bnd % lpj_discharge / lpj_state%lad_warea
elsewhere
bnd % lpj_discharge = 0
endwhere
#endif
! dont use land mask here, because runoff into ocean is just outside land boundaries
!if (maxval(lpj_state%runoff_out) == 0.0) __NOTE__('update_lpj_bnd_slow: runoff_out is 0')
if (id_runoff_out_lad > 0) &
used = send_data(id_runoff_out_lad, lpj_state%runoff_out, lpj_state%current_time)
if (id_discharge_out_lad > 0) &
used = send_data(id_discharge_out_lad, lpj_state%discharge, lpj_state%current_time)
!!
! ============================================================================
!
!
! Prepare coastal runoff routing map from LPJmL river routing map
!
!
! Parts of this code was taken from / inspired by land_lad/soil/river.F90
!
! The input arrays contain coordinates, and the LPJ cell index
! of the next cell in the river routing chain.
! -1 marks dead ends, i.e. inland sinks, or coastal river mouths.
! -9 means no data, i.e. ocean cells. Should not occur in this data.
! TODO: check that!!
! Follow the method used in the land_lad code to prepare
! the coastal runoff routing map:
! Iterate over the map.
! First look at dead ends of runoff routes, i.e. cells with
! value -1. For each such sink cell, search its vicinity for an
! ocean cell. If there are ocean cells nearby, store the i,j indices
! of the nearest ocean cell as runoff destination in a
! global coast cell map.
! Further, check for cells where the runoff destination already is
! an ocean cell. Copy the values of those cells also into the
! coast cell map.
! The computed routing information is stored in module-global arrays.
!
subroutine prepare_coastal_runoff_routing(global_number_lpj_cells, &
lpj_cell_lons, lpj_cell_lats, lpj_cell_next, &
lad_gfrac, &
id_runoff_map, id_runoff_coastcells, &
id_runoff_dest_lon, id_runoff_dest_lat, &
id_runoff_dcount, id_runoff_basin, &
id_runin_coastcells, id_inland_sinks)
integer, intent(in) :: global_number_lpj_cells
real, intent(inout), dimension(global_number_lpj_cells) :: lpj_cell_lons, lpj_cell_lats
integer, intent(inout), dimension(global_number_lpj_cells) :: lpj_cell_next
real, dimension(:,:), intent(in) :: lad_gfrac ! global fraction of the cell covered by land
integer, intent(in) :: id_runoff_map, id_runoff_coastcells
integer, intent(in) :: id_runoff_dest_lon, id_runoff_dest_lat
integer, intent(in) :: id_runoff_dcount, id_runoff_basin
integer, intent(in) :: id_runin_coastcells
integer, intent(in) :: id_inland_sinks
character(len=*), parameter :: module_name = 'prepare_coastal_runoff_routing'
integer :: i, j, cell, ii, this_i, this_j, next_i, next_j, used
real :: dlon, dlat, maxbasin
! count statistics of obtained mappings
integer :: n_10=0, n_9=0, n_1=0, n_other=0, n0=0, nvalid=0
integer :: ncoastcells = 0, nsinks_already = 0, nsinks_redirected = 0, nsinks_remain = 0
! Global runoff routing maps.
! For each coast cell, the i and j indices of the target ocean cell are stored,
! in Fortran numbering 1..nlons, 1..nlats
! Ocean cells are -9, inland sinks are -1.
integer, dimension(1:nlons_global,1:nlats_global) :: lpj_runoff_coast_map_i
integer, dimension(1:nlons_global,1:nlats_global) :: lpj_runoff_coast_map_j
! for static diagnostic output (dcount, basin)
real, dimension(nlons_global, nlats_global) :: map, diag
!real, dimension(is:ie,js:je) :: inland_sinks
! path to find redirection points in the vicinity
integer, parameter :: n_didj = 81
integer, parameter :: di(n_didj) = &
(/ -4, 4,-4, 4,-4, 4,-3, 3, 4,-4, 3,-3,-2, 2,-4, 4,-4, 4,-2, 2, &
3,-3, 3,-3,-1, 1, 4,-4,-4, 4,-1, 1, 0, 0,-4, 4, &
2,-3, 3,-3, 3,-3, 2,-2, &
1,-1, 3,-3,-3, 3, 1,-1, 3,-3, 0, 0, &
-2, 2,-2, 2,-1, 1, 2,-2,-2, 2,-1, 1, 0, 0,-2, 2, &
1,-1, 1,-1, 1,-1, 0, 0, 0/)
integer, parameter :: dj(n_didj) = &
(/ 4,-4,-4, 4, 3,-3,-4, 4, 3,-3,-4, 4, 4,-4,-2, 2,-2,-2,-4, 4, &
3,-3,-3, 3,-4, 4,-1, 1,-1, 1, 4,-4,-4, 4, 0, 0, &
3,-2,-2, 2, 2,-2,-3, 3, &
-3, 3, 1,-1, 1,-1, 3,-3, 0, 0, 3,-3, &
2,-2,-2, 2,-2, 2,-1, 1,-1, 1, 2,-2,-2, 2, 0, 0, &
1,-1,-1, 1, 0, 0, 1,-1, 0/)
character(len=nlons_global) :: buffer
logical :: changed
integer :: riverlength
integer :: riveradd ! how much dcount to add to new cells in the path
integer :: riverincr ! increment from cell cell, 1 or 0
__NOTE__("starting")
! first print map before the calls to mpp_sum.
! LPJ counts cells from 0, thus subtract 1 in the printout
!write(*,*) 'discharge map as obtained from LPJ, only local domain filled in:'
!do cell = 1, global_number_lpj_cells
! write(*,'(a,i5,f8.2,f8.2,a,i5)') 'discharge map A ',cell-1,lpj_cell_lons(cell),lpj_cell_lats(cell),' -> ',lpj_cell_next(cell)
!end do
dlon = 360.0/nlons_global
dlat = 180.0/nlats_global
write(*,'(a,f5.2,a,f5.2)') 'dlon ',dlon,' dlat ',dlat
! we (ab)use global sum. This is safe, because each task
! contributes non-zero values only for its local cells.
! This was verified using the above and below write statements.
! It works, despite strange warning from compiler and from
! valgrind about using uninitialized values. Sigh.
call mpp_sum(lpj_cell_lons, global_number_lpj_cells)
call mpp_sum(lpj_cell_lats, global_number_lpj_cells)
call mpp_sum(lpj_cell_next, global_number_lpj_cells)
! print map resulting from mpp_sum.
!write(*,*) 'discharge map after mpp_sum, now global domain should be filled in:'
!do cell = 1, global_number_lpj_cells
! write(*,'(a,i5,f8.2,f8.2,a,i5)') 'discharge map B ',cell-1,lpj_cell_lons(cell),lpj_cell_lats(cell),' -> ',lpj_cell_next(cell)
!end do
lpj_runoff_coast_map_i = -10
lpj_runoff_coast_map_j = -10
write(*,*) 'convert the LPJ next cell indices to FMS global grid indices i,j'
do cell = 1, global_number_lpj_cells
! LPJ counts cells from 0, thus subtract 1 in the printout
this_i = lon2i(lpj_cell_lons(cell), dlon)
this_j = lat2j(lpj_cell_lats(cell), dlat)
if (is_coastal(this_i,this_j,lad_gfrac)) ncoastcells = ncoastcells+1
!write(*,'(a,i5,f8.2,f8.2,i5,i5,a,i5$)') ' discharge map ',cell-1,lpj_cell_lons(cell),lpj_cell_lats(cell),this_i,this_j,' -> ',lpj_cell_next(cell)
__CHECK__(this_i >= 1 .and. this_i <= nlons_global, "this_i out of range")
__CHECK__(this_j >= 1 .and. this_j <= nlats_global, "this_j out of range")
if (lpj_cell_next(cell) >= 0) then
__CHECK__(lpj_cell_next(cell) < global_number_lpj_cells, "next out of range")
nvalid = nvalid+1
next_i = lon2i(lpj_cell_lons(lpj_cell_next(cell)+1), dlon)
next_j = lat2j(lpj_cell_lats(lpj_cell_next(cell)+1), dlat)
!write(*,'(a,f8.2,f8.2,i5,i5$)') ' ',lpj_cell_lons(lpj_cell_next(cell)+1),lpj_cell_lats(lpj_cell_next(cell)+1),next_i,next_j
__CHECK__(next_i >= 1 .and. next_i <= nlons_global, "next_i out of range")
__CHECK__(next_j >= 1 .and. next_j <= nlats_global, "next_j out of range")
!if (this_i == next_i .and. this_j == next_j) write(*,'(a$)') ' loop to self'
else if (lpj_cell_next(cell) == 0) then
n0 = n0+1
next_i = 0 ! ?
next_j = 0 ! ?
else if (lpj_cell_next(cell) == -1) then ! dead end
next_i = -1
next_j = -1
n_1 = n_1+1
else if (lpj_cell_next(cell) == -9) then ! ocean cell in land?
n_9 = n_9+1
next_i = lpj_cell_next(cell)
next_j = lpj_cell_next(cell)
else ! some other negative value
n_other = n_other+1
next_i = lpj_cell_next(cell)
next_j = lpj_cell_next(cell)
!write(*,'(a,i7,a,i7,a,f7.2,a,f7.2)') &
! 'Warning: '//module_name//' found river routing next == ', &
! lpj_cell_next(cell), &
! ' in cell ', cell, ' at ', &
! lpj_cell_lons(cell),',',lpj_cell_lats(cell)
end if
!write(*,*) ' ' ! just end the line
lpj_runoff_coast_map_i(this_i, this_j) = next_i
lpj_runoff_coast_map_j(this_i, this_j) = next_j
end do
write(*,*) 'got from LPJ -10/-9/-1/-X/0/valid/other', &
n_10, n_9, n_1, n_other, n0, nvalid, &
global_number_lpj_cells-(n_10+n_9+n_1+n_other+n0+nvalid)
write(*,*) 'river route map has ', ncoastcells, ' coast cells'
! Now the remaining -9 entries in lpj_runoff_coast_map_i (and _j)
! should match the ocean part of the lad_land_mask. Check that.
! (this could be done in above loop already,
! but the code would look uglier)
do j=1,nlats_global
do i=1,nlons_global
if (lad_gfrac(i,j) > 0) then
if (-9 == lpj_runoff_coast_map_i(i,j)) then
!write(*,'(a,i5,a,i5)') &
! 'Warning: '//module_name &
! //' found land cell in lad_gfrac' &
! //' but LPJ ocean cell in river routing next' &
! //' at index ', i,',',j
end if
if (-10 == lpj_runoff_coast_map_i(i,j)) then
!write(*,'(a,i5,a,i5)') &
! 'Warning: '//module_name &
! //' found land cell in lad_gfrac' &
! //' but lad ocean cell in river routing next' &
! //' at index ', i,',',j
end if
else
if (-9 /= lpj_runoff_coast_map_i(i,j)) then
!write(*,'(a,i5,a,i5)') &
! 'Warning: '//module_name &
! //' found LPJ ocean cell in lad_gfrac,' &
! //' but land cell in river routing next' &
! //' at index ', i,',',j
end if
if (-10 /= lpj_runoff_coast_map_i(i,j)) then
!write(*,'(a,i5,a,i5)') &
! 'Warning: '//module_name &
! //' found lad ocean cell in lad_gfrac,' &
! //' but land cell in river routing next' &
! //' at index ', i,',',j
end if
end if
end do
end do
!write(*,*) 'print a nice map'
!do j=nlats_global,1,-1 ! print N to S
! do i=1,nlons_global
! select case (lpj_runoff_coast_map_i(i,j))
! case (-10)
! buffer(i:i) = '.'
! case (-9)
! buffer(i:i) = '+'
! case (-1)
! buffer(i:i) = 'X'
! case (0)
! buffer(i:i) = '0'
! case default
! if (lpj_runoff_coast_map_i(i,j) == i .and. &
! lpj_runoff_coast_map_j(i,j) == j) then
! buffer(i:i) = '@' ! route into self
! else
! buffer(i:i) = '*'
! end if
! end select
! end do
! write(*,'(a)') buffer
!end do
__NOTE__("Now try to re-route coastal sink cells")
! In the first pass, just count the number of successes,
! so that we know how much space to allocate for the final map.
! TODO: Optimization: in each task, store and process only those parts
! of the coast runoff map, that pours into the local compute domain.
nsinks_redirected = 0
nsinks_already = 0
do j=1,nlats_global
do i=1,nlons_global
next_i = lpj_runoff_coast_map_i(i,j)
next_j = lpj_runoff_coast_map_j(i,j)
if (next_i > 0) then
if (is_coastal(next_i, next_j, lad_gfrac)) then
! cell i,j already routes into coast cell.
! Just count it.
nsinks_already = nsinks_already + 1
cycle
end if
end if
if (-1 == lpj_runoff_coast_map_i(i,j)) then
! try to find a nearby ocean cell to redirect runoff to.
! to avoid branching out of loop, currently search all
! possibilities, over-writing from most distant to least distant.
! backward search direction, from closest to farthest,
! and break out of loop.
do ii = size(di(:)), 1, -1
next_i = modulo(i+di(ii)-1, nlons_global) + 1
next_j = max(min(j+dj(ii), nlats_global), 1)
if (is_coastal(next_i, next_j, lad_gfrac)) then
! found a new destination that is OK
nsinks_redirected = nsinks_redirected + 1
exit ! from loop over ii
end if
end do
end if
end do
end do
nsinks_remain = n_1 - nsinks_redirected
write(*,*) nsinks_already, ' cells already routing into ocean cell'
write(*,*) 'redirected routing for ',nsinks_redirected,' coast cells.'
write(*,*) 'remaining inland sink cells: ',nsinks_remain
! Prepare a compressed coastal runoff map in the shape of 4 1-D arrays,
! on the assumption that this is cheaper to handle than full global maps.
!lpj_state%ncoastcells = nsinks_already+nsinks_redirected
lpj_state%ncoastcells = nsinks_redirected+nsinks_already
allocate( &
lpj_state%coast_src_i(nsinks_redirected+nsinks_already), &
lpj_state%coast_src_j(nsinks_redirected+nsinks_already), &
lpj_state%coast_dst_i(nsinks_redirected+nsinks_already), &
lpj_state%coast_dst_j(nsinks_redirected+nsinks_already))
__NOTE__("2nd pass")
! 2nd pass
map = 0
diag = -1
allocate(lpj_state%inland_sinks(is:ie,js:je))
lpj_state%inland_sinks = 0
cell = 0
do j=1,nlats_global
do i=1,nlons_global
next_i = lpj_runoff_coast_map_i(i,j)
next_j = lpj_runoff_coast_map_j(i,j)
if (next_i > 0) then
map(i, j) = (next_j-1)*nlons+next_i+1
if (is_coastal(next_i, next_j, lad_gfrac)) then
!write(*,'(a,i5,i5,a,i5,i5,a,f6.4)') &
! 'cell ',i,j,' already routes into coast cell ',next_i, next_j,' gfrac ',lad_gfrac(i,j)
cell = cell + 1
lpj_state%coast_src_i(cell) = i
lpj_state%coast_src_j(cell) = j
lpj_state%coast_dst_i(cell) = next_i
lpj_state%coast_dst_j(cell) = next_j
diag(i, j) = 0
end if
else
map(i, j) = next_i ! copy missing value
end if
if (-1 == next_i) then
! try to find a nearby ocean cell to redirect runoff to.
! to avoid branching out of loop, currently search all
! possibilities, over-writing from most distant to least distant.
! backward search direction, from closest to farthest,
! and break out of loop.
do ii = size(di(:)), 1, -1
next_i = modulo(i+di(ii)-1, nlons_global) + 1
next_j = max(min(j+dj(ii), nlats_global), 1)
if (is_coastal(next_i, next_j, lad_gfrac)) then
! found a new destination that is OK
lpj_runoff_coast_map_i(i,j) = next_i
lpj_runoff_coast_map_j(i,j) = next_j
!write(*,'(a,i5,a,i5,a,i5,a,i5,a,f6.4$)') 'redirect coastal runoff from ', &
! i,',',j,' to ', next_i,',',next_j,' gfrac ',lad_gfrac(i,j)
!if (i == next_i .and. j == next_j) write(*,'(a$)') ' redirected loop to self'
!write(*,*) ' '
cell = cell + 1
lpj_state%coast_src_i(cell) = i
lpj_state%coast_src_j(cell) = j
lpj_state%coast_dst_i(cell) = next_i
lpj_state%coast_dst_j(cell) = next_j
diag(i, j) = ii
map(i, j) = (next_j-1)*nlons+next_i+1
exit ! from loop over ii
end if
end do ! loop over ii
if (0 == ii) then
write(*,'(a,i5,a,i5)') 'no coastal redirection found from ',i,',',j
if (i>=is.and.i<=ie.and.j>=js.and.j<=je) lpj_state%inland_sinks(i,j) = 1 ! prepare diagnostic output of inland sink locations
end if
end if
end do
end do
if (id_runoff_map > 0) &
used = send_data(id_runoff_map, map(is:ie,js:je))
if (id_runoff_coastcells > 0) &
used = send_data(id_runoff_coastcells, diag(is:ie,js:je))
if (id_inland_sinks > 0) &
used = send_data(id_inland_sinks, lpj_state%inland_sinks(is:ie,js:je))
__CHECK__(cell == nsinks_redirected+nsinks_already, "2nd redirection pass got inconsistent result")
if (id_runoff_dest_lon > 0) then
diag = -1000.0
do j=js,je
do i=is,ie
if (lpj_runoff_coast_map_i(i,j) > 0) &
diag(i,j) = i2lon(lpj_runoff_coast_map_i(i,j), dlon)
end do
end do
used = send_data(id_runoff_dest_lon, diag(is:ie,js:je))
end if
if (id_runoff_dest_lat > 0) then
diag = -1000.0
do j=js,je
do i=is,ie
if (lpj_runoff_coast_map_j(i,j) > 0) &
diag(i,j) = j2lat(lpj_runoff_coast_map_j(i,j), dlat)
end do
end do
used = send_data(id_runoff_dest_lat, diag(is:ie,js:je))
end if
! static diagnostics:
if (id_runoff_dcount > 0 .or. id_runin_coastcells > 0) then
__NOTE__("dcount diagnostics")
!TODO: use array map to improve cycle detection by coloring it with
! the number of the currently followed river flow. If we encounter a map cell
! which already is colored with the current number, we have a cycle.
! That saves counting up to the limit.
diag = -10 ! ocean cells
map = 0
ii = 0 ! use here for river number
do j=1,nlats_global
do i=1,nlons_global
riverlength = 0
if (diag(i,j) /= -10) &
cycle ! We have been here before. No sense to add 0 all along downstream path
if (lpj_runoff_coast_map_i(i,j) <= 0) cycle
! follow down-river
!write(*,'(a,i5,i5,a,i5,i5)') 'follow from ',i,j,' to ',lpj_runoff_coast_map_i(i,j), lpj_runoff_coast_map_j(i,j)
this_i = i
this_j = j
ii = ii+1
riveradd = 1 ! how much to add to new cells in the path
riverincr = 1 ! increment from cell to cell, might be set to 0 below
riverlength = 1
!write(*,'(a,i5,i5$)') ' start ',i,j
do while (.true.) ! sigh. F90 does not have && shortcut evaluation like C
next_i = lpj_runoff_coast_map_i(this_i,this_j)
next_j = lpj_runoff_coast_map_j(this_i,this_j)
!write(*,'(a,i5,i5$)') ' -> ',next_i, next_j
if (next_i <= 0 .or. next_j <= 0) then ! sink
!write(*,*) ' sink nodata len ',riverlength
exit
end if
if (next_i == this_i .and. next_j == this_j) then ! sink to self
!write(*,*) ' sink self len ',riverlength
exit
end if
if (map(this_i,this_j) == ii & ! routing loop, or bug
.or. riverlength == 1000 &
.or. diag(next_i, next_j) >= global_number_lpj_cells &
) then
!write(*,*) ' loop len ',riverlength
!write(*,'(a,i5,a,i5,a,f7.2,f7.2,a,i7)') 'Warning: river routing loop at ',next_i,',',next_j,' = ',i2lon(next_i,dlon),j2lat(next_j,dlat),' len ',riverlength
exit
end if
map(this_i,this_j) = ii ! color it
!if (diag(next_i, next_j) >= global_number_lpj_cells) exit ! loop or bug
if (diag(next_i, next_j) == -10) then
diag(next_i, next_j) = riveradd
else
diag(next_i, next_j) = diag(next_i, next_j) + riveradd
! We have been here before. Stop increasing riveradd,
! in order to not count downstream cells multiple times.
riverincr = 0
end if
this_i = next_i ! and move on
this_j = next_j
riveradd = riveradd + riverincr
riverlength = riverlength + 1
end do
end do
end do
if (id_runoff_dcount > 0) &
used = send_data(id_runoff_dcount, diag(is:ie,js:je))
!write(*,*) 'print a nice dcount map'
!do j=nlats_global,1,-1 ! print N to S
! write(*,'(*(i6))') int(diag(1:nlons_global,j))
!end do
if (id_runin_coastcells > 0) then
__NOTE__("beach cell diagnostics")
! beach cells := ocean cells with number of river cells routing into them
map = -1.0
do cell=1,nsinks_redirected
this_i = lpj_state%coast_src_i(cell)
this_j = lpj_state%coast_src_j(cell)
next_i = lpj_state%coast_dst_i(cell)
next_j = lpj_state%coast_dst_j(cell)
if (map(next_i, next_j) < 0) map(next_i, next_j) = 0
map(next_i, next_j) = map(next_i, next_j) + diag(this_i, this_j)
end do
used = send_data(id_runin_coastcells, map(is:ie,js:je))
end if ! id_runin_coastcells > 0
end if ! id_runoff_dcount > 0 .or. id_runin_coastcells > 0
if (id_runoff_basin > 0) then
__NOTE__("basin diagnostics")
! basins of size 0, i.e. cells without outflow and without inflow get 0.
! All non-trivial basins get a unique number.
! Ocean cells get -10
diag = -10
maxbasin = 1
changed = .true. ! Termination flag. If anything was changed, we need another loop over the cells.
do while (changed)
changed = .false.
do j=1,nlats_global
do i=1,nlons_global
!write(*,*) 'looking at ',i,j, lpj_runoff_coast_map_i(i,j), lpj_runoff_coast_map_j(i,j)
if (lpj_runoff_coast_map_i(i,j) > 0) then
this_i = i
this_j = j
next_i = lpj_runoff_coast_map_i(this_i,this_j)
next_j = lpj_runoff_coast_map_j(this_i,this_j)
if (diag(this_i, this_j) == -10 & ! new cell
.and. (next_i <= 0 & ! no outflow
.or. (next_i == this_i .and. next_j == this_j))) then ! outflow to self
!write(*,*) ' new cell without outflow -> basin size 0'
diag(this_i, this_j) = 0 ! might be a basin of size 0
changed = .true.
cycle
end if
if (diag(this_i,this_j) == -10) then ! found a new cell with outflow
!write(*,*) ' new cell with outflow -> basin',maxbasin
diag(this_i,this_j) = maxbasin
maxbasin = maxbasin+1
changed = .true.
end if
! follow down-river
riverlength = 1
!write(*,'(a,i5,i5,a,i5,i5)') 'follow from ',i,j,' to ',lpj_runoff_coast_map_i(i,j), lpj_runoff_coast_map_j(i,j)
do while (.true.) ! sigh. F90 does not have && shortcut evaluation like C
if (next_i <= 0 .or. next_j <= 0) exit
if (next_i == this_i .and. next_j == this_j) exit
if (riverlength >= global_number_lpj_cells) exit
if (diag(next_i,next_j) == -10 & ! new cell
.or. diag(next_i,next_j) == 0) then ! cell without outflow
diag(next_i,next_j) = diag(this_i,this_j) ! now we found an inflow
!write(*,*) ' found inflow to basin 0 from', diag(this_i,this_j)
changed = .true.
else if (diag(next_i,next_j) == diag(this_i,this_j)) then
!write(*,*) ' stop follow downriver'
exit ! apparently we have been along this path before, no more changes to do downstream
else
! Found inflow into a previously found different basin.
! Propagate the lower basin number.
diag(next_i,next_j) = min(diag(next_i,next_j), diag(this_i,this_j))
diag(this_i,this_j) = diag(next_i,next_j)
!write(*,*) ' connecting basins',diag(next_i,next_j), diag(this_i,this_j)
changed = .true.
end if
riverlength = riverlength + 1
!if (riverlength == global_number_lpj_cells) &
! write(*,*) 'Warning: river routing loop at ',next_i,',',next_j
this_i = next_i ! and move on
this_j = next_j
next_i = lpj_runoff_coast_map_i(this_i,this_j)
next_j = lpj_runoff_coast_map_j(this_i,this_j)
end do
end if
end do
end do
end do
write(*,*) 'maxbasin ', maxbasin
used = send_data(id_runoff_basin, diag(is:ie,js:je))
end if ! id_runoff_basin > 0
end subroutine prepare_coastal_runoff_routing
!
!
! Convert from longitude to grid index
integer function lon2i(lon_in, dlon)
real, intent(in) :: lon_in, dlon
real :: lon
lon = lon_in
if (lon_in < 0.0) lon = lon + 360.0 ! transform -180 .. 180 to 0 .. 360
lon2i = lon / dlon + 1
end function lon2i
!
!
! Convert from latitude to grid index
integer function lat2j(lat, dlat)
real, intent(in) :: lat, dlat
lat2j = (lat+90.0) / dlat + 1
end function lat2j
!
!
! Convert from grid index to longitude
real function i2lon(i_in, dlon)
integer, intent(in) :: i_in
real, intent(in) :: dlon
i2lon = dlon/2.0 + (i_in-1) * dlon
end function i2lon
!
!
! Convert from grid index to longitude
real function j2lat(j_in, dlat)
integer, intent(in) :: j_in
real, intent(in) :: dlat
j2lat = -90.0 + dlat/2.0 + (j_in-1) * dlat
end function j2lat
!
!
!
! Tests whether a point is a coastal point.
!
!
!
! Returns true if and only if point (i,j) is a coastal point on the map gfrac.
!
!
!
! value=is_coastal(i, j, gfrac)
!
logical function is_coastal(i, j, gfrac)
integer, intent(in) :: i, j ! coordinates of the point in question
real, intent(in) :: gfrac(:,:) ! fractional area of the land
integer ip,in, jp,jn ! coordinates of surrounding points
! find coordinates of surrounding points
ip = i-1; if (ip<1) ip = size(gfrac, 1)
in = i+1; if (in>size(gfrac,1)) in = 1
jp = max(1, j-1);
jn = min(size(gfrac,2), j+1)
!write(*,*) 'is_coastal bounds(gfrac) ',lbound(gfrac), ubound(gfrac)
!write(*,'(a,i4,i4,i4,i4,i4,i4)') ' i,j,ip,jp,in,jn ',i,j,ip,jp,in,jn
__ASSERT__(i >= lbound(gfrac,1) .and. i <= ubound(gfrac,1), "i out of range")
__ASSERT__(j >= lbound(gfrac,2) .and. j <= ubound(gfrac,2), "j out of range")
__ASSERT__(ip >= lbound(gfrac,1) .and. ip <= ubound(gfrac,1), "ip out of range")
__ASSERT__(jp >= lbound(gfrac,2) .and. jp <= ubound(gfrac,2), "jp out of range")
__ASSERT__(in >= lbound(gfrac,1) .and. in <= ubound(gfrac,1), "in out of range")
__ASSERT__(jn >= lbound(gfrac,2) .and. jn <= ubound(gfrac,2), "jn out of range")
is_coastal = gfrac(i,j) < (1.0-min_water_frac) .and.( &
gfrac(i ,j ) > min_land_frac .or. &
gfrac(in,j ) > min_land_frac .or. &
gfrac(ip,j ) > min_land_frac .or. &
gfrac(i ,jn) > min_land_frac .or. &
gfrac(i ,jp) > min_land_frac )
end function is_coastal
!
! ==============================================================================
! Reports error, including file name and line.
subroutine my_error(mod_name, message, mode, file, line)
character(len=*), intent(in) :: mod_name
character(len=*), intent(in) :: message
integer, intent(in) :: mode
character(len=*), intent(in), optional :: file
integer, intent(in), optional :: line
! ---- local vars ----------------------------------------------------------
character(len=512) :: mesg
if(present(file)) then ! assume that file and line are either both present or not
write(mesg,'("File ",a," Line ",i4.4," :: ",a)')&
file, line, trim(message)
else
mesg = trim(message)
endif
call error_mesg(mod_name, mesg, mode)
end subroutine
end module lpj_driver_mod
#endif