Changeset 599 for codes/icosagcm/trunk/src/transport/advect.F90

Timestamp:

10/19/17 17:04:26 (7 years ago)

Author:

dubos

Message:

trunk : backported commits r582-r598 (transport diagnostics)

File:

• codes/icosagcm/trunk/src/transport/advect.F90 (moved) (moved from codes/icosagcm/trunk/src/transport/advect.f90) (3 diffs)

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

@@ -425 +425 @@
   ! Horizontal transport (S. Dubey, T. Dubos)
   ! Slope-limited van Leer approach with hexagons
-  SUBROUTINE compute_advect_horiz(update_mass,hfluxt,cc,gradq3d, mass,qi) 
+  SUBROUTINE compute_advect_horiz(update_mass,diagflux_on, hfluxt,cc,gradq3d, mass,qi,qfluxt)
     USE trace
     USE omp_para
     IMPLICIT NONE
-    LOGICAL, INTENT(IN)       :: update_mass
+    LOGICAL, INTENT(IN)       :: update_mass, diagflux_on
     REAL(rstd), INTENT(IN)    :: gradq3d(iim*jjm,llm,3) 
     REAL(rstd), INTENT(IN)    :: hfluxt(3*iim*jjm,llm) ! mass flux
@@ -435 +435 @@
     REAL(rstd), INTENT(INOUT) :: mass(iim*jjm,llm)
     REAL(rstd), INTENT(INOUT) :: qi(iim*jjm,llm)
+    REAL(rstd), INTENT(INOUT) :: qfluxt(3*iim*jjm,llm) ! time-integrated tracer flux
 
     REAL(rstd) :: dq,dmass,qe,ed(3), newmass
@@ -441 +442 @@
 
     CALL trace_start("compute_advect_horiz")
-
-    ! evaluate tracer field at point cc using piecewise linear reconstruction
-    ! q(cc)= q0 + gradq.(cc-xyz_i)  with xi centroid of hexagon
-    ! ne*hfluxt>0 iff outgoing
-    DO l = ll_begin,ll_end
-
-!DIR$ SIMD
-      DO ij=ij_begin_ext,ij_end_ext
-
-             IF (ne(ij,right)*hfluxt(ij+u_right,l)>0) THEN 
-                ed = cc(ij+u_right,l,:) - xyz_i(ij,:)
-!                qe = qi(ij,l)+sum2(gradq3d(ij,l,:),ed) 
-                 qe = qi(ij,l)+gradq3d(ij,l,1)*ed(1)+gradq3d(ij,l,2)*ed(2)+gradq3d(ij,l,3)*ed(3) 
-             ELSE
-                ed = cc(ij+u_right,l,:) - xyz_i(ij+t_right,:)
-!                qe = qi(ij+t_right,l)+sum2(gradq3d(ij+t_right,l,:),ed) 
-                qe = qi(ij+t_right,l) + gradq3d(ij+t_right,l,1)*ed(1)+gradq3d(ij+t_right,l,2)*ed(2)+gradq3d(ij+t_right,l,3)*ed(3) 
-             ENDIF
-             qflux(ij+u_right,l) = hfluxt(ij+u_right,l)*qe
-             
-             IF (ne(ij,lup)*hfluxt(ij+u_lup,l)>0) THEN 
-                ed = cc(ij+u_lup,l,:) - xyz_i(ij,:)
-!                qe = qi(ij,l)+sum2(gradq3d(ij,l,:),ed)
-                qe = qi(ij,l)  + gradq3d(ij,l,1)*ed(1)+gradq3d(ij,l,2)*ed(2)+gradq3d(ij,l,3)*ed(3) 
-             ELSE
-                ed = cc(ij+u_lup,l,:) - xyz_i(ij+t_lup,:)
-!                qe = qi(ij+t_lup,l)+sum2(gradq3d(ij+t_lup,l,:),ed) 
-                qe = qi(ij+t_lup,l) + gradq3d(ij+t_lup,l,1)*ed(1)+gradq3d(ij+t_lup,l,2)*ed(2)+gradq3d(ij+t_lup,l,3)*ed(3) 
-             ENDIF
-             qflux(ij+u_lup,l) = hfluxt(ij+u_lup,l)*qe 
-
-             IF (ne(ij,ldown)*hfluxt(ij+u_ldown,l)>0) THEN 
-                ed = cc(ij+u_ldown,l,:) - xyz_i(ij,:)
-!                qe = qi(ij,l)+sum2(gradq3d(ij,l,:),ed) 
-                qe = qi(ij,l) + gradq3d(ij,l,1)*ed(1)+gradq3d(ij,l,2)*ed(2)+gradq3d(ij,l,3)*ed(3) 
-             ELSE
-                ed = cc(ij+u_ldown,l,:) - xyz_i(ij+t_ldown,:)
-!                qe = qi(ij+t_ldown,l)+sum2(gradq3d(ij+t_ldown,l,:),ed) 
-                qe = qi(ij+t_ldown,l)+gradq3d(ij+t_ldown,l,1)*ed(1)+gradq3d(ij+t_ldown,l,2)*ed(2)+gradq3d(ij+t_ldown,l,3)*ed(3) 
-             ENDIF
-             qflux(ij+u_ldown,l) = hfluxt(ij+u_ldown,l)*qe
-       END DO
-    END DO
-
-    ! update q and, if update_mass, update mass
-    DO l = ll_begin,ll_end
-!DIR$ SIMD
-      DO ij=ij_begin,ij_end
-             ! sign convention as Ringler et al. (2010) eq. 21
-             dmass =   hfluxt(ij+u_right,l)*ne(ij,right)   &
-                     + hfluxt(ij+u_lup,l)  *ne(ij,lup)     &
-                     + hfluxt(ij+u_ldown,l)*ne(ij,ldown)   &
-                     + hfluxt(ij+u_rup,l)  *ne(ij,rup)     &
-                     + hfluxt(ij+u_left,l) *ne(ij,left)    &
-                     + hfluxt(ij+u_rdown,l)*ne(ij,rdown)
-
-             dq    =   qflux(ij+u_right,l) *ne(ij,right)   &
-                     + qflux(ij+u_lup,l)   *ne(ij,lup)     &
-                     + qflux(ij+u_ldown,l) *ne(ij,ldown)   &
-                     + qflux(ij+u_rup,l)   *ne(ij,rup)     &
-                     + qflux(ij+u_left,l)  *ne(ij,left)    &
-                     + qflux(ij+u_rdown,l) *ne(ij,rdown)
-
-             
-             newmass = mass(ij,l) - dmass/Ai(ij)
-             qi(ij,l) = (qi(ij,l)*mass(ij,l) - dq/Ai(ij) ) / newmass
-             IF(update_mass) mass(ij,l) = newmass
-
-       END DO
-    END DO
-
+#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