adding common module; command line arguments in main

compile.sh

 #!/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

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)
 

srcF/sfx_common.f90

+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

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

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
 

srcF/sfx_main.f90

 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))