diff --git a/compile.sh b/compile.sh
index f0c5bf9cb7bfee4f97dd848fe4fd32af4e1881a5..107662793561f45fcded8574cf3e2f2f41fa07fb 100755
--- a/compile.sh
+++ b/compile.sh
@@ -1,8 +1,9 @@
#!/bin/bash
gfortran -c -cpp -Wuninitialized srcF/sfx_phys_const.f90
+gfortran -c -cpp -Wuninitialized srcF/sfx_common.f90
gfortran -c -cpp -Wuninitialized srcF/sfx_esm_param.f90
gfortran -c -cpp -Wuninitialized srcF/sfx_esm.f90
gfortran -c -cpp -Wuninitialized srcF/sfx_main.f90
-gfortran -o sfx.exe sfx_main.o sfx_esm.o sfx_esm_param.o sfx_phys_const.o
+gfortran -o sfx.exe sfx_main.o sfx_esm.o sfx_esm_param.o sfx_common.o sfx_phys_const.o
rm *.o *.mod
diff --git a/makefile b/makefile
index 4a0e47b2bdf5c9eda1ea57e19134ce17b5dfd3a5..05e2436e4687416b83d0c1a259264bc626949fed 100644
--- a/makefile
+++ b/makefile
@@ -10,7 +10,7 @@ ifeq ($(COMPILER),gnu)
FC = gfortran
endif
-OBJ_F90 = sfx_phys_const.o sfx_esm_param.o sfx_esm.o sfx_main.o
+OBJ_F90 = sfx_phys_const.o sfx_common.o sfx_esm_param.o sfx_esm.o sfx_main.o
OBJ_F =
OBJ = $(OBJ_F90) $(OBJ_F)
diff --git a/srcF/sfx_common.f90 b/srcF/sfx_common.f90
new file mode 100644
index 0000000000000000000000000000000000000000..e7852ff79e402d809c5804c6167ff7dbdedf0d63
--- /dev/null
+++ b/srcF/sfx_common.f90
@@ -0,0 +1,31 @@
+module sfx_common
+ !> @brief surface flux code common subroutines
+ public
+
+contains
+
+ ! ----------------------------------------------------------------------------
+ elemental subroutine str2int(int, str, stat)
+ !> @brief string to int conversion
+ ! ----------------------------------------------------------------------------
+ implicit none
+ integer, intent(out) :: int
+ integer, intent(out) :: stat
+
+ character(len = *), intent(in) :: str
+ ! ----------------------------------------------------------------------------
+
+ read(str, * , iostat=stat) int
+ end subroutine str2int
+ ! ----------------------------------------------------------------------------
+
+ elemental function is_finite(value)
+ use ieee_arithmetic
+ implicit none
+ logical :: is_finite
+ real, intent(in) :: value
+
+ is_finite = ieee_is_finite(value)
+ end function is_finite
+
+end module sfx_common
\ No newline at end of file
diff --git a/srcF/sfx_esm.f90 b/srcF/sfx_esm.f90
index 188ab3845d1719ab4dea361674c8c6f8b64a1a17..a611bd2e4a159e88e4b914e0b6062e31b4db2833 100644
--- a/srcF/sfx_esm.f90
+++ b/srcF/sfx_esm.f90
@@ -3,11 +3,15 @@ module sfx_esm
! macro defs.
! --------------------------------------------------------------------------------
-#define SFX_FORCE_DEPRECATED_CODE
+!#define SFX_FORCE_DEPRECATED_CODE
+!#define SFX_CHECK_NAN
! --------------------------------------------------------------------------------
! modules used
! --------------------------------------------------------------------------------
+#ifdef SFX_CHECK_NAN
+ use sfx_common
+#endif
use sfx_esm_param
! --------------------------------------------------------------------------------
@@ -111,7 +115,7 @@ contains
! ----------------------------------------------------------------------------
do i = 1, n
-#ifdef USE_DEPRECATED_CODE
+#ifdef SFX_FORCE_DEPRECATED_CODE
#else
meteo_cell = meteoDataType(&
h = meteo%h(i), &
@@ -156,6 +160,10 @@ contains
!> @brief surface flux calculation for single cell
!> @details contains C/C++ interface
! ----------------------------------------------------------------------------
+#ifdef SFX_CHECK_NAN
+ use ieee_arithmetic
+#endif
+
type (sfxDataType), intent(out) :: sfx
type (meteoDataType), intent(in) :: meteo
@@ -201,8 +209,26 @@ contains
integer surface_type !> surface type = (ocean || land)
real fval !> just a shortcut for partial calculations
+
+#ifdef SFX_CHECK_NAN
+ real NaN
+#endif
! ----------------------------------------------------------------------------
+#ifdef SFX_CHECK_NAN
+ ! --- checking if arguments are finite
+ if (.not.(is_finite(meteo%U).and.is_finite(meteo%Tsemi).and.is_finite(meteo%dT).and.is_finite(meteo%dQ) &
+ .and.is_finite(meteo%z0_m).and.is_finite(meteo%h))) then
+
+ NaN = ieee_value(0.0, ieee_quiet_nan) ! setting NaN
+ sfx = sfxDataType(zeta = NaN, Rib = NaN, &
+ Re = NaN, B = NaN, z0_m = NaN, z0_t = NaN, &
+ Rib_conv_lim = NaN, &
+ Cm = NaN, Ct = NaN, Km = NaN, Pr_t_inv = NaN)
+ return
+ end if
+#endif
+
! --- shadowing names for clarity
U = meteo%U
Tsemi = meteo%Tsemi
diff --git a/srcF/sfx_esm_param.f90 b/srcF/sfx_esm_param.f90
index fad0d7588e2816aa039b475063330174c6619b78..3aec59df790a95e63a61dec30cc353c1a5ae24cd 100644
--- a/srcF/sfx_esm_param.f90
+++ b/srcF/sfx_esm_param.f90
@@ -12,6 +12,7 @@ module sfx_esm_param
implicit none
! --------------------------------------------------------------------------------
+ ! *: add lake surface & add surface type as input value
integer, public, parameter :: surface_ocean = 0 !> ocean surface
integer, public, parameter :: surface_land = 1 !> land surface
diff --git a/srcF/sfx_main.f90 b/srcF/sfx_main.f90
index ad5f42a2de26639300e9e34ad94fe38f60ef08e7..96a435cf954d1dd8e4535822c0d6b716c994fe0c 100644
--- a/srcF/sfx_main.f90
+++ b/srcF/sfx_main.f90
@@ -1,118 +1,259 @@
program sfx_main
+ ! modules used
+ ! --------------------------------------------------------------------------------
use sfx_phys_const
+ use sfx_common
+
use sfx_esm_param
use sfx_esm
+ ! --------------------------------------------------------------------------------
+
+ ! directives list
+ ! --------------------------------------------------------------------------------
+ implicit none
+ ! --------------------------------------------------------------------------------
+
+ !> dataset ID:
+ !> = 0, USER
+ !> = 1, MOSAiC dataset
+ !> = 2, IRGASON dataset
+ !> = 3, SHEBA dataset
+ integer :: dataset_id
+ character(len = 256) :: dataset_name
+ integer, parameter :: dataset_MOSAiC = 1
+ integer, parameter :: dataset_IRGASON = 2
+ integer, parameter :: dataset_SHEBA = 3 !> please spell 'SHIBA'
+ integer, parameter :: dataset_USER = 4 !> used defined dataset
+
+ ! input/output data
+ ! --------------------------------------------------------------------------------
+ type(meteoDataVecType) :: meteo !> meteorological data (input)
+ type(meteoDataType) :: meteo_cell
+
+ type(sfxDataVecType) :: sfx !> surface fluxes (output)
+
+ type(numericsType) :: numerics !> surface flux module numerics parameters
+
+ integer :: num !> number of 'cells' in input
+
+
+ ! --- input/output filenames
+ character(len = 256) :: filename_in_common
+ character(len = 256) :: filename_in
+ character(len = 256) :: filename_out
+ ! --------------------------------------------------------------------------------
+
+ ! command line arguments
+ ! --------------------------------------------------------------------------------
+ integer :: num_args
+ character(len = 128) :: arg
+
+ character(len = 128), parameter :: arg_key_dataset = '--dataset'
+ character(len = 128), parameter :: arg_key_output = '--output'
+ character(len = 128), parameter :: arg_key_nmax = '--nmax'
+ character(len = 128), parameter :: arg_key_help = '--help'
+
+ character(len = 128), parameter :: arg_key_mosaic = 'mosaic'
+ character(len = 128), parameter :: arg_key_irgason = 'irgason'
+ character(len = 128), parameter :: arg_key_sheba = 'sheba'
+ character(len = 128), parameter :: arg_key_user = 'user'
+
+ integer :: is_output_set
+ integer :: nmax
+ ! --------------------------------------------------------------------------------
+
+ ! local variables
+ ! --------------------------------------------------------------------------------
+ integer :: i
+ integer :: status
+ ! --------------------------------------------------------------------------------
+
+
+ !> @brief define dataset
+ dataset_id = dataset_MOSAiC !> default = MOSAiC
+
+ is_output_set = 0
+ nmax = 0
+ num_args = command_argument_count()
+ do i = 1, num_args
+ call get_command_argument(i, arg)
+ if (trim(arg) == trim(arg_key_help)) then
+ write(*, *) ' sfx model, usage:'
+ write(*, *) ' --help '
+ write(*, *) ' print usage options '
+ write(*, *) ' --dataset [key]'
+ write(*, *) ' key = mosaic || irgason || sheba || user [files]'
+ write(*, *) ' files = in-common-file in-file out-file'
+ write(*, *) ' --output [file]'
+ write(*, *) ' set output filename '
+ write(*, *) ' --nmax [value]'
+ write(*, *) ' max number of data points > 0 '
+ stop
+ end if
+ if (trim(arg) == trim(arg_key_dataset)) then
+ if (i == num_args) then
+ write(*, *) ' FAILURE! > missing dataset [key] argument'
+ stop
+ end if
+ call get_command_argument(i + 1, arg)
+ if (trim(arg) == trim(arg_key_mosaic)) then
+ dataset_id = dataset_MOSAiC
+ else if (trim(arg) == trim(arg_key_irgason)) then
+ dataset_id = dataset_IRGASON
+ else if (trim(arg) == trim(arg_key_sheba)) then
+ dataset_id = dataset_SHEBA
+ else if (trim(arg) == trim(arg_key_user)) then
+ dataset_id = dataset_USER
+ if (i + 4 > num_args) then
+ write(*, *) ' FAILURE! > incorrect arguments for [user] dataset'
+ stop
+ end if
+ call get_command_argument(i + 2, filename_in_common)
+ call get_command_argument(i + 3, filename_in)
+ call get_command_argument(i + 4, filename_out)
+ else
+ write(*, *) ' FAILURE! > unknown dataset [key]: ', trim(arg)
+ stop
+ end if
+ end if
+ if (trim(arg) == trim(arg_key_output)) then
+ if (i == num_args) then
+ write(*, *) ' FAILURE! > missing dataset [key] argument'
+ stop
+ end if
+ is_output_set = 1
+ call get_command_argument(i + 1, filename_out)
+ end if
+ if (trim(arg) == trim(arg_key_nmax)) then
+ if (i == num_args) then
+ write(*, *) ' FAILURE! > missing nmax [key] argument'
+ stop
+ end if
+ call get_command_argument(i + 1, arg)
+ call str2int(nmax, arg, status)
+ if (status /= 0) then
+ write(*, *) ' FAILURE! > expecting int nmax [value]'
+ stop
+ end if
+ if (nmax <= 0) then
+ write(*, *) ' FAILURE! > nmax [value] should be positive'
+ stop
+ end if
+ end if
+ end do
+
+
+ !> @brief set filenames & data for specific dataset
+ if (dataset_id == dataset_MOSAiC) then
+ dataset_name = 'MOSAiC'
- integer, parameter :: test = 1
-
-
- type(meteoDataType):: data_in1
- type(meteoDataVecType) :: meteo
- type(sfxDataType) :: data_outdef1
- type(sfxDataVecType) :: data_outMAS
- type(numericsType) :: data_par1
-
-
- integer :: numst, i
- real :: cflh, z0in
- character(len = 50) :: filename_in
- character(len = 50) :: filename_out
- character(len = 50) :: filename_in2
-
- !input
- ! mas_w - abs(wind velocity) at constant flux layer (cfl) hight (m/s)
- ! mas_dt - difference between potential temperature at cfl hight and at surface ( deg. k)
- ! mas_st - semi-sum of potential temperature at cfl hight and and at surface ( deg. k)
- ! mas_dq - difference between humidity at cfl hight and a surface ( gr/gr )
- ! cflh - - cfl hight ( m )
- ! z0in=0.01 - roughness of surface ( m );
- ! it - number of iterations
- ! lu_indx - 1 for land, 2 for sea, 3 for lake
- ! test - file input
-
-
- !output
- !masout_zl - non-dimensional cfl hight
- !masout_ri - richardson number
- !masout_re - reynods number
- !masout_lnzuzt - ln(zu/zt)
- !masout_zu - dynamical roughness zu (m)
- !masout_ztout - thermal roughness zt (m)
- !masout_rith - critical richardson number
- !masout_cm - transfer coefficient for momentum
- !masout_ch - transfer coefficient fr heat
- !masout_ct - coefficient of turbulence (km) at cfl hight (m**2/s)
- !masout_ckt - alft=kt/km ( kt-coefficient of turbulence for heat)
-
- !> @brief Test - file selection for test
-
- write(*,*) 'running code'
-
- if (TEST==1) then
- filename_in='data/MOSAiC.txt'
- filename_out='out_MOSAiC.txt'
- filename_in2='data/MOSAiC_zh.txt'
- elseif (TEST==2) then
- filename_in='data/Irgason1.txt'
- filename_out='out_IRGASON1.txt'
- filename_in2='data/IRGASON_zh.txt'
- endif
-
- open (1, file= filename_in, status ='old')
- open (2, file=filename_out)
- numst=0
- do WHILE (ioer.eq.0)
- read (1,*, iostat=ioer) data_in1%U, data_in1%dT, data_in1%Tsemi, data_in1%dQ
- numst=numst+1
+ filename_in_common = 'data/MOSAiC_zh.txt'
+ filename_in = 'data/MOSAiC.txt'
+ if (is_output_set == 0) filename_out = 'out_MOSAiC.txt'
+ else if (dataset_id == dataset_IRGASON) then
+ dataset_name = 'IRGASON'
+
+ filename_in_common = 'data/IRGASON_zh.txt'
+ filename_in = 'data/Irgason1.txt'
+ if (is_output_set == 0) filename_out = 'out_IRGASON1.txt'
+ else if (dataset_id == dataset_SHEBA) then
+ dataset_name = 'SHEBA'
+
+ write(*, *) ' FAILURE! > SHEBA dataset is not supported yet:( '
+ stop
+ else if (dataset_id == dataset_USER) then
+ dataset_name = 'USER'
+ else
+ write(*, *) ' FAILURE! > unknown dataset id: ', dataset_id
+ end if
+
+ write(*, *) ' Running SFX model'
+ write(*, *) ' dataset = ', trim(dataset_name)
+ write(*, *) ' filename[IN-COMMON] = ', trim(filename_in_common)
+ write(*, *) ' filename[IN] = ', trim(filename_in)
+ write(*, *) ' filename[OUT] = ', trim(filename_out)
+
+
+ !> @brief define number of cells
+ open(1, file= filename_in, status ='old')
+ status = 0
+
+ num = 0
+ do while (status.eq.0)
+ read (1, *, iostat = status) meteo_cell%U, meteo_cell%dT, meteo_cell%Tsemi, meteo_cell%dQ
+ num = num + 1
enddo
- close (1)
- numst=numst-1
-
-
- allocate(meteo%h(numst))
- allocate(meteo%U(numst))
- allocate(meteo%dT(numst))
- allocate(meteo%Tsemi(numst))
- allocate(meteo%dQ(numst))
- allocate(meteo%z0_m(numst))
-
-
- allocate(data_outMAS%zeta(numst))
- allocate(data_outMAS%Rib(numst))
- allocate(data_outMAS%Re(numst))
- allocate(data_outMAS%B(numst))
- allocate(data_outMAS%z0_m(numst))
- allocate(data_outMAS%z0_t(numst))
- allocate(data_outMAS%Rib_conv_lim(numst))
- allocate(data_outMAS%Cm(numst))
- allocate(data_outMAS%Ct(numst))
- allocate(data_outMAS%Km(numst))
- allocate(data_outMAS%Pr_t_inv(numst))
-
- open (11, file=filename_in2, status ='old')
- open (1, file= filename_in, status ='old')
- read (11, *) cflh, z0in
- do i=1,numst
- read (1,*) data_in1%U, data_in1%dT, data_in1%Tsemi, data_in1%dQ
-
- meteo%h(i)=cflh
- meteo%U(i) = meteo%U(i)+data_in1%U
- meteo%dT(i) = meteo%dT(i)+data_in1%dT
- meteo%Tsemi(i) = meteo%Tsemi(i)+data_in1%Tsemi
- meteo%dQ(i) = meteo%dQ(i)+data_in1%dQ
- meteo%z0_m(i)=z0in
+ num = num - 1
+
+ close(1)
+
+ ! --- print number of elements in dataset
+ write(*, *) ' size = ', num
+ if (nmax > 0) then
+ write(*, *) ' nmax = ', nmax
+ num = min(num, nmax)
+ end if
+
+
+ !> @brief allocate input & output data
+ allocate(meteo%h(num))
+ allocate(meteo%U(num))
+ allocate(meteo%dT(num))
+ allocate(meteo%Tsemi(num))
+ allocate(meteo%dQ(num))
+ allocate(meteo%z0_m(num))
+
+ allocate(sfx%zeta(num))
+ allocate(sfx%Rib(num))
+ allocate(sfx%Re(num))
+ allocate(sfx%B(num))
+ allocate(sfx%z0_m(num))
+ allocate(sfx%z0_t(num))
+ allocate(sfx%Rib_conv_lim(num))
+ allocate(sfx%Cm(num))
+ allocate(sfx%Ct(num))
+ allocate(sfx%Km(num))
+ allocate(sfx%Pr_t_inv(num))
+
+
+ !> @brief read input data common parameters
+ open(1, file = filename_in_common, status = 'old')
+ read(1, *) meteo_cell%h, meteo_cell%z0_m
+ close(1)
+
+
+ !> @brief read input data
+ open(1, file = filename_in, status = 'old')
+ do i = 1, num
+ read(1, *) meteo_cell%U, meteo_cell%dT, meteo_cell%Tsemi, meteo_cell%dQ
+
+ meteo%h(i) = meteo_cell%h
+ meteo%U(i) = meteo_cell%U
+ meteo%dT(i) = meteo_cell%dT
+ meteo%Tsemi(i) = meteo_cell%Tsemi
+ meteo%dQ(i) = meteo_cell%dQ
+ meteo%z0_m(i) = meteo_cell%z0_m
enddo
+ close(1)
+
- CALL get_surface_fluxes_vec(data_outMAS, meteo, &
- data_par1, numst)
+ !> @brief calling flux module
+ call get_surface_fluxes_vec(sfx, meteo, numerics, num)
- do i=1,numst
- write (2,20) data_outMAS%zeta(i), data_outMAS%Rib(i), data_outMAS%Re(i), data_outMAS%B(i),&
- data_outMAS%z0_m(i), data_outMAS%z0_t(i), data_outMAS%Rib_conv_lim(i), data_outMAS%Cm(i),&
- data_outMAS%Ct(i), data_outMAS%Km(i), data_outMAS%Pr_t_inv(i)
+
+ !> @brief write output data
+ open(2, file = filename_out)
+ do i = 1, num
+ write(2, 20) sfx%zeta(i), sfx%Rib(i), &
+ sfx%Re(i), sfx%B(i), sfx%z0_m(i), sfx%z0_t(i), &
+ sfx%Rib_conv_lim(i), &
+ sfx%Cm(i),sfx%Ct(i), sfx%Km(i), sfx%Pr_t_inv(i)
enddo
+ close(2)
+
+ !> @brief deallocate input & output data
deallocate(meteo%h)
deallocate(meteo%U)
deallocate(meteo%dT)
@@ -120,18 +261,19 @@ program sfx_main
deallocate(meteo%dQ)
deallocate(meteo%z0_m)
- deallocate(data_outMAS%zeta)
- deallocate(data_outMAS%Rib)
- deallocate(data_outMAS%Re)
- deallocate(data_outMAS%B)
- deallocate(data_outMAS%z0_m)
- deallocate(data_outMAS%z0_t)
- deallocate(data_outMAS%Rib_conv_lim)
- deallocate(data_outMAS%Cm)
- deallocate(data_outMAS%Ct)
- deallocate(data_outMAS%Km)
- deallocate(data_outMAS%Pr_t_inv)
+ deallocate(sfx%zeta)
+ deallocate(sfx%Rib)
+ deallocate(sfx%Re)
+ deallocate(sfx%B)
+ deallocate(sfx%z0_m)
+ deallocate(sfx%z0_t)
+ deallocate(sfx%Rib_conv_lim)
+ deallocate(sfx%Cm)
+ deallocate(sfx%Ct)
+ deallocate(sfx%Km)
+ deallocate(sfx%Pr_t_inv)
+ ! *: remove format(10) if not needed
10 format (f8.4,2x,f8.4)
20 format (11(f10.4,3x))