Changeset 373 for codes/icosagcm

Timestamp:

04/08/16 09:34:17 (8 years ago)

Author:

dubos

Message:

Mass-based NH

Location:

codes/icosagcm/trunk/src

Files:

• caldyn_hevi.f90 (modified) (2 diffs)
• caldyn_kernels_base.f90 (modified) (5 diffs)

codes/icosagcm/trunk/src/caldyn_hevi.f90

@@ -118 +118 @@
           dW = f_dW_fast(ind)
           dPhi = f_dPhi_fast(ind)
-          CALL compute_caldyn_solver(tau,mass,theta,pk,geopot,W,dPhi,dW,du) ! computes d(Phi,W)_fast and updates Phi,W
+          CALL compute_caldyn_solver(tau,mass,theta,pk,geopot,W,dPhi,dW,du) ! computes d(Phi,W,du)_fast and updates Phi,W
        END IF
        u=f_u(ind)
-       CALL compute_caldyn_fast(tau,u,mass,theta,pk,geopot,du) ! computes du_fast and updates du
+       CALL compute_caldyn_fast(tau,u,mass,theta,pk,geopot,du) ! computes du_fast and updates u
     ENDDO
     
@@ -167 +167 @@
           wwuu=f_wwuu(ind)
           dps=f_dps(ind)
-          IF(hydrostatic) THEN
-             CALL compute_caldyn_vert(u,theta,mass,convm, wflux,wwuu, dps, dtheta_rhodz, du)
-          ELSE
-             PRINT *, 'NH dynamics limited to Lagrangian vertical coordinate so far !' 
-             STOP
+          CALL compute_caldyn_vert(u,theta,mass,convm, wflux,wwuu, dps, dtheta_rhodz, du)
+          IF(.NOT.hydrostatic) THEN
+             CALL compute_caldyn_vert_NH(mass,geopot,W,wflux, du,dPhi,dW)
           END IF
        END IF

codes/icosagcm/trunk/src/caldyn_kernels_base.f90

@@ -2 +2 @@
   USE icosa
   USE transfert_mod
+  USE disvert_mod
+  USE omp_para
+  USE trace
   IMPLICIT NONE
   PRIVATE
@@ -20 +23 @@
   TYPE(t_message),PUBLIC :: req_ps, req_mass, req_theta_rhodz, req_u, req_qu, req_geopot, req_w
 
-  PUBLIC :: compute_geopot, compute_caldyn_vert
+  PUBLIC :: compute_geopot, compute_caldyn_vert, compute_caldyn_vert_nh
 
 CONTAINS
@@ -27 +30 @@
 
   SUBROUTINE compute_geopot(ps,rhodz,theta, pk,geopot) 
-    USE icosa
-    USE disvert_mod
-    USE exner_mod
-    USE trace
-    USE omp_para
-    IMPLICIT NONE
     REAL(rstd),INTENT(INOUT) :: ps(iim*jjm)
     REAL(rstd),INTENT(IN)    :: rhodz(iim*jjm,llm)
@@ -134 +131 @@
 
   SUBROUTINE compute_caldyn_vert(u,theta,rhodz,convm, wflux,wwuu, dps,dtheta_rhodz,du)
-    USE icosa
-    USE disvert_mod
-    USE exner_mod
-    USE trace
-    USE omp_para
-    IMPLICIT NONE
     REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)
     REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm)
@@ -264 +255 @@
   END SUBROUTINE compute_caldyn_vert
 
+  SUBROUTINE compute_caldyn_vert_NH(mass,geopot,W,wflux, du,dPhi,dW)
+    REAL(rstd),INTENT(IN) :: mass(iim*jjm,llm)
+    REAL(rstd),INTENT(IN) :: geopot(iim*jjm,llm+1)
+    REAL(rstd),INTENT(IN) :: W(iim*jjm,llm+1)
+    REAL(rstd),INTENT(IN) :: wflux(iim*jjm,llm+1)
+    REAL(rstd),INTENT(INOUT) :: du(iim*3*jjm,llm)
+    REAL(rstd),INTENT(INOUT) :: dPhi(iim*jjm,llm+1)
+    REAL(rstd),INTENT(INOUT) :: dW(iim*jjm,llm+1)
+    ! local arrays
+    REAL(rstd) :: eta_dot(iim*jjm) ! eta_dot in full layers
+    REAL(rstd) :: wcov(iim*jjm) ! covariant vertical momentum in full layers
+    REAL(rstd) :: W_etadot(iim*jjm,llm) ! vertical flux of vertical momentum
+    ! indices and temporary values
+    INTEGER    :: ij, l
+    REAL(rstd) :: wflux_ij, w_ij
+
+    CALL trace_start("compute_caldyn_vert_nh")
+
+    DO l=ll_begin,ll_end
+       ! compute the local arrays
+       !DIR$ SIMD
+       DO ij=ij_begin_ext,ij_end_ext
+          wflux_ij = .5*(wflux(ij,l)+wflux(ij,l+1))
+          w_ij = .5*(W(ij,l)+W(ij,l+1))/mass(ij,l)
+          W_etadot(ij,l) = wflux_ij*w_ij
+          eta_dot(ij) = wflux_ij / mass(ij,l)
+          wcov(ij) = w_ij*(geopot(ij,l+1)-geopot(ij,l))
+       ENDDO
+       ! add NH term to du
+      !DIR$ SIMD
+      DO ij=ij_begin,ij_end
+          du(ij+u_right,l) = du(ij+u_right,l) &
+               - .5*(wcov(ij+t_right)+wcov(ij)) &
+               *ne_right*(eta_dot(ij+t_right)-eta_dot(ij))
+          du(ij+u_lup,l) = du(ij+u_lup,l) &
+               - .5*(wcov(ij+t_lup)+wcov(ij)) &
+               *ne_lup*(eta_dot(ij+t_lup)-eta_dot(ij))
+          du(ij+u_ldown,l) = du(ij+u_ldown,l) &
+               - .5*(wcov(ij+t_ldown)+wcov(ij)) &
+               *ne_ldown*(eta_dot(ij+t_ldown)-eta_dot(ij))
+       END DO
+    ENDDO
+    ! add NH terms to dW, dPhi
+    ! FIXME : TODO top and bottom
+    DO l=ll_beginp1,ll_end ! inner interfaces only
+       !DIR$ SIMD
+       DO ij=ij_begin,ij_end
+          dW(ij,l) = dW(ij,l) + W_etadot(ij,l-1) - W_etadot(ij,l)
+          dPhi(ij,l) = dPhi(ij,l) - wflux(ij,l) &
+               * (geopot(ij,l+1)-geopot(ij,l-1))/(mass(ij,l-1)+mass(ij,l))
+       END DO
+    END DO
+    CALL trace_end("compute_caldyn_vert_nh")
+
+  END SUBROUTINE compute_caldyn_vert_NH
 END MODULE caldyn_kernels_base_mod