Changeset 628 for trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90

Timestamp:

Apr 20, 2012, 12:11:47 PM (13 years ago)

Author:

acolaitis

Message:

Correction to advection problems in thermals + made the thermals model faster by limiting the vertical extension of loops with the height reached by thermals.

File:

• trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90 (modified) (12 diffs)

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

@@ -99 +99 @@
       REAL f(ngridmx)
 
+      REAL detrmod(ngridmx,nlayermx)
+
       REAL teta_th_int(ngridmx,nlayermx)
       REAL teta_env_int(ngridmx,nlayermx)
       REAL teta_down_int(ngridmx,nlayermx)
 
-      CHARACTER (LEN=20) :: modname
       CHARACTER (LEN=80) :: abort_message
+      INTEGER ndt
 
 ! ============= PLUME VARIABLES ============
@@ -123 +125 @@
       REAL denom(ngridmx)
       REAL zlevinter(ngridmx)
+      INTEGER zlmax
 
 ! =========================================
@@ -186 +189 @@
       detr(:,:)=0.
       fm(:,:)=0.
-      zu(:,:)=pu(:,:)
-      zv(:,:)=pv(:,:)
+!      zu(:,:)=pu(:,:)
+!      zv(:,:)=pv(:,:)
       zhc(:,:)=pt(:,:)/zpopsk(:,:)
+      ndt=1
 
 ! **********************************************************************
@@ -888 +892 @@
 ! ===========================================================================
 
+      zlmax=MAXVAL(lmax(:))+2
+      if (zlmax .ge. nlayermx) then
+        print*,'thermals have reached last layer of the model'
+        print*,'this is not good !'
+      endif
+
 ! Choix de la fonction d'alimentation utilisee pour la fermeture.
 
@@ -968 +978 @@
 !-------------------------------------------------------------------------
 
-      do l=1,nlayermx
+      do l=1,zlmax
          entr(:,l)=f(:)*(entr_star(:,l)+alim_star(:,l))
          detr(:,l)=f(:)*detr_star(:,l)
       enddo
 
-      do l=1,nlayermx
+      do l=1,zlmax
          do ig=1,ngridmx
             if (l.lt.lmax(ig)) then
@@ -1007 +1017 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-      do l=1,nlayermx
+      do l=1,zlmax
 
          do ig=1,ngridmx
@@ -1156 +1166 @@
 !-----------------------------------------------------------------------
 
-      do l=1,nlayermx-1
+      do l=1,zlmax
          do ig=1,ngridmx
             eee0=entr(ig,l)
@@ -1347 +1357 @@
 ! Calcul de la fraction de l'ascendance
 !------------------------------------------------------------------
-      do ig=1,ngridmx
-         fraca(ig,1)=0.
-         fraca(ig,nlayermx+1)=0.
-      enddo
-      do l=2,nlayermx
+      fraca(:,:)=0.
+      do l=2,zlmax
          do ig=1,ngridmx
             if (zw2(ig,l).gt.1.e-10) then
@@ -1520 +1527 @@
 
       else
-
-      call thermcell_dqup(ngridmx,nlayermx,ptimestep                &
-     &      ,fm,entr,  &
-     &     masse,zu,pduadj)
-
-      call thermcell_dqup(ngridmx,nlayermx,ptimestep    &
-     &       ,fm,entr,  &
-     &     masse,zv,pdvadj)
+!      detrmod(:,:)=0.
+!      do k=1,zlmax
+!         do ig=1,ngridmx
+!            detrmod(ig,k)=fm(ig,k)-fm(ig,k+1) &
+!     &      +entr(ig,k)
+!            if (detrmod(ig,k).lt.0.) then
+!               entr(ig,k)=entr(ig,k)-detrmod(ig,k)
+!               detrmod(ig,k)=0.
+!            endif
+!         enddo
+!      enddo
+!
+!
+!      call thermcell_dqup(ngridmx,nlayermx,ptimestep                &
+!     &      ,fm,entr,detrmod,  &
+!     &     masse,zu,pduadj,ndt,zlmax)
+!
+!      call thermcell_dqup(ngridmx,nlayermx,ptimestep    &
+!     &       ,fm,entr,detrmod,  &
+!     &     masse,zv,pdvadj,ndt,zlmax)
 
       endif
@@ -1538 +1557 @@
 ! The rest is transported outside the sub-timestep loop
 
+      ratiom(:,:)=1.
+
       if (ico2.ne.0) then
-!      if (nqmx .ne. 0.) then
+      detrmod(:,:)=0.
+      do k=1,zlmax
+         do ig=1,ngridmx
+            detrmod(ig,k)=fm(ig,k)-fm(ig,k+1) &
+     &      +entr(ig,k)
+            if (detrmod(ig,k).lt.0.) then
+               entr(ig,k)=entr(ig,k)-detrmod(ig,k)
+               detrmod(ig,k)=0.
+            endif
+         enddo
+      enddo
+
       call thermcell_dqup(ngridmx,nlayermx,ptimestep     &
-     &     ,fm,entr,  &
-     &    masse,pq(:,:,ico2),pdqadj(:,:,ico2))
-!      endif
+     &     ,fm,entr,detrmod,  &
+     &    masse,pq(:,:,ico2),pdqadj(:,:,ico2),ndt,zlmax)
 
 ! Compute the ratio between theta and theta_m
       
-       do l=1,nlayermx
+       do l=1,zlmax
           do ig=1,ngridmx
              ratiom(ig,l)=1./(A*(pq(ig,l,ico2)+pdqadj(ig,l,ico2)*ptimestep)+B)
           enddo
        enddo
-       else
-          ratiom(:,:)=1.
+
        endif
-
 
 !------------------------------------------------------------------
@@ -1561 +1590 @@
 !------------------------------------------------------------------
 
-      do l=1,nlayermx
+      pdtadj(:,:)=0.
+      do l=1,zlmax
          do ig=1,ngridmx
          pdtadj(ig,l)=(zdthladj(ig,l)+zdthladj_down(ig,l))*zpopsk(ig,l)*ratiom(ig,l)
@@ -1601 +1631 @@
        teta_down_int(ig,nlayermx)=teta_down_int(ig,nlayermx-1)
       enddo
-      do l=1,nlayermx
+        heatFlux(:,:)=0.
+        buoyancyOut(:,:)=0.
+        buoyancyEst(:,:)=0.
+        heatFlux_down(:,:)=0.
+      do l=1,zlmax
        do ig=1,ngridmx
         heatFlux(ig,l)=fm(ig,l)*(teta_th_int(ig,l)-teta_env_int(ig,l))/(rhobarz(ig,l))