Addition of Water Vapor Age Distribution Moments Tracers

.gitignore

@@ -110,4 +110,4 @@ ENV/
 .spyproject
 
 # mkdocs documentation
-/site
+/site

bin/mkmf.template.nyu

@@ -3,7 +3,7 @@
 # mkmf -t template.ifc -c"-Duse_libMPI -Duse_netCDF" path_names /usr/local/include
 CPPFLAGS = -I/usr/local/include 
 CPPFLAGS += -I${NETCDF_INC}
-FFLAGS = $(CPPFLAGS) -fpp -stack_temps -safe_cray_ptr -ftz -i-dynamic -assume byterecl -i4 -r8
+FFLAGS = $(CPPFLAGS) -fpp -stack_temps -safe_cray_ptr -ftz -i-dynamic -assume byterecl -i4 -r8 -flexiblas
 FFLAGS += -O2
 FFLAGS += ${DEBUG}
 FFLAGS += -I${NETCDF_INC} -I${MPI_INC}

exp/test_cases/MiMA/MiMA_test_case.py

@@ -19,8 +19,6 @@
 # is used to load the correct compilers.  The env file is always loaded from
 # $GFDL_BASE and not the checked out git repo.
 
-cb.compile()  # compile the source code to working directory $GFDL_WORK/codebase
-
 # create an Experiment object to handle the configuration of model parameters
 # and output diagnostics
 exp = Experiment('mima_test_experiment', codebase=cb)
@@ -174,6 +172,9 @@
 })
 #Lets do a run!
 if __name__=="__main__":
+
+    cb.compile()  # compile the source code to working directory $GFDL_WORK/codebase
+
     exp.run(1, use_restart=False, num_cores=NCORES)
     for i in range(2,121):
         exp.run(i, num_cores=NCORES)

exp/test_cases/ape_aquaplanet/socrates_ape_aquaplanet_T42.py

@@ -161,7 +161,7 @@
         'evaporation':True,  
         'depth': 10.0,                          #Depth of mixed layer used
         'albedo_value': 0.38,                   #Albedo value used      
