Changeset 1403 for LMDZ4/trunk/libf/phylmd/thermcell_env.F90

Timestamp:

Jul 1, 2010, 11:02:53 AM (14 years ago)

Author:

Laurent Fairhead

Message:

Merged LMDZ4V5.0-dev branch changes r1292:r1399 to trunk.

Validation:
Validation consisted in compiling the HEAD revision of the trunk,
LMDZ4V5.0-dev branch and the merged sources and running different
configurations on local and SX8 machines comparing results.

Local machine: bench configuration, 32x24x11, gfortran

• IPSLCM5A configuration (comparison between trunk and merged sources):
  • numerical convergence on dynamical fields over 3 days
  • start files are equivalent (except for RN and PB fields)
  • daily history files equivalent

• MH07 configuration, new physics package (comparison between LMDZ4V5.0-dev branch and merged sources):
  • numerical convergence on dynamical fields over 3 days
  • start files are equivalent (except for RN and PB fields)
  • daily history files equivalent

SX8 machine (brodie), 96x95x39 on 4 processors:

• IPSLCM5A configuration:
  • start files are equivalent (except for RN and PB fields)
  • monthly history files equivalent

• MH07 configuration:
  • start files are equivalent (except for RN and PB fields)
  • monthly history files equivalent

Changes to the makegcm and create_make_gcm scripts to take into account
main programs in F90 files

---

Fusion de la branche LMDZ4V5.0-dev (r1292:r1399) au tronc principal

Validation:
La validation a consisté à compiler la HEAD de le trunk et de la banche
LMDZ4V5.0-dev et les sources fusionnées et de faire tourner le modéle selon
différentes configurations en local et sur SX8 et de comparer les résultats

En local: 32x24x11, config bench/gfortran

• pour une config IPSLCM5A (comparaison tronc/fusion):
  • convergence numérique sur les champs dynamiques après 3 jours
  • restart et restartphy égaux (à part sur RN et Pb)
  • fichiers histoire égaux

• pour une config nlle physique (MH07) (comparaison LMDZ4v5.0-dev/fusion):
  • convergence numérique sur les champs dynamiques après 3 jours
  • restart et restartphy égaux
  • fichiers histoire équivalents

Sur brodie, 96x95x39 sur 4 proc:

• pour une config IPSLCM5A:
  • restart et restartphy égaux (à part sur RN et PB)
  • pas de différence dans les fichiers histmth.nc

• pour une config MH07
  • restart et restartphy égaux (à part sur RN et PB)
  • pas de différence dans les fichiers histmth.nc

Changement sur makegcm et create_make-gcm pour pouvoir prendre en compte des
programmes principaux en *F90

Location:

LMDZ4/trunk

Files:

• . (modified) (1 prop)
• libf/phylmd/thermcell_env.F90 (modified) (3 diffs)

LMDZ4/trunk/libf/phylmd/thermcell_env.F90

@@ -1 +1 @@
       SUBROUTINE thermcell_env(ngrid,nlay,po,pt,pu,pv,pplay,  &
-     &           pplev,zo,zh,zl,ztv,zthl,zu,zv,zpspsk,zqsat,lev_out)
+     &           pplev,zo,zh,zl,ztv,zthl,zu,zv,zpspsk,pqsat,lev_out)
 
 !--------------------------------------------------------------
@@ -31 +31 @@
       REAL zu(ngrid,nlay)
       REAL zv(ngrid,nlay)
-      REAL zqsat(ngrid,nlay)
+      REAL pqsat(ngrid,nlay)
 
-      INTEGER ig,l,ll
+      INTEGER ig,ll
 
-      real zcor,zdelta,zcvm5,qlbef
-      real Tbef,qsatbef
-      real dqsat_dT,DT,num,denom
-      REAL RLvCp,DDT0
-      PARAMETER (DDT0=.01)
-      LOGICAL Zsat
+      real dqsat_dT
+      real RLvCp
 
-      Zsat=.false.
-      RLvCp = RLVTT/RCPD
+logical mask(ngrid,nlay)
 
