!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! !! !! 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 !! !! !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! $Id: fv_pack.F90,v 17.0 2009/07/21 02:53:12 fms Exp $ module fv_pack use timingModule #ifdef USE_LIMA use shr_kind_mod, only: r8 => shr_kind_r8 #endif !----------------------------------------------------------------------- ! "SPMD" stands for Single Program Multiple Data. It activates ! domain decomposition with inter-process communication done ! by MPI-1, MPI-2, or SHMEM -- this is the first level of the ! parallelism. The second level parallelism uses the standard ! "OpenMP". To activate OpenMP, add "-mp" to the compiler and the load ! flags on the SGI. !----------------------------------------------------------------------- #ifdef SPMD use mod_comm, only: gid, mp_init, mp_exit, y_decomp, & mp_send3d, mp_recv3d, & mp_scatter3d, mp_send4d_ns, mp_recv4d_ns #endif ! --- FMS ---- use fms_mod, only: file_exist, open_namelist_file, & error_mesg, FATAL, & close_file, check_nml_error, & write_version_number, stdlog, stdout, stderr use mpp_mod, only: mpp_pe, mpp_root_pe, mpp_npes, MPP_VERBOSE, & mpp_get_current_pelist, mpp_sync, mpp_error use mpp_domains_mod, only: mpp_define_domains, mpp_domains_init, domain2D, & mpp_get_compute_domain use constants_mod, only: radius, omega, grav, kappa, rdgas, cp_air, rvgas ! --- tracer manager --- use tracer_manager_mod, only : tm_get_number_tracers => get_number_tracers, & tm_get_tracer_index => get_tracer_index, & tm_get_tracer_indices => get_tracer_indices, & tm_set_tracer_profile => set_tracer_profile, & tm_get_tracer_names => get_tracer_names, & tm_check_if_prognostic=> check_if_prognostic,& tm_register_tracers => register_tracers use field_manager_mod, only : MODEL_ATMOS use fv_arrays_mod, only: fv_arrays_init, fv_arrays_exit, fv_array_check, & fv_stack_push, fv_array_sync, isp, iep, jsp, jep, ksp, kep, layout, & nlev, ncnst, n_zonal, fv_print_chksums, fv_pset_stack #ifndef USE_LIMA use pft_module, only: pft_init, pft_end #else use pft_module, only: pft_init #endif ! --- FMS ---- implicit none private ! --- FMS ---- !================== ! fv_core namelist: !================== !---------------- ! Dimensionality: !---------------- integer nlon ! E-W (total number of longitudes) integer mlat ! N-S (total number of latitudes) ! integer nlev ! Vertical dimension ! integer ncnst ! Total number of tracers integer :: nt_phys = 4 ! Total number of tracers needed by physics integer :: nt_prog = 0 ! Number of prognostic tracers integer :: pnats = 0 ! Number of non-advected tracers; default to 0 logical :: do_fms_tracer_manager = .FALSE. ! get list of tracers from field_table ! name of fms tracers restart file character(len=*), parameter :: fms_tracers_file='atmos_tracers.res' !------------------------------------------ ! Parameters that control the FV algorithms: !------------------------------------------ ! Momentum transport options: integer :: iord_mt = 4 ! E-W integer :: jord_mt = 4 ! N-S integer :: kord_mt = 4 ! vertical mapping option ! Heat transport options: integer :: iord_tm = 4 ! E-W integer :: jord_tm = 4 ! N-S integer :: kord_tm = 7 ! vertical mapping option ! Tracer transport options: integer :: iord_tr = 4 ! E-W integer :: jord_tr = 4 ! N-S integer :: kord_tr = 7 ! vertical mapping option for tracers integer :: n_split = 0 ! Number of time splits for the lagrangian dynamics ! Default = 0 (automatic computation of best value) integer :: n_spong = 1 ! Number of sponge layers from the top integer :: n_diffu = 0 ! Number of diffusive layers from the bottom ! integer :: n_zonal = 0 ! loop decomposition in East-West direction #ifndef USE_LIMA integer :: dyn_step= 0 ! Multi-Dimensional transport algorithms for dynamics integer :: tra_step= 1 ! Multi-Dimensional transport algorithms for tracers ! 0: NO directional split ! 1: start with meridional transport ! 2: start with zonal transport #endif logical :: age_tracer= .false. ! Carry the age tracer as the last tracer logical :: adiabatic = .false. ! Run without physics (full or idealized). ! this is mainly to turn off the Held-Suarez ! forcing in the "solo" mode logical :: full_phys = .true. ! Run with full physics with virtual effects ! for R (gas constant) and Cp real:: consv_te = 0. ! Total energy fixer. Range[0,1.]. For example 0.75 ! would provide 75% correction ! Divergence damping parameters: D = a_div + b_div * cos(lat) real:: a_div = 0.6 ! real:: b_div = 0.6 ! #ifndef USE_LIMA !!!1 real:: dry_mass = 98290. #endif #ifdef MARS_GCM real :: p_ref = 700. real :: fv_sponge_damp= 10.0 real :: fv_sponge_lev= 0.5 ! level (mb) for enhanced divergence damping real :: dry_mass = 7.01E2 #else real :: p_ref = 1.E5 real :: dry_mass = 98290. #endif logical :: cold_start = .false. ! read in initial condition logical :: change_time = .false. ! ignore time stamp in fv_rst.res logical :: adjust_dry_mass = .false. ! adjust initial total dry air mass logical :: use_set_eta = .false. ! use set_eta(), instead of restart, for (ak,bk) logical :: pft_phys = .true. ! polar filter physical tendencies logical :: use_tendency = .false. ! use the tendency approach logical :: do_ch4_chem = .false. ! Relax H2o to Haloe data between 1 and 10 mb integer :: print_freq = 0 ! Print max/min of selected fields ! 0: off ! positive n: every n hours ! negative n: every time step #ifdef SW_DYN logical :: fill = .false. #else logical :: fill = .true. ! Activate a simple filling algorithm for tracers #endif integer :: map_dt ! Remapping time step (seconds); default to model ! time step ! Default case for the shallow water test: Rossby wave-4 integer :: icase = 6 ! Test case number if "SW_DYN" is defined ! see tools/nit_sw_ic.f90 for details character(len=24) :: restart_format = 'native' ! native or netcdf character(len=128) :: version = '$Id: fv_pack.F90 1.1.2.15.2.2' character(len=128) :: tagname = '$Name: mom4p1_pubrel_dec2009_nnz $' #ifndef USE_LIMA namelist /fv_core_nml/nlon, mlat, nlev, & ncnst, nt_phys, pnats, & iord_mt, iord_tm, iord_tr, & jord_mt, jord_tm, jord_tr, & kord_mt, kord_tm, kord_tr, & n_split, n_zonal, map_dt, & consv_te, restart_format, adiabatic, full_phys, & fill, adjust_dry_mass, n_spong, n_diffu, & change_time, a_div, b_div, use_set_eta, & use_tendency, do_ch4_chem, pft_phys, age_tracer, & tra_step, dyn_step, dry_mass, print_freq, icase, & #ifdef MARS_GCM p_ref, fv_sponge_lev, fv_sponge_damp, & #endif MARS_GCM layout, fv_pset_stack #else namelist /fv_core_nml/nlon, mlat, nlev, & ncnst, nt_phys, pnats, & iord_mt, iord_tm, iord_tr, & jord_mt, jord_tm, jord_tr, & kord_mt, kord_tm, kord_tr, & n_split, n_zonal, map_dt, & consv_te, restart_format, adiabatic, full_phys, & fill, adjust_dry_mass, n_spong, n_diffu, & change_time, a_div, b_div, use_set_eta, & use_tendency, do_ch4_chem, pft_phys, age_tracer, & #ifdef MARS_GCM p_ref, fv_sponge_lev, fv_sponge_damp, & #endif MARS_GCM print_freq, icase, layout #endif logical mountain ! If terrain is present integer nqrst ! Total number of tracers in fv_rst file !----------------------------------------------------------------- ! Domain decomposition parameters (not to be changed by the user): !----------------------------------------------------------------- integer beglon ! beginning lon index integer endlon ! ending lon index integer beglat ! beginning lat index integer endlat ! ending lat index integer ng_d ! algorithm dependent ghost zones integer ng_s ! algorithm dependent ghost zones for u-wind logical master ! If "master process"? used mostly for print ! and read/write restart ! FMS specific: type(domain2D), save :: fv_domain !----------------------------------------------------------------------- ! Eulerian vertical coordinate: ! pe(k) = ak(k) + bk(k) * ps !----------------------------------------------------------------------- real, allocatable :: ak(:) real, allocatable :: bk(:) integer ks ! number of layers in pure p region real ptop ! pressure at model top (pascal) ! Horizontal grid geometry #ifndef USE_LIMA real, allocatable:: f_d(:,:) ! Coriolis parm at cell center #else real, allocatable:: f_d(:) ! Coriolis parm at cell center #endif real, allocatable:: dx(:) ! grid distance (m) along latitude circles real, allocatable:: one_by_dx(:) ! 1 / dx real one_by_dy ! 1 / dy real acap ! scaled polar area = nlon * cosp(1) real rcap ! 1 / acap !----------------------------------- ! Geometric arrays !----------------------------------- real, allocatable :: lon(:), lonb(:), rlonb(:) real, allocatable :: lat(:), latb(:), rlatb(:) real, allocatable :: rlat(:,:) real, allocatable :: rlon(:,:) real, allocatable :: area(:,:) ! grid cell area: m**2 real, allocatable :: sine(:) real, allocatable :: cose(:) real, allocatable :: sinp(:) real, allocatable :: cosp(:) real, allocatable :: acosp(:) real, allocatable :: acosu(:) ! for u-wind averaging real, allocatable :: grid_weight(:) real, allocatable :: coslon(:) real, allocatable :: sinlon(:) real, allocatable :: cosl5(:) real, allocatable :: sinl5(:) real, parameter::cv_wv = 1463. ! specific heat of water vapor at constant ! volume [J/(K KG)] real, parameter::cp_wv = 1952. ! specific heat of water vapor at constant ! pressure [J/(K KG)] real cp_vir real pi #ifdef USE_LIMA ! Polar filter related integer ifax(13) !ECMWF fft real(r8), allocatable :: trigs(:) real(r8), allocatable :: dc(:,:), de(:,:), sc(:), se(:) integer:: ipft1, ipft2 #endif !----------------------------------------------------------------------- ! Others: !----------------------------------------------------------------------- ! real, allocatable :: pesouth(:,:) ! "pe at the south" for SPMD logical :: u_ghosted = .false. ! u-wind is ghosted at the first lat to the N real, parameter::ptop_min = 1.E-6 ! minimum pressure (pa) at model top to avoid ! floating point exception; this is not needed ! if model top is not at zero real, parameter :: too_big = 1.E30 ! use this for missing values real :: age_time ! Time (sec) since init of age_tracer #ifdef USE_LIMA integer:: js2g0, jn1g1 #endif integer:: ighost !PUBLIC VARIABLES public :: nlon, mlat, nlev, beglon, endlon, beglat, endlat, & age_tracer, age_time, do_ch4_chem, u_ghosted, ng_d, ng_s, & use_tendency, pft_phys, ak, bk, get_eta_level, ptop, full_phys, & coslon, sinlon, ighost, cosp, sinp, cose, sine, master, & nt_phys, nt_prog, rlat, rlon, area, rlonb, rlatb, rcap, & one_by_dy, one_by_dx, f_d, acosp, ks, lon, lat, lonb, latb, & fv_domain, print_freq, do_fms_tracer_manager, nqrst, mountain, & adjust_dry_mass, use_set_eta, fms_tracers_file, restart_format, & cold_start, gid, acap, acosu, ptop_min, pi, n_spong, n_diffu, & a_div, b_div, cosl5, sinl5, dx, n_split, iord_mt, iord_tm, & iord_tr, jord_mt, jord_tm, jord_tr, kord_mt, kord_tm, kord_tr, & fill, n_zonal, cp_vir, icase, ncnst, pnats, consv_te, & change_time, grid_weight #ifdef USE_LIMA public :: sc, se, dc, de, ipft1, ipft2, ifax, trigs, too_big #else public :: dry_mass, tra_step, dyn_step, p_ref #endif #ifdef MARS_GCM public :: fv_sponge_lev, fv_sponge_damp #endif !Balaji: needed by solo driver public :: map_dt, adiabatic !PUBLIC ROUTINES public fv_init, fv_end, d2a3d, p_var, polavg, drymadj, tracer_mass ! !REVISION HISTORY: ! ! 03.10.02 SJ Lin Modifications from fvGCM ! 03.11.12 SJ Lin Inline documentation ! 03.11.21 SJ Lin Change pt from potential temperature to temperature ! 04.11.29 SJ Lin Pre-2D decomposition modifications !----------------------------------------------------------------------- contains subroutine fv_init(ndt) use fms_io_mod, only : get_instance_filename #ifdef USE_LIMA use pft_module, only: pft_init, fftfax #endif !#include "fv_arrays.h" ! SJL note: ! This routine initializes various parameters for FV and touch the arrays ! that should promote better memory placement on share memory machines (e.g., SGI) integer, intent(in) :: ndt ! model time step (seconds) ! LOCAL VARIABLES: #ifdef USE_LIMA real rat, ycrit #endif real:: dim0 = 180. ! base dimension real:: dt0 = 1800. ! base time step real:: ns0 = 4. ! base nsplit for base dimension real dimx integer :: i, j, k, m, n integer :: unit, ierr, io integer :: kfirst, klast ! useless (but required) input to y_decomp at this point integer, allocatable:: ysize(:) ! tracers integer :: num_family, index ! output of register_tracers integer, allocatable :: tracer_indices(:) logical :: is_in_nonadv_tracer_section character(len=128) :: msg, tname, restart_file_name integer :: log_unit, out_unit map_dt = ndt ! mapping time step default to model time step call get_instance_filename('INPUT/fv_rst.res',restart_file_name ) cold_start = (.not.file_exist('INPUT/fv_rst.res') & & .and. .not.file_exist('INPUT/fv_rst.res.nc')) & & .and. .not.file_exist(restart_file_name) & & .and. .not.file_exist('INPUT/'//fms_tracers_file) !----- read namelist ----- if (file_exist('input.nml')) then unit = open_namelist_file ( ) ierr=1; do while (ierr /= 0) ! Every PE reading the same file? Ok on the SGI read (unit, nml=fv_core_nml, iostat=io, end=10) ierr = check_nml_error (io, 'fv_core_nml') enddo 10 call close_file (unit) endif ! ! override number of tracers by reading field_table if (file_exist('field_table')) do_fms_tracer_manager = .TRUE. if (do_fms_tracer_manager) then !{ #ifndef MARS_GCM #ifndef USE_LIMA if(full_phys .and. nt_phys /= 4) then !{ call error_mesg('FV_INIT', & & 'Wrong number of physics tracers nt_phys, should be 4! ', & FATAL) endif !} #else if(nt_phys /= 4) then !{ call error_mesg('FV_INIT', & & 'Wrong number of physics tracers nt_phys, should be 4! ', & FATAL) endif !} #endif #endif call tm_register_tracers (MODEL_ATMOS, ncnst, nt_prog, pnats, num_family) if(mpp_pe() == mpp_root_pe()) then out_unit = stdout() write(out_unit, *)'ncnst=', ncnst,' num_prog=',nt_prog,' pnats=',pnats,' num_family=',num_family endif #ifndef MARS_GCM #ifndef USE_LIMA if(full_phys .and. nt_prog < nt_phys) then !{ call error_mesg('FV_INIT', & & 'Number of prognostic tracers in field_table should be >= number of physics tracers', & FATAL) endif !} #else if(nt_prog < nt_phys) then !{ call error_mesg('FV_INIT', & & 'Number of prognostic tracers in field_table should be >= number of physics tracers', & FATAL) endif !} #endif #endif allocate(tracer_indices(ncnst)) call tm_get_tracer_indices(MODEL_ATMOS, tracer_indices) #ifndef USE_LIMA if(full_phys) then !{ #endif #ifndef MARS_GCM ! Check tracer positions in field_table ! ! sphum should the first entry ! followed by liq_wat, ice_wat, cld_amt (in any order) ! followed by all other advected tracers (if any) ! followed by non-advected tracers (if any) index = tm_get_tracer_index(MODEL_ATMOS, 'sphum', tracer_indices) if(index /= 1) then !{ call error_mesg('FV_INIT', & & 'sphum should be the first tracer entry in field table ', & FATAL) endif !} index = tm_get_tracer_index(MODEL_ATMOS, 'liq_wat', tracer_indices) if(index < 2 .or. index > 4) then !{ call error_mesg('FV_INIT', & & 'liq_wat should be tracer entry no 2 to 4 in field table ', & FATAL) endif !} index = tm_get_tracer_index(MODEL_ATMOS, 'ice_wat', tracer_indices) if(index < 2 .or. index > 4) then !{ call error_mesg('FV_INIT', & & 'ice_wat should be tracer entry no 2 to 4 in field table ', & FATAL) endif !} index = tm_get_tracer_index(MODEL_ATMOS, 'cld_amt', tracer_indices) if(index < 2 .or. index > 4) then !{ call error_mesg('FV_INIT', & & 'cld_amt should be tracer entry no 2 to 4 in field table ', & FATAL) endif #endif #ifndef USE_LIMA endif #endif if(age_tracer) then !{ ! age of air tracer should be present in field_table index = tm_get_tracer_index(MODEL_ATMOS, 'aoa', tracer_indices) if(index < 1 .or. index > ncnst) then !{ call error_mesg('FV_INIT', & & 'age_tracer=.TRUE.: entry "aoa" is required in field table but could not find it', & FATAL) endif !} endif !} ! check that non-advected tracers are defined last in field table is_in_nonadv_tracer_section = .FALSE. do i = 1, ncnst call tm_get_tracer_names(MODEL_ATMOS, tracer_indices(i), tname) if( is_in_nonadv_tracer_section .AND. & tm_check_if_prognostic(MODEL_ATMOS, tracer_indices(i)) ) then call tm_get_tracer_names(MODEL_ATMOS, tracer_indices(i), tname) write(msg, '(a,a,a)') 'Advected tracer ', tname, & ' entry found after non-advected tracer entry in field table.' call error_mesg('FV_INIT', msg, FATAL) endif if( .NOT. tm_check_if_prognostic(MODEL_ATMOS, tracer_indices(i)) ) then is_in_nonadv_tracer_section = .TRUE. endif enddo deallocate(tracer_indices) endif !} call write_version_number ( version, tagname ) if ( mpp_pe() == mpp_root_pe() ) then log_unit = stdlog() write (log_unit, nml=fv_core_nml) endif #ifdef SPMD ng_d = max( 2, min(3, max(jord_mt, jord_tm, jord_tr)) ) ng_s = 3 ighost = 0 gid = -999 !the PEs not participating in mod_comm will show this value call fv_arrays_init( nlon, mlat, nlev, ncnst, ng_d, ng_s, ighost, & fv_domain ) master = (gid == 0) !set some variables from SPMD correctly... call mpp_get_compute_domain( fv_domain, beglon, endlon, beglat, endlat ) !gid.NE.mpp_pe when npex.NE.1! call y_decomp( mlat, nlev, beglat, endlat, kfirst, klast, gid ) #else beglat = 1; endlat = mlat beglon = 1; endlon = nlon ng_d = 0; ng_s = 0 ighost = 0 master = .true. if( n_zonal == 0 ) n_zonal = 1 #endif if(master) then write(6,*) 'SMP pseudo domain decomposition in E-W =', n_zonal write(6,*) 'Cold_start=', cold_start endif call timing_init() call timing_on('Total') ! call timing_on('FV_init') !---------------------------------- ! Determine nsplit for the dynamics !---------------------------------- if ( n_split == 0 ) then dimx = max ( nlon, 2*(mlat-1) ) n_split = nint ( ns0*abs(ndt)*dimx/(dt0*dim0) + 0.49 ) n_split = max ( 1, n_split ) if(master) write(6,*) 'n_split is set to ', n_split, ' for (im,jm,dt)=',nlon,mlat,ndt else if(master) write(6,*) 'Using n_split from the namelist: ', n_split endif if(master) write(6,*) 'Total number of sponge layers= ', n_spong if ( adiabatic ) full_phys = .false. ! Set up some constants: if ( full_phys ) then cp_vir = (cv_wv + rvgas) / cp_air - 1. else cp_vir = 0. endif #ifdef MARS_GCM cp_vir = 0.0 #endif !-------------------------------------------- ! Set up horizontal grid geometrical factors: !-------------------------------------------- call set_fv_geom #ifndef USE_LIMA call pft_init(nlon, mlat, beglat, endlat, ighost, cosp, cose, master) #else ipft1 = max(2, beglat-ighost) js2g0 = max(2, beglat) ipft2 = min(mlat-1, endlat+ighost) jn1g1 = min(mlat, endlat+1) allocate( sc(ipft1:ipft2), se(js2g0:jn1g1) ) allocate( dc(nlon,ipft1:ipft2), de(nlon,js2g0:jn1g1) ) allocate( trigs(3*nlon/2+1) ) call fftfax(nlon, ifax, trigs) !--------------------------------------------- ! Determine ycrit such that effective DX >= DY !--------------------------------------------- rat = float(nlon)/float(2*(mlat-1)) ycrit = acos( min(0.81, rat) ) * (180./pi) if ( master ) write(6,*) 'Critical latitude for pft=', ycrit write (6,*) 'nlon, mlat, js2g0, jn1g1 ', nlon, mlat, js2g0, jn1g1 write (6,*) ' size(sc,1), size(dc,1), size(dc,2) ', size(sc,1), size(dc,1), size(dc,2) call pft_init(nlon, mlat, js2g0, jn1g1, & sc, se, dc, de, & cosp, cose, ycrit, ipft1, ipft2) #endif ! ! Vertical reference Eulerian coordinate for remapping: ! ks = -1000 ! something large enough to cause trouble if not initialized ptop = too_big allocate ( ak(nlev+1) ) allocate ( bk(nlev+1) ) call initialize_arrays call fv_print_chksums( 'Exiting fv_init' ) end subroutine fv_init subroutine initialize_arrays #include "fv_arrays.h" ! real, parameter :: too_big=huge(1.) integer :: i, j, k, m #include "fv_point.inc" !--------------------------- ! Initialize prognostic vars: !--------------------------- ! Note: the use of omp parallel is for better memory placement ! (the "first touch" principle) !$omp parallel do private (i, j, k) do k=ksp,kep !1,nlev do j=beglat-ng_d,endlat+ng_s do i=1,nlon u(i,j,k) = too_big enddo enddo do j=beglat-ng_d,endlat+ng_d do i=1,nlon v(i,j,k) = too_big enddo enddo do j=beglat-ng_d,endlat+ng_d do i=1,nlon pt(i,j,k) = too_big enddo enddo do j=beglat,endlat do i=1,nlon delp(i,j,k) = too_big enddo enddo do j=beglat,endlat do i=1,nlon pkz(i,j,k) = too_big enddo enddo enddo !------------------ ! Initialize tracers !------------------ if(do_fms_tracer_manager) then !{ ! use tracer_manager profile spec do m = 1, ncnst !!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! should use 1D decomp here !!!!!!!!!!!!!!!!!!!!!!!!!!!! call tm_set_tracer_profile (MODEL_ATMOS, m, q(:,:,:,m)) enddo else !}{ ! default do m=1,ncnst !$omp parallel do private (i, j, k) ! do k=1,nlev do k=ksp,kep do j=beglat-ng_d, endlat+ng_d do i=1,nlon q(i,j,k,m) = too_big enddo enddo enddo enddo endif !} !$omp parallel do private (i, j, k) ! do k=1,nlev+1 do k=ksp,kep+1 !NOTE that kep+1 overlaps between PEs! do j=beglat,endlat do i=1,nlon pk(i,j,k) = too_big enddo enddo do i = 1,nlon pesouth(i,k) = too_big end do enddo call fv_array_sync() ! The follwoing var are mostly used by physics !$omp parallel do private (i, j, k) ! do j=beglat, endlat do j=jsp,jep do i=1,nlon ps(i,j) = too_big enddo do k=1,nlev do i=1,nlon omga(i,j,k) = too_big enddo enddo do k=1, nlev+1 do i=1,nlon pe (i,k,j) = too_big peln(i,k,j) = too_big enddo enddo ! required to solve a denorm problem for HS (pletzer) do i = 1,nlon ps_bp(i,j) = too_big phis (i,j) = too_big u_srf(i,j) = too_big v_srf(i,j) = too_big end do enddo call fv_array_sync() end subroutine initialize_arrays !----------------------------------------------------------------------- subroutine get_eta_level(km, p_s, pf, ph, pscale) integer, intent(in) :: km real, intent(in) :: p_s ! unit: pascal real, intent(out) :: pf(km) real, intent(out) :: ph(km+1) real, intent(in), optional :: pscale integer k real ak1 if(nlev /= km) & call error_mesg('get_eta_level:','dimensionally inconsistent', FATAL) ph(1) = ak (1) do k=2,nlev+1 ph(k) = ak (k) + bk(k)*p_s enddo if ( present(pscale) ) then do k=1,nlev+1 ph(k) = pscale*ph(k) enddo endif if( ak(1) > ptop_min ) then pf(1) = (ph(2) - ph(1)) / log(ph(2)/ph(1)) else ak1 = (kappa + 1.) / kappa pf(1) = (ph(2) - ph(1)) / ak1 ! if(master) write(*,*) 'Modified p_full(1)=',pf(1), 0.5*(ph(1)+ph(2)) endif do k=2,nlev pf(k) = (ph(k+1) - ph(k)) / log(ph(k+1)/ph(k)) enddo end subroutine get_eta_level !----------------------------------------------------------------------- subroutine set_fv_geom ! Sets the geometric factors used in FV core for the lat-lon grid ! This routine should only be called once per run integer i, j, imh real dp, dl real ph5, zamda, zam5 real dg2r, asum, gmean allocate ( lon(nlon) ) allocate ( lonb(nlon+1) ) allocate ( rlonb(nlon+1) ) allocate ( lat(mlat) ) allocate ( latb(mlat+1) ) allocate ( rlatb(mlat+1) ) allocate ( rlat(nlon,beglat:endlat) ) allocate ( rlon(nlon,beglat:endlat) ) allocate ( area(nlon,beglat:endlat) ) allocate ( cose(mlat) ) allocate ( sine(mlat) ) allocate ( cosp(mlat) ) allocate ( sinp(mlat) ) allocate (acosp(mlat) ) allocate (acosu(mlat) ) allocate ( grid_weight(mlat) ) allocate ( coslon(nlon) ) allocate ( sinlon(nlon) ) allocate ( cosl5(nlon) ) allocate ( sinl5(nlon) ) allocate ( dx(mlat) ) allocate ( one_by_dx(mlat) ) #ifndef USE_LIMA allocate ( f_d(nlon,beglat-ng_d:endlat+ng_d) ) #else allocate ( f_d(beglat-ng_d:endlat+ng_d) ) #endif pi = 4.d0 * atan(1.d0) dl = (pi+pi) / nlon dp = pi/(mlat-1) sine(1) = -1. do j=2,mlat ph5 = -0.5d0*pi + (dble(j-1)-0.5d0)*(pi/dble(mlat-1)) sine(j) = sin(ph5) enddo do j=1,mlat-1 grid_weight(j) = sine(j+1) - sine(j) enddo grid_weight(mlat) = 1. - sine(mlat) ! Define area-conservative cosine at "cell center" do j=1,mlat cosp(j) = grid_weight(j) / dp acosp(j) = dp / grid_weight(j) grid_weight(j) = grid_weight(j) / (2*nlon) enddo ! Define cosine at edges.. cose(1) = 0. cose(2) = 0.5 * cosp(2) do j=3,mlat-1 cose(j) = 0.5 * (cosp(j-1) + cosp(j)) enddo cose(mlat) = cose(2) do j=2,mlat-1 acosu(j) = 2. / (cose(j) + cose(j+1)) enddo do j=1,mlat-1 sinp(j) = 0.5 * (sine(j) + sine(j+1)) enddo sinp(mlat) = 0.5 * (sine(mlat) + 1.) acap = nlon*(1.+sine(2))/dp ! acap = nlon * cosp(1) rcap = 1.d0 / acap ! Coriolis parameter at cell center #ifndef USE_LIMA do j = max(1,beglat-ng_d), min(mlat,endlat+ng_d) do i=1,nlon f_d(i,j) = (omega+omega) * sinp(j) enddo enddo #else do j = max(1,beglat-ng_d), min(mlat,endlat+ng_d) f_d(j) = (omega+omega) * sinp(j) enddo #endif !------------- ! Bounding box !------------- do j = 2, mlat latb(j) = -90. + (real(j-1)-0.5) * 180. / real(mlat-1) enddo latb( 1 ) = -90. latb( mlat+1 ) = 90. do j=1, mlat lat(j) = 0.5 * ( latb(j) + latb(j+1) ) enddo one_by_dy = 1. / ( radius*dp ) do j=2, mlat-1 dx(j) = dl*radius*cosp(j) one_by_dx(j) = 1. / dx(j) enddo dg2r = pi/180. do j=1,mlat+1 rlatb(j) = dg2r * latb(j) enddo !------------------------------- ! Set longitude dependent param !------------------------------- do i=1,nlon+1 #ifdef OLD_WAY lonb(i) = (real(i-1)-0.5) * 360. / real(nlon) #else !-------------------------- ! FMS B grid way: !-------------------------- lonb(i) = real(i-1)*360./real(nlon) #endif enddo do i=1,nlon lon(i) = 0.5 * ( lonb(i) + lonb(i+1) ) enddo do i=1,nlon+1 rlonb(i) = dg2r * lonb(i) enddo imh = nlon/2 if(nlon /= 4*(imh/2)) then call error_mesg('FV_INIT','nlon must be dividable by 4', FATAL) endif do i=1,imh #ifdef OLD_WAY zam5 = ((i-1)-0.5d0) * dl zamda = (i-1)*dl #else zam5 = (i-1) * dl zamda = (0.5d0+(i-1))*dl #endif cosl5(i) = cos(zam5) cosl5(i+imh) = -cosl5(i) sinl5(i) = sin(zam5) sinl5(i+imh) = -sinl5(i) coslon(i) = cos(zamda) coslon(i+imh) = -coslon(i) sinlon(i) = sin(zamda) sinlon(i+imh) = -sinlon(i) enddo do j=beglat, endlat do i=1,nlon rlat(i,j) = dg2r * lat(j) rlon(i,j) = dg2r * lon(i) ! Cubed sphere ready area(i,j) = (radius*dl*cosp(j)) * (dp*radius) enddo enddo asum = gmean(nlon, mlat, beglat, endlat, area) if(master) then write(6,*) 'Mean cell width (km)=', 1.e-3*sqrt(asum) endif end subroutine set_fv_geom subroutine fv_end(days, seconds) integer, intent(in):: days integer, intent(in):: seconds deallocate ( lon ) deallocate ( lat ) deallocate ( rlat ) deallocate ( rlon ) deallocate ( area ) deallocate ( lonb ) deallocate ( latb ) deallocate ( rlonb ) deallocate ( rlatb ) deallocate ( cose ) deallocate ( sine ) deallocate ( cosp ) deallocate ( sinp ) deallocate ( grid_weight ) deallocate ( acosp ) deallocate ( acosu ) deallocate ( coslon ) deallocate ( sinlon ) deallocate ( cosl5 ) deallocate ( sinl5 ) deallocate ( ak ) deallocate ( bk ) deallocate ( dx ) deallocate ( one_by_dx ) call fv_arrays_exit ! Clean up polar filter module: #ifndef USE_LIMA call pft_end #endif call timing_off('Total') #ifdef SPMD call timing_prt(gid) call mp_exit( .false. ) !Call mp_exit to cleanly exit mod_comm #else call timing_prt(0) #endif end subroutine fv_end subroutine d2a3d(u, v, ua, va, im, jm, km, & jfirst, jlast, ng_d, ng_s, coslon, sinlon) ! !ROUTINE: d2a3d -- Second order D-to-A grid transformation (3D) ! !INPUT PARAMETERS: integer, intent(in):: im, jm, km ! Dimensions integer, intent(in):: jfirst, jlast ! Latitude strip integer, intent(in):: ng_d, ng_s ! as defined in fvgcm.F real, intent(in):: coslon(im),sinlon(im) ! Sine and cosine in longitude real, intent(in):: v(im,jfirst-ng_d:jlast+ng_d,km) real, intent(inout):: u(im,jfirst-ng_d:jlast+ng_s,km) ! !OUTPUT PARAMETERS: real, intent(out):: ua(im,jfirst:jlast,km) ! U-Wind on A Grid real, intent(out):: va(im,jfirst:jlast,km) ! V-Wind on A Grid ! local integer i, j, k integer imh, js2g0 real un, vn, us, vs integer :: ks, ke imh = im/2 js2g0 = max(jfirst, 2) call fv_array_check( LOC(u ) ) call fv_array_check( LOC(v ) ) call fv_array_check( LOC(ua) ) call fv_array_check( LOC(va) ) #ifdef SPMD if ( .not. u_ghosted ) then call mp_send3d(gid-1, gid+1, im, jm, km, & 1, im, jfirst-ng_d, jlast+ng_s, 1, km, & 1, im, jfirst, jfirst, 1, km, u) endif #endif if( km.EQ.1 )then !this means processing the last level km=nlev, surface ks=1; ke=1 else if( km.EQ.nlev )then ks=ksp; ke=kep !shared arrays else call mpp_error( FATAL, 'D2A3D: shared pointers requires km.EQ.1 or km.EQ.nlev' ) end if !$omp parallel do private(i,j,k,us,vs,un,vn) do k=ks,ke !1,km do j=js2g0,jlast-1 do i=1,im ua(i,j,k) = 0.5*(u(i,j,k) + u(i,j+1,k)) enddo enddo do j=js2g0,min(jlast,jm-1) do i=1,im-1 va(i,j,k) = 0.5*(v(i,j,k) + v(i+1,j,k)) enddo va(im,j,k) = 0.5*(v(im,j,k) + v(1,j,k)) enddo if( jfirst == 1 ) then ! Projection at SP us = 0. vs = 0. do i=1,imh us = us + (ua(i+imh,2,k)-ua(i,2,k))*sinlon(i) & + (va(i,2,k)-va(i+imh,2,k))*coslon(i) vs = vs + (ua(i+imh,2,k)-ua(i,2,k))*coslon(i) & + (va(i+imh,2,k)-va(i,2,k))*sinlon(i) enddo us = us/im vs = vs/im do i=1,imh ua(i,1,k) = -us*sinlon(i) - vs*coslon(i) va(i,1,k) = us*coslon(i) - vs*sinlon(i) ua(i+imh,1,k) = -ua(i,1,k) va(i+imh,1,k) = -va(i,1,k) enddo endif if ( jlast == jm ) then ! Projection at NP un = 0. vn = 0. do i=1,imh un = un + (ua(i+imh,jm-1,k)-ua(i,jm-1,k))*sinlon(i) & + (va(i+imh,jm-1,k)-va(i,jm-1,k))*coslon(i) vn = vn + (ua(i,jm-1,k)-ua(i+imh,jm-1,k))*coslon(i) & + (va(i+imh,jm-1,k)-va(i,jm-1,k))*sinlon(i) enddo un = un/im vn = vn/im do i=1,imh ua(i,jm,k) = -un*sinlon(i) + vn*coslon(i) va(i,jm,k) = -un*coslon(i) - vn*sinlon(i) ua(i+imh,jm,k) = -ua(i,jm,k) va(i+imh,jm,k) = -va(i,jm,k) enddo endif enddo #ifdef SPMD if ( .not. u_ghosted ) then call mp_recv3d(gid+1, im, jm, km, & 1, im, jfirst-ng_d, jlast+ng_s, 1, km, & 1, im, jlast+1, jlast+1, 1, km, u) u_ghosted = .true. endif call fv_array_sync() if ( jlast < jm ) then !balaji, switched parallelization to i, since km might be 1 !$omp parallel do private(i, k) do k=1,km do i=isp,iep ua(i,jlast,k) = 0.5 * ( u(i,jlast,k) + u(i,jlast+1,k) ) enddo enddo endif #endif call fv_array_sync() end subroutine d2a3d subroutine p_var(im, jm, km, jfirst, jlast, ptop, & delp, ps, pe, peln, pk, pkz, & cappa, q, ng, nq, adjust_dry_mass ) ! Given (ptop, delp) computes (ps, pk, pe, peln, pkz) ! ! Input: integer, intent(in):: im, jm, km ! Dimensions integer, intent(in):: jfirst, jlast ! Latitude strip integer, intent(in):: nq integer, intent(in):: ng logical, intent(in):: adjust_dry_mass real, intent(in):: ptop real, intent(in):: cappa real, intent(inout):: delp(im,jfirst:jlast, km) real, intent(inout):: q(im,jfirst-ng:jlast+ng, km, nq) ! Output: real, intent(out) :: ps(im, jfirst:jlast) real, intent(out) :: pk(im, jfirst:jlast, km+1) real, intent(out) :: pe(im, km+1, jfirst:jlast) ! Edge pressure real, intent(out) :: peln(im, km+1, jfirst:jlast) ! Edge pressure real, intent(out) :: pkz(im, jfirst:jlast, km) ! Local real pek real lnp real ak1 integer i, j, k pek = ptop ** cappa ak1 = (cappa + 1.) / cappa call fv_array_check( LOC(delp) ) call fv_array_check( LOC(q) ) call fv_array_check( LOC(ps) ) call fv_array_check( LOC(pk) ) call fv_array_check( LOC(pe) ) call fv_array_check( LOC(peln) ) call fv_array_check( LOC(pkz) ) !$omp parallel do private(i,j,k) do j=jsp,jep !jfirst,jlast do i=1,im pe(i,1,j) = ptop pk(i,j,1) = pek enddo do k=2,km+1 do i=1,im pe(i,k,j) = pe(i,k-1,j) + delp(i,j,k-1) peln(i,k,j) = log(pe(i,k,j)) pk(i,j,k) = pe(i,k,j)**cappa enddo enddo do i=1,im ps(i,j) = pe(i,km+1,j) enddo !---- GFDL modification if( ptop < ptop_min ) then do i=1,im peln(i,1,j) = peln(i,2,j) - ak1 enddo else lnp = log( ptop ) do i=1,im peln(i,1,j) = lnp enddo endif !---- GFDL modification do k=1,km do i=1,im pkz(i,j,k) = (pk(i,j,k+1)-pk(i,j,k))/(cappa*(peln(i,k+1,j)-peln(i,k,j))) enddo enddo enddo call fv_array_sync() ! Check dry air mass: call drymadj( im, jm, km, jfirst, jlast, ng, mountain, cappa, & ptop, ps, delp, pe, pk, peln, pkz, nq, q, adjust_dry_mass ) end subroutine p_var subroutine polavg(p, im, jm, jfirst, jlast) integer, intent(in):: im, jm, jfirst, jlast real, intent(inout):: p(im,jfirst:jlast) ! local: integer i real sum1 if ( jfirst == 1 ) then sum1 = 0. do i=1,im sum1 = sum1 + p(i,1) enddo sum1 = sum1/im do i=1,im p(i,1) = sum1 enddo endif if ( jlast == jm ) then sum1 = 0. do i=1,im sum1 = sum1 + p(i,jm) enddo sum1 = sum1/im do i=1,im p(i,jm) = sum1 enddo endif end subroutine polavg subroutine drymadj( im, jm, km, jfirst, jlast, ng, & moun, cappa, ptop, ps, delp, pe, & pk, peln, pkz, nq, q, adjust_dry_mass ) ! !INPUT PARAMETERS: integer im, jm, km ! Dimensions integer jfirst, jlast ! Latitude strip integer nq ! Number of tracers integer ng logical moun logical, intent(in):: adjust_dry_mass real, intent(in):: ptop real, intent(in):: cappa ! !INPUT/OUTPUT PARAMETERS: real delp(im,jfirst:jlast,km) ! real pe(im,km+1,jfirst:jlast) ! real q(im,jfirst-ng:jlast+ng,km,nq) real ps(im,jfirst:jlast) ! surface pressure real, intent(inout):: pk(im, jfirst:jlast, km+1) real, intent(inout):: peln(im, km+1, jfirst:jlast) ! Edge pressure real, intent(inout):: pkz(im, jfirst:jlast, km) #ifdef USE_LIMA real drym ! parameter ( drym = 98222. ) ! setting for US NAVY 10 min data parameter ( drym = 98288. ) ! setting for USGS ! parameter ( drym = 98290. ) ! New setting for USGS #endif integer i, j, k real psmo real psdry, qtot real gmean real dpd integer ic, ip ! Local real psq(im,jfirst:jlast) real psd(im,jfirst:jlast) ! surface pressure due to dry air mass #ifdef use_shared_pointers pointer( p_psq, psq ) pointer( p_psd, psd ) call fv_stack_push( p_psq, im*(jlast-jfirst+1) ) call fv_stack_push( p_psd, im*(jlast-jfirst+1) ) #endif call fv_array_check( LOC(delp) ) call fv_array_check( LOC(pe) ) call fv_array_check( LOC(q) ) call fv_array_check( LOC(ps) ) call fv_array_check( LOC(pk) ) call fv_array_check( LOC(peln) ) call fv_array_check( LOC(pkz) ) ip = min(3,size(q,4)) !$omp parallel do & !$omp default(shared) & !$omp private(i,j,k) do 1000 j=jsp,jep !jfirst,jlast do i=1,im psd(i,j) = ptop psq(i,j) = 0. enddo #ifdef MARS_GCM do k=1,km do i=1,im psd(i,j) = psd(i,j) + delp(i,j,k) enddo enddo #else if( full_phys ) then do k=1,km do i=1,im psd(i,j) = psd(i,j) + delp(i,j,k)*(1.-q(i,j,k,1)) ! psq(i,j) = psq(i,j) + delp(i,j,k)*(q(i,j,k,1)+q(i,j,k,2)+q(i,j,k,3)) psq(i,j) = psq(i,j) + delp(i,j,k)*sum( q(i,j,k,1:ip) ) enddo enddo else do k=1,km do i=1,im psd(i,j) = psd(i,j) + delp(i,j,k) enddo enddo endif #endif MARS_GCM 1000 continue call fv_array_sync() !probably needed before gmean ! Check global maximum/minimum psmo = gmean(im, jm, jfirst, jlast, ps (1,jfirst)) psdry = gmean(im, jm, jfirst, jlast, psd(1,jfirst)) qtot = gmean(im, jm, jfirst, jlast, psq(1,jfirst)) if(master) then write(6,*) 'Total surface pressure (mb) = ', 0.01*psmo if ( full_phys ) then write(6,*) 'mean dry surface pressure = ', 0.01*psdry write(6,*) 'TPW-vapor (kg/m**2) =', (psmo-psdry)/GRAV write(6,*) 'TPW-total (kg/m**2) =', qtot/GRAV endif endif if( .not. adjust_dry_mass ) return #ifndef USE_LIMA if(moun) then dpd = dry_mass - psdry else dpd = 1000.*100. - psdry endif #else if(moun) then dpd = drym - psdry else dpd = 1000.*100. - psdry endif #endif if(master) write(6,*) 'dry mass to be added (pascals) =', dpd !$omp parallel do & !$omp default(shared) & !$omp private(i,j, ic) do 2000 j=jsp,jep !jfirst,jlast do ic=1,nq do i=1,im q(i,j,km,ic) = q(i,j,km,ic)*delp(i,j,km) / (delp(i,j,km)+dpd) enddo enddo ! Adjust the lowest Lagrangian layer do i=1,im delp(i,j,km) = delp(i,j,km) + dpd pe(i,km+1,j) = pe(i,km,j) + delp(i,j,km) ps(i,j) = pe(i,km+1,j) ! adjust pk, peln, pkz pk(i,j,km+1) = pe(i,km+1,j) ** cappa peln(i,km+1,j) = log(pe(i,km+1,j)) pkz(i,j,km) = (pk(i,j,km+1)-pk(i,j,km)) / & (cappa*(peln(i,km+1,j)-peln(i,km,j))) enddo 2000 continue call fv_array_sync() !probably needed before gmean psmo = gmean(im, jm, jfirst, jlast, ps(1,jfirst) ) if( master ) write(6,*) & 'Total (moist) surface pressure after adjustment = ', & 0.01*psmo end subroutine drymadj subroutine tracer_mass(im, jm, km, jfirst, jlast, ng, nq, q, delp) ! !INPUT PARAMETERS: integer im, jm, km ! Dimensions integer jfirst, jlast ! Latitude strip integer nq ! Number of tracers integer ng real delp(im,jfirst:jlast,km) ! real q(im,jfirst-ng:jlast+ng,km,nq) ! Local real gmean real qz(im,jfirst:jlast,nq) real qmass integer i, j, k integer iq #ifdef use_shared_pointers pointer( p_qz, qz ) call fv_stack_push( p_qz, im*(jlast-jfirst+1)*nq ) #endif call fv_array_check( LOC(delp) ) call fv_array_check( LOC(q) ) do iq=1,nq !$omp parallel do & !$omp default(shared) & !$omp private(i,j,k) do 1000 j=jsp,jep do i=1,im qz(i,j,iq) = 0. enddo do k=1,km do i=1,im qz(i,j,iq) = qz(i,j,iq) + delp(i,j,k)*q(i,j,k,iq) enddo enddo 1000 continue call fv_array_sync() qmass = gmean(im, jm, jfirst, jlast, qz(1,jfirst,iq)) if(master) write(6,*) iq, 'tracer mass = ', qmass enddo end subroutine tracer_mass end module fv_pack