Changeset 1146 for LMDZ4/trunk/libf/phylmd/cv3_routines.F

Timestamp:

Apr 9, 2009, 12:11:35 PM (15 years ago)

Author:

Laurent Fairhead

Message:

Réintegration dans le tronc des modifications issues de la branche LMDZ-dev
comprises entre la révision 1074 et 1145
Validation: une simulation de 1 jour en séquentiel sur PC donne les mêmes
résultats entre la trunk et la dev
LF

Location:

LMDZ4/trunk

Files:

• . (modified) (1 prop)
• libf/phylmd/cv3_routines.F (modified) (9 diffs)

LMDZ4/trunk/libf/phylmd/cv3_routines.F

@@ -1 +1 @@
 !
-! $Header$
+! $Header: /home/cvsroot/LMDZ4/libf/phylmd/cv3_routines.F,v 1.16 2008-11-06 16:29:35 lmdzadmin Exp $
 !
 c
@@ -120 +120 @@
 
 c ori      do 110 k=1,nlp
-      do 110 k=1,nl ! convect3
+! abderr     do 110 k=1,nl ! convect3
+       do 110 k=1,nlp
+      
         do 100 i=1,len
 cdebug          lv(i,k)= lv0-clmcpv*(t(i,k)-t0)
@@ -2256 +2258 @@
       SUBROUTINE cv3_yield(nloc,ncum,nd,na,ntra
      :                    ,icb,inb,delt
-     :                    ,t,rr,t_wake,rr_wake,u,v,tra
+     :                    ,t,rr,t_wake,rr_wake,s_wake,u,v,tra
      :                    ,gz,p,ph,h,hp,lv,cpn,th,th_wake
      :                    ,ep,clw,m,tp,mp,rp,up,vp,trap
@@ -2283 +2285 @@
       real t(nloc,nd), rr(nloc,nd), u(nloc,nd), v(nloc,nd)
       real t_wake(nloc,nd), rr_wake(nloc,nd)
+      real s_wake(nloc)
       real tra(nloc,nd,ntra), sig(nloc,nd)
       real gz(nloc,na), ph(nloc,nd+1), h(nloc,na), hp(nloc,na)
@@ -2327 +2330 @@
       real esum(nloc), fsum(nloc), gsum(nloc), hsum(nloc)
       real th_wake(nloc,nd)
+      real alpha_qpos(nloc)
       real qcond(nloc,nd), nqcond(nloc,nd), wa(nloc,nd)  ! cld
       real siga(nloc,nd), sax(nloc,nd), mac(nloc,nd)      ! cld
@@ -2961 +2965 @@
 c   ***          integrated enthalpy and water tendencies         ***
 c
+c Correction bug le 18-03-09
       do 503 il=1,ncum
       IF (iflag(il) .le. 1) THEN
-      ax=0.1*ment(il,inb(il),inb(il))*(hp(il,inb(il))-h(il,inb(il))
-     : +t(il,inb(il))*(cpv-cpd)
+        if (cvflag_grav) then
+      ax=0.01*grav*ment(il,inb(il),inb(il))*(hp(il,inb(il))
+     : -h(il,inb(il))+t(il,inb(il))*(cpv-cpd)
      : *(rr(il,inb(il))-qent(il,inb(il),inb(il))))
      :  /(cpn(il,inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1)))
@@ -2972 +2978 @@
      :    *(ph(il,inb(il)-1)-ph(il,inb(il))))
 
