Changeset 620 for trunk/LMDZ.MARS/libf

Timestamp:

Apr 13, 2012, 6:29:11 PM (13 years ago)

Author:

acolaitis

Message:

THERMALS: modified tracer advection scheme with a homemade one, coded in implicit. Tracer conservation is now better than compared to convadj, with constant and sharp tracer profiles

File:

• trunk/LMDZ.MARS/libf/phymars/thermcell_dqup.F90 (modified) (2 diffs)

trunk/LMDZ.MARS/libf/phymars/thermcell_dqup.F90

@@ -41 +41 @@
       REAL gammae(ngridmx,nlayermx)
       REAL alphamid(ngridmx,nlayermx)
-
+      REAL Z(2,2,ngridmx,nlayermx)
+      REAL X(2,1,ngridmx,nlayermx)
 ! =========== Init ==============================================
 
@@ -51 +52 @@
             if (fm0(ig,k)+entr0(ig,k) .gt. 0.) then
               gammaf(ig,k)=fm0(ig,k)/(fm0(ig,k)+entr0(ig,k))
-              gammae(ig,k)=fm0(ig,k+1)/(fm0(ig,k)+entr0(ig,k))
             else
               gammaf(ig,k)=0.
+            endif
+            if (fm0(ig,k+1) .gt. 0.) then
+              gammae(ig,k)=(fm0(ig,k)+entr0(ig,k))/fm0(ig,k+1)
+            else
               gammae(ig,k)=1.
             endif
          enddo
       enddo
+      do k=1,nlayermx
+        do ig=1, ngridmx
+          Z(1,1,ig,k)=1.
+          Z(1,2,ig,k)=0.
+          Z(2,1,ig,k)=0.
+          Z(2,2,ig,k)=1.          
+        enddo
+      enddo
+      do k=2,nlayermx
+        do ig=1, ngridmx
+          if (gammaf(ig,k) .gt. 0.) then
+          Z(1,1,ig,k)=1./gammae(ig,k-1) + (gammaf(ig,k)-1.)*            &
+     &               (1.-gammae(ig,k-1))/(gammae(ig,k-1)*gammaf(ig,k))
+          Z(1,2,ig,k)=-(1.-gammae(ig,k-1))                              &
+     &                        /(gammae(ig,k-1)*gammaf(ig,k))
+          Z(2,1,ig,k)=(gammaf(ig,k)-1.)/gammaf(ig,k)
+          Z(2,2,ig,k)=1./gammaf(ig,k)
+          else
+          Z(1,1,ig,k)=1./gammae(ig,k-1) -1.
+          Z(1,2,ig,k)=1.
+          Z(2,:,ig,k)=0.
+          endif
 
-
-! =========== Updraft ============================================
-
-!      qa(:,1)=q_therm(:,1)
-
-
-      do ig=1,ngridmx
-         do k=2, nlayermx
-
-             qa(ig,k)=gammaf(ig,k)*qa(ig,k-1)+(1.-gammaf(ig,k))*q(ig,k)
-
-         enddo
+        enddo
       enddo
 
-! =========== Environnement =======================================
+      X(1,1,:,:)=q_therm(:,:)
+      X(2,1,:,:)=q_therm(:,:)
 
-!     q(:,nlayermx)=q_therm(:,nlayermx)
+      do k=nlayermx,2,-1
+        do ig=1, ngridmx
+       X(:,1,ig,k-1)=MATMUL(Z(:,:,ig,k),X(:,1,ig,k))
 
-      do ig=1,ngridmx
-         do k=nlayermx-1, 1,-1
- 
-             q(ig,k)=gammae(ig,k)*q(ig,k+1)+(1.-gammae(ig,k))*qa(ig,k)
-
-         enddo
+        enddo
       enddo
 
-! ====== dq ======================================================
+      q(:,:)=X(1,1,:,:)
+      qa(:,:)=X(2,1,:,:)
 
       dq_therm(:,:)=0.