Changeset 955

Timestamp:

07/15/19 23:21:31 (5 years ago)

Author:

dubos

Message:

devel : fixed diagnosis of hydrostatic pressure for Lagrangian vertical coordinate

Location:

codes/icosagcm/devel/src

Files:

• base/init_grid_param.f90 (modified) (2 diffs)
• diagnostics/compute_diagnostics.f90 (modified) (2 diffs)
• diagnostics/compute_pression.F90 (modified) (5 diffs)
• diagnostics/observable.f90 (modified) (2 diffs)
• dynamics/compute_geopot.F90 (modified) (2 diffs)

codes/icosagcm/devel/src/base/init_grid_param.f90

@@ -70 +70 @@
     compute_pression_mid => compute_pression_mid_hex
     compute_temperature  => compute_temperature_hex
+    compute_hydrostatic_pressure => compute_hydrostatic_pressure_hex
     ! dynamics
     compute_pvort_only        => compute_pvort_only_hex
@@ -85 +86 @@
     compute_pression_mid => compute_pression_mid_unst
     compute_temperature  => compute_temperature_unst
+    compute_hydrostatic_pressure => compute_hydrostatic_pressure_unst
     ! dynamics
     compute_pvort_only         => compute_pvort_only_unst

codes/icosagcm/devel/src/diagnostics/compute_diagnostics.f90

@@ -31 +31 @@
     END SUBROUTINE comp_temperature
 
+    SUBROUTINE comp_hydro_press(rhodz, theta_rhodz, ps, p)
+      IMPORT
+      REAL(rstd),INTENT(IN)  :: rhodz(iim_jjm_i, llm_)
+      REAL(rstd),INTENT(IN)  :: theta_rhodz(iim_jjm_i, llm_, nqdyn_)
+      REAL(rstd),INTENT(OUT) :: ps(iim_jjm_i)
+      REAL(rstd),INTENT(OUT) :: p(iim_jjm_i, llm_)
+    END SUBROUTINE comp_hydro_press
+
   END INTERFACE
 
@@ -36 +44 @@
   PROCEDURE(comp_pression),    POINTER :: compute_pression => NULL(), compute_pression_mid => NULL()
   PROCEDURE(comp_temperature), POINTER :: compute_temperature => NULL()
+  PROCEDURE(comp_hydro_press), POINTER :: compute_hydrostatic_pressure => NULL()
 
-  PUBLIC :: compute_rhodz, compute_pression, compute_pression_mid, compute_temperature
+  PUBLIC :: compute_rhodz, compute_pression, compute_pression_mid, compute_temperature, compute_hydrostatic_pressure
 
 END MODULE compute_diagnostics_mod

codes/icosagcm/devel/src/diagnostics/compute_pression.F90

@@ -3 +3 @@
   USE icosa
   USE omp_para
-  USE disvert_mod, ONLY : ap, bp, ap_bp_present
+  USE disvert_mod, ONLY : ap, bp, ap_bp_present, ptop, caldyn_eta, eta_lag
   IMPLICIT NONE
   PRIVATE
@@ -10 +10 @@
 
   PUBLIC :: pression, compute_pression_hex, compute_pression_unst, &
-       pression_mid, compute_pression_mid_hex, compute_pression_mid_unst
+       pression_mid, compute_pression_mid_hex, compute_pression_mid_unst, &
+       hydrostatic_pressure, compute_hydrostatic_pressure_unst, compute_hydrostatic_pressure_hex
 
 CONTAINS
@@ -92 +93 @@
   END SUBROUTINE pression_mid
 
+  SUBROUTINE hydrostatic_pressure(f_rhodz, f_theta_rhodz, f_ps, f_p)
+    TYPE(t_field), POINTER :: f_rhodz(:), f_theta_rhodz(:), f_ps(:), f_p(:)
+    REAL(rstd), POINTER :: ps(:), rhodz(:,:), p(:,:), theta_rhodz(:,:,:)
+    INTEGER :: ind
+    DO ind=1,ndomain
+      IF (.NOT. assigned_domain(ind)) CYCLE
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+      rhodz=f_rhodz(ind)
+      theta_rhodz=f_theta_rhodz(ind)
+      ps=f_ps(ind)
+      p=f_p(ind)
+      CALL compute_hydrostatic_pressure(rhodz, theta_rhodz, ps, p)
+    ENDDO
+  END SUBROUTINE hydrostatic_pressure
+
 !------------- hexagonal-mesh compute kernels --------
 
