diagnistics module created

CMakeLists.txt

@@ -22,6 +22,7 @@ set(SOURCES
     obl_forcing_and_boundary.f90
     obl_rhs.f90
     obl_richardson.f90
+    obl_diag.f90
     obl_k_epsilon.f90
     obl_pacanowski_philander.f90
     obl_pacanowski_philander_plus.f90

obl_diag.f90

+#include "obl_def.fi"
+module obl_diag
+    !< @brief diagnostics module
+    ! --------------------------------------------------------------------------------
+
+    ! modules used
+    ! --------------------------------------------------------------------------------
+
+    ! modules used
+    use obl_grid
+    use obl_math
+    use obl_common_phys_parameters
+    use obl_state_eq 
+
+    ! directives list
+    implicit none 
+    private
+
+    ! public interface
+    ! --------------------------------------------------------------------------------
+    public :: get_hml, get_lab_hml
+    ! --------------------------------------------------------------------------------
+    
+    contains
+
+    subroutine get_hml(hml, maxcellnumber, N2, dz, nz, z)
+        !< @brief calculate mixed layer depth 
+        ! ----------------------------------------------------------------------------
+        real, intent(out) :: hml                    !< mixed layer depth, [m]
+        integer, intent(out) :: maxcellnumber       !< index of cell with maximum N2
+
+        integer, intent(in) :: nz                   !< number of z-steps
+        real, intent(in), dimension(nz) :: N2       !< square of buoyancy frequency, [s**(-2)]
+        real, intent(in), dimension(nz) :: dz       !< depth(z) step [m] 
+        real, intent(in) :: z                       !< depth, [m]
+        
+        real :: max_N2, hml_tmp !< temporary
+        real :: hml_min, hml_max !< hml limits
+        integer :: k !< counter
+        ! ----------------------------------------------------------------------------
+        hml_min = 0.5 * dz(nz)
+        hml_max = z
+
+        max_N2 = N2(nz)
+        maxcellnumber = 1
+        hml = 0.0
+        hml_tmp = 0.0
+
+        do k=nz - 1, 1, -1
+            if (N2(k) > max_N2) then
+                max_N2 = N2(k)
+                maxcellnumber = k
+            end if
+        end do 
+
+        do k=nz, maxcellnumber + 1, -1
+            hml_tmp = hml_tmp + dz(k)
+        enddo
+        hml_tmp = hml_tmp + 0.5 * dz(maxcellnumber)
+        
+        hml = hml_tmp 
+
+        if (hml < hml_min) hml = hml_min
+        if (hml > hml_max) hml = hml_max
+
+    end subroutine get_hml
+
+    subroutine get_lab_hml(lab_hml, flux_u_surf, flux_v_surf, time, z)
+        !< @brief calculate theoretical mixed layer depth(Price, 1970)
+        ! ----------------------------------------------------------------------------
+
+        real, intent(out) :: lab_hml                    !< theoretical mixed layer depth, [m]
+
+        real, intent(in) :: flux_u_surf, flux_v_surf    !< u_dyn**2 & v_dyn**2, [(m/s)**2]
+        real, intent(in) :: time                        !< time, [s]
+        real, intent(in) :: z                           !<depth, [m]
+        ! ----------------------------------------------------------------------------
+        lab_hml = 1.05  * (flux_u_surf**(2.0)+flux_v_surf**(2.0))**(1.0/4.0) * sqrt(time) / ((0.0004)**(1.0/4.0))
+
+        if (lab_hml > z) lab_hml = z
+
+    end subroutine
+
+endmodule
\ No newline at end of file

obl_main.f90

@@ -19,6 +19,7 @@ program obl_main
     use obl_rhs
     use obl_equation
     use obl_richardson
+    use obl_diag
     use obl_pacanowski_philander
     use obl_pacanowski_philander_plus
     use obl_k_epsilon

obl_richardson.f90

@@ -184,57 +184,4 @@ module obl_richardson
         enddo
     end subroutine
 
-    subroutine get_hml(hml, maxcellnumber, N2, dz, nz, z)
-        !< @brief calculate mixed layer depth 
-        integer, intent(in) :: nz !< number of z-steps
-        real, intent(out) :: hml !< mixed layer depth, [m]
-        real, intent(in) :: N2(nz) !< square of buoyancy frequency, [s**(-2)]
-        real, intent(in) :: z  !< depth, [m]
-        integer :: k !< counter
-        real :: max_N2, hml_tmp !< temporary
-        integer, intent(out) :: maxcellnumber !< index of cell with maximum N2
-        real, intent(in) :: dz(nz) !< depth(z) step [m] 
-        real :: hml_min, hml_max !< hml limits
-
-        hml_min = 0.5 * dz(nz)
-        hml_max = z
-
-        max_N2 = N2(nz)
-        maxcellnumber = 1
-        hml = 0.0
-        hml_tmp = 0.0
-
-        do k=nz - 1, 1, -1
-            if (N2(k) > max_N2) then
-                max_N2 = N2(k)
-                maxcellnumber = k
-            end if
-        end do 
-
-        do k=nz, maxcellnumber + 1, -1
-            hml_tmp = hml_tmp + dz(k)
-        enddo
-        hml_tmp = hml_tmp + 0.5 * dz(maxcellnumber)
-        
-        hml = hml_tmp!z - hml_tmp
-
-        if (hml < hml_min) hml = hml_min
-        if (hml > hml_max) hml = hml_max
-
-    end subroutine
-
-    subroutine get_lab_hml(lab_hml, flux_u_surf, flux_v_surf, time, z)
-        !< @brief calculate theoretical mixed layer depth(Price, 1970)
-        real, intent(out) :: lab_hml !< theoretical mixed layer depth, [m]
-        real :: flux_u_surf, flux_v_surf !< u_dyn**2 & v_dyn**2, [(m/s)**2]
-        real :: time !< time, [s]
-        real :: z !<depth, [m]
-
-        lab_hml = 1.05  * (flux_u_surf**(2.0)+flux_v_surf**(2.0))**(1.0/4.0) * sqrt(time) / ((0.0004)**(1.0/4.0))
-
-        if (lab_hml > z) lab_hml = z
-            
-
-    end subroutine
-
 endmodule
\ No newline at end of file