+      bx=0.01*grav*ment(il,inb(il),inb(il))*(qent(il,inb(il),inb(il))
+     :    -rr(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
+      fr(il,inb(il))=fr(il,inb(il))-bx
+      fr(il,inb(il)-1)=fr(il,inb(il)-1)
+     :   +bx*(ph(il,inb(il))-ph(il,inb(il)+1))
+     :      /(ph(il,inb(il)-1)-ph(il,inb(il)))
+
+      cx=0.01*grav*ment(il,inb(il),inb(il))*(uent(il,inb(il),inb(il))
+     :       -u(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
+      fu(il,inb(il))=fu(il,inb(il))-cx
+      fu(il,inb(il)-1)=fu(il,inb(il)-1)
+     :     +cx*(ph(il,inb(il))-ph(il,inb(il)+1))
+     :        /(ph(il,inb(il)-1)-ph(il,inb(il)))
+
+      dx=0.01*grav*ment(il,inb(il),inb(il))*(vent(il,inb(il),inb(il))
+     :      -v(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
+      fv(il,inb(il))=fv(il,inb(il))-dx
+      fv(il,inb(il)-1)=fv(il,inb(il)-1)
+     :    +dx*(ph(il,inb(il))-ph(il,inb(il)+1))
+     :       /(ph(il,inb(il)-1)-ph(il,inb(il)))
+       else
+       ax=0.1*ment(il,inb(il),inb(il))*(hp(il,inb(il))
+     : -h(il,inb(il))+t(il,inb(il))*(cpv-cpd)
+     : *(rr(il,inb(il))-qent(il,inb(il),inb(il))))
+     :  /(cpn(il,inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1)))
+      ft(il,inb(il))=ft(il,inb(il))-ax
+      ft(il,inb(il)-1)=ft(il,inb(il)-1)+ax*cpn(il,inb(il))
+     :    *(ph(il,inb(il))-ph(il,inb(il)+1))/(cpn(il,inb(il)-1)
+     :    *(ph(il,inb(il)-1)-ph(il,inb(il))))
+
       bx=0.1*ment(il,inb(il),inb(il))*(qent(il,inb(il),inb(il))
      :    -rr(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
@@ -2992 +3028 @@
      :    +dx*(ph(il,inb(il))-ph(il,inb(il)+1))
      :       /(ph(il,inb(il)-1)-ph(il,inb(il)))
+       endif
       ENDIF    !iflag
 503   continue
@@ -3065 +3102 @@
        enddo
       enddo
+
+c
+c   ***   Check that moisture stays positive. If not, scale tendencies
+c        in order to ensure moisture positivity
+      DO il = 1,ncum
+       IF (iflag(il) .le. 1) THEN
+        alpha_qpos(il) = max(1. , -delt*fr(il,1)/
+     :     (s_wake(il)*rr_wake(il,1)+(1.-s_wake(il))*rr(il,1)))
+       ENDIF
+      ENDDO
+      DO i = 2,nl
+       DO il = 1,ncum
+        IF (iflag(il) .le. 1) THEN
+        alpha_qpos(il) = max(alpha_qpos(il) , -delt*fr(il,i)/
+     :     (s_wake(il)*rr_wake(il,i)+(1.-s_wake(il))*rr(il,i)))
+        ENDIF
+       ENDDO
+      ENDDO
+      DO il = 1,ncum
+       IF (iflag(il) .le. 1 .and. alpha_qpos(il) .gt. 1.001) THEN
+        alpha_qpos(il) = alpha_qpos(il)*1.1
+       ENDIF
+      ENDDO
+      DO il = 1,ncum
+       IF (iflag(il) .le. 1) THEN
+        sigd(il) = sigd(il)/alpha_qpos(il)
+        precip(il) = precip(il)/alpha_qpos(il)
+       ENDIF
+      ENDDO
+      DO i = 1,nl
+       DO il = 1,ncum
+        IF (iflag(il) .le. 1) THEN
+         fr(il,i) = fr(il,i)/alpha_qpos(il)
+         ft(il,i) = ft(il,i)/alpha_qpos(il)
+         fqd(il,i) = fqd(il,i)/alpha_qpos(il)
+         ftd(il,i) = ftd(il,i)/alpha_qpos(il)
+         fu(il,i) = fu(il,i)/alpha_qpos(il)
+         fv(il,i) = fv(il,i)/alpha_qpos(il)
+         m(il,i) = m(il,i)/alpha_qpos(il)
+         mp(il,i) = mp(il,i)/alpha_qpos(il)
+         Vprecip(il,i) = Vprecip(il,i)/alpha_qpos(il)
+        ENDIF
+       ENDDO
+      ENDDO
+      DO i = 1,nl
+      DO j = 1,nl
+       DO il = 1,ncum
+        IF (iflag(il) .le. 1) THEN
+         ment(il,i,j) = ment(il,i,j)/alpha_qpos(il)
+        ENDIF
+       ENDDO
+      ENDDO
+      ENDDO
+      DO j = 1,ntra
+      DO i = 1,nl
+       DO il = 1,ncum
+        IF (iflag(il) .le. 1) THEN
+         ftra(il,i,j) = ftra(il,i,j)/alpha_qpos(il)
+        ENDIF
+       ENDDO
+      ENDDO
+      ENDDO
+
 c
 c   ***           reset counter and return           ***