Changeset 1614 for LMDZ5/trunk/libf/dyn3d

Timestamp:

Feb 3, 2012, 11:07:08 AM (12 years ago)

Author:

Ehouarn Millour

Message:

Correction to enforce having 1+1=2 in the dynamics.
Still not sure why changing surface pressure during dissipation step leads to having 1+1!=2. But clearly there is no reason to recompute polar surface pressure during a dissipation step.
EM

File:

• LMDZ5/trunk/libf/dyn3d/leapfrog.F (modified) (5 diffs)

LMDZ5/trunk/libf/dyn3d/leapfrog.F

@@ -221 +221 @@
 c   ----------------------------------
 
-   1  CONTINUE
+   1  CONTINUE ! Matsuno Forward step begins here
 
       jD_cur = jD_ref + day_ini - day_ref +                             &
@@ -261 +261 @@
       CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
 
-   2  CONTINUE
+   2  CONTINUE ! Matsuno backward or leapfrog step begins here
 
 c-----------------------------------------------------------------------
@@ -304 +304 @@
 c   --------------------------------
 
+      ! compute geopotential phi()
       CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi   )
 
@@ -537 +538 @@
         ENDDO
 
-        DO ij =  1,iim
-          tppn(ij)  = aire(  ij    ) * ps (  ij    )
-          tpps(ij)  = aire(ij+ip1jm) * ps (ij+ip1jm)
-        ENDDO
-          tpn  = SSUM(iim,tppn,1)/apoln
-          tps  = SSUM(iim,tpps,1)/apols
-
-        DO ij = 1, iip1
-          ps(  ij    ) = tpn
-          ps(ij+ip1jm) = tps
-        ENDDO
-
+        if (1 == 0) then
+!!! Ehouarn: lines here 1) kill 1+1=2 in the dynamics
+!!!                     2) should probably not be here anyway
+!!! but are kept for those who would want to revert to previous behaviour
+           DO ij =  1,iim
+             tppn(ij)  = aire(  ij    ) * ps (  ij    )
+             tpps(ij)  = aire(ij+ip1jm) * ps (ij+ip1jm)
+           ENDDO
+             tpn  = SSUM(iim,tppn,1)/apoln
+             tps  = SSUM(iim,tpps,1)/apols
+
+           DO ij = 1, iip1
+             ps(  ij    ) = tpn
+             ps(ij+ip1jm) = tps
+           ENDDO
+        endif ! of if (1 == 0)
 
       END IF ! of IF(apdiss)
@@ -664 +669 @@
 
               CLOSE(99)
+              !!! Ehouarn: Why not stop here and now?
             ENDIF ! of IF (itau.EQ.itaufin)