!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! !! !! 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 !! !! under the terms of the GNU General Public License as published by !! !! the Free Software Foundation, either version 3 of the License, or !! !! (at your option) any later version. !! !! !! !! 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, see: http://www.gnu.org/licenses/gpl.txt !! !! !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! module atmosphere_mod use fms_mod, only: set_domain, write_version_number, field_size, file_exist, & mpp_pe, mpp_root_pe, error_mesg, FATAL, read_data, write_data, nullify_domain use constants_mod, only: grav, pi use transforms_mod, only: trans_grid_to_spherical, trans_spherical_to_grid, & get_deg_lon, get_deg_lat, get_grid_boundaries, grid_domain, & spectral_domain, get_grid_domain, get_lon_max, get_lat_max, atmosphere_domain use time_manager_mod, only: time_type, get_time, operator( + ) use press_and_geopot_mod, only: compute_pressures_and_heights use spectral_dynamics_mod, only: spectral_dynamics_init, spectral_dynamics, spectral_dynamics_end, & get_num_levels, get_axis_id, spectral_diagnostics, get_initial_fields, & complete_robert_filter, get_surf_geopotential use tracer_type_mod, only: tracer_type use field_manager_mod, only: MODEL_ATMOS use tracer_manager_mod, only: get_number_tracers #ifdef IDEALIZED_MOIST_MODEL use idealized_moist_phys_mod, only: idealized_moist_phys_init , idealized_moist_phys , idealized_moist_phys_end #else use hs_forcing_mod, only: hs_forcing_init, hs_forcing #endif implicit none private !================================================================================================================================= character(len=128) :: version= & '$Id: atmosphere.F90,v 20.0 2013/12/13 23:25:37 fms Exp $' character(len=128) :: tagname= & '$Name: tikal $' character(len=10), parameter :: mod_name='atmosphere' !================================================================================================================================= public :: atmosphere_init, atmosphere, atmosphere_end, atmosphere_domain !================================================================================================================================= integer, parameter :: num_time_levels = 2 integer :: is, ie, js, je, num_levels, num_tracers, nhum logical :: dry_model real, allocatable, dimension(:,:,:,:) :: p_half, p_full real, allocatable, dimension(:,:,:,:) :: z_half, z_full type(tracer_type), allocatable, dimension(:) :: tracer_attributes real, allocatable, dimension(:,:,:,:,:) :: grid_tracers real, allocatable, dimension(:,:,: ) :: psg, wg_full real, allocatable, dimension(:,:,:,: ) :: ug, vg, tg real, allocatable, dimension(:,: ) :: dt_psg real, allocatable, dimension(:,:,: ) :: dt_ug, dt_vg, dt_tg real, allocatable, dimension(:,:,:,:) :: dt_tracers real, allocatable, dimension(:) :: deg_lon, deg_lat real, allocatable, dimension(:,:) :: rad_lon_2d, rad_lat_2d real, allocatable, dimension(:,:) :: surf_geopotential real, allocatable, dimension(:,:) :: rad_lonb_2d, rad_latb_2d real, allocatable, dimension(:) :: rad_lonb, rad_latb integer :: previous, current, future logical :: module_is_initialized =.false. integer :: dt_integer real :: dt_real type(time_type) :: Time_step !================================================================================================================================= contains !================================================================================================================================= subroutine atmosphere_init(Time_init, Time, Time_step_in) type (time_type), intent(in) :: Time_init, Time, Time_step_in integer :: seconds, days, lon_max, lat_max, ntr, nt, i, j integer, dimension(4) :: siz real, dimension(2) :: time_pointers character(len=64) :: file, tr_name character(len=256) :: message if(module_is_initialized) return call write_version_number(version, tagname) Time_step = Time_step_in call get_time(Time_step, seconds, days) dt_integer = 86400*days + seconds dt_real = float(dt_integer) call get_number_tracers(MODEL_ATMOS, num_prog=num_tracers) allocate (tracer_attributes(num_tracers)) call spectral_dynamics_init(Time, Time_step, tracer_attributes, dry_model, nhum) call get_grid_domain(is, ie, js, je) call get_num_levels(num_levels) allocate (p_half (is:ie, js:je, num_levels+1, num_time_levels)) allocate (z_half (is:ie, js:je, num_levels+1, num_time_levels)) allocate (p_full (is:ie, js:je, num_levels, num_time_levels)) allocate (z_full (is:ie, js:je, num_levels, num_time_levels)) allocate (wg_full (is:ie, js:je, num_levels)) allocate (psg (is:ie, js:je, num_time_levels)) allocate (ug (is:ie, js:je, num_levels, num_time_levels)) allocate (vg (is:ie, js:je, num_levels, num_time_levels)) allocate (tg (is:ie, js:je, num_levels, num_time_levels)) allocate (grid_tracers (is:ie, js:je, num_levels, num_time_levels, num_tracers )) allocate (dt_psg (is:ie, js:je)) allocate (dt_ug (is:ie, js:je, num_levels)) allocate (dt_vg (is:ie, js:je, num_levels)) allocate (dt_tg (is:ie, js:je, num_levels)) allocate (dt_tracers (is:ie, js:je, num_levels, num_tracers)) allocate (deg_lon (is:ie )) allocate (rad_lon_2d (is:ie, js:je)) allocate (deg_lat ( js:je)) allocate (rad_lat_2d (is:ie, js:je)) allocate (rad_lonb_2d(is:ie+1, js:je+1)) allocate (rad_latb_2d(is:ie+1, js:je+1)) allocate (rad_lonb (is:ie+1)) allocate (rad_latb (js:je+1)) p_half = 0.; z_half = 0.; p_full = 0.; z_full = 0. wg_full = 0.; psg = 0.; ug = 0.; vg = 0.; tg = 0.; grid_tracers = 0. dt_psg = 0.; dt_ug = 0.; dt_vg = 0.; dt_tg = 0.; dt_tracers = 0. allocate (surf_geopotential(is:ie, js:je)) call get_surf_geopotential(surf_geopotential) !-------------------------------------------------------------------------------------------------------------------------------- file = 'INPUT/atmosphere.res.nc' if(file_exist(trim(file))) then call get_lon_max(lon_max) call get_lat_max(lat_max) call field_size(trim(file), 'ug', siz) if(lon_max /= siz(1) .or. lat_max /= siz(2)) then write(message,*) 'Resolution of restart data does not match resolution specified on namelist. Restart data: lon_max=', & siz(1),', lat_max=',siz(2),' Namelist: lon_max=',lon_max,', lat_max=',lat_max call error_mesg('atmosphere_init', message, FATAL) endif call nullify_domain() ! call read_data(trim(file), 'previous', previous) ! No interface of read_data exists to read integer scalars ! call read_data(trim(file), 'current', current) ! No interface of read_data exists to read integer scalars call read_data(trim(file), 'time_pointers', time_pointers) ! Getaround for no interface to read integer scalars previous = int(time_pointers(1)) ! Getaround for no interface to read integer scalars current = int(time_pointers(2)) ! Getaround for no interface to read integer scalars do nt=1,num_time_levels call read_data(trim(file), 'ug', ug(:,:,:,nt), grid_domain, timelevel=nt) call read_data(trim(file), 'vg', vg(:,:,:,nt), grid_domain, timelevel=nt) call read_data(trim(file), 'tg', tg(:,:,:,nt), grid_domain, timelevel=nt) call read_data(trim(file), 'psg', psg(:,:, nt), grid_domain, timelevel=nt) do ntr = 1,num_tracers tr_name = trim(tracer_attributes(ntr)%name) call read_data(trim(file), trim(tr_name), grid_tracers(:,:,:,nt,ntr), grid_domain, timelevel=nt) enddo ! end loop over tracers enddo ! end loop over time levels call read_data(trim(file), 'wg_full', wg_full, grid_domain) else previous = 1; current = 1 call get_initial_fields(ug(:,:,:,1), vg(:,:,:,1), tg(:,:,:,1), psg(:,:,1), grid_tracers(:,:,:,1,:)) endif !-------------------------------------------------------------------------------------------------------------------------------- if(dry_model) then call compute_pressures_and_heights(tg(:,:,:,current), psg(:,:,current), surf_geopotential, & z_full(:,:,:,current), z_half(:,:,:,current), p_full(:,:,:,current), p_half(:,:,:,current)) call compute_pressures_and_heights(tg(:,:,:,previous), psg(:,:,previous), surf_geopotential, & z_full(:,:,:,previous), z_half(:,:,:,previous), p_full(:,:,:,previous), p_half(:,:,:,previous)) else call compute_pressures_and_heights(tg(:,:,:,current), psg(:,:,current), surf_geopotential, & z_full(:,:,:,current), z_half(:,:,:,current), p_full(:,:,:,current), p_half(:,:,:,current), & grid_tracers(:,:,:,current,nhum)) call compute_pressures_and_heights(tg(:,:,:,previous), psg(:,:,previous), surf_geopotential, & z_full(:,:,:,previous), z_half(:,:,:,previous), p_full(:,:,:,previous), p_half(:,:,:,previous), & grid_tracers(:,:,:,previous,nhum)) endif call get_deg_lon(deg_lon) do i=is,ie rad_lon_2d(i,:) = deg_lon(i)*pi/180. enddo call get_deg_lat(deg_lat) do j=js,je rad_lat_2d(:,j) = deg_lat(j)*pi/180. enddo call get_grid_boundaries(rad_lonb, rad_latb) do i = is,ie+1 rad_lonb_2d(i,:) = rad_lonb(i) enddo do j = js,je+1 rad_latb_2d(:,j) = rad_latb(j) enddo #ifdef IDEALIZED_MOIST_MODEL call idealized_moist_phys_init(Time, Time_step, nhum, rad_lat_2d, rad_lonb_2d, rad_latb_2d, tg(:,:,num_levels,current)) #else call hs_forcing_init(get_axis_id(), Time, rad_lonb_2d, rad_latb_2d) #endif module_is_initialized = .true. return end subroutine atmosphere_init !================================================================================================================================= subroutine atmosphere(Time) type(time_type), intent(in) :: Time real :: delta_t type(time_type) :: Time_next if(.not.module_is_initialized) then call error_mesg('atmosphere','atmosphere module is not initialized',FATAL) endif dt_ug = 0.0 dt_vg = 0.0 dt_tg = 0.0 dt_psg = 0.0 dt_tracers = 0.0 if(current == previous) then delta_t = dt_real else delta_t = 2*dt_real endif Time_next = Time + Time_step #ifdef IDEALIZED_MOIST_MODEL call idealized_moist_phys(Time, p_half, p_full, z_half, z_full, ug, vg, tg, grid_tracers, & previous, current, dt_ug, dt_vg, dt_tg, dt_tracers) #else call hs_forcing(1, ie-is+1, 1, je-js+1, delta_t, Time_next, rad_lon_2d, rad_lat_2d, & p_half(:,:,:,current ), p_full(:,:,:,current ), & ug(:,:,:,previous), vg(:,:,:,previous ), & tg(:,:,:,previous), grid_tracers(:,:,:,previous,:), & ug(:,:,:,previous), vg(:,:,:,previous ), & tg(:,:,:,previous), grid_tracers(:,:,:,previous,:), & dt_ug(:,:,: ), dt_vg(:,:,: ), & dt_tg(:,:,: ), dt_tracers(:,:,:,:)) #endif if(previous == current) then future = num_time_levels + 1 - current else future = previous endif call spectral_dynamics(Time, psg(:,:,future), ug(:,:,:,future), vg(:,:,:,future), & tg(:,:,:,future), tracer_attributes, grid_tracers(:,:,:,:,:), future, & dt_psg, dt_ug, dt_vg, dt_tg, dt_tracers, wg_full, & p_full(:,:,:,current), p_half(:,:,:,current), z_full(:,:,:,current)) call complete_robert_filter(tracer_attributes) if(dry_model) then call compute_pressures_and_heights(tg(:,:,:,future), psg(:,:,future), surf_geopotential, & z_full(:,:,:,future), z_half(:,:,:,future), p_full(:,:,:,future), p_half(:,:,:,future)) else call compute_pressures_and_heights(tg(:,:,:,future), psg(:,:,future), surf_geopotential, & z_full(:,:,:,future), z_half(:,:,:,future), p_full(:,:,:,future), p_half(:,:,:,future), & grid_tracers(:,:,:,future,nhum)) endif call spectral_diagnostics(Time_next, psg(:,:,future), ug(:,:,:,future), vg(:,:,:,future), & tg(:,:,:,future), wg_full, grid_tracers(:,:,:,:,:), future) previous = current current = future return end subroutine atmosphere !================================================================================================================================= subroutine atmosphere_end integer :: ntr, nt character(len=64) :: file, tr_name if(.not.module_is_initialized) return file='RESTART/atmosphere.res' call nullify_domain() call write_data(trim(file), 'time_pointers', (/real(previous),real(current)/)) do nt=1,num_time_levels call write_data(trim(file), 'ug', ug(:,:,:,nt), grid_domain) call write_data(trim(file), 'vg', vg(:,:,:,nt), grid_domain) call write_data(trim(file), 'tg', tg(:,:,:,nt), grid_domain) call write_data(trim(file), 'psg', psg(:,:, nt), grid_domain) do ntr = 1,num_tracers tr_name = trim(tracer_attributes(ntr)%name) call write_data(trim(file), tr_name, grid_tracers(:,:,:,nt,ntr), grid_domain) enddo enddo call write_data(trim(file), 'wg_full', wg_full, grid_domain) deallocate (p_half, z_half, p_full, z_full, wg_full, psg, ug, vg, tg, grid_tracers) deallocate (dt_psg, dt_ug, dt_vg, dt_tg, dt_tracers) deallocate (deg_lon, rad_lon_2d, deg_lat, rad_lat_2d) call set_domain(grid_domain) #ifdef IDEALIZED_MOIST_MODEL call idealized_moist_phys_end #endif call spectral_dynamics_end(tracer_attributes) deallocate(tracer_attributes) module_is_initialized = .false. end subroutine atmosphere_end !================================================================================================================================= end module atmosphere_mod