! $Id: fv_pack.F90,v 19.0 2012/01/06 20:00:50 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, input_nml_file 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: tikal $' #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 ----- #ifdef INTERNAL_FILE_NML read (input_nml_file,fv_core_nml,iostat=io) ierr = check_nml_error (io, 'fv_core_nml') #else 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 #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