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

Timestamp:

06/08/16 22:24:19 (8 years ago)

Author:

dubos

Message:

Moist dynamics (to be tested)

File:

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

codes/icosagcm/trunk/src/caldyn_kernels_hevi.f90

@@ -382 +382 @@
     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(IN)    :: theta(iim*jjm,llm,nqdyn)
     REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm)
     REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1)
     REAL(rstd),INTENT(INOUT)   :: du(iim*3*jjm,llm)
     REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function
+    REAL(rstd) :: berniv(iim*jjm,llm)  ! moist Bernoulli function
 
     INTEGER :: i,j,ij,l
-    REAL(rstd) :: due_right, due_lup, due_ldown
+    REAL(rstd) :: Rd, qv, temp, chi, nu, due_right, due_lup, due_ldown
 
     CALL trace_start("compute_caldyn_fast")
+
+    Rd=cpp*kappa
 
     ! Compute Bernoulli term
     IF(boussinesq) THEN
-
        DO l=ll_begin,ll_end
           !DIR$ SIMD
@@ -417 +419 @@
              DO ij=ij_begin,ij_end         
                 berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
-                     + pk(ij,l)*(cpp-theta(ij,l)) ! pk=temperature, theta=entropy
+                     + pk(ij,l)*(cpp-theta(ij,l,1)) ! pk=temperature, theta=entropy
              END DO
+          CASE(thermo_moist) 
+             !DIR$ SIMD
+             DO ij=ij_begin,ij_end         
+                ! du/dt = grad(Bd)+rv.grad(Bv)+s.grad(T)
+                ! Bd = Phi + gibbs_d
+                ! Bv = Phi + gibbs_v
+                ! pk=temperature, theta=entropy
+                qv = theta(ij,l,2)
+                temp = pk(ij,l)
+                chi = log(temp/Treff)
+                nu = (chi*(cpp+qv*cppv)-theta(ij,l,1))/(Rd+qv*Rv) ! log(p/preff)
+                berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
+                     + temp*(cpp*(1.-chi)+Rd*nu)
+                berniv(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
+                     + temp*(cppv*(1.-chi)+Rv*nu)
+            END DO
           END SELECT
        END DO
@@ -426 +444 @@
 !!! u:=u+tau*du, du = -grad(B)-theta.grad(pi)
     DO l=ll_begin,ll_end
-       !DIR$ SIMD
-       DO ij=ij_begin,ij_end         
-          due_right =  0.5*(theta(ij,l)+theta(ij+t_right,l))  &
-                          *(pk(ij+t_right,l)-pk(ij,l))        &
-                         +  berni(ij+t_right,l)-berni(ij,l)
-          due_lup = 0.5*(theta(ij,l)+theta(ij+t_lup,l))    &
-                       *(pk(ij+t_lup,l)-pk(ij,l))          &
-                         +  berni(ij+t_lup,l)-berni(ij,l)
-          due_ldown = 0.5*(theta(ij,l)+theta(ij+t_ldown,l)) &
-                         *(pk(ij+t_ldown,l)-pk(ij,l))      &
-                        +   berni(ij+t_ldown,l)-berni(ij,l)
-          du(ij+u_right,l) = du(ij+u_right,l) - ne_right*due_right
-          du(ij+u_lup,l)   = du(ij+u_lup,l)   - ne_lup*due_lup
-          du(ij+u_ldown,l) = du(ij+u_ldown,l) - ne_ldown*due_ldown
-          u(ij+u_right,l) = u(ij+u_right,l) + tau*du(ij+u_right,l)
-          u(ij+u_lup,l)   = u(ij+u_lup,l)   + tau*du(ij+u_lup,l)
-          u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*du(ij+u_ldown,l)
-       ENDDO
-    ENDDO
-
+       IF(caldyn_thermo == thermo_moist) THEN
+          !DIR$ SIMD
+          DO ij=ij_begin,ij_end         
+             due_right =  berni(ij+t_right,l)-berni(ij,l)       &
+                  + 0.5*(theta(ij,l,1)+theta(ij+t_right,l,1))   &
+                       *(pk(ij+t_right,l)-pk(ij,l))             &
+                  + 0.5*(theta(ij,l,2)+theta(ij+t_right,l,2))   &
+                       *(berniv(ij+t_right,l)-berniv(ij,l))
+             
+             due_lup = berni(ij+t_lup,l)-berni(ij,l)            &
+                  + 0.5*(theta(ij,l,1)+theta(ij+t_lup,l,1))     &
+                       *(pk(ij+t_lup,l)-pk(ij,l))               &
+                  + 0.5*(theta(ij,l,2)+theta(ij+t_lup,l,2))     &
+                       *(berniv(ij+t_lup,l)-berniv(ij,l))
+             
+             due_ldown = berni(ij+t_ldown,l)-berni(ij,l)        &
+                  + 0.5*(theta(ij,l,1)+theta(ij+t_ldown,l,1))   &
+                       *(pk(ij+t_ldown,l)-pk(ij,l))             &
+                  + 0.5*(theta(ij,l,2)+theta(ij+t_ldown,l,2))   &
+                       *(berniv(ij+t_ldown,l)-berniv(ij,l))
+             
+             du(ij+u_right,l) = du(ij+u_right,l) - ne_right*due_right
+             du(ij+u_lup,l)   = du(ij+u_lup,l)   - ne_lup*due_lup
+             du(ij+u_ldown,l) = du(ij+u_ldown,l) - ne_ldown*due_ldown
+             u(ij+u_right,l) = u(ij+u_right,l) + tau*du(ij+u_right,l)
+             u(ij+u_lup,l)   = u(ij+u_lup,l)   + tau*du(ij+u_lup,l)
+             u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*du(ij+u_ldown,l)
+          END DO
+       ELSE
+          !DIR$ SIMD
+          DO ij=ij_begin,ij_end         
+             due_right =  0.5*(theta(ij,l,1)+theta(ij+t_right,l,1))  &
+                  *(pk(ij+t_right,l)-pk(ij,l))        &
+                  +  berni(ij+t_right,l)-berni(ij,l)
+             due_lup = 0.5*(theta(ij,l,1)+theta(ij+t_lup,l,1))    &
+                  *(pk(ij+t_lup,l)-pk(ij,l))          &
+                  +  berni(ij+t_lup,l)-berni(ij,l)
+             due_ldown = 0.5*(theta(ij,l,1)+theta(ij+t_ldown,l,1)) &
+                  *(pk(ij+t_ldown,l)-pk(ij,l))      &
+                  +   berni(ij+t_ldown,l)-berni(ij,l)
+             du(ij+u_right,l) = du(ij+u_right,l) - ne_right*due_right
+             du(ij+u_lup,l)   = du(ij+u_lup,l)   - ne_lup*due_lup
+             du(ij+u_ldown,l) = du(ij+u_ldown,l) - ne_ldown*due_ldown
+             u(ij+u_right,l) = u(ij+u_right,l) + tau*du(ij+u_right,l)
+             u(ij+u_lup,l)   = u(ij+u_lup,l)   + tau*du(ij+u_lup,l)
+             u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*du(ij+u_ldown,l)
+          END DO
+       END IF
+    END DO
+       
     CALL trace_end("compute_caldyn_fast")