Changeset 375 for codes/icosagcm/trunk/src

Timestamp:

04/08/16 09:48:36 (8 years ago)

Author:

dubos

Message:

Cosmetic rewrite of gradient operators

File:

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

codes/icosagcm/trunk/src/caldyn_kernels_hevi.f90

@@ -357 +357 @@
     !    ENDDO
     
-    ! Compute Exner function (needed by compute_caldyn_fast) and du=g^-2.grad(w^2)
+    ! Compute Exner function (needed by compute_caldyn_fast) and du=-g^2.grad(w^2)
     DO l=1,llm
        !DIR$ SIMD
        DO ij=ij_begin_ext,ij_end_ext
           pk(ij,l) = cpp*((pk(ij,l)/preff)**kappa) ! other formulae possible if exponentiation is slow
-          berni(ij) = (-.5*g*g)*( (W(ij,l)/m_il(ij,l))**2 &
+          berni(ij) = (-.5*g*g)*(              &
+                 (W(ij,l)/m_il(ij,l))**2       &
                + (W(ij,l+1)/m_il(ij,l+1))**2 )
        ENDDO
        DO ij=ij_begin,ij_end         
-          du(ij+u_right,l) = ne_right*berni(ij)+ne_left *berni(ij+t_right)
-          du(ij+u_lup,l)   = ne_lup  *berni(ij)+ne_rdown*berni(ij+t_lup)
-          du(ij+u_ldown,l) = ne_ldown*berni(ij)+ne_rup  *berni(ij+t_ldown)
+          du(ij+u_right,l) = ne_right*(berni(ij)-berni(ij+t_right))
+          du(ij+u_lup,l)   = ne_lup  *(berni(ij)-berni(ij+t_lup))
+          du(ij+u_ldown,l) = ne_ldown*(berni(ij)-berni(ij+t_ldown))
        ENDDO
     ENDDO
@@ -418 +419 @@
        !DIR$ SIMD
        DO ij=ij_begin,ij_end         
-          due_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)
-          due_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)
-          due_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)
-          du(ij+u_right,l) = du(ij+u_right,l) + due_right
-          du(ij+u_lup,l)   = du(ij+u_lup,l)   + due_lup
-          du(ij+u_ldown,l) = du(ij+u_ldown,l) + due_ldown
+          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)
@@ -623 +624 @@
                le_de(ij+u_rdown)*u(ij+u_rdown,l)**2 )  
        ENDDO
-       ! Compute du=grad(Bernoulli)
+       ! Compute du=-grad(Bernoulli)
        !DIR$ SIMD
        DO ij=ij_begin,ij_end
-             du(ij+u_right,l) = ne_right*berni(ij)+ne_left*berni(ij+t_right)
-             du(ij+u_lup,l)   = ne_lup*berni(ij)+ne_rdown*berni(ij+t_lup)
-             du(ij+u_ldown,l) = ne_ldown*berni(ij)+ne_rup*berni(ij+t_ldown)
+             du(ij+u_right,l) = ne_right*(berni(ij)-berni(ij+t_right))
+             du(ij+u_lup,l)   = ne_lup*(berni(ij)-berni(ij+t_lup))
+             du(ij+u_ldown,l) = ne_ldown*(berni(ij)-berni(ij+t_ldown))
        END DO
     END DO
@@ -751 +752 @@
                     -0.5*(F_el(ij+u_right,l)+F_el(ij+u_right,l+1))
           hflux(ij+u_right,l) = uu_right*le_de(ij+u_right)
-          du(ij+u_right,l) = ne_right*berni(ij)+ne_left*berni(ij+t_right)
+          du(ij+u_right,l) = ne_right*(berni(ij)-berni(ij+t_right))
           ! Compute on edge 'lup'
           uu_lup = 0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l) &
                   -0.5*(F_el(ij+u_lup,l)+F_el(ij+u_lup,l+1))
           hflux(ij+u_lup,l) = uu_lup*le_de(ij+u_lup)
-          du(ij+u_lup,l) = ne_lup*berni(ij)+ne_rdown*berni(ij+t_lup)
+          du(ij+u_lup,l) = ne_lup*(berni(ij)-berni(ij+t_lup))
           ! Compute on edge 'ldown'
           uu_ldown = 0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l) &
                     -0.5*(F_el(ij+u_ldown,l)+F_el(ij+u_ldown,l+1))
           hflux(ij+u_ldown,l) = uu_ldown*le_de(ij+u_ldown)
-          du(ij+u_ldown,l) = ne_ldown*berni(ij)+ne_rup*berni(ij+t_ldown)
+          du(ij+u_ldown,l) = ne_ldown*(berni(ij)-berni(ij+t_ldown))
        END DO
     END DO