Changeset 183

Timestamp:

12/16/13 00:48:24 (11 years ago)

Author:

dubos

Message:

moved Boussniesq calculation of pressure into compute_caldyn_horiz

File:

• codes/icosagcm/trunk/src/caldyn_gcm.f90 (modified) (11 diffs)

codes/icosagcm/trunk/src/caldyn_gcm.f90

@@ -221 +221 @@
           pk = f_pk(ind)
           geopot = f_geopot(ind)  
-          CALL compute_geopot(ps,mass,theta,u, pk,geopot)
+          CALL compute_geopot(ps,mass,theta, pk,geopot)
           hflux=f_hflux(ind)
           convm = f_dmass(ind)
@@ -249 +249 @@
           pk = f_pk(ind)
           geopot = f_geopot(ind)  
-          CALL compute_geopot(ps,mass,theta,u, pk,geopot)
+          CALL compute_geopot(ps,mass,theta, pk,geopot)
           hflux=f_hflux(ind)
           convm = f_dmass(ind)
@@ -405 +405 @@
   END SUBROUTINE compute_pvort
   
-  SUBROUTINE compute_geopot(ps,rhodz,theta,u, pk,geopot) 
+  SUBROUTINE compute_geopot(ps,rhodz,theta, pk,geopot) 
   USE icosa
   USE disvert_mod
@@ -415 +415 @@
     REAL(rstd),INTENT(IN)    :: rhodz(iim*jjm,llm)
     REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm)    ! potential temperature
-    REAL(rstd),INTENT(IN)    :: u(3*iim*jjm,llm)      ! prognostic velocity
     REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm)       ! Exner function
     REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1) ! geopotential
     
-    REAL(rstd) :: urel(3*iim*jjm,llm)  ! relative velocity
-    REAL(rstd) :: cor_NT(iim*jjm,llm)  ! NT coriolis force u.(du/dPhi)
     INTEGER :: i,j,ij,l
     REAL(rstd) :: p_ik, exner_ik
@@ -448 +445 @@
        ! Notice that the computation below should work also when caldyn_eta=eta_mass          
 
-       IF(boussinesq) THEN 
-          ! use theta*rhodz in hydrostatic balance
-          ! uppermost layer
-          !DIR$ SIMD
-          DO ij=ij_begin_ext,ij_end_ext         
-             pk(ij,llm) = ptop + (.5*g)*theta(ij,llm)*rhodz(ij,llm)
-          END DO
-          ! other layers
-          DO l = llm-1, 1, -1
-             !$OMP DO SCHEDULE(STATIC) 
-             !DIR$ SIMD
-             DO ij=ij_begin_ext,ij_end_ext         
-                pk(ij,l) = pk(ij,l+1) + (.5*g)*(theta(ij,l)*rhodz(ij,l)+theta(ij,l+1)*rhodz(ij,l+1))
-             END DO
-          END DO
-          ! surface pressure (for diagnostics)
-          DO ij=ij_begin_ext,ij_end_ext         
-             ps(ij) = pk(ij,1) + (.5*g)*theta(ij,1)*rhodz(ij,1)
-          END DO
-          ! now pk contains the Lagrange multiplier (pressure)
-
+       IF(boussinesq) THEN ! compute only geopotential : pressure pk will be computed in compute_caldyn_horiz
           ! specific volume 1 = dphi/g/rhodz
           DO l = 1,llm
@@ -477 +454 @@
              ENDDO
           ENDDO
-       ELSE ! non-Boussinesq
+       ELSE ! non-Boussinesq, compute geopotential and Exner pressure
           ! uppermost layer
           !DIR$ SIMD
@@ -496 +473 @@
           END DO
 
-
           ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz
           DO l = 1,llm
@@ -523 +499 @@
   USE omp_para
   IMPLICIT NONE
-    REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)
+    REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)    ! prognostic "velocity"
     REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)
     REAL(rstd),INTENT(IN)  :: qu(iim*3*jjm,llm)
@@ -535 +511 @@
     REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm)
 
+    REAL(rstd) :: cor_NT(iim*jjm,llm)  ! NT coriolis force u.(du/dPhi)
     REAL(rstd) :: urel(3*iim*jjm,llm)  ! relative velocity
     REAL(rstd) :: Ftheta(3*iim*jjm,llm) ! theta flux
@@ -686 +663 @@
 !!! Compute bernouilli term = Kinetic Energy + geopotential
     IF(boussinesq) THEN
+       ! first use hydrostatic balance with theta*rhodz to find pk (Lagrange multiplier=pressure) 
+       ! uppermost layer
+       !DIR$ SIMD
+       DO ij=ij_begin_ext,ij_end_ext         
+          pk(ij,llm) = ptop + (.5*g)*theta(ij,llm)*rhodz(ij,llm)
+       END DO
+       ! other layers
+       DO l = llm-1, 1, -1
+          !$OMP DO SCHEDULE(STATIC) 
+          !DIR$ SIMD
+          DO ij=ij_begin_ext,ij_end_ext         
+             pk(ij,l) = pk(ij,l+1) + (.5*g)*(theta(ij,l)*rhodz(ij,l)+theta(ij,l+1)*rhodz(ij,l+1))
+          END DO
+       END DO
+       ! surface pressure (for diagnostics) FIXME
+       ! DO ij=ij_begin_ext,ij_end_ext         
+       !   ps(ij) = pk(ij,1) + (.5*g)*theta(ij,1)*rhodz(ij,1)
+       ! END DO
+       ! now pk contains the Lagrange multiplier (pressure)
 
        DO l=ll_begin,ll_end
 !DIR$ SIMD
           DO ij=ij_begin,ij_end         
-                ! until this point pk contains the Lagrange multiplier (pressure)
                 berni(ij,l) = pk(ij,l) &
                      + 1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +    &
@@ -703 +698 @@
        ENDDO
 
-    ELSE
+    ELSE ! compressible
 
        DO l=ll_begin,ll_end