diff --git a/.gitignore b/.gitignore index 0b3138728d..e534738ed7 100644 --- a/.gitignore +++ b/.gitignore @@ -6,16 +6,3 @@ html MOM6 build/ deps/ -#.testing/*/available_diags.* -#.testing/*/CPU_stats -#.testing/*/chksum_diag -#.testing/*/exitcode -#.testing/*/logfile.*.out -#.testing/*/MOM_parameter_doc.* -#.testing/*/ocean_geometry.nc -#.testing/*/ocean.stats -#.testing/*/ocean.stats.nc -#.testing/*/RESTART/ -#.testing/*/time_stamp.out -#.testing/*/Vertical_coordinate.nc -#.testing/*/GOLD_IC.nc diff --git a/.testing/.gitignore b/.testing/.gitignore index f119a40591..a096823fcd 100644 --- a/.testing/.gitignore +++ b/.testing/.gitignore @@ -5,9 +5,11 @@ exitcode logfile.*.out MOM_parameter_doc.* ocean_geometry.nc -ocean.stats -ocean.stats.nc +ocean.stats* RESTART/ time_stamp.out Vertical_coordinate.nc GOLD_IC.nc +debug.out +chksum_diag.* +config.mk diff --git a/.testing/Makefile b/.testing/Makefile index 1dee0e2100..d3093bf523 100644 --- a/.testing/Makefile +++ b/.testing/Makefile @@ -37,7 +37,7 @@ MKMF_TEMPLATE ?= $(DEPS)/mkmf/templates/ncrc-gnu.mk # Executables BUILDS = symmetric asymmetric repro -CONFIGS := $(foreach n,$(shell seq 0 3),tc$(n)) +CONFIGS := $(wildcard tc*) TESTS = grids layouts restarts repros nans dims # The following variables are configured by Travis: diff --git a/.testing/tc0/diag_table b/.testing/tc0/diag_table index 1527de166b..091dc46933 100644 --- a/.testing/tc0/diag_table +++ b/.testing/tc0/diag_table @@ -1,2 +1,2 @@ -"Unit tests" +"MOM test configuration 0" 1 1 1 0 0 0 diff --git a/.testing/tc1.a/MOM_input b/.testing/tc1.a/MOM_input new file mode 120000 index 0000000000..dca928737e --- /dev/null +++ b/.testing/tc1.a/MOM_input @@ -0,0 +1 @@ +../tc1/MOM_input \ No newline at end of file diff --git a/.testing/tc1.a/MOM_override b/.testing/tc1.a/MOM_override new file mode 100644 index 0000000000..e69de29bb2 diff --git a/.testing/tc1.a/MOM_tc_variant b/.testing/tc1.a/MOM_tc_variant new file mode 100644 index 0000000000..8032901a82 --- /dev/null +++ b/.testing/tc1.a/MOM_tc_variant @@ -0,0 +1 @@ +#override SPLIT=False diff --git a/.testing/tc1.a/diag_table b/.testing/tc1.a/diag_table new file mode 120000 index 0000000000..bf2ad677b6 --- /dev/null +++ b/.testing/tc1.a/diag_table @@ -0,0 +1 @@ +../tc1/diag_table \ No newline at end of file diff --git a/.testing/tc1.a/input.nml b/.testing/tc1.a/input.nml new file mode 100644 index 0000000000..3c7dcf7bea --- /dev/null +++ b/.testing/tc1.a/input.nml @@ -0,0 +1,20 @@ +&mom_input_nml + output_directory = './' + input_filename = 'n' + restart_input_dir = 'INPUT/' + restart_output_dir = 'RESTART/' + parameter_filename = + 'MOM_input', + 'MOM_tc_variant', + 'MOM_override', +/ + +&diag_manager_nml +/ + +&fms_nml + clock_grain = 'ROUTINE' + clock_flags = 'SYNC' + domains_stack_size = 955296 + stack_size = 0 +/ diff --git a/.testing/tc1/diag_table b/.testing/tc1/diag_table index 19d6a32e1e..220d65d34f 100644 --- a/.testing/tc1/diag_table +++ b/.testing/tc1/diag_table @@ -1,86 +1,10 @@ -"MOM benchmark Experiment" +"MOM test configuration 1" 1 1 1 0 0 0 "prog", 1,"days",1,"days","time", -#"ave_prog", 5,"days",1,"days","Time",365,"days" -#"cont", 5,"days",1,"days","Time",365,"days" - -#This is the field section of the diag_table. # Prognostic Ocean fields: -#========================= - "ocean_model","u","u","prog","all",.false.,"none",2 "ocean_model","v","v","prog","all",.false.,"none",2 "ocean_model","h","h","prog","all",.false.,"none",1 "ocean_model","e","e","prog","all",.false.,"none",2 "ocean_model","temp","temp","prog","all",.false.,"none",2 -#"ocean_model","salt","salt","prog","all",.false.,"none",2 - -#"ocean_model","u","u","ave_prog_%4yr_%3dy","all",.true.,"none",2 -#"ocean_model","v","v","ave_prog_%4yr_%3dy","all",.true.,"none",2 -#"ocean_model","h","h","ave_prog_%4yr_%3dy","all",.true.,"none",1 -#"ocean_model","e","e","ave_prog_%4yr_%3dy","all",.true.,"none",2 - -# Auxilary Tracers: -#================== -#"ocean_model","vintage","vintage","prog_%4yr_%3dy","all",.false.,"none",2 -#"ocean_model","age","age","prog_%4yr_%3dy","all",.false.,"none",2 - -# Continuity Equation Terms: -#=========================== -#"ocean_model","dhdt","dhdt","cont_%4yr_%3dy","all",.true.,"none",2 -#"ocean_model","wd","wd","cont_%4yr_%3dy","all",.true.,"none",2 -#"ocean_model","uh","uh","cont_%4yr_%3dy","all",.true.,"none",2 -#"ocean_model","vh","vh","cont_%4yr_%3dy","all",.true.,"none",2 -#"ocean_model","h_rho","h_rho","cont_%4yr_%3dy","all",.true.,"none",2 -#"ocean_model","uh_rho","uh_rho","cont_%4yr_%3dy","all",.true.,"none",2 -#"ocean_model","vh_rho","vh_rho","cont_%4yr_%3dy","all",.true.,"none",2 -#"ocean_model","uhGM_rho","uhGM_rho","cont_%4yr_%3dy","all",.true.,"none",2 -#"ocean_model","vhGM_rho","vhGM_rho","cont_%4yr_%3dy","all",.true.,"none",2 - -# -# Tracer Fluxes: -#================== -#"ocean_model","T_adx", "T_adx", "ave_prog_%4yr_%3dy","all",.true.,"none",2 -#"ocean_model","T_ady", "T_ady", "ave_prog_%4yr_%3dy","all",.true.,"none",2 -#"ocean_model","T_diffx","T_diffx","ave_prog_%4yr_%3dy","all",.true.,"none",2 -#"ocean_model","T_diffy","T_diffy","ave_prog_%4yr_%3dy","all",.true.,"none",2 -#"ocean_model","S_adx", "S_adx", "ave_prog_%4yr_%3dy","all",.true.,"none",2 -#"ocean_model","S_ady", "S_ady", "ave_prog_%4yr_%3dy","all",.true.,"none",2 -#"ocean_model","S_diffx","S_diffx","ave_prog_%4yr_%3dy","all",.true.,"none",2 -#"ocean_model","S_diffy","S_diffy","ave_prog_%4yr_%3dy","all",.true.,"none",2 - -#============================================================================================= -# -#===- This file can be used with diag_manager/v2.0a (or higher) ==== -# -# -# FORMATS FOR FILE ENTRIES (not all input values are used) -# ------------------------ -# -#"file_name", output_freq, "output_units", format, "time_units", "time_long_name", ... -# (opt) new_file_frequecy, (opt) "new_file_freq_units", "new_file_start_date" -# -# -#output_freq: > 0 output frequency in "output_units" -# = 0 output frequency every time step -# =-1 output frequency at end of run -# -#output_units = units used for output frequency -# (years, months, days, minutes, hours, seconds) -# -#time_units = units used to label the time axis -# (days, minutes, hours, seconds) -# -# -# FORMAT FOR FIELD ENTRIES (not all input values are used) -# ------------------------ -# -#"module_name", "field_name", "output_name", "file_name" "time_sampling", time_avg, "other_opts", packing -# -#time_avg = .true. or .false. -# -#packing = 1 double precision -# = 2 float -# = 4 packed 16-bit integers -# = 8 packed 1-byte (not tested?) diff --git a/.testing/tc2/diag_table b/.testing/tc2/diag_table index 19d6a32e1e..941b9c0c15 100644 --- a/.testing/tc2/diag_table +++ b/.testing/tc2/diag_table @@ -1,86 +1,2 @@ -"MOM benchmark Experiment" +"MOM test configuration 2" 1 1 1 0 0 0 -"prog", 1,"days",1,"days","time", -#"ave_prog", 5,"days",1,"days","Time",365,"days" -#"cont", 5,"days",1,"days","Time",365,"days" - -#This is the field section of the diag_table. - -# Prognostic Ocean fields: -#========================= - -"ocean_model","u","u","prog","all",.false.,"none",2 -"ocean_model","v","v","prog","all",.false.,"none",2 -"ocean_model","h","h","prog","all",.false.,"none",1 -"ocean_model","e","e","prog","all",.false.,"none",2 -"ocean_model","temp","temp","prog","all",.false.,"none",2 -#"ocean_model","salt","salt","prog","all",.false.,"none",2 - -#"ocean_model","u","u","ave_prog_%4yr_%3dy","all",.true.,"none",2 -#"ocean_model","v","v","ave_prog_%4yr_%3dy","all",.true.,"none",2 -#"ocean_model","h","h","ave_prog_%4yr_%3dy","all",.true.,"none",1 -#"ocean_model","e","e","ave_prog_%4yr_%3dy","all",.true.,"none",2 - -# Auxilary Tracers: -#================== -#"ocean_model","vintage","vintage","prog_%4yr_%3dy","all",.false.,"none",2 -#"ocean_model","age","age","prog_%4yr_%3dy","all",.false.,"none",2 - -# Continuity Equation Terms: -#=========================== -#"ocean_model","dhdt","dhdt","cont_%4yr_%3dy","all",.true.,"none",2 -#"ocean_model","wd","wd","cont_%4yr_%3dy","all",.true.,"none",2 -#"ocean_model","uh","uh","cont_%4yr_%3dy","all",.true.,"none",2 -#"ocean_model","vh","vh","cont_%4yr_%3dy","all",.true.,"none",2 -#"ocean_model","h_rho","h_rho","cont_%4yr_%3dy","all",.true.,"none",2 -#"ocean_model","uh_rho","uh_rho","cont_%4yr_%3dy","all",.true.,"none",2 -#"ocean_model","vh_rho","vh_rho","cont_%4yr_%3dy","all",.true.,"none",2 -#"ocean_model","uhGM_rho","uhGM_rho","cont_%4yr_%3dy","all",.true.,"none",2 -#"ocean_model","vhGM_rho","vhGM_rho","cont_%4yr_%3dy","all",.true.,"none",2 - -# -# Tracer Fluxes: -#================== -#"ocean_model","T_adx", "T_adx", "ave_prog_%4yr_%3dy","all",.true.,"none",2 -#"ocean_model","T_ady", "T_ady", "ave_prog_%4yr_%3dy","all",.true.,"none",2 -#"ocean_model","T_diffx","T_diffx","ave_prog_%4yr_%3dy","all",.true.,"none",2 -#"ocean_model","T_diffy","T_diffy","ave_prog_%4yr_%3dy","all",.true.,"none",2 -#"ocean_model","S_adx", "S_adx", "ave_prog_%4yr_%3dy","all",.true.,"none",2 -#"ocean_model","S_ady", "S_ady", "ave_prog_%4yr_%3dy","all",.true.,"none",2 -#"ocean_model","S_diffx","S_diffx","ave_prog_%4yr_%3dy","all",.true.,"none",2 -#"ocean_model","S_diffy","S_diffy","ave_prog_%4yr_%3dy","all",.true.,"none",2 - -#============================================================================================= -# -#===- This file can be used with diag_manager/v2.0a (or higher) ==== -# -# -# FORMATS FOR FILE ENTRIES (not all input values are used) -# ------------------------ -# -#"file_name", output_freq, "output_units", format, "time_units", "time_long_name", ... -# (opt) new_file_frequecy, (opt) "new_file_freq_units", "new_file_start_date" -# -# -#output_freq: > 0 output frequency in "output_units" -# = 0 output frequency every time step -# =-1 output frequency at end of run -# -#output_units = units used for output frequency -# (years, months, days, minutes, hours, seconds) -# -#time_units = units used to label the time axis -# (days, minutes, hours, seconds) -# -# -# FORMAT FOR FIELD ENTRIES (not all input values are used) -# ------------------------ -# -#"module_name", "field_name", "output_name", "file_name" "time_sampling", time_avg, "other_opts", packing -# -#time_avg = .true. or .false. -# -#packing = 1 double precision -# = 2 float -# = 4 packed 16-bit integers -# = 8 packed 1-byte (not tested?) diff --git a/.testing/tc3/MOM_input b/.testing/tc3/MOM_input index 1689ef993e..430ce24b61 100644 --- a/.testing/tc3/MOM_input +++ b/.testing/tc3/MOM_input @@ -35,11 +35,11 @@ NJHALO = 4 ! default = 2 ! y-direction. With STATIC_MEMORY_ this is set as NJHALO_ ! in MOM_memory.h at compile time; without STATIC_MEMORY_ ! the default is NJHALO_ in MOM_memory.h (if defined) or 2. -NIGLOBAL = 25 ! +NIGLOBAL = 13 ! ! The total number of thickness grid points in the ! x-direction in the physical domain. With STATIC_MEMORY_ ! this is set in MOM_memory.h at compile time. -NJGLOBAL = 25 ! +NJGLOBAL = 13 ! ! The total number of thickness grid points in the ! y-direction in the physical domain. With STATIC_MEMORY_ ! this is set in MOM_memory.h at compile time. diff --git a/.testing/tc3/diag_table b/.testing/tc3/diag_table index e31244cbd4..64043b6e0d 100644 --- a/.testing/tc3/diag_table +++ b/.testing/tc3/diag_table @@ -1,207 +1,2 @@ -"MOM Experiment" +"MOM test configuration 3" 1 1 1 0 0 0 -"prog", 2,"minutes",1,"days","Time", -#"ave_prog", 1,"hours",1,"days","Time", -#"cont", 1,"hours",1,"days","Time", -#"trac", 5,"days",1,"days","Time", -#"mom", 5,"days",1,"days","Time", -#"bt_mom", 5,"days",1,"days","Time", -#"visc", 5,"days",1,"days","Time", -#"energy", 5,"days",1,"days","Time", -#"ML_TKE", 5,"days",1,"days","Time", -#"forcing", 5,"days",1,"days","Time", - -#This is the field section of the diag_table. - -# Prognostic Ocean fields: -#========================= - -"ocean_model","u","u","prog","all",.false.,"none",2 -"ocean_model","v","v","prog","all",.false.,"none",2 -"ocean_model","h","h","prog","all",.false.,"none",1 -"ocean_model","e","e","prog","all",.false.,"none",2 -#"ocean_model","SSH","SSH","prog","all",.false.,"none",2 -#"ocean_model","temp","temp","prog","all",.false.,"none",2 -#"ocean_model","salt","salt","prog","all",.false.,"none",2 -#"ocean_model","Rml","Rml","prog","all",.false.,"none",2 -#"ocean_model","tr_D1","tr1","prog","all",.false.,"none",2 - -#"ocean_model","RV","RV","prog","all",.false.,"none",2 -#"ocean_model","PV","PV","prog","all",.false.,"none",2 -#"ocean_model","e_D","e_D","prog","all",.false.,"none",2 - -#"ocean_model","u","u","ave_prog","all",.true.,"none",2 -#"ocean_model","v","v","ave_prog","all",.true.,"none",2 -#"ocean_model","h","h","ave_prog","all",.true.,"none",1 -#"ocean_model","e","e","ave_prog","all",.true.,"none",2 -#"ocean_model","temp","temp","ave_prog","all",.true.,"none",2 -#"ocean_model","salt","salt","ave_prog","all",.true.,"none",2 -#"ocean_model","Rml","Rml","ave_prog","all",.true.,"none",2 - -# Auxilary Tracers: -#================== -#"ocean_model","vintage","vintage","prog","all",.false.,"none",2 -#"ocean_model","age","age","prog","all",.false.,"none",2 - -# Tracers: -#========= -#"ocean_model","tr_D1","tr1","trac","all",.false.,"none",2 -#"ocean_model","tr_D2","tr2","trac","all",.false.,"none",2 -#"ocean_model","tr_D3","tr3","trac","all",.false.,"none",2 -#"ocean_model","tr_D4","tr4","trac","all",.false.,"none",2 -#"ocean_model","tr_D5","tr5","trac","all",.false.,"none",2 -#"ocean_model","tr_D6","tr6","trac","all",.false.,"none",2 -#"ocean_model","tr_D7","tr7","trac","all",.false.,"none",2 -#"ocean_model","tr_D8","tr8","trac","all",.false.,"none",2 -#"ocean_model","tr_D9","tr9","trac","all",.false.,"none",2 -#"ocean_model","tr_D10","tr10","trac","all",.false.,"none",2 -#"ocean_model","tr_D11","tr11","trac","all",.false.,"none",2 - -# Continuity Equation Terms: -#=========================== -#"ocean_model","dhdt","dhdt","cont","all",.true.,"none",2 -#"ocean_model","wd","wd","cont","all",.true.,"none",2 -#"ocean_model","uh","uh","cont","all",.true.,"none",2 -#"ocean_model","vh","vh","cont","all",.true.,"none",2 -#"ocean_model","uhGM","uhGM","cont","all",.true.,"none",2 -#"ocean_model","vhGM","vhGM","cont","all",.true.,"none",2 -#"ocean_model","uhbt","uhbt","cont","all",.true.,"none",2 -#"ocean_model","vhbt","vhbt","cont","all",.true.,"none",2 - -# Continuity Equation Terms In Pure Potential Density Coordiantes: -#================================================================= -#"ocean_model","h_rho","h_rho","cont","all",.true.,"none",2 -#"ocean_model","uh_rho","uh_rho","cont","all",.true.,"none",2 -#"ocean_model","vh_rho","vh_rho","cont","all",.true.,"none",2 -#"ocean_model","uhGM_rho","uhGM_rho","cont","all",.true.,"none",2 -#"ocean_model","vhGM_rho","vhGM_rho","cont","all",.true.,"none",2 - -# -# Tracer Fluxes: -#================== -#"ocean_model","T_adx", "T_adx", "ave_prog","all",.true.,"none",2 -#"ocean_model","T_ady", "T_ady", "ave_prog","all",.true.,"none",2 -#"ocean_model","T_diffx","T_diffx","ave_prog","all",.true.,"none",2 -#"ocean_model","T_diffy","T_diffy","ave_prog","all",.true.,"none",2 -#"ocean_model","S_adx", "S_adx", "ave_prog","all",.true.,"none",2 -#"ocean_model","S_ady", "S_ady", "ave_prog","all",.true.,"none",2 -#"ocean_model","S_diffx","S_diffx","ave_prog","all",.true.,"none",2 -#"ocean_model","S_diffy","S_diffy","ave_prog","all",.true.,"none",2 - - -# Momentum Balance Terms: -#======================= -#"ocean_model","dudt","dudt","mom","all",.true.,"none",2 -#"ocean_model","dvdt","dvdt","mom","all",.true.,"none",2 -#"ocean_model","CAu","CAu","mom","all",.true.,"none",2 -#"ocean_model","CAv","CAv","mom","all",.true.,"none",2 -#"ocean_model","PFu","PFu","mom","all",.true.,"none",2 -#"ocean_model","PFv","PFv","mom","all",.true.,"none",2 -#"ocean_model","du_dt_visc","du_dt_visc","mom","all",.true.,"none",2 -#"ocean_model","dv_dt_visc","dv_dt_visc","mom","all",.true.,"none",2 -#"ocean_model","diffu","diffu","mom","all",.true.,"none",2 -#"ocean_model","diffv","diffv","mom","all",.true.,"none",2 -#"ocean_model","dudt_dia","dudt_dia","mom","all",.true.,"none",2 -#"ocean_model","dvdt_dia","dvdt_dia","mom","all",.true.,"none",2 -# Subterms that should not be added to a closed budget. -#"ocean_model","gKEu","gKEu","mom","all",.true.,"none",2 -#"ocean_model","gKEv","gKEv","mom","all",.true.,"none",2 -#"ocean_model","rvxu","rvxu","mom","all",.true.,"none",2 -#"ocean_model","rvxv","rvxv","mom","all",.true.,"none",2 -#"ocean_model","PFu_bc","PFu_bc","mom","all",.true.,"none",2 -#"ocean_model","PFv_bc","PFv_bc","mom","all",.true.,"none",2 - -# Barotropic Momentum Balance Terms: -# (only available with split time stepping.) -#=========================================== -#"ocean_model","PFuBT","PFuBT","bt_mom","all",.true.,"none",2 -#"ocean_model","PFvBT","PFvBT","bt_mom","all",.true.,"none",2 -#"ocean_model","CoruBT","CoruBT","bt_mom","all",.true.,"none",2 -#"ocean_model","CorvBT","CorvBT","bt_mom","all",.true.,"none",2 -#"ocean_model","ubtforce","ubtforce","bt_mom","all",.true.,"none",2 -#"ocean_model","vbtforce","vbtforce","bt_mom","all",.true.,"none",2 -#"ocean_model","u_accel_bt","u_accel_bt","bt_mom","all",.true.,"none",2 -#"ocean_model","v_accel_bt","v_accel_bt","bt_mom","all",.true.,"none",2 -# -# Viscosities and diffusivities: -#=============================== -#"ocean_model","Kd_effective","Kd_effective","visc","all",.true.,"none",2 -#"ocean_model","Ahh","Ahh","visc","all",.true.,"none",2 -#"ocean_model","Ahq","Ahq","visc","all",.true.,"none",2 -#"ocean_model","Khh","Khh","visc","all",.true.,"none",2 -#"ocean_model","Khq","Khq","visc","all",.true.,"none",2 -#"ocean_model","bbl_thick_u","bbl_thick_u","visc","all",.true.,"none",2 -#"ocean_model","kv_bbl_u","kv_bbl_u","visc","all",.true.,"none",2 -#"ocean_model","bbl_thick_v","bbl_thick_v","visc","all",.true.,"none",2 -#"ocean_model","kv_bbl_v","kv_bbl_v","visc","all",.true.,"none",2 -#"ocean_model","av_visc","av_visc","visc","all",.true.,"none",2 -#"ocean_model","au_visc","au_visc","visc","all",.true.,"none",2 -# -# Kinetic Energy Balance Terms: -#============================= -#"ocean_model","KE","KE","energy","all",.true.,"none",2 -#"ocean_model","dKE_dt","dKE_dt","energy","all",.true.,"none",2 -#"ocean_model","PE_to_KE","PE_to_KE","energy","all",.true.,"none",2 -#"ocean_model","KE_Coradv","KE_Coradv","energy","all",.true.,"none",2 -#"ocean_model","KE_adv","KE_adv","energy","all",.true.,"none",2 -#"ocean_model","KE_visc","KE_visc","energy","all",.true.,"none",2 -#"ocean_model","KE_horvisc","KE_horvisc","energy","all",.true.,"none",2 -#"ocean_model","KE_dia","KE_dia","energy","all",.true.,"none",2 -# -# Mixed Layer TKE Budget Terms: -#=========================== -#"ocean_model","TKE_wind","TKE_wind","ML_TKE","all",.true.,"none",2 -#"ocean_model","TKE_RiBulk","TKE_RiBulk","ML_TKE","all",.true.,"none",2 -#"ocean_model","TKE_conv","TKE_conv","ML_TKE","all",.true.,"none",2 -#"ocean_model","TKE_pen_SW","TKE_pen_SW","ML_TKE","all",.true.,"none",2 -#"ocean_model","TKE_mixing","TKE_mixing","ML_TKE","all",.true.,"none",2 -#"ocean_model","TKE_mech_decay","TKE_mech_decay","ML_TKE","all",.true.,"none",2 -#"ocean_model","TKE_conv_decay","TKE_conv_decay","ML_TKE","all",.true.,"none",2 - -# Surface Forcing: -#================= -#"ocean_model","taux","taux","forcing","all",.true.,"none",2 -#"ocean_model","tauy","tauy","forcing","all",.true.,"none",2 -#"ocean_model","ustar","ustar","forcing","all",.true.,"none",2 -#"ocean_model","PRCmE","PRCmE","forcing","all",.true.,"none",2 -#"ocean_model","SW","SW","forcing","all",.true.,"none",2 -#"ocean_model","LwLatSens","LwLatSens","forcing","all",.true.,"none",2 -#"ocean_model","p_surf","p_surf","forcing","all",.true.,"none",2 -#"ocean_model","salt_flux","salt_flux","forcing","all",.true.,"none",2 -# - - -#============================================================================================= -# -#====> This file can be used with diag_manager/v2.0a (or higher) <==== -# -# -# FORMATS FOR FILE ENTRIES (not all input values are used) -# ------------------------ -# -#"file_name", output_freq, "output_units", format, "time_units", "time_long_name", ... -# (opt) new_file_frequecy, (opt) "new_file_freq_units", "new_file_start_date" -# -# -#output_freq: > 0 output frequency in "output_units" -# = 0 output frequency every time step -# =-1 output frequency at end of run -# -#output_units = units used for output frequency -# (years, months, days, minutes, hours, seconds) -# -#time_units = units used to label the time axis -# (days, minutes, hours, seconds) -# -# -# FORMAT FOR FIELD ENTRIES (not all input values are used) -# ------------------------ -# -#"module_name", "field_name", "output_name", "file_name" "time_sampling", time_avg, "other_opts", packing -# -#time_avg = .true. or .false. -# -#packing = 1 double precision -# = 2 float -# = 4 packed 16-bit integers -# = 8 packed 1-byte (not tested?) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 06634753e6..b8f88c43ff 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -1677,8 +1677,8 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & "faster by eliminating subroutine calls.", default=.false.) call get_param(param_file, "MOM", "DO_DYNAMICS", CS%do_dynamics, & "If False, skips the dynamics calls that update u & v, as well as "//& - "the gravity wave adjustment to h. This is a fragile feature and "//& - "thus undocumented.", default=.true., do_not_log=.true. ) + "the gravity wave adjustment to h. This may be a fragile feature, "//& + "but can be useful during development", default=.true.) call get_param(param_file, "MOM", "ADVECT_TS", advect_TS, & "If True, advect temperature and salinity horizontally "//& "If False, T/S are registered for advection. "//& @@ -2364,7 +2364,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & endif call tracer_advect_init(Time, G, param_file, diag, CS%tracer_adv_CSp) - call tracer_hor_diff_init(Time, G, US, param_file, diag, CS%tv%eqn_of_state, & + call tracer_hor_diff_init(Time, G, US, param_file, diag, CS%tv%eqn_of_state, CS%diabatic_CSp, & CS%tracer_diff_CSp) call lock_tracer_registry(CS%tracer_Reg) diff --git a/src/core/MOM_unit_tests.F90 b/src/core/MOM_unit_tests.F90 index ff5a93a62c..4197cfea3f 100644 --- a/src/core/MOM_unit_tests.F90 +++ b/src/core/MOM_unit_tests.F90 @@ -5,10 +5,11 @@ module MOM_unit_tests use MOM_error_handler, only : MOM_error, FATAL, is_root_pe -use MOM_string_functions, only : string_functions_unit_tests -use MOM_remapping, only : remapping_unit_tests -use MOM_neutral_diffusion, only : neutral_diffusion_unit_tests -use MOM_diag_vkernels, only : diag_vkernels_unit_tests +use MOM_string_functions, only : string_functions_unit_tests +use MOM_remapping, only : remapping_unit_tests +use MOM_neutral_diffusion, only : neutral_diffusion_unit_tests +use MOM_diag_vkernels, only : diag_vkernels_unit_tests +use MOM_lateral_boundary_diffusion, only : near_boundary_unit_tests implicit none ; private @@ -35,6 +36,9 @@ subroutine unit_tests(verbosity) "MOM_unit_tests: neutralDiffusionUnitTests FAILED") if (diag_vkernels_unit_tests(verbose)) call MOM_error(FATAL, & "MOM_unit_tests: diag_vkernels_unit_tests FAILED") + if (near_boundary_unit_tests(verbose)) call MOM_error(FATAL, & + "MOM_unit_tests: near_boundary_unit_tests FAILED") + endif end subroutine unit_tests diff --git a/src/equation_of_state/MOM_EOS.F90 b/src/equation_of_state/MOM_EOS.F90 index d3b056827b..0b966e8549 100644 --- a/src/equation_of_state/MOM_EOS.F90 +++ b/src/equation_of_state/MOM_EOS.F90 @@ -82,6 +82,11 @@ module MOM_EOS module procedure calculate_TFreeze_scalar, calculate_TFreeze_array end interface calculate_TFreeze +!> Calculates the compressibility of water from T, S, and P +interface calculate_compress + module procedure calculate_compress_scalar, calculate_compress_array +end interface calculate_compress + !> A control structure for the equation of state type, public :: EOS_type ; private integer :: form_of_EOS = 0 !< The equation of state to use. @@ -523,16 +528,16 @@ subroutine calculate_specific_vol_derivs(T, S, pressure, dSV_dT, dSV_dS, start, end subroutine calculate_specific_vol_derivs !> Calls the appropriate subroutine to calculate the density and compressibility for 1-D array inputs. -subroutine calculate_compress(T, S, pressure, rho, drho_dp, start, npts, EOS) - real, dimension(:), intent(in) :: T !< Potential temperature referenced to the surface [degC] - real, dimension(:), intent(in) :: S !< Salinity [ppt] +subroutine calculate_compress_array(T, S, pressure, rho, drho_dp, start, npts, EOS) + real, dimension(:), intent(in) :: T !< Potential temperature referenced to the surface [degC] + real, dimension(:), intent(in) :: S !< Salinity [PSU] real, dimension(:), intent(in) :: pressure !< Pressure [Pa] - real, dimension(:), intent(out) :: rho !< In situ density [kg m-3]. - real, dimension(:), intent(out) :: drho_dp !< The partial derivative of density with pressure - !! (also the inverse of the square of sound speed) in s2 m-2. - integer, intent(in) :: start !< Starting index within the array - integer, intent(in) :: npts !< The number of values to calculate - type(EOS_type), pointer :: EOS !< Equation of state structure + real, dimension(:), intent(out) :: rho !< In situ density [kg m-3]. + real, dimension(:), intent(out) :: drho_dp !< The partial derivative of density with pressure + !! (also the inverse of the square of sound speed) [s2 m-2]. + integer, intent(in) :: start !< Starting index within the array + integer, intent(in) :: npts !< The number of values to calculate + type(EOS_type), pointer :: EOS !< Equation of state structure if (.not.associated(EOS)) call MOM_error(FATAL, & "calculate_compress called with an unassociated EOS_type EOS.") @@ -554,8 +559,29 @@ subroutine calculate_compress(T, S, pressure, rho, drho_dp, start, npts, EOS) "calculate_compress: EOS%form_of_EOS is not valid.") end select -end subroutine calculate_compress +end subroutine calculate_compress_array + +!> Calculate density and compressibility for a scalar. This just promotes the scalar to an array with a singleton +!! dimension and calls calculate_compress_array +subroutine calculate_compress_scalar(T, S, pressure, rho, drho_dp, EOS) + real, intent(in) :: T !< Potential temperature referenced to the surface (degC) + real, intent(in) :: S !< Salinity (PSU) + real, intent(in) :: pressure !< Pressure (Pa) + real, intent(out) :: rho !< In situ density in kg m-3. + real, intent(out) :: drho_dp !< The partial derivative of density with pressure + !! (also the inverse of the square of sound speed) in s2 m-2. + type(EOS_type), pointer :: EOS !< Equation of state structure + + real, dimension(1) :: Ta, Sa, pa, rhoa, drho_dpa + + if (.not.associated(EOS)) call MOM_error(FATAL, & + "calculate_compress called with an unassociated EOS_type EOS.") + Ta(1) = T ; Sa(1) = S; pa(1) = pressure + + call calculate_compress_array(Ta, Sa, pa, rhoa, drho_dpa, 1, 1, EOS) + rho = rhoa(1) ; drho_dp = drho_dpa(1) +end subroutine calculate_compress_scalar !> Calls the appropriate subroutine to alculate analytical and nearly-analytical !! integrals in pressure across layers of geopotential anomalies, which are !! required for calculating the finite-volume form pressure accelerations in a diff --git a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 index 2ff0b3efe1..71b2df1e49 100644 --- a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 @@ -196,7 +196,6 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive, Waves) !! False => compute G'(1) as in LMD94 if (associated(CS)) call MOM_error(FATAL, 'MOM_CVMix_KPP, KPP_init: '// & 'Control structure has already been initialized') - allocate(CS) ! Read parameters call log_version(paramFile, mdl, version, 'This is the MOM wrapper to CVMix:KPP\n' // & @@ -207,6 +206,7 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive, Waves) default=.false.) ! Forego remainder of initialization if not using this scheme if (.not. KPP_init) return + allocate(CS) call openParameterBlock(paramFile,'KPP') call get_param(paramFile, mdl, 'PASSIVE', CS%passiveMode, & diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index e6f644d210..4a99bb9b2b 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -2864,12 +2864,14 @@ end subroutine layered_diabatic !> Returns pointers or values of members within the diabatic_CS type. For extensibility, !! each returned argument is an optional argument -subroutine extract_diabatic_member(CS, opacity_CSp, optics_CSp, & - evap_CFL_limit, minimum_forcing_depth, diabatic_aux_CSp) - type(diabatic_CS), intent(in ) :: CS !< module control structure +subroutine extract_diabatic_member(CS, opacity_CSp, optics_CSp, evap_CFL_limit, & + minimum_forcing_depth, KPP_CSp, energetic_PBL_CSp, diabatic_aux_CSp) + type(diabatic_CS), intent(in ) :: CS !< module control structure ! All output arguments are optional type(opacity_CS), optional, pointer :: opacity_CSp !< A pointer to be set to the opacity control structure type(optics_type), optional, pointer :: optics_CSp !< A pointer to be set to the optics control structure + type(KPP_CS), optional, pointer :: KPP_CSp !< A pointer to be set to the KPP CS + type(energetic_PBL_CS), optional, pointer :: energetic_PBL_CSp !< A pointer to be set to the ePBL CS real, optional, intent( out) :: evap_CFL_limit ! CS%opacity_CSp - if (present(optics_CSp)) optics_CSp => CS%optics - if (present(diabatic_aux_CSp)) diabatic_aux_CSp => CS%diabatic_aux_CSp + if (present(opacity_CSp)) opacity_CSp => CS%opacity_CSp + if (present(optics_CSp)) optics_CSp => CS%optics + if (present(KPP_CSp)) KPP_CSp => CS%KPP_CSp + if (present(energetic_PBL_CSp)) energetic_PBL_CSp => CS%energetic_PBL_CSp ! Constants within diabatic_CS if (present(evap_CFL_limit)) evap_CFL_limit = CS%evap_CFL_limit diff --git a/src/tracer/MOM_lateral_boundary_diffusion.F90 b/src/tracer/MOM_lateral_boundary_diffusion.F90 new file mode 100644 index 0000000000..07f062d1d2 --- /dev/null +++ b/src/tracer/MOM_lateral_boundary_diffusion.F90 @@ -0,0 +1,1218 @@ +!> Calculates and applies diffusive fluxes as a parameterization of lateral mixing (non-neutral) by +!! mesoscale eddies near the top and bottom (to be implemented) boundary layers of the ocean. + +module MOM_lateral_boundary_diffusion + +! This file is part of MOM6. See LICENSE.md for the license. + +use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end +use MOM_cpu_clock, only : CLOCK_MODULE, CLOCK_ROUTINE +use MOM_domains, only : pass_var +use MOM_diag_mediator, only : diag_ctrl, time_type +use MOM_diag_mediator, only : post_data, register_diag_field +use MOM_error_handler, only : MOM_error, FATAL, WARNING, MOM_mesg, is_root_pe +use MOM_file_parser, only : get_param, log_version, param_file_type +use MOM_file_parser, only : openParameterBlock, closeParameterBlock +use MOM_grid, only : ocean_grid_type +use MOM_remapping, only : remapping_CS, initialize_remapping +use MOM_remapping, only : extract_member_remapping_CS, build_reconstructions_1d +use MOM_remapping, only : average_value_ppoly, remappingSchemesDoc, remappingDefaultScheme +use MOM_tracer_registry, only : tracer_registry_type, tracer_type +use MOM_unit_scaling, only : unit_scale_type +use MOM_verticalGrid, only : verticalGrid_type +use MOM_CVMix_KPP, only : KPP_get_BLD, KPP_CS +use MOM_energetic_PBL, only : energetic_PBL_get_MLD, energetic_PBL_CS +use MOM_diabatic_driver, only : diabatic_CS, extract_diabatic_member + +implicit none ; private + +public near_boundary_unit_tests, lateral_boundary_diffusion, lateral_boundary_diffusion_init +public boundary_k_range + +! Private parameters to avoid doing string comparisons for bottom or top boundary layer +integer, public, parameter :: SURFACE = -1 !< Set a value that corresponds to the surface bopundary +integer, public, parameter :: BOTTOM = 1 !< Set a value that corresponds to the bottom boundary +#include + +!> Sets parameters for lateral boundary mixing module. +type, public :: lateral_boundary_diffusion_CS ; private + integer :: method !< Determine which of the three methods calculate + !! and apply near boundary layer fluxes + !! 1. Bulk-layer approach + !! 2. Along layer + !! 3. Decomposition onto pressure levels + integer :: deg !< Degree of polynomial reconstruction + integer :: surface_boundary_scheme !< Which boundary layer scheme to use + !! 1. ePBL; 2. KPP + type(remapping_CS) :: remap_CS !< Control structure to hold remapping configuration + type(KPP_CS), pointer :: KPP_CSp => NULL() !< KPP control structure needed to get BLD + type(energetic_PBL_CS), pointer :: energetic_PBL_CSp => NULL() !< ePBL control structure needed to get BLD + type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to + !! regulate the timing of diagnostic output. +end type lateral_boundary_diffusion_CS + +! This include declares and sets the variable "version". +#include "version_variable.h" +character(len=40) :: mdl = "MOM_lateral_boundary_diffusion" !< Name of this module + +contains + +!> Initialization routine that reads runtime parameters and sets up pointers to other control structures that might be +!! needed for lateral boundary diffusion. +logical function lateral_boundary_diffusion_init(Time, G, param_file, diag, diabatic_CSp, CS) + type(time_type), target, intent(in) :: Time !< Time structure + type(ocean_grid_type), intent(in) :: G !< Grid structure + type(param_file_type), intent(in) :: param_file !< Parameter file structure + type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure + type(diabatic_CS), pointer :: diabatic_CSp !< KPP control structure needed to get BLD + type(lateral_boundary_diffusion_CS), pointer :: CS !< Lateral boundary mixing control structure + + ! local variables + character(len=80) :: string ! Temporary strings + logical :: boundary_extrap + + if (ASSOCIATED(CS)) then + call MOM_error(FATAL, "lateral_boundary_diffusion_init called with associated control structure.") + return + endif + + ! Log this module and master switch for turning it on/off + call log_version(param_file, mdl, version, & + "This module implements lateral diffusion of tracers near boundaries") + call get_param(param_file, mdl, "USE_LATERAL_BOUNDARY_DIFFUSION", lateral_boundary_diffusion_init, & + "If true, enables the lateral boundary tracer's diffusion module.", & + default=.false.) + + if (.not. lateral_boundary_diffusion_init) then + return + endif + + allocate(CS) + CS%diag => diag + call extract_diabatic_member(diabatic_CSp, KPP_CSp=CS%KPP_CSp) + call extract_diabatic_member(diabatic_CSp, energetic_PBL_CSp=CS%energetic_PBL_CSp) + + CS%surface_boundary_scheme = -1 + if ( .not. ASSOCIATED(CS%energetic_PBL_CSp) .and. .not. ASSOCIATED(CS%KPP_CSp) ) then + call MOM_error(FATAL,"Lateral boundary diffusion is true, but no valid boundary layer scheme was found") + endif + + ! Read all relevant parameters and write them to the model log. + call get_param(param_file, mdl, "LATERAL_BOUNDARY_METHOD", CS%method, & + "Determine how to apply boundary lateral diffusion of tracers: \n"//& + "1. Bulk layer approach \n"//& + "2. Along layer approach \n"//& + "3. Decomposition on to pressure levels", default=1) + call get_param(param_file, mdl, "LBD_BOUNDARY_EXTRAP", boundary_extrap, & + "Use boundary extrapolation in LBD code", & + default=.false.) + call get_param(param_file, mdl, "LBD_REMAPPING_SCHEME", string, & + "This sets the reconstruction scheme used "//& + "for vertical remapping for all variables. "//& + "It can be one of the following schemes: "//& + trim(remappingSchemesDoc), default=remappingDefaultScheme) + call initialize_remapping( CS%remap_CS, string, boundary_extrapolation = boundary_extrap ) + call extract_member_remapping_CS(CS%remap_CS, degree=CS%deg) + +end function lateral_boundary_diffusion_init + +!> Driver routine for calculating lateral diffusive fluxes near the top and bottom boundaries. Two different methods +!! Method 1: Calculate fluxes from bulk layer integrated quantities +subroutine lateral_boundary_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS) + type(ocean_grid_type), intent(inout) :: G !< Grid type + type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & + intent(in) :: h !< Layer thickness [H ~> m or kg m-2] + real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: Coef_x !< dt * Kh * dy / dx at u-points [m2] + real, dimension(SZI_(G),SZJB_(G)), intent(in) :: Coef_y !< dt * Kh * dx / dy at v-points [m2] + real, intent(in) :: dt !< Tracer time step * I_numitts + !! (I_numitts in tracer_hordiff) + type(tracer_registry_type), pointer :: Reg !< Tracer registry + type(lateral_boundary_diffusion_CS), intent(in) :: CS !< Control structure for this module + + ! Local variables + real, dimension(SZI_(G),SZJ_(G)) :: hbl !< bnd. layer depth [m] + real, dimension(SZI_(G),SZJ_(G),SZK_(G),CS%deg+1) :: ppoly0_coefs !< Coefficients of polynomial + real, dimension(SZI_(G),SZJ_(G),SZK_(G),2) :: ppoly0_E !< Edge values from reconstructions + real, dimension(SZK_(G),CS%deg+1) :: ppoly_S !< Slopes from reconstruction (placeholder) + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: uFlx !< Zonal flux of tracer [H conc ~> m conc or conc kg m-2] + real, dimension(SZIB_(G),SZJ_(G)) :: uFLx_bulk !< Total calculated bulk-layer u-flux for the tracer + real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: vFlx !< Meridional flux of tracer + real, dimension(SZI_(G),SZJB_(G)) :: vFlx_bulk !< Total calculated bulk-layer v-flux for the tracer + real, dimension(SZIB_(G),SZJ_(G)) :: uwork_2d !< Layer summed u-flux transport + real, dimension(SZI_(G),SZJB_(G)) :: vwork_2d !< Layer summed v-flux transport + type(tracer_type), pointer :: Tracer => NULL() !< Pointer to the current tracer + integer :: remap_method !< Reconstruction method + integer :: i,j,k,m !< indices to loop over + real :: Idt !< inverse of the time step [s-1] + + Idt = 1./dt + hbl(:,:) = 0. + if (ASSOCIATED(CS%KPP_CSp)) call KPP_get_BLD(CS%KPP_CSp, hbl, G) + if (ASSOCIATED(CS%energetic_PBL_CSp)) call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, hbl, G, US) + + call pass_var(hbl,G%Domain) + + do m = 1,Reg%ntr + tracer => Reg%tr(m) + do j = G%jsc-1, G%jec+1 + ! Interpolate state to interface + do i = G%isc-1, G%iec+1 + call build_reconstructions_1d( CS%remap_CS, G%ke, h(i,j,:), tracer%t(i,j,:), ppoly0_coefs(i,j,:,:), & + ppoly0_E(i,j,:,:), ppoly_S, remap_method, GV%H_subroundoff, GV%H_subroundoff) + enddo + enddo + ! Diffusive fluxes in the i-direction + uFlx(:,:,:) = 0. + vFlx(:,:,:) = 0. + uFlx_bulk(:,:) = 0. + vFlx_bulk(:,:) = 0. + + ! Method #1 + if ( CS%method == 1 ) then + do j=G%jsc,G%jec + do i=G%isc-1,G%iec + if (G%mask2dCu(I,j)>0.) then + call fluxes_bulk_method(SURFACE, GV%ke, CS%deg, h(I,j,:), h(I+1,j,:), hbl(I,j), hbl(I+1,j), & + G%areaT(I,j), G%areaT(I+1,j), tracer%t(I,j,:), tracer%t(I+1,j,:), & + ppoly0_coefs(I,j,:,:), ppoly0_coefs(I+1,j,:,:), ppoly0_E(I,j,:,:), & + ppoly0_E(I+1,j,:,:), remap_method, Coef_x(I,j), uFlx_bulk(I,j), uFlx(I,j,:)) + endif + enddo + enddo + do J=G%jsc-1,G%jec + do i=G%isc,G%iec + if (G%mask2dCv(i,J)>0.) then + call fluxes_bulk_method(SURFACE, GV%ke, CS%deg, h(i,J,:), h(i,J+1,:), hbl(i,J), hbl(i,J+1), & + G%areaT(i,J), G%areaT(i,J+1), tracer%t(i,J,:), tracer%t(i,J+1,:), & + ppoly0_coefs(i,J,:,:), ppoly0_coefs(i,J+1,:,:), ppoly0_E(i,J,:,:), & + ppoly0_E(i,J+1,:,:), remap_method, Coef_y(i,J), vFlx_bulk(i,J), vFlx(i,J,:)) + endif + enddo + enddo + ! Post tracer bulk diags + if (tracer%id_lbd_bulk_dfx>0) call post_data(tracer%id_lbd_bulk_dfx, uFlx_bulk*Idt, CS%diag) + if (tracer%id_lbd_bulk_dfy>0) call post_data(tracer%id_lbd_bulk_dfy, vFlx_bulk*Idt, CS%diag) + + ! TODO: this is where we would filter vFlx and uFlux to get rid of checkerboard noise + + ! Method #2 + elseif (CS%method == 2) then + do j=G%jsc,G%jec + do i=G%isc-1,G%iec + if (G%mask2dCu(I,j)>0.) then + call fluxes_layer_method(SURFACE, GV%ke, CS%deg, h(i,j,:), h(i+1,j,:), hbl(i,j), hbl(i+1,j), & + tracer%t(i,j,:), tracer%t(i+1,j,:), ppoly0_coefs(i,j,:,:), ppoly0_coefs(i+1,j,:,:), ppoly0_E(i,j,:,:), & + ppoly0_E(i+1,j,:,:), remap_method, Coef_x(I,j), uFlx(I,j,:)) + endif + enddo + enddo + do J=G%jsc-1,G%jec + do i=G%isc,G%iec + if (G%mask2dCv(i,J)>0.) then + call fluxes_layer_method(SURFACE, GV%ke, CS%deg, h(i,J,:), h(i,J+1,:), hbl(i,J), hbl(i,J+1), & + tracer%t(i,J,:), tracer%t(i,J+1,:), ppoly0_coefs(i,J,:,:), ppoly0_coefs(i,J+1,:,:), ppoly0_E(i,J,:,:), & + ppoly0_E(i,J+1,:,:), remap_method, Coef_y(i,J), vFlx(i,J,:)) + endif + enddo + enddo + endif + + ! Update the tracer fluxes + do k=1,GV%ke ; do j=G%jsc,G%jec ; do i=G%isc,G%iec + if (G%mask2dT(i,j)>0.) then + tracer%t(i,j,k) = tracer%t(i,j,k) + (( (uFlx(I-1,j,k)-uFlx(I,j,k)) ) + ( (vFlx(i,J-1,k)-vFlx(i,J,k) ) ))* & + (G%IareaT(i,j)/( h(i,j,k) + GV%H_subroundoff)) + endif + enddo ; enddo ; enddo + + ! Post the tracer diagnostics + if (tracer%id_lbd_dfx>0) call post_data(tracer%id_lbd_dfx, uFlx*Idt, CS%diag) + if (tracer%id_lbd_dfy>0) call post_data(tracer%id_lbd_dfy, vFlx*Idt, CS%diag) + if (tracer%id_lbd_dfx_2d>0) then + uwork_2d(:,:) = 0. + do k=1,GV%ke; do j=G%jsc,G%jec; do I=G%isc-1,G%iec + uwork_2d(I,j) = uwork_2d(I,j) + (uFlx(I,j,k) * Idt) + enddo; enddo; enddo + call post_data(tracer%id_lbd_dfx_2d, uwork_2d, CS%diag) + endif + + if (tracer%id_lbd_dfy_2d>0) then + vwork_2d(:,:) = 0. + do k=1,GV%ke; do J=G%jsc-1,G%jec; do i=G%isc,G%iec + vwork_2d(i,J) = vwork_2d(i,J) + (vFlx(i,J,k) * Idt) + enddo; enddo; enddo + call post_data(tracer%id_lbd_dfy_2d, vwork_2d, CS%diag) + endif + enddo + +end subroutine lateral_boundary_diffusion + +!< Calculate bulk layer value of a scalar quantity as the thickness weighted average +real function bulk_average(boundary, nk, deg, h, hBLT, phi, ppoly0_E, ppoly0_coefs, method, k_top, zeta_top, k_bot, & + zeta_bot) + integer :: boundary !< SURFACE or BOTTOM [nondim] + integer :: nk !< Number of layers [nondim] + integer :: deg !< Degree of polynomial [nondim] + real, dimension(nk) :: h !< Layer thicknesses [m] + real :: hBLT !< Depth of the boundary layer [m] + real, dimension(nk) :: phi !< Scalar quantity + real, dimension(nk,2) :: ppoly0_E(:,:) !< Edge value of polynomial + real, dimension(nk,deg+1) :: ppoly0_coefs(:,:) !< Coefficients of polynomial + integer :: method !< Remapping scheme to use + + integer :: k_top !< Index of the first layer within the boundary + real :: zeta_top !< Fraction of the layer encompassed by the bottom boundary layer + !! (0 if none, 1. if all). For the surface, this is always 0. because + !! integration starts at the surface [nondim] + integer :: k_bot !< Index of the last layer within the boundary + real :: zeta_bot !< Fraction of the layer encompassed by the surface boundary layer + !! (0 if none, 1. if all). For the bottom boundary layer, this is always 1. + !! because integration starts at the bottom [nondim] + ! Local variables + real :: htot ! Running sum of the thicknesses (top to bottom) + integer :: k ! k indice + + + htot = 0. + bulk_average = 0. + if (hblt == 0.) return + if (boundary == SURFACE) then + htot = (h(k_bot) * zeta_bot) + bulk_average = average_value_ppoly( nk, phi, ppoly0_E, ppoly0_coefs, method, k_bot, 0., zeta_bot) * htot + do k = k_bot-1,1,-1 + bulk_average = bulk_average + phi(k)*h(k) + htot = htot + h(k) + enddo + elseif (boundary == BOTTOM) then + htot = (h(k_top) * zeta_top) + ! (note 1-zeta_top because zeta_top is the fraction of the layer) + bulk_average = average_value_ppoly( nk, phi, ppoly0_E, ppoly0_coefs, method, k_top, (1.-zeta_top), 1.) * htot + do k = k_top+1,nk + bulk_average = bulk_average + phi(k)*h(k) + htot = htot + h(k) + enddo + else + call MOM_error(FATAL, "bulk_average: a valid boundary type must be provided.") + endif + + bulk_average = bulk_average / hBLT + +end function bulk_average + +!> Calculate the harmonic mean of two quantities +!! See \ref section_harmonic_mean. +real function harmonic_mean(h1,h2) + real :: h1 !< Scalar quantity + real :: h2 !< Scalar quantity + if (h1 + h2 == 0.) then + harmonic_mean = 0. + else + harmonic_mean = 2.*(h1*h2)/(h1+h2) + endif +end function harmonic_mean + +!> Find the k-index range corresponding to the layers that are within the boundary-layer region +subroutine boundary_k_range(boundary, nk, h, hbl, k_top, zeta_top, k_bot, zeta_bot) + integer, intent(in ) :: boundary !< SURFACE or BOTTOM [nondim] + integer, intent(in ) :: nk !< Number of layers [nondim] + real, dimension(nk), intent(in ) :: h !< Layer thicknesses of the coluymn [m] + real, intent(in ) :: hbl !< Thickness of the boundary layer [m] + !! If surface, with respect to zbl_ref = 0. + !! If bottom, with respect to zbl_ref = SUM(h) + integer, intent( out) :: k_top !< Index of the first layer within the boundary + real, intent( out) :: zeta_top !< Distance from the top of a layer to the intersection of the + !! top extent of the boundary layer (0 at top, 1 at bottom) [nondim] + integer, intent( out) :: k_bot !< Index of the last layer within the boundary + real, intent( out) :: zeta_bot !< Distance of the lower layer to the boundary layer depth + !! (0 at top, 1 at bottom) [nondim] + ! Local variables + real :: htot + integer :: k + ! Surface boundary layer + if ( boundary == SURFACE ) then + k_top = 1 + zeta_top = 0. + htot = 0. + k_bot = 1 + zeta_bot = 0. + if (hbl == 0.) return + do k=1,nk + htot = htot + h(k) + if ( htot >= hbl) then + k_bot = k + zeta_bot = 1 - (htot - hbl)/h(k) + return + endif + enddo + ! Bottom boundary layer + elseif ( boundary == BOTTOM ) then + k_top = nk + zeta_top = 1. + k_bot = nk + zeta_bot = 1. + htot = 0. + if (hbl == 0.) return + do k=nk,1,-1 + htot = htot + h(k) + if (htot >= hbl) then + k_top = k + zeta_top = 1 - (htot - hbl)/h(k) + return + endif + enddo + else + call MOM_error(FATAL,"Houston, we've had a problem in boundary_k_range") + endif + +end subroutine boundary_k_range + + +!> Calculate the lateral boundary diffusive fluxes using the layer by layer method. +!! See \ref section_method2 +subroutine fluxes_layer_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, phi_L, phi_R, ppoly0_coefs_L, & + ppoly0_coefs_R, ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_layer) + + integer, intent(in ) :: boundary !< Which boundary layer SURFACE or BOTTOM [nondim] + integer, intent(in ) :: nk !< Number of layers [nondim] + integer, intent(in ) :: deg !< order of the polynomial reconstruction [nondim] + real, dimension(nk), intent(in ) :: h_L !< Layer thickness (left) [m] + real, dimension(nk), intent(in ) :: h_R !< Layer thickness (right) [m] + real, intent(in ) :: hbl_L !< Thickness of the boundary boundary + !! layer (left) [m] + real, intent(in ) :: hbl_R !< Thickness of the boundary boundary + !! layer (left) [m] + real, dimension(nk), intent(in ) :: phi_L !< Tracer values (left) [ nondim m^-3 ] + real, dimension(nk), intent(in ) :: phi_R !< Tracer values (right) [ nondim m^-3 ] + real, dimension(nk,deg+1), intent(in ) :: ppoly0_coefs_L !< Tracer reconstruction (left) [ nondim m^-3 ] + real, dimension(nk,deg+1), intent(in ) :: ppoly0_coefs_R !< Tracer reconstruction (right) [ nondim m^-3 ] + real, dimension(nk,2), intent(in ) :: ppoly0_E_L !< Polynomial edge values (left) [ nondim ] + real, dimension(nk,2), intent(in ) :: ppoly0_E_R !< Polynomial edge values (right) [ nondim ] + integer, intent(in ) :: method !< Method of polynomial integration [ nondim ] + real, intent(in ) :: khtr_u !< Horizontal diffusivities times delta t at U-point [m^2] + real, dimension(nk), intent( out) :: F_layer !< Layerwise diffusive flux at U-point [m^2 conc] + ! Local variables + real, dimension(nk) :: h_means ! Calculate the layer-wise harmonic means [m] + real, dimension(nk) :: h_u ! Thickness at the u-point [m] + real :: hbl_u ! Boundary layer Thickness at the u-point [m] + real :: khtr_avg ! Thickness-weighted diffusivity at the u-point [m^2 s^-1] + real :: heff ! Harmonic mean of layer thicknesses [m] + real :: inv_heff ! Inverse of the harmonic mean of layer thicknesses [m^[-1] + real :: phi_L_avg, phi_R_avg ! Bulk, thickness-weighted tracer averages (left and right column) + ! [conc m^-3 ] + real :: htot ! Total column thickness [m] + integer :: k, k_bot_min, k_top_max + integer :: k_top_L, k_bot_L, k_top_u + integer :: k_top_R, k_bot_R, k_bot_u + real :: zeta_top_L, zeta_top_R, zeta_top_u + real :: zeta_bot_L, zeta_bot_R, zeta_bot_u + real :: h_work_L, h_work_R ! dummy variables + real :: hbl_min ! minimum BLD (left and right) + + F_layer(:) = 0.0 + if (hbl_L == 0. .or. hbl_R == 0.) then + return + endif + + ! Calculate vertical indices containing the boundary layer + call boundary_k_range(boundary, nk, h_L, hbl_L, k_top_L, zeta_top_L, k_bot_L, zeta_bot_L) + call boundary_k_range(boundary, nk, h_R, hbl_R, k_top_R, zeta_top_R, k_bot_R, zeta_bot_R) + + if (boundary == SURFACE) then + k_bot_min = MIN(k_bot_L, k_bot_R) + ! make sure left and right k indices span same range + if (k_bot_min .ne. k_bot_L) then + k_bot_L = k_bot_min + zeta_bot_L = 1.0 + endif + if (k_bot_min .ne. k_bot_R) then + k_bot_R= k_bot_min + zeta_bot_R = 1.0 + endif + + h_work_L = (h_L(k_bot_L) * zeta_bot_L) + h_work_R = (h_R(k_bot_R) * zeta_bot_R) + + phi_L_avg = average_value_ppoly( nk, phi_L, ppoly0_E_L, ppoly0_coefs_L, method, k_bot_L, 0., zeta_bot_L) + phi_R_avg = average_value_ppoly( nk, phi_R, ppoly0_E_R, ppoly0_coefs_R, method, k_bot_R, 0., zeta_bot_R) + heff = harmonic_mean(h_work_L, h_work_R) + ! tracer flux where the minimum BLD intersets layer + F_layer(k_bot_min) = -heff * (phi_R_avg - phi_L_avg) + do k = k_bot_min-1,1,-1 + heff = harmonic_mean(h_L(k), h_R(k)) + F_layer(k) = -heff * (phi_R(k) - phi_L(k)) + enddo + endif + + if (boundary == BOTTOM) then + k_top_max = MAX(k_top_L, k_top_R) + ! make sure left and right k indices span same range + if (k_top_max .ne. k_top_L) then + k_top_L = k_top_max + zeta_top_L = 1.0 + endif + if (k_top_max .ne. k_top_R) then + k_top_R= k_top_max + zeta_top_R = 1.0 + endif + + h_work_L = (h_L(k_top_L) * zeta_top_L) + h_work_R = (h_R(k_top_R) * zeta_top_R) + + phi_L_avg = average_value_ppoly( nk, phi_L, ppoly0_E_L, ppoly0_coefs_L, method, k_top_L, 1.0-zeta_top_L, 1.0) + phi_R_avg = average_value_ppoly( nk, phi_R, ppoly0_E_R, ppoly0_coefs_R, method, k_top_R, 1.0-zeta_top_R, 1.0) + heff = harmonic_mean(h_work_L, h_work_R) + + ! tracer flux where the minimum BLD intersets layer + F_layer(k_top_max) = (-heff * khtr_u) * (phi_R_avg - phi_L_avg) + do k = k_top_max+1,nk + heff = harmonic_mean(h_L(k), h_R(k)) + F_layer(k) = -(heff * khtr_u) * (phi_R(k) - phi_L(k)) + enddo + endif + +end subroutine fluxes_layer_method + +!> Apply the lateral boundary diffusive fluxes calculated from a 'bulk model' +!! See \ref section_method1 +subroutine fluxes_bulk_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, ppoly0_coefs_L, & + ppoly0_coefs_R, ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_bulk, F_layer, F_limit) + + integer, intent(in ) :: boundary !< Which boundary layer SURFACE or BOTTOM [nondim] + integer, intent(in ) :: nk !< Number of layers [nondim] + integer, intent(in ) :: deg !< order of the polynomial reconstruction [nondim] + real, dimension(nk), intent(in ) :: h_L !< Layer thickness (left) [m] + real, dimension(nk), intent(in ) :: h_R !< Layer thickness (right) [m] + real, intent(in ) :: hbl_L !< Thickness of the boundary boundary + !! layer (left) [m] + real, intent(in ) :: hbl_R !< Thickness of the boundary boundary + !! layer (left) [m] + real, intent(in ) :: area_L !< Area of the horizontal grid (left) [m^2] + real, intent(in ) :: area_R !< Area of the horizontal grid (right) [m^2] + real, dimension(nk), intent(in ) :: phi_L !< Tracer values (left) [ nondim m^-3 ] + real, dimension(nk), intent(in ) :: phi_R !< Tracer values (right) [ nondim m^-3 ] + real, dimension(nk,deg+1), intent(in ) :: ppoly0_coefs_L !< Tracer reconstruction (left) [ nondim m^-3 ] + real, dimension(nk,deg+1), intent(in ) :: ppoly0_coefs_R !< Tracer reconstruction (right) [ nondim m^-3 ] + real, dimension(nk,2), intent(in ) :: ppoly0_E_L !< Polynomial edge values (left) [ nondim ] + real, dimension(nk,2), intent(in ) :: ppoly0_E_R !< Polynomial edge values (right) [ nondim ] + integer, intent(in ) :: method !< Method of polynomial integration [ nondim ] + real, intent(in ) :: khtr_u !< Horizontal diffusivities times delta t at U-point [m^2] + real, intent( out) :: F_bulk !< The bulk mixed layer lateral flux [m^2 conc] + real, dimension(nk), intent( out) :: F_layer !< Layerwise diffusive flux at U-point [m^2 conc] + real, optional, dimension(nk), intent( out) :: F_limit !< The amount of flux not applied due to limiter + !! F_layer(k) - F_max [m^2 conc] + ! Local variables + real, dimension(nk) :: h_means ! Calculate the layer-wise harmonic means [m] + real, dimension(nk) :: h_u ! Thickness at the u-point [m] + real :: hbl_u ! Boundary layer Thickness at the u-point [m] + real :: khtr_avg ! Thickness-weighted diffusivity at the u-point [m^2 s^-1] + real :: heff ! Harmonic mean of layer thicknesses [m] + real :: inv_heff ! Inverse of the harmonic mean of layer thicknesses [m^[-1] + real :: phi_L_avg, phi_R_avg ! Bulk, thickness-weighted tracer averages (left and right column) + ! [conc m^-3 ] + real :: htot ! Total column thickness [m] + integer :: k, k_min, k_max + integer :: k_top_L, k_bot_L, k_top_u + integer :: k_top_R, k_bot_R, k_bot_u + real :: zeta_top_L, zeta_top_R, zeta_top_u + real :: zeta_bot_L, zeta_bot_R, zeta_bot_u + real :: h_work_L, h_work_R ! dummy variables + real :: F_max !< The maximum amount of flux that can leave a cell + logical :: limited !< True if the flux limiter was applied + real :: hfrac, F_bulk_remain + + if (hbl_L == 0. .or. hbl_R == 0.) then + F_bulk = 0. + F_layer(:) = 0. + return + endif + + ! Calculate vertical indices containing the boundary layer + call boundary_k_range(boundary, nk, h_L, hbl_L, k_top_L, zeta_top_L, k_bot_L, zeta_bot_L) + call boundary_k_range(boundary, nk, h_R, hbl_R, k_top_R, zeta_top_R, k_bot_R, zeta_bot_R) + + ! Calculate bulk averages of various quantities + phi_L_avg = bulk_average(boundary, nk, deg, h_L, hbl_L, phi_L, ppoly0_E_L, ppoly0_coefs_L, method, k_top_L, & + zeta_top_L, k_bot_L, zeta_bot_L) + phi_R_avg = bulk_average(boundary, nk, deg, h_R, hbl_R, phi_R, ppoly0_E_R, ppoly0_coefs_R, method, k_top_R, & + zeta_top_R, k_bot_R, zeta_bot_R) + do k=1,nk + h_u(k) = 0.5 * (h_L(k) + h_R(k)) + enddo + hbl_u = 0.5*(hbl_L + hbl_R) + call boundary_k_range(boundary, nk, h_u, hbl_u, k_top_u, zeta_top_u, k_bot_u, zeta_bot_u) + + ! Calculate the 'bulk' diffusive flux from the bulk averaged quantities + heff = harmonic_mean(hbl_L, hbl_R) + F_bulk = -(khtr_u * heff) * (phi_R_avg - phi_L_avg) + F_bulk_remain = F_bulk + ! Calculate the layerwise sum of the vertical effective thickness. This is different than the heff calculated + ! above, but is used as a way to decompose the fluxes onto the individual layers + h_means(:) = 0. + + if (boundary == SURFACE) then + k_min = MIN(k_bot_L, k_bot_R) + + ! left hand side + if (k_bot_L == k_min) then + h_work_L = h_L(k_min) * zeta_bot_L + else + h_work_L = h_L(k_min) + endif + + ! right hand side + if (k_bot_R == k_min) then + h_work_R = h_R(k_min) * zeta_bot_R + else + h_work_R = h_R(k_min) + endif + + h_means(k_min) = harmonic_mean(h_work_L,h_work_R) + + do k=1,k_min-1 + h_means(k) = harmonic_mean(h_L(k),h_R(k)) + enddo + + elseif (boundary == BOTTOM) then + k_max = MAX(k_top_L, k_top_R) + ! left hand side + if (k_top_L == k_max) then + h_work_L = h_L(k_max) * zeta_top_L + else + h_work_L = h_L(k_max) + endif + + ! right hand side + if (k_top_R == k_max) then + h_work_R = h_R(k_max) * zeta_top_R + else + h_work_R = h_R(k_max) + endif + + h_means(k_max) = harmonic_mean(h_work_L,h_work_R) + + do k=nk,k_max+1,-1 + h_means(k) = harmonic_mean(h_L(k),h_R(k)) + enddo + endif + + if ( SUM(h_means) == 0. ) then + return + else + ! Initialize remaining thickness + inv_heff = 1./SUM(h_means) + ! Decompose the bulk flux onto the individual layers + do k=1,nk + if (h_means(k) > 0.) then + ! Limit the tracer flux so that the donor cell with positive concentration can't go negative + ! If a tracer can go negative, it is unclear what the limiter should be. BOB ALISTAIR?! + hfrac = h_means(k)*inv_heff + F_layer(k) = F_bulk * hfrac + if ( SIGN(1.,F_bulk) == SIGN(1., F_layer(k))) then + ! limit the flux to 0.25 of the total tracer in the cell + if (F_bulk < 0. .and. phi_R(k) >= 0.) then + F_max = 0.25 * (area_R*(phi_R(k)*h_R(k))) + elseif (F_bulk > 0. .and. phi_L(k) >= 0.) then + F_max = 0.25 * (area_L*(phi_L(k)*h_L(k))) + else ! The above quantities are always positive, so we can use F_max < -1 to see if we don't need to limit + F_max = -1. + endif + ! Distribute bulk flux onto layers + if ( ((boundary == SURFACE) .and. (k == k_min)) .or. ((boundary == BOTTOM) .and. (k == nk)) ) then + F_layer(k) = F_bulk_remain + endif + F_bulk_remain = F_bulk_remain - F_layer(k) + + ! Apply flux limiter calculated above + if (F_max >= 0.) then + if (F_layer(k) > 0.) then + limited = F_layer(k) > F_max + F_layer(k) = MIN(F_layer(k),F_max) + elseif (F_layer(k) < 0.) then + limited = F_layer(k) < -F_max + F_layer(k) = MAX(F_layer(k),-F_max) ! Note negative to make the sign of flux consistent + endif + endif + + if (PRESENT(F_limit)) then + if (limited) then + F_limit(k) = F_layer(k) - F_max + else + F_limit(k) = 0. + endif + endif + else + F_bulk_remain = F_bulk_remain - F_layer(k) + F_layer(k) = 0. + endif + else + F_layer(k) = 0. + endif + enddo + endif + +end subroutine fluxes_bulk_method + +! TODO: GMM, this is a placeholder for the pressure reconstruction. +! get rid of all the T/S related variables below. We need to use the +! continuous version since pressure will be continuous. However, +! for tracer we will need to use a discontinuous reconstruction. +! Mimic the neutral diffusion driver to calculate and apply sub-layer +! fluxes. + +!> Returns positions within left/right columns of combined interfaces using continuous reconstructions of T/S +!subroutine find_neutral_surface_positions_continuous(nk, Pl, Pr, PoL, PoR, KoL, KoR, hEff) +! integer, intent(in) :: nk !< Number of levels +! real, dimension(nk+1), intent(in) :: Pl !< Left-column interface pressure [Pa] +! real, dimension(2*nk+2), intent(inout) :: PoL !< Fractional position of neutral surface within +! !! layer KoL of left column +! real, dimension(2*nk+2), intent(inout) :: PoR !< Fractional position of neutral surface within +! !! layer KoR of right column +! integer, dimension(2*nk+2), intent(inout) :: KoL !< Index of first left interface above neutral surface +! integer, dimension(2*nk+2), intent(inout) :: KoR !< Index of first right interface above neutral surface +! real, dimension(2*nk+1), intent(inout) :: hEff !< Effective thickness between two neutral surfaces [Pa] +! +! ! Local variables +! integer :: ns ! Number of neutral surfaces +! integer :: k_surface ! Index of neutral surface +! integer :: kl ! Index of left interface +! integer :: kr ! Index of right interface +! real :: dRdT, dRdS ! dRho/dT and dRho/dS for the neutral surface +! logical :: searching_left_column ! True if searching for the position of a right interface in the left column +! logical :: searching_right_column ! True if searching for the position of a left interface in the right column +! logical :: reached_bottom ! True if one of the bottom-most interfaces has been used as the target +! integer :: krm1, klm1 +! real :: dRho, dRhoTop, dRhoBot, hL, hR +! integer :: lastK_left, lastK_right +! real :: lastP_left, lastP_right +! +! ns = 2*nk+2 +! ! Initialize variables for the search +! kr = 1 ; lastK_right = 1 ; lastP_right = 0. +! kl = 1 ; lastK_left = 1 ; lastP_left = 0. +! reached_bottom = .false. +! +! ! Loop over each neutral surface, working from top to bottom +! neutral_surfaces: do k_surface = 1, ns +! klm1 = max(kl-1, 1) +! if (klm1>nk) stop 'find_neutral_surface_positions(): klm1 went out of bounds!' +! krm1 = max(kr-1, 1) +! if (krm1>nk) stop 'find_neutral_surface_positions(): krm1 went out of bounds!' +! +! ! TODO: GMM, instead of dRho we need dP (pressure at right - pressure at left) +! +! ! Potential density difference, rho(kr) - rho(kl) +! dRho = 0.5 * ( ( dRdTr(kr) + dRdTl(kl) ) * ( Tr(kr) - Tl(kl) ) & +! + ( dRdSr(kr) + dRdSl(kl) ) * ( Sr(kr) - Sl(kl) ) ) +! ! Which column has the lighter surface for the current indexes, kr and kl +! if (.not. reached_bottom) then +! if (dRho < 0.) then +! searching_left_column = .true. +! searching_right_column = .false. +! elseif (dRho > 0.) then +! searching_right_column = .true. +! searching_left_column = .false. +! else ! dRho == 0. +! if (kl + kr == 2) then ! Still at surface +! searching_left_column = .true. +! searching_right_column = .false. +! else ! Not the surface so we simply change direction +! searching_left_column = .not. searching_left_column +! searching_right_column = .not. searching_right_column +! endif +! endif +! endif +! +! if (searching_left_column) then +! ! Interpolate for the neutral surface position within the left column, layer klm1 +! ! Potential density difference, rho(kl-1) - rho(kr) (should be negative) +! dRhoTop = 0.5 * ( ( dRdTl(klm1) + dRdTr(kr) ) * ( Tl(klm1) - Tr(kr) ) & +! + ( dRdSl(klm1) + dRdSr(kr) ) * ( Sl(klm1) - Sr(kr) ) ) +! ! Potential density difference, rho(kl) - rho(kr) (will be positive) +! dRhoBot = 0.5 * ( ( dRdTl(klm1+1) + dRdTr(kr) ) * ( Tl(klm1+1) - Tr(kr) ) & +! + ( dRdSl(klm1+1) + dRdSr(kr) ) * ( Sl(klm1+1) - Sr(kr) ) ) +! +! ! Because we are looking left, the right surface, kr, is lighter than klm1+1 and should be denser than klm1 +! ! unless we are still at the top of the left column (kl=1) +! if (dRhoTop > 0. .or. kr+kl==2) then +! PoL(k_surface) = 0. ! The right surface is lighter than anything in layer klm1 +! elseif (dRhoTop >= dRhoBot) then ! Left layer is unstratified +! PoL(k_surface) = 1. +! else +! ! Linearly interpolate for the position between Pl(kl-1) and Pl(kl) where the density difference +! ! between right and left is zero. +! +! ! TODO: GMM, write the linear solution instead of using interpolate_for_nondim_position +! PoL(k_surface) = interpolate_for_nondim_position( dRhoTop, Pl(klm1), dRhoBot, Pl(klm1+1) ) +! endif +! if (PoL(k_surface)>=1. .and. klm1= is really ==, when PoL==1 we point to the bottom of the cell +! klm1 = klm1 + 1 +! PoL(k_surface) = PoL(k_surface) - 1. +! endif +! if (real(klm1-lastK_left)+(PoL(k_surface)-lastP_left)<0.) then +! PoL(k_surface) = lastP_left +! klm1 = lastK_left +! endif +! KoL(k_surface) = klm1 +! if (kr <= nk) then +! PoR(k_surface) = 0. +! KoR(k_surface) = kr +! else +! PoR(k_surface) = 1. +! KoR(k_surface) = nk +! endif +! if (kr <= nk) then +! kr = kr + 1 +! else +! reached_bottom = .true. +! searching_right_column = .true. +! searching_left_column = .false. +! endif +! elseif (searching_right_column) then +! ! Interpolate for the neutral surface position within the right column, layer krm1 +! ! Potential density difference, rho(kr-1) - rho(kl) (should be negative) +! dRhoTop = 0.5 * ( ( dRdTr(krm1) + dRdTl(kl) ) * ( Tr(krm1) - Tl(kl) ) & +! + ( dRdSr(krm1) + dRdSl(kl) ) * ( Sr(krm1) - Sl(kl) ) ) +! ! Potential density difference, rho(kr) - rho(kl) (will be positive) +! dRhoBot = 0.5 * ( ( dRdTr(krm1+1) + dRdTl(kl) ) * ( Tr(krm1+1) - Tl(kl) ) & +! + ( dRdSr(krm1+1) + dRdSl(kl) ) * ( Sr(krm1+1) - Sl(kl) ) ) +! +! ! Because we are looking right, the left surface, kl, is lighter than krm1+1 and should be denser than krm1 +! ! unless we are still at the top of the right column (kr=1) +! if (dRhoTop >= 0. .or. kr+kl==2) then +! PoR(k_surface) = 0. ! The left surface is lighter than anything in layer krm1 +! elseif (dRhoTop >= dRhoBot) then ! Right layer is unstratified +! PoR(k_surface) = 1. +! else +! ! Linearly interpolate for the position between Pr(kr-1) and Pr(kr) where the density difference +! ! between right and left is zero. +! PoR(k_surface) = interpolate_for_nondim_position( dRhoTop, Pr(krm1), dRhoBot, Pr(krm1+1) ) +! endif +! if (PoR(k_surface)>=1. .and. krm1= is really ==, when PoR==1 we point to the bottom of the cell +! krm1 = krm1 + 1 +! PoR(k_surface) = PoR(k_surface) - 1. +! endif +! if (real(krm1-lastK_right)+(PoR(k_surface)-lastP_right)<0.) then +! PoR(k_surface) = lastP_right +! krm1 = lastK_right +! endif +! KoR(k_surface) = krm1 +! if (kl <= nk) then +! PoL(k_surface) = 0. +! KoL(k_surface) = kl +! else +! PoL(k_surface) = 1. +! KoL(k_surface) = nk +! endif +! if (kl <= nk) then +! kl = kl + 1 +! else +! reached_bottom = .true. +! searching_right_column = .false. +! searching_left_column = .true. +! endif +! else +! stop 'Else what?' +! endif +! +! lastK_left = KoL(k_surface) ; lastP_left = PoL(k_surface) +! lastK_right = KoR(k_surface) ; lastP_right = PoR(k_surface) +! +! ! Effective thickness +! ! NOTE: This would be better expressed in terms of the layers thicknesses rather +! ! than as differences of position - AJA +! +! ! TODO: GMM, we need to import absolute_position from neutral diffusion. This gives us +! !! the depth of the interface on the left and right side. +! +! if (k_surface>1) then +! hL = absolute_position(nk,ns,Pl,KoL,PoL,k_surface) - absolute_position(nk,ns,Pl,KoL,PoL,k_surface-1) +! hR = absolute_position(nk,ns,Pr,KoR,PoR,k_surface) - absolute_position(nk,ns,Pr,KoR,PoR,k_surface-1) +! if ( hL + hR > 0.) then +! hEff(k_surface-1) = 2. * hL * hR / ( hL + hR ) ! Harmonic mean of layer thicknesses +! else +! hEff(k_surface-1) = 0. +! endif +! endif +! +! enddo neutral_surfaces +!end subroutine find_neutral_surface_positions_continuous + +!> Unit tests for near-boundary horizontal mixing +logical function near_boundary_unit_tests( verbose ) + logical, intent(in) :: verbose !< If true, output additional information for debugging unit tests + + ! Local variables + integer, parameter :: nk = 2 ! Number of layers + integer, parameter :: deg = 1 ! Degree of reconstruction (linear here) + integer, parameter :: method = 1 ! Method used for integrating polynomials + real, dimension(nk) :: phi_L, phi_R ! Tracer values (left and right column) [ nondim m^-3 ] + real, dimension(nk) :: phi_L_avg, phi_R_avg ! Bulk, thickness-weighted tracer averages (left and right column) + real, dimension(nk,deg+1) :: phi_pp_L, phi_pp_R ! Coefficients for the linear pseudo-reconstructions + ! [ nondim m^-3 ] + + real, dimension(nk,2) :: ppoly0_E_L, ppoly0_E_R! Polynomial edge values (left and right) [concentration] + real, dimension(nk) :: h_L, h_R ! Layer thickness (left and right) [m] + real :: khtr_u ! Horizontal diffusivities at U-point [m^2 s^-1] + real :: hbl_L, hbl_R ! Depth of the boundary layer (left and right) [m] + real :: F_bulk ! Total diffusive flux across the U point [nondim s^-1] + real, dimension(nk) :: F_layer ! Diffusive flux within each layer at U-point [nondim s^-1] + real :: h_u, hblt_u ! Thickness at the u-point [m] + real :: khtr_avg ! Thickness-weighted diffusivity at the u-point [m^2 s^-1] + real :: heff ! Harmonic mean of layer thicknesses [m] + real :: inv_heff ! Inverse of the harmonic mean of layer thicknesses [m^[-1] + character(len=120) :: test_name ! Title of the unit test + integer :: k_top ! Index of cell containing top of boundary + real :: zeta_top ! Nondimension position + integer :: k_bot ! Index of cell containing bottom of boundary + real :: zeta_bot ! Nondimension position + real :: area_L,area_R ! Area of grid cell [m^2] + area_L = 1.; area_R = 1. ! Set to unity for all unit tests + + near_boundary_unit_tests = .false. + + ! Unit tests for boundary_k_range + test_name = 'Surface boundary spans the entire top cell' + h_L = (/5.,5./) + call boundary_k_range(SURFACE, nk, h_L, 5., k_top, zeta_top, k_bot, zeta_bot) + near_boundary_unit_tests = test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 0., 1, 1., test_name, verbose) + + test_name = 'Surface boundary spans the entire column' + h_L = (/5.,5./) + call boundary_k_range(SURFACE, nk, h_L, 10., k_top, zeta_top, k_bot, zeta_bot) + near_boundary_unit_tests = test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 0., 2, 1., test_name, verbose) + + test_name = 'Bottom boundary spans the entire bottom cell' + h_L = (/5.,5./) + call boundary_k_range(BOTTOM, nk, h_L, 5., k_top, zeta_top, k_bot, zeta_bot) + near_boundary_unit_tests = test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 2, 0., 2, 1., test_name, verbose) + + test_name = 'Bottom boundary spans the entire column' + h_L = (/5.,5./) + call boundary_k_range(BOTTOM, nk, h_L, 10., k_top, zeta_top, k_bot, zeta_bot) + near_boundary_unit_tests = test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 0., 2, 1., test_name, verbose) + + test_name = 'Surface boundary intersects second layer' + h_L = (/10.,10./) + call boundary_k_range(SURFACE, nk, h_L, 17.5, k_top, zeta_top, k_bot, zeta_bot) + near_boundary_unit_tests = test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 0., 2, 0.75, test_name, verbose) + + test_name = 'Surface boundary intersects first layer' + h_L = (/10.,10./) + call boundary_k_range(SURFACE, nk, h_L, 2.5, k_top, zeta_top, k_bot, zeta_bot) + near_boundary_unit_tests = test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 0., 1, 0.25, test_name, verbose) + + test_name = 'Bottom boundary intersects first layer' + h_L = (/10.,10./) + call boundary_k_range(BOTTOM, nk, h_L, 17.5, k_top, zeta_top, k_bot, zeta_bot) + near_boundary_unit_tests = test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 0.75, 2, 1., test_name, verbose) + + test_name = 'Bottom boundary intersects second layer' + h_L = (/10.,10./) + call boundary_k_range(BOTTOM, nk, h_L, 2.5, k_top, zeta_top, k_bot, zeta_bot) + near_boundary_unit_tests = test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 2, 0.25, 2, 1., test_name, verbose) + + ! All cases in this section have hbl which are equal to the column thicknesses + test_name = 'Equal hbl and same layer thicknesses (gradient from right to left)' + hbl_L = 10; hbl_R = 10 + h_L = (/5.,5./) ; h_R = (/5.,5./) + phi_L = (/0.,0./) ; phi_R = (/1.,1./) + phi_pp_L(1,1) = 0.; phi_pp_L(1,2) = 0. + phi_pp_L(2,1) = 0.; phi_pp_L(2,2) = 0. + phi_pp_R(1,1) = 1.; phi_pp_R(1,2) = 0. + phi_pp_R(2,1) = 1.; phi_pp_R(2,2) = 0. + ppoly0_E_L(1,1) = 0.; ppoly0_E_L(1,2) = 0. + ppoly0_E_L(2,1) = 0.; ppoly0_E_L(2,2) = 0. + ppoly0_E_R(1,1) = 1.; ppoly0_E_R(1,2) = 1. + ppoly0_E_R(2,1) = 1.; ppoly0_E_R(2,2) = 1. + khtr_u = 1. + call fluxes_bulk_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, phi_pp_L, phi_pp_R, & + ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_bulk, F_layer) + near_boundary_unit_tests = test_layer_fluxes( verbose, nk, test_name, F_layer, (/-5.0,-5.0/) ) + + test_name = 'Equal hbl and same layer thicknesses (gradient from left to right)' + hbl_L = 10.; hbl_R = 10. + h_L = (/5.,5./) ; h_R = (/5.,5./) + phi_L = (/1.,1./) ; phi_R = (/0.,0./) + phi_pp_L(1,1) = 1.; phi_pp_L(1,2) = 0. + phi_pp_L(2,1) = 1.; phi_pp_L(2,2) = 0. + phi_pp_R(1,1) = 0.; phi_pp_R(1,2) = 0. + phi_pp_R(2,1) = 0.; phi_pp_R(2,2) = 0. + ppoly0_E_L(1,1) = 1.; ppoly0_E_L(1,2) = 1. + ppoly0_E_L(2,1) = 1.; ppoly0_E_L(2,2) = 1. + ppoly0_E_R(1,1) = 0.; ppoly0_E_R(1,2) = 0. + ppoly0_E_R(2,1) = 0.; ppoly0_E_R(2,2) = 0. + khtr_u = 1. + call fluxes_bulk_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, phi_pp_L, phi_pp_R,& + ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_bulk, F_layer) + near_boundary_unit_tests = test_layer_fluxes( verbose, nk, test_name, F_layer, (/5.0,5.0/) ) + + test_name = 'Equal hbl and same layer thicknesses (no gradient)' + hbl_L = 10; hbl_R = 10 + h_L = (/5.,5./) ; h_R = (/5.,5./) + phi_L = (/1.,1./) ; phi_R = (/1.,1./) + phi_pp_L(1,1) = 1.; phi_pp_L(1,2) = 0. + phi_pp_L(2,1) = 1.; phi_pp_L(2,2) = 0. + phi_pp_R(1,1) = 1.; phi_pp_R(1,2) = 0. + phi_pp_R(2,1) = 1.; phi_pp_R(2,2) = 0. + ppoly0_E_L(1,1) = 1.; ppoly0_E_L(1,2) = 0. + ppoly0_E_L(2,1) = 1.; ppoly0_E_L(2,2) = 0. + ppoly0_E_R(1,1) = 1.; ppoly0_E_R(1,2) = 1. + ppoly0_E_R(2,1) = 1.; ppoly0_E_R(2,2) = 1. + khtr_u = 1. + call fluxes_bulk_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, phi_pp_L, phi_pp_R,& + ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_bulk, F_layer) + near_boundary_unit_tests = test_layer_fluxes( verbose, nk, test_name, F_layer, (/0.0,0.0/) ) + + test_name = 'Equal hbl and different layer thicknesses (gradient right to left)' + hbl_L = 16.; hbl_R = 16. + h_L = (/10.,6./) ; h_R = (/6.,10./) + phi_L = (/0.,0./) ; phi_R = (/1.,1./) + phi_pp_L(1,1) = 0.; phi_pp_L(1,2) = 0. + phi_pp_L(2,1) = 0.; phi_pp_L(2,2) = 0. + phi_pp_R(1,1) = 1.; phi_pp_R(1,2) = 0. + phi_pp_R(2,1) = 1.; phi_pp_R(2,2) = 0. + ppoly0_E_L(1,1) = 0.; ppoly0_E_L(1,2) = 0. + ppoly0_E_L(2,1) = 0.; ppoly0_E_L(2,2) = 0. + ppoly0_E_R(1,1) = 1.; ppoly0_E_R(1,2) = 1. + ppoly0_E_R(2,1) = 1.; ppoly0_E_R(2,2) = 1. + khtr_u = 1. + call fluxes_bulk_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, phi_pp_L, phi_pp_R,& + ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_bulk, F_layer) + near_boundary_unit_tests = test_layer_fluxes( verbose, nk, test_name, F_layer, (/-8.0,-8.0/) ) + + test_name = 'Equal hbl and same layer thicknesses (diagonal tracer values)' + hbl_L = 10.; hbl_R = 10. + h_L = (/5.,5./) ; h_R = (/5.,5./) + phi_L = (/1.,0./) ; phi_R = (/0.,1./) + phi_pp_L(1,1) = 1.; phi_pp_L(1,2) = 0. + phi_pp_L(2,1) = 0.; phi_pp_L(2,2) = 0. + phi_pp_R(1,1) = 0.; phi_pp_R(1,2) = 0. + phi_pp_R(2,1) = 1.; phi_pp_R(2,2) = 0. + ppoly0_E_L(1,1) = 1.; ppoly0_E_L(1,2) = 1. + ppoly0_E_L(2,1) = 0.; ppoly0_E_L(2,2) = 0. + ppoly0_E_R(1,1) = 0.; ppoly0_E_R(1,2) = 0. + ppoly0_E_R(2,1) = 1.; ppoly0_E_R(2,2) = 1. + khtr_u = 1. + call fluxes_bulk_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, phi_pp_L, phi_pp_R,& + ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_bulk, F_layer) + near_boundary_unit_tests = test_layer_fluxes( verbose, nk, test_name, F_layer, (/0.0,0.0/) ) + + test_name = 'Different hbl and different layer thicknesses (gradient from right to left)' + hbl_L = 12; hbl_R = 20 + h_L = (/6.,6./) ; h_R = (/10.,10./) + phi_L = (/0.,0./) ; phi_R = (/1.,1./) + phi_pp_L(1,1) = 0.; phi_pp_L(1,2) = 0. + phi_pp_L(2,1) = 0.; phi_pp_L(2,2) = 0. + phi_pp_R(1,1) = 1.; phi_pp_R(1,2) = 0. + phi_pp_R(2,1) = 1.; phi_pp_R(2,2) = 0. + ppoly0_E_L(1,1) = 0.; ppoly0_E_L(1,2) = 0. + ppoly0_E_L(2,1) = 0.; ppoly0_E_L(2,2) = 0. + ppoly0_E_R(1,1) = 1.; ppoly0_E_R(1,2) = 1. + ppoly0_E_R(2,1) = 1.; ppoly0_E_R(2,2) = 1. + khtr_u = 1. + call fluxes_bulk_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, phi_pp_L, phi_pp_R,& + ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_bulk, F_layer) + near_boundary_unit_tests = test_layer_fluxes( verbose, nk, test_name, F_layer, (/-7.5,-7.5/) ) + + ! Cases where hbl < column thickness (polynomial coefficients specified for pseudo-linear reconstruction) + + test_name = 'hbl < column thickness, hbl same, constant concentration each column' + hbl_L = 2; hbl_R = 2 + h_L = (/1.,2./) ; h_R = (/1.,2./) + phi_L = (/0.,0./) ; phi_R = (/1.,1./) + phi_pp_L(1,1) = 0.; phi_pp_L(1,2) = 0. + phi_pp_L(2,1) = 0.; phi_pp_L(2,2) = 0. + phi_pp_R(1,1) = 1.; phi_pp_R(1,2) = 0. + phi_pp_R(2,1) = 1.; phi_pp_R(2,2) = 0. + khtr_u = 1. + call fluxes_bulk_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, phi_pp_L, phi_pp_R,& + ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_bulk, F_layer) + near_boundary_unit_tests = test_layer_fluxes( verbose, nk, test_name, F_layer, (/-1.,-1./) ) + + test_name = 'hbl < column thickness, hbl same, linear profile right' + hbl_L = 2; hbl_R = 2 + h_L = (/1.,2./) ; h_R = (/1.,2./) + phi_L = (/0.,0./) ; phi_R = (/0.5,2./) + phi_pp_L(1,1) = 0.; phi_pp_L(1,2) = 0. + phi_pp_L(2,1) = 0.; phi_pp_L(2,2) = 0. + phi_pp_R(1,1) = 0.; phi_pp_R(1,2) = 1. + phi_pp_R(2,1) = 1.; phi_pp_R(2,2) = 2. + khtr_u = 1. + ppoly0_E_L(1,1) = 0.; ppoly0_E_L(1,2) = 0. + ppoly0_E_L(2,1) = 0.; ppoly0_E_L(2,2) = 0. + ppoly0_E_R(1,1) = 0.; ppoly0_E_R(1,2) = 1. + ppoly0_E_R(2,1) = 1.; ppoly0_E_R(2,2) = 3. + call fluxes_bulk_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, phi_pp_L, phi_pp_R,& + ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_bulk, F_layer) + near_boundary_unit_tests = test_layer_fluxes( verbose, nk, test_name, F_layer, (/-1.,-1./) ) + + test_name = 'hbl < column thickness, hbl same, linear profile right, khtr=2' + hbl_L = 2; hbl_R = 2 + h_L = (/1.,2./) ; h_R = (/1.,2./) + phi_L = (/0.,0./) ; phi_R = (/0.5,2./) + phi_pp_L(1,1) = 0.; phi_pp_L(1,2) = 0. + phi_pp_L(2,1) = 0.; phi_pp_L(2,2) = 0. + phi_pp_R(1,1) = 0.; phi_pp_R(1,2) = 1. + phi_pp_R(2,1) = 1.; phi_pp_R(2,2) = 2. + khtr_u = 2. + ppoly0_E_L(1,1) = 0.; ppoly0_E_L(1,2) = 0. + ppoly0_E_L(2,1) = 0.; ppoly0_E_L(2,2) = 0. + ppoly0_E_R(1,1) = 0.; ppoly0_E_R(1,2) = 1. + ppoly0_E_R(2,1) = 1.; ppoly0_E_R(2,2) = 3. + call fluxes_bulk_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, phi_pp_L, phi_pp_R,& + ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_bulk, F_layer) + near_boundary_unit_tests = test_layer_fluxes( verbose, nk, test_name, F_layer, (/-2.,-2./) ) + + ! unit tests for layer by layer method + test_name = 'Different hbl and different column thicknesses (gradient from right to left)' + hbl_L = 12; hbl_R = 20 + h_L = (/6.,6./) ; h_R = (/10.,10./) + phi_L = (/0.,0./) ; phi_R = (/1.,1./) + phi_pp_L(1,1) = 0.; phi_pp_L(1,2) = 0. + phi_pp_L(2,1) = 0.; phi_pp_L(2,2) = 0. + phi_pp_R(1,1) = 1.; phi_pp_R(1,2) = 0. + phi_pp_R(2,1) = 1.; phi_pp_R(2,2) = 0. + ppoly0_E_L(1,1) = 0.; ppoly0_E_L(1,2) = 0. + ppoly0_E_L(2,1) = 0.; ppoly0_E_L(2,2) = 0. + ppoly0_E_R(1,1) = 1.; ppoly0_E_R(1,2) = 1. + ppoly0_E_R(2,1) = 1.; ppoly0_E_R(2,2) = 1. + khtr_u = 1. + call fluxes_layer_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, phi_L, phi_R, phi_pp_L, phi_pp_R,& + ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_layer) + near_boundary_unit_tests = test_layer_fluxes( verbose, nk, test_name, F_layer, (/-7.5,0.0/) ) + + test_name = 'Different hbl and different column thicknesses (linear profile right)' + + hbl_L = 15; hbl_R = 6 + h_L = (/10.,10./) ; h_R = (/12.,10./) + phi_L = (/0.,0./) ; phi_R = (/1.,3./) + phi_pp_L(1,1) = 0.; phi_pp_L(1,2) = 0. + phi_pp_L(2,1) = 0.; phi_pp_L(2,2) = 0. + phi_pp_R(1,1) = 0.; phi_pp_R(1,2) = 2. + phi_pp_R(2,1) = 2.; phi_pp_R(2,2) = 2. + ppoly0_E_L(1,1) = 0.; ppoly0_E_L(1,2) = 0. + ppoly0_E_L(2,1) = 0.; ppoly0_E_L(2,2) = 0. + ppoly0_E_R(1,1) = 0.; ppoly0_E_R(1,2) = 2. + ppoly0_E_R(2,1) = 2.; ppoly0_E_R(2,2) = 4. + khtr_u = 1. + call fluxes_layer_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, phi_L, phi_R, phi_pp_L, phi_pp_R,& + ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_layer) + near_boundary_unit_tests = test_layer_fluxes( verbose, nk, test_name, F_layer, (/-3.75,0.0/) ) +end function near_boundary_unit_tests + +!> Returns true if output of near-boundary unit tests does not match correct computed values +!! and conditionally writes results to stream +logical function test_layer_fluxes(verbose, nk, test_name, F_calc, F_ans) + logical, intent(in) :: verbose !< If true, write results to stdout + character(len=80), intent(in) :: test_name !< Brief description of the unit test + integer, intent(in) :: nk !< Number of layers + real, dimension(nk), intent(in) :: F_calc !< Fluxes of the unitless tracer from the algorithm [s^-1] + real, dimension(nk), intent(in) :: F_ans !< Fluxes of the unitless tracer calculated by hand [s^-1] + ! Local variables + integer :: k + integer, parameter :: stdunit = 6 + + test_layer_fluxes = .false. + do k=1,nk + if ( F_calc(k) /= F_ans(k) ) then + test_layer_fluxes = .true. + write(stdunit,*) "UNIT TEST FAILED: ", test_name + write(stdunit,10) k, F_calc(k), F_ans(k) + elseif (verbose) then + write(stdunit,10) k, F_calc(k), F_ans(k) + endif + enddo + +10 format("Layer=",i3," F_calc=",f20.16," F_ans",f20.16) +end function test_layer_fluxes + +!> Return true if output of unit tests for boundary_k_range does not match answers +logical function test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, k_top_ans, zeta_top_ans,& + k_bot_ans, zeta_bot_ans, test_name, verbose) + integer :: k_top !< Index of cell containing top of boundary + real :: zeta_top !< Nondimension position + integer :: k_bot !< Index of cell containing bottom of boundary + real :: zeta_bot !< Nondimension position + integer :: k_top_ans !< Index of cell containing top of boundary + real :: zeta_top_ans !< Nondimension position + integer :: k_bot_ans !< Index of cell containing bottom of boundary + real :: zeta_bot_ans !< Nondimension position + character(len=80) :: test_name !< Name of the unit test + logical :: verbose !< If true always print output + + integer, parameter :: stdunit = 6 + + test_boundary_k_range = k_top .ne. k_top_ans + test_boundary_k_range = test_boundary_k_range .or. (zeta_top .ne. zeta_top_ans) + test_boundary_k_range = test_boundary_k_range .or. (k_bot .ne. k_bot_ans) + test_boundary_k_range = test_boundary_k_range .or. (zeta_bot .ne. zeta_bot_ans) + + if (test_boundary_k_range) write(stdunit,*) "UNIT TEST FAILED: ", test_name + if (test_boundary_k_range .or. verbose) then + write(stdunit,20) "k_top", k_top, "k_top_ans", k_top_ans + write(stdunit,20) "k_bot", k_bot, "k_bot_ans", k_bot_ans + write(stdunit,30) "zeta_top", zeta_top, "zeta_top_ans", zeta_top_ans + write(stdunit,30) "zeta_bot", zeta_bot, "zeta_bot_ans", zeta_bot_ans + endif + + 20 format(A,"=",i3,X,A,"=",i3) + 30 format(A,"=",f20.16,X,A,"=",f20.16) + + +end function test_boundary_k_range + +!> \namespace mom_lbd +!! +!! \section section_LBD The Lateral Boundary Diffusion (LBD) framework +!! +!! The LBD framework accounts for the effects of diabatic mesoscale fluxes +!! within surface and bottom boundary layers. Unlike the equivalent adiabatic +!! fluxes, which is applied along neutral density surfaces, LBD is purely +!! horizontal. +!! +!! The bottom boundary layer fluxes remain to be implemented, although most +!! of the steps needed to do so have already been added and tested. +!! +!! Boundary lateral diffusion can be applied using one of the three methods: +!! +!! * [Method #1: Bulk layer](@ref section_method1) (default); +!! * [Method #2: Along layer](ref section_method2); +!! * [Method #3: Decomposition on to pressure levels](@ref section_method3). +!! +!! A brief summary of these methods is provided below. +!! +!! \subsection section_method1 Bulk layer approach (Method #1) +!! +!! Apply the lateral boundary diffusive fluxes calculated from a 'bulk model' +!! +!! Step #1: get vertical indices containing the boundary layer depth. These are +!! k_top, k_bot, zeta_top, zeta_bot +!! +!! Step #2: compute bulk averages (thickness weighted). phi_L and phi_R +!! +!! Step #3: compute a diffusive bulk flux +!! \f[ F_{bulk} = -(KHTR \times heff) \times (\phi_R - \phi_L), \f] +!! where heff is the harmonic mean of the boundary layer depth in the left and +!! right columns (\f[ HBL_L \f] and \f[ HBL_R \f], respectively). +!! +!! Step #4: limit the tracer flux so that the donor cell, with positive +!! concentration, cannot go negative. If a tracer can go negative (e.g., +!! temperature at high latitudes) it is unclear what limiter should be used. +!! (TODO: ask Bob and Alistair). +!! +!! Step #5: decompose the bulk flux into individual layers and keep track of +!! the remaining flux. The limiter described above is also applied during +!! this step. +!! +!! \subsection section_method2 Along layer approach (Method #2) +!! +!! \subsection section_method3 Decomposition on to pressure levels (Method #3) +!! +!! To be implemented +!! +!! \subsection section_harmonic_mean Harmonic Mean +!! +!! +end module MOM_lateral_boundary_diffusion diff --git a/src/tracer/MOM_neutral_diffusion.F90 b/src/tracer/MOM_neutral_diffusion.F90 index a13eace934..841ab56981 100644 --- a/src/tracer/MOM_neutral_diffusion.F90 +++ b/src/tracer/MOM_neutral_diffusion.F90 @@ -5,19 +5,16 @@ module MOM_neutral_diffusion use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end use MOM_cpu_clock, only : CLOCK_MODULE, CLOCK_ROUTINE +use MOM_domains, only : pass_var use MOM_diag_mediator, only : diag_ctrl, time_type use MOM_diag_mediator, only : post_data, register_diag_field use MOM_EOS, only : EOS_type, EOS_manual_init, calculate_compress, calculate_density_derivs -use MOM_EOS, only : calculate_density_second_derivs +use MOM_EOS, only : calculate_density, calculate_density_second_derivs use MOM_EOS, only : extract_member_EOS, EOS_LINEAR, EOS_TEOS10, EOS_WRIGHT use MOM_error_handler, only : MOM_error, FATAL, WARNING, MOM_mesg, is_root_pe use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_file_parser, only : openParameterBlock, closeParameterBlock use MOM_grid, only : ocean_grid_type -use MOM_neutral_diffusion_aux, only : ndiff_aux_CS_type, set_ndiff_aux_params -use MOM_neutral_diffusion_aux, only : mark_unstable_cells, increment_interface, calc_drho, drho_at_pos -use MOM_neutral_diffusion_aux, only : search_other_column, interpolate_for_nondim_position, refine_nondim_position -use MOM_neutral_diffusion_aux, only : check_neutral_positions use MOM_remapping, only : remapping_CS, initialize_remapping use MOM_remapping, only : extract_member_remapping_CS, build_reconstructions_1d use MOM_remapping, only : average_value_ppoly, remappingSchemesDoc, remappingDefaultScheme @@ -27,7 +24,10 @@ module MOM_neutral_diffusion use polynomial_functions, only : evaluation_polynomial, first_derivative_polynomial use PPM_functions, only : PPM_reconstruction, PPM_boundary_extrapolation use regrid_edge_values, only : edge_values_implicit_h4 - +use MOM_CVMix_KPP, only : KPP_get_BLD, KPP_CS +use MOM_energetic_PBL, only : energetic_PBL_get_MLD, energetic_PBL_CS +use MOM_diabatic_driver, only : diabatic_CS, extract_diabatic_member +use MOM_lateral_boundary_diffusion, only : boundary_k_range, SURFACE, BOTTOM implicit none ; private #include @@ -38,17 +38,18 @@ module MOM_neutral_diffusion !> The control structure for the MOM_neutral_diffusion module type, public :: neutral_diffusion_CS ; private - integer :: nkp1 !< Number of interfaces for a column = nk + 1 - integer :: nsurf !< Number of neutral surfaces - integer :: deg = 2 !< Degree of polynomial used for reconstructions + integer :: nkp1 !< Number of interfaces for a column = nk + 1 + integer :: nsurf !< Number of neutral surfaces + integer :: deg = 2 !< Degree of polynomial used for reconstructions logical :: continuous_reconstruction = .true. !< True if using continuous PPM reconstruction at interfaces - logical :: refine_position = .false. !< If true, iterate to refine the corresponding positions - !! in neighboring columns logical :: debug = .false. !< If true, write verbose debugging messages + logical :: hard_fail_heff !< Bring down the model if a problem with heff is detected integer :: max_iter !< Maximum number of iterations if refine_position is defined - real :: tolerance !< Convergence criterion representing difference from true neutrality + real :: drho_tol !< Convergence criterion representing difference from true neutrality + real :: x_tol !< Convergence criterion for how small an update of the position can be real :: ref_pres !< Reference pressure, negative if using locally referenced neutral density - + logical :: interior_only !< If true, only applies neutral diffusion in the ocean interior. + !! That is, the algorithm will exclude the surface and bottom boundary layers. ! Positions of neutral surfaces in both the u, v directions real, allocatable, dimension(:,:,:) :: uPoL !< Non-dimensional position with left layer uKoL-1, u-point real, allocatable, dimension(:,:,:) :: uPoR !< Non-dimensional position with right layer uKoR-1, u-point @@ -74,21 +75,27 @@ module MOM_neutral_diffusion real, allocatable, dimension(:,:,:) :: Sint !< Interface S [ppt] real, allocatable, dimension(:,:,:) :: Pint !< Interface pressure [Pa] ! Variables needed for discontinuous reconstructions - real, allocatable, dimension(:,:,:,:) :: T_i !< Top edge reconstruction of temperature [degC] - real, allocatable, dimension(:,:,:,:) :: S_i !< Top edge reconstruction of salinity [ppt] - real, allocatable, dimension(:,:,:,:) :: dRdT_i !< dRho/dT [kg m-3 degC-1] at top edge - real, allocatable, dimension(:,:,:,:) :: dRdS_i !< dRho/dS [kg m-3 ppt-1] at top edge + real, allocatable, dimension(:,:,:,:) :: T_i !< Top edge reconstruction of temperature (degC) + real, allocatable, dimension(:,:,:,:) :: S_i !< Top edge reconstruction of salinity (ppt) + real, allocatable, dimension(:,:,:,:) :: P_i !< Interface pressure (Pa) + real, allocatable, dimension(:,:,:,:) :: dRdT_i !< dRho/dT (kg/m3/degC) at top edge + real, allocatable, dimension(:,:,:,:) :: dRdS_i !< dRho/dS (kg/m3/ppt) at top edge integer, allocatable, dimension(:,:) :: ns !< Number of interfacs in a column logical, allocatable, dimension(:,:,:) :: stable_cell !< True if the cell is stably stratified wrt to the next cell type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to !! regulate the timing of diagnostic output. + integer :: neutral_pos_method !< Method to find the position of a neutral surface within the layer + character(len=40) :: delta_rho_form !< Determine which (if any) approximation is made to the + !! equation describing the difference in density + integer :: id_uhEff_2d = -1 !< Diagnostic IDs integer :: id_vhEff_2d = -1 !< Diagnostic IDs real :: C_p !< heat capacity of seawater (J kg-1 K-1) type(EOS_type), pointer :: EOS !< Equation of state parameters type(remapping_CS) :: remap_CS !< Remapping control structure used to create sublayers - type(ndiff_aux_CS_type), pointer :: ndiff_aux_CS !< Store parameters for iteratively finding neutral surface + type(KPP_CS), pointer :: KPP_CSp => NULL() !< KPP control structure needed to get BLD + type(energetic_PBL_CS), pointer :: energetic_PBL_CSp => NULL() !< ePBL control structure needed to get MLD end type neutral_diffusion_CS ! This include declares and sets the variable "version". @@ -98,27 +105,26 @@ module MOM_neutral_diffusion contains !> Read parameters and allocate control structure for neutral_diffusion module. -logical function neutral_diffusion_init(Time, G, param_file, diag, EOS, CS) +logical function neutral_diffusion_init(Time, G, param_file, diag, EOS, diabatic_CSp, CS) type(time_type), target, intent(in) :: Time !< Time structure type(ocean_grid_type), intent(in) :: G !< Grid structure type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure type(param_file_type), intent(in) :: param_file !< Parameter file structure type(EOS_type), target, intent(in) :: EOS !< Equation of state + type(diabatic_CS), pointer :: diabatic_CSp!< KPP control structure needed to get BLD type(neutral_diffusion_CS), pointer :: CS !< Neutral diffusion control structure ! Local variables character(len=256) :: mesg ! Message for error messages. character(len=80) :: string ! Temporary strings logical :: boundary_extrap - ! For refine_pos - integer :: max_iter - real :: drho_tol, xtol, ref_pres if (associated(CS)) then call MOM_error(FATAL, "neutral_diffusion_init called with associated control structure.") return endif + ! Log this module and master switch for turning it on/off call log_version(param_file, mdl, version, & "This module implements neutral diffusion of tracers") @@ -131,7 +137,6 @@ logical function neutral_diffusion_init(Time, G, param_file, diag, EOS, CS) endif allocate(CS) - allocate(CS%ndiff_aux_CS) CS%diag => diag CS%EOS => EOS ! call openParameterBlock(param_file,'NEUTRAL_DIFF') @@ -148,13 +153,16 @@ logical function neutral_diffusion_init(Time, G, param_file, diag, EOS, CS) "the equation of state. If negative (default), local "//& "pressure is used.", & default = -1.) + call get_param(param_file, mdl, "NDIFF_INTERIOR_ONLY", CS%interior_only, & + "If true, only applies neutral diffusion in the ocean interior."//& + "That is, the algorithm will exclude the surface and bottom"//& + "boundary layers.",default = .false.) + ! Initialize and configure remapping if (CS%continuous_reconstruction .eqv. .false.) then call get_param(param_file, mdl, "NDIFF_BOUNDARY_EXTRAP", boundary_extrap, & - "Uses a rootfinding approach to find the position of a "//& - "neutral surface within a layer taking into account the "//& - "nonlinearity of the equation of state and the "//& - "polynomial reconstructions of T/S.", & + "Extrapolate at the top and bottommost cells, otherwise \n"// & + "assume boundaries are piecewise constant", & default=.false.) call get_param(param_file, mdl, "NDIFF_REMAPPING_SCHEME", string, & "This sets the reconstruction scheme used "//& @@ -163,32 +171,52 @@ logical function neutral_diffusion_init(Time, G, param_file, diag, EOS, CS) trim(remappingSchemesDoc), default=remappingDefaultScheme) call initialize_remapping( CS%remap_CS, string, boundary_extrapolation = boundary_extrap ) call extract_member_remapping_CS(CS%remap_CS, degree=CS%deg) - call get_param(param_file, mdl, "NDIFF_REFINE_POSITION", CS%refine_position, & - "Uses a rootfinding approach to find the position of a "//& - "neutral surface within a layer taking into account the "//& - "nonlinearity of the equation of state and the "//& - "polynomial reconstructions of T/S.", & - default=.false.) - if (CS%refine_position) then - call get_param(param_file, mdl, "NDIFF_DRHO_TOL", drho_tol, & - "Sets the convergence criterion for finding the neutral "//& + call get_param(param_file, mdl, "NEUTRAL_POS_METHOD", CS%neutral_pos_method, & + "Method used to find the neutral position \n"// & + "1. Delta_rho varies linearly, find 0 crossing \n"// & + "2. Alpha and beta vary linearly from top to bottom, \n"// & + " Newton's method for neutral position \n"// & + "3. Full nonlinear equation of state, use regula falsi \n"// & + " for neutral position", default=3) + if (CS%neutral_pos_method > 4 .or. CS%neutral_pos_method < 0) then + call MOM_error(FATAL,"Invalid option for NEUTRAL_POS_METHOD") + endif + + call get_param(param_file, mdl, "DELTA_RHO_FORM", CS%delta_rho_form, & + "Determine how the difference in density is calculated \n"// & + " full : Difference of in-situ densities \n"// & + " no_pressure: Calculated from dRdT, dRdS, but no \n"// & + " pressure dependence", & + default="mid_pressure") + if (CS%neutral_pos_method > 1) then + call get_param(param_file, mdl, "NDIFF_DRHO_TOL", CS%drho_tol, & + "Sets the convergence criterion for finding the neutral\n"// & "position within a layer in kg m-3.", & default=1.e-10) - call get_param(param_file, mdl, "NDIFF_X_TOL", xtol, & - "Sets the convergence criterion for a change in nondim "//& + call get_param(param_file, mdl, "NDIFF_X_TOL", CS%x_tol, & + "Sets the convergence criterion for a change in nondim\n"// & "position within a layer.", & default=0.) - call get_param(param_file, mdl, "NDIFF_MAX_ITER", max_iter, & - "The maximum number of iterations to be done before "//& + call get_param(param_file, mdl, "NDIFF_MAX_ITER", CS%max_iter, & + "The maximum number of iterations to be done before \n"// & "exiting the iterative loop to find the neutral surface", & default=10) - call set_ndiff_aux_params(CS%ndiff_aux_CS, max_iter = max_iter, drho_tol = drho_tol, xtol = xtol) endif call get_param(param_file, mdl, "NDIFF_DEBUG", CS%debug, & "Turns on verbose output for discontinuous neutral "//& "diffusion routines.", & default = .false.) - call set_ndiff_aux_params(CS%ndiff_aux_CS, deg=CS%deg, ref_pres = CS%ref_pres, EOS = EOS, debug = CS%debug) + call get_param(param_file, mdl, "HARD_FAIL_HEFF", CS%hard_fail_heff, & + "Bring down the model if a problem with heff is detected",& + default = .true.) + endif + + if (CS%interior_only) then + call extract_diabatic_member(diabatic_CSp, KPP_CSp=CS%KPP_CSp) + call extract_diabatic_member(diabatic_CSp, energetic_PBL_CSp=CS%energetic_PBL_CSp) + if ( .not. ASSOCIATED(CS%energetic_PBL_CSp) .and. .not. ASSOCIATED(CS%KPP_CSp) ) then + call MOM_error(FATAL,"NDIFF_INTERIOR_ONLY is true, but no valid boundary layer scheme was found") + endif endif ! call get_param(param_file, mdl, "KHTR", CS%KhTr, & @@ -203,6 +231,7 @@ logical function neutral_diffusion_init(Time, G, param_file, diag, EOS, CS) CS%nsurf = 4*G%ke ! Discontinuous means that every interface has four connections allocate(CS%T_i(SZI_(G),SZJ_(G),SZK_(G),2)) ; CS%T_i(:,:,:,:) = 0. allocate(CS%S_i(SZI_(G),SZJ_(G),SZK_(G),2)) ; CS%S_i(:,:,:,:) = 0. + allocate(CS%P_i(SZI_(G),SZJ_(G),SZK_(G),2)) ; CS%P_i(:,:,:,:) = 0. allocate(CS%dRdT_i(SZI_(G),SZJ_(G),SZK_(G),2)) ; CS%dRdT_i(:,:,:,:) = 0. allocate(CS%dRdS_i(SZI_(G),SZJ_(G),SZK_(G),2)) ; CS%dRdS_i(:,:,:,:) = 0. allocate(CS%ppoly_coeffs_T(SZI_(G),SZJ_(G),SZK_(G),CS%deg+1)) ; CS%ppoly_coeffs_T(:,:,:,:) = 0. @@ -231,9 +260,10 @@ end function neutral_diffusion_init !> Calculate remapping factors for u/v columns used to map adjoining columns to !! a shared coordinate space. -subroutine neutral_diffusion_calc_coeffs(G, GV, h, T, S, CS) +subroutine neutral_diffusion_calc_coeffs(G, GV, US, h, T, S, CS) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2] real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: T !< Potential temperature [degC] real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: S !< Salinity [ppt] @@ -244,13 +274,36 @@ subroutine neutral_diffusion_calc_coeffs(G, GV, h, T, S, CS) ! Variables used for reconstructions real, dimension(SZK_(G),2) :: ppoly_r_S ! Reconstruction slopes real, dimension(SZI_(G), SZJ_(G)) :: hEff_sum + real, dimension(SZI_(G),SZJ_(G)) :: hbl !< bnd. layer depth [m] integer :: iMethod real, dimension(SZI_(G)) :: ref_pres ! Reference pressure used to calculate alpha/beta + real, dimension(SZI_(G)) :: rho_tmp ! Routine to calculate drho_dp, returns density which is not used real :: h_neglect, h_neglect_edge + integer, dimension(SZI_(G), SZJ_(G)) :: k_top !< Index of the first layer within the boundary + real, dimension(SZI_(G), SZJ_(G)) :: zeta_top !< Distance from the top of a layer to the intersection of the + !! top extent of the boundary layer (0 at top, 1 at bottom) [nondim] + integer, dimension(SZI_(G), SZJ_(G)) :: k_bot !< Index of the last layer within the boundary + real, dimension(SZI_(G), SZJ_(G)) :: zeta_bot !< Distance of the lower layer to the boundary layer depth real :: pa_to_H pa_to_H = 1. / GV%H_to_pa + k_top(:,:) = 1 ; k_bot(:,:) = 1 + zeta_top(:,:) = 0. ; zeta_bot(:,:) = 1. + + ! check if hbl needs to be extracted + if (CS%interior_only) then + hbl(:,:) = 0. + if (ASSOCIATED(CS%KPP_CSp)) call KPP_get_BLD(CS%KPP_CSp, hbl, G) + if (ASSOCIATED(CS%energetic_PBL_CSp)) call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, hbl, G, US) + call pass_var(hbl,G%Domain) + ! get k-indices and zeta + do j=G%jsc-1, G%jec+1 ; do i=G%isc-1,G%iec+1 + call boundary_k_range(SURFACE, G%ke, h(i,j,:), hbl(i,j), k_top(i,j), zeta_top(i,j), k_bot(i,j), zeta_bot(i,j)) + enddo; enddo + ! TODO: add similar code for BOTTOM boundary layer + endif + !### Try replacing both of these with GV%H_subroundoff if (GV%Boussinesq) then h_neglect = GV%m_to_H*1.0e-30 ; h_neglect_edge = GV%m_to_H*1.0e-10 @@ -281,6 +334,19 @@ subroutine neutral_diffusion_calc_coeffs(G, GV, h, T, S, CS) CS%Pint(i,j,k+1) = CS%Pint(i,j,k) + h(i,j,k)*GV%H_to_Pa enddo ; enddo ; enddo + ! Pressures at the interfaces, this is redundant as P_i(k,1) = P_i(k-1,2) however retain tis + ! for now ensure consitency of indexing for diiscontinuous reconstructions + if (.not. CS%continuous_reconstruction) then + do j=G%jsc-1, G%jec+1 ; do i=G%isc-1,G%iec+1 + CS%P_i(i,j,1,1) = 0. + CS%P_i(i,j,1,2) = h(i,j,1)*GV%H_to_Pa + enddo ; enddo + do k=2,G%ke ; do j=G%jsc-1, G%jec+1 ; do i=G%isc-1,G%iec+1 + CS%P_i(i,j,k,1) = CS%P_i(i,j,k-1,2) + CS%P_i(i,j,k,2) = CS%P_i(i,j,k-1,2) + h(i,j,k)*GV%H_to_Pa + enddo ; enddo ; enddo + endif + do j = G%jsc-1, G%jec+1 ! Interpolate state to interface do i = G%isc-1, G%iec+1 @@ -292,6 +358,14 @@ subroutine neutral_diffusion_calc_coeffs(G, GV, h, T, S, CS) CS%T_i(i,j,:,:), ppoly_r_S, iMethod, h_neglect, h_neglect_edge ) call build_reconstructions_1d( CS%remap_CS, G%ke, h(i,j,:), S(i,j,:), CS%ppoly_coeffs_S(i,j,:,:), & CS%S_i(i,j,:,:), ppoly_r_S, iMethod, h_neglect, h_neglect_edge ) + ! In the current ALE formulation, interface values are not exactly at the 0. or 1. of the + ! polynomial reconstructions + do k=1,G%ke + CS%T_i(i,j,k,1) = evaluation_polynomial( CS%ppoly_coeffs_T(i,j,k,:), CS%deg+1, 0. ) + CS%T_i(i,j,k,2) = evaluation_polynomial( CS%ppoly_coeffs_T(i,j,k,:), CS%deg+1, 1. ) + CS%S_i(i,j,k,1) = evaluation_polynomial( CS%ppoly_coeffs_S(i,j,k,:), CS%deg+1, 0. ) + CS%S_i(i,j,k,2) = evaluation_polynomial( CS%ppoly_coeffs_S(i,j,k,:), CS%deg+1, 1. ) + enddo endif enddo @@ -305,11 +379,13 @@ subroutine neutral_diffusion_calc_coeffs(G, GV, h, T, S, CS) else ! Discontinuous reconstruction do k = 1, G%ke if (CS%ref_pres<0) ref_pres(:) = CS%Pint(:,j,k) + ! Calculate derivatives for the top interface call calculate_density_derivs(CS%T_i(:,j,k,1), CS%S_i(:,j,k,1), ref_pres, & CS%dRdT_i(:,j,k,1), CS%dRdS_i(:,j,k,1), G%isc-1, G%iec-G%isc+3, CS%EOS) if (CS%ref_pres<0) then ref_pres(:) = CS%Pint(:,j,k+1) endif + ! Calcualte derivatives at the bottom interface call calculate_density_derivs(CS%T_i(:,j,k,2), CS%S_i(:,j,k,2), ref_pres, & CS%dRdT_i(:,j,k,2), CS%dRdS_i(:,j,k,2), G%isc-1, G%iec-G%isc+3, CS%EOS) enddo @@ -318,8 +394,14 @@ subroutine neutral_diffusion_calc_coeffs(G, GV, h, T, S, CS) if (.not. CS%continuous_reconstruction) then do j = G%jsc-1, G%jec+1 ; do i = G%isc-1, G%iec+1 - call mark_unstable_cells( G%ke, CS%dRdT_i(i,j,:,:), CS%dRdS_i(i,j,:,:), CS%T_i(i,j,:,:), CS%S_i(i,j,:,:), & - CS%stable_cell(i,j,:), CS%ns(i,j) ) + call mark_unstable_cells( CS, G%ke, CS%T_i(i,j,:,:), CS%S_i(i,j,:,:), CS%P_i(i,j,:,:), CS%stable_cell(i,j,:) ) + if (CS%interior_only) then + if (.not. CS%stable_cell(i,j,k_bot(i,j))) zeta_bot(i,j) = -1. + ! set values in the surface and bottom boundary layer to false. + do k = 1, k_bot(i,j)-1 + CS%stable_cell(i,j,k) = .false. + enddo + endif enddo ; enddo endif @@ -341,16 +423,16 @@ subroutine neutral_diffusion_calc_coeffs(G, GV, h, T, S, CS) call find_neutral_surface_positions_continuous(G%ke, & CS%Pint(i,j,:), CS%Tint(i,j,:), CS%Sint(i,j,:), CS%dRdT(i,j,:), CS%dRdS(i,j,:), & CS%Pint(i+1,j,:), CS%Tint(i+1,j,:), CS%Sint(i+1,j,:), CS%dRdT(i+1,j,:), CS%dRdS(i+1,j,:), & - CS%uPoL(I,j,:), CS%uPoR(I,j,:), CS%uKoL(I,j,:), CS%uKoR(I,j,:), CS%uhEff(I,j,:) ) + CS%uPoL(I,j,:), CS%uPoR(I,j,:), CS%uKoL(I,j,:), CS%uKoR(I,j,:), CS%uhEff(I,j,:), & + k_bot(I,j), k_bot(I+1,j), 1.-zeta_bot(I,j), 1.-zeta_bot(I+1,j)) else - call find_neutral_surface_positions_discontinuous(CS, G%ke, CS%ns(i,j)+CS%ns(i+1,j), & - CS%Pint(i,j,:), h(i,j,:), CS%T_i(i,j,:,:), CS%S_i(i,j,:,:), & - CS%dRdT_i(i,j,:,:), CS%dRdS_i(i,j,:,:), CS%stable_cell(i,j,:), & - CS%Pint(i+1,j,:), h(i+1,j,:), CS%T_i(i+1,j,:,:), CS%S_i(i+1,j,:,:), & - CS%dRdT_i(i+1,j,:,:), CS%dRdS_i(i+1,j,:,:), CS%stable_cell(i+1,j,:), & - CS%uPoL(I,j,:), CS%uPoR(I,j,:), CS%uKoL(I,j,:), CS%uKoR(I,j,:), CS%uhEff(I,j,:), & - CS%ppoly_coeffs_T(i,j,:,:), CS%ppoly_coeffs_S(i,j,:,:), & - CS%ppoly_coeffs_T(i+1,j,:,:), CS%ppoly_coeffs_S(i+1,j,:,:)) + call find_neutral_surface_positions_discontinuous(CS, G%ke, & + CS%P_i(i,j,:,:), h(i,j,:), CS%T_i(i,j,:,:), CS%S_i(i,j,:,:), CS%ppoly_coeffs_T(i,j,:,:), & + CS%ppoly_coeffs_S(i,j,:,:),CS%stable_cell(i,j,:), & + CS%P_i(i+1,j,:,:), h(i+1,j,:), CS%T_i(i+1,j,:,:), CS%S_i(i+1,j,:,:), CS%ppoly_coeffs_T(i+1,j,:,:), & + CS%ppoly_coeffs_S(i+1,j,:,:), CS%stable_cell(i+1,j,:), & + CS%uPoL(I,j,:), CS%uPoR(I,j,:), CS%uKoL(I,j,:), CS%uKoR(I,j,:), CS%uhEff(I,j,:), & + hard_fail_heff = CS%hard_fail_heff) endif endif enddo ; enddo @@ -359,28 +441,27 @@ subroutine neutral_diffusion_calc_coeffs(G, GV, h, T, S, CS) do J = G%jsc-1, G%jec ; do i = G%isc, G%iec if (G%mask2dCv(i,J) > 0.) then if (CS%continuous_reconstruction) then - call find_neutral_surface_positions_continuous(G%ke, & + call find_neutral_surface_positions_continuous(G%ke, & CS%Pint(i,j,:), CS%Tint(i,j,:), CS%Sint(i,j,:), CS%dRdT(i,j,:), CS%dRdS(i,j,:), & CS%Pint(i,j+1,:), CS%Tint(i,j+1,:), CS%Sint(i,j+1,:), CS%dRdT(i,j+1,:), CS%dRdS(i,j+1,:), & - CS%vPoL(i,J,:), CS%vPoR(i,J,:), CS%vKoL(i,J,:), CS%vKoR(i,J,:), CS%vhEff(i,J,:) ) + CS%vPoL(i,J,:), CS%vPoR(i,J,:), CS%vKoL(i,J,:), CS%vKoR(i,J,:), CS%vhEff(i,J,:), & + k_bot(i,J), k_bot(i,J+1), 1.-zeta_bot(i,J), 1.-zeta_bot(i,J+1)) else - call find_neutral_surface_positions_discontinuous(CS, G%ke, CS%ns(i,j)+CS%ns(i,j+1), & - CS%Pint(i,j,:), h(i,j,:), CS%T_i(i,j,:,:), CS%S_i(i,j,:,:), & - CS%dRdT_i(i,j,:,:), CS%dRdS_i(i,j,:,:), CS%stable_cell(i,j,:), & - CS%Pint(i,j+1,:), h(i,j+1,:), CS%T_i(i,j+1,:,:), CS%S_i(i,j+1,:,:), & - CS%dRdT_i(i,j+1,:,:), CS%dRdS_i(i,j+1,:,:), CS%stable_cell(i,j+1,:), & - CS%vPoL(I,j,:), CS%vPoR(I,j,:), CS%vKoL(I,j,:), CS%vKoR(I,j,:), CS%vhEff(I,j,:), & - CS%ppoly_coeffs_T(i,j,:,:), CS%ppoly_coeffs_S(i,j,:,:), & - CS%ppoly_coeffs_T(i,j+1,:,:), CS%ppoly_coeffs_S(i,j+1,:,:)) - + call find_neutral_surface_positions_discontinuous(CS, G%ke, & + CS%P_i(i,j,:,:), h(i,j,:), CS%T_i(i,j,:,:), CS%S_i(i,j,:,:), CS%ppoly_coeffs_T(i,j,:,:), & + CS%ppoly_coeffs_S(i,j,:,:),CS%stable_cell(i,j,:), & + CS%P_i(i,j+1,:,:), h(i,j+1,:), CS%T_i(i,j+1,:,:), CS%S_i(i,j+1,:,:), CS%ppoly_coeffs_T(i,j+1,:,:), & + CS%ppoly_coeffs_S(i,j+1,:,:), CS%stable_cell(i,j+1,:), & + CS%vPoL(I,j,:), CS%vPoR(I,j,:), CS%vKoL(I,j,:), CS%vKoR(I,j,:), CS%vhEff(I,j,:), & + hard_fail_heff = CS%hard_fail_heff) endif endif enddo ; enddo ! Continuous reconstructions calculate hEff as the difference between the pressures of the ! neutral surfaces which need to be reconverted to thickness units. The discontinuous version - ! calculates hEff from the fraction of the nondimensional fraction of the layer occupied by - ! the... (Please finish this thought. -RWH) + ! calculates hEff from the fraction of the nondimensional fraction of the layer spanned by + ! adjacent neutral surfaces. if (CS%continuous_reconstruction) then do k = 1, CS%nsurf-1 ; do j = G%jsc, G%jec ; do I = G%isc-1, G%iec if (G%mask2dCu(I,j) > 0.) CS%uhEff(I,j,k) = CS%uhEff(I,j,k) * pa_to_H @@ -498,6 +579,11 @@ subroutine neutral_diffusion(G, GV, h, Coef_x, Coef_y, dt, Reg, US, CS) do k = 1, GV%ke tracer%t(i,j,k) = tracer%t(i,j,k) + dTracer(k) * & ( G%IareaT(i,j) / ( h(i,j,k) + GV%H_subroundoff ) ) +! if (tracer%t(i,j,k) < 0.) then +! do ks = 1,CS%nsurf-1 +! print *, uFlx(I,j,ks), uFlx(I-1,j,ks), vFlx(i,J,ks), vFlx(i,J-1,ks) +! enddo +! endif enddo if (tracer%id_dfxy_conc > 0 .or. tracer%id_dfxy_cont > 0 .or. tracer%id_dfxy_cont_2d > 0 ) then @@ -809,7 +895,7 @@ end function fvlsq_slope !> Returns positions within left/right columns of combined interfaces using continuous reconstructions of T/S subroutine find_neutral_surface_positions_continuous(nk, Pl, Tl, Sl, dRdTl, dRdSl, Pr, Tr, Sr, & - dRdTr, dRdSr, PoL, PoR, KoL, KoR, hEff) + dRdTr, dRdSr, PoL, PoR, KoL, KoR, hEff, bl_kl, bl_kr, bl_zl, bl_zr) integer, intent(in) :: nk !< Number of levels real, dimension(nk+1), intent(in) :: Pl !< Left-column interface pressure [Pa] real, dimension(nk+1), intent(in) :: Tl !< Left-column interface potential temperature [degC] @@ -828,6 +914,10 @@ subroutine find_neutral_surface_positions_continuous(nk, Pl, Tl, Sl, dRdTl, dRdS integer, dimension(2*nk+2), intent(inout) :: KoL !< Index of first left interface above neutral surface integer, dimension(2*nk+2), intent(inout) :: KoR !< Index of first right interface above neutral surface real, dimension(2*nk+1), intent(inout) :: hEff !< Effective thickness between two neutral surfaces [Pa] + integer, optional, intent(in) :: bl_kl !< Layer index of the boundary layer (left) + integer, optional, intent(in) :: bl_kr !< Layer index of the boundary layer (right) + real, optional, intent(in) :: bl_zl !< Nondimensional position of the boundary layer (left) + real, optional, intent(in) :: bl_zr !< Nondimensional position of the boundary layer (right) ! Local variables integer :: ns ! Number of neutral surfaces @@ -842,13 +932,22 @@ subroutine find_neutral_surface_positions_continuous(nk, Pl, Tl, Sl, dRdTl, dRdS real :: dRho, dRhoTop, dRhoBot, hL, hR integer :: lastK_left, lastK_right real :: lastP_left, lastP_right + logical :: interior_limit ns = 2*nk+2 + ! Initialize variables for the search - kr = 1 ; lastK_right = 1 ; lastP_right = 0. - kl = 1 ; lastK_left = 1 ; lastP_left = 0. + kr = 1 ; + kl = 1 ; + lastP_right = 0. + lastP_left = 0. + lastK_right = 1 + lastK_left = 1 reached_bottom = .false. + ! Check to see if we should limit the diffusion to the interior + interior_limit = PRESENT(bl_kl) .and. PRESENT(bl_kr) .and. PRESENT(bl_zr) .and. PRESENT(bl_zl) + ! Loop over each neutral surface, working from top to bottom neutral_surfaces: do k_surface = 1, ns klm1 = max(kl-1, 1) @@ -967,10 +1066,23 @@ subroutine find_neutral_surface_positions_continuous(nk, Pl, Tl, Sl, dRdTl, dRdS else stop 'Else what?' endif + if (interior_limit) then + if (KoL(k_surface)<=bl_kl) then + KoL(k_surface) = bl_kl + if (PoL(k_surface) Returns the non-dimensional position between Pneg and Ppos where the +!! interpolated density difference equals zero. +!! The result is always bounded to be between 0 and 1. +real function interpolate_for_nondim_position(dRhoNeg, Pneg, dRhoPos, Ppos) + real, intent(in) :: dRhoNeg !< Negative density difference + real, intent(in) :: Pneg !< Position of negative density difference + real, intent(in) :: dRhoPos !< Positive density difference + real, intent(in) :: Ppos !< Position of positive density difference + + if (PposdRhoPos) then + write(0,*) 'dRhoNeg, Pneg, dRhoPos, Ppos=',dRhoNeg, Pneg, dRhoPos, Ppos + elseif (dRhoNeg>dRhoPos) then + stop 'interpolate_for_nondim_position: Houston, we have a problem! dRhoNeg>dRhoPos' + endif + if (Ppos<=Pneg) then ! Handle vanished or inverted layers + interpolate_for_nondim_position = 0.5 + elseif ( dRhoPos - dRhoNeg > 0. ) then + interpolate_for_nondim_position = min( 1., max( 0., -dRhoNeg / ( dRhoPos - dRhoNeg ) ) ) + elseif ( dRhoPos - dRhoNeg == 0) then + if (dRhoNeg>0.) then + interpolate_for_nondim_position = 0. + elseif (dRhoNeg<0.) then + interpolate_for_nondim_position = 1. + else ! dRhoPos = dRhoNeg = 0 + interpolate_for_nondim_position = 0.5 + endif + else ! dRhoPos - dRhoNeg < 0 + interpolate_for_nondim_position = 0.5 + endif + if ( interpolate_for_nondim_position < 0. ) & + stop 'interpolate_for_nondim_position: Houston, we have a problem! Pint < Pneg' + if ( interpolate_for_nondim_position > 1. ) & + stop 'interpolate_for_nondim_position: Houston, we have a problem! Pint > Ppos' +end function interpolate_for_nondim_position !> Higher order version of find_neutral_surface_positions. Returns positions within left/right columns -!! of combined interfaces using intracell reconstructions of T/S -subroutine find_neutral_surface_positions_discontinuous(CS, nk, ns, Pres_l, hcol_l, Tl, Sl, & - dRdT_l, dRdS_l, stable_l, Pres_r, hcol_r, Tr, Sr, dRdT_r, dRdS_r, stable_r, & - PoL, PoR, KoL, KoR, hEff, ppoly_T_l, ppoly_S_l, ppoly_T_r, ppoly_S_r) - type(neutral_diffusion_CS), intent(inout) :: CS !< Neutral diffusion control structure - integer, intent(in) :: nk !< Number of levels - integer, intent(in) :: ns !< Number of neutral surfaces - real, dimension(nk+1), intent(in) :: Pres_l !< Left-column interface pressure [Pa] - real, dimension(nk), intent(in) :: hcol_l !< Left-column layer thicknesses - real, dimension(nk,2), intent(in) :: Tl !< Left-column top interface potential temperature [degC] - real, dimension(nk,2), intent(in) :: Sl !< Left-column top interface salinity [ppt] - real, dimension(nk,2), intent(in) :: dRdT_l !< Left-column, top interface dRho/dT [kg m-3 degC-1] - real, dimension(nk,2), intent(in) :: dRdS_l !< Left-column, top interface dRho/dS [kg m-3 ppt-1] - logical, dimension(nk), intent(in) :: stable_l !< Left-column, top interface is stable - real, dimension(nk+1), intent(in) :: Pres_r !< Right-column interface pressure [Pa] - real, dimension(nk), intent(in) :: hcol_r !< Left-column layer thicknesses - real, dimension(nk,2), intent(in) :: Tr !< Right-column top interface potential temperature [degC] - real, dimension(nk,2), intent(in) :: Sr !< Right-column top interface salinity [ppt] - real, dimension(nk,2), intent(in) :: dRdT_r !< Right-column, top interface dRho/dT [kg m-3 degC-1] - real, dimension(nk,2), intent(in) :: dRdS_r !< Right-column, top interface dRho/dS [kg m-3 ppt-1] - logical, dimension(nk), intent(in) :: stable_r !< Right-column, top interface is stable - real, dimension(4*nk), intent(inout) :: PoL !< Fractional position of neutral surface within - !! layer KoL of left column - real, dimension(4*nk), intent(inout) :: PoR !< Fractional position of neutral surface within - !! layer KoR of right column - integer, dimension(4*nk), intent(inout) :: KoL !< Index of first left interface above neutral surface - integer, dimension(4*nk), intent(inout) :: KoR !< Index of first right interface above neutral surface - real, dimension(4*nk-1), intent(inout) :: hEff !< Effective thickness between two neutral surfaces [Pa] - real, dimension(nk,CS%deg+1), & - optional, intent(in) :: ppoly_T_l !< Left-column coefficients of T reconstruction - real, dimension(nk,CS%deg+1), & - optional, intent(in) :: ppoly_S_l !< Left-column coefficients of S reconstruction - real, dimension(nk,CS%deg+1), & - optional, intent(in) :: ppoly_T_r !< Right-column coefficients of T reconstruction - real, dimension(nk,CS%deg+1), & - optional, intent(in) :: ppoly_S_r !< Right-column coefficients of S reconstruction - +!! of combined interfaces using intracell reconstructions of T/S. Note that the polynomial reconstrcutions +!! of T and S are optional to aid with unit testing, but will always be passed otherwise +subroutine find_neutral_surface_positions_discontinuous(CS, nk, Pres_l, hcol_l, Tl, Sl, ppoly_T_l, ppoly_S_l, stable_l,& + Pres_r, hcol_r, Tr, Sr, ppoly_T_r, ppoly_S_r, stable_r,& + PoL, PoR, KoL, KoR, hEff, zeta_bot_L, zeta_bot_R, & + k_bot_L, k_bot_R, hard_fail_heff) + + type(neutral_diffusion_CS), intent(inout) :: CS !< Neutral diffusion control structure + integer, intent(in) :: nk !< Number of levels + real, dimension(nk,2), intent(in) :: Pres_l !< Left-column interface pressure (Pa) + real, dimension(nk), intent(in) :: hcol_l !< Left-column layer thicknesses + real, dimension(nk,2), intent(in) :: Tl !< Left-column top interface potential temperature (degC) + real, dimension(nk,2), intent(in) :: Sl !< Left-column top interface salinity (ppt) + real, dimension(:,:), intent(in) :: ppoly_T_l !< Left-column coefficients of T reconstruction + real, dimension(:,:), intent(in) :: ppoly_S_l !< Left-column coefficients of S reconstruction + logical, dimension(nk), intent(in) :: stable_l !< Left-column, top interface dRho/dS (kg/m3/ppt) + real, dimension(nk,2), intent(in) :: Pres_r !< Right-column interface pressure (Pa) + real, dimension(nk), intent(in) :: hcol_r !< Left-column layer thicknesses + real, dimension(nk,2), intent(in) :: Tr !< Right-column top interface potential temperature (degC) + real, dimension(nk,2), intent(in) :: Sr !< Right-column top interface salinity (ppt) + real, dimension(:,:), intent(in) :: ppoly_T_r !< Right-column coefficients of T reconstruction + real, dimension(:,:), intent(in) :: ppoly_S_r !< Right-column coefficients of S reconstruction + logical, dimension(nk), intent(in) :: stable_r !< Left-column, top interface dRho/dS (kg/m3/ppt) + real, dimension(4*nk), intent(inout) :: PoL !< Fractional position of neutral surface within + !! layer KoL of left column + real, dimension(4*nk), intent(inout) :: PoR !< Fractional position of neutral surface within + !! layer KoR of right column + integer, dimension(4*nk), intent(inout) :: KoL !< Index of first left interface above neutral surface + integer, dimension(4*nk), intent(inout) :: KoR !< Index of first right interface above neutral surface + real, dimension(4*nk-1), intent(inout) :: hEff !< Effective thickness between two neutral surfaces (Pa) + real, optional, intent(in) :: zeta_bot_L!< Non-dimensional distance to where the boundary layer + !! intersetcs the cell (left) [nondim] + real, optional, intent(in) :: zeta_bot_R!< Non-dimensional distance to where the boundary layer + !! intersetcs the cell (right) [nondim] + + integer, optional, intent(in) :: k_bot_L !< k-index for the boundary layer (left) [nondim] + integer, optional, intent(in) :: k_bot_R !< k-index for the boundary layer (right) [nondim] + logical, optional, intent(in) :: hard_fail_heff !< If true (default) bring down the model if the + !! neutral surfaces ever cross [logical] ! Local variables + integer :: ns ! Number of neutral surfaces integer :: k_surface ! Index of neutral surface integer :: kl_left, kl_right ! Index of layers on the left/right integer :: ki_left, ki_right ! Index of interfaces on the left/right @@ -1034,235 +1185,614 @@ subroutine find_neutral_surface_positions_discontinuous(CS, nk, ns, Pres_l, hcol logical :: searching_right_column ! True if searching for the position of a left interface in the right column logical :: reached_bottom ! True if one of the bottom-most interfaces has been used as the target logical :: search_layer - integer :: k, kl_left_0, kl_right_0 + logical :: fail_heff ! By default, real :: dRho, dRhoTop, dRhoBot, hL, hR - integer :: lastK_left, lastK_right + real :: z0, pos + real :: dRdT_from_top, dRdS_from_top ! Density derivatives at the searched from interface + real :: dRdT_from_bot, dRdS_from_bot ! Density derivatives at the searched from interface + real :: dRdT_to_top, dRdS_to_top ! Density derivatives at the interfaces being searched + real :: dRdT_to_bot, dRdS_to_bot ! Density derivatives at the interfaces being searched + real :: T_ref, S_ref, P_ref, P_top, P_bot real :: lastP_left, lastP_right - real :: min_bound - real :: T_other, S_other, P_other, dRdT_other, dRdS_other - logical, dimension(nk) :: top_connected_l, top_connected_r - logical, dimension(nk) :: bot_connected_l, bot_connected_r - - top_connected_l(:) = .false. ; top_connected_r(:) = .false. - bot_connected_l(:) = .false. ; bot_connected_r(:) = .false. - - ! Check to make sure that polynomial reconstructions were passed if refine_pos defined) - if (CS%refine_position) then - if (.not. ( present(ppoly_T_l) .and. present(ppoly_S_l) .and. & - present(ppoly_T_r) .and. present(ppoly_S_r) ) ) & - call MOM_error(FATAL, "fine_neutral_surface_positions_discontinuous: refine_pos is requested, but " //& - "polynomial coefficients not available for T and S") - endif - do k = 1,nk - if (stable_l(k)) then - kl_left = k - kl_left_0 = k - exit - endif - enddo - do k = 1,nk - if (stable_r(k)) then - kl_right = k - kl_right_0 = k - exit - endif - enddo - + integer :: k_init_L, k_init_R ! Starting indices layers for left and right + real :: p_init_L, p_init_R ! Starting positions for left and right ! Initialize variables for the search - ki_right = 1 ; lastK_right = 1 ; lastP_right = -1. - ki_left = 1 ; lastK_left = 1 ; lastP_left = -1. - + ns = 4*nk + ki_right = 1 + ki_left = 1 + kl_left = 1 + kl_right = 1 + lastP_left = 0. + lastP_right = 0. reached_bottom = .false. searching_left_column = .false. searching_right_column = .false. + fail_heff = .true. + if (PRESENT(hard_fail_heff)) fail_heff = hard_fail_heff + + if (PRESENT(k_bot_L) .and. PRESENT(k_bot_R) .and. PRESENT(zeta_bot_L) .and. PRESENT(zeta_bot_R)) then + k_init_L = k_bot_L; k_init_R = k_bot_R + p_init_L = zeta_bot_L; p_init_R = zeta_bot_R + lastP_left = zeta_bot_L; lastP_right = zeta_bot_R + kl_left = k_bot_L; kl_right = k_bot_R + else + k_init_L = 1 ; k_init_R = 1 + p_init_L = 0. ; p_init_R = 0. + endif ! Loop over each neutral surface, working from top to bottom neutral_surfaces: do k_surface = 1, ns - ! Potential density difference, rho(kr) - rho(kl) - dRho = 0.5 * ( ( dRdT_r(kl_right,ki_right) + dRdT_l(kl_left,ki_left) ) * & - ( Tr(kl_right,ki_right) - Tl(kl_left,ki_left) ) & - + ( dRdS_r(kl_right,ki_right) + dRdS_l(kl_left,ki_left) ) * & - ( Sr(kl_right,ki_right) - Sl(kl_left,ki_left) ) ) - if (CS%debug) write(*,'(A,I2,A,E12.4,A,I2,A,I2,A,I2,A,I2)') "k_surface=",k_surface," dRho=",dRho, & - "kl_left=",kl_left, " ki_left=",ki_left," kl_right=",kl_right, " ki_right=",ki_right - ! Which column has the lighter surface for the current indexes, kr and kl - if (.not. reached_bottom) then - if (dRho < 0.) then - searching_left_column = .true. - searching_right_column = .false. - elseif (dRho > 0.) then - searching_right_column = .true. - searching_left_column = .false. - else ! dRho == 0. - if ( ( kl_left == kl_left_0) .and. ( kl_right == kl_right_0 ) .and. & - (ki_left + ki_right == 2) ) then ! Still at surface - searching_left_column = .true. - searching_right_column = .false. - else ! Not the surface so we simply change direction - searching_left_column = .not. searching_left_column - searching_right_column = .not. searching_right_column - endif - endif - endif - if (searching_left_column) then - ! delta_rho is referenced to the right interface T, S, and P - if (CS%ref_pres>=0.) then - P_other = CS%ref_pres + if (k_surface == ns) then + PoL(k_surface) = 1. + PoR(k_surface) = 1. + KoL(k_surface) = nk + KoR(k_surface) = nk + ! If the layers are unstable, then simply point the surface to the previous location + elseif (.not. stable_l(kl_left)) then + if (k_surface > 1) then + PoL(k_surface) = ki_left - 1 ! Top interface is at position = 0., Bottom is at position = 1 + KoL(k_surface) = kl_left + PoR(k_surface) = PoR(k_surface-1) + KoR(k_surface) = KoR(k_surface-1) else - if (ki_right == 1) P_other = Pres_r(kl_right) - if (ki_right == 2) P_other = Pres_r(kl_right+1) + PoR(k_surface) = p_init_R + KoR(k_surface) = k_init_R + PoL(k_surface) = p_init_L + KoL(k_Surface) = k_init_L endif - T_other = Tr(kl_right, ki_right) - S_other = Sr(kl_right, ki_right) - dRdT_other = dRdT_r(kl_right, ki_right) - dRdS_other = dRdS_r(kl_right, ki_right) - if (CS%refine_position .and. (lastK_left == kl_left)) then - call drho_at_pos(CS%ndiff_aux_CS, T_other, S_other, dRdT_other, dRdS_other, Pres_l(kl_left), & - Pres_l(kl_left+1), ppoly_T_l(kl_left,:), ppoly_S_l(kl_left,:), lastP_left, dRhoTop) + call increment_interface(nk, kl_left, ki_left, reached_bottom, searching_left_column, searching_right_column) + searching_left_column = .true. + searching_right_column = .false. + elseif (.not. stable_r(kl_right)) then ! Check the right layer for stability + if (k_surface > 1) then + PoR(k_surface) = ki_right - 1 ! Top interface is at position = 0., Bottom is at position = 1 + KoR(k_surface) = kl_right + PoL(k_surface) = PoL(k_surface-1) + KoL(k_surface) = KoL(k_surface-1) else - dRhoTop = calc_drho(Tl(kl_left,1), Sl(kl_left,1), dRdT_l(kl_left,1), dRdS_l(kl_left,1), T_other, S_other, & - dRdT_other, dRdS_other) - endif - ! Potential density difference, rho(kl) - rho(kl_right,ki_right) (will be positive) - dRhoBot = calc_drho(Tl(kl_left,2), Sl(kl_left,2), dRdT_l(kl_left,2), dRdS_l(kl_left,2), & - T_other, S_other, dRdT_other, dRdS_other) - if (CS%debug) then - write(*,'(A,I2,A,E12.4,A,E12.4,A,E12.4)') "Searching left layer ", kl_left, & - " dRhoTop=", dRhoTop, " dRhoBot=", dRhoBot - write(*,'(A,I2,X,I2)') "Searching from right: ", kl_right, ki_right - write(*,*) "Temp/Salt Reference: ", T_other, S_other - write(*,*) "Temp/Salt Top L: ", Tl(kl_left,1), Sl(kl_left,1) - write(*,*) "Temp/Salt Bot L: ", Tl(kl_left,2), Sl(kl_left,2) + PoR(k_surface) = 0. + KoR(k_surface) = 1 + PoL(k_surface) = 0. + KoL(k_surface) = 1 endif - - ! Set the position within the starting column - PoR(k_surface) = REAL(ki_right-1) - KoR(k_surface) = kl_right - - ! Set position within the searched column - call search_other_column(dRhoTop, dRhoBot, Pres_l(kl_left), Pres_l(kl_left+1), & - lastP_left, lastK_left, kl_left, kl_left_0, ki_left, & - top_connected_l, bot_connected_l, PoL(k_surface), KoL(k_surface), search_layer) - - if ( CS%refine_position .and. search_layer ) then - min_bound = 0. - if (k_surface > 1) then - if ( KoL(k_surface) == KoL(k_surface-1) ) min_bound = PoL(k_surface-1) + call increment_interface(nk, kl_right, ki_right, reached_bottom, searching_right_column, searching_left_column) + searching_left_column = .false. + searching_right_column = .true. + else ! Layers are stable so need to figure out whether we need to search right or left + ! For convenience, the left column uses the searched "from" interface variables, and the right column + ! uses the searched 'to'. These will get reset in subsequent calc_delta_rho calls + + call calc_delta_rho_and_derivs(CS, & + Tr(kl_right, ki_right), Sr(kl_right, ki_right), Pres_r(kl_right,ki_right), & + Tl(kl_left, ki_left), Sl(kl_left, ki_left) , Pres_l(kl_left,ki_left), & + dRho) + if (CS%debug) write(*,'(A,I2,A,E12.4,A,I2,A,I2,A,I2,A,I2)') "k_surface=",k_surface," dRho=",dRho, & + "kl_left=",kl_left, " ki_left=",ki_left," kl_right=",kl_right, " ki_right=",ki_right + ! Which column has the lighter surface for the current indexes, kr and kl + if (.not. reached_bottom) then + if (dRho < 0.) then + searching_left_column = .true. + searching_right_column = .false. + elseif (dRho > 0.) then + searching_left_column = .false. + searching_right_column = .true. + else ! dRho == 0. + if ( ( kl_left + kl_right == 2 ) .and. (ki_left + ki_right == 2) ) then ! Still at surface + searching_left_column = .true. + searching_right_column = .false. + else ! Not the surface so we simply change direction + searching_left_column = .not. searching_left_column + searching_right_column = .not. searching_right_column + endif endif - PoL(k_surface) = refine_nondim_position( CS%ndiff_aux_CS, T_other, S_other, dRdT_other, dRdS_other, & - Pres_l(kl_left), Pres_l(kl_left+1), ppoly_T_l(kl_left,:), ppoly_S_l(kl_left,:), & - dRhoTop, dRhoBot, min_bound ) endif - if (PoL(k_surface) == 0.) top_connected_l(KoL(k_surface)) = .true. - if (PoL(k_surface) == 1.) bot_connected_l(KoL(k_surface)) = .true. - call increment_interface(nk, kl_right, ki_right, stable_r, reached_bottom, & - searching_right_column, searching_left_column) - - elseif (searching_right_column) then - if (CS%ref_pres>=0.) then - P_other = CS%ref_pres - else - if (ki_left == 1) P_other = Pres_l(kl_left) - if (ki_left == 2) P_other = Pres_l(kl_left+1) - endif - T_other = Tl(kl_left, ki_left) - S_other = Sl(kl_left, ki_left) - dRdT_other = dRdT_l(kl_left, ki_left) - dRdS_other = dRdS_l(kl_left, ki_left) - ! Interpolate for the neutral surface position within the right column, layer krm1 - ! Potential density difference, rho(kr-1) - rho(kl) (should be negative) - - if (CS%refine_position .and. (lastK_right == kl_right)) then - call drho_at_pos(CS%ndiff_aux_CS, T_other, S_other, dRdT_other, dRdS_other, Pres_r(kl_right), & - Pres_l(kl_right+1), ppoly_T_r(kl_right,:), ppoly_S_r(kl_right,:), lastP_right, dRhoTop) - else - dRhoTop = calc_drho(Tr(kl_right,1), Sr(kl_right,1), dRdT_r(kl_right,1), dRdS_r(kl_right,1), & - T_other, S_other, dRdT_other, dRdS_other) - endif - dRhoBot = calc_drho(Tr(kl_right,2), Sr(kl_right,2), dRdT_r(kl_right,2), dRdS_r(kl_right,2), & - T_other, S_other, dRdT_other, dRdS_other) - if (CS%debug) then - write(*,'(A,I2,A,E12.4,A,E12.4,A,E12.4)') "Searching right layer ", kl_right, & - " dRhoTop=", dRhoTop, " dRhoBot=", dRhoBot - write(*,'(A,I2,X,I2)') "Searching from left: ", kl_left, ki_left - write(*,*) "Temp/Salt Reference: ", T_other, S_other - write(*,*) "Temp/Salt Top R: ", Tr(kl_right,1), Sr(kl_right,1) - write(*,*) "Temp/Salt Bot R: ", Tr(kl_right,2), Sr(kl_right,2) - endif - ! Set the position within the starting column - PoL(k_surface) = REAL(ki_left-1) - KoL(k_surface) = kl_left - - ! Set position within the searched column - call search_other_column(dRhoTop, dRhoBot, Pres_r(kl_right), Pres_r(kl_right+1), lastP_right, lastK_right, & - kl_right, kl_right_0, ki_right, top_connected_r, bot_connected_r, PoR(k_surface), & - KoR(k_surface), search_layer) - if ( CS%refine_position .and. search_layer) then - min_bound = 0. - if (k_surface > 1) then - if ( KoR(k_surface) == KoR(k_surface-1) ) min_bound = PoR(k_surface-1) + if (searching_left_column) then + ! Position of the right interface is known and all quantities are fixed + PoR(k_surface) = ki_right - 1. + KoR(k_surface) = kl_right + PoL(k_surface) = search_other_column(CS, k_surface, lastP_left, & + Tr(kl_right, ki_right), Sr(kl_right, ki_right), Pres_r(kl_right, ki_right), & + Tl(kl_left,1), Sl(kl_left,1), Pres_l(kl_left,1), & + Tl(kl_left,2), Sl(kl_left,2), Pres_l(kl_left,2), & + ppoly_T_l(kl_left,:), ppoly_S_l(kl_left,:)) + KoL(k_surface) = kl_left + + if (CS%debug) then + write(*,'(A,I2)') "Searching left layer ", kl_left + write(*,'(A,I2,X,I2)') "Searching from right: ", kl_right, ki_right + write(*,*) "Temp/Salt Reference: ", Tr(kl_right,ki_right), Sr(kl_right,ki_right) + write(*,*) "Temp/Salt Top L: ", Tl(kl_left,1), Sl(kl_left,1) + write(*,*) "Temp/Salt Bot L: ", Tl(kl_left,2), Sl(kl_left,2) + endif + call increment_interface(nk, kl_right, ki_right, reached_bottom, searching_right_column, searching_left_column) + lastP_left = PoL(k_surface) + ! If the right layer increments, then we need to reset the last position on the right + if ( kl_right == (KoR(k_surface) + 1) ) lastP_right = 0. + + elseif (searching_right_column) then + ! Position of the right interface is known and all quantities are fixed + PoL(k_surface) = ki_left - 1. + KoL(k_surface) = kl_left + PoR(k_surface) = search_other_column(CS, k_surface, lastP_right, & + Tl(kl_left, ki_left), Sl(kl_left, ki_left), Pres_l(kl_left, ki_left), & + Tr(kl_right,1), Sr(kl_right,1), Pres_r(kl_right,1), & + Tr(kl_right,2), Sr(kl_right,2), Pres_r(kl_right,2), & + ppoly_T_r(kl_right,:), ppoly_S_r(kl_right,:)) + KoR(k_surface) = kl_right + + if (CS%debug) then + write(*,'(A,I2)') "Searching right layer ", kl_right + write(*,'(A,I2,X,I2)') "Searching from left: ", kl_left, ki_left + write(*,*) "Temp/Salt Reference: ", Tl(kl_left,ki_left), Sl(kl_left,ki_left) + write(*,*) "Temp/Salt Top L: ", Tr(kl_right,1), Sr(kl_right,1) + write(*,*) "Temp/Salt Bot L: ", Tr(kl_right,2), Sr(kl_right,2) endif - PoR(k_surface) = refine_nondim_position(CS%ndiff_aux_CS, T_other, S_other, dRdT_other, dRdS_other, & - Pres_r(kl_right), Pres_r(kl_right+1), ppoly_T_r(kl_right,:), ppoly_S_r(kl_right,:), & - dRhoTop, dRhoBot, min_bound ) + call increment_interface(nk, kl_left, ki_left, reached_bottom, searching_left_column, searching_right_column) + lastP_right = PoR(k_surface) + ! If the right layer increments, then we need to reset the last position on the right + if ( kl_left == (KoL(k_surface) + 1) ) lastP_left = 0. + else + stop 'Else what?' endif - if (PoR(k_surface) == 0.) top_connected_r(KoR(k_surface)) = .true. - if (PoR(k_surface) == 1.) bot_connected_r(KoR(k_surface)) = .true. - call increment_interface(nk, kl_left, ki_left, stable_l, reached_bottom, & - searching_left_column, searching_right_column) - - else - stop 'Else what?' + if (CS%debug) write(*,'(A,I3,A,ES16.6,A,I2,A,ES16.6)') "KoL:", KoL(k_surface), " PoL:", PoL(k_surface), & + " KoR:", KoR(k_surface), " PoR:", PoR(k_surface) endif - lastK_left = KoL(k_surface) ; lastP_left = PoL(k_surface) - lastK_right = KoR(k_surface) ; lastP_right = PoR(k_surface) - - if (CS%debug) write(*,'(A,I3,A,ES16.6,A,I2,A,ES16.6)') "KoL:", KoL(k_surface), " PoL:", PoL(k_surface), & - " KoR:", KoR(k_surface), " PoR:", PoR(k_surface) ! Effective thickness if (k_surface>1) then - ! This is useful as a check to make sure that positions are monotonically increasing - hL = absolute_position(nk,ns,Pres_l,KoL,PoL,k_surface) - absolute_position(nk,ns,Pres_l,KoL,PoL,k_surface-1) - hR = absolute_position(nk,ns,Pres_r,KoR,PoR,k_surface) - absolute_position(nk,ns,Pres_r,KoR,PoR,k_surface-1) - ! In the case of a layer being unstably stratified, may get a negative thickness. Set the previous position - ! to the current location - if ( hL<0. .or. hR<0. ) then - hEff(k_surface-1) = 0. - call MOM_error(WARNING, "hL or hR is negative") - elseif ( hL > 0. .and. hR > 0.) then + if ( KoL(k_surface) == KoL(k_surface-1) .and. KoR(k_surface) == KoR(k_surface-1) ) then hL = (PoL(k_surface) - PoL(k_surface-1))*hcol_l(KoL(k_surface)) hR = (PoR(k_surface) - PoR(k_surface-1))*hcol_r(KoR(k_surface)) - hEff(k_surface-1) = 2. * ( (hL * hR) / ( hL + hR ) )! Harmonic mean + if (hL < 0. .or. hR < 0.) then + if (fail_heff) then + call MOM_error(FATAL,"Negative thicknesses in neutral diffusion") + else + if (searching_left_column) then + PoL(k_surface) = PoL(k_surface-1) + KoL(k_surface) = KoL(k_surface-1) + elseif (searching_right_column) then + PoR(k_surface) = PoR(k_surface-1) + KoR(k_surface) = KoR(k_surface-1) + endif + endif + elseif ( hL + hR == 0. ) then + hEff(k_surface-1) = 0. + else + hEff(k_surface-1) = 2. * ( (hL * hR) / ( hL + hR ) )! Harmonic mean + if ( KoL(k_surface) /= KoL(k_surface-1) ) then + call MOM_error(FATAL,"Neutral sublayer spans multiple layers") + endif + if ( KoR(k_surface) /= KoR(k_surface-1) ) then + call MOM_error(FATAL,"Neutral sublayer spans multiple layers") + endif + endif else hEff(k_surface-1) = 0. endif endif enddo neutral_surfaces - if (CS%debug) then - write (*,*) "==========Start Neutral Surfaces==========" - do k = 1,ns-1 - if (hEff(k)>0.) then - kl_left = KoL(k) - kl_right = KoR(k) - write (*,'(A,I3,X,ES16.6,X,I3,X,ES16.6)') "Top surface KoL, PoL, KoR, PoR: ", kl_left, PoL(k), kl_right, PoR(k) - call check_neutral_positions(CS%ndiff_aux_CS, Pres_l(kl_left), Pres_l(kl_left+1), Pres_r(kl_right), & - Pres_r(kl_right+1), PoL(k), PoR(k), ppoly_T_l(kl_left,:), ppoly_T_r(kl_right,:), & - ppoly_S_l(kl_left,:), ppoly_S_r(kl_right,:)) - kl_left = KoL(k+1) - kl_right = KoR(k+1) - write (*,'(A,I3,X,ES16.6,X,I3,X,ES16.6)') "Bot surface KoL, PoL, KoR, PoR: ", kl_left, PoL(k+1), kl_right, PoR(k) - call check_neutral_positions(CS%ndiff_aux_CS, Pres_l(kl_left), Pres_l(kl_left+1), Pres_r(kl_right), & - Pres_r(kl_right+1), PoL(k), PoR(k), ppoly_T_l(kl_left,:), ppoly_T_r(kl_right,:), & - ppoly_S_l(kl_left,:), ppoly_S_r(kl_right,:)) +end subroutine find_neutral_surface_positions_discontinuous + +!> Sweep down through the column and mark as stable if the bottom interface of a cell is denser than the top +subroutine mark_unstable_cells(CS, nk, T, S, P, stable_cell) + type(neutral_diffusion_CS), intent(inout) :: CS !< Neutral diffusion control structure + integer, intent(in) :: nk !< Number of levels in a column + real, dimension(nk,2), intent(in) :: T !< Temperature at interfaces + real, dimension(nk,2), intent(in) :: S !< Salinity at interfaces + real, dimension(nk,2), intent(in) :: P !< Pressure at interfaces + logical, dimension(nk), intent( out) :: stable_cell !< True if this cell is unstably stratified + + integer :: k, first_stable, prev_stable + real :: delta_rho + + do k = 1,nk + call calc_delta_rho_and_derivs( CS, T(k,2), S(k,2), max(P(k,2),CS%ref_pres), & + T(k,1), S(k,1), max(P(k,1),CS%ref_pres), delta_rho ) + stable_cell(k) = delta_rho > 0. + enddo +end subroutine mark_unstable_cells + +!> Searches the "other" (searched) column for the position of the neutral surface +real function search_other_column(CS, ksurf, pos_last, T_from, S_from, P_from, T_top, S_top, P_top, & + T_bot, S_bot, P_bot, T_poly, S_poly ) result(pos) + type(neutral_diffusion_CS), intent(in ) :: CS !< Neutral diffusion control structure + integer, intent(in ) :: ksurf !< Current index of neutral surface + real, intent(in ) :: pos_last !< Last position within the current layer, used as the lower + !! bound in the rootfinding algorithm + real, intent(in ) :: T_from !< Temperature at the searched from interface + real, intent(in ) :: S_from !< Salinity at the searched from interface + real, intent(in ) :: P_from !< Pressure at the searched from interface + real, intent(in ) :: T_top !< Temperature at the searched to top interface + real, intent(in ) :: S_top !< Salinity at the searched to top interface + real, intent(in ) :: P_top !< Pressure at the searched to top interface + real, intent(in ) :: T_bot !< Temperature at the searched to bottom interface + real, intent(in ) :: S_bot !< Salinity at the searched to bottom interface + real, intent(in ) :: P_bot !< Pressure at the searched to bottom interface + real, dimension(:), intent(in ) :: T_poly !< Temperature polynomial reconstruction coefficients + real, dimension(:), intent(in ) :: S_poly !< Salinity polynomial reconstruction coefficients + ! Local variables + real :: dRhotop, dRhobot + real :: dRdT_top, dRdS_top, dRdT_bot, dRdS_bot + real :: dRdT_from, dRdS_from + real :: P_mid + + ! Calculate the differencei in density at the tops or the bottom + if (CS%neutral_pos_method == 1 .or. CS%neutral_pos_method == 3) then + call calc_delta_rho_and_derivs(CS, T_top, S_top, P_top, T_from, S_from, P_from, dRhoTop) + call calc_delta_rho_and_derivs(CS, T_bot, S_bot, P_bot, T_from, S_from, P_from, dRhoBot) + elseif (CS%neutral_pos_method == 2) then + call calc_delta_rho_and_derivs(CS, T_top, S_top, P_top, T_from, S_from, P_from, dRhoTop, & + dRdT_top, dRdS_top, dRdT_from, dRdS_from) + call calc_delta_rho_and_derivs(CS, T_bot, S_bot, P_bot, T_from, S_from, P_from, dRhoBot, & + dRdT_bot, dRdS_bot, dRdT_from, dRdS_from) + endif + + ! Handle all the special cases EXCEPT if it connects within the layer + if ( (dRhoTop > 0.) .or. (ksurf == 1) ) then ! First interface or lighter than anything in layer + pos = pos_last + if (CS%debug) print *, "Lighter" + elseif ( dRhoTop > dRhoBot ) then ! Unstably stratified + pos = 1. + if (CS%debug) print *, "Unstable" + elseif ( dRhoTop < 0. .and. dRhoBot < 0.) then ! Denser than anything in layer + pos = 1. + if (CS%debug) print *, "Denser" + elseif ( dRhoTop == 0. .and. dRhoBot == 0. ) then ! Perfectly unstratified + pos = 1. + if (CS%debug) print *, "Unstratified" + elseif ( dRhoBot == 0. ) then ! Matches perfectly at the Top + pos = 1. + if (CS%debug) print *, "Bottom" + elseif ( dRhoTop == 0. ) then ! Matches perfectly at the Bottom + pos = pos_last + if (CS%debug) print *, "Top" + else ! Neutral surface within layer + pos = -1 + if (CS%debug) print *, "Interpolate" + endif + + ! Can safely return if position is >= 0 otherwise will need to find the position within the layer + if (pos>=0) return + + if (CS%neutral_pos_method==1) then + pos = interpolate_for_nondim_position( dRhoTop, P_top, dRhoBot, P_bot ) + ! For the 'Linear' case of finding the neutral position, the fromerence pressure to use is the average + ! of the midpoint of the layer being searched and the interface being searched from + elseif (CS%neutral_pos_method == 2) then + pos = find_neutral_pos_linear( CS, pos_last, T_from, S_from, P_from, dRdT_from, dRdS_from, & + P_top, dRdT_top, dRdS_top, & + P_bot, dRdT_bot, dRdS_bot, T_poly, S_poly ) + elseif (CS%neutral_pos_method == 3) then + pos = find_neutral_pos_full( CS, pos_last, T_from, S_from, P_from, P_top, P_bot, T_poly, S_poly) + endif + +end function search_other_column + +!> Increments the interface which was just connected and also set flags if the bottom is reached +subroutine increment_interface(nk, kl, ki, reached_bottom, searching_this_column, searching_other_column) + integer, intent(in ) :: nk !< Number of vertical levels + integer, intent(inout) :: kl !< Current layer (potentially updated) + integer, intent(inout) :: ki !< Current interface + logical, intent(inout) :: reached_bottom !< Updated when kl == nk and ki == 2 + logical, intent(inout) :: searching_this_column !< Updated when kl == nk and ki == 2 + logical, intent(inout) :: searching_other_column !< Updated when kl == nk and ki == 2 + integer :: k + + reached_bottom = .false. + if (ki == 2) then ! At the bottom interface + if ((ki == 2) .and. (kl < nk) ) then ! Not at the bottom so just go to the next layer + kl = kl+1 + ki = 1 + elseif ((kl == nk) .and. (ki==2)) then + reached_bottom = .true. + searching_this_column = .false. + searching_other_column = .true. + endif + elseif (ki==1) then ! At the top interface + ki = 2 ! Next interface is same layer, but bottom interface + else + call MOM_error(FATAL,"Unanticipated eventuality in increment_interface") + endif +end subroutine increment_interface + +!> Search a layer to find where delta_rho = 0 based on a linear interpolation of alpha and beta of the top and bottom +!! being searched and polynomial reconstructions of T and S. Compressibility is not needed because either, we are +!! assuming incompressibility in the equation of state for this module or alpha and beta are calculated having been +!! displaced to the average pressures of the two pressures We need Newton's method because the T and S reconstructions +!! make delta_rho a polynomial function of z if using PPM or higher. If Newton's method would search fall out of the +!! interval [0,1], a bisection step would be taken instead. Also this linearization of alpha, beta means that second +!! derivatives of the EOS are not needed. Note that delta in variable names below refers to horizontal differences and +!! 'd' refers to vertical differences +function find_neutral_pos_linear( CS, z0, T_ref, S_ref, P_ref, dRdT_ref, dRdS_ref, & + P_top, dRdT_top, dRdS_top, & + P_bot, dRdT_bot, dRdS_bot, ppoly_T, ppoly_S ) result( z ) + type(neutral_diffusion_CS),intent(in) :: CS !< Control structure with parameters for this module + real, intent(in) :: z0 !< Lower bound of position, also serves as the initial guess + real, intent(in) :: T_ref !< Temperature at the searched from interface + real, intent(in) :: S_ref !< Salinity at the searched from interface + real, intent(in) :: P_ref !< Pressure at the searched from interface + real, intent(in) :: dRdT_ref !< dRho/dT at the searched from interface + real, intent(in) :: dRdS_ref !< dRho/dS at the searched from interface + real, intent(in) :: P_top !< Pressure at top of layer being searched + real, intent(in) :: dRdT_top !< dRho/dT at top of layer being searched + real, intent(in) :: dRdS_top !< dRho/dS at top of layer being searched + real, intent(in) :: P_bot !< Pressure at bottom of layer being searched + real, intent(in) :: dRdT_bot !< dRho/dT at bottom of layer being searched + real, intent(in) :: dRdS_bot !< dRho/dS at bottom of layer being searched + real, dimension(:), intent(in) :: ppoly_T !< Coefficients of the polynomial reconstruction of T within + !! the layer to be searched. + real, dimension(:), intent(in) :: ppoly_S !< Coefficients of the polynomial reconstruction of T within + !! the layer to be searched. + real :: z !< Position where drho = 0 + ! Local variables + real :: dRdT_diff, dRdS_diff + real :: drho, drho_dz, dRdT_z, dRdS_z, T_z, S_z, deltaT, deltaS, deltaP, dT_dz, dS_dz + real :: drho_min, drho_max, ztest, zmin, zmax, dRdT_sum, dRdS_sum, dz, P_z, dP_dz + real :: a1, a2 + integer :: iter + integer :: nterm + real :: T_top, T_bot, S_top, S_bot + + nterm = SIZE(ppoly_T) + + ! Position independent quantities + dRdT_diff = dRdT_bot - dRdT_top + dRdS_diff = dRdS_bot - dRdS_top + ! Assume a linear increase in pressure from top and bottom of the cell + dP_dz = P_bot - P_top + ! Initial starting drho (used for bisection) + zmin = z0 ! Lower bounding interval + zmax = 1. ! Maximum bounding interval (bottom of layer) + a1 = 1. - zmin + a2 = zmin + T_z = evaluation_polynomial( ppoly_T, nterm, zmin ) + S_z = evaluation_polynomial( ppoly_S, nterm, zmin ) + dRdT_z = a1*dRdT_top + a2*dRdT_bot + dRdS_z = a1*dRdS_top + a2*dRdS_bot + P_z = a1*P_top + a2*P_bot + drho_min = delta_rho_from_derivs(T_z, S_z, P_z, dRdT_z, dRdS_z, & + T_ref, S_ref, P_ref, dRdT_ref, dRdS_ref) + + T_z = evaluation_polynomial( ppoly_T, nterm, 1. ) + S_z = evaluation_polynomial( ppoly_S, nterm, 1. ) + drho_max = delta_rho_from_derivs(T_z, S_z, P_bot, dRdT_bot, dRdS_bot, & + T_ref, S_ref, P_ref, dRdT_ref, dRdS_ref) + + if (drho_min >= 0.) then + z = z0 + return + elseif (drho_max == 0.) then + z = 1. + return + endif + if ( SIGN(1.,drho_min) == SIGN(1.,drho_max) ) then + print *, drho_min, drho_max + call MOM_error(FATAL, "drho_min is the same sign as dhro_max") + endif + + z = z0 + ztest = z0 + do iter = 1, CS%max_iter + ! Calculate quantities at the current nondimensional position + a1 = 1.-z + a2 = z + dRdT_z = a1*dRdT_top + a2*dRdT_bot + dRdS_z = a1*dRdS_top + a2*dRdS_bot + T_z = evaluation_polynomial( ppoly_T, nterm, z ) + S_z = evaluation_polynomial( ppoly_S, nterm, z ) + P_z = a1*P_top + a2*P_bot + deltaT = T_z - T_ref + deltaS = S_z - S_ref + deltaP = P_z - P_ref + dRdT_sum = dRdT_ref + dRdT_z + dRdS_sum = dRdS_ref + dRdS_z + drho = delta_rho_from_derivs(T_z, S_z, P_z, dRdT_z, dRdS_z, & + T_ref, S_ref, P_ref, dRdT_ref, dRdS_ref) + + ! Check for convergence + if (ABS(drho) <= CS%drho_tol) exit + ! Update bisection bracketing intervals + if (drho < 0. .and. drho > drho_min) then + drho_min = drho + zmin = z + elseif (drho > 0. .and. drho < drho_max) then + drho_max = drho + zmax = z + endif + + ! Calculate a Newton step + dT_dz = first_derivative_polynomial( ppoly_T, nterm, z ) + dS_dz = first_derivative_polynomial( ppoly_S, nterm, z ) + drho_dz = 0.5*( (dRdT_diff*deltaT + dRdT_sum*dT_dz) + (dRdS_diff*deltaS + dRdS_sum*dS_dz) ) + + ztest = z - drho/drho_dz + ! Take a bisection if z falls out of [zmin,zmax] + if (ztest < zmin .or. ztest > zmax) then + if ( drho < 0. ) then + ztest = 0.5*(z + zmax) + else + ztest = 0.5*(zmin + z) endif - enddo - write(*,'(A,E16.6)') "Total thickness of sublayers: ", SUM(hEff) - write(*,*) "==========End Neutral Surfaces==========" + endif + + ! Test to ensure we haven't stalled out + if ( abs(z-ztest) <= CS%x_tol ) exit + ! Reset for next iteration + z = ztest + enddo + +end function find_neutral_pos_linear + +!> Use the full equation of state to calculate the difference in locally referenced potential density. The derivatives +!! in this case are not trivial to calculate, so instead we use a regula falsi method +function find_neutral_pos_full( CS, z0, T_ref, S_ref, P_ref, P_top, P_bot, ppoly_T, ppoly_S ) result( z ) + type(neutral_diffusion_CS),intent(in) :: CS !< Control structure with parameters for this module + real, intent(in) :: z0 !< Lower bound of position, also serves as the initial guess + real, intent(in) :: T_ref !< Temperature at the searched from interface + real, intent(in) :: S_ref !< Salinity at the searched from interface + real, intent(in) :: P_ref !< Pressure at the searched from interface + real, intent(in) :: P_top !< Pressure at top of layer being searched + real, intent(in) :: P_bot !< Pressure at bottom of layer being searched + real, dimension(:), intent(in) :: ppoly_T !< Coefficients of the polynomial reconstruction of T within + !! the layer to be searched. + real, dimension(:), intent(in) :: ppoly_S !< Coefficients of the polynomial reconstruction of T within + !! the layer to be searched. + real :: z !< Position where drho = 0 + ! Local variables + integer :: iter + integer :: nterm + + real :: drho_a, drho_b, drho_c + real :: a, b, c, Ta, Tb, Tc, Sa, Sb, Sc, Pa, Pb, Pc + integer :: side + + side = 0 + ! Set the first two evaluation to the endpoints of the interval + b = z0; c = 1 + nterm = SIZE(ppoly_T) + + ! Calculate drho at the minimum bound + Tb = evaluation_polynomial( ppoly_T, nterm, b ) + Sb = evaluation_polynomial( ppoly_S, nterm, b ) + Pb = P_top*(1.-b) + P_bot*b + call calc_delta_rho_and_derivs(CS, Tb, Sb, Pb, T_ref, S_ref, P_ref, drho_b) + + ! Calculate drho at the maximum bound + Tc = evaluation_polynomial( ppoly_T, nterm, 1. ) + Sc = evaluation_polynomial( ppoly_S, nterm, 1. ) + Pc = P_Bot + call calc_delta_rho_and_derivs(CS, Tc, Sc, Pc, T_ref, S_ref, P_ref, drho_c) + + if (drho_b >= 0.) then + z = z0 + return + elseif (drho_c == 0.) then + z = 1. + return + endif + if ( SIGN(1.,drho_b) == SIGN(1.,drho_c) ) then +! print *, drho_b, drho_c +! call MOM_error(WARNING, "drho_b is the same sign as dhro_c") + z = z0 + return endif -end subroutine find_neutral_surface_positions_discontinuous + do iter = 1, CS%max_iter + ! Calculate new position and evaluate if we have converged + a = (drho_b*c - drho_c*b)/(drho_b-drho_c) + Ta = evaluation_polynomial( ppoly_T, nterm, a ) + Sa = evaluation_polynomial( ppoly_S, nterm, a ) + Pa = P_top*(1.-a) + P_bot*a + call calc_delta_rho_and_derivs(CS, Ta, Sa, Pa, T_ref, S_ref, P_ref, drho_a) + if (ABS(drho_a) < CS%drho_tol) then + z = a + return + endif + + if (drho_a*drho_c > 0.) then + if ( ABS(a-c) 0 ) then + if ( ABS(a-b) Calculate the difference in density between two points in a variety of ways +subroutine calc_delta_rho_and_derivs(CS, T1, S1, p1_in, T2, S2, p2_in, drho, & + drdt1_out, drds1_out, drdt2_out, drds2_out ) + type(neutral_diffusion_CS) :: CS !< Neutral diffusion control structure + real, intent(in ) :: T1 !< Temperature at point 1 + real, intent(in ) :: S1 !< Salinity at point 1 + real, intent(in ) :: p1_in !< Pressure at point 1 + real, intent(in ) :: T2 !< Temperature at point 2 + real, intent(in ) :: S2 !< Salinity at point 2 + real, intent(in ) :: p2_in !< Pressure at point 2 + real, intent( out) :: drho !< Difference in density between the two points + real, optional, intent( out) :: dRdT1_out !< drho_dt at point 1 + real, optional, intent( out) :: dRdS1_out !< drho_ds at point 1 + real, optional, intent( out) :: dRdT2_out !< drho_dt at point 2 + real, optional, intent( out) :: dRdS2_out !< drho_ds at point 2 + ! Local variables + real :: rho1, rho2, p1, p2, pmid + real :: drdt1, drdt2, drds1, drds2, drdp1, drdp2, rho_dummy + + ! Use the same reference pressure or the in-situ pressure + if (CS%ref_pres > 0.) then + p1 = CS%ref_pres + p2 = CS%ref_pres + else + p1 = p1_in + p2 = p2_in + endif + + ! Use the full linear equation of state to calculate the difference in density (expensive!) + if (TRIM(CS%delta_rho_form) == 'full') then + pmid = 0.5 * (p1 + p2) + call calculate_density( T1, S1, pmid, rho1, CS%EOS ) + call calculate_density( T2, S2, pmid, rho2, CS%EOS ) + drho = rho1 - rho2 + ! Use the density derivatives at the average of pressures and the differentces int temperature + elseif (TRIM(CS%delta_rho_form) == 'mid_pressure') then + pmid = 0.5 * (p1 + p2) + if (CS%ref_pres>=0) pmid = CS%ref_pres + call calculate_density_derivs(T1, S1, pmid, drdt1, drds1, CS%EOS) + call calculate_density_derivs(T2, S2, pmid, drdt2, drds2, CS%EOS) + drho = delta_rho_from_derivs( T1, S1, P1, drdt1, drds1, T2, S2, P2, drdt2, drds2) + elseif (TRIM(CS%delta_rho_form) == 'local_pressure') then + call calculate_density_derivs(T1, S1, p1, drdt1, drds1, CS%EOS) + call calculate_density_derivs(T2, S2, p2, drdt2, drds2, CS%EOS) + drho = delta_rho_from_derivs( T1, S1, P1, drdt1, drds1, T2, S2, P2, drdt2, drds2) + else + call MOM_error(FATAL, "delta_rho_form is not recognized") + endif + if (PRESENT(drdt1_out)) drdt1_out = drdt1 + if (PRESENT(drds1_out)) drds1_out = drds1 + if (PRESENT(drdt2_out)) drdt2_out = drdt2 + if (PRESENT(drds2_out)) drds2_out = drds2 + +end subroutine calc_delta_rho_and_derivs + +!> Calculate delta rho from derivatives and gradients of properties +!! \f$ \Delta \rho$ = \frac{1}{2}\left[ (\alpha_1 + \alpha_2)*(T_1-T_2) + +!! (\beta_1 + \beta_2)*(S_1-S_2) + +!! (\gamma^{-1}_1 + \gamma%{-1}_2)*(P_1-P_2) \right] \f$ +function delta_rho_from_derivs( T1, S1, P1, dRdT1, dRdS1, & + T2, S2, P2, dRdT2, dRdS2 ) result (drho) + real :: T1 !< Temperature at point 1 + real :: S1 !< Salinity at point 1 + real :: P1 !< Pressure at point 1 + real :: dRdT1 !< Pressure at point 1 + real :: dRdS1 !< Pressure at point 1 + real :: T2 !< Temperature at point 2 + real :: S2 !< Salinity at point 2 + real :: P2 !< Pressure at point 2 + real :: dRdT2 !< Pressure at point 2 + real :: dRdS2 !< Pressure at point 2 + ! Local variables + real :: drho + + drho = 0.5 * ( (dRdT1+dRdT2)*(T1-T2) + (dRdS1+dRdS2)*(S1-S2)) + +end function delta_rho_from_derivs !> Converts non-dimensional position within a layer to absolute position (for debugging) real function absolute_position(n,ns,Pint,Karr,NParr,k_surface) integer, intent(in) :: n !< Number of levels @@ -1731,7 +2261,7 @@ logical function ndiff_unit_tests_continuous(verbose) (/0.,0.,0.,0.,0.,0.,1.,1./), & ! pL (/0.,0.,0.,0.,0.,0.,1.,1./), & ! pR (/0.,10.,0.,10.,0.,10.,0./), & ! hEff - 'Indentical columns with mixed layer') + 'Identical columns with mixed layer') ! Right column with unstable mixed layer call find_neutral_surface_positions_continuous(3, & @@ -1789,14 +2319,15 @@ logical function ndiff_unit_tests_discontinuous(verbose) integer, parameter :: ns = nk*4 real, dimension(nk) :: Sl, Sr, Tl, Tr, hl, hr real, dimension(nk,2) :: TiL, SiL, TiR, SiR - real, dimension(nk+1) :: Pres_l, Pres_R + real, dimension(nk,2) :: Pres_l, Pres_r integer, dimension(ns) :: KoL, KoR real, dimension(ns) :: PoL, PoR real, dimension(ns-1) :: hEff, Flx type(neutral_diffusion_CS) :: CS !< Neutral diffusion control structure type(EOS_type), pointer :: EOS !< Structure for linear equation of state type(remapping_CS), pointer :: remap_CS !< Remapping control structure (PLM) - real, dimension(nk,2) :: poly_T_l, poly_T_r, poly_S, poly_slope ! Linear reconstruction for T + real, dimension(nk,2) :: ppoly_T_l, ppoly_T_r ! Linear reconstruction for T + real, dimension(nk,2) :: ppoly_S_l, ppoly_S_r ! Linear reconstruction for S real, dimension(nk,2) :: dRdT, dRdS logical, dimension(nk) :: stable_l, stable_r integer :: iMethod @@ -1808,180 +2339,226 @@ logical function ndiff_unit_tests_discontinuous(verbose) v = verbose ndiff_unit_tests_discontinuous = .false. ! Normally return false write(*,*) '==== MOM_neutral_diffusion: ndiff_unit_tests_discontinuous =' - +! h_neglect = 1.0e-30 ; h_neglect_edge = 1.0e-10 ! Unit tests for find_neutral_surface_positions_discontinuous ! Salinity is 0 for all these tests - Sl(:) = 0. ; Sr(:) = 0. ; poly_S(:,:) = 0. ; SiL(:,:) = 0. ; SiR(:,:) = 0. - dRdT(:,:) = -1. ; dRdS(:,:) = 0. - + allocate(CS%EOS) + call EOS_manual_init(CS%EOS, form_of_EOS = EOS_LINEAR, dRho_dT = -1., dRho_dS = 0.) + Sl(:) = 0. ; Sr(:) = 0. ; ; SiL(:,:) = 0. ; SiR(:,:) = 0. + ppoly_T_l(:,:) = 0.; ppoly_T_r(:,:) = 0. + ppoly_S_l(:,:) = 0.; ppoly_S_r(:,:) = 0. ! Intialize any control structures needed for unit tests - CS%refine_position = .false. CS%ref_pres = -1. - allocate(remap_CS) - call initialize_remapping( remap_CS, "PLM", boundary_extrapolation = .true. ) - hL = (/10.,10.,10./) ; hR = (/10.,10.,10./) ; Pres_l(1) = 0. ; Pres_r(1) = 0. - do k = 1,nk ; Pres_l(k+1) = Pres_l(k) + hL(k) ; Pres_r(k+1) = Pres_r(k) + hR(k) ; enddo - ! Identical columns - Tl = (/20.,16.,12./) ; Tr = (/20.,16.,12./) - call build_reconstructions_1d( remap_CS, nk, hL, Tl, poly_T_l, TiL, poly_slope, iMethod, h_neglect, h_neglect_edge ) - call build_reconstructions_1d( remap_CS, nk, hR, Tr, poly_T_r, TiR, poly_slope, iMethod, h_neglect, h_neglect_edge ) - call mark_unstable_cells( nk, dRdT, dRdS, Til, Sil, stable_l, ns_l ) - call mark_unstable_cells( nk, dRdT, dRdS, Tir, Sir, stable_r, ns_r ) - call find_neutral_surface_positions_discontinuous(CS, nk, ns_l+ns_r, Pres_l, hL, TiL, SiL, dRdT, dRdS, stable_l, & - Pres_r, hR, TiR, SiR, dRdT, dRdS, stable_r, PoL, PoR, KoL, KoR, hEff) + hL = (/10.,10.,10./) ; hR = (/10.,10.,10./) + Pres_l(1,1) = 0. ; Pres_l(1,2) = hL(1) ; Pres_r(1,1) = 0. ; Pres_r(1,2) = hR(1) + do k = 2,nk + Pres_l(k,1) = Pres_l(k-1,2) + Pres_l(k,2) = Pres_l(k,1) + hL(k) + Pres_r(k,1) = Pres_r(k-1,2) + Pres_r(k,2) = Pres_r(k,1) + hR(k) + enddo + CS%delta_rho_form = 'mid_pressure' + CS%neutral_pos_method = 1 + + TiL(1,:) = (/ 22.00, 18.00 /); TiL(2,:) = (/ 18.00, 14.00 /); TiL(3,:) = (/ 14.00, 10.00 /); + TiR(1,:) = (/ 22.00, 18.00 /); TiR(2,:) = (/ 18.00, 14.00 /); TiR(3,:) = (/ 14.00, 10.00 /); + call mark_unstable_cells( CS, nk, Til, Sil, Pres_l, stable_l ) + call mark_unstable_cells( CS, nk, Tir, Sir, Pres_r, stable_r ) + call find_neutral_surface_positions_discontinuous(CS, nk, Pres_l, hL, TiL, SiL, ppoly_T_l, ppoly_S_l, stable_l, & + Pres_r, hR, TiR, SiR, ppoly_T_r, ppoly_S_r, stable_r, PoL, PoR, KoL, KoR, hEff) ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, KoL, KoR, PoL, PoR, hEff, & - (/1,1,1,1,2,2,2,2,3,3,3,3/), & ! KoL - (/1,1,1,1,2,2,2,2,3,3,3,3/), & ! KoR - (/0.,0.,1.,1.,0.,0.,1.,1.,0.,0.,1.,1./), & ! pL - (/0.,0.,1.,1.,0.,0.,1.,1.,0.,0.,1.,1./), & ! pR - (/0.,10.,0.,0.,0.,10.,0.,0.,0.,10.,0.,0./), & ! hEff - 'Identical columns') - Tl = (/20.,16.,12./) ; Tr = (/18.,14.,10./) - call build_reconstructions_1d( remap_CS, nk, hL, Tl, poly_T_l, TiL, poly_slope, iMethod, h_neglect, h_neglect_edge ) - call build_reconstructions_1d( remap_CS, nk, hR, Tr, poly_T_r, TiR, poly_slope, iMethod, h_neglect, h_neglect_edge ) - call mark_unstable_cells( nk, dRdT, dRdS, Til, Sil, stable_l, ns_l ) - call mark_unstable_cells( nk, dRdT, dRdS, Tir, Sir, stable_r, ns_r ) - call find_neutral_surface_positions_discontinuous(CS, nk, ns_l+ns_r, Pres_l, hL, TiL, SiL, dRdT, dRdS, stable_l, & - Pres_r, hR, TiR, SiR, dRdT, dRdS, stable_r, PoL, PoR, KoL, KoR, hEff) + (/ 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3 /), & ! KoL + (/ 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3 /), & ! KoR + (/ 0.00, 0.00, 1.00, 1.00, 0.00, 0.00, 1.00, 1.00, 0.00, 0.00, 1.00, 1.00 /), & ! PoL + (/ 0.00, 0.00, 1.00, 0.00, 0.00, 0.00, 1.00, 0.00, 0.00, 0.00, 1.00, 1.00 /), & ! PoR + (/ 0.00, 10.00, 0.00, 0.00, 0.00, 10.00, 0.00, 0.00, 0.00, 10.00, 0.00 /), & ! hEff + 'Identical Columns') + + TiL(1,:) = (/ 22.00, 18.00 /); TiL(2,:) = (/ 18.00, 14.00 /); TiL(3,:) = (/ 14.00, 10.00 /); + TiR(1,:) = (/ 20.00, 16.00 /); TiR(2,:) = (/ 16.00, 12.00 /); TiR(3,:) = (/ 12.00, 8.00 /); + call mark_unstable_cells( CS, nk, Til, Sil, Pres_l, stable_l ) + call mark_unstable_cells( CS, nk, Tir, Sir, Pres_r, stable_r ) + call find_neutral_surface_positions_discontinuous(CS, nk, Pres_l, hL, TiL, SiL, ppoly_T_l, ppoly_S_l, stable_l, & + Pres_r, hR, TiR, SiR, ppoly_T_r, ppoly_S_r, stable_r, PoL, PoR, KoL, KoR, hEff) ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, KoL, KoR, PoL, PoR, hEff, & - (/1,1,1,2,2,2,2,3,3,3,3,3/), & ! KoL - (/1,1,1,1,1,2,2,2,2,3,3,3/), & ! KoR - (/0.0, 0.5, 1.0, 0.0, 0.5, 0.5, 1.0, 0.0, 0.5, 0.5, 1.0, 1.0/), & ! pL - (/0.0, 0.0, 0.5, 0.5, 1.0, 0.0, 0.5, 0.5, 1.0, 0.0, 0.5, 1.0/), & ! pR - (/0.0, 5.0, 0.0, 5.0, 0.0, 5.0, 0.0, 5.0, 0.0, 5.0, 0.0/), & ! hEff - 'Right column slightly cooler') - Tl = (/18.,14.,10./) ; Tr = (/20.,16.,12./) - call build_reconstructions_1d( remap_CS, nk, hL, Tl, poly_T_l, TiL, poly_slope, iMethod, h_neglect, h_neglect_edge ) - call build_reconstructions_1d( remap_CS, nk, hR, Tr, poly_T_r, TiR, poly_slope, iMethod, h_neglect, h_neglect_edge ) - call mark_unstable_cells( nk, dRdT, dRdS, Til, Sil, stable_l, ns_l ) - call mark_unstable_cells( nk, dRdT, dRdS, Tir, Sir, stable_r, ns_r ) - call find_neutral_surface_positions_discontinuous(CS, nk, ns_l+ns_r, Pres_l, hL, TiL, SiL, dRdT, dRdS, stable_l, & - Pres_r, hR, TiR, SiR, dRdT, dRdS, stable_r, PoL, PoR, KoL, KoR, hEff) + (/ 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3 /), & ! KoL + (/ 1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3 /), & ! KoR + (/ 0.00, 0.50, 1.00, 0.00, 0.50, 0.50, 1.00, 0.00, 0.50, 0.50, 1.00, 1.00 /), & ! PoL + (/ 0.00, 0.00, 0.50, 0.50, 1.00, 0.00, 0.50, 0.50, 1.00, 0.00, 0.50, 1.00 /), & ! PoR + (/ 0.00, 5.00, 0.00, 5.00, 0.00, 5.00, 0.00, 5.00, 0.00, 5.00, 0.00 /), & ! hEff + 'Right slightly cooler') + + TiL(1,:) = (/ 20.00, 16.00 /); TiL(2,:) = (/ 16.00, 12.00 /); TiL(3,:) = (/ 12.00, 8.00 /); + TiR(1,:) = (/ 22.00, 18.00 /); TiR(2,:) = (/ 18.00, 14.00 /); TiR(3,:) = (/ 14.00, 10.00 /); + call mark_unstable_cells( CS, nk, Til, Sil, Pres_l, stable_l ) + call mark_unstable_cells( CS, nk, Tir, Sir, Pres_r, stable_r ) + call find_neutral_surface_positions_discontinuous(CS, nk, Pres_l, hL, TiL, SiL, ppoly_T_l, ppoly_S_l, stable_l, & + Pres_r, hR, TiR, SiR, ppoly_T_r, ppoly_S_r, stable_r, PoL, PoR, KoL, KoR, hEff) ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, KoL, KoR, PoL, PoR, hEff, & - (/1,1,1,1,1,2,2,2,2,3,3,3/), & ! KoL - (/1,1,1,2,2,2,2,3,3,3,3,3/), & ! KoR - (/0.0, 0.0, 0.5, 0.5, 1.0, 0.0, 0.5, 0.5, 1.0, 0.0, 0.5, 1.0/), & ! pL - (/0.0, 0.5, 1.0, 0.0, 0.5, 0.5, 1.0, 0.0, 0.5, 0.5, 1.0, 1.0/), & ! pR - (/0.0, 5.0, 0.0, 5.0, 0.0, 5.0, 0.0, 5.0, 0.0, 5.0, 0.0/), & ! hEff - 'Left column slightly cooler') - Tl = (/20.,16.,12./) ; Tr = (/14.,10.,6./) - call build_reconstructions_1d( remap_CS, nk, hL, Tl, poly_T_l, TiL, poly_slope, iMethod, h_neglect, h_neglect_edge ) - call build_reconstructions_1d( remap_CS, nk, hR, Tr, poly_T_r, TiR, poly_slope, iMethod, h_neglect, h_neglect_edge ) - call mark_unstable_cells( nk, dRdT, dRdS, Til, Sil, stable_l, ns_l ) - call mark_unstable_cells( nk, dRdT, dRdS, Tir, Sir, stable_r, ns_r ) - call find_neutral_surface_positions_discontinuous(CS, nk, ns_l+ns_r, Pres_l, hL, TiL, SiL, dRdT, dRdS, stable_l, & - Pres_r, hR, TiR, SiR, dRdT, dRdS, stable_r, PoL, PoR, KoL, KoR, hEff) + (/ 1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3 /), & ! KoL + (/ 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3 /), & ! KoR + (/ 0.00, 0.00, 0.50, 0.50, 1.00, 0.00, 0.50, 0.50, 1.00, 0.00, 0.50, 1.00 /), & ! PoL + (/ 0.00, 0.50, 1.00, 0.00, 0.50, 0.50, 1.00, 0.00, 0.50, 0.50, 1.00, 1.00 /), & ! PoR + (/ 0.00, 5.00, 0.00, 5.00, 0.00, 5.00, 0.00, 5.00, 0.00, 5.00, 0.00 /), & ! hEff + 'Left slightly cooler') + + TiL(1,:) = (/ 22.00, 20.00 /); TiL(2,:) = (/ 18.00, 16.00 /); TiL(3,:) = (/ 14.00, 12.00 /); + TiR(1,:) = (/ 32.00, 24.00 /); TiR(2,:) = (/ 22.00, 14.00 /); TiR(3,:) = (/ 12.00, 4.00 /); + call mark_unstable_cells( CS, nk, Til, Sil, Pres_l, stable_l ) + call mark_unstable_cells( CS, nk, Tir, Sir, Pres_r, stable_r ) + call find_neutral_surface_positions_discontinuous(CS, nk, Pres_l, hL, TiL, SiL, ppoly_T_l, ppoly_S_l, stable_l, & + Pres_r, hR, TiR, SiR, ppoly_T_r, ppoly_S_r, stable_r, PoL, PoR, KoL, KoR, hEff) ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, KoL, KoR, PoL, PoR, hEff, & - (/1,1,2,2,2,3,3,3,3,3,3,3/), & ! KoL - (/1,1,1,1,1,1,1,2,2,2,3,3/), & ! KoR - (/0.0, 1.0, 0.0, 0.5, 1.0, 0.0, 0.5, 0.5, 1.0, 1.0, 1.0, 1.0/), & ! pL - (/0.0, 0.0, 0.0, 0.0, 0.5, 0.5, 1.0, 0.0, 0.5, 1.0, 0.0, 1.0/), & ! pR - (/0.0, 0.0, 0.0, 5.0, 0.0, 5.0, 0.0, 5.0, 0.0, 0.0, 0.0/), & ! hEff - 'Right column somewhat cooler') - Tl = (/20.,16.,12./) ; Tr = (/8.,6.,4./) - call build_reconstructions_1d( remap_CS, nk, hL, Tl, poly_T_l, TiL, poly_slope, iMethod, h_neglect, h_neglect_edge ) - call build_reconstructions_1d( remap_CS, nk, hR, Tr, poly_T_r, TiR, poly_slope, iMethod, h_neglect, h_neglect_edge ) - call mark_unstable_cells( nk, dRdT, dRdS, Til, Sil, stable_l, ns_l ) - call mark_unstable_cells( nk, dRdT, dRdS, Tir, Sir, stable_r, ns_r ) - call find_neutral_surface_positions_discontinuous(CS, nk, ns_l+ns_r, Pres_l, hL, TiL, SiL, dRdT, dRdS, stable_l, & - Pres_r, hR, TiR, SiR, dRdT, dRdS, stable_r, PoL, PoR, KoL, KoR, hEff) + (/ 1, 1, 1, 1, 1, 2, 2, 3, 3, 3, 3, 3 /), & ! KoL + (/ 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3 /), & ! KoR + (/ 0.00, 0.00, 0.00, 0.00, 1.00, 0.00, 1.00, 0.00, 0.00, 1.00, 1.00, 1.00 /), & ! PoL + (/ 0.00, 1.00, 0.00, 0.00, 0.25, 0.50, 0.75, 1.00, 0.00, 0.00, 0.00, 1.00 /), & ! PoR + (/ 0.00, 0.00, 0.00, 4.00, 0.00, 4.00, 0.00, 0.00, 0.00, 0.00, 0.00 /), & ! hEff + 'Right more strongly stratified') + + TiL(1,:) = (/ 22.00, 18.00 /); TiL(2,:) = (/ 18.00, 14.00 /); TiL(3,:) = (/ 14.00, 10.00 /); + TiR(1,:) = (/ 14.00, 14.00 /); TiR(2,:) = (/ 14.00, 14.00 /); TiR(3,:) = (/ 12.00, 8.00 /); + call mark_unstable_cells( CS, nk, Til, Sil, Pres_l, stable_l ) + call mark_unstable_cells( CS, nk, Tir, Sir, Pres_r, stable_r ) + call find_neutral_surface_positions_discontinuous(CS, nk, Pres_l, hL, TiL, SiL, ppoly_T_l, ppoly_S_l, stable_l, & + Pres_r, hR, TiR, SiR, ppoly_T_r, ppoly_S_r, stable_r, PoL, PoR, KoL, KoR, hEff) ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, KoL, KoR, PoL, PoR, hEff, & - (/1,1,2,2,3,3,3,3,3,3,3,3/), & ! KoL - (/1,1,1,1,1,1,1,1,2,2,3,3/), & ! KoR - (/0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0/), & ! pL - (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0/), & ! pR - (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/), & ! hEff - 'Right column much cooler') - Tl = (/14.,14.,10./) ; Tr = (/14.,14.,10./) - call build_reconstructions_1d( remap_CS, nk, hL, Tl, poly_T_l, TiL, poly_slope, iMethod, h_neglect, h_neglect_edge ) - call build_reconstructions_1d( remap_CS, nk, hR, Tr, poly_T_r, TiR, poly_slope, iMethod, h_neglect, h_neglect_edge ) - call mark_unstable_cells( nk, dRdT, dRdS, Til, Sil, stable_l, ns_l ) - call mark_unstable_cells( nk, dRdT, dRdS, Tir, Sir, stable_r, ns_r ) - call find_neutral_surface_positions_discontinuous(CS, nk, ns_l+ns_r, Pres_l, hL, TiL, SiL, dRdT, dRdS, stable_l, & - Pres_r, hR, TiR, SiR, dRdT, dRdS, stable_r, PoL, PoR, KoL, KoR, hEff) + (/ 1, 1, 1, 1, 1, 1, 2, 2, 3, 3, 3, 3 /), & ! KoL + (/ 1, 1, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3 /), & ! KoR + (/ 0.00, 0.00, 0.00, 0.00, 0.00, 1.00, 0.00, 1.00, 0.00, 0.50, 1.00, 1.00 /), & ! PoL + (/ 0.00, 1.00, 0.00, 1.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.50, 1.00 /), & ! PoR + (/ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 5.00, 0.00 /), & ! hEff + 'Deep Mixed layer on the right') + + TiL(1,:) = (/ 14.00, 14.00 /); TiL(2,:) = (/ 14.00, 12.00 /); TiL(3,:) = (/ 10.00, 8.00 /); + TiR(1,:) = (/ 14.00, 14.00 /); TiR(2,:) = (/ 14.00, 14.00 /); TiR(3,:) = (/ 14.00, 14.00 /); + call mark_unstable_cells( CS, nk, Til, Sil, Pres_l, stable_l ) + call mark_unstable_cells( CS, nk, Tir, Sir, Pres_r, stable_r ) + call find_neutral_surface_positions_discontinuous(CS, nk, Pres_l, hL, TiL, SiL, ppoly_T_l, ppoly_S_l, stable_l, & + Pres_r, hR, TiR, SiR, ppoly_T_r, ppoly_S_r, stable_r, PoL, PoR, KoL, KoR, hEff) ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, KoL, KoR, PoL, PoR, hEff, & - (/1,1,1,1,2,2,2,2,3,3,3,3/), & ! KoL - (/1,1,1,1,2,2,2,2,3,3,3,3/), & ! KoR - (/0.,0.,1.,1.,0.,0.,1.,1.,0.,0.,1.,1./), & ! pL - (/0.,0.,1.,1.,0.,0.,1.,1.,0.,0.,1.,1./), & ! pR - (/0.,10.,0.,0.,0.,10.,0.,0.,0.,10.,0.,0./), & ! hEff - 'Identical columns with mixed layer') - Tl = (/20.,16.,12./) ; Tr = (/14.,14.,10./) - call build_reconstructions_1d( remap_CS, nk, hL, Tl, poly_T_l, TiL, poly_slope, iMethod, h_neglect, h_neglect_edge ) - call build_reconstructions_1d( remap_CS, nk, hR, Tr, poly_T_r, TiR, poly_slope, iMethod, h_neglect, h_neglect_edge ) - call mark_unstable_cells( nk, dRdT, dRdS, Til, Sil, stable_l, ns_l ) - call mark_unstable_cells( nk, dRdT, dRdS, Tir, Sir, stable_r, ns_r ) - call find_neutral_surface_positions_discontinuous(CS, nk, ns_l+ns_r, Pres_l, hL, TiL, SiL, dRdT, dRdS, stable_l, & - Pres_r, hR, TiR, SiR, dRdT, dRdS, stable_r, PoL, PoR, KoL, KoR, hEff) + (/ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3 /), & ! KoL + (/ 1, 1, 1, 1, 2, 2, 3, 3, 3, 3, 3, 3 /), & ! KoR + (/ 0.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00 /), & ! PoL + (/ 0.00, 0.00, 0.00, 1.00, 0.00, 1.00, 0.00, 1.00, 1.00, 1.00, 1.00, 1.00 /), & ! PoR + (/ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 /), & ! hEff + 'Right unstratified column') + + TiL(1,:) = (/ 14.00, 14.00 /); TiL(2,:) = (/ 14.00, 12.00 /); TiL(3,:) = (/ 10.00, 8.00 /); + TiR(1,:) = (/ 14.00, 14.00 /); TiR(2,:) = (/ 14.00, 14.00 /); TiR(3,:) = (/ 12.00, 4.00 /); + call mark_unstable_cells( CS, nk, Til, Sil, Pres_l, stable_l ) + call mark_unstable_cells( CS, nk, Tir, Sir, Pres_r, stable_r ) + call find_neutral_surface_positions_discontinuous(CS, nk, Pres_l, hL, TiL, SiL, ppoly_T_l, ppoly_S_l, stable_l, & + Pres_r, hR, TiR, SiR, ppoly_T_r, ppoly_S_r, stable_r, PoL, PoR, KoL, KoR, hEff) ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, KoL, KoR, PoL, PoR, hEff, & - (/1,1,2,2,2,3,3,3,3,3,3,3/), & ! KoL - (/1,1,1,1,1,1,2,2,2,3,3,3/), & ! KoR - (/0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 1.0/), & ! pL - (/0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.5, 1.0/), & ! pR - (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 5.0, 0.0/), & ! hEff - 'Right column with mixed layer') - Tl = (/14.,14.,6./) ; Tr = (/12.,16.,8./) - call build_reconstructions_1d( remap_CS, nk, hL, Tl, poly_T_l, TiL, poly_slope, iMethod, h_neglect, h_neglect_edge ) - call build_reconstructions_1d( remap_CS, nk, hR, Tr, poly_T_r, TiR, poly_slope, iMethod, h_neglect, h_neglect_edge ) - call mark_unstable_cells( nk, dRdT, dRdS, Til, Sil, stable_l, ns_l ) - call mark_unstable_cells( nk, dRdT, dRdS, Tir, Sir, stable_r, ns_r ) - call find_neutral_surface_positions_discontinuous(CS, nk, ns_l+ns_r, Pres_l, hL, TiL, SiL, dRdT, dRdS, stable_l, & - Pres_r, hR, TiR, SiR, dRdT, dRdS, stable_r, PoL, PoR, KoL, KoR, hEff) - ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 10, KoL, KoR, PoL, PoR, hEff, & - (/1,1,1,1,2,2,2,3,3,3/), & ! KoL - (/2,2,2,3,3,3,3,3,3,3/), & ! KoR - (/0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, .75, 1.0/), & ! pL - (/0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, .25, 1.0, 1.0/), & ! pR - (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 7.5, 0.0/), & ! hEff - 'Left mixed layer, right unstable mixed layer') - - Tl = (/10.,11.,6./) ; Tr = (/12.,13.,8./) - Til(:,1) = (/8.,12.,10./) ; Til(:,2) = (/12.,10.,2./) - Tir(:,1) = (/10.,14.,12./) ; TiR(:,2) = (/14.,12.,4./) - call mark_unstable_cells( nk, dRdT, dRdS, Til, Sil, stable_l, ns_l ) - call mark_unstable_cells( nk, dRdT, dRdS, Tir, Sir, stable_r, ns_r ) - call find_neutral_surface_positions_discontinuous(CS, nk, ns_l+ns_r, Pres_l, hL, TiL, SiL, dRdT, dRdS, stable_l, & - Pres_r, hR, TiR, SiR, dRdT, dRdS, stable_r, PoL, PoR, KoL, KoR, hEff) - ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 8, KoL, KoR, PoL, PoR, hEff, & - (/2,2,2,2,2,3,3,3/), & ! KoL - (/2,2,2,3,3,3,3,3/), & ! KoR - (/0.0, 0.0, 0.0, 0.0, 1.0, 0.0, .75, 1.0/), & ! pL - (/0.0, 1.0, 1.0, 0.0, .25, .25, 1.0, 1.0/), & ! pR - (/0.0, 0.0, 0.0, 4.0, 0.0, 7.5, 0.0/), & ! hEff - 'Two unstable mixed layers') - deallocate(remap_CS) - - allocate(EOS) - call EOS_manual_init(EOS, form_of_EOS = EOS_LINEAR, dRho_dT = -1., dRho_dS = 2.) - ! Unit tests for refine_nondim_position - allocate(CS%ndiff_aux_CS) - call set_ndiff_aux_params(CS%ndiff_aux_CS, deg = 1, max_iter = 10, drho_tol = 0., xtol = 0., EOS = EOS) - ! Tests using Newton's method - ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. (test_rnp(0.5,refine_nondim_position( & - CS%ndiff_aux_CS, 20., 35., -1., 2., 0., 1., (/21., -2./), (/35., 0./), -1., 1., 0.), & - "Temperature stratified (Newton) ")) - ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. (test_rnp(0.5,refine_nondim_position( & - CS%ndiff_aux_CS, 20., 35., -1., 2., 0., 1., (/20., 0./), (/34., 2./), -2., 2., 0.), & - "Salinity stratified (Newton) ")) - ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. (test_rnp(0.5,refine_nondim_position( & - CS%ndiff_aux_CS, 20., 35., -1., 2., 0., 1., (/21., -2./), (/34., 2./), -1., 1., 0.), & - "Temp/Salt stratified (Newton) ")) - call set_ndiff_aux_params(CS%ndiff_aux_CS, force_brent = .true.) - ! Tests using Brent's method - ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. (test_rnp(0.5,refine_nondim_position( & - CS%ndiff_aux_CS, 20., 35., -1., 2., 0., 1., (/21., -2./), (/35., 0./), -1., 1., 0.), & - "Temperature stratified (Brent) ")) - ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. (test_rnp(0.5,refine_nondim_position( & - CS%ndiff_aux_CS, 20., 35., -1., 2., 0., 1., (/20., 0./), (/34., 2./), -2., 2., 0.), & - "Salinity stratified (Brent) ")) - ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. (test_rnp(0.5,refine_nondim_position( & - CS%ndiff_aux_CS, 20., 35., -1., 2., 0., 1., (/21., -2./), (/34., 2./), -1., 1., 0.), & - "Temp/Salt stratified (Brent) ")) - deallocate(EOS) - + (/ 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3 /), & ! KoL + (/ 1, 1, 1, 1, 2, 2, 3, 3, 3, 3, 3, 3 /), & ! KoR + (/ 0.00, 1.00, 1.00, 1.00, 1.00, 1.00, 0.00, 1.00, 1.00, 0.00, 1.00, 1.00 /), & ! PoL + (/ 0.00, 0.00, 0.00, 1.00, 0.00, 1.00, 0.00, 0.00, 0.00, 0.25, 0.50, 1.00 /), & ! PoR + (/ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 4.00, 0.00 /), & ! hEff + 'Right unstratified column') + + TiL(1,:) = (/ 14.00, 14.00 /); TiL(2,:) = (/ 14.00, 10.00 /); TiL(3,:) = (/ 10.00, 2.00 /); + TiR(1,:) = (/ 14.00, 14.00 /); TiR(2,:) = (/ 14.00, 10.00 /); TiR(3,:) = (/ 10.00, 2.00 /); + call mark_unstable_cells( CS, nk, Til, Sil, Pres_l, stable_l ) + call mark_unstable_cells( CS, nk, Tir, Sir, Pres_r, stable_r ) + call find_neutral_surface_positions_discontinuous(CS, nk, Pres_l, hL, TiL, SiL, ppoly_T_l, ppoly_S_l, stable_l, & + Pres_r, hR, TiR, SiR, ppoly_T_r, ppoly_S_r, stable_r, PoL, PoR, KoL, KoR, hEff) + ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, KoL, KoR, PoL, PoR, hEff, & + (/ 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3 /), & ! KoL + (/ 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 3 /), & ! KoR + (/ 0.00, 1.00, 1.00, 1.00, 0.00, 0.00, 1.00, 1.00, 0.00, 0.00, 1.00, 1.00 /), & ! PoL + (/ 0.00, 0.00, 0.00, 1.00, 0.00, 0.00, 1.00, 0.00, 0.00, 0.00, 1.00, 1.00 /), & ! PoR + (/ 0.00, 0.00, 0.00, 0.00, 0.00, 10.00, 0.00, 0.00, 0.00, 10.00, 0.00 /), & ! hEff + 'Identical columns with mixed layer') + + TiL(1,:) = (/ 14.00, 12.00 /); TiL(2,:) = (/ 10.00, 10.00 /); TiL(3,:) = (/ 8.00, 2.00 /); + TiR(1,:) = (/ 14.00, 12.00 /); TiR(2,:) = (/ 12.00, 8.00 /); TiR(3,:) = (/ 8.00, 2.00 /); + call mark_unstable_cells( CS, nk, Til, Sil, Pres_l, stable_l ) + call mark_unstable_cells( CS, nk, Tir, Sir, Pres_r, stable_r ) + call find_neutral_surface_positions_discontinuous(CS, nk, Pres_l, hL, TiL, SiL, ppoly_T_l, ppoly_S_l, stable_l, & + Pres_r, hR, TiR, SiR, ppoly_T_r, ppoly_S_r, stable_r, PoL, PoR, KoL, KoR, hEff) + ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, KoL, KoR, PoL, PoR, hEff, & + (/ 1, 1, 1, 1, 2, 2, 3, 3, 3, 3, 3, 3 /), & ! KoL + (/ 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3 /), & ! KoR + (/ 0.00, 0.00, 1.00, 1.00, 0.00, 1.00, 0.00, 0.00, 0.00, 0.00, 1.00, 1.00 /), & ! PoL + (/ 0.00, 0.00, 1.00, 0.00, 0.00, 0.00, 0.00, 1.00, 1.00, 0.00, 1.00, 1.00 /), & ! PoR + (/ 0.00, 10.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 10.00, 0.00 /), & ! hEff + 'Left interior unstratified') + + TiL(1,:) = (/ 12.00, 12.00 /); TiL(2,:) = (/ 12.00, 10.00 /); TiL(3,:) = (/ 10.00, 6.00 /); + TiR(1,:) = (/ 12.00, 10.00 /); TiR(2,:) = (/ 10.00, 12.00 /); TiR(3,:) = (/ 8.00, 4.00 /); + call mark_unstable_cells( CS, nk, Til, Sil, Pres_l, stable_l ) + call mark_unstable_cells( CS, nk, Tir, Sir, Pres_r, stable_r ) + call find_neutral_surface_positions_discontinuous(CS, nk, Pres_l, hL, TiL, SiL, ppoly_T_l, ppoly_S_l, stable_l, & + Pres_r, hR, TiR, SiR, ppoly_T_r, ppoly_S_r, stable_r, PoL, PoR, KoL, KoR, hEff) + ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, KoL, KoR, PoL, PoR, hEff, & + (/ 1, 1, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3 /), & ! KoL + (/ 1, 1, 1, 1, 1, 1, 2, 2, 3, 3, 3, 3 /), & ! KoR + (/ 0.00, 1.00, 0.00, 0.00, 1.00, 0.00, 0.00, 0.00, 0.00, 0.50, 1.00, 1.00 /), & ! PoL + (/ 0.00, 0.00, 0.00, 0.00, 1.00, 1.00, 0.00, 1.00, 0.00, 0.00, 0.50, 1.00 /), & ! PoR + (/ 0.00, 0.00, 0.00, 10.00, 0.00, 0.00, 0.00, 0.00, 0.00, 5.00, 0.00 /), & ! hEff + 'Left mixed layer, Right unstable interior') + + TiL(1,:) = (/ 14.00, 14.00 /); TiL(2,:) = (/ 10.00, 10.00 /); TiL(3,:) = (/ 8.00, 6.00 /); + TiR(1,:) = (/ 10.00, 14.00 /); TiR(2,:) = (/ 16.00, 16.00 /); TiR(3,:) = (/ 12.00, 4.00 /); + call mark_unstable_cells( CS, nk, Til, Sil, Pres_l, stable_l ) + call mark_unstable_cells( CS, nk, Tir, Sir, Pres_r, stable_r ) + call find_neutral_surface_positions_discontinuous(CS, nk, Pres_l, hL, TiL, SiL, ppoly_T_l, ppoly_S_l, stable_l, & + Pres_r, hR, TiR, SiR, ppoly_T_r, ppoly_S_r, stable_r, PoL, PoR, KoL, KoR, hEff) + ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, KoL, KoR, PoL, PoR, hEff, & + (/ 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3 /), & ! KoL + (/ 1, 1, 1, 1, 1, 1, 2, 2, 3, 3, 3, 3 /), & ! KoR + (/ 0.00, 1.00, 0.00, 1.00, 1.00, 1.00, 1.00, 1.00, 0.00, 0.00, 1.00, 1.00 /), & ! PoL + (/ 0.00, 0.00, 0.00, 0.00, 0.00, 1.00, 0.00, 1.00, 0.00, 0.50, 0.75, 1.00 /), & ! PoR + (/ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 4.00, 0.00 /), & ! hEff + 'Left thick mixed layer, Right unstable mixed') + + TiL(1,:) = (/ 8.00, 12.00 /); TiL(2,:) = (/ 12.00, 10.00 /); TiL(3,:) = (/ 8.00, 4.00 /); + TiR(1,:) = (/ 10.00, 14.00 /); TiR(2,:) = (/ 14.00, 12.00 /); TiR(3,:) = (/ 10.00, 6.00 /); + call mark_unstable_cells( CS, nk, Til, Sil, Pres_l, stable_l ) + call mark_unstable_cells( CS, nk, Tir, Sir, Pres_r, stable_r ) + call find_neutral_surface_positions_discontinuous(CS, nk, Pres_l, hL, TiL, SiL, ppoly_T_l, ppoly_S_l, stable_l, & + Pres_r, hR, TiR, SiR, ppoly_T_r, ppoly_S_r, stable_r, PoL, PoR, KoL, KoR, hEff) + ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, KoL, KoR, PoL, PoR, hEff, & + (/ 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3 /), & ! KoL + (/ 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 3 /), & ! KoR + (/ 0.00, 1.00, 1.00, 1.00, 0.00, 0.00, 0.00, 1.00, 0.00, 0.00, 0.50, 1.00 /), & ! PoL + (/ 0.00, 0.00, 0.00, 1.00, 0.00, 1.00, 1.00, 0.00, 0.00, 0.50, 1.00, 1.00 /), & ! PoR + (/ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 5.00, 0.00 /), & ! hEff + 'Unstable mixed layers, left cooler') + + call EOS_manual_init(CS%EOS, form_of_EOS = EOS_LINEAR, dRho_dT = -1., dRho_dS = 2.) + ! Tests for linearized version of searching the layer for neutral surface position + ! EOS linear in T, uniform alpha + CS%max_iter = 10 + ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. (test_rnp(0.5, & + find_neutral_pos_linear(CS, 0., 10., 35., 0., -0.2, 0., & + 0., -0.2, 0., 10., -0.2, 0., & + (/12.,-4./), (/34.,0./)), "Temp Uniform Linearized Alpha/Beta")) + ! EOS linear in S, uniform beta + ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. (test_rnp(0.5, & + find_neutral_pos_linear(CS, 0., 10., 35., 0., 0., 0.8, & + 0., 0., 0.8, 10., 0., 0.8, & + (/12.,0./), (/34.,2./)), "Salt Uniform Linearized Alpha/Beta")) + ! EOS linear in T/S, uniform alpha/beta + ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. (test_rnp(0.5, & + find_neutral_pos_linear(CS, 0., 10., 35., 0., -0.5, 0.5, & + 0., -0.5, 0.5, 10., -0.5, 0.5, & + (/12.,-4./), (/34.,2./)), "Temp/salt Uniform Linearized Alpha/Beta")) + ! EOS linear in T, insensitive to So + ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. (test_rnp(0.5, & + find_neutral_pos_linear(CS, 0., 10., 35., 0., -0.2, 0., & + 0., -0.4, 0., 10., -0.6, 0., & + (/12.,-4./), (/34.,0./)), "Temp stratified Linearized Alpha/Beta")) +! ! EOS linear in S, insensitive to T + ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. (test_rnp(0.5, & + find_neutral_pos_linear(CS, 0., 10., 35., 0., 0., 0.8, & + 0., 0., 1.0, 10., 0., 0.5, & + (/12.,0./), (/34.,2./)), "Salt stratified Linearized Alpha/Beta")) if (.not. ndiff_unit_tests_discontinuous) write(*,*) 'Pass' end function ndiff_unit_tests_discontinuous @@ -2230,7 +2807,7 @@ logical function test_rnp(expected_pos, test_pos, title) character(len=*), intent(in) :: title !< A label for this test ! Local variables integer :: stdunit = 6 ! Output to standard error - test_rnp = expected_pos /= test_pos + test_rnp = ABS(expected_pos - test_pos) > 2*EPSILON(test_pos) if (test_rnp) then write(stdunit,'(A, f20.16, " .neq. ", f20.16, " <-- WRONG")') title, expected_pos, test_pos else diff --git a/src/tracer/MOM_neutral_diffusion_aux.F90 b/src/tracer/MOM_neutral_diffusion_aux.F90 deleted file mode 100644 index 88df1ddbc5..0000000000 --- a/src/tracer/MOM_neutral_diffusion_aux.F90 +++ /dev/null @@ -1,669 +0,0 @@ -!> A column-wise toolbox for implementing neutral diffusion -module MOM_neutral_diffusion_aux - -use MOM_EOS, only : EOS_type, extract_member_EOS, EOS_LINEAR, EOS_TEOS10, EOS_WRIGHT -use MOM_EOS, only : calculate_density_derivs, calculate_density_second_derivs -use MOM_error_handler, only : MOM_error, FATAL, WARNING -use polynomial_functions, only : evaluation_polynomial, first_derivative_polynomial - -! This file is part of MOM6. See LICENSE.md for the license. -implicit none ; private - -public set_ndiff_aux_params -public mark_unstable_cells -public increment_interface -public calc_drho -public drho_at_pos -public search_other_column -public interpolate_for_nondim_position -public refine_nondim_position -public check_neutral_positions -public kahan_sum - -!> The control structure for this module -type, public :: ndiff_aux_CS_type ; private - integer :: nterm !< Number of terms in polynomial (deg+1) - integer :: max_iter !< Maximum number of iterations - real :: drho_tol !< Tolerance criterion for difference in density [kg m-3] - real :: xtol !< Criterion for how much position changes [nondim] - real :: ref_pres !< Determines whether a constant reference pressure is used everywhere or locally referenced - !< density is done. ref_pres <-1 is the latter, ref_pres >= 0. otherwise - logical :: force_brent = .false. !< Use Brent's method instead of Newton even when second derivatives are available - logical :: debug !< If true, write verbose debugging messages and checksusm - type(EOS_type), pointer :: EOS !< Pointer to equation of state used in the model -end type ndiff_aux_CS_type - -contains - -!> Initialize the parameters used to iteratively find the neutral direction -subroutine set_ndiff_aux_params( CS, deg, max_iter, drho_tol, xtol, ref_pres, force_brent, EOS, debug) - type(ndiff_aux_CS_type), intent(inout) :: CS !< Control structure for refine_pos - integer, optional, intent(in ) :: deg !< Degree of polynommial used in reconstruction - integer, optional, intent(in ) :: max_iter !< Maximum number of iterations - real, optional, intent(in ) :: drho_tol !< Tolerance for function convergence - real, optional, intent(in ) :: xtol !< Tolerance for change in position - real, optional, intent(in ) :: ref_pres !< Reference pressure to use - logical, optional, intent(in ) :: force_brent !< Force Brent method for linear, TEOS-10, and WRIGHT - logical, optional, intent(in ) :: debug !< If true, print output use to help debug neutral diffusion - type(EOS_type), target, optional, intent(in ) :: EOS !< Equation of state - - if (present( deg )) CS%nterm = deg + 1 - if (present( max_iter )) CS%max_iter = max_iter - if (present( drho_tol )) CS%drho_tol = drho_tol - if (present( xtol )) CS%xtol = xtol - if (present( ref_pres )) CS%ref_pres = ref_pres - if (present( force_brent )) CS%force_brent = force_brent - if (present( EOS )) CS%EOS => EOS - if (present( debug )) CS%debug = debug - -end subroutine set_ndiff_aux_params - -!> Given the reconsturcitons of dRdT, dRdS, T, S mark the cells which are stably stratified parts of the water column -!! For an layer to be unstable the top interface must be denser than the bottom or the bottom interface of the layer -subroutine mark_unstable_cells(nk, dRdT, dRdS,T, S, stable_cell, ns) - integer, intent(in) :: nk !< Number of levels in a column - real, dimension(nk,2), intent(in) :: dRdT !< drho/dT [kg m-3 degC-1] at interfaces - real, dimension(nk,2), intent(in) :: dRdS !< drho/dS [kg m-3 ppt-1] at interfaces - real, dimension(nk,2), intent(in) :: T !< Temperature [degC] at interfaces - real, dimension(nk,2), intent(in) :: S !< Salinity [ppt] at interfaces - logical, dimension(nk), intent( out) :: stable_cell !< True if this cell is unstably stratified - integer, intent( out) :: ns !< Number of neutral surfaces in unmasked part of the column - - integer :: k, first_stable, prev_stable - real :: delta_rho - - ! First check to make sure that density profile between the two interfaces of the cell are stable - ! Note that we neglect a factor of 0.5 because we only care about the sign of delta_rho not magnitude - do k = 1,nk - ! Compare density of bottom interface to top interface, should be positive (or zero) if stable - delta_rho = (dRdT(k,2) + dRdT(k,1))*(T(k,2) - T(k,1)) + (dRdS(k,2) + dRdS(k,1))*(S(k,2) - S(k,1)) - stable_cell(k) = delta_rho >= 0. - enddo - - first_stable = 1 - ! Check to see that bottom interface of upper cell is lighter than the upper interface of the lower cell - do k=1,nk - if (stable_cell(k)) then - first_stable = k - exit - endif - enddo - prev_stable = first_stable - - ! Start either with the first stable cell or the layer just below the surface - do k = prev_stable+1, nk - ! Don't do anything if the cell has already been marked as unstable - if (.not. stable_cell(k)) cycle - ! Otherwise, we need to check to see if this cell's upper interface is denser than the previous stable_cell - ! Compare top interface of lower cell to bottom interface of upper cell, positive or zero if bottom cell is stable - delta_rho = (dRdT(k,1) + dRdT(prev_stable,2))*(T(k,1) - T(prev_stable,2)) + & - (dRdS(k,1) + dRdS(prev_stable,2))*(S(k,1) - S(prev_stable,2)) - stable_cell(k) = delta_rho >= 0. - ! If the lower cell is marked as stable, then it should be the next reference cell - if (stable_cell(k)) prev_stable = k - enddo - - ! Number of interfaces is the 2 times number of stable cells in the water column - ns = 0 - do k = 1,nk - if (stable_cell(k)) ns = ns + 2 - enddo - -end subroutine mark_unstable_cells - -!> Increments the interface which was just connected and also set flags if the bottom is reached -subroutine increment_interface(nk, kl, ki, stable, reached_bottom, searching_this_column, searching_other_column) - integer, intent(in ) :: nk !< Number of vertical levels - integer, intent(inout) :: kl !< Current layer (potentially updated) - integer, intent(inout) :: ki !< Current interface - logical, dimension(nk), intent(in ) :: stable !< True if the cell is stably stratified - logical, intent(inout) :: reached_bottom !< Updated when kl == nk and ki == 2 - logical, intent(inout) :: searching_this_column !< Updated when kl == nk and ki == 2 - logical, intent(inout) :: searching_other_column !< Updated when kl == nk and ki == 2 - integer :: k - - if (ki == 1) then - ki = 2 - elseif ((ki == 2) .and. (kl < nk) ) then - do k = kl+1,nk - if (stable(kl)) then - kl = k - ki = 1 - exit - endif - ! If we did not find another stable cell, then the current cell is essentially the bottom - ki = 2 - reached_bottom = .true. - searching_this_column = .true. - searching_other_column = .false. - enddo - elseif ((kl == nk) .and. (ki==2)) then - reached_bottom = .true. - searching_this_column = .true. - searching_other_column = .false. - else - call MOM_error(FATAL,"Unanticipated eventuality in increment_interface") - endif -end subroutine increment_interface - -!> Calculates difference in density at two points (rho1-rho2) with known density derivatives, T, and S -real function calc_drho(T1, S1, dRdT1, dRdS1, T2, S2, dRdT2, dRdS2) - real, intent(in ) :: T1 !< Temperature at point 1 - real, intent(in ) :: S1 !< Salinity at point 1 - real, intent(in ) :: dRdT1 !< dRhodT at point 1 - real, intent(in ) :: dRdS1 !< dRhodS at point 1 - real, intent(in ) :: T2 !< Temperature at point 2 - real, intent(in ) :: S2 !< Salinity at point 2 - real, intent(in ) :: dRdT2 !< dRhodT at point 2 - real, intent(in ) :: dRdS2 !< dRhodS at point - - calc_drho = 0.5*( (dRdT1+dRdT2)*(T1-T2) + (dRdS1+dRdS2)*(S1-S2) ) -end function calc_drho - -!> Calculate the difference in neutral density between a reference T, S, alpha, and beta -!! at a point on the polynomial reconstructions of T, S -subroutine drho_at_pos(CS, T_ref, S_ref, alpha_ref, beta_ref, P_top, P_bot, ppoly_T, ppoly_S, x0, & - delta_rho, P_out, T_out, S_out, alpha_avg_out, beta_avg_out, delta_T_out, delta_S_out) - type(ndiff_aux_CS_type), intent(in) :: CS !< Control structure with parameters for this module - real, intent(in) :: T_ref !< Temperature at reference surface - real, intent(in) :: S_ref !< Salinity at reference surface - real, intent(in) :: alpha_ref !< dRho/dT at reference surface - real, intent(in) :: beta_ref !< dRho/dS at reference surface - real, intent(in) :: P_top !< Pressure (Pa) at top interface of layer to be searched - real, intent(in) :: P_bot !< Pressure (Pa) at bottom interface - real, dimension(CS%nterm), intent(in) :: ppoly_T !< Coefficients of T reconstruction - real, dimension(CS%nterm), intent(in) :: ppoly_S !< Coefficients of S reconstruciton - real, intent(in) :: x0 !< Nondimensional position to evaluate - real, intent(out) :: delta_rho !< The density difference from a reference value - real, optional, intent(out) :: P_out !< Pressure at point x0 - real, optional, intent(out) :: T_out !< Temperature at point x0 - real, optional, intent(out) :: S_out !< Salinity at point x0 - real, optional, intent(out) :: alpha_avg_out !< Average of alpha between reference and x0 - real, optional, intent(out) :: beta_avg_out !< Average of beta between reference and x0 - real, optional, intent(out) :: delta_T_out !< Difference in temperature between reference and x0 - real, optional, intent(out) :: delta_S_out !< Difference in salinity between reference and x0 - - real :: alpha, beta, alpha_avg, beta_avg, P_int, T, S, delta_T, delta_S - - P_int = (1. - x0)*P_top + x0*P_bot - T = evaluation_polynomial( ppoly_T, CS%nterm, x0 ) - S = evaluation_polynomial( ppoly_S, CS%nterm, x0 ) - ! Interpolated pressure if using locally referenced neutral density - if (CS%ref_pres<0.) then - call calculate_density_derivs( T, S, P_int, alpha, beta, CS%EOS ) - else - ! Constant reference pressure (isopycnal) - call calculate_density_derivs( T, S, CS%ref_pres, alpha, beta, CS%EOS ) - endif - - ! Calculate the f(P) term for Newton's method - alpha_avg = 0.5*( alpha + alpha_ref ) - beta_avg = 0.5*( beta + beta_ref ) - delta_T = T - T_ref - delta_S = S - S_ref - delta_rho = alpha_avg*delta_T + beta_avg*delta_S - - ! If doing a Newton step, these quantities are needed, otherwise they can just be optional - if (present(P_out)) P_out = P_int - if (present(T_out)) T_out = T - if (present(S_out)) S_out = S - if (present(alpha_avg_out)) alpha_avg_out = alpha_avg - if (present(beta_avg_out)) beta_avg_out = beta_avg - if (present(delta_T_out)) delta_T_out = delta_T - if (present(delta_S_out)) delta_S_out = delta_S - -end subroutine drho_at_pos - -!> Searches the "other" (searched) column for the position of the neutral surface -subroutine search_other_column(dRhoTop, dRhoBot, Ptop, Pbot, lastP, lastK, kl, kl_0, ki, & - top_connected, bot_connected, out_P, out_K, search_layer) - real, intent(in ) :: dRhoTop !< Density difference across top interface - real, intent(in ) :: dRhoBot !< Density difference across top interface - real, intent(in ) :: Ptop !< Pressure at top interface - real, intent(in ) :: Pbot !< Pressure at bottom interface - real, intent(in ) :: lastP !< Last position connected in the searched column - integer, intent(in ) :: lastK !< Last layer connected in the searched column - integer, intent(in ) :: kl !< Layer in the searched column - integer, intent(in ) :: kl_0 !< Layer in the searched column - integer, intent(in ) :: ki !< Interface of the searched column - logical, dimension(:), intent(inout) :: top_connected !< True if the top interface was pointed to - logical, dimension(:), intent(inout) :: bot_connected !< True if the top interface was pointed to - real, intent( out) :: out_P !< Position within searched column - integer, intent( out) :: out_K !< Layer within searched column - logical, intent( out) :: search_layer !< Neutral surface within cell - - search_layer = .false. - if (kl > kl_0) then ! Away from top cell - if (kl == lastK) then ! Searching in the same layer - if (dRhoTop > 0.) then - if (lastK == kl) then - out_P = lastP - else - out_P = 0. - endif - out_K = kl -! out_P = max(0.,lastP) ; out_K = kl - elseif ( dRhoTop == dRhoBot ) then - if (top_connected(kl)) then - out_P = 1. ; out_K = kl - else - out_P = max(0.,lastP) ; out_K = kl - endif - elseif (dRhoTop >= dRhoBot) then - out_P = 1. ; out_K = kl - elseif (dRhoTop < 0. .and. dRhoBot < 0.)then - out_P = 1. ; out_K = kl - else - out_K = kl - out_P = max(interpolate_for_nondim_position( dRhoTop, Ptop, dRhoBot, Pbot ),lastP) - search_layer = .true. - endif - else ! Searching across the interface - if (.not. bot_connected(kl-1) ) then - out_K = kl-1 - out_P = 1. - else - out_K = kl - out_P = 0. - endif - endif - else ! At the top cell - if (ki == 1) then - out_P = 0. ; out_K = kl - elseif (dRhoTop > 0.) then - if (lastK == kl) then - out_P = lastP - else - out_P = 0. - endif - out_K = kl -! out_P = max(0.,lastP) ; out_K = kl - elseif ( dRhoTop == dRhoBot ) then - if (top_connected(kl)) then - out_P = 1. ; out_K = kl - else - out_P = max(0.,lastP) ; out_K = kl - endif - elseif (dRhoTop >= dRhoBot) then - out_P = 1. ; out_K = kl - elseif (dRhoTop < 0. .and. dRhoBot < 0.)then - out_P = 1. ; out_K = kl - else - out_K = kl - out_P = max(interpolate_for_nondim_position( dRhoTop, Ptop, dRhoBot, Pbot ),lastP) - search_layer = .true. - endif - endif - -end subroutine search_other_column - -!> Returns the non-dimensional position between Pneg and Ppos where the -!! interpolated density difference equals zero. -!! The result is always bounded to be between 0 and 1. -real function interpolate_for_nondim_position(dRhoNeg, Pneg, dRhoPos, Ppos) - real, intent(in) :: dRhoNeg !< Negative density difference - real, intent(in) :: Pneg !< Position of negative density difference - real, intent(in) :: dRhoPos !< Positive density difference - real, intent(in) :: Ppos !< Position of positive density difference - - if (PposdRhoPos) then - write(0,*) 'dRhoNeg, Pneg, dRhoPos, Ppos=',dRhoNeg, Pneg, dRhoPos, Ppos - elseif (dRhoNeg>dRhoPos) then - stop 'interpolate_for_nondim_position: Houston, we have a problem! dRhoNeg>dRhoPos' - endif - if (Ppos<=Pneg) then ! Handle vanished or inverted layers - interpolate_for_nondim_position = 0.5 - elseif ( dRhoPos - dRhoNeg > 0. ) then - interpolate_for_nondim_position = min( 1., max( 0., -dRhoNeg / ( dRhoPos - dRhoNeg ) ) ) - elseif ( dRhoPos - dRhoNeg == 0) then - if (dRhoNeg>0.) then - interpolate_for_nondim_position = 0. - elseif (dRhoNeg<0.) then - interpolate_for_nondim_position = 1. - else ! dRhoPos = dRhoNeg = 0 - interpolate_for_nondim_position = 0.5 - endif - else ! dRhoPos - dRhoNeg < 0 - interpolate_for_nondim_position = 0.5 - endif - if ( interpolate_for_nondim_position < 0. ) & - stop 'interpolate_for_nondim_position: Houston, we have a problem! Pint < Pneg' - if ( interpolate_for_nondim_position > 1. ) & - stop 'interpolate_for_nondim_position: Houston, we have a problem! Pint > Ppos' -end function interpolate_for_nondim_position - -!> Use root-finding methods to find where dRho = 0, based on the equation of state and the polynomial -!! reconstructions of temperature, salinity. Initial guess is based on the zero crossing of based on linear -!! profiles of dRho, T, and S, between the top and bottom interface. If second derivatives of the EOS are available, -!! it starts with a Newton's method. However, Newton's method is not guaranteed to be bracketed, a check is performed -!! to see if it it diverges outside the interval. In that case (or in the case that second derivatives are not -!! available), Brent's method is used following the implementation found at -!! https://people.sc.fsu.edu/~jburkardt/f_src/brent/brent.f90 -real function refine_nondim_position(CS, T_ref, S_ref, alpha_ref, beta_ref, P_top, P_bot, & - ppoly_T, ppoly_S, drho_top, drho_bot, min_bound) - type(ndiff_aux_CS_type), intent(in) :: CS !< Control structure with parameters for this module - real, intent(in) :: T_ref !< Temperature of the neutral surface at the searched from interface - real, intent(in) :: S_ref !< Salinity of the neutral surface at the searched from interface - real, intent(in) :: alpha_ref !< dRho/dT of the neutral surface at the searched from interface - real, intent(in) :: beta_ref !< dRho/dS of the neutral surface at the searched from interface - real, intent(in) :: P_top !< Pressure at the top interface in the layer to be searched - real, intent(in) :: P_bot !< Pressure at the bottom interface in the layer to be searched - real, dimension(:), intent(in) :: ppoly_T !< Coefficients of the order N polynomial reconstruction of T within - !! the layer to be searched. - real, dimension(:), intent(in) :: ppoly_S !< Coefficients of the order N polynomial reconstruction of T within - !! the layer to be searched. - real, intent(in) :: drho_top !< Delta rho at top interface (or previous position in cell - real, intent(in) :: drho_bot !< Delta rho at bottom interface - real, intent(in) :: min_bound !< Lower bound of position, also serves as the initial guess - - ! Local variables - integer :: form_of_EOS - integer :: iter - logical :: do_newton, do_brent - - real :: delta_rho, d_delta_rho_dP ! Terms for the Newton iteration - real :: P_int, P_min, P_ref ! Interpolated pressure - real :: delta_rho_init, delta_rho_final - real :: neg_x, neg_fun - real :: T, S, alpha, beta, alpha_avg, beta_avg - ! Newton's Method with variables - real :: dT_dP, dS_dP, delta_T, delta_S, delta_P - real :: dbeta_dS, dbeta_dT, dalpha_dT, dalpha_dS, dbeta_dP, dalpha_dP - real :: a, b, c, b_last - ! Extra Brent's Method variables - real :: d, e, f, fa, fb, fc, m, p, q, r, s0, sa, sb, tol, machep - - real :: P_last - - machep = EPSILON(machep) - if (CS%ref_pres>=0.) P_ref = CS%ref_pres - delta_P = P_bot-P_top - refine_nondim_position = min_bound - - call extract_member_EOS(CS%EOS, form_of_EOS = form_of_EOS) - do_newton = (form_of_EOS == EOS_LINEAR) .or. (form_of_EOS == EOS_TEOS10) .or. (form_of_EOS == EOS_WRIGHT) - do_brent = .not. do_newton - if (CS%force_brent) then - do_newton = .not. CS%force_brent - do_brent = CS%force_brent - endif - - ! Calculate the initial values - call drho_at_pos(CS, T_ref, S_ref, alpha_ref, beta_ref, P_top, P_bot, ppoly_T, ppoly_S, min_bound, & - delta_rho, P_int, T, S, alpha_avg, beta_avg, delta_T, delta_S) - delta_rho_init = delta_rho - if ( ABS(delta_rho_init) <= CS%drho_tol ) then - refine_nondim_position = min_bound - return - endif - if (ABS(drho_bot) <= CS%drho_tol) then - refine_nondim_position = 1. - return - endif - - ! Set the initial values to ensure that the algorithm returns a 'negative' value - neg_fun = delta_rho - neg_x = min_bound - - if (CS%debug) then - write (*,*) "------" - write (*,*) "Starting x0, delta_rho: ", min_bound, delta_rho - endif - - ! For now only linear, Wright, and TEOS-10 equations of state have functions providing second derivatives and - ! thus can use Newton's method for the equation of state - if (do_newton) then - refine_nondim_position = min_bound - ! Set lower bound of pressure - P_min = P_top*(1.-min_bound) + P_bot*(min_bound) - fa = delta_rho_init ; a = min_bound - fb = delta_rho_init ; b = min_bound - fc = drho_bot ; c = 1. - ! Iterate over Newton's method for the function: x0 = x0 - delta_rho/d_delta_rho_dP - do iter = 1, CS%max_iter - P_int = P_top*(1. - b) + P_bot*b - ! Evaluate total derivative of delta_rho - if (CS%ref_pres<0.) P_ref = P_int - call calculate_density_second_derivs( T, S, P_ref, dbeta_dS, dbeta_dT, dalpha_dT, dbeta_dP, dalpha_dP, CS%EOS ) - ! In the case of a constant reference pressure, no dependence on neutral direction with pressure - if (CS%ref_pres>=0.) then - dalpha_dP = 0. ; dbeta_dP = 0. - endif - dalpha_dS = dbeta_dT ! Cross derivatives are identicial - ! By chain rule dT_dP= (dT_dz)*(dz/dP) = dT_dz / (Pbot-Ptop) - dT_dP = first_derivative_polynomial( ppoly_T, CS%nterm, b ) / delta_P - dS_dP = first_derivative_polynomial( ppoly_S, CS%nterm, b ) / delta_P - ! Total derivative of d_delta_rho wrt P - d_delta_rho_dP = 0.5*( delta_S*(dS_dP*dbeta_dS + dT_dP*dbeta_dT + dbeta_dP) + & - ( delta_T*(dS_dP*dalpha_dS + dT_dP*dalpha_dT + dalpha_dP))) + & - dS_dP*beta_avg + dT_dP*alpha_avg - ! This probably won't happen, but if it does take a bisection - if (d_delta_rho_dP == 0.) then - b = 0.5*(a+c) - else - ! Newton step update - P_int = P_int - (fb / d_delta_rho_dP) - ! This line is equivalent to the next - ! refine_nondim_position = (P_top-P_int)/(P_top-P_bot) - b_last = b - b = (P_int-P_top)/delta_P - ! Test to see if it fell out of the bracketing interval. If so, take a bisection step - if (b < a .or. b > c) then - b = 0.5*(a + c) - endif - endif - call drho_at_pos(CS, T_ref, S_ref, alpha_ref, beta_ref, P_top, P_bot, ppoly_T, ppoly_S, & - b, fb, P_int, T, S, alpha_avg, beta_avg, delta_T, delta_S) - if (CS%debug) write(*,'(A,I3.3,X,ES23.15,X,ES23.15)') "Iteration, b, fb: ", iter, b, fb - - if (fb < 0. .and. fb > neg_fun) then - neg_fun = fb - neg_x = b - endif - - ! For the logic to find neutral surfaces to work properly, the function needs to converge to zero - ! or a small negative value - if ((fb <= 0.) .and. (fb >= -CS%drho_tol)) then - refine_nondim_position = b - exit - endif - ! Exit if method has stalled out - if ( ABS(b-b_last)<=CS%xtol ) then - refine_nondim_position = b - exit - endif - - ! Update the bracket - if (SIGN(1.,fa)*SIGN(1.,fb)<0.) then - c = b - fc = delta_rho - else - a = b - fa = delta_rho - endif - enddo - refine_nondim_position = b - delta_rho = fb - endif - if (delta_rho > 0.) then - refine_nondim_position = neg_x - delta_rho = neg_fun - endif - ! Do Brent if analytic second derivatives don't exist - if (do_brent) then - sa = max(refine_nondim_position,min_bound) ; fa = delta_rho - sb = 1. ; fb = drho_bot - c = sa ; fc = fa ; e = sb - sa; d = e - - - ! This is from https://people.sc.fsu.edu/~jburkardt/f_src/brent/brent.f90 - do iter = 1,CS%max_iter - if ( abs ( fc ) < abs ( fb ) ) then - sa = sb - sb = c - c = sa - fa = fb - fb = fc - fc = fa - endif - tol = 2. * machep * abs ( sb ) + CS%xtol - m = 0.5 * ( c - sb ) - if ( abs ( m ) <= tol .or. fb == 0. ) then - exit - endif - if ( abs ( e ) < tol .or. abs ( fa ) <= abs ( fb ) ) then - e = m - d = e - else - s0 = fb / fa - if ( sa == c ) then - p = 2. * m * s0 - q = 1. - s0 - else - q = fa / fc - r = fb / fc - p = s0 * ( 2. * m * q * ( q - r ) - ( sb - sa ) * ( r - 1. ) ) - q = ( q - 1. ) * ( r - 1. ) * ( s0 - 1. ) - endif - if ( 0. < p ) then - q = - q - else - p = - p - endif - s0 = e - e = d - if ( 2. * p < 3. * m * q - abs ( tol * q ) .and. & - p < abs ( 0.5 * s0 * q ) ) then - d = p / q - else - e = m - d = e - endif - endif - sa = sb - fa = fb - if ( tol < abs ( d ) ) then - sb = sb + d - elseif ( 0. < m ) then - sb = sb + tol - else - sb = sb - tol - endif - call drho_at_pos(CS, T_ref, S_ref, alpha_ref, beta_ref, P_top, P_bot, ppoly_T, ppoly_S, & - sb, fb) - if ( ( 0. < fb .and. 0. < fc ) .or. & - ( fb <= 0. .and. fc <= 0. ) ) then - c = sa - fc = fa - e = sb - sa - d = e - endif - enddo - ! Modified from original to ensure that the minimum is found - fa = ABS(fa) ; fb = ABS(fb) ; fc = ABS(fc) - delta_rho = MIN(fa, fb, fc) - - if (fb==delta_rho) then - refine_nondim_position = max(sb,min_bound) - elseif (fa==delta_rho) then - refine_nondim_position = max(sa,min_bound) - elseif (fc==delta_rho) then - refine_nondim_position = max(c, min_bound) - endif - endif - - ! Make sure that the result is bounded between 0 and 1 - if (refine_nondim_position>1.) then - if (CS%debug) then - write (*,*) "T, T Poly Coeffs: ", T, ppoly_T - write (*,*) "S, S Poly Coeffs: ", S, ppoly_S - write (*,*) "T_ref, alpha_ref: ", T_ref, alpha_ref - write (*,*) "S_ref, beta_ref : ", S_ref, beta_ref - write (*,*) "x0: ", min_bound - write (*,*) "refine_nondim_position: ", refine_nondim_position - endif - call MOM_error(WARNING, "refine_nondim_position>1.") - refine_nondim_position = 1. - endif - - if (refine_nondim_position Do a compensated sum to account for roundoff level -subroutine kahan_sum(sum, summand, c) - real, intent(inout) :: sum !< Running sum - real, intent(in ) :: summand !< Term to be added - real ,intent(inout) :: c !< Keep track of roundoff - real :: y, t - y = summand - c - t = sum + y - c = (t-sum) - y - sum = t - -end subroutine kahan_sum - -end module MOM_neutral_diffusion_aux diff --git a/src/tracer/MOM_tracer_hor_diff.F90 b/src/tracer/MOM_tracer_hor_diff.F90 index 098a647ec8..848841caf6 100644 --- a/src/tracer/MOM_tracer_hor_diff.F90 +++ b/src/tracer/MOM_tracer_hor_diff.F90 @@ -11,6 +11,7 @@ module MOM_tracer_hor_diff use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type use MOM_domains, only : pass_vector use MOM_debugging, only : hchksum, uvchksum +use MOM_diabatic_driver, only : diabatic_CS use MOM_EOS, only : calculate_density, EOS_type use MOM_error_handler, only : MOM_error, FATAL, WARNING, MOM_mesg, is_root_pe use MOM_error_handler, only : MOM_set_verbosity, callTree_showQuery @@ -22,6 +23,8 @@ module MOM_tracer_hor_diff use MOM_neutral_diffusion, only : neutral_diffusion_init, neutral_diffusion_end use MOM_neutral_diffusion, only : neutral_diffusion_CS use MOM_neutral_diffusion, only : neutral_diffusion_calc_coeffs, neutral_diffusion +use MOM_lateral_boundary_diffusion, only : lateral_boundary_diffusion_CS, lateral_boundary_diffusion_init +use MOM_lateral_boundary_diffusion, only : lateral_boundary_diffusion use MOM_tracer_registry, only : tracer_registry_type, tracer_type, MOM_tracer_chksum use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : thermo_var_ptrs @@ -57,7 +60,13 @@ module MOM_tracer_hor_diff !! the CFL limit is not violated. logical :: use_neutral_diffusion !< If true, use the neutral_diffusion module from within !! tracer_hor_diff. + logical :: use_lateral_boundary_diffusion !< If true, use the lateral_boundary_diffusion module from within + !! tracer_hor_diff. + logical :: recalc_neutral_surf !< If true, recalculate the neutral surfaces if CFL has been + !! exceeded type(neutral_diffusion_CS), pointer :: neutral_diffusion_CSp => NULL() !< Control structure for neutral diffusion. + type(lateral_boundary_diffusion_CS), pointer :: lateral_boundary_diffusion_CSp => NULL() !< Control structure for + !! lateral boundary mixing. type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to !! regulate the timing of diagnostic output. logical :: debug !< If true, write verbose checksums for debugging purposes. @@ -381,6 +390,31 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, G, GV, US, CS, Reg, tv, do_online endif enddo + if (CS%use_lateral_boundary_diffusion) then + + if (CS%show_call_tree) call callTree_waypoint("Calling lateral boundary mixing (tracer_hordiff)") + + call do_group_pass(CS%pass_t, G%Domain, clock=id_clock_pass) + + do J=js-1,je ; do i=is,ie + Coef_y(i,J) = I_numitts * khdt_y(i,J) + enddo ; enddo + do j=js,je + do I=is-1,ie + Coef_x(I,j) = I_numitts * khdt_x(I,j) + enddo + enddo + + do itt=1,num_itts + if (CS%show_call_tree) call callTree_waypoint("Calling lateral boundary diffusion (tracer_hordiff)",itt) + if (itt>1) then ! Update halos for subsequent iterations + call do_group_pass(CS%pass_t, G%Domain, clock=id_clock_pass) + endif + call lateral_boundary_diffusion(G, GV, US, h, Coef_x, Coef_y, I_numitts*dt, Reg, & + CS%lateral_boundary_diffusion_CSp) + enddo ! itt + endif + if (CS%use_neutral_diffusion) then if (CS%show_call_tree) call callTree_waypoint("Calling neutral diffusion coeffs (tracer_hordiff)") @@ -390,7 +424,7 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, G, GV, US, CS, Reg, tv, do_online ! lateral diffusion iterations. Otherwise the call to neutral_diffusion_calc_coeffs() ! would be inside the itt-loop. -AJA - call neutral_diffusion_calc_coeffs(G, GV, h, tv%T, tv%S, CS%neutral_diffusion_CSp) + call neutral_diffusion_calc_coeffs(G, GV, US, h, tv%T, tv%S, CS%neutral_diffusion_CSp) do J=js-1,je ; do i=is,ie Coef_y(i,J) = I_numitts * khdt_y(i,J) enddo ; enddo @@ -404,6 +438,9 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, G, GV, US, CS, Reg, tv, do_online if (CS%show_call_tree) call callTree_waypoint("Calling neutral diffusion (tracer_hordiff)",itt) if (itt>1) then ! Update halos for subsequent iterations call do_group_pass(CS%pass_t, G%Domain, clock=id_clock_pass) + if (CS%recalc_neutral_surf) then + call neutral_diffusion_calc_coeffs(G, GV, US, h, tv%T, tv%S, CS%neutral_diffusion_CSp) + endif endif call neutral_diffusion(G, GV, h, Coef_x, Coef_y, I_numitts*dt, Reg, US, CS%neutral_diffusion_CSp) enddo ! itt @@ -1383,12 +1420,13 @@ end subroutine tracer_epipycnal_ML_diff !> Initialize lateral tracer diffusion module -subroutine tracer_hor_diff_init(Time, G, US, param_file, diag, EOS, CS) +subroutine tracer_hor_diff_init(Time, G, US, param_file, diag, EOS, diabatic_CSp, CS) type(time_type), target, intent(in) :: Time !< current model time type(ocean_grid_type), intent(in) :: G !< ocean grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(diag_ctrl), target, intent(inout) :: diag !< diagnostic control type(EOS_type), target, intent(in) :: EOS !< Equation of state CS + type(diabatic_CS), pointer, intent(in) :: diabatic_CSp !< Equation of state CS type(param_file_type), intent(in) :: param_file !< parameter file type(tracer_hor_diff_CS), pointer :: CS !< horz diffusion control structure @@ -1446,6 +1484,10 @@ subroutine tracer_hor_diff_init(Time, G, US, param_file, diag, EOS, CS) "below this value. The number of diffusive iterations "//& "is often this value or the next greater integer.", & units="nondim", default=-1.0) + call get_param(param_File, mdl, "RECALC_NEUTRAL_SURF", CS%recalc_neutral_surf, & + "If true, then recalculate the neutral surfaces if the \n"//& + "diffusive CFL is exceeded. If false, assume that the \n"//& + "positions of the surfaces do not change \n", default = .false.) CS%ML_KhTR_scale = 1.0 if (CS%Diffuse_ML_interior) then call get_param(param_file, mdl, "ML_KHTR_SCALE", CS%ML_KhTR_scale, & @@ -1455,9 +1497,14 @@ subroutine tracer_hor_diff_init(Time, G, US, param_file, diag, EOS, CS) units="nondim", default=1.0) endif - CS%use_neutral_diffusion = neutral_diffusion_init(Time, G, param_file, diag, EOS, CS%neutral_diffusion_CSp) + CS%use_neutral_diffusion = neutral_diffusion_init(Time, G, param_file, diag, EOS, diabatic_CSp, & + CS%neutral_diffusion_CSp ) if (CS%use_neutral_diffusion .and. CS%Diffuse_ML_interior) call MOM_error(FATAL, "MOM_tracer_hor_diff: "// & "USE_NEUTRAL_DIFFUSION and DIFFUSE_ML_TO_INTERIOR are mutually exclusive!") + CS%use_lateral_boundary_diffusion = lateral_boundary_diffusion_init(Time, G, param_file, diag, diabatic_CSp, & + CS%lateral_boundary_diffusion_CSp) + if (CS%use_neutral_diffusion .and. CS%Diffuse_ML_interior) call MOM_error(FATAL, "MOM_tracer_hor_diff: "// & + "USE_LATERAL_BOUNDARY_DIFFUSION and DIFFUSE_ML_TO_INTERIOR are mutually exclusive!") call get_param(param_file, mdl, "DEBUG", CS%debug, default=.false.) diff --git a/src/tracer/MOM_tracer_registry.F90 b/src/tracer/MOM_tracer_registry.F90 index 4680c058b4..f71760d361 100644 --- a/src/tracer/MOM_tracer_registry.F90 +++ b/src/tracer/MOM_tracer_registry.F90 @@ -56,6 +56,18 @@ module MOM_tracer_registry !! [conc H m2 s-1 ~> conc m3 s-1 or conc kg s-1] real, dimension(:,:,:), pointer :: df_y => NULL() !< diagnostic array for y-diffusive tracer flux !! [conc H m2 s-1 ~> conc m3 s-1 or conc kg s-1] + real, dimension(:,:,:), pointer :: lbd_dfx => NULL() !< diagnostic array for x-diffusive tracer flux + !! [conc H m2 s-1 ~> conc m3 s-1 or conc kg s-1] + real, dimension(:,:,:), pointer :: lbd_dfy => NULL() !< diagnostic array for y-diffusive tracer flux + !! [conc H m2 s-1 ~> conc m3 s-1 or conc kg s-1] + real, dimension(:,:), pointer :: lbd_dfx_2d => NULL() !< diagnostic array for x-diffusive tracer flux + !! [conc H m2 s-1 ~> conc m3 s-1 or conc kg s-1] + real, dimension(:,:), pointer :: lbd_dfy_2d => NULL() !< diagnostic array for y-diffusive tracer flux + !! [conc H m2 s-1 ~> conc m3 s-1 or conc kg s-1] + real, dimension(:,:), pointer :: lbd_bulk_df_x => NULL() !< diagnostic array for x-diffusive tracer flux + !! [conc H m2 s-1 ~> conc m3 s-1 or conc kg s-1] + real, dimension(:,:), pointer :: lbd_bulk_df_y => NULL() !< diagnostic array for y-diffusive tracer flux + !! [conc H m2 s-1 ~> conc m3 s-1 or conc kg s-1] real, dimension(:,:), pointer :: df2d_x => NULL() !< diagnostic vertical sum x-diffusive flux !! [conc H m2 s-1 ~> conc m3 s-1 or conc kg s-1] real, dimension(:,:), pointer :: df2d_y => NULL() !< diagnostic vertical sum y-diffusive flux @@ -109,6 +121,8 @@ module MOM_tracer_registry !>@{ Diagnostic IDs integer :: id_tr = -1 integer :: id_adx = -1, id_ady = -1, id_dfx = -1, id_dfy = -1 + integer :: id_lbd_bulk_dfx = -1, id_lbd_bulk_dfy = -1, id_lbd_dfx = -1, id_lbd_dfy = -1 + integer :: id_lbd_dfx_2d = -1 , id_lbd_dfy_2d = -1 integer :: id_adx_2d = -1, id_ady_2d = -1, id_dfx_2d = -1, id_dfy_2d = -1 integer :: id_adv_xy = -1, id_adv_xy_2d = -1 integer :: id_dfxy_cont = -1, id_dfxy_cont_2d = -1, id_dfxy_conc = -1 @@ -398,8 +412,14 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, use_ALE) diag%axesCuL, Time, trim(flux_longname)//" diffusive zonal flux" , & trim(flux_units), v_extensive = .true., y_cell_method = 'sum') Tr%id_dfy = register_diag_field("ocean_model", trim(shortnm)//"_dfy", & - diag%axesCvL, Time, trim(flux_longname)//" diffusive zonal flux" , & + diag%axesCvL, Time, trim(flux_longname)//" diffusive merdional flux" , & trim(flux_units), v_extensive = .true., x_cell_method = 'sum') + Tr%id_lbd_dfx = register_diag_field("ocean_model", trim(shortnm)//"_lbd_diffx", & + diag%axesCuL, Time, trim(flux_longname)//" diffusive zonal flux from the lateral boundary diffusion "& + "scheme", trim(flux_units), v_extensive = .true., y_cell_method = 'sum') + Tr%id_lbd_dfy = register_diag_field("ocean_model", trim(shortnm)//"_lbd_diffy", & + diag%axesCvL, Time, trim(flux_longname)//" diffusive meridional flux from the lateral boundary diffusion"& + " scheme", trim(flux_units), v_extensive = .true., x_cell_method = 'sum') else Tr%id_adx = register_diag_field("ocean_model", trim(shortnm)//"_adx", & diag%axesCuL, Time, "Advective (by residual mean) Zonal Flux of "//trim(flux_longname), & @@ -413,11 +433,19 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, use_ALE) Tr%id_dfy = register_diag_field("ocean_model", trim(shortnm)//"_diffy", & diag%axesCvL, Time, "Diffusive Meridional Flux of "//trim(flux_longname), & flux_units, v_extensive=.true., conversion=Tr%flux_scale, x_cell_method = 'sum') + Tr%id_lbd_dfx = register_diag_field("ocean_model", trim(shortnm)//"_lbd_diffx", & + diag%axesCuL, Time, "Lateral Boundary Diffusive Zonal Flux of "//trim(flux_longname), & + flux_units, v_extensive=.true., conversion=Tr%flux_scale, y_cell_method = 'sum') + Tr%id_lbd_dfy = register_diag_field("ocean_model", trim(shortnm)//"_lbd_diffy", & + diag%axesCvL, Time, "Lateral Boundary Diffusive Meridional Flux of "//trim(flux_longname), & + flux_units, v_extensive=.true., conversion=Tr%flux_scale, x_cell_method = 'sum') endif if (Tr%id_adx > 0) call safe_alloc_ptr(Tr%ad_x,IsdB,IedB,jsd,jed,nz) if (Tr%id_ady > 0) call safe_alloc_ptr(Tr%ad_y,isd,ied,JsdB,JedB,nz) if (Tr%id_dfx > 0) call safe_alloc_ptr(Tr%df_x,IsdB,IedB,jsd,jed,nz) if (Tr%id_dfy > 0) call safe_alloc_ptr(Tr%df_y,isd,ied,JsdB,JedB,nz) + if (Tr%id_lbd_dfx > 0) call safe_alloc_ptr(Tr%lbd_dfx,IsdB,IedB,jsd,jed,nz) + if (Tr%id_lbd_dfy > 0) call safe_alloc_ptr(Tr%lbd_dfy,isd,ied,JsdB,JedB,nz) Tr%id_adx_2d = register_diag_field("ocean_model", trim(shortnm)//"_adx_2d", & diag%axesCu1, Time, & @@ -435,11 +463,29 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, use_ALE) diag%axesCv1, Time, & "Vertically Integrated Diffusive Meridional Flux of "//trim(flux_longname), & flux_units, conversion=Tr%flux_scale, x_cell_method = 'sum') + Tr%id_lbd_bulk_dfx = register_diag_field("ocean_model", trim(shortnm)//"_lbd_bulk_diffx", & + diag%axesCu1, Time, & + "Total Bulk Diffusive Zonal Flux of "//trim(flux_longname), & + flux_units, conversion=Tr%flux_scale, y_cell_method = 'sum') + Tr%id_lbd_bulk_dfy = register_diag_field("ocean_model", trim(shortnm)//"_lbd_bulk_diffy", & + diag%axesCv1, Time, & + "Total Bulk Diffusive Meridional Flux of "//trim(flux_longname), & + flux_units, conversion=Tr%flux_scale, x_cell_method = 'sum') + Tr%id_lbd_dfx_2d = register_diag_field("ocean_model", trim(shortnm)//"_lbd_diffx_2d", & + diag%axesCu1, Time, "Vertically-integrated zonal diffusive flux from the lateral boundary diffusion "//& + "scheme for "//trim(flux_longname), flux_units, conversion=Tr%flux_scale, y_cell_method = 'sum') + Tr%id_lbd_dfy_2d = register_diag_field("ocean_model", trim(shortnm)//"_lbd_diffy_2d", & + diag%axesCv1, Time, "Vertically-integrated meridional diffusive flux from the lateral boundary diffusion "//& + "scheme for "//trim(flux_longname), flux_units, conversion=Tr%flux_scale, x_cell_method = 'sum') if (Tr%id_adx_2d > 0) call safe_alloc_ptr(Tr%ad2d_x,IsdB,IedB,jsd,jed) if (Tr%id_ady_2d > 0) call safe_alloc_ptr(Tr%ad2d_y,isd,ied,JsdB,JedB) if (Tr%id_dfx_2d > 0) call safe_alloc_ptr(Tr%df2d_x,IsdB,IedB,jsd,jed) if (Tr%id_dfy_2d > 0) call safe_alloc_ptr(Tr%df2d_y,isd,ied,JsdB,JedB) + if (Tr%id_lbd_bulk_dfx > 0) call safe_alloc_ptr(Tr%lbd_bulk_df_x,IsdB,IedB,jsd,jed) + if (Tr%id_lbd_bulk_dfy > 0) call safe_alloc_ptr(Tr%lbd_bulk_df_y,isd,ied,JsdB,JedB) + if (Tr%id_lbd_dfx_2d > 0) call safe_alloc_ptr(Tr%lbd_dfx_2d,IsdB,IedB,jsd,jed) + if (Tr%id_lbd_dfy_2d > 0) call safe_alloc_ptr(Tr%lbd_dfy_2d,isd,ied,JsdB,JedB) Tr%id_adv_xy = register_diag_field('ocean_model', trim(shortnm)//"_advection_xy", & diag%axesTL, Time, &