@@ -115 +132 @@
 #undef AP
 #undef BP
+
+  SUBROUTINE compute_hydrostatic_pressure_hex(rhodz, theta_rhodz, ps, pk)
+    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm) ! mass per unit surface in each model level
+    REAL(rstd),INTENT(IN)  :: theta_rhodz(iim*jjm,llm, nqdyn) ! dynamical tracers (theta/entropy)
+    REAL(rstd),INTENT(OUT) :: ps(iim*jjm)        ! surface pressure, diagnosed if Lagrangian vertical coordinate 
+    REAL(rstd),INTENT(OUT) :: pk(iim*jjm,llm)    ! pressure at full levels
+    INTEGER :: ij,l, ij_omp_begin_ext, ij_omp_end_ext
+    !$OMP BARRIER
+    CALL distrib_level(ij_begin_ext,ij_end_ext, ij_omp_begin_ext,ij_omp_end_ext)
+#include "../kernels_hex/compute_hydrostatic_pressure.k90"
+    !$OMP BARRIER
+  END SUBROUTINE compute_hydrostatic_pressure_hex
 
 !----------- unstructured-mesh compute kernels --------
@@ -140 +169 @@
 #undef BP
 
+  SUBROUTINE compute_hydrostatic_pressure_unst(rhodz, theta_rhodz, ps, pk)
+    FIELD_MASS,  INTENT(IN)  :: rhodz ! mass per unit surface in each model level
+    FIELD_THETA, INTENT(IN)  :: theta_rhodz(iim*jjm,llm, nqdyn) ! dynamical tracers (theta/entropy)
+    FIELD_PS,    INTENT(OUT) :: ps(iim*jjm)        ! surface pressure, diagnosed if Lagrangian vertical coordinate 
+    FIELD_MASS,  INTENT(OUT) :: pk(iim*jjm,llm)    ! pressure at full levels
+    DECLARE_INDICES
+#include "../kernels_unst/compute_hydrostatic_pressure.k90"
+  END SUBROUTINE compute_hydrostatic_pressure_unst
+  
+
 END MODULE compute_pression_mod

codes/icosagcm/devel/src/diagnostics/observable.f90

@@ -43 +43 @@
     USE xios_mod
     USE earth_const
-    USE compute_pression_mod, ONLY : pression_mid
+    USE compute_pression_mod, ONLY : pression_mid, hydrostatic_pressure
     USE compute_temperature_mod
     USE compute_velocity_mod
@@ -98 +98 @@
     END IF
 
-    CALL pression_mid(f_ps, f_pmid)
+    CALL hydrostatic_pressure(f_mass, f_theta_rhodz, f_ps, f_pmid) 
     CALL temperature(f_pmid, f_q, f_buf_i) ! f_buf_i : IN = theta, out = T
 

codes/icosagcm/devel/src/dynamics/compute_geopot.F90

@@ -17 +17 @@
 
 #ifdef BEGIN_DYSL
-
-#define THECELL {{ thecell }}
 
 {# ---------------- macro to generate code computing pressure top-down ---------------
@@ -59 +57 @@
 
 #define END_GEOPOT {% endcall %}
+
+#define THECELL {{ thecell }}
+
+KERNEL(compute_hydrostatic_pressure)
+  SELECT CASE(caldyn_thermo)
+  CASE(thermo_boussinesq)
+    ! use hydrostatic balance with theta*rhodz to find pk (=Lagrange multiplier=pressure)
+    BALANCE( theta_rhodz(THECELL,1) )
+  CASE(thermo_theta, thermo_entropy, thermo_variable_Cp)
+    BALANCE( rhodz(THECELL) )
+  CASE(thermo_moist)
+    BALANCE( (rhodz(THECELL)+theta_rhodz(THECELL,2)) )
+  END SELECT
+END_BLOCK
 
 KERNEL(compute_geopot)