diff --git a/CMakeLists.txt b/CMakeLists.txt index 2e079bf38..adf6e812d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -114,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_processmodel.F90) list(APPEND fv3_srcs ${model_srcs} ${tools_srcs}) diff --git a/docs/ideal_processmodel_doc/README.md b/docs/ideal_processmodel_doc/README.md new file mode 100644 index 000000000..e33a4cc94 --- /dev/null +++ b/docs/ideal_processmodel_doc/README.md @@ -0,0 +1,83 @@ +Some notes on the ideal "Process Model" setup: + +The idealized "Process Model" provides the option to initialize FV3 with horizontally homogeneous environment in a doubly-periodic domain with any physics suite. The environmental setup is done through the fv_processmodel_nml namelist. This option is separate and distinct from the test cases in test_case.F90 that provide specific dycore tests. Here we provide the ability for the user to choose their 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_processmodel.F90 in the atmos_cubed_sphere directory for code specifics. + +A self-contained experiment directory (working as of initial release) can be downloaded from + https://github.com/MicroTed/ufs_ideal_run.git + +Example namelist settings: + + &fv_processmodel_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.) NOTE: These random perts are not currently reproducible under different domain decompositions (i.e., different number of MPI threads) + 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 physics suite 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_processmodel_doc/create_ideal_sfc_oro_input.py b/docs/ideal_processmodel_doc/create_ideal_sfc_oro_input.py new file mode 100755 index 000000000..9f31b00cf --- /dev/null +++ b/docs/ideal_processmodel_doc/create_ideal_sfc_oro_input.py @@ -0,0 +1,210 @@ +#!/usr/bin/env python3 +# NOTE: requires the netcdf4 package, which is not standard on system installs. +# On hera, try /home/Ted.Mansell/miniconda3/bin/python3 which has netcdf installed, or +# otherwise install netcdf in your own miniconda +# +# Purpose: Creates sfc_data.nc and oro_data.nc for UFS ideal setup with user-specified +# values (affecting land surface, radiation, etc.) +# 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 = 'sfc_data.nc' +orofilename = 'oro_data.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/docs/ideal_processmodel_doc/input_sounding b/docs/ideal_processmodel_doc/input_sounding new file mode 100644 index 000000000..b309ffb4f --- /dev/null +++ b/docs/ideal_processmodel_doc/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/driver/UFS/atmosphere.F90 b/driver/UFS/atmosphere.F90 index 61eb52d28..01206d078 100644 --- a/driver/UFS/atmosphere.F90 +++ b/driver/UFS/atmosphere.F90 @@ -973,13 +973,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 @@ -993,6 +994,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_processmodel.F90 b/driver/UFS/fv_processmodel.F90 new file mode 100644 index 000000000..9cfb91ad8 --- /dev/null +++ b/driver/UFS/fv_processmodel.F90 @@ -0,0 +1,912 @@ +!*********************************************************************** +!* 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 . +!*********************************************************************** + +!*********************************************************************** +! Initialization module for "process model" capability. These routines can read or +! generate a single sounding for a horizontally homogeneous environment (doubly- +! periodic) and introduce thermal perturbations (bubbles) to initiate convection, +! if desired. Eventually can add other capabilities. +! Namelist of input parameters is fv_processmodel_nml +!*********************************************************************** + + module fv_processmodel_mod +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! > +! +! +! +! +!
Module NameFunctions Included
constants_modcnst_radius=>radius, pi=>pi_8, omega, grav, kappa, rdgas, cp_air, rvgas
init_hydro_modp_var
fms_modcheck_nml_error
fv_arrays_modfv_grid_type, fv_flags_type, fv_grid_bounds_type, R_GRID
fv_grid_utils_modptop_min
fv_mp_modis_master
mpp_domains_moddomain2d
mpp_modmpp_error, FATAL, stdlog, input_nml_file
+ +#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 fms_mod, only: check_nml_error + use fv_arrays_mod, only: fv_grid_type, fv_flags_type, fv_grid_bounds_type, R_GRID + use fv_grid_utils_mod, only: ptop_min + use fv_mp_mod, only: is_master + use mpp_domains_mod, only: domain2d + use mpp_mod, only: mpp_error, FATAL, stdlog, input_nml_file + 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_processmodel + public :: read_namelist_fv_processmodel + contains + + subroutine fv_init_processmodel(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 + zrad = Abs(zm - zc)/zc + + 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_processmodel + + subroutine read_namelist_fv_processmodel(nml_filename) + + character(*), intent(IN) :: nml_filename + integer :: ierr, f_unit, unit, ios + namelist /fv_processmodel_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_processmodel_nml,iostat=ios) + ierr = check_nml_error(ios,'fv_processmodel_nml') + write(unit, nml=fv_processmodel_nml) + + + end subroutine read_namelist_fv_processmodel + + 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_processmodel_mod diff --git a/model/fv_arrays.F90 b/model/fv_arrays.F90 index 15f4daa8c..ab9d0ccd7 100644 --- a/model/fv_arrays.F90 +++ b/model/fv_arrays.F90 @@ -918,6 +918,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 50d60157f..3813433d5 100644 --- a/model/fv_control.F90 +++ b/model/fv_control.F90 @@ -141,6 +141,9 @@ 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 +#ifdef GFS_PHYS + use fv_processmodel_mod, only: read_namelist_fv_processmodel +#endif 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 @@ -403,6 +406,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 @@ -582,7 +587,12 @@ 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 + +#ifdef GFS_PHYS + call read_namelist_fv_processmodel(Atm(this_grid)%nml_filename) +#endif 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) @@ -987,6 +997,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 @@ -1093,7 +1105,7 @@ subroutine read_namelist_fv_core_nml(Atm) sa3dtke_dyco, & 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 063589282..10124fb8e 100644 --- a/model/fv_dynamics.F90 +++ b/model/fv_dynamics.F90 @@ -1430,6 +1430,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 @@ -1443,6 +1445,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/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_restart.F90 b/tools/fv_restart.F90 index d87d93582..6ff582c20 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_processmodel_mod, only: fv_init_processmodel +#endif implicit none private @@ -430,6 +433,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_processmodel" + call fv_init_processmodel(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, & @@ -443,6 +462,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") @@ -645,6 +665,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 @@ -657,11 +678,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