-        'do_ape_sst': True,
+        'do_ape_sst' : True,
     },
 
     'qe_moist_convection_nml': {

exp/test_cases/bucket_hydrology/bucket_model_test_case.py

@@ -19,8 +19,6 @@
 # is used to load the correct compilers.  The env file is always loaded from
 # $GFDL_BASE and not the checked out git repo.
 
-cb.compile()  # compile the source code to working directory $GFDL_WORK/codebase
-
 # create an Experiment object to handle the configuration of model parameters
 # and output diagnostics
 exp = Experiment('bucket_test_experiment', codebase=cb)
@@ -179,6 +177,8 @@
 
 #Lets do a run!
 if __name__=="__main__":
+    cb.compile()  # compile the source code to working directory $GFDL_WORK/codebase
+
     exp.run(1, use_restart=False, num_cores=NCORES)
     for i in range(2,121):
         exp.run(i, num_cores=NCORES)

exp/test_cases/frierson/frierson_test_case.py

@@ -19,8 +19,6 @@
 # is used to load the correct compilers.  The env file is always loaded from
 # $GFDL_BASE and not the checked out git repo.
 
-cb.compile()  # compile the source code to working directory $GFDL_WORK/codebase
-
 # create an Experiment object to handle the configuration of model parameters
 # and output diagnostics
 exp = Experiment('frierson_test_experiment', codebase=cb)
@@ -173,6 +171,8 @@
 
 #Lets do a run!
 if __name__=="__main__":
+    cb.compile()  # compile the source code to working directory $GFDL_WORK/codebase
+
     exp.run(1, use_restart=False, num_cores=NCORES)
     for i in range(2,121):
         exp.run(i, num_cores=NCORES)

exp/test_cases/realistic_continents/realistic_continents_fixed_sst_test_case.py

@@ -21,8 +21,6 @@
 # is used to load the correct compilers.  The env file is always loaded from
 # $GFDL_BASE and not the checked out git repo.
 
-cb.compile()  # compile the source code to working directory $GFDL_WORK/codebase
-
 # create an Experiment object to handle the configuration of model parameters
 # and output diagnostics
 exp = Experiment('realistic_continents_fixed_sst_test_experiment', codebase=cb)
@@ -72,6 +70,8 @@
 
 #Lets do a run!
 if __name__=="__main__":
+    cb.compile()  # compile the source code to working directory $GFDL_WORK/codebase
+
     exp.run(1, use_restart=False, num_cores=NCORES)
     for i in range(2,121):
         exp.run(i, num_cores=NCORES)

exp/test_cases/realistic_continents/realistic_continents_fixed_sst_wv_age_test_case.py

@@ -0,0 +1,90 @@
+import os
+
+import numpy as np
+
+import f90nml
+
+from isca import IscaCodeBase, DiagTable, Experiment, Namelist, GFDL_BASE
+
+
+NCORES = 4
+base_dir = os.path.dirname(os.path.realpath(__file__))
+# a CodeBase can be a directory on the computer,
+# useful for iterative development
+cb = IscaCodeBase.from_directory(GFDL_BASE)
+
+# or it can point to a specific git repo and commit id.
+# This method should ensure future, independent, reproducibility of results.
+# cb = DryCodeBase.from_repo(repo='https://github.com/isca/isca', commit='isca1.1')
+
+# compilation depends on computer specific settings.  The $GFDL_ENV
+# environment variable is used to determine which `$GFDL_BASE/src/extra/env` file
+# is used to load the correct compilers.  The env file is always loaded from
+# $GFDL_BASE and not the checked out git repo.
+
+# create an Experiment object to handle the configuration of model parameters
+# and output diagnostics
+
+# Select field table for age moment tracers
+n_moments = 2
+field_table_name = f"field_table_age_{n_moments}"
+
+
+exp = Experiment('realistic_continents_fixed_sst_wv_age_test_experiment', codebase=cb,field_table_name = field_table_name)
+
+#Add any input files that are necessary for a particular experiment.
+exp.inputfiles = [os.path.join(GFDL_BASE,'input/land_masks/era_land_t42.nc'),os.path.join(GFDL_BASE,'input/rrtm_input_files/ozone_1990.nc'),
+                      os.path.join(base_dir,'input/sst_clim_amip.nc'), os.path.join(base_dir,'input/siconc_clim_amip.nc')]
+
+#Tell model how to write diagnostics
+diag = DiagTable()
+diag.add_file('atmos_monthly', 30, 'days', time_units='days')
+
+
+#Tell model which diagnostics to write
+diag.add_field('dynamics', 'ps', time_avg=True)
+diag.add_field('dynamics', 'bk')
+diag.add_field('dynamics', 'pk')
+diag.add_field('dynamics', 'zsurf') #need at least ps, pk, bk and zsurf to do vertical interpolation onto plevels from sigma
+diag.add_field('atmosphere', 'precipitation', time_avg=True)
+diag.add_field('mixed_layer', 't_surf', time_avg=True)
+diag.add_field('dynamics', 'sphum', time_avg=True)
+diag.add_field('dynamics', 'ucomp', time_avg=True)
+diag.add_field('dynamics', 'vcomp', time_avg=True)
+diag.add_field('dynamics', 'temp', time_avg=True)
+diag.add_field('dynamics', 'vor', time_avg=True)
+diag.add_field('dynamics', 'div', time_avg=True)
+
+# WV age moments diagnostics
+for ind in range(n_moments):
+    name = f"sphum_age_{ind+1}"
+    diag.add_field('dynamics', name, time_avg=True)
+
+exp.diag_table = diag
+
+
+#Empty the run directory ready to run
+exp.clear_rundir()
+
+#Define values for the 'core' namelist
+namelist_name = os.path.join(GFDL_BASE, 'exp/test_cases/realistic_continents/namelist_basefile.nml')
+nml = f90nml.read(namelist_name)
+exp.namelist = nml
+
+exp.update_namelist({
+    'mixed_layer_nml': {  
+        'do_qflux' : False, #Don't use the prescribed analytical formula for q-fluxes
+        'do_read_sst' : True, #Read in sst values from input file
+        'do_sc_sst' : True, #Do specified ssts (need both to be true)
+        'sst_file' : 'sst_clim_amip', #Set name of sst input file
+        'specify_sst_over_ocean_only' : True, #Make sure sst only specified in regions of ocean.
+    }
+})
+
+#Lets do a run!
+if __name__=="__main__":
+    cb.compile()  # compile the source code to working directory $GFDL_WORK/codebase
+
+    exp.run(1, use_restart=False, num_cores=NCORES)
+    for i in range(2,121):
+        exp.run(i, num_cores=NCORES)

exp/test_cases/realistic_continents/realistic_continents_variable_qflux_test_case.py

@@ -21,8 +21,6 @@
 # is used to load the correct compilers.  The env file is always loaded from
 # $GFDL_BASE and not the checked out git repo.
 
-cb.compile()  # compile the source code to working directory $GFDL_WORK/codebase
-
 # create an Experiment object to handle the configuration of model parameters
 # and output diagnostics
 exp = Experiment('realistic_continents_qflux_test_experiment', codebase=cb)

exp/test_cases/trip_test/trip_test_functions.py

@@ -193,7 +193,7 @@ def conduct_comparison_on_test_case(base_commit, later_commit, test_case_name, r
         else:
             cb = IscaCodeBase(repo=repo_to_use, commit=s)
         try:
-            cb.compile()
+            cb.compile(debug = False)
             exp = Experiment(exp_name, codebase=cb)
             exp.namelist = nml_use.copy()
             exp.diag_table = diag_use

src/atmos_column/column.F90

@@ -58,7 +58,7 @@ module column_mod
 !real, allocatable, dimension(:) :: lat_boundaries_global, lon_boundaries_global
 
 real :: virtual_factor, dt_real
-integer :: pe, npes, num_tracers, nhum, step_number
+integer :: pe, npes, num_tracers, nhum,nhum_age,nhum_sink, step_number
 integer :: is, ie, js, je
 integer :: previous, current, future 
 
@@ -109,14 +109,16 @@ module column_mod
 
 contains 
 
-subroutine column_init(Time, Time_step_in, tracer_attributes, dry_model_out, nhum_out)
+subroutine column_init(Time, Time_step_in, tracer_attributes, dry_model_out, nhum_out,nhum_out_age,nhum_out_sink)
 
     type(time_type), intent(in) :: Time, Time_step_in
     type(tracer_type), intent(inout), dimension(:) :: tracer_attributes
     logical, intent(out) :: dry_model_out
     integer, intent(out) :: nhum_out
+    integer, intent(out) :: nhum_out_age
+    integer, intent(out) :: nhum_out_sink    
 
-    integer :: unit, ierr, io, ntr, nsphum, nmix_rat, seconds, days 
+    integer :: unit, ierr, io, ntr, nsphum,nsphum_age,nsphum_sink, nmix_rat, seconds, days 
     !real :: del_lon, del_lat !!! NTL START HERE
     !real :: longitude_origin_local = 0.0
 
@@ -172,11 +174,15 @@ subroutine column_init(Time, Time_step_in, tracer_attributes, dry_model_out, nhu
 
     enddo
     nsphum   = get_tracer_index(MODEL_ATMOS, 'sphum')
+    nsphum_age   = get_tracer_index(MODEL_ATMOS, 'sphum_age')
+    nsphum_sink   = get_tracer_index(MODEL_ATMOS, 'sphum_sink')
     nmix_rat = get_tracer_index(MODEL_ATMOS, 'mix_rat')
 
     if(nsphum == NO_TRACER) then
         if(nmix_rat == NO_TRACER) then
           nhum = 0
+          nhum_age = 0
+          nhum_sink = 0
           dry_model = .true.
         else
           nhum = nmix_rat
@@ -185,13 +191,17 @@ subroutine column_init(Time, Time_step_in, tracer_attributes, dry_model_out, nhu
     else
         if(nmix_rat == NO_TRACER) then
           nhum = nsphum
+          nhum_age = nsphum_age
+          nhum_sink = nsphum_sink
           dry_model = .false.
         else
           call error_mesg('column_init','sphum and mix_rat cannot both be specified as tracers at the same time', FATAL)
         endif
     endif
     dry_model_out = dry_model
     nhum_out = nhum
+    nhum_out_age = nhum_age
+    nhum_out_sink = nhum_sink
 
 
     call read_restart_or_do_coldstart(tracer_attributes)

src/atmos_param/qe_moist_convection/qe_moist_convection.F90

@@ -336,6 +336,7 @@ subroutine SBM_convection_scheme(dt, Tin, qin, p_full, p_half, rain,  snow, &
              ! Calculate the precipitation rate Pq 
              call Pq_calculation(kLZB, k_surface, qref_parcel, qin(i,j,:),    &
                   p_half(i,j,:), deltaq_parcel, Pq_parcel, dt)
+
              ! Calculate the humidity change that would be necessary
              ! to balance temperature change by latent heat release 
              call Pt_calculation(kLZB, k_surface, Tref_parcel, Tin(i,j,:),    &

src/atmos_solo/atmos_model.F90

@@ -115,15 +115,11 @@ program atmos_model
  call constants_init
  call fms_init ( )
  call atmos_model_init
-
 !   ------ atmosphere integration loop -------
 
     call mpp_clock_begin (id_loop)
-
     do na = 1, num_atmos_calls
-
        call atmosphere (Time)
-
        Time = Time + Time_step_atmos
 
        if(modulo(na,memuse_interval) == 0 .and. print_memuse) then
@@ -132,7 +128,6 @@ program atmos_model
        endif
 
     enddo
-
     call mpp_clock_end (id_loop)
 
 !   ------ end of atmospheric time step loop -----

src/atmos_spectral/driver/solo/age_moments.F90

@@ -0,0 +1,29 @@
+module age_moments_mod
+    implicit none
+    public :: get_age_moments
+contains
+
+subroutine get_age_moments(nsphum, nsphum_age, previous, grid_tracers, sink, dt_tracers)
+    implicit none
+
+    integer,                    intent(in)    :: nsphum
+    integer,                    intent(in)    :: nsphum_age
+    integer,                    intent(in)    :: previous
+    real, dimension(:,:,:,:,:), intent(in)    :: grid_tracers
+    real, dimension(:,:,:),   intent(in)      :: sink
+    real, dimension(:,:,:,:),   intent(inout)   :: dt_tracers
+
+    real :: eps_blowup
+    integer :: i
+
+    eps_blowup = 1e-10
+
+    ! Calculates the RHS of age-moment evolution equation
+    do i = 2, nsphum_age
+        dt_tracers(:,:,:,i)  = dt_tracers(:,:,:,i) +  (i-1) * grid_tracers(:,:,:,previous,i-1) + sink * (grid_tracers(:,:,:,previous,i)/(eps_blowup+grid_tracers(:,:,:,previous,nsphum)))
+    end do
+
+end subroutine get_age_moments
+
+end module age_moments_mod
+