Changeset 594

Timestamp:

10/18/17 18:46:11 (7 years ago)

Author:

dubos

Message:

devel : diagnose energy fluxes

Location:

codes/icosagcm/devel/src

Files:

• diagnostics/diagflux.F90 (modified) (7 diffs)
• diagnostics/observable.f90 (modified) (1 diff)
• kernels/energy_fluxes.k90 (added)
• time/timeloop_gcm.f90 (modified) (2 diffs)
• transport/advect.F90 (modified) (1 diff)

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

@@ -1 +1 @@
 MODULE diagflux_mod
   USE icosa
+  USE omp_para
   IMPLICIT NONE
   SAVE
-  
-  TYPE(t_field),POINTER :: &
+  PRIVATE
+
+  TYPE(t_field), POINTER, PUBLIC :: &
        f_masst(:), f_qmasst(:), & ! time-averaged mass, tracer mass, 
        f_massfluxt(:), f_qfluxt(:), & ! time-integrated mass flux and tracer flux
@@ -13 +15 @@
   !$OMP THREADPRIVATE(diagflux_on)
 
+  PUBLIC :: diagflux_on, init_diagflux, zero_qfluxt, qflux_centered_lonlat, flux_centered_lonlat, diagflux_energy
+
 CONTAINS
 
@@ -20 +24 @@
     diagflux_on = .FALSE.
     CALL getin("diagflux", diagflux_on)
-    ll = MERGE(ll,0,diagflux_on)
+    ll = MERGE(llm,1,diagflux_on)
     CALL allocate_field(f_masst,         field_t,type_real,ll,       name="masst")
     CALL allocate_field(f_epot,          field_t,type_real,ll,       name="epot")
@@ -36 +40 @@
   END SUBROUTINE init_diagflux
 
+#define ZERO2(field) buf2=field(ind) ; buf2(:,ll_begin:ll_end)=0.
+#define ZERO3(field) buf3=field(ind) ; buf3(:,ll_begin:ll_end,:)=0.
+
   SUBROUTINE zero_qfluxt
-    USE omp_para
     INTEGER :: ind
     REAL(rstd), POINTER :: buf2(:,:),buf3(:,:,:)
@@ -43 +49 @@
        IF (.NOT. assigned_domain(ind)) CYCLE
        CALL swap_dimensions(ind)
-       buf2=f_masst(ind)
-       buf2(:,ll_begin:ll_end)=0.
-       buf2=f_massfluxt(ind)
-       buf2(:,ll_begin:ll_end)=0.
-       buf3=f_qmasst(ind)
-       buf3(:,ll_begin:ll_end,:)=0.
-       buf3=f_qfluxt(ind)
-       buf3(:,ll_begin:ll_end,:)=0.
+       ZERO2(f_masst)
+       ZERO2(f_epot)
+       ZERO2(f_ekin)
+       ZERO2(f_enthalpy)
+       ZERO3(f_qmasst)
+       ZERO2(f_massfluxt)
+       ZERO2(f_epotfluxt)
+       ZERO2(f_ekinfluxt)
+       ZERO2(f_enthalpyfluxt)
+       ZERO3(f_qfluxt)
     END DO
   END SUBROUTINE zero_qfluxt
 
-  SUBROUTINE flux_centered_lonlat(scale, f_massflux, f_flux, f_massflux_lon, f_massflux_lat, f_flux_lon, f_flux_lat)
+!------------------------------------ Reconstruct fluxes at cell centers ---------------------------------------
+
+  SUBROUTINE qflux_centered_lonlat(scale, f_flux, f_flux_lon, f_flux_lat)
     REAL(rstd), INTENT(IN) :: scale
-    TYPE(t_field),POINTER :: f_flux(:), f_flux_lon(:), f_flux_lat(:), &
-         f_massflux(:), f_massflux_lon(:), f_massflux_lat(:)
-    REAL(rstd), POINTER :: flux(:,:,:), flux_lon(:,:,:), flux_lat(:,:,:), &
-         massflux(:,:), massflux_lon(:,:), massflux_lat(:,:)
+    TYPE(t_field),POINTER :: f_flux(:), f_flux_lon(:), f_flux_lat(:)
+    REAL(rstd), POINTER :: flux(:,:,:), flux_lon(:,:,:), flux_lat(:,:,:)
     INTEGER :: ind, itrac
     DO ind=1,ndomain
