! Stefan Petri ! ! ! ! !----------------------------------------------------------------------- ! ! Standalone main program for testing the interface between FMS/F90 ! and AtmosphereIII/C++ ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! !! !! 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 !! !! !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! program testaeolus use sinsol_mod, only: sinsol_mod_init !----------------------------------------------------------------------- implicit none !----------------------------------------------------------------------- ! type surf_diff_type real, pointer, dimension(:,:) :: dtmass => NULL() ! dt/mass, where dt = atmospheric time step (sec) ! mass = mass of lowest atmospheric layer (Kg/m2) real, pointer, dimension(:,:) :: dflux_t => NULL() ! derivative of the temperature flux at the top of the lowest ! atmospheric layer with respect to the temperature ! of that layer (J/(m2 K)) real, pointer, dimension(:,:,:) :: dflux_tr => NULL() ! derivative of the flux of specific humidity ! at the top of the lowest atmospheric layer with respect to ! the specific humidity of that layer (--/(m2 K)) real, pointer, dimension(:,:) :: delta_t => NULL() ! the increment in temperature in the lowest atmospheric ! layer (((i+1)-(i-1) if atmos model is leapfrog) (K) ! (defined in gcm_vert_diff_down as the increment computed up ! to this point in model, including effect of vertical ! diffusive flux at top of lowest model level, presumed ! to be modified to include effects of surface fluxes ! outside of this module, then used to start upward ! tridiagonal sweep, real, pointer, dimension(:,:,:):: delta_tr => NULL() ! similarly for the increment in specific humidity ! (non-dimensional = Kg/Kg) real, pointer, dimension(:,:) :: delta_u => NULL() real, pointer, dimension(:,:) :: delta_v => NULL() real, pointer, dimension(:,:) :: sst_miz => NULL() end type surf_diff_type ! !----------------------------------------------------------------------- character(len=128) :: version = '$Id$' character(len=128) :: tagname = '$Name: mom4p1_pubrel_dec2009_nnz $' !---- namelist (saved in file input.nml) ---- ! ! ! Longitudinal resolution of atmosphere grid (at the equator). ! Must match the data in the grid specification file. ! ! ! Latitudinal resolution of atmosphere grid. ! Must match the data in the grid specification file. ! ! ! Number of vertical layers of atmosphere grid. ! ! ! If true, longitudinal adjacent grid cells are merged towards the poles, ! so that cell sizes, and with them the time steps, dont get too small. ! If false, simply use a regular lat/lon grid. ! ! namelist /atmosphere_nml/ & NoNx, NoNy, NoNz, NoNz_Strato, & ! Number of grid cells in lon, lat, height (Tropo+Strato) polar_merging, & ! merge cells longitudinal towards pole, or regular lat/lon grid? lon_0_360 ! whether longitudes run 0 to 360 in the grid spec file !----------------------------------------------------------------------- !---- private data ---- integer :: NoNx = 96, NoNy = 48, NoNz = 5, NoNz_Strato = 0 ! values from namelist logical :: polar_merging = .true. , lon_0_360 = .true. integer :: nr_tracers = 1 ! total number of tracers in surf_diff_type fields integer :: q_ind = 1 ! index of humidity in tracers array - for now hardcoded! integer :: is, ie, js, je ! local domain boundary indices !------------------------- axis names ---------------------------------- character(len=8) :: axiset = 'aeolus' character(len=8) :: mod_name = 'aeolus' !----------------------------------------------------------------------- type(surf_diff_type) :: Surf_diff integer :: Time_init_days, Time_init_seconds integer :: Time_days, Time_seconds integer :: Time_step_days, Time_step_seconds integer :: unit, ierr, io integer :: log_unit integer :: pe, root_pe, npes, layout(2) real(kind=8), dimension(:), allocatable :: lonb, latb, lon, lat ! for the axes integer :: id_lonb, id_latb, id_lon, id_lat integer, dimension(1:2) :: axes_2d character(len=64) :: filename ! 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 :: deadbeef, fse; real(8) :: onetwothree, pitest logical :: true = .true., false = .false. integer :: kind_i = KIND(deadbeef), kind_r = KIND(pitest), kind_l = KIND(true) real(8), dimension(:,:,:), allocatable :: testarray3d real(8), dimension(:,:), allocatable :: & solarm, & ! daily solar radiation at the top of the atmosphere coszm ! daily averaged solar zenit angle integer :: DaysPerYear integer :: i, j, k ! start of program code deadbeef = X'eadbeef' ! X'deadbeef' yields overflow error on 32bit machines fse = 4711 onetwothree = 1234567809.0987654321D0 pitest = 3.1415926535897932384626433832795029D0 ! from linux is = 1 ie = NoNx js = 1 je = NoNy !----- read namelist ----- open(file='input.nml', unit=42) ierr=1; do while (ierr /= 0) read (unit=42, nml=atmosphere_nml, iostat=io, end=10) !ierr = check_nml_error (io, 'atmosphere_nml') enddo 10 close(42) !----- write version and namelist to log file ----- write (*, nml=atmosphere_nml) !---- compute physics/atmos time step in seconds ---- ! initialize domain decomposition pe = 0 root_pe = 0 npes = 1 !call mpp_get_compute_domain(grid_domain, is, ie, js, je) ! --- now we can allocate variables --- ! hack: in case Aeolus uses more than 1 tracer other than water vapor, revamp this code. allocate (Surf_diff%dtmass (is:ie, js:je) ) allocate (Surf_diff%dflux_t (is:ie, js:je) ) allocate (Surf_diff%dflux_tr (is:ie, js:je, q_ind) ) allocate (Surf_diff%delta_t (is:ie, js:je) ) allocate (Surf_diff%delta_tr (is:ie, js:je, q_ind) ) allocate (Surf_diff%delta_u (is:ie, js:je) ) allocate (Surf_diff%delta_v (is:ie, js:je) ) Surf_diff%dtmass = 0.0 Surf_diff%dflux_t = 0.0 Surf_diff%dflux_tr = 0.0 Surf_diff%delta_t = 0.0 Surf_diff%delta_tr = 0.0 Surf_diff%delta_u = 0.0 Surf_diff%delta_v = 0.0 Time_init_seconds = 0 Time_init_days = 0 Time_seconds = 7200 Time_days = 42 Time_step_seconds = 3600 Time_step_days = 0 allocate(testarray3d(2,3,4)) do i=1,2 do j=1,3 do k=1,4 testarray3d(i,j,k) = i*100.0+j+k/10.0 write(*,*) 'testarray3d', i, j, k, testarray3d(i, j, k) end do end do end do ! now call out into the C++ - implemented heart of Aeolus write(*,*) 'MAXEXPONENT ', MAXEXPONENT(pitest) call aeolus_init( & deadbeef, pitest, & true, false, & testarray3d, & kind_i, kind_r, kind_l, & Time_init_days, Time_init_seconds, & Time_days, Time_seconds, & Time_step_days, Time_step_seconds, & NoNx, NoNy, NoNz, NoNz_Strato, & polar_merging, & lon_0_360, & pe, root_pe, npes, & is, ie, js, je, & nr_tracers, q_ind, & onetwothree, fse & ) !----- initialize atmos_axes and diagnostics allocate(lonb(1:NoNx+1), latb(1:NoNy+1), lon(1:NoNx), lat(1:NoNy)) call aeolus_get_axes_coords(lonb, latb, lon, lat) write (*,*) 'lonb', lonb write (*,*) 'latb', latb write (*,*) 'lon', lon write (*,*) 'lat', lat ! now init the tables of cosine zenith angle and daily insolation DaysPerYear = 365 allocate(solarm(NoNy, DaysPerYear)) allocate(coszm(NoNy, DaysPerYear)) call sinsol_mod_init(real(Time_init_days / DaysPerYear, 8), DaysPerYear, NoNy, lat, solarm, coszm) deallocate(lonb) deallocate(latb) deallocate(lon) deallocate(lat) call test_inner(is,ie,js,je, Surf_diff) call aeolus_finish(Time_seconds, Time_days) contains subroutine test_inner(is, ie, js, je, Surf_diff) integer, intent(in) :: is, ie, js, je type(surf_diff_type), intent(inout) :: Surf_diff integer :: i, j ! for ``down'' test real(8), dimension(is:ie,js:je) :: frac_land, t_surf, albedo real(8), dimension(is:ie,js:je) :: albedo_vis_dir, albedo_nir_dir real(8), dimension(is:ie,js:je) :: albedo_vis_dif, albedo_nir_dif real(8), dimension(is:ie,js:je) :: rough_mom real(8), dimension(is:ie,js:je) :: u_star, b_star, q_star, dtau_du, dtau_dv real(8), dimension(is:ie,js:je) :: u_flux, v_flux real(8), dimension(is:ie,js:je) :: gust, coszen real(8), dimension(is:ie,js:je) :: flux_sw real(8), dimension(is:ie,js:je) :: flux_sw_dir, flux_sw_dif, flux_sw_down_vis_dir real(8), dimension(is:ie,js:je) :: flux_sw_down_vis_dif, flux_sw_down_total_dir real(8), dimension(is:ie,js:je) :: flux_sw_down_total_dif, flux_sw_vis real(8), dimension(is:ie,js:je) :: flux_sw_vis_dir, flux_sw_vis_dif real(8), dimension(is:ie,js:je) :: flux_lw ! for ``up'' test !real, dimension(is:ie,js:je) :: frac_land !type(surf_diff_type) :: Surf_diff real(8), dimension(is:ie,js:je) :: lprec, fprec !, gust !real, dimension(:,:) :: u_star, b_star, q_star ! for cell_area test real(8), dimension(is:ie,js:je) :: area_out ! for boundary test real(8) :: lon_boundaries(is:ie+1,js:je+1), lat_boundaries(is:ie+1,js:je+1) ! radians !! ! for bottom_mass test real(8), dimension(is:ie,js:je) :: t_bot, p_bot, z_bot, p_surf, slp real(8), dimension(is:ie,js:je,3) :: tr_bot ! for wind test real(8), dimension(is:ie,js:je) :: u_bot, v_bot ! for stock test integer :: idx ! stock type index real(8) :: val do i=is,ie do j=js,je frac_land(i,j) = i+j/100.0 end do end do ! get info from SphericalGrid call aeolus_get_boundary(lon_boundaries, lat_boundaries, .true., & lbound(lon_boundaries, dim=1), ubound(lon_boundaries, dim=1), & lbound(lon_boundaries, dim=2), ubound(lon_boundaries, dim=2), & lbound(lat_boundaries, dim=1), ubound(lat_boundaries, dim=1), & lbound(lat_boundaries, dim=2), ubound(lat_boundaries, dim=2)) write(*,*) 'lon_boundaries ', lon_boundaries write(*,*) 'lat_boundaries ', lat_boundaries write(*,*) 'end boundaries' ! get area information from SphericalGrid call aeolus_cell_area(area_out) !do xsize=is,ie ! do ysize=js,je ! write(*,*) 'i ', xsize, ' j ', ysize, ' area ', area_out(xsize, ysize) ! end do !end do call aeolus_test_var_access( & frac_land, & t_bot, tr_bot, p_bot) do i=is,ie do j=js,je write(*,*) 't_bot ', i, j, t_bot(i,j) write(*,*) 'p_bot ', i, j, p_bot(i,j) write(*,*) 'tr_bot', i, j, tr_bot(i,j,1) end do end do call flush(6) !call aeolus_atmosphere_down ( & ! Time_seconds, Time_days, & ! frac_land, t_surf, albedo, & ! albedo_vis_dir, albedo_nir_dir, & ! albedo_vis_dif, albedo_nir_dif, & ! rough_mom, & ! u_star, b_star, q_star, & ! dtau_du, dtau_dv, u_flux, v_flux, & ! gust, coszen, flux_sw, & ! flux_sw_dir, flux_sw_dif, & ! flux_sw_down_vis_dir, & ! flux_sw_down_vis_dif, & ! flux_sw_down_total_dir, & ! flux_sw_down_total_dif, & ! flux_sw_vis, & ! flux_sw_vis_dir, & ! flux_sw_vis_dif, & ! flux_lw, & ! Surf_diff%delta_tr, & ! Surf_diff%dflux_tr, & ! Surf_diff%delta_t, & ! Surf_diff%dflux_t, & ! Surf_diff%dtmass, & ! Surf_diff%delta_u, & ! Surf_diff%delta_v, & ! ) !call aeolus_atmosphere_up ( & ! Time_seconds, Time_days, & ! frac_land, & ! Surf_diff%delta_tr, & ! Surf_diff%dflux_tr, & ! lprec, fprec, gust, & ! u_star, b_star, q_star ) !call aeolus_get_bottom_mass(t_bot, tr_bot, p_bot, z_bot, p_surf, slp) !call aeolus_get_bottom_wind(u_bot, v_bot) !call aeolus_get_stock_pe(idx, val) end subroutine test_inner end program testaeolus