Changeset 167 for codes/icosagcm/trunk/src/caldyn_gcm.f90

Timestamp:

06/29/13 02:35:03 (11 years ago)

Author:

dubos

Message:

Multi-layer Saint-Venant / Ripa equations - tested with Williamson91.5 & JW06

File:

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

codes/icosagcm/trunk/src/caldyn_gcm.f90

@@ -22 +22 @@
 
   PUBLIC init_caldyn, caldyn_BC, caldyn, write_output_fields, &
-       req_ps, req_mass, caldyn_eta, eta_mass, eta_lag
+       req_ps, req_mass
 
 CONTAINS
@@ -122 +122 @@
        f_hflux, f_wflux, f_dps, f_dmass, f_dtheta_rhodz, f_du)
     USE icosa
-    USE disvert_mod
+    USE disvert_mod, ONLY : caldyn_eta, eta_mass
     USE vorticity_mod
     USE kinetic_mod
@@ -285 +285 @@
 SUBROUTINE compute_pvort(ps,u,theta_rhodz, rhodz,theta,qu,qv)
   USE icosa
-  USE disvert_mod
+  USE disvert_mod, ONLY : mass_dak, mass_dbk, caldyn_eta, eta_mass
   USE exner_mod
   USE trace
@@ -318 +318 @@
            DO i=ii_begin-1,ii_end+1
               ij=(j-1)*iim+i
-              m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g 
-!              m = ( mass_dak(l)+ps(ij)*mass_dbk(l) )/g
+              m = ( mass_dak(l)+ps(ij)*mass_dbk(l) )/g
               rhodz(ij,l) = m
               theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
@@ -420 +419 @@
        ENDDO
 
-    ELSE
-
+    ELSE 
        ! We are using a Lagrangian vertical coordinate
        ! Pressure must be computed first top-down (temporarily stored in pk)
        ! Then Exner pressure and geopotential are computed bottom-up
-       ! Notice that the computation below should work also when caldyn_eta=eta_mass
+       ! Notice that the computation below should work also when caldyn_eta=eta_mass          
 
        ! uppermost layer
@@ -451 +449 @@
           END DO
        END DO
-       DO l = 1,llm
-          !$OMP DO SCHEDULE(STATIC) 
-          DO j=jj_begin-1,jj_end+1
-             DO i=ii_begin-1,ii_end+1
-                ij=(j-1)*iim+i
-                p_ik = pk(ij,l)
-                exner_ik = cpp * (p_ik/preff) ** kappa
-                ! 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 
-                pk(ij,l) = exner_ik
+
+       IF(boussinesq) THEN ! specific volume 1 = dphi/g/rhodz
+          DO l = 1,llm
+             !$OMP DO SCHEDULE(STATIC) 
+             DO j=jj_begin-1,jj_end+1
+                DO i=ii_begin-1,ii_end+1
+                   ij=(j-1)*iim+i
+                   geopot(ij,l+1) = geopot(ij,l) + g*rhodz(ij,l)
+                ENDDO
              ENDDO
           ENDDO
-       ENDDO
+       ELSE ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz
+          DO l = 1,llm
+             !$OMP DO SCHEDULE(STATIC) 
+             DO j=jj_begin-1,jj_end+1
+                DO i=ii_begin-1,ii_end+1
+                   ij=(j-1)*iim+i
+                   p_ik = pk(ij,l)
+                   exner_ik = cpp * (p_ik/preff) ** kappa
+                   geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik 
+                   pk(ij,l) = exner_ik
+                ENDDO
+             ENDDO
+          ENDDO
+       END IF
 
     END IF
@@ -482 +492 @@
     REAL(rstd),INTENT(IN)  :: qu(iim*3*jjm,llm)
     REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm)  ! potential temperature
-    REAL(rstd),INTENT(IN)  :: pk(iim*jjm,llm) ! Exner function
+    REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm) ! Exner function
     REAL(rstd),INTENT(IN)  :: geopot(iim*jjm,llm+1)    ! geopotential
 
@@ -647 +657 @@
   
 !!! Compute bernouilli term = Kinetic Energy + geopotential