@@ -71 +79 @@
           CALL compute_flux_centered_lonlat(scale, flux(:,:,itrac), flux_lon(:,:,itrac), flux_lat(:,:,itrac))
        END DO
-       massflux=f_massflux(ind)
-       massflux_lon=f_massflux_lon(ind)
-       massflux_lat=f_massflux_lat(ind)
-       CALL compute_flux_centered_lonlat(scale, massflux, massflux_lon, massflux_lat)
+    END DO
+  END SUBROUTINE qflux_centered_lonlat
+  
+  SUBROUTINE flux_centered_lonlat(scale, f_flux, f_flux_lon, f_flux_lat)
+    REAL(rstd), INTENT(IN) :: scale
+    TYPE(t_field),POINTER :: f_flux(:), f_flux_lon(:), f_flux_lat(:)
+    REAL(rstd), POINTER :: flux(:,:), flux_lon(:,:), flux_lat(:,:)
+    INTEGER :: ind
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       flux=f_flux(ind)
+       flux_lon=f_flux_lon(ind)
+       flux_lat=f_flux_lat(ind)
+       CALL compute_flux_centered_lonlat(scale, flux, flux_lon, flux_lat)
     END DO
   END SUBROUTINE flux_centered_lonlat
@@ -88 +108 @@
   END SUBROUTINE compute_flux_centered_lonlat
 
+!------------------------------------ Compute energy fluxes ---------------------------------------
+
+  SUBROUTINE diagflux_energy(frac, f_phis,f_rhodz,f_theta_rhodz,f_u, f_geopot,f_theta, f_hfluxt)
+    REAL(rstd), INTENT(IN) :: frac
+    TYPE(t_field),POINTER :: f_phis(:),f_rhodz(:),f_theta_rhodz(:),f_u(:), f_geopot(:), f_theta(:), f_hfluxt(:)
+    REAL(rstd), POINTER :: phis(:), rhodz(:,:), theta_rhodz(:,:,:), u(:,:), &
+         geopot(:,:), pk(:,:,:), hfluxt(:,:), &
+         epot(:,:), ekin(:,:), enthalpy(:,:), &
+         epotflux(:,:), ekinflux(:,:), enthalpyflux(:,:)
+    INTEGER :: ind
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       hfluxt = f_hfluxt(ind)
+       phis = f_phis(ind)
+       rhodz = f_rhodz(ind)
+       theta_rhodz = f_theta_rhodz(ind)
+       u = f_u(ind)
+       geopot = f_geopot(ind)
+       pk = f_theta(ind) ! buffer
+       epot = f_epot(ind)
+       ekin = f_ekin(ind)
+       enthalpy = f_enthalpy(ind)
+       epotflux = f_epotfluxt(ind)
+       ekinflux = f_ekinfluxt(ind)
+       enthalpyflux = f_enthalpyfluxt(ind)
+       CALL compute_diagflux_energy(frac,hfluxt, phis,rhodz,theta_rhodz,u, geopot,pk, epot,ekin,enthalpy, epotflux, ekinflux, enthalpyflux)
+    END DO
+  END SUBROUTINE diagflux_energy
+
+  SUBROUTINE compute_diagflux_energy(frac, massflux, phis,rhodz,theta_rhodz,u, geopot,pk, epot,ekin,enthalpy, epot_flux, ekin_flux, enthalpy_flux)
+    USE disvert_mod, ONLY : ptop
+    REAL(rstd), INTENT(IN) :: frac
+    REAL(rstd), INTENT(IN) :: massflux(3*iim*jjm,llm), u(3*iim*jjm,llm),&
+                              phis(iim*jjm), rhodz(iim*jjm,llm), theta_rhodz(iim*jjm,llm,nqtot)
+    REAL(rstd), INTENT(INOUT) :: geopot(iim*jjm,llm+1), pk(iim*jjm,llm) ! pk = buffer
+    REAL(rstd), INTENT(INOUT), DIMENSION(iim*jjm, llm)   ::  epot, ekin, enthalpy
+    REAL(rstd), INTENT(INOUT), DIMENSION(3*iim*jjm, llm) ::  epot_flux, ekin_flux, enthalpy_flux    
+    REAL(rstd) :: energy, p_ik, theta_ik, temp_ik, gv, Rd 
+    INTEGER :: ij, l, ij_omp_begin_ext, ij_omp_end_ext
+    Rd = kappa*cpp
+    CALL distrib_level(ij_end_ext-ij_begin_ext+1,ij_omp_begin_ext,ij_omp_end_ext)
+    ij_omp_begin_ext = ij_omp_begin_ext+ij_begin_ext-1
+    ij_omp_end_ext = ij_omp_end_ext+ij_begin_ext-1
+#include "../kernels/energy_fluxes.k90"
+  END SUBROUTINE compute_diagflux_energy
+
 END MODULE diagflux_mod

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

