Changeset 1326 for LMDZ4/branches

Timestamp:

Mar 15, 2010, 5:37:02 PM (14 years ago)

Author:

idelkadi

Message:

Optimisation de thermcell_dv2.F90

File:

• LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_dv2.F90 (modified) (6 diffs)

LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_dv2.F90

@@ -10 +10 @@
 !   de "thermiques" explicitement representes
 !   calcul du dq/dt une fois qu'on connait les ascendances
+!
+! Vectorisation, FH : 2010/03/08
 !
 !=======================================================================
@@ -31 +33 @@
       real qa(ngrid,nlay),detr(ngrid,nlay),zf,zf2
       real wvd(ngrid,nlay+1),wud(ngrid,nlay+1)
-      real gamma0,gamma(ngrid,nlay+1)
+      real gamma0(ngrid,nlay+1),gamma(ngrid,nlay+1)
       real ue(ngrid,nlay),ve(ngrid,nlay)
-      real dua,dva
+      LOGICAL ltherm(ngrid,nlay)
+      real dua(ngrid,nlay),dva(ngrid,nlay)
       integer iter
 
-      integer ig,k
+      integer ig,k,nlarga0
+
+!-------------------------------------------------------------------------
 
 !   calcul du detrainement
+!---------------------------
+
+      print*,'THERMCELL DV2 OPTIMISE 3'
+
+      nlarga0=0.
 
       do k=1,nlay
@@ -59 +69 @@
       do k=2,nlay
          do ig=1,ngrid
-            if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.  &
-     &         1.e-5*masse(ig,k)) then
-!   On itère sur la valeur du coeff de freinage.
-!              gamma0=rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k))
-!IM 060508 beg
-!             if(0.5*(fraca(ig,k+1)+fraca(ig,k)).LT.0.) THEN
-!              print*,'th_dv2 ig k fraca(:,k) fraca(:k+1)', &
-!    &         ig,k,fraca(ig,k),fraca(ig,k+1)
-!             endif
-!             if(larga(ig).EQ.0.) THEN
-!              print*,'th_dv2 ig larga=0.',ig
-!             endif
-              if(larga(ig).GT.0.) THEN
-!IM 060508 end
-               gamma0=masse(ig,k)  &
+            ltherm(ig,k)=(fm(ig,k+1)+detr(ig,k))*ptimestep > 1.e-5*masse(ig,k)
+            if(ltherm(ig,k).and.larga(ig)>0.) then
+               gamma0(ig,k)=masse(ig,k)  &
      &         *sqrt( 0.5*(fraca(ig,k+1)+fraca(ig,k)) )  &
      &         *0.5/larga(ig)  &
      &         *1.
-!IM 060508 beg
-              else 
-               if(prt_level.GE.10) print*,'WARNING cas ELSE on initialise gamma0=0.'
-               gamma0=0.
-              endif !(larga(ig).GT.0.) THEN
-!IM 060508 end
-!    s         *0.5
-!              gamma0=0.
-               zf=0.5*(fraca(ig,k)+fraca(ig,k+1))
+            else
+               gamma0(ig,k)=0.
+            endif
+            if (ltherm(ig,k).and.larga(ig)<=0.) nlarga0=nlarga0+1
+         enddo
+      enddo
+
+      gamma(:,:)=0.
+
+      do k=2,nlay
+
+         do ig=1,ngrid
+            if (ltherm(ig,k)) then
+               dua(ig,k)=ua(ig,k-1)-u(ig,k-1)
+               dva(ig,k)=va(ig,k-1)-v(ig,k-1)
+            else
+               ua(ig,k)=u(ig,k)
+               va(ig,k)=v(ig,k)
+               ue(ig,k)=u(ig,k)
+               ve(ig,k)=v(ig,k)
+            endif
+         enddo
+
+
+! Debut des iterations
+!----------------------
+do iter=1,5
+         do ig=1,ngrid
+! Pour memoire : calcul prenant en compte la fraction reelle
+!              zf=0.5*(fraca(ig,k)+fraca(ig,k+1))
+!              zf2=1./(1.-zf)
+! Calcul avec fraction infiniement petite
                zf=0.
