Changeset 776 for trunk/LMDZ.COMMON/libf/dyn3d/leapfrog.F

Timestamp:

Sep 7, 2012, 2:49:58 PM (12 years ago)

Author:

emillour

Message:

Common dynamics: updates to keep up with LMDZ5 Earth (rev 1649)
See file "DOC/chantiers/commit_importants.log" for details.
EM

File:

• trunk/LMDZ.COMMON/libf/dyn3d/leapfrog.F (modified) (11 diffs)

trunk/LMDZ.COMMON/libf/dyn3d/leapfrog.F

@@ -124 +124 @@
 
       REAL  SSUM
-      REAL time_0 , finvmaold(ip1jmp1,llm)
+      REAL time_0
+!     REAL finvmaold(ip1jmp1,llm)
 
 cym      LOGICAL  lafin
@@ -243 +244 @@
       dq(:,:,:)=0.
       CALL pression ( ip1jmp1, ap, bp, ps, p       )
-      if (disvert_type==1) then
+      if (pressure_exner) then
         CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
-      else ! we assume that we are in the disvert_type==2 case
+      else
         CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
       endif
@@ -271 +272 @@
 c   ----------------------------------
 
-   1  CONTINUE
+   1  CONTINUE ! Matsuno Forward step begins here
 
       jD_cur = jD_ref + day_ini - day_ref +                             &
-     &          int (itau * dtvr / daysec) 
+     &          itau/day_step 
       jH_cur = jH_ref + start_time +                                    &
-     &          (itau * dtvr / daysec - int(itau * dtvr / daysec)) 
+     &          mod(itau,day_step)/float(day_step) 
       jD_cur = jD_cur + int(jH_cur)
       jH_cur = jH_cur - int(jH_cur)
@@ -307 +308 @@
 
 c   ...    P.Le Van .26/04/94  ....
-
-      CALL SCOPY   ( ijp1llm,   masse, 1, finvmaold,     1 )
-      CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
-
-   2  CONTINUE
+! Ehouarn: finvmaold is actually not used
+!      CALL SCOPY   ( ijp1llm,   masse, 1, finvmaold,     1 )
+!      CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
+
+   2  CONTINUE ! Matsuno backward or leapfrog step begins here
 
 c-----------------------------------------------------------------------
@@ -357 +358 @@
       call tpot2t(ijp1llm,teta,temp,pk)
       tsurpk = cpp*temp/pk
+      ! compute geopotential phi()
       CALL geopot  ( ip1jmp1, tsurpk  , pk , pks,  phis  , phi   )
 
@@ -372 +374 @@
 
 !      IF( forward. OR . leapf )  THEN
+! Ehouarn: NB: at this point p with ps are not synchronized
+!              (whereas mass and ps are...)
       IF((.not.forward).OR. leapf )  THEN
         ! Ehouarn: gather mass fluxes during backward Matsuno or LF step
@@ -398 +402 @@
 
        CALL integrd ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
-     $         dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ,
-     $              finvmaold                                    )
+     $         dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis )
+!     $              finvmaold                                    )
 
        IF ((planet_type.eq."titan").and.(tidal)) then
@@ -431 +435 @@
 
          CALL pression (  ip1jmp1, ap, bp, ps,  p      )
-         if (disvert_type==1) then
+         if (pressure_exner) then
            CALL exner_hyb(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
-         else ! we assume that we are in the disvert_type==2 case
+         else
            CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
          endif
 
            jD_cur = jD_ref + day_ini - day_ref +                        &
-     &          int (itau * dtvr / daysec) 
+     &          itau/day_step 
            jH_cur = jH_ref + start_time +                               &
-     &          (itau * dtvr / daysec - int(itau * dtvr / daysec)) 
+     &          mod(itau,day_step)/float(day_step) 
            jD_cur = jD_cur + int(jH_cur)
            jH_cur = jH_cur - int(jH_cur)
@@ -545 +549 @@
 
         CALL pression ( ip1jmp1, ap, bp, ps, p                  )
-        if (disvert_type==1) then
+        if (pressure_exner) then
           CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
-        else ! we assume that we are in the disvert_type==2 case
+        else
           CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
         endif
@@ -613 +617 @@
         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)
@@ -749 +757 @@
 
               CLOSE(99)
+              !!! Ehouarn: Why not stop here and now?
             ENDIF ! of IF (itau.EQ.itaufin)