@@ -163 +163 @@
     IF(.NOT. init) THEN
        IF(diagflux_on) THEN
+          CALL flux_centered_lonlat(1./(itau_out*dt) , f_massfluxt, f_buf_ulon, f_buf_ulat)
           CALL output_field("mass_t", f_masst)
-          CALL output_field("qmass_t", f_qmasst)
-          CALL flux_centered_lonlat(1./(itau_out*dt) , f_massfluxt, f_qfluxt, &
-               f_buf_ulon, f_buf_ulat, f_qfluxt_lon, f_qfluxt_lat)
           CALL output_field("massflux_lon",f_buf_ulon)
           CALL output_field("massflux_lat",f_buf_ulat)
+          CALL flux_centered_lonlat(1./(itau_out*dt) , f_epotfluxt, f_buf_ulon, f_buf_ulat)
+          CALL output_field("epot_t", f_epot)
+          CALL output_field("epotflux_lon",f_buf_ulon)
+          CALL output_field("epotflux_lat",f_buf_ulat)
+          CALL flux_centered_lonlat(1./(itau_out*dt) , f_ekinfluxt, f_buf_ulon, f_buf_ulat)
+          CALL output_field("ekin_t", f_ekin)
+          CALL output_field("ekinflux_lon",f_buf_ulon)
+          CALL output_field("ekinflux_lat",f_buf_ulat)
+          CALL flux_centered_lonlat(1./(itau_out*dt) , f_enthalpyfluxt, f_buf_ulon, f_buf_ulat)
+          CALL output_field("enthalpy_t", f_enthalpy)
+          CALL output_field("enthalpyflux_lon",f_buf_ulon)
+          CALL output_field("enthalpyflux_lat",f_buf_ulat)
+          CALL qflux_centered_lonlat(1./(itau_out*dt) , f_qfluxt, f_qfluxt_lon, f_qfluxt_lat)
+          CALL output_field("qmass_t", f_qmasst)
           CALL output_field("qflux_lon",f_qfluxt_lon)
           CALL output_field("qflux_lat",f_qfluxt_lat)

codes/icosagcm/devel/src/time/timeloop_gcm.f90

@@ -273 +273 @@
           CALL HEVI_scheme(it, fluxt_zero)
        END SELECT
-       
+
        IF (MOD(it,itau_dissip)==0) THEN
           
@@ -304 +304 @@
        
        IF(MOD(it,itau_adv)==0) THEN
+          IF(diagflux_on) CALL diagflux_energy(adv_over_out, f_phis,f_rhodz,f_theta_rhodz,f_u, f_geopot,f_theta, f_hfluxt)
           CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz, & ! update q and rhodz after RK step
                adv_over_out, f_masst,f_qmasst,f_massfluxt, f_qfluxt)  ! accumulate mass and fluxes if diagflux_on
-          fluxt_zero=.TRUE.
-          ! FIXME : check that rhodz is consistent with ps
-          IF (check_rhodz) THEN
+          fluxt_zero=.TRUE. ! restart accumulation of hfluxt and qfluxt at next call to explicit_scheme / HEVI_scheme
+          IF (check_rhodz) THEN ! check that rhodz is consistent with ps
              DO ind=1,ndomain
                 IF (.NOT. assigned_domain(ind)) CYCLE

codes/icosagcm/devel/src/transport/advect.F90

@@ -444 +444 @@
 #include "../kernels/advect_horiz.k90"
     CALL trace_end("compute_advect_horiz")
-  CONTAINS
-
-    !====================================================================================
-    PURE REAL(rstd) FUNCTION sum2(a1,a2)
-      IMPLICIT NONE
-      REAL(rstd),INTENT(IN):: a1(3), a2(3)
-      sum2 = a1(1)*a2(1)+a1(2)*a2(2)+a1(3)*a2(3)
-      ! sum2 = 0. ! Godunov scheme
-    END FUNCTION sum2
-
   END SUBROUTINE compute_advect_horiz