Changeset 366 for codes/icosagcm/trunk/src/caldyn_kernels_hevi.f90

Timestamp:

10/30/15 15:41:06 (9 years ago)

Author:

dubos

Message:

Progress towards NH

File:

• codes/icosagcm/trunk/src/caldyn_kernels_hevi.f90 (modified) (4 diffs)

codes/icosagcm/trunk/src/caldyn_kernels_hevi.f90

@@ -6 +6 @@
   PRIVATE
 
-  PUBLIC :: compute_theta, compute_pvort_only, compute_caldyn_fast, compute_caldyn_slow
+  PUBLIC :: compute_theta, compute_pvort_only, compute_caldyn_slow, &
+       compute_caldyn_solver, compute_caldyn_fast
 
 CONTAINS
@@ -120 +121 @@
   END SUBROUTINE compute_pvort_only
 
-  SUBROUTINE compute_caldyn_fast(tau,u,rhodz,theta,pk,geopot, du)
+  SUBROUTINE compute_caldyn_solver(tau,rhodz,theta,pk, geopot,W, dPhi,dW)
+    USE icosa
+    USE disvert_mod
+    USE exner_mod
+    USE trace
+    USE omp_para
+    USE disvert_mod
+    IMPLICIT NONE
+    REAL(rstd),INTENT(IN) :: tau ! "solve" Phi-tau*dPhi/dt = Phi_rhs
+    REAL(rstd),INTENT(IN)    :: rhodz(iim*jjm,llm)
+    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm)
+    REAL(rstd),INTENT(OUT)   :: pk(iim*jjm,llm)
+    REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1)
+    REAL(rstd),INTENT(INOUT) :: W(iim*jjm,llm+1) ! OUT if tau>0
+    REAL(rstd),INTENT(OUT)   :: dW(iim*jjm,llm+1)
+    REAL(rstd),INTENT(OUT)   :: dPhi(iim*jjm,llm+1)
+
+    INTEGER :: ij,l
+    INTEGER    :: ij_omp_begin_ext, ij_omp_end_ext
+    REAL(rstd) :: gamma, rho_ij, X_ij, Y_ij, m_il
+
+    CALL trace_start("compute_caldyn_solver")
+
+    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
+
+    IF(tau>0) THEN
+       PRINT *, 'No implicit solver for NH yet'
+       STOP
+    END IF
+
+    ! Compute pressure, stored temporarily in pk
+    ! kappa = R/Cp
+    ! 1-kappa = Cv/Cp
+    ! Cp/Cv = 1/(1-kappa)
+    gamma = 1./(1.-kappa)
+    DO l=ll_begin,ll_end
+       !DIR$ SIMD
+       DO ij=ij_omp_begin_ext,ij_omp_end_ext
+          rho_ij = (g*rhodz(ij,l))/(geopot(ij,l+1)-geopot(ij,l))
+          X_ij = (cpp/preff)*kappa*theta(ij,l)*rho_ij
+          ! kappa.theta.rho = p/exner 
+          ! => X = (p/p0)/(exner/Cp)
+          !      = (p/p0)^(1-kappa)
+          pk(ij,l) = preff*(X_ij**gamma)
+       ENDDO
+    ENDDO
+
+    ! Compute tendencies for geopotential and vertical momentum
+    DO l=ll_beginp1,ll_end
+       !DIR$ SIMD
+       DO ij=ij_omp_begin_ext,ij_omp_end_ext
+          m_il = .5*(rhodz(ij,l)+rhodz(ij,l-1))
+          dW(ij,l) = (1./g)*(pk(ij,l-1)-pk(ij,l)) - m_il
+          dPhi(ij,l) = g*g*W(ij,l)/m_il
+       ENDDO
+       !       PRINT *,'Max dPhi', l,ij_begin,ij_end, MAXVAL(abs(dPhi(ij_begin:ij_end,l)))
+       !       PRINT *,'Max dW', l,ij_begin,ij_end, MAXVAL(abs(dW(ij_begin:ij_end,l))) 
+    ENDDO
+    ! Lower BC (FIXME : no orography yet !)
+    DO ij=ij_begin,ij_end         
+       dPhi(ij,1)=0
+       W(ij,1)=0
+       dW(ij,1)=0
+       dPhi(ij,llm+1)=0 ! rigid lid
+       W(ij,llm+1)=0
+       dW(ij,llm+1)=0
+    ENDDO
+    ! Upper BC p=ptop
+!    DO ij=ij_omp_begin_ext,ij_omp_end_ext
+!       dPhi(ij,llm+1) = W(ij,llm+1)/rhodz(ij,llm)
+!       dW(ij,llm+1) = (1./g)*(pk(ij,llm)-ptop) - .5*rhodz(ij,llm)
+!    ENDDO
+    
+    ! Compute Exner function (needed by compute_caldyn_fast)
+    DO l=ll_begin,ll_end
+       !DIR$ SIMD
+       DO ij=ij_omp_begin_ext,ij_omp_end_ext
+          pk(ij,l) = cpp*((pk(ij,l)/preff)**kappa) ! other formulae possible if exponentiation is slow
+       ENDDO
+    ENDDO
+    
+    CALL trace_end("compute_caldyn_solver")
+    
+  END SUBROUTINE compute_caldyn_solver
+  
+  SUBROUTINE compute_caldyn_fast(tau,u,rhodz,theta,pk,geopot,du)
     USE icosa
     USE disvert_mod
@@ -127 +215 @@
     USE omp_para
     IMPLICIT NONE
-    REAL(rstd), INTENT(IN) :: tau ! "solve" u-tau*du/dt = rhs
-    REAL(rstd),INTENT(INOUT)  :: u(iim*3*jjm,llm)    ! prognostic "velocity"
-    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)
-    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm)  ! potential temperature
-    REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm) ! Exner function
-    REAL(rstd),INTENT(IN)  :: geopot(iim*jjm,llm+1)    ! geopotential
-    REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm)
+    REAL(rstd),INTENT(IN)    :: tau                ! "solve" u-tau*du/dt = rhs
+    REAL(rstd),INTENT(INOUT) :: u(iim*3*jjm,llm)   ! OUT if tau>0
+    REAL(rstd),INTENT(IN)    :: rhodz(iim*jjm,llm)
+    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm)
+    REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm)
+    REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1)
+    REAL(rstd),INTENT(OUT)   :: du(iim*3*jjm,llm)
     REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function
 
@@ -140 +228 @@
 
     CALL trace_start("compute_caldyn_fast")
-    
-    ! Compute bernouilli term
+
+    ! Compute Bernoulli term
     IF(boussinesq) THEN
+
        DO l=ll_begin,ll_end
           !DIR$ SIMD