-   DO l=ll_begin,ll_end
-    DO j=jj_begin,jj_end
-      DO i=ii_begin,ii_end
-        ij=(j-1)*iim+i
-
-        berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
+    IF(boussinesq) THEN
+
+       DO l=ll_begin,ll_end
+          DO j=jj_begin,jj_end
+             DO i=ii_begin,ii_end
+                ij=(j-1)*iim+i
+                
+                berni(ij,l) = pk(ij,l) + &
+                     + 1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +    &
+                                     le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +          &
+                                     le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +          &
+                                     le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +       &
+                                     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 )  
+                pk(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
+             ENDDO
+          ENDDO
+       ENDDO
+
+    ELSE
+
+       DO l=ll_begin,ll_end
+          DO j=jj_begin,jj_end
+             DO i=ii_begin,ii_end
+                ij=(j-1)*iim+i
+                
+                berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
                      + 1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +    &
                                      le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +          &
@@ -660 +691 @@
                                      le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 )  
        
-      ENDDO
-    ENDDO
-  ENDDO
-
- 
+             ENDDO
+          ENDDO
+       ENDDO
+
+    END IF ! Boussinesq/compressible
+
 !!! Add gradients of Bernoulli and Exner functions to du
-
-   DO l=ll_begin,ll_end
-    DO j=jj_begin,jj_end
-      DO i=ii_begin,ii_end
-        ij=(j-1)*iim+i
+    DO l=ll_begin,ll_end
+       DO j=jj_begin,jj_end
+          DO i=ii_begin,ii_end
+             ij=(j-1)*iim+i
         
-        du(ij+u_right,l) = du(ij+u_right,l) + 1/de(ij+u_right) * (             &
-                          0.5*(theta(ij,l)+theta(ij+t_right,l))                &
-                             *( ne_right*pk(ij,l)+ne_left*pk(ij+t_right,l))    &
-                        + ne_right*berni(ij,l)+ne_left*berni(ij+t_right,l) )
+             du(ij+u_right,l) = du(ij+u_right,l) + 1/de(ij+u_right) * (             &
+                               0.5*(theta(ij,l)+theta(ij+t_right,l))                &
+                                  *( ne_right*pk(ij,l)+ne_left*pk(ij+t_right,l))    &
+                                   + ne_right*berni(ij,l)+ne_left*berni(ij+t_right,l) )
         
         
-        du(ij+u_lup,l) = du(ij+u_lup,l) +  1/de(ij+u_lup) * (                  &
-                          0.5*(theta(ij,l)+theta(ij+t_lup,l))                  &
-                             *( ne_lup*pk(ij,l)+ne_rdown*pk(ij+t_lup,l))       &
-                        + ne_lup*berni(ij,l)+ne_rdown*berni(ij+t_lup,l) )
+             du(ij+u_lup,l) = du(ij+u_lup,l) +  1/de(ij+u_lup) * (                  &
+                               0.5*(theta(ij,l)+theta(ij+t_lup,l))                  &
+                                  *( ne_lup*pk(ij,l)+ne_rdown*pk(ij+t_lup,l))       &
+                                   + ne_lup*berni(ij,l)+ne_rdown*berni(ij+t_lup,l) )
                 
-        du(ij+u_ldown,l) = du(ij+u_ldown,l) + 1/de(ij+u_ldown) * (             &
-                          0.5*(theta(ij,l)+theta(ij+t_ldown,l))                &
-                             *( ne_ldown*pk(ij,l)+ne_rup*pk(ij+t_ldown,l))     &
-                        + ne_ldown*berni(ij,l)+ne_rup*berni(ij+t_ldown,l) )
-
-      ENDDO
+             du(ij+u_ldown,l) = du(ij+u_ldown,l) + 1/de(ij+u_ldown) * (             &
+                               0.5*(theta(ij,l)+theta(ij+t_ldown,l))                &
+                                  *( ne_ldown*pk(ij,l)+ne_rup*pk(ij+t_ldown,l))     &
+                                   + ne_ldown*berni(ij,l)+ne_rup*berni(ij+t_ldown,l) )
+
+          ENDDO
+       ENDDO
     ENDDO
-  ENDDO
-
-  CALL trace_end("compute_caldyn_horiz")
+ 
+    CALL trace_end("compute_caldyn_horiz")
 
 END SUBROUTINE compute_caldyn_horiz