Changeset 181

Timestamp:

12/15/13 01:54:31 (11 years ago)

Author:

dubos

Message:

fixed shallow-water hydrostatic balance

File:

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

codes/icosagcm/trunk/src/caldyn_gcm.f90

@@ -426 +426 @@
 
     IF(caldyn_eta==eta_mass) THEN
-    
+
 !!! Compute exner function and geopotential
        DO l = 1,llm
           !$OMP DO SCHEDULE(STATIC) 
-!DIR$ SIMD
-         DO ij=ij_begin_ext,ij_end_ext         
-            p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij) ! FIXME : leave ps for the moment ; change ps to Ms later
-            !         p_ik = ptop + g*(mass_ak(l)+ mass_bk(l)*ps(i,j))
-            exner_ik = cpp * (p_ik/preff) ** kappa
-            pk(ij,l) = exner_ik
-            ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz
-            geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik
+          !DIR$ SIMD
+          DO ij=ij_begin_ext,ij_end_ext         
+             p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij) ! FIXME : leave ps for the moment ; change ps to Ms later
+             !         p_ik = ptop + g*(mass_ak(l)+ mass_bk(l)*ps(i,j))
+             exner_ik = cpp * (p_ik/preff) ** kappa
+             pk(ij,l) = exner_ik
+             ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz
+             geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik
           ENDDO
        ENDDO
@@ -447 +447 @@
        ! Notice that the computation below should work also when caldyn_eta=eta_mass          
 
-       ! uppermost layer
-!DIR$ SIMD
-       DO ij=ij_begin_ext,ij_end_ext         
-         pk(ij,llm) = ptop + (.5*g)*rhodz(ij,llm)
-       END DO
-       ! other layers
-       DO l = llm-1, 1, -1
-          !$OMP DO SCHEDULE(STATIC) 
-!DIR$ SIMD
+       IF(boussinesq) THEN 
+          ! use theta*rhodz in hydrostatic balance
+          ! uppermost layer
+          !DIR$ SIMD
           DO ij=ij_begin_ext,ij_end_ext         
-             pk(ij,l) = pk(ij,l+1) + (.5*g)*(rhodz(ij,l)+rhodz(ij,l+1))
+             pk(ij,llm) = ptop + (.5*g)*theta(ij,llm)*rhodz(ij,llm)
           END DO
-       END DO
-       ! surface pressure (for diagnostics)
-       DO ij=ij_begin_ext,ij_end_ext         
-         ps(ij) = pk(ij,1) + (.5*g)*rhodz(ij,1)
-       END DO
-
-       IF(boussinesq) THEN ! specific volume 1 = dphi/g/rhodz
+          ! 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)
+
+          ! specific volume 1 = dphi/g/rhodz
           DO l = 1,llm
              !$OMP DO SCHEDULE(STATIC) 
-!DIR$ SIMD
+             !DIR$ SIMD
              DO ij=ij_begin_ext,ij_end_ext         
-               geopot(ij,l+1) = geopot(ij,l) + g*rhodz(ij,l)
+                geopot(ij,l+1) = geopot(ij,l) + g*rhodz(ij,l)
              ENDDO
           ENDDO
-       ELSE ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz
+       ELSE ! non-Boussinesq
+          ! uppermost layer
+          !DIR$ SIMD
+          DO ij=ij_begin_ext,ij_end_ext         
+             pk(ij,llm) = ptop + (.5*g)*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)*(rhodz(ij,l)+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)*rhodz(ij,1)
+          END DO
+
+
+          ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz
           DO l = 1,llm
              !$OMP DO SCHEDULE(STATIC) 
-!DIR$ SIMD
-              DO ij=ij_begin_ext,ij_end_ext         
+             !DIR$ SIMD
+             DO ij=ij_begin_ext,ij_end_ext         
                 p_ik = pk(ij,l)
                 exner_ik = cpp * (p_ik/preff) ** kappa
@@ -666 +689 @@
 !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 +    &
@@ -674 +697 @@
                                      le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    &
                                      le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 )  
+                ! from now on pk contains the vertically-averaged geopotential
                 pk(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
           ENDDO