-!
-! Pr Tprec=Tl calcul de qsat 
-! Si qsat>qT T=Tl, q=qT
-! Sinon DDT=(-Tprec+Tl+RLVCP (qT-qsat(T')) / (1+RLVCP dqsat/dt) 
-! On cherche DDT < DDT0
+
+!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+! Initialisations :
+!------------------
+
+mask(:,:)=.true.
+RLvCp = RLVTT/RCPD
+
 !
 ! calcul des caracteristiques de l environnement
@@ -57 +55 @@
             zl(ig,ll)=0.
             zh(ig,ll)=pt(ig,ll)
-            zqsat(ig,ll)=0.
          EndDO
        EndDO
 !
 !
-!recherche de saturation dans l environnement
-       DO ll=1,nlay
-! les points insatures sont definitifs
-         DO ig=1,ngrid
-            Tbef=pt(ig,ll)
-            zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
-            qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,ll)
-            qsatbef=MIN(0.5,qsatbef)
-            zcor=1./(1.-retv*qsatbef)
-            qsatbef=qsatbef*zcor
-            Zsat = (max(0.,po(ig,ll)-qsatbef) .gt. 1.e-10)
-            if (Zsat) then
-            qlbef=max(0.,po(ig,ll)-qsatbef)
-! si sature: ql est surestime, d'ou la sous-relax
-            DT = 0.5*RLvCp*qlbef
-! on pourra enchainer 2 ou 3 calculs sans Do while
-            do while (abs(DT).gt.DDT0)
-! il faut verifier si c,a conserve quand on repasse en insature ...
-              Tbef=Tbef+DT
-              zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
-              qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,ll)
-              qsatbef=MIN(0.5,qsatbef)
-              zcor=1./(1.-retv*qsatbef)
-              qsatbef=qsatbef*zcor
-! on veut le signe de qlbef
-              qlbef=po(ig,ll)-qsatbef
-!          dqsat_dT
-              zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
-              zcvm5=R5LES*(1.-zdelta) + R5IES*zdelta
-              zcor=1./(1.-retv*qsatbef)
-              dqsat_dT=FOEDE(Tbef,zdelta,zcvm5,qsatbef,zcor)
-              num=-Tbef+pt(ig,ll)+RLvCp*qlbef
-              denom=1.+RLvCp*dqsat_dT
-              if (denom.lt.1.e-10) then
-                  print*,'pb denom'
-              endif
-              DT=num/denom
-            enddo
-! on ecrit de maniere conservative (sat ou non)
-            zl(ig,ll) = max(0.,qlbef)
-!          T = Tl +Lv/Cp ql
-            zh(ig,ll) = pt(ig,ll)+RLvCp*zl(ig,ll)
-            zo(ig,ll) = po(ig,ll)-zl(ig,ll)
-           endif
-!on ecrit zqsat 
-            zqsat(ig,ll)=qsatbef     
-         EndDO
-       EndDO
+! Condensation :
+!---------------
+! Calcul de l'humidite a saturation et de la condensation
+
+call thermcell_qsat(ngrid*nlay,mask,pplev,pt,po,pqsat)
+DO ll=1,nlay
+   DO ig=1,ngrid
+      zl(ig,ll) = max(0.,po(ig,ll)-pqsat(ig,ll))
+      zh(ig,ll) = pt(ig,ll)+RLvCp*zl(ig,ll)         !   T = Tl + Lv/Cp ql
+      zo(ig,ll) = po(ig,ll)-zl(ig,ll)
+   ENDDO
+ENDDO
 !
 !
 !-----------------------------------------------------------------------
-!   incrementation eventuelle de tendances precedentes:
-!   ---------------------------------------------------
 
       if (prt_level.ge.1) print*,'0 OK convect8'
 
-      DO 1010 l=1,nlay
-         DO 1015 ig=1,ngrid
-             zpspsk(ig,l)=(pplay(ig,l)/100000.)**RKAPPA
-             zu(ig,l)=pu(ig,l)
-             zv(ig,l)=pv(ig,l)
+      DO ll=1,nlay
+         DO ig=1,ngrid
+             zpspsk(ig,ll)=(pplay(ig,ll)/100000.)**RKAPPA
+             zu(ig,ll)=pu(ig,ll)
+             zv(ig,ll)=pv(ig,ll)
 !attention zh est maintenant le profil de T et plus le profil de theta !
+! Quelle horreur ! A eviter.
 !
 !   T-> Theta
-            ztv(ig,l)=zh(ig,l)/zpspsk(ig,l)
+            ztv(ig,ll)=zh(ig,ll)/zpspsk(ig,ll)
 !Theta_v
-            ztv(ig,l)=ztv(ig,l)*(1.+RETV*(zo(ig,l))  &
-     &           -zl(ig,l))
+            ztv(ig,ll)=ztv(ig,ll)*(1.+RETV*(zo(ig,ll))-zl(ig,ll))
 !Thetal
-            zthl(ig,l)=pt(ig,l)/zpspsk(ig,l)
+            zthl(ig,ll)=pt(ig,ll)/zpspsk(ig,ll)
 !            
-1015     CONTINUE
-1010  CONTINUE
+         ENDDO
+      ENDDO
  
       RETURN