Changeset 1311 for LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/yamada4.F

Timestamp:

Feb 18, 2010, 2:14:02 PM (14 years ago)

Author:

Laurent Fairhead

Message:

Modifications to thermals for TKE transport

---

Modifications aux thermiques pour le transport de la TKE

pbl_surface_mode.F90 : ok_flux_surf=.false. seulement pour klon>1
physiq.F : option iflag_pbl=10 pour transporter la TKE avec les thermiques.
calltherm.F90 : passage de iflag_thermals_ed en argument pour thermcell_plume
thermcell_main.F90 : Appel a plusieurs version de thermcell_plume en option
thermcell_plume.F90 : plusieurs versions dans le meme fichier (temporaire)
thermcell_height.F90 : verrue pour les cas ou les thermiques montent tout

en haut

yamada4 : inclusion de la diffusion verticale en option iflag_pbl=9

+ variables anciennement common, puis save/allocatable, remises en local

File:

• LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/yamada4.F (modified) (8 diffs)

LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/yamada4.F

@@ -1 +1 @@
 !
-! $Id$
+! $Header$
 !
       SUBROUTINE yamada4(ngrid,dt,g,rconst,plev,temp
@@ -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.
-
-      character (len=20) :: modname='yamada4'
-      character (len=80) :: abort_message
-
 c$OMP THREADPRIVATE(firstcall)       
       frif(ri)=0.6588*(ri+0.1776-sqrt(ri*ri-0.3221*ri+0.03156))
@@ -123 +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.
@@ -131 +118 @@
 
 
-      if (.not.(iflag_pbl.ge.6.and.iflag_pbl.le.9)) then
-              abort_message = 'probleme de coherence dans appel a MY'
-              CALL abort_gcm (modname,abort_message,1)
+      if (.not.(iflag_pbl.ge.6.and.iflag_pbl.le.10)) then
+           stop'probleme de coherence dans appel a MY'
       endif
 
@@ -417 +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
+      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
@@ -497 +484 @@
       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