Changeset 972 for LMDZ4/trunk/libf/phylmd/thermcell_dq.F90

Timestamp:

Jun 19, 2008, 12:24:22 PM (16 years ago)

Author:

lmdzadmin

Message:

Version thermique FH/CRio
Ajout tests cas physiques non pris en comptes et ajout/enleve prints
Nouvelle routine thermcell_flux2.F90
IM

File:

• LMDZ4/trunk/libf/phylmd/thermcell_dq.F90 (modified) (5 diffs)

LMDZ4/trunk/libf/phylmd/thermcell_dq.F90

@@ -23 +23 @@
       real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
 
+      real zzm
+
       integer ig,k
+      real cfl
+
+      real qold(ngrid,nlay)
+      real ztimestep
+      integer niter,iter
+
+
+
+! Calcul du critere CFL pour l'advection dans la subsidence
+      do k=1,nlay
+         do ig=1,ngrid
+            zzm=masse(ig,k)/ptimestep
+            cfl=max(cfl,fm(ig,k)/zzm)
+            if (entr(ig,k).gt.zzm) then
+               print*,'entr dt > m ',entr(ig,k)*ptimestep,masse(ig,k)
+               stop
+            endif
+         enddo
+      enddo
+
+!IM 090508     print*,'CFL CFL CFL CFL ',cfl
+
+#undef CFL
+#ifdef CFL
+! On subdivise le calcul en niter pas de temps.
+      niter=int(cfl)+1
+#else
+      niter=1
+#endif
+
+      ztimestep=ptimestep/niter
+      qold=q
+
+
+do iter=1,niter
+      if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
 
 !   calcul du detrainement
-
-      if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
-
       do k=1,nlay
          do ig=1,ngrid
@@ -56 +91 @@
       do k=2,nlay
          do ig=1,ngrid
-            if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.  &
+            if ((fm(ig,k+1)+detr(ig,k))*ztimestep.gt.  &
      &         1.e-5*masse(ig,k)) then
          qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k))  &
@@ -72 +107 @@
       enddo
 
+! Calcul du flux subsident
+
       do k=2,nlay
          do ig=1,ngrid
-!             wqd(ig,k)=fm(ig,k)*0.5*(q(ig,k-1)+q(ig,k))
+#undef centre
+#ifdef centre
+             wqd(ig,k)=fm(ig,k)*0.5*(q(ig,k-1)+q(ig,k))
+#else
+
+#define plusqueun
+#ifdef plusqueun
+! Schema avec advection sur plus qu'une maille.
+            zzm=masse(ig,k)/ztimestep
+            if (fm(ig,k)>zzm) then
+               wqd(ig,k)=zzm*q(ig,k)+(fm(ig,k)-zzm)*q(ig,k+1)
+            else
+               wqd(ig,k)=fm(ig,k)*q(ig,k)
+            endif
+#else
             wqd(ig,k)=fm(ig,k)*q(ig,k)
+#endif
+#endif
+
             if (wqd(ig,k).lt.0.) then
 !               print*,'wqd<0!!!'
@@ -86 +140 @@
       enddo
      
+
+! Calcul des tendances
       do k=1,nlay
          do ig=1,ngrid
-            dq(ig,k)=(detr(ig,k)*qa(ig,k)-entr(ig,k)*q(ig,k)  &
+            q(ig,k)=q(ig,k)+(detr(ig,k)*qa(ig,k)-entr(ig,k)*q(ig,k)  &
      &               -wqd(ig,k)+wqd(ig,k+1))  &
-     &               /masse(ig,k)
+     &               *ztimestep/masse(ig,k)
 !            if (dq(ig,k).lt.0.) then
 !               print*,'dq<0!!!'
@@ -97 +153 @@
       enddo
 
+
+enddo
+
+
+! Calcul des tendances
+      do k=1,nlay
+         do ig=1,ngrid
+            dq(ig,k)=(q(ig,k)-qold(ig,k))/ptimestep
+            q(ig,k)=qold(ig,k)
+         enddo
+      enddo
+
       return
       end