Changeset 1403 for LMDZ4/trunk/libf/phylmd/yamada4.F

Timestamp:

Jul 1, 2010, 11:02:53 AM (15 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/yamada4.F (modified) (7 diffs)

LMDZ4/trunk/libf/phylmd/yamada4.F

@@ -38 +38 @@
 c    iflag_pbl=7 : MY 2.0.Fournier
 c    iflag_pbl=8 : MY 2.5
-c    iflag_pbl=9 : un test ?
+c    iflag_pbl>=9 : MY 2.5 avec diffusion verticale
 
 c.......................................................................
@@ -66 +66 @@
 
       integer nlay,nlev
-cym      PARAMETER (nlay=klev)
-cym      PARAMETER (nlev=klev+1)
 
       logical first
@@ -98 +96 @@
       real fl,zzz,zl0,zq2,zn2
 
-cym      real rino(klon,klev+1),smyam(klon,klev),styam(klon,klev)
-cym     s  ,lyam(klon,klev),knyam(klon,klev)
-cym     s  ,w2yam(klon,klev),t2yam(klon,klev)
-      real,allocatable,save,dimension(:,:) :: rino,smyam,styam,lyam,
-     s                                        knyam,w2yam,t2yam
-cym      common/pbldiag/rino,smyam,styam,lyam,knyam,w2yam,t2yam
-c$OMP THREADPRIVATE(rino,smyam,styam,lyam,knyam,w2yam,t2yam)
+      real rino(klon,klev+1),smyam(klon,klev),styam(klon,klev)
+     s  ,lyam(klon,klev),knyam(klon,klev)
+     s  ,w2yam(klon,klev),t2yam(klon,klev)
       logical,save :: firstcall=.true.
 c$OMP THREADPRIVATE(firstcall)       
@@ -119 +113 @@
       
       if (firstcall) then
-        allocate(rino(klon,klev+1),smyam(klon,klev),styam(klon,klev))
-        allocate(lyam(klon,klev),knyam(klon,klev))
-        allocate(w2yam(klon,klev),t2yam(klon,klev))
         allocate(l0(klon))
         firstcall=.false.
@@ -127 +118 @@
 
 
-      if (.not.(iflag_pbl.ge.6.and.iflag_pbl.le.9)) then
+      if (.not.(iflag_pbl.ge.6.and.iflag_pbl.le.10)) then
            stop'probleme de coherence dans appel a MY'
       endif
@@ -412 +403 @@
       enddo
 
-c     if (iflag_pbl.ge.7..and.0.eq.1) then
-c        q2(:,1)=q2(:,2)
-c        call vdif_q2(dt,g,rconst,plev,temp,kq,q2)
-c     endif
+! Transport diffusif vertical de la TKE.
+      if (iflag_pbl.ge.9) then
+!       print*,'YAMADA VDIF'
+        q2(:,1)=q2(:,2)
+        call vdif_q2(dt,g,rconst,ngrid,plev,temp,kq,q2)
+      endif
 
 c   Traitement des cas noctrunes avec l'introduction d'une longueur
@@ -492 +485 @@
       return
       end
+      SUBROUTINE vdif_q2(timestep,gravity,rconst,ngrid,plev,temp,
+     &  kmy,q2)
+      use dimphy
+      IMPLICIT NONE
+c.......................................................................
+#include "dimensions.h"
+cccc#include "dimphy.h"
+c.......................................................................
+c
+c dt : pas de temps
+
+      real plev(klon,klev+1)
+      real temp(klon,klev)
+      real timestep
+      real gravity,rconst
+      real kstar(klon,klev+1),zz
+      real kmy(klon,klev+1)
+      real q2(klon,klev+1)
+      real deltap(klon,klev+1)
+      real denom(klon,klev+1),alpha(klon,klev+1),beta(klon,klev+1)
+      integer ngrid
+
+      integer i,k
+
+!       print*,'RD=',rconst
+      do k=1,klev
+         do i=1,ngrid
+c test
+!       print*,'i,k',i,k
+!       print*,'temp(i,k)=',temp(i,k)
+!       print*,'(plev(i,k)-plev(i,k+1))=',plev(i,k),plev(i,k+1)
+            zz=(plev(i,k)+plev(i,k+1))*gravity/(rconst*temp(i,k))
+            kstar(i,k)=0.125*(kmy(i,k+1)+kmy(i,k))*zz*zz
+     s      /(plev(i,k)-plev(i,k+1))*timestep
+         enddo
+      enddo
+
+      do k=2,klev
+         do i=1,ngrid
+            deltap(i,k)=0.5*(plev(i,k-1)-plev(i,k+1))
+         enddo
+      enddo
+      do i=1,ngrid
+         deltap(i,1)=0.5*(plev(i,1)-plev(i,2))
+         deltap(i,klev+1)=0.5*(plev(i,klev)-plev(i,klev+1))
+         denom(i,klev+1)=deltap(i,klev+1)+kstar(i,klev)
+         alpha(i,klev+1)=deltap(i,klev+1)*q2(i,klev+1)/denom(i,klev+1)
+         beta(i,klev+1)=kstar(i,klev)/denom(i,klev+1)
+      enddo
+
+      do k=klev,2,-1
+         do i=1,ngrid
+            denom(i,k)=deltap(i,k)+(1.-beta(i,k+1))*
+     s      kstar(i,k)+kstar(i,k-1)
+c   correction d'un bug 10 01 2001
+            alpha(i,k)=(q2(i,k)*deltap(i,k)
+     s      +kstar(i,k)*alpha(i,k+1))/denom(i,k)
+            beta(i,k)=kstar(i,k-1)/denom(i,k)
+         enddo
+      enddo
+
+c  Si on recalcule q2(1)
+      if(1.eq.0) then
+      do i=1,ngrid
+         denom(i,1)=deltap(i,1)+(1-beta(i,2))*kstar(i,1)
+         q2(i,1)=(q2(i,1)*deltap(i,1)
+     s      +kstar(i,1)*alpha(i,2))/denom(i,1)
+      enddo
+      endif
+c   sinon, on peut sauter cette boucle...
+
+      do k=2,klev+1
+         do i=1,ngrid
+            q2(i,k)=alpha(i,k)+beta(i,k)*q2(i,k-1)
+         enddo
+      enddo
+
+      return
+      end
+      SUBROUTINE vdif_q2e(timestep,gravity,rconst,ngrid,
+     &   plev,temp,kmy,q2)
+      use dimphy
+      IMPLICIT NONE
+c.......................................................................
+#include "dimensions.h"
+cccc#include "dimphy.h"
+c.......................................................................
+c
+c dt : pas de temps
+
+      real plev(klon,klev+1)
+      real temp(klon,klev)
+      real timestep
+      real gravity,rconst
+      real kstar(klon,klev+1),zz
+      real kmy(klon,klev+1)
+      real q2(klon,klev+1)
+      real deltap(klon,klev+1)
+      real denom(klon,klev+1),alpha(klon,klev+1),beta(klon,klev+1)
+      integer ngrid
+
+      integer i,k
+
+      do k=1,klev
+         do i=1,ngrid
+            zz=(plev(i,k)+plev(i,k+1))*gravity/(rconst*temp(i,k))
+            kstar(i,k)=0.125*(kmy(i,k+1)+kmy(i,k))*zz*zz
+     s      /(plev(i,k)-plev(i,k+1))*timestep
+         enddo
+      enddo
+
+      do k=2,klev
+         do i=1,ngrid
+            deltap(i,k)=0.5*(plev(i,k-1)-plev(i,k+1))
+         enddo
+      enddo
+      do i=1,ngrid
+         deltap(i,1)=0.5*(plev(i,1)-plev(i,2))
+         deltap(i,klev+1)=0.5*(plev(i,klev)-plev(i,klev+1))
+      enddo
+
+      do k=klev,2,-1
+         do i=1,ngrid
+            q2(i,k)=q2(i,k)+
+     s      ( kstar(i,k)*(q2(i,k+1)-q2(i,k))
+     s       -kstar(i,k-1)*(q2(i,k)-q2(i,k-1)) )
+     s      /deltap(i,k)
+         enddo
+      enddo
+
+      do i=1,ngrid
+         q2(i,1)=q2(i,1)+
+     s   ( kstar(i,1)*(q2(i,2)-q2(i,1))
+     s                                      )
+     s   /deltap(i,1)
+         q2(i,klev+1)=q2(i,klev+1)+
+     s   ( 
+     s    -kstar(i,klev)*(q2(i,klev+1)-q2(i,klev)) )
+     s   /deltap(i,klev+1)
+      enddo
+
+      return
+      end