-               zf2=1./(1.-zf)
-!   la première fois on multiplie le coefficient de freinage
-!   par le module du vent dans la couche en dessous.
-               dua=ua(ig,k-1)-u(ig,k-1)
-               dva=va(ig,k-1)-v(ig,k-1)
-               do iter=1,5
+               zf2=1.
+
+!  la première fois on multiplie le coefficient de freinage
+!  par le module du vent dans la couche en dessous.
+!  Mais pourquoi donc ???
+               if (ltherm(ig,k)) then
 !   On choisit une relaxation lineaire.
-                  gamma(ig,k)=gamma0
+!                 gamma(ig,k)=gamma0(ig,k)
 !   On choisit une relaxation quadratique.
-                  gamma(ig,k)=gamma0*sqrt(dua**2+dva**2)
+                  gamma(ig,k)=gamma0(ig,k)*sqrt(dua(ig,k)**2+dva(ig,k)**2)
                   ua(ig,k)=(fm(ig,k)*ua(ig,k-1)  &
      &               +(zf2*entr(ig,k)+gamma(ig,k))*u(ig,k))  &
@@ -105 +126 @@
      &               /(fm(ig,k+1)+detr(ig,k)+entr(ig,k)*zf*zf2  &
      &                 +gamma(ig,k))
-!                 print*,k,ua(ig,k),va(ig,k),u(ig,k),v(ig,k),dua,dva
-                  dua=ua(ig,k)-u(ig,k)
-                  dva=va(ig,k)-v(ig,k)
+!                 print*,k,ua(ig,k),va(ig,k),u(ig,k),v(ig,k),dua(ig,k),dva(ig,k)
+                  dua(ig,k)=ua(ig,k)-u(ig,k)
+                  dva(ig,k)=va(ig,k)-v(ig,k)
                   ue(ig,k)=(u(ig,k)-zf*ua(ig,k))*zf2
                   ve(ig,k)=(v(ig,k)-zf*va(ig,k))*zf2
-               enddo
-            else
-               ua(ig,k)=u(ig,k)
-               va(ig,k)=v(ig,k)
-               ue(ig,k)=u(ig,k)
-               ve(ig,k)=v(ig,k)
-               gamma(ig,k)=0.
-            endif
+               endif
          enddo
-      enddo
+! Fin des iterations
+!--------------------
+enddo
 
+      enddo ! k=2,nlay
+
+
+! Calcul du flux vertical de moment dans l'environnement.
+!---------------------------------------------------------
       do k=2,nlay
          do ig=1,ngrid
@@ -134 +155 @@
       enddo
 
+! calcul des tendances.
+!-----------------------
       do k=1,nlay
          do ig=1,ngrid
-!IM
-         if(prt_level.GE.10) print*,'th_dv2 ig k gamma entr detr ua ue va ve wud wvd masse',ig,k,gamma(ig,k), &
-     &   entr(ig,k),detr(ig,k),ua(ig,k),ue(ig,k),va(ig,k),ve(ig,k),wud(ig,k),wvd(ig,k),wud(ig,k+1),wvd(ig,k+1), &
-     &   masse(ig,k)
-!
             du(ig,k)=((detr(ig,k)+gamma(ig,k))*ua(ig,k)  &
      &               -(entr(ig,k)+gamma(ig,k))*ue(ig,k)  &
@@ -152 +170 @@
       enddo
 
+
+! Sorties eventuelles.
+!----------------------
+
+   if(prt_level.GE.10) then
+      do k=1,nlay
+         do ig=1,ngrid
+           print*,'th_dv2 ig k gamma entr detr ua ue va ve wud wvd masse',ig,k,gamma(ig,k), &
+     &   entr(ig,k),detr(ig,k),ua(ig,k),ue(ig,k),va(ig,k),ve(ig,k),wud(ig,k),wvd(ig,k),wud(ig,k+1),wvd(ig,k+1), &
+     &   masse(ig,k)
+         enddo
+      enddo
+   endif
+!
+     if (nlarga0>0) then
+          print*,'WARNING !!!!!! DANS THERMCELL_DV2 '
+          print*,nlarga0,' points pour lesquels laraga=0. dans un thermique'
+          print*,'Il faudrait decortiquer ces points'
+     endif
+
       return
       end