diff --git a/CMakeLists.txt b/CMakeLists.txt index afd224ed6..585cfb042 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -35,6 +35,8 @@ option(use_WRTCOMP "Enable compiler definition -Duse_WRTCOMP" OFF) option(INTERNAL_FILE_NML "Enable compiler definition -DINTERNAL_FILE_NML" ON) option(ENABLE_QUAD_PRECISION "Enable compiler definition -DENABLE_QUAD_PRECISION" ON) option(HYDRO "Enable compiler definition -DHYDRO" OFF) +option(ENABLE_RRFS_WAR "Enable independent parallel I/O when io_layout=1,1 -DENABLE_RRFS_WAR" OFF) +option(ENABLE_PARALLELRESTART "Enable collective parallel reads -DENABLE_PARALLELRESTART" OFF) find_package(MPI REQUIRED) if(OPENMP) @@ -112,7 +114,8 @@ list(APPEND driver_srcs driver/UFS/DYCORE_typedefs.F90 driver/UFS/fv_nggps_diag.F90 driver/UFS/fv_ufs_restart_io.F90 - driver/UFS/atmosphere.F90) + driver/UFS/atmosphere.F90 + driver/UFS/fv_ideal.F90) list(APPEND fv3_srcs ${model_srcs} ${tools_srcs}) @@ -167,6 +170,14 @@ if(ENABLE_QUAD_PRECISION) list(APPEND fv3_defs ENABLE_QUAD_PRECISION) endif() +if(ENABLE_PARALLELRESTART) + list(APPEND fv3_defs ENABLE_PARALLELRESTART) +endif() + +if(ENABLE_RRFS_WAR) + list(APPEND fv3_defs ENABLE_RRFS_WAR) +endif() + if(OPENMP) list(APPEND fv3_defs OPENMP) endif() diff --git a/docs/ideal_doc/README.md b/docs/ideal_doc/README.md new file mode 100644 index 000000000..2b39bc818 --- /dev/null +++ b/docs/ideal_doc/README.md @@ -0,0 +1,80 @@ +Some notes on ideal setup: + +The idealized option runs FV3 in a doubly-periodic domain. The environmental setup is done through the test_case_nml namelist. There are many hard-coded test cases in test_case.F90 that are only intended to work in FV3 SOLO, but here we provide the ability for the user to choose their idealized set up at run time. Options include several different hard-coded thermal and wind profile soundings along with the ability to read in a provided sounding file. Users can also add initial thermal bubbles at their chosen locations and intensities. + +Use of this capability require sfc_data.nc and oro_data.nc even if surface/PBL physics are not utilized. The grids in these files must be at least 1 grid point larger, in both x & y directions, than npx & npy. A configurable python script (create_ideal_sfc_oro_input.py) is provided here to generate idealized versions of these files. + +See driver/UFS/fv_ideal.F90 in the atmos_cubed_sphere directory for code specifics. + +Example namelist settings: + + &test_case_nml + bubble_type = 3 + ! 0 (default) : no bubble ; 1 : single bubble in domain center ; 2 : Line of N-S bubbles on left side of domain to initialize a squall line ; 3 : User-entered x,y locations of n_nun bubble centers + n_bub = 3 + ! number of bubbles (valid only if bubble_type = 3) + jcenters = 60, 60, 60 + ! j index of bubble centers; must be n_bub in length (valid only if bubble_type = 3) + icenters = 25, 50, 75 + ! i index of bubble centers; must be n_bub in length (valid only if bubble_type = 3) + t_profile = -1 + ! 0 prescribed Weisman-Klemp sounding (default) + ! 1 adiabatic sounding with adi_t K value + ! 2 isothermal sounding with iso_t K value + ! -1 read in sounding from external file input_sounding + q_profile = -1 + ! 0 prescribed Weisman-Klemp sounding (default) + ! 1 dry sounding + ! -1 read in sounding from external file input_sounding + ws_profile = -1 + ! 0 quarter circle hodograph (Harris Approximation) scaled by us0 (default) + ! 1 unidirectional Weisman-Klemp shear profile scaled by us0 + ! 2 constant u = us0, v = 0 + ! 3 constant v=us0, u=0 + ! 4 no wind + ! -1 read in sounding from external file input_sounding + bubble_t = 2.0 ! max center T (K) perturbation of bubbles + bubble_q = 0. ! max center qv (g/kg) perturbation + bubble_rad_x = 4000. ! bubble radius (m) in the x-direction (m; 10000 default) + bubble_rad_y = 4000. ! bubble radius (m) in the x-direction (m; 10000 default) + bubble_zc = 1500. ! bubble radius (m) in the z-direction (m; 1400 default) + do_coriolis = .false. ! set true for f-plane based using value of deglat + us0 = 12. ! scaling for wind profiles (m/s ; default 30) + adi_th = 300. ! adiabatic temperature (K) for t_profile=1 (default 300 K) + iso_th = 300. ! isothermal temperature (K) for t_profile=2 (default 300 K) + do_rand_perts = .false. ! if n_bub > 0, applies small amplitude(0.2k; 1E-7 kg/kg), random perturbations to T and qv values inside bubbles (default=.false.) + umove = 0.0 ! 'grid motion' u-component (m/s) subtracted from sounding U-wind to help keep storm centered in the domain + vmove = 0.0 ! 'grid motion' v-component (m/s) subtracted from sounding V-wind to help keep storm centered in the domain + / + +Other settings: + + &fv_core_nml + external_eta = .false. !required .false. + external_ic = .false. !required .false. + grid_type = 4 ! selects Cartesian periodic grid; required 4 + npz_type = 'meso' ! options here are buried in the code; ‘meso’ provides a lower top than !operational configuration (~20 hPa) + dx_const = 3000. ! set to desired value of dx (meters) + dy_const = 3000. ! set to desired value of dy (meters) + regional = .false. ! Required .false. + deglat = 15 ! grid latitude (default 15 degrees) Affects radiation and coriolis (f-plane) + deglon = 0 ! grid longitude (default 0) Only affects radiation (time of day relative to UTC) + / + +suite FV3_mp_nssl_ideal: +Simplest case is microphysics only, with this one set up to run the NSSL MP scheme. + +example input.nml: input.nml.idealbubble.test.3m (NSSL 3-moment mp) + +suite FV3_ideal_pbl_mp_nssl: + +This suite adds PBL/Surface physics. It needs input files + aerosol.dat + sfc_emissivity_idx.txt + solarconstant_noaa_an.txt + co2historicaldata_* + + Example input.nml: input.nml.idealpblbubble.test.3m (NSSL 3-moment mp) + + + diff --git a/docs/ideal_doc/create_ideal_sfc_oro_input.py b/docs/ideal_doc/create_ideal_sfc_oro_input.py new file mode 100755 index 000000000..6feda7ca7 --- /dev/null +++ b/docs/ideal_doc/create_ideal_sfc_oro_input.py @@ -0,0 +1,208 @@ +#!/usr/bin/env /home/Ted.Mansell/miniconda3/bin/python3 +# #!/usr/bin/env /home/Larissa.Reames/scratch/miniconda3/bin/python3 +# #!/usr/bin/env /Users/Shared/miniconda3/bin/python3 + +# Creates sfcdat_ideal.nc and orogdat_ideal.nc for UFS ideal setup +# THIS WILL NOT OVERWRITE EXISITNG FILES (Error code will be issued) +# nx and ny can be any value greater than or equal to the npy, npy in the fv_core_nml +# +# The values of lat and lon here don't matter (yet) for ideal cases. Use deglat +# and deglon in the fv_core_nml namelist to set these (probably only important for +# radiation and Coriolis, if enabled) +# +from netCDF4 import Dataset # Note: python is case-sensitive! +import numpy as np + +# filenames can be anything, but then must link as INPUT/sfc_data.nc and INPUT/oro_data.nc +sfcfilename = 'sfcdat_ideal.nc' +orofilename = 'orogdat_ideal.nc' + +# dimensions: +nx = 300 # should have nx >= npx +ny = 300 # should have ny >= npy +nz = 4 # number of soil layer; should be 4 for noahmp; Can set to 9 for ruc lsm and set namelist appropriately + +# The following are used to set various array values +tmpslmsk = 1 # sea/land mask (0=sea, 1=land) +tmpstype = 1 # land surface type; if set to 0 or 14, will get water surface +# soil types: 1, 'SAND' : 2, 'LOAMY SAND' : 3, 'SANDY LOAM' : 4, 'SILT LOAM' : 5, 'SILT' : 6, 'LOAM' : 7, 'SANDY CLAY LOAM' : 8, 'SILTY CLAY LOAM' : 9, 'CLAY LOAM' : 10,'SANDY CLAY' : 11,'SILTY CLAY' : 12,'CLAY' : 13,'ORGANIC MATERIAL' : 14,'WATER' +tmpvtype = 7 # vegetation type; see noahmptable.tbl; 7 = grassland +tmpvfrac = 0.9 # vegetation fraction +tmpsmc = 0.2 # soil moisture content; set to 1 for water surface; sets all levels to the same value +tmpslc = 0.2 # soil liquid content; set to 0 for water surface; sets all levels to the same value +tmpstc = 300 # soil temperature (all levels); sets all levels to the same value +lat1 = 35 # value is ignored by FV3, which will use deglat from input namelist +lon1 = 270 # value is ignored by FV3, which will use deglon from input namelist +tem2m = 300 # t2m = 2-meter temperature +q2mtmp = 0.009 # 2-meter specific humidity +ust = 0.2 # uustar : for surface layer scheme; cannot be zero + +# Note that the variable names (key) are always the same as the long names, so longName is not actually used +sfc_vars = { 'slmsk': {'longName':'slmsk','units':'none','ndims':'2','value':tmpslmsk}, # sea/land mask array (sea:0,land:1,sea-ice:2) + 'tsea': {'longName':'tsea','units':'none','ndims':'2','value':tmpstc}, # tsea is read into variable 'tsfc' : surface air temperature + 'sheleg': {'longName':'sheleg','units':'none','ndims':'2','value':0}, # read into variable 'weasd' : water equiv of accumulated snow depth (kg/m**2) + 'tg3': {'longName':'tg3','units':'none','ndims':'2','value':300}, # deep soil temperature + 'zorl': {'longName':'zorl','units':'none','ndims':'2','value':0.0001}, # composite surface roughness in cm + 'alvsf': {'longName':'alvsf','units':'none','ndims':'2','value':0.06}, # mean vis albedo with strong cosz dependency + 'alvwf': {'longName':'alvwf','units':'none','ndims':'2','value':0.06}, # mean vis albedo with weak cosz dependency + 'alnsf': {'longName':'alnsf','units':'none','ndims':'2','value':0.06}, # mean nir albedo with strong cosz dependency + 'alnwf': {'longName':'alnwf','units':'none','ndims':'2','value':0.06}, # mean nir albedo with weak cosz dependency + 'facsf': {'longName':'facsf','units':'none','ndims':'2','value':0}, # fractional coverage with strong cosz dependency + 'facwf': {'longName':'facwf','units':'none','ndims':'2','value':0}, # fractional coverage with weak cosz dependency + 'vfrac': {'longName':'vfrac','units':'none','ndims':'2','value':tmpvfrac}, # vegetation fraction + 'canopy': {'longName':'canopy','units':'none','ndims':'2','value':0}, # canopy water + 'f10m': {'longName':'f10m','units':'none','ndims':'2','value':0.997}, # fm at 10 m : Ratio of sigma level 1 wind and 10m wind + 't2m': {'longName':'t2m','units':'none','ndims':'2','value':tem2m}, # 2 meter temperature + 'q2m': {'longName':'q2m','units':'none','ndims':'2','value':q2mtmp}, # 2 meter humidity + 'vtype': {'longName':'vtype','units':'none','ndims':'2','value':tmpvtype}, # vegetation type + 'stype': {'longName':'stype','units':'none','ndims':'2','value':tmpstype}, # soil type + 'uustar': {'longName':'uustar','units':'none','ndims':'2','value':ust}, # boundary layer parameter (must be nonzero positive) + 'ffmm': {'longName':'ffmm','units':'none','ndims':'2','value':10}, # fm parameter from PBL scheme + 'ffhh': {'longName':'ffhh','units':'none','ndims':'2','value':0}, # fh parameter from PBL scheme + 'hice': {'longName':'hice','units':'none','ndims':'2','value':0}, # sea ice thickness + 'fice': {'longName':'fice','units':'none','ndims':'2','value':0}, # ice fraction over open water grid + 'tisfc': {'longName':'tisfc','units':'none','ndims':'2','value':290},# surface temperature over ice fraction + 'tprcp': {'longName':'tprcp','units':'none','ndims':'2','value':0}, # sfc_fld%tprcp - total precipitation + 'srflag': {'longName':'srflag','units':'none','ndims':'2','value':0},# sfc_fld%srflag - snow/rain flag for precipitation + 'snwdph': {'longName':'snwdph','units':'none','ndims':'2','value':0}, # snow depth water equivalent in mm ; same as snwdph + 'shdmin': {'longName':'shdmin','units':'none','ndims':'2','value':0}, # min fractional coverage of green veg + 'shdmax': {'longName':'shdmax','units':'none','ndims':'2','value':0}, # max fractnl cover of green veg (not used) + 'slope': {'longName':'slope','units':'none','ndims':'2','value':0}, # sfc slope type for lsm + 'snoalb': {'longName':'snoalb','units':'none','ndims':'2','value':0}, # maximum snow albedo in fraction + 'stc': {'longName':'stc','units':'none','ndims':'3','value':tmpstc}, # soil temperature (noah) + 'smc': {'longName':'stc','units':'none','ndims':'3','value':tmpsmc}, # total soil moisture content (noah) + 'slc': {'longName':'stc','units':'none','ndims':'3','value':tmpslc}, # liquid soil moisture content (noah) + 'tref': {'longName':'tref','units':'none','ndims':'2','value':290}, # stuff for near surface sea temperature model (NSSTM) + 'z_c': {'longName':'z_c','units':'none','ndims':'2','value':0.001}, # NSSTM + 'c_0': {'longName':'c_0','units':'none','ndims':'2','value':0.05}, # NSSTM + 'c_d': {'longName':'c_d','units':'none','ndims':'2','value':-50}, # NSSTM + 'w_0': {'longName':'w_0','units':'none','ndims':'2','value':0}, # NSSTM + 'w_d': {'longName':'w_d','units':'none','ndims':'2','value':0}, # NSSTM + 'xt': {'longName':'xt','units':'none','ndims':'2','value':0}, # NSSTM + 'xs': {'longName':'xs','units':'none','ndims':'2','value':0}, # NSSTM + 'xu': {'longName':'xu','units':'none','ndims':'2','value':0}, # NSSTM + 'xv': {'longName':'xv','units':'none','ndims':'2','value':0}, # NSSTM + 'xz': {'longName':'xz','units':'none','ndims':'2','value':30}, # NSSTM + 'zm': {'longName':'zm','units':'none','ndims':'2','value':0}, # NSSTM + 'xtts': {'longName':'xtts','units':'none','ndims':'2','value':0}, # NSSTM + 'xzts': {'longName':'xzts','units':'none','ndims':'2','value':0}, # NSSTM + 'd_conv': {'longName':'d_conv','units':'none','ndims':'2','value':0}, # NSSTM + 'ifd': {'longName':'ifd','units':'none','ndims':'2','value':0}, # NSSTM + 'dt_cool':{'longName':'dt_cool','units':'none','ndims':'2','value':0.3}, # NSSTM + 'qrain': {'longName':'qrain','units':'none','ndims':'2','value':0} # NSSTM + } + +oro_vars = { # keys +'slmsk': {'longName':'slmsk','value':tmpslmsk}, # sea/land mask array (sea:0,land:1,sea-ice:2) +'land_frac': {'longName':'land_frac','value':1}, # land fraction [0:1] +'orog_raw': {'longName':'orog_raw','value':0}, # oro_uf (:) => null() !< unfiltered orography +'orog_filt': {'longName':'orog_filt','value':0}, # oro (:) => null() !< orography +'stddev': {'longName':'stddev','value':0}, # hprime (:,1) ; orographic metrics +'convexity': {'longName':'convexity','value':0}, # hprime (:,2) ; orographic metrics +'oa1': {'longName':'oa1','value':0}, # hprime (:,3) ; orographic metrics +'oa2': {'longName':'oa2','value':0}, # hprime (:,4) ; orographic metrics +'oa3': {'longName':'oa3','value':0}, # hprime (:,5) ; orographic metrics +'oa4': {'longName':'oa4','value':0}, # hprime (:,6) ; orographic metrics +'ol1': {'longName':'ol1','value':0}, # hprime (:,7) ; orographic metrics +'ol2': {'longName':'ol2','value':0}, # hprime (:,8) ; orographic metrics +'ol3': {'longName':'ol3','value':0}, # hprime (:,9) ; orographic metrics +'ol4': {'longName':'ol4','value':0}, # hprime (:,10) ; orographic metrics +'theta': {'longName':'theta','value':0}, # hprime (:,11) ; orographic metrics +'gamma': {'longName':'gamma','value':0}, # hprime (:,12) ; orographic metrics +'sigma': {'longName':'sigma','value':0}, # hprime (:,13) ; orographic metrics +'elvmax': {'longName':'elvmax','value':0} # hprime (:,14) ; orographic metrics +} + +ncfile = Dataset(sfcfilename,mode='w',clobber=False,format='NETCDF4_CLASSIC') + + +x_dim = ncfile.createDimension('xaxis_1', nx) # latitude axis +y_dim = ncfile.createDimension('yaxis_1', ny) # longitude axis +z_dim = ncfile.createDimension('zaxis_1', nz) # longitude axis +time_dim = ncfile.createDimension('Time', 1) # unlimited axis (can be appended to). + +xdim = ncfile.createVariable('xaxis_1', np.float32, ('xaxis_1',)) +xdim.units = 'none' +xdim.longName = 'xaxis_1' +xdim.cartesian_axis = 'X' +xdim[:] = 1 + np.arange(nx) + +ydim = ncfile.createVariable('yaxis_1', np.float32, ('yaxis_1',)) +ydim.units = 'none' +ydim.longName = 'yaxis_1' +ydim.cartesian_axis = 'Y' +ydim[:] = 1 + np.arange(ny) + +zdim = ncfile.createVariable('zaxis_1', np.float32, ('zaxis_1',)) +zdim.units = 'none' +zdim.longName = 'zaxis_1' +zdim.cartesian_axis = 'Z' +zdim[:] = 1 + np.arange(nz) + +timedim = ncfile.createVariable('Time', np.float32, ('Time',)) +timedim.units = 'Time level' +timedim.longName = 'Time' +timedim.cartesian_axis = 'T' +timedim[:] = 1 + +geolat = ncfile.createVariable('geolat', np.float64, ('yaxis_1','xaxis_1'),zlib=True) +geolat.units = 'degrees_north' +geolat.longName = 'Latitude' +geolat[:] = lat1 + +geolon = ncfile.createVariable('geolon', np.float64, ('yaxis_1','xaxis_1'),zlib=True) +geolon.units = 'degrees_east' +geolon.longName = 'Longitude' +geolon[:] = lon1 + +for k, key in enumerate(sfc_vars): + if sfc_vars[key]["ndims"] == '2': + var = ncfile.createVariable(key, np.float64, ('yaxis_1','xaxis_1'),zlib=True) + else: + var = ncfile.createVariable(key, np.float64, ('zaxis_1','yaxis_1','xaxis_1'),zlib=True) + #print('key = ',key) + print('sfcvar = ',sfc_vars[key]) + #print('sfcvar = ',sfc_vars[key]["longName"]) + var.longName = key + var.units = sfc_vars[key]["units"] + var.coordinates = "geolon geolat" + var[:] = sfc_vars[key].get("value") + + +# first print the Dataset object to see what we've got +print(ncfile) +# close the Dataset. +ncfile.close(); print('Sfc Dataset is closed!') + + +# Ideal oro_data.nc +ncfile = Dataset(orofilename,mode='w',clobber=False,format='NETCDF4_CLASSIC') + +lat = ncfile.createDimension('lat', nx) # latitude axis +lon = ncfile.createDimension('lon', ny) # longitude axis + +geolat = ncfile.createVariable('geolat', np.float32, ('lat','lon'),zlib=True) +geolat.units = 'degrees_north' +geolat.long_name = 'Latitude' +geolat[:] = lat1 + +geolon = ncfile.createVariable('geolon', np.float32, ('lat','lon'),zlib=True) +geolon.units = 'degrees_east' +geolon.long_name = 'Longitude' +geolon[:] = lon1 + + +for k, key in enumerate(oro_vars): + var = ncfile.createVariable(key, np.float32, ('lat','lon'),zlib=True) + #print('key = ',key) + print('orovar = ',oro_vars[key]) + #print('sfcvar = ',sfc_vars[key]["longName"]) + var.coordinates = "geolon geolat" + var[:] = oro_vars[key].get("value") + + +# first print the Dataset object to see what we've got +print(ncfile) +# close the Dataset. +ncfile.close(); print('Orog Dataset is closed!') + diff --git a/driver/UFS/atmosphere.F90 b/driver/UFS/atmosphere.F90 index b9abdb28d..b6a37c373 100644 --- a/driver/UFS/atmosphere.F90 +++ b/driver/UFS/atmosphere.F90 @@ -968,13 +968,14 @@ end subroutine get_nth_domain_info !! decomposition for the current cubed-sphere tile. !>@detail Coupling is done using the mass/temperature grid with no halos. subroutine atmosphere_domain ( fv_domain, rd_domain, layout, regional, nested, & - ngrids_atmos, mygrid_atmos, pelist ) + ngrids_atmos, mygrid_atmos, pelist, grid_type) type(domain2d), intent(out) :: fv_domain, rd_domain integer, intent(out) :: layout(2) logical, intent(out) :: regional logical, intent(out) :: nested integer, intent(out) :: ngrids_atmos integer, intent(out) :: mygrid_atmos + integer, intent(out) :: grid_type integer, pointer, intent(out) :: pelist(:) integer :: n @@ -988,6 +989,7 @@ subroutine atmosphere_domain ( fv_domain, rd_domain, layout, regional, nested, & mygrid_atmos = mygrid call set_atmosphere_pelist() pelist => Atm(mygrid)%pelist + grid_type = Atm(mygrid)%flagstruct%grid_type end subroutine atmosphere_domain diff --git a/driver/UFS/fv_ideal.F90 b/driver/UFS/fv_ideal.F90 new file mode 100644 index 000000000..d2d232a93 --- /dev/null +++ b/driver/UFS/fv_ideal.F90 @@ -0,0 +1,904 @@ +!*********************************************************************** +!* GNU Lesser General Public License +!* +!* This file is part of the FV3 dynamical core. +!* +!* The FV3 dynamical core is free software: you can redistribute it +!* and/or modify it under the terms of the +!* GNU Lesser General Public License as published by the +!* Free Software Foundation, either version 3 of the License, or +!* (at your option) any later version. +!* +!* The FV3 dynamical core is distributed in the hope that it will be +!* useful, but WITHOUT ANYWARRANTY; 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 Lesser General Public +!* License along with the FV3 dynamical core. +!* If not, see . +!*********************************************************************** + + module fv_ideal_mod + +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! > +! +! +! +! +!
Module NameFunctions Included
constants_modcnst_radius=>radius, pi=>pi_8, omega, grav, rdgas, +! cp_air, rvgas
init_hydro_modp_var
fv_mp_modis_master
fv_grid_utils_modptop_min
mpp_modmpp_error, FATAL, stdlog, input_nml_file
fms_modcheck_nml_error
mpp_domains_moddomain2d
fv_arrays_modfv_grid_type, fv_flags_type, fv_grid_bounds_type, R_GRID
+ +#ifdef OVERLOAD_R4 + use constantsR4_mod, only: pi=>pi_8, omega, grav, kappa, rdgas, cp_air, rvgas +#else + use constants_mod, only: pi=>pi_8, omega, grav, kappa, rdgas, cp_air, rvgas +#endif + use init_hydro_mod, only: p_var + use fv_mp_mod, only: is_master + use fv_grid_utils_mod, only: ptop_min + use mpp_mod, only: mpp_error, FATAL, stdlog, input_nml_file + use fms_mod, only: check_nml_error + use mpp_domains_mod, only: domain2d + use fv_arrays_mod, only: fv_grid_type, fv_flags_type, fv_grid_bounds_type, R_GRID + implicit none + private + + integer :: sphum, theta_d + integer, parameter :: max_bub=100 + real(kind=R_GRID), parameter :: one = 1.d0 + integer :: test_case = 11 + logical :: bubble_do = .false. + logical :: do_rand_perts = .false. + real :: alpha = 0.0 + integer :: Nsolitons = 1, n_bub = 1 + real :: soliton_size = 750.e3, soliton_Umax = 50. + integer :: t_profile = 0, q_profile = 0, ws_profile = 0 + integer :: do_coriolis = 0, bubble_type = 0 + real :: bubble_t = 2., bubble_q = 0., bubble_rad_x = 10.0E3 + real :: umove = 0.0, vmove = 0.0 + real :: p00_in = 1.e5 + real :: bubble_rad_y = 10.0E3, bubble_zc = 1.4E3 + real :: iso_t = 300., adi_th = 300., us0 = 30. + real,dimension(max_bub) :: icenters, jcenters + + public :: fv_init_ideal + public :: read_namelist_fv_ideal + contains + + subroutine fv_init_ideal(u,v,w,pt,delp,q,phis, ps,pe,peln,pk,pkz, uc,vc, ua,va, ak, bk, & + gridstruct, flagstruct, npx, npy, npz, ng, ncnst, nwat, ndims, nregions, dry_mass, & + mountain, moist_phys, hydrostatic, hybrid_z, delz, ze0, ks, ptop, domain_in, tile_in, bd) + + type(fv_grid_bounds_type), intent(IN) :: bd + real , intent(INOUT) :: u(bd%isd:bd%ied ,bd%jsd:bd%jed+1,npz) + real , intent(INOUT) :: v(bd%isd:bd%ied+1,bd%jsd:bd%jed ,npz) + real , intent(INOUT) :: w(bd%isd: ,bd%jsd: ,1:) + real , intent(INOUT) :: pt(bd%isd:bd%ied ,bd%jsd:bd%jed ,npz) + real , intent(INOUT) :: delp(bd%isd:bd%ied ,bd%jsd:bd%jed ,npz) + real , intent(INOUT) :: q(bd%isd:bd%ied ,bd%jsd:bd%jed ,npz, ncnst) + + real , intent(INOUT) :: phis(bd%isd:bd%ied ,bd%jsd:bd%jed ) + + real , intent(INOUT) :: ps(bd%isd:bd%ied ,bd%jsd:bd%jed ) + real , intent(INOUT) :: pe(bd%is-1:bd%ie+1,npz+1,bd%js-1:bd%je+1) + real , intent(INOUT) :: pk(bd%is:bd%ie ,bd%js:bd%je ,npz+1) + real , intent(INOUT) :: peln(bd%is :bd%ie ,npz+1 ,bd%js:bd%je) + real , intent(INOUT) :: pkz(bd%is:bd%ie ,bd%js:bd%je ,npz ) + real , intent(INOUT) :: uc(bd%isd:bd%ied+1,bd%jsd:bd%jed ,npz) + real , intent(INOUT) :: vc(bd%isd:bd%ied ,bd%jsd:bd%jed+1,npz) + real , intent(INOUT) :: ua(bd%isd:bd%ied ,bd%jsd:bd%jed ,npz) + real , intent(INOUT) :: va(bd%isd:bd%ied ,bd%jsd:bd%jed ,npz) + real , intent(inout) :: delz(bd%is:,bd%js:,1:) + real , intent(inout) :: ze0(bd%is:,bd%js:,1:) + + real , intent(inout) :: ak(npz+1) + real , intent(inout) :: bk(npz+1) + + integer, intent(IN) :: npx, npy, npz + integer, intent(IN) :: ng, ncnst, nwat + integer, intent(IN) :: ndims + integer, intent(IN) :: nregions + + real, intent(IN) :: dry_mass + logical, intent(IN) :: mountain + logical, intent(IN) :: moist_phys + logical, intent(IN) :: hydrostatic, hybrid_z + integer, intent(INOUT) :: ks + integer, intent(INOUT), target :: tile_in + real, intent(INOUT) :: ptop + + type(domain2d), intent(IN), target :: domain_in + + type(fv_grid_type), intent(IN), target :: gridstruct + type(fv_flags_type), intent(IN), target :: flagstruct + + integer, parameter :: nl_max = 2500 + real, dimension(nl_max) :: z_snd, p_snd, t_snd, rho_snd, u_snd, v_snd, qv_snd + real, dimension(bd%is:bd%ie):: pm, qs + real, dimension(1:npz):: pk1, pe1, ts1, qs1, dummy + !real :: us0 = 30. + real :: dist,r0, f0_const, prf, rgrav, xrad, yrad, zrad,RAD + real :: xradbub, yradbub + real :: ptmp, ze, zc, zm, utmp, vtmp + real :: t00, p00, xmax, xc, xx, yy, pk0, pturb, ztop + real :: ze1(npz+1) + real:: dz1(npz) + real:: zvir, rand1, rand2 + real:: amplitude = 0.2 + integer :: i, j, k, m, icenter, jcenter, nl_snd + + real, pointer, dimension(:,:,:) :: agrid, grid + real(kind=R_GRID), pointer, dimension(:,:) :: area + real, pointer, dimension(:,:) :: rarea, fC, f0 + real, pointer, dimension(:,:,:) :: ee1, ee2, en1, en2 + real, pointer, dimension(:,:,:,:) :: ew, es + real, pointer, dimension(:,:) :: dx,dy, dxa,dya, rdxa, rdya, dxc,dyc + + logical, pointer :: cubed_sphere, latlon + + type(domain2d), pointer :: domain + integer, pointer :: tile + + logical, pointer :: have_south_pole, have_north_pole + + integer, pointer :: ntiles_g + real, pointer :: acapN, acapS, globalarea + + real(kind=R_GRID), pointer :: dx_const, dy_const + + integer :: is, ie, js, je + integer :: isd, ied, jsd, jed + + integer :: b + + is = bd%is + ie = bd%ie + js = bd%js + je = bd%je + isd = bd%isd + ied = bd%ied + jsd = bd%jsd + jed = bd%jed + + agrid => gridstruct%agrid + grid => gridstruct%grid + + area => gridstruct%area_64 + + dx => gridstruct%dx + dy => gridstruct%dy + dxa => gridstruct%dxa + dya => gridstruct%dya + rdxa => gridstruct%rdxa + rdya => gridstruct%rdya + dxc => gridstruct%dxc + dyc => gridstruct%dyc + + fC => gridstruct%fC + f0 => gridstruct%f0 + + !These are frequently used and so have pointers set up for them + dx_const => flagstruct%dx_const + dy_const => flagstruct%dy_const + + domain => domain_in + tile => tile_in + + have_south_pole => gridstruct%have_south_pole + have_north_pole => gridstruct%have_north_pole + + ntiles_g => gridstruct%ntiles_g + acapN => gridstruct%acapN + acapS => gridstruct%acapS + globalarea => gridstruct%globalarea + + if (do_coriolis>=1) then + f0_const = 2.*omega*sin(flagstruct%deglat/180.*pi) + else + f0_const = 0.0 + endif + + f0(:,:) = f0_const + fC(:,:) = f0_const + + q = 0. + + zvir = rvgas/rdgas - 1. + p00 = p00_in ! 1000.E2 + print*, do_coriolis, t_profile, q_profile, ws_profile + if (t_profile == -1 .or. q_profile == -1 .or. ws_profile == -1) then + call get_sounding( z_snd, p_snd, t_snd, rho_snd, u_snd, v_snd, qv_snd, nl_max, nl_snd, p00) + IF ( is_master() ) write(*,*) 'Using p00 from sounding file: p00 = ',p00 + endif + + ! Set up pressure arrays + ps(:,:) = p00 + phis(:,:) = 0. + do j=js,je + do i=is,ie + pk(i,j,1) = ptop**kappa + pe(i,1,j) = ptop + peln(i,1,j) = log(ptop) + enddo + enddo + + do k=1,npz + do j=js,je + do i=is,ie + delp(i,j,k) = ak(k+1)-ak(k) + ps(i,j)*(bk(k+1)-bk(k)) + pe(i,k+1,j) = ak(k+1) + ps(i,j)*bk(k+1) + peln(i,k+1,j) = log(pe(i,k+1,j)) + pk(i,j,k+1) = exp( kappa*peln(i,k+1,j) ) + enddo + enddo + enddo + + i = is + j = js + do k=1,npz + pk1(k) = (pk(i,j,k+1)-pk(i,j,k))/(kappa*(peln(i,k+1,j)-peln(i,k,j))) + pe1(k) = (pe(i,k+1,j)-pe(i,k,j))/(peln(i,k+1,j)-peln(i,k,j)) + enddo + + !Set model thermodynamic profile based on namelist inputs + ! Read profile from sounding + if (t_profile == -1) then + do k = 1,npz + ts1(k) = interp_log( t_snd, p_snd, pe1(k), nl_max, nl_snd ) + enddo + elseif ( t_profile == 0 ) then + !Generate GFDL's supercell sounding + call SuperCell_Sounding(npz, p00, pk1, ts1, qs1) + elseif (t_profile == 1 ) then + !adiabatic + print*, "kappa, adi_th = ", kappa, adi_th + do k=1,npz + ts1(k) = adi_th * ( pe1(k) / 1E5) ** kappa + enddo + elseif (t_profile == 2) then + !isothermal + ts1(:) = iso_t + else + call mpp_error(FATAL, " t_profile ", t_profile ," not defined" ) + endif + + if (q_profile == -1) then + ! Read profile from sounding + do k = 1,npz + qs1(k) = interp_log( qv_snd, p_snd, pe1(k), nl_max,nl_snd ) + enddo + elseif ( q_profile==0 ) then + if (t_profile == 0) then + ! qs1 already computed prior, move along + else + ! Generate GFDL's supercell sounding + call SuperCell_Sounding(npz, p00, pk1, dummy, qs1) + endif + elseif (q_profile==1 ) then + ! dry environment + qs1(:) = 1E-9 + else + call mpp_error(FATAL, " q_profile ", q_profile ," not defined" ) + endif + + ! Compute delz from ts1 and qs1 + w(:,:,:) = 0. + q(:,:,:,:) = 0. + + do k=1,npz + do j=js,je + do i=is,ie + pt(i,j,k) = ts1(k) + q(i,j,k,1) = qs1(k) + delz(i,j,k)=rdgas/grav*ts1(k)*(1.+zvir*qs1(k))*(peln(i,k,j)-peln(i,k+1,j)) + enddo + enddo + enddo + + ze1(npz+1) = 0. + do k=npz,1,-1 + ze1(k) = ze1(k+1) - delz(is,js,k) + enddo + + ! Set up model winds + !Read winds from sounding + if (ws_profile == -1) then + do k = 1,npz + zm = 0.5*(ze1(k)+ze1(k+1)) + utmp = interp_lin( u_snd, z_snd, zm, nl_max, nl_snd ) + vtmp = interp_lin( v_snd, z_snd,zm, nl_max, nl_snd ) + do j=js,je+1 + do i=is,ie + u(i,j,k) = utmp + enddo + enddo + do j=js,je + do i=is,ie+1 + v(i,j,k) = vtmp + enddo + enddo + enddo + + elseif ( ws_profile==0 ) then + ! Quarter-circle hodograph (Harris approximation) + do k=1,npz + zm = 0.5*(ze1(k)+ze1(k+1)) + if ( zm .le. 2.e3 ) then + utmp = 8.*(1.-cos(pi*zm/4.e3)) + vtmp = 8.*sin(pi*zm/4.e3) + elseif (zm .le. 6.e3 ) then + utmp = 8. + (us0-8.)*(zm-2.e3)/4.e3 + vtmp = 8. + else + utmp = us0 + vtmp = 8. + endif +! u-wind + do j=js,je+1 + do i=is,ie + u(i,j,k) = utmp - 8. + enddo + enddo +! v-wind + do j=js,je + do i=is,ie+1 + v(i,j,k) = vtmp - 4. + enddo + enddo + enddo + elseif (ws_profile==1 ) then + ! Unidirectional WK shear + v(:,:,:) = 0. + do k=1,npz + zm = 0.5*(ze1(k)+ze1(k+1)) + if ( zm .le. 6.e3 ) then + u(:,:,k) = us0 * tanh(zm/3.e3) + else + u(:,:,k) = us0 + endif + enddo + elseif (ws_profile==2 ) then + ! constant u, v=0 + u(:,:,:) = us0 + v(:,:,:) = 0. + elseif (ws_profile==3 ) then + ! constant v, u=0 + u(:,:,:) = 0. + v(:,:,:) = us0 + elseif (ws_profile==4 ) then + u(:,:,:) = 0. + v(:,:,:) = 0. + elseif (ws_profile==5) then + ! Linear WK shear below 2.5km + v(:,:,:) = 0. + do k=1,npz + zm = 0.5*(ze1(k)+ze1(k+1)) + if ( zm .le. 2.5e3 ) then + u(:,:,k) = us0 * tanh(zm/3.e3) + else + u(:,:,k) = us0 + endif + enddo + else + call mpp_error(FATAL, " ws_profile ", ws_profile ," not defined" ) + endif + + IF ( is_master() ) THEN + write(*,*) 'Final sounding: k, z, t, q, u, v, p' + do k = 1,npz + write(*,*) k,ze1(k),ts1(k),qs1(k),u(is,js,k),v(is,js,k),pe1(k) + enddo + ENDIF + + call p_var(npz, is, ie, js, je, ptop, ptop_min, delp, delz, pt, ps,& + pe, peln, pk, pkz, kappa, q, ng, ncnst, area, dry_mass,.false.,.false., & + .true., hydrostatic, nwat, domain, flagstruct%adiabatic) + + ! Add in (or don't) bubble(s) to initialize storms + if (bubble_type > 0) then + if (is_master()) print*, "ADDING BUBBLE" +! *** Add Initial perturbation *** + pturb = bubble_t + xradbub = bubble_rad_x + yradbub = bubble_rad_y + zc = bubble_zc ! center of bubble from surface + if (bubble_type == 1) then + ! One bubble in domain center + icenter = (npx-1)/2 + 1 + jcenter = (npy-1)/2 + 1 + n_bub = 1 + elseif (bubble_type == 2) then + !Line of centered N-S bubbles for squall line + !n_bub = floor(float(npy)*dy_const/r0) + if ( is_master() ) print*, "initializing ", n_bub , " bubbles" + icenter = 0 + elseif (bubble_type == 3) then + ! User entry of i/j bubble locations + if ( is_master() ) print*, "initializing ", n_bub , " bubbles" + if ( is_master() ) print*, "at locations i = ", icenters(1:n_bub), "j = ", jcenters(1:n_bub) + endif + do j = js, je + do i = is, ie + do k=1, npz + do b = 1, n_bub + call random_number(rand1) + call random_number(rand2) + if (bubble_type == 2) then + jcenter=((npy-1)/2+1)-((n_bub+1)/2-b)*30000.0/dy_const + elseif (bubble_type == 3) then + icenter = icenters(b) + jcenter = jcenters(b) + endif + zm = 0.5*(ze1(k)+ze1(k+1)) + yrad = dy_const*float(j-jcenter)/yradbub + xrad = dx_const*float(i-icenter)/xradbub + + RAD=SQRT(xrad*xrad+yrad*yrad+zrad*zrad) + IF(RAD <= 1.) THEN + if (do_rand_perts) then + ! Add in random, small amplitude perturbations to bubble thermodynamic state + pt(i,j,k) = pt(i,j,k) + pturb*COS(.5*pi*RAD)**2 + 0.2 * (2.0*rand1-1.0) + q(i,j,k,1) = q(i,j,k,1) + bubble_q *COS(.5*pi*RAD)**2 + 1.0E-7 *(2.0*rand2-1.0) + else + pt(i,j,k) = pt(i,j,k) + pturb*COS(.5*pi*RAD)**2 + q(i,j,k,1) = q(i,j,k,1) + bubble_q *COS(.5*pi*RAD)**2 + endif + ENDIF + enddo !nbub + enddo!npz + enddo!i + enddo!j + endif !bubbletype + + uc(isd:ied,:,:) = u(:,jsd:jed,:) + uc(ied+1,:,:) = u(ied,jsd:jed,:) + ua(:,:,:) = u(:,jsd:jed,:) + + vc(:,jsd:jed,:) = v(isd:ied,:,:) + vc(:,jed+1,:) = v(isd:ied,jed,:) + va(:,:,:) = v(isd:ied,:,:) + + nullify(grid) + nullify(agrid) + + nullify(area) + + nullify(fC) + nullify(f0) + + nullify(ee1) + nullify(ee2) + nullify(ew) + nullify(es) + nullify(en1) + nullify(en2) + + nullify(dx) + nullify(dy) + nullify(dxa) + nullify(dya) + nullify(rdxa) + nullify(rdya) + nullify(dxc) + nullify(dyc) + + nullify(dx_const) + nullify(dy_const) + + nullify(domain) + nullify(tile) + + nullify(have_south_pole) + nullify(have_north_pole) + + nullify(ntiles_g) + nullify(acapN) + nullify(globalarea) + + end subroutine fv_init_ideal + + subroutine read_namelist_fv_ideal(nml_filename) + + character(*), intent(IN) :: nml_filename + integer :: ierr, f_unit, unit, ios + namelist /fv_ideal_nml/bubble_do, & + t_profile, q_profile, ws_profile, bubble_t, bubble_q, & + bubble_zc, do_coriolis, iso_t, adi_th, us0, bubble_type,n_bub, & + icenters,jcenters, bubble_rad_x, bubble_rad_y, do_rand_perts, & + p00_in, umove, vmove + +#include + + unit = stdlog() + + ! Make alpha = 0 the default: + alpha = 0. + bubble_do = .false. + + ! Read Test_Case namelist + read (input_nml_file,fv_ideal_nml,iostat=ios) + ierr = check_nml_error(ios,'fv_ideal_nml') + write(unit, nml=fv_ideal_nml) + + + end subroutine read_namelist_fv_ideal + + subroutine get_sounding( zk, p, t, rho, u, v, qv, nl_max, nl_in, p00 ) + implicit none + + integer nl_max, nl_in + real p00 + real zk(nl_max), p(nl_max), theta(nl_max), rho(nl_max), & + u(nl_max), v(nl_max), qv(nl_max), p_dry(nl_max), t(nl_max) + + integer n + parameter(n=1000) + logical debug + parameter( debug = .true.) + character*256 message + +! input sounding data + + real p_surf, th_surf, qv_surf + real pi_surf, pi(n) + real h_input(n), th_input(n), qv_input(n), u_input(n), v_input(n) + +! diagnostics + + real rho_surf, p_input(n), rho_input(n) + real pm_input(n) ! this are for full moist sounding + +! local data + + real r, g,cp + parameter (r = rdgas) + parameter (g = grav) + parameter (cp = cp_air) + real, parameter :: p1000mb = 100000.0, cvpm = -718./cp_air + real, parameter :: rvovrd = rvgas/rdgas + integer k, it, nl, nl_file, istat + real qvf, qvf1, dz + character*256 line + +! first, read the sounding + print*,"OPEN SOUNDING FILE: ", trim('input_sounding') + open(14, file=trim('input_sounding'), form='formatted', iostat=istat) + if (istat /= 0) then + call mpp_error(FATAL,"ERROR OPENING VARIABLE MAPPING FILE") + endif + + nl = 0 + nl_file = 0 + + !Loop over lines of file to count the number of levels + do + read(14, '(A)', iostat=istat) line + if (istat/=0) exit + if ( trim(line) .eq. '' ) cycle + nl_file = nl_file + 1 + enddo + if ( nl_file == 0) call mpp_error(FATAL,"VARMAP FILE IS EMPTY.") + nl = nl_file -1 + + nl_in = nl + if(nl_in .gt. nl_max ) then + print*, ' too many levels for input arrays ',nl_in,nl_max + call mpp_error (FATAL, 'Too many levels for input arrays ' ) + end if + + rewind(14) + read(14,*,iostat=istat) p_surf, th_surf, qv_surf + do k = 2,nl_file + read(14, *, iostat=istat) h_input(nl_in-k+2), th_input(nl_in-k+2), qv_input(nl_in-k+2), u_input(nl_in-k+2), v_input(nl_in-k+2) + if (istat /= 0) call mpp_error(FATAL,"READING VARIABLE MAPPING FILE") + enddo + close(14) + + ! compute diagnostics, + ! first, convert qv(g/kg) to qv(g/g) + + if (qv_surf > 0.1) THEN + do k=1,nl + qv_input(k) = 0.001*qv_input(k) + enddo + endif + + do k=1,nl + IF ( qv_input(k) < 1E-9 ) THEN + write(*,*) 'Warning Input sounding has qv = 0, resetting to 1E-9 kg/kg' + qv_input(k) = 1E-9 ! Max(0.00001,qv_input(k)) + ENDIF + enddo + + IF ( p_surf < 2000. ) p_surf = 100.*p_surf ! convert to pascals + p00 = p_surf + qvf = 1. + rvovrd*qv_input(1) + rho_surf = 1./((r/p1000mb)*th_surf*qvf*((p_surf/p1000mb)**cvpm)) + pi_surf = (p_surf/p1000mb)**(r/cp) + + + + ! integrate moist sounding hydrostatically, starting from the + ! specified surface pressure + ! -> first, integrate from surface to lowest level + + qvf = 1. + rvovrd*qv_input(nl) + qvf1 = 1. + qv_input(nl) + rho_input(nl) = rho_surf + dz = h_input(nl) + do it=1,10 + pm_input(nl) = p_surf & + - 0.5*dz*(rho_surf+rho_input(nl))*g*qvf1 + rho_input(nl) = 1./((r/p1000mb)*th_input(nl)*qvf*((pm_input(nl)/p1000mb)**cvpm)) + enddo + do k=nl-1,1,-1 + rho_input(k) = rho_input(k+1) + dz = h_input(k)-h_input(k+1) + qvf1 = 0.5*(2.+(qv_input(k+1)+qv_input(k))) + qvf = 1. + rvovrd*qv_input(k) ! qv is in g/kg here + + do it=1,10 + pm_input(k) = pm_input(k+1) & + - 0.5*dz*(rho_input(k)+rho_input(k+1))*g*qvf1 + IF(pm_input(k) .LE. 0. )THEN + print*, "Integrated pressure has gone negative - toocold for chosen height" + WRITE(message,*)'k,pm_input(k),h_input(k),th_input(k) =',k,pm_input(k),h_input(k),th_input(k) + CALL mpp_error (FATAL, message ) + ENDIF + rho_input(k) =1./((r/p1000mb)*th_input(k)*qvf*((pm_input(k)/p1000mb)**cvpm)) + enddo + enddo + + do k=1,nl + + zk(k) = h_input(k) + p(k) = pm_input(k) + t(k) = th_input(k) * (pm_input(k) / p1000mb)**(rdgas/cp_air) + u(k) = u_input(k) - umove + v(k) = v_input(k) - vmove + if(is_master()) print*, zk(k) + qv(k) = qv_input(k) + + enddo + end subroutine get_sounding + + real function interp_log( v_in, p_in, p_out, nzmax,nz_in ) + implicit none + integer, intent(in) :: nz_in, nzmax + real, intent(in) :: v_in(nzmax), p_in(nzmax) + real, intent(in) :: p_out + + integer kp, k, im, ip, nz_out + logical interp, increasing_z + real pres + double precision :: w1, w2, v1, v2 + logical debug + parameter ( debug = .false. ) + + pres = p_out + + IF (pres > p_in(nz_in)) then + ! if (is_master()) print*,"interp_log 1: p_in(nz_in), p_in(1) = " , p_in(nz_in), p_in(1) + w2 = (log(p_in(nz_in))-log(pres))/(log(p_in(nz_in))-log(p_in(nz_in-1))) + w1 = 1.-w2 + interp_log = v_in(nz_in)**w1 * v_in(nz_in-1)**w2 + ELSE IF (pres < p_in(1)) then ! extrapolate to lower pressure + w2 = (log(p_in(2))-log(pres))/(log(p_in(2))-log(p_in(1))) + w1 = 1.-w2 + v1 = v_in(1) + v2 = v_in(2) + ! interp_log = w1*v_in(2) + w2*v_in(1) ! + interp_log = dble(v_in(2))**w1 * dble(v_in(1))**w2 + ! if (is_master()) print*,"interp_log 2: p_in(nz_in), p_in(1) = " , p_in(nz_in), p_in(1),pres, & + ! log(p_in(2))-log(pres), log(p_in(2))-log(p_in(1)), w2, w1, v_in(2) , v_in(1), & + ! v_in(2)**w1, v_in(1)**w2, v2**w1 * v1**w2, w1*v_in(2) + w2*v_in(1), interp_log + ELSE + interp = .false. + kp = nz_in + DO WHILE ( (interp .eqv. .false.) .and. (kp .ge. 2) ) + IF( ((p_in(kp) .ge. pres) .and. & + (p_in(kp-1) .le. pres)) ) THEN + w2 = (log(pres)-log(p_in(kp)))/(log(p_in(kp-1))-log(p_in(kp))) + w1 = 1.-w2 + interp_log = v_in(kp)**w1 * v_in(kp-1)**w2 + interp = .true. + ! if (is_master()) print*,"interp_log 3: pres,p(kp),p(kp-1),w2,w1 = " , pres, p_in(kp), p_in(kp-1), & + ! w2, w1, v_in(kp) , v_in(kp-1), interp_log, w1*v_in(kp) + w2*v_in(kp-1) + + END IF + kp = kp-1 + ENDDO + ENDIF + + end function interp_log + + real function interp_lin( v_in, z_in, z_out, nzmax,nz_in ) + implicit none + integer, intent(in) :: nz_in, nzmax + real, intent(in) :: v_in(nzmax), z_in(nzmax) + real, intent(in) :: z_out + + integer kp, k, im, ip, nz_out + logical interp, increasing_z + real height, w1, w2 + logical debug + parameter ( debug = .true. ) + + height = z_out + IF (height < z_in(nz_in)) then + if(debug .and. is_master()) print*, ' point 1 in interp_lin ', height + w2 = (z_in(nz_in)-height)/(z_in(nz_in)-z_in(nz_in-1)) + w1 = 1.-w2 + interp_lin = w1*v_in(nz_in) + w2*v_in(nz_in-1) + ELSE IF (height > z_in(1)) then + if(debug .and. is_master()) print*, ' point 2 in interp_lin ', height + w2 = (z_in(2)-height)/(z_in(2)-z_in(1)) + w1 = 1.-w2 + interp_lin = w1*v_in(2) + w2*v_in(1) + ELSE + if(debug .and. is_master()) print*, ' point 3 in interp_lin ', height + interp = .false. + kp = nz_in + height = z_out + DO WHILE ( (interp .eqv. .false.) .and. (kp .ge. 2) ) + IF( ((z_in(kp) .le. height) .and. & + (z_in(kp-1) .ge. height)) ) THEN + w2 = (height-z_in(kp))/(z_in(kp-1)-z_in(kp)) + w1 = 1.-w2 + interp_lin = w1*v_in(kp) + w2*v_in(kp-1) + interp = .true. + END IF + kp = kp-1 + ENDDO + ENDIF + end function interp_lin + + subroutine SuperCell_Sounding(km, ps, pk1, tp, qp) + ! Morris Weisman & J. Klemp 2002 sounding + ! Output sounding on pressure levels: + integer, intent(in):: km + real, intent(in):: ps ! surface pressure (Pa) + real, intent(in), dimension(km):: pk1 + real, intent(out), dimension(km):: tp, qp + ! Local: + integer, parameter:: ns = 401 + integer, parameter:: nx = 3 + real, dimension(ns):: zs, pt, qs, us, rh, pp, pk, dpk, dqdt + real, parameter:: Tmin = 175. + real, parameter:: p00 = 1.0e5 + real, parameter:: qst = 3.0e-6 + real, parameter:: qv0 = 1.4e-2 + real, parameter:: ztr = 12.E3 + real, parameter:: ttr = 213. + real, parameter:: ptr = 343. ! Tropopause potential temp. + real, parameter:: pt0 = 300. ! surface potential temperature + real:: dz0, zvir, fac_z, pk0, temp1, p2 + integer:: k, n, kk + + ! tp=0. + ! qp=0. + + !!#else + + zvir = rvgas/rdgas - 1. + pk0 = p00**kappa + pp(ns) = ps + pk(ns) = ps**kappa + if ( (is_master()) ) then + write(*,*) 'Computing sounding for super-cell test' + endif + + dz0 = 50. + zs(ns) = 0. + qs(:) = qst + rh(:) = 0.25 + + do k=ns-1, 1, -1 + zs(k) = zs(k+1) + dz0 + enddo + + do k=1,ns + ! Potential temperature + if ( zs(k) .gt. ztr ) then + ! Stratosphere: + pt(k) = ptr*exp(grav*(zs(k)-ztr)/(cp_air*ttr)) + else + ! Troposphere: + fac_z = (zs(k)/ztr)**1.25 + pt(k) = pt0 + (ptr-pt0)* fac_z + rh(k) = 1. - 0.75 * fac_z + ! First guess on q: + qs(k) = qv0 - (qv0-qst)*fac_z + endif + pt(k) = pt(k) / pk0 + enddo + + !-------------------------------------- + ! Iterate nx times with virtual effect: + !-------------------------------------- + do n=1, nx + do k=1,ns-1 + temp1 = 0.5*(pt(k)*(1.+zvir*qs(k)) + pt(k+1)*(1.+zvir*qs(k+1))) + dpk(k) = grav*(zs(k)-zs(k+1))/(cp_air*temp1) ! DPK > 0 + enddo + + do k=ns-1,1,-1 + pk(k) = pk(k+1) - dpk(k) + enddo + + do k=1, ns + temp1 = pt(k)*pk(k) + ! if ( (is_master()) ) write(*,*) k, temp1, rh(k) + if ( pk(k) > 0. ) then + pp(k) = exp(log(pk(k))/kappa) + qs(k) = 380./pp(k)*exp(17.27*(temp1-273.)/(temp1-36.)) + qs(k) = min( qv0, rh(k)*qs(k) ) + if ( (is_master()) ) write(*,*) 0.01*pp(k), qs(k) + else + if ( (is_master()) ) write(*,*) n, k, pk(k) + call mpp_error(FATAL, 'Super-Cell case: pk < 0') + endif + enddo + enddo + + ! Interpolate to p levels using pk1: p**kappa + do 555 k=1, km + if ( pk1(k) .le. pk(1) ) then + tp(k) = pt(1)*pk(1)/pk1(k) ! isothermal above + qp(k) = qst ! set to stratosphere value + elseif ( pk1(k) .ge. pk(ns) ) then + tp(k) = pt(ns) + qp(k) = qs(ns) + else + do kk=1,ns-1 + if( (pk1(k).le.pk(kk+1)) .and. (pk1(k).ge.pk(kk)) ) then + fac_z = (pk1(k)-pk(kk))/(pk(kk+1)-pk(kk)) + tp(k) = pt(kk) + (pt(kk+1)-pt(kk))*fac_z + qp(k) = qs(kk) + (qs(kk+1)-qs(kk))*fac_z + goto 555 + endif + enddo + endif + 555 continue + + do k=1,km + tp(k) = tp(k)*pk1(k) ! temperature + tp(k) = max(Tmin, tp(k)) + enddo + + + end subroutine SuperCell_Sounding + + end module fv_ideal_mod diff --git a/driver/UFS/fv_nggps_diag.F90 b/driver/UFS/fv_nggps_diag.F90 index a8b8aa41d..2a257bd25 100644 --- a/driver/UFS/fv_nggps_diag.F90 +++ b/driver/UFS/fv_nggps_diag.F90 @@ -651,7 +651,6 @@ subroutine fv_nggps_diag(Atm, zvir, Time) enddo call store_data(id_delp, wk, Time, kstt_delp, kend_delp) endif - !--- Surface Pressure (PS) ! Re-compute pressure (dry_mass + water_vapor) surface pressure if(id_ps > 0) then diff --git a/input_sounding b/input_sounding new file mode 100644 index 000000000..b309ffb4f --- /dev/null +++ b/input_sounding @@ -0,0 +1,81 @@ +1000.0000 300.5000 14.0000 +125.0000 300.5000 14.0000 0.6250 0.0000000E+00 +375.0000 300.5650 14.0000 1.8750 0.0000000E+00 +625.0000 301.0699 14.0000 3.1250 0.0000000E+00 +875.0000 301.6293 14.0000 4.3750 0.0000000E+00 +1125.0000 302.2307 14.0000 5.6250 0.0000000E+00 +1375.0000 302.8666 13.1191 6.8750 0.0000000E+00 +1625.0000 303.5323 11.8756 8.1250 0.0000000E+00 +1875.0000 304.2242 10.7420 9.3750 0.0000000E+00 +2125.0000 304.9395 9.7079 10.6250 0.0000000E+00 +2375.0000 305.6764 8.7643 11.8750 0.0000000E+00 +2625.0000 306.4329 7.9033 13.1250 0.0000000E+00 +2875.0000 307.2076 7.1179 14.3750 0.0000000E+00 +3125.0000 307.9994 6.4018 15.6250 0.0000000E+00 +3375.0000 308.8071 5.7493 16.8750 0.0000000E+00 +3625.0000 309.6300 5.1553 18.1250 0.0000000E+00 +3875.0000 310.4672 4.6150 19.3750 0.0000000E+00 +4125.0000 311.3181 4.1241 20.6250 0.0000000E+00 +4375.0000 312.1819 3.6788 21.8750 0.0000000E+00 +4625.0000 313.0581 3.2752 23.1250 0.0000000E+00 +4875.0000 313.9464 2.9100 24.3750 0.0000000E+00 +5125.0000 314.8460 2.5800 25.0000 0.0000000E+00 +5375.0000 315.7567 2.2825 25.0000 0.0000000E+00 +5625.0000 316.6780 2.0145 25.0000 0.0000000E+00 +5875.0000 317.6097 1.7738 25.0000 0.0000000E+00 +6125.0000 318.5513 1.5579 25.0000 0.0000000E+00 +6375.0000 319.5026 1.3647 25.0000 0.0000000E+00 +6625.0000 320.4632 1.1922 25.0000 0.0000000E+00 +6875.0000 321.4330 1.0385 25.0000 0.0000000E+00 +7125.0000 322.4116 0.9019 25.0000 0.0000000E+00 +7375.0000 323.3988 0.7809 25.0000 0.0000000E+00 +7625.0000 324.3945 0.6739 25.0000 0.0000000E+00 +7875.0000 325.3983 0.5796 25.0000 0.0000000E+00 +8125.0000 326.4102 0.4967 25.0000 0.0000000E+00 +8375.0000 327.4298 0.4241 25.0000 0.0000000E+00 +8625.0000 328.4571 0.3606 25.0000 0.0000000E+00 +8875.0000 329.4919 0.3054 25.0000 0.0000000E+00 +9125.0000 330.5339 0.2576 25.0000 0.0000000E+00 +9375.0000 331.5832 0.2162 25.0000 0.0000000E+00 +9625.0000 332.6394 0.1806 25.0000 0.0000000E+00 +9875.0000 333.7025 0.1502 25.0000 0.0000000E+00 +10125.0000 334.7725 0.1241 25.0000 0.0000000E+00 +10375.0000 335.8490 0.1020 25.0000 0.0000000E+00 +10625.0000 336.9320 0.0834 25.0000 0.0000000E+00 +10875.0000 338.0214 0.0677 25.0000 0.0000000E+00 +11125.0000 339.1171 0.0545 25.0000 0.0000000E+00 +11375.0000 340.2190 0.0436 25.0000 0.0000000E+00 +11625.0000 341.3269 0.0346 25.0000 0.0000000E+00 +11875.0000 342.4408 0.0272 25.0000 0.0000000E+00 +12125.0000 344.9724 0.0246 25.0000 0.0000000E+00 +12375.0000 348.9513 0.0257 25.0000 0.0000000E+00 +12625.0000 352.9761 0.0269 25.0000 0.0000000E+00 +12875.0000 357.0473 0.0282 25.0000 0.0000000E+00 +13125.0000 361.1655 0.0295 25.0000 0.0000000E+00 +13375.0000 365.3311 0.0309 25.0000 0.0000000E+00 +13625.0000 369.5448 0.0323 25.0000 0.0000000E+00 +13875.0000 373.8072 0.0338 25.0000 0.0000000E+00 +14125.0000 378.1187 0.0354 25.0000 0.0000000E+00 +14375.0000 382.4798 0.0371 25.0000 0.0000000E+00 +14625.0000 386.8913 0.0388 25.0000 0.0000000E+00 +14875.0000 391.3537 0.0406 25.0000 0.0000000E+00 +15125.0000 395.8676 0.0425 25.0000 0.0000000E+00 +15375.0000 400.4335 0.0445 25.0000 0.0000000E+00 +15625.0000 405.0521 0.0467 25.0000 0.0000000E+00 +15875.0000 409.7239 0.0489 25.0000 0.0000000E+00 +16125.0000 414.4497 0.0512 25.0000 0.0000000E+00 +16375.0000 419.2299 0.0536 25.0000 0.0000000E+00 +16625.0000 424.0653 0.0562 25.0000 0.0000000E+00 +16875.0000 428.9565 0.0588 25.0000 0.0000000E+00 +17125.0000 433.9040 0.0616 25.0000 0.0000000E+00 +17375.0000 438.9086 0.0646 25.0000 0.0000000E+00 +17625.0000 443.9710 0.0677 25.0000 0.0000000E+00 +17875.0000 449.0917 0.0709 25.0000 0.0000000E+00 +18125.0000 454.2715 0.0743 25.0000 0.0000000E+00 +18375.0000 459.5110 0.0779 25.0000 0.0000000E+00 +18625.0000 464.8110 0.0816 25.0000 0.0000000E+00 +18875.0000 470.1721 0.0856 25.0000 0.0000000E+00 +19125.0000 475.5951 0.0897 25.0000 0.0000000E+00 +19375.0000 481.0806 0.0940 25.0000 0.0000000E+00 +19625.0000 486.6293 0.0940 25.0000 0.0000000E+00 +19875.0000 492.2421 0.0940 25.0000 0.0000000E+00 diff --git a/model/fv_arrays.F90 b/model/fv_arrays.F90 index 26a9cf153..6a91ed09b 100644 --- a/model/fv_arrays.F90 +++ b/model/fv_arrays.F90 @@ -916,6 +916,12 @@ module fv_arrays_mod real(kind=R_GRID) :: deglat=15. !< Latitude (in degrees) used to compute the uniform f-plane !< Coriolis parameter for doubly-periodic simulations !< (grid_type = 4). The default value is 15. + real(kind=R_GRID) :: deglon=0. !< Longitude (in degrees) used for radiation + !< calculations for doubly-periodic simulations. + !< (grid_type = 4). Default value is 0. + real(kind=R_GRID) :: deg_domain=10. !< Domain size (north-south & east-east) in + !< degrees used to calculate radiation. + !< (grid_type = 4). Default value is 10. !The following deglat_*, deglon_* options are not used. real(kind=R_GRID) :: deglon_start = -30., deglon_stop = 30., & !< boundaries of latlon patch deglat_start = -30., deglat_stop = 30. diff --git a/model/fv_control.F90 b/model/fv_control.F90 index 2ead0274c..99355e4c1 100644 --- a/model/fv_control.F90 +++ b/model/fv_control.F90 @@ -141,6 +141,7 @@ module fv_control_mod use fv_mp_mod, only: mp_start, domain_decomp, mp_assign_gid, global_nest_domain use fv_mp_mod, only: broadcast_domains, mp_barrier, is_master, setup_master, grids_master_procs, tile_fine use fv_mp_mod, only: MAX_NNEST, MAX_NTILE + use fv_ideal_mod, only: read_namelist_fv_ideal use test_cases_mod, only: read_namelist_test_case_nml use fv_timing_mod, only: timing_on, timing_off, timing_init, timing_prt use mpp_domains_mod, only: domain2D @@ -399,6 +400,8 @@ subroutine fv_control_init(Atm, dt_atmos, this_grid, grids_on_this_pe, p_split, real(kind=R_GRID), pointer :: deglon_start, deglon_stop, & ! boundaries of latlon patch deglat_start, deglat_stop real(kind=R_GRID), pointer :: deglat + real(kind=R_GRID), pointer :: deglon + real(kind=R_GRID), pointer :: deg_domain logical, pointer :: nested, twowaynest logical, pointer :: regional @@ -578,7 +581,11 @@ subroutine fv_control_init(Atm, dt_atmos, this_grid, grids_on_this_pe, p_split, call read_namelist_molecular_diffusion_nml(Atm(this_grid)%nml_filename, & Atm(this_grid)%flagstruct%ncnst, Atm(this_grid)%flagstruct%nwat) endif + + call read_namelist_fv_ideal(Atm(this_grid)%nml_filename) + call read_namelist_test_case_nml(Atm(this_grid)%nml_filename) + call mpp_get_current_pelist(Atm(this_grid)%pelist, commID=commID) ! for commID call mp_start(commID,halo_update_type) @@ -978,6 +985,8 @@ subroutine set_namelist_pointers(Atm) deglat_stop => Atm%flagstruct%deglat_stop deglat => Atm%flagstruct%deglat + deglon => Atm%flagstruct%deglon + deg_domain => Atm%flagstruct%deg_domain nested => Atm%neststruct%nested twowaynest => Atm%neststruct%twowaynest @@ -1083,7 +1092,7 @@ subroutine read_namelist_fv_core_nml(Atm) tau, tau_w, fast_tau_w_sec, tau_h2o, rf_cutoff, nf_omega, hydrostatic, fv_sg_adj, sg_cutoff, breed_vortex_inline, & na_init, nudge_dz, hybrid_z, Make_NH, n_zs_filter, nord_zs_filter, full_zs_filter, reset_eta, & pnats, dnats, dnrts, a2b_ord, remap_t, p_ref, d2_bg_k1, d2_bg_k2, & - c2l_ord, dx_const, dy_const, umax, deglat, & + c2l_ord, dx_const, dy_const, umax, deglat, deglon, deg_domain, & deglon_start, deglon_stop, deglat_start, deglat_stop, & phys_hydrostatic, use_hydro_pressure, make_hybrid_z, old_divg_damp, add_noise, butterfly_effect, & molecular_diffusion, dz_min, psm_bc, nested, twowaynest, nudge_qv, & diff --git a/model/fv_dynamics.F90 b/model/fv_dynamics.F90 index 36b048136..251a0aff4 100644 --- a/model/fv_dynamics.F90 +++ b/model/fv_dynamics.F90 @@ -1420,6 +1420,8 @@ subroutine Rayleigh_Friction(dt, npx, npy, npz, ks, pm, tau, u, v, w, pt, & if ( .not. RF_initialized ) then allocate( rf(npz) ) + rf(:) = 0.0 + kmax = 0 if( is_master() ) write(6,*) 'Rayleigh friction E-folding time (days):' do k=1, npz if ( pm(k) < rf_cutoff ) then @@ -1433,6 +1435,8 @@ subroutine Rayleigh_Friction(dt, npx, npy, npz, ks, pm, tau, u, v, w, pt, & RF_initialized = .true. endif + if ( kmax == 0 ) return + allocate( u2f(isd:ied,jsd:jed,kmax) ) call c2l_ord2(u, v, ua, va, gridstruct, npz, gridstruct%grid_type, bd, gridstruct%bounded_domain) diff --git a/tools/external_ic.F90 b/tools/external_ic.F90 index 157f0467f..14ff72ca2 100644 --- a/tools/external_ic.F90 +++ b/tools/external_ic.F90 @@ -154,7 +154,17 @@ module external_ic_mod use mpp_mod, only: mpp_error, FATAL, NOTE, mpp_pe, mpp_root_pe use mpp_mod, only: stdlog, input_nml_file, mpp_npes, mpp_get_current_pelist use mpp_parameter_mod, only: AGRID_PARAM=>AGRID - use mpp_domains_mod, only: mpp_get_tile_id, domain2d, mpp_update_domains, NORTH, EAST +#ifdef ENABLE_PARALLELRESTART + use mpp_domains_mod, only: mpp_get_tile_id, domain2d, mpp_update_domains, & + NORTH, EAST, mpp_get_domain_tile_commid, & + mpp_get_io_domain_layout, mpp_copy_domain, & + mpp_define_io_domain, mpp_get_layout +#else + use mpp_domains_mod, only: mpp_get_tile_id, domain2d, mpp_update_domains, & + NORTH, EAST, & + mpp_get_io_domain_layout, mpp_copy_domain, & + mpp_define_io_domain, mpp_get_layout +#endif use tracer_manager_mod, only: get_tracer_names, get_number_tracers, get_tracer_index use tracer_manager_mod, only: set_tracer_profile use field_manager_mod, only: MODEL_ATMOS @@ -875,6 +885,8 @@ subroutine get_nggps_ic (Atm) contains subroutine read_gfs_ic() + integer :: read_layout(2) + type(domain2D) :: domain_for_read ! !--- read in ak and bk from the gfs control file using fms_io read_data --- ! @@ -916,7 +928,17 @@ subroutine read_gfs_ic() dim_names_3d4(1) = "levp" ! surface pressure (Pa) +#ifdef ENABLE_PARALLELRESTART + call mpp_get_layout(Atm%domain,read_layout) + call mpp_copy_domain(Atm%domain, domain_for_read) + call mpp_define_io_domain(domain_for_read, read_layout) + + GFS_restart%use_collective = .true. + GFS_restart%tile_comm = mpp_get_domain_tile_commid(Atm%domain) + if( open_file(GFS_restart, fn_gfs_ics, "read", domain_for_read, is_restart=.true., dont_add_res_to_filename=.true.) ) then +#else if( open_file(GFS_restart, fn_gfs_ics, "read", Atm%domain_for_read, is_restart=.true., dont_add_res_to_filename=.true.) ) then +#endif call register_axis(GFS_restart, "lat", "y") call register_axis(GFS_restart, "lon", "x") call register_axis(GFS_restart, "lonp", "x", domain_position=east) diff --git a/tools/fv_grid_tools.F90 b/tools/fv_grid_tools.F90 index 67b33ccf2..9fcd8f0c8 100644 --- a/tools/fv_grid_tools.F90 +++ b/tools/fv_grid_tools.F90 @@ -730,7 +730,7 @@ subroutine init_grid(Atm, grid_name, grid_file, npx, npy, npz, ndims, nregions, if (Atm%flagstruct%grid_type == 4) then call setup_cartesian(npx, npy, Atm%flagstruct%dx_const, Atm%flagstruct%dy_const, & - Atm%flagstruct%deglat, Atm%bd) + Atm%flagstruct%deglat, Atm%flagstruct%deglon, Atm%flagstruct%deg_domain, Atm%bd) elseif (Atm%flagstruct%grid_type == 5) then call setup_orthogonal_grid(npx, npy, Atm%bd, grid_file) else @@ -1472,11 +1472,11 @@ subroutine init_grid(Atm, grid_name, grid_file, npx, npy, npz, ndims, nregions, contains - subroutine setup_cartesian(npx, npy, dx_const, dy_const, deglat, bd) + subroutine setup_cartesian(npx, npy, dx_const, dy_const, deglat, deglon, deg_domain, bd) type(fv_grid_bounds_type), intent(IN) :: bd integer, intent(in):: npx, npy - real(kind=R_GRID), intent(IN) :: dx_const, dy_const, deglat + real(kind=R_GRID), intent(IN) :: dx_const, dy_const, deglat, deglon, deg_domain real(kind=R_GRID) lat_rad, lon_rad, domain_rad integer i,j integer :: is, ie, js, je @@ -1491,9 +1491,9 @@ subroutine setup_cartesian(npx, npy, dx_const, dy_const, deglat, bd) jsd = bd%jsd jed = bd%jed - domain_rad = pi/16. ! arbitrary + domain_rad = deg_domain * pi/180. !pi/16. ! arbitrary lat_rad = deglat * pi/180. - lon_rad = 0. ! arbitrary + lon_rad = deglon * pi/180. !0. ! arbitrary dx(:,:) = dx_const rdx(:,:) = 1./dx_const diff --git a/tools/fv_io.F90 b/tools/fv_io.F90 index 74d6e1ce4..e050db220 100644 --- a/tools/fv_io.F90 +++ b/tools/fv_io.F90 @@ -94,11 +94,19 @@ module fv_io_mod variable_exists, read_data, set_filename_appendix, get_dimension_size use mpp_mod, only: mpp_error, FATAL, NOTE, WARNING, mpp_root_pe, & mpp_sync, mpp_pe, mpp_declare_pelist, mpp_get_current_pelist, & +#ifdef ENABLE_PARALLELRESTART + mpp_npes, MPP_COMM_NULL +#else mpp_npes +#endif use mpp_domains_mod, only: domain2d, EAST, WEST, NORTH, CENTER, SOUTH, CORNER, & mpp_get_compute_domain, mpp_get_data_domain, & - mpp_get_layout, mpp_get_ntile_count, & + mpp_get_layout, mpp_get_ntile_count, mpp_copy_domain, & +#ifdef ENABLE_PARALLELRESTART + mpp_get_global_domain, mpp_get_domain_tile_commid, mpp_define_io_domain +#else mpp_get_global_domain +#endif use tracer_manager_mod, only: tr_get_tracer_names=>get_tracer_names, & get_tracer_names, get_number_tracers, & set_tracer_profile, & @@ -438,6 +446,8 @@ subroutine fv_io_read_restart(fv_domain,Atm) character(len=20) :: suffix character(len=1) :: tile_num integer, allocatable, dimension(:) :: pes !< Array of the pes in the current pelist + type(domain2D) :: domain_for_read + integer :: read_layout(2), layout(2) allocate(pes(mpp_npes())) call mpp_get_current_pelist(pes) @@ -451,7 +461,7 @@ subroutine fv_io_read_restart(fv_domain,Atm) call close_file(Atm(1)%Fv_restart) Atm(1)%Fv_restart_is_open = .false. endif - deallocate(pes) + !deallocate(pes) if (Atm(1)%flagstruct%external_eta) then call set_external_eta(Atm(1)%ak, Atm(1)%bk, Atm(1)%ptop, Atm(1)%ks) @@ -469,17 +479,38 @@ subroutine fv_io_read_restart(fv_domain,Atm) endif fname = 'INPUT/fv_core.res'//trim(suffix)//'.nc' +#ifdef ENABLE_PARALLELRESTART + Atm(1)%Fv_restart_tile%use_collective = .true. + + call mpp_get_layout(Atm(1)%domain,read_layout) + call mpp_copy_domain(Atm(1)%domain, domain_for_read) + call mpp_define_io_domain(domain_for_read, read_layout) + Atm(1)%Fv_restart_tile%tile_comm = mpp_get_domain_tile_commid(Atm(1)%domain) + Atm(1)%Fv_restart_tile_is_open = open_file(Atm(1)%Fv_restart_tile, fname, "read", domain_for_read, is_restart=.true.) +#else Atm(1)%Fv_restart_tile_is_open = open_file(Atm(1)%Fv_restart_tile, fname, "read", fv_domain, is_restart=.true.) +#endif if (Atm(1)%Fv_restart_tile_is_open) then call fv_io_register_restart(Atm(1)) call read_restart(Atm(1)%Fv_restart_tile, ignore_checksum=Atm(1)%flagstruct%ignore_rst_cksum) call close_file(Atm(1)%Fv_restart_tile) Atm(1)%Fv_restart_tile_is_open = .false. endif +#ifdef ENABLE_PARALLELRESTART + ! Disable use_collective so the reuse of this Fv_restart_tile during checkpointing doesn't fail + Atm(1)%Fv_restart_tile%use_collective = .false. + Atm(1)%Fv_restart_tile%tile_comm = MPP_COMM_NULL +#endif !--- restore data for fv_tracer - if it exists fname = 'INPUT/fv_tracer.res'//trim(suffix)//'.nc' +#ifdef ENABLE_PARALLELRESTART + Atm(1)%Tra_restart%use_collective = .true. + Atm(1)%Tra_restart%tile_comm = mpp_get_domain_tile_commid(Atm(1)%domain) + Atm(1)%Tra_restart_is_open = open_file(Atm(1)%Tra_restart, fname, "read", domain_for_read, is_restart=.true.) +#else Atm(1)%Tra_restart_is_open = open_file(Atm(1)%Tra_restart, fname, "read", fv_domain, is_restart=.true.) +#endif if (Atm(1)%Tra_restart_is_open) then call fv_io_register_restart(Atm(1)) call read_restart(Atm(1)%Tra_restart, ignore_checksum=Atm(1)%flagstruct%ignore_rst_cksum) @@ -488,6 +519,11 @@ subroutine fv_io_read_restart(fv_domain,Atm) else call mpp_error(NOTE,'==> Warning from fv_read_restart: Expected file '//trim(fname)//' does not exist') endif +#ifdef ENABLE_PARALLELRESTART + ! Disable use_collective so the reuse of this Tra_restart during checkpointing doesn't fail + Atm(1)%Tra_restart%use_collective = .false. + Atm(1)%Fv_restart_tile%tile_comm = MPP_COMM_NULL +#endif !--- restore data for surface winds - if it exists fname = 'INPUT/fv_srf_wnd.res'//trim(suffix)//'.nc' @@ -528,6 +564,7 @@ subroutine fv_io_read_restart(fv_domain,Atm) endif endif + deallocate(pes) return end subroutine fv_io_read_restart diff --git a/tools/fv_mp_mod.F90 b/tools/fv_mp_mod.F90 index 4f2c0f2c2..bf6b1ede4 100644 --- a/tools/fv_mp_mod.F90 +++ b/tools/fv_mp_mod.F90 @@ -658,12 +658,16 @@ subroutine domain_decomp(grid_num,npx,npy,nregions,grid_type,nested,layout,io_la !--- if io_layout\=(1,1) then read io_layout=io_layout (no change) l_layout = mpp_get_io_domain_layout(domain) call mpp_copy_domain(domain, domain_for_read) +#ifdef ENABLE_RRFS_WAR if (ALL(l_layout == 1)) then call mpp_get_layout(domain, l_layout) call mpp_define_io_domain(domain_for_read, l_layout) else call mpp_define_io_domain(domain_for_read, l_layout) endif +#else + call mpp_define_io_domain(domain_for_read, l_layout) +#endif endif deallocate(pe_start,pe_end) diff --git a/tools/fv_restart.F90 b/tools/fv_restart.F90 index dae3d0496..2b356f78a 100644 --- a/tools/fv_restart.F90 +++ b/tools/fv_restart.F90 @@ -173,6 +173,9 @@ module fv_restart_mod #ifdef MULTI_GASES use multi_gases_mod, only: virq #endif +#ifdef GFS_PHYS + use fv_ideal_mod, only: fv_init_ideal +#endif implicit none private @@ -434,6 +437,22 @@ subroutine fv_restart(fv_domain, Atm, seconds, days, cold_start, grid_type, this Atm(n)%ptop, Atm(n)%domain, Atm(n)%tile_of_mosaic, Atm(n)%bd) endif elseif (grid_type == 4) then +#ifdef GFS_PHYS + print*, "CALLING FV_INIT_IDEAL" + call fv_init_ideal(Atm(n)%u,Atm(n)%v,Atm(n)%w,Atm(n)%pt, & + Atm(n)%delp,Atm(n)%q,Atm(n)%phis, Atm(n)%ps,Atm(n)%pe, & + Atm(n)%peln,Atm(n)%pk,Atm(n)%pkz, & + Atm(n)%uc,Atm(n)%vc, Atm(n)%ua,Atm(n)%va, & + Atm(n)%ak, Atm(n)%bk, & + Atm(n)%gridstruct, Atm(n)%flagstruct, & + Atm(n)%npx, Atm(n)%npy, npz, Atm(n)%ng, & + ncnst, Atm(n)%flagstruct%nwat, & + Atm(n)%flagstruct%ndims, Atm(n)%flagstruct%ntiles, & + Atm(n)%flagstruct%dry_mass, Atm(n)%flagstruct%mountain, & + Atm(n)%flagstruct%moist_phys, Atm(n)%flagstruct%hydrostatic, & + hybrid, Atm(n)%delz, Atm(n)%ze0, Atm(n)%ks, Atm(n)%ptop, & + Atm(n)%domain, Atm(n)%tile_of_mosaic, Atm(n)%bd) +#else call init_double_periodic(Atm(n)%u,Atm(n)%v,Atm(n)%w,Atm(n)%pt, & Atm(n)%delp,Atm(n)%q,Atm(n)%phis, Atm(n)%ps,Atm(n)%pe, & Atm(n)%peln,Atm(n)%pk,Atm(n)%pkz, & @@ -447,6 +466,7 @@ subroutine fv_restart(fv_domain, Atm, seconds, days, cold_start, grid_type, this Atm(n)%flagstruct%moist_phys, Atm(n)%flagstruct%hydrostatic, & hybrid, Atm(n)%delz, Atm(n)%ze0, Atm(n)%ks, Atm(n)%ptop, & Atm(n)%domain, Atm(n)%tile_of_mosaic, Atm(n)%bd) +#endif if( is_master() ) write(*,*) 'Doubly Periodic IC generated' elseif (grid_type == 5 .or. grid_type == 6) then call mpp_error(FATAL, "Idealized test cases for grid_type == 5,6 (global lat-lon) grid not supported") @@ -649,6 +669,7 @@ subroutine fv_restart(fv_domain, Atm, seconds, days, cold_start, grid_type, this sin(Atm(n)%gridstruct%agrid(i,j,2))*cos(alpha) ) enddo enddo +#ifndef GFS_PHYS else f00 = 2.*omega*sin(Atm(n)%flagstruct%deglat/180.*pi) do j=jsd,jed+1 @@ -661,11 +682,13 @@ subroutine fv_restart(fv_domain, Atm, seconds, days, cold_start, grid_type, this Atm(n)%gridstruct%f0(i,j) = f00 enddo enddo +#endif endif call mpp_update_domains( Atm(n)%gridstruct%f0, Atm(n)%domain ) if ( Atm(n)%gridstruct%cubed_sphere .and. (.not. Atm(n)%gridstruct%bounded_domain))then call fill_corners(Atm(n)%gridstruct%f0, Atm(n)%npx, Atm(n)%npy, Corners_YDir) endif + endif