Changeset 1403 for LMDZ4/trunk/libf/phylmd/thermcell_dv2.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_dv2.F90 (modified) (6 diffs)

LMDZ4/trunk/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