Changeset 1294 for LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_closure.F90

Timestamp:

Jan 14, 2010, 3:35:30 PM (14 years ago)

Author:

Laurent Fairhead

Message:

Modifications pour la nouvelle version des thermiques (2009/2010) CR et FH

File:

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

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

@@ +1 @@
+!
+! $Header$
+!
       SUBROUTINE thermcell_closure(ngrid,nlay,r_aspect,ptimestep,rho,  &
-     &   zlev,lalim,alim_star,alim_star_tot,zmax_sec,wmax_sec,zmax,wmax,f,lev_out)
+     &   zlev,lalim,alim_star,f_star,zmax,wmax,f,lev_out)
 
 !-------------------------------------------------------------------------
 !thermcell_closure: fermeture, determination de f
+!
+! Modification 7 septembre 2009
+! 1. On enleve alim_star_tot des arguments pour le recalculer et etre ainis
+! coherent avec l'integrale au numerateur.
+! 2. On ne garde qu'une version des couples wmax,zmax et wmax_sec,zmax_sec
+! l'idee etant que le choix se fasse a l'appel de thermcell_closure
+! 3. Vectorisation en mettant les boucles en l l'exterieur avec des if
 !-------------------------------------------------------------------------
       IMPLICIT NONE
@@ -9 +19 @@
 #include "iniprint.h"
 #include "thermcell.h"
-      INTEGER ngrid,nlay
-      INTEGER ig,k       
-      REAL r_aspect,ptimestep
-      integer lev_out                           ! niveau pour les print
+INTEGER ngrid,nlay
+INTEGER ig,k       
+REAL r_aspect,ptimestep
+integer lev_out                           ! niveau pour les print
 
-      INTEGER lalim(ngrid)
-      REAL alim_star(ngrid,nlay)
-      REAL alim_star_tot(ngrid)
-      REAL rho(ngrid,nlay)
-      REAL zlev(ngrid,nlay)
-      REAL zmax(ngrid),zmax_sec(ngrid)
-      REAL wmax(ngrid),wmax_sec(ngrid)
-      real zdenom
+INTEGER lalim(ngrid)
+REAL alim_star(ngrid,nlay)
+REAL f_star(ngrid,nlay+1)
+REAL rho(ngrid,nlay)
+REAL zlev(ngrid,nlay)
+REAL zmax(ngrid)
+REAL wmax(ngrid)
+REAL zdenom(ngrid)
+REAL alim_star2(ngrid)
+REAL f(ngrid)
 
-      REAL alim_star2(ngrid)
+REAL alim_star_tot(ngrid)
+INTEGER llmax
 
-      REAL f(ngrid)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+print*,'THERMCELL CLOSURE 26E'
 
-      do ig=1,ngrid
-         alim_star2(ig)=0.
-      enddo
-      do ig=1,ngrid
-         if (alim_star(ig,1).LT.1.e-10) then
-            f(ig)=0.
-         else   
-             do k=1,lalim(ig)
-                alim_star2(ig)=alim_star2(ig)+alim_star(ig,k)**2  &
-     &                    /(rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k)))
-             enddo
-             zdenom=max(500.,zmax(ig))*r_aspect*alim_star2(ig)
-             if (zdenom<1.e-14) then
-                print*,'ig=',ig
-                print*,'alim_star2',alim_star2(ig)
-                print*,'zmax',zmax(ig)
-                print*,'r_aspect',r_aspect
-                print*,'zdenom',zdenom
-                print*,'alim_star',alim_star(ig,:)
-                print*,'zmax_sec',zmax_sec(ig)
-                print*,'wmax_sec',wmax_sec(ig)
-                stop
-             endif
-             if ((zmax_sec(ig).gt.1.e-10).and.(iflag_thermals_ed.eq.0)) then 
-             f(ig)=wmax_sec(ig)*alim_star_tot(ig)/(max(500.,zmax_sec(ig))*r_aspect  &
-     &             *alim_star2(ig))
-!            f(ig)=f(ig)+(f0(ig)-f(ig))*exp((-ptimestep/  &
-!    &                     zmax_sec(ig))*wmax_sec(ig))
-             if(prt_level.GE.10) write(lunout,*)'closure dry',f(ig),wmax_sec(ig),alim_star_tot(ig),zmax_sec(ig)
-             else
-             f(ig)=wmax(ig)*alim_star_tot(ig)/zdenom
-!            f(ig)=f(ig)+(f0(ig)-f(ig))*exp((-ptimestep/  &
-!     &                     zmax(ig))*wmax(ig))
-             if(prt_level.GE.10) print*,'closure moist',f(ig),wmax(ig),alim_star_tot(ig),zmax(ig)
-             endif
-         endif
-!         f0(ig)=f(ig)
-      enddo
-      if (prt_level.ge.1) print*,'apres fermeture'
+alim_star2(:)=0.
+alim_star_tot(:)=0.
+f(:)=0.
 
-!
+! Indice vertical max (max de lalim) atteint par les thermiques sur le domaine
+llmax=1
+do ig=1,ngrid
+   if (lalim(ig)>llmax) llmax=lalim(ig)
+enddo
+
+
+! Calcul des integrales sur la verticale de alim_star et de
+!   alim_star^2/(rho dz)
+do k=1,llmax-1
+   do ig=1,ngrid
+      if (k<lalim(ig)) then
+         alim_star2(ig)=alim_star2(ig)+alim_star(ig,k)**2  &
+&                    /(rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k)))
+         alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,k)
+      endif
+   enddo
+enddo
+
+
+do ig=1,ngrid
+   if (alim_star2(ig)>1.e-10) then
+      f(ig)=wmax(ig)*alim_star_tot(ig)/  &
+&     (max(500.,zmax(ig))*r_aspect*alim_star2(ig))
+   endif
+enddo
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! TESTS POUR UNE NOUVELLE FERMETURE DANS LAQUELLE ALIM_STAR NE SERAIT
+! PAS NORMALISE
+!           f(ig)=f(ig)*f_star(ig,2)/(f_star(ig,lalim(ig)))
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
       return
       end