Changeset 1301 for trunk/LMDZ.VENUS/libf/phyvenus/physiq.F

Timestamp:

Jun 26, 2014, 12:25:31 PM (10 years ago)

Author:

slebonnois

Message:

SL: many bug corrections in phyvenus, some cleaning, and a new ksi matrix format for Venus IR

File:

• trunk/LMDZ.VENUS/libf/phyvenus/physiq.F (modified) (15 diffs)

trunk/LMDZ.VENUS/libf/phyvenus/physiq.F

@@ -7 +7 @@
      .            paprs,pplay,ppk,pphi,pphis,presnivs,
      .            u,v,t,qx,
-     .            omega,
+     .            flxmw,
      .            d_u, d_v, d_t, d_qx, d_ps)
 
@@ -46 +46 @@
 c qx------input-R-mass mixing ratio traceurs (kg/kg) 
 c d_t_dyn-input-R-tendance dynamique pour "t" (K/s)
-c omega---input-R-vitesse verticale en Pa/s
+c flxmw---input-R-flux de masse vertical en kg/s
 c
 c d_u-----output-R-tendance physique de "u" (m/s/s)
@@ -125 +125 @@
       REAL d_t_dyn(klon,klev)
 
-      REAL omega(klon,klev)
+      REAL flxmw(klon,klev)
 
       REAL d_u(klon,klev)
@@ -146 +146 @@
 c Variables propres a la physique
 c
-      REAL,save,allocatable :: rlev(:,:) ! altitude a chaque niveau (interface inferieure de la couche) 
       INTEGER,save :: itap        ! compteur pour la physique
       REAL delp(klon,klev)        ! epaisseur d'une couche
+      REAL omega(klon,klev)       ! vitesse verticale en Pa/s
+
       
       INTEGER igwd,idx(klon),itest(klon)
@@ -242 +243 @@
 c
       REAL zphi(klon,klev)
-c
+      REAL zzlev(klon,klev+1),zzlay(klon,klev),z1,z2
+
 c Variables du changement
 c
@@ -362 +364 @@
 c========================
       IF (debut) THEN
-         allocate(rlev(klon,klevp1))
          allocate(source(klon,nqmax))
 
@@ -621 +622 @@
       ENDIF
 c====================================================================
+
+c Calcule de vitesse verticale a partir de flux de masse verticale
+      DO k = 1, klev
+       DO i = 1, klon
+        omega(i,k) = RG*flxmw(i,k) / airephy(i)
+       END DO
+      END DO
+
 c
 c Ajouter le geopotentiel du sol:
@@ -629 +638 @@
       ENDDO
       ENDDO
+
+c   calcul du geopotentiel aux niveaux intercouches
+c   ponderation des altitudes au niveau des couches en dp/p
+
+      DO k=1,klev
+         DO i=1,klon
+            zzlay(i,k)=zphi(i,k)/RG
+         ENDDO
+      ENDDO
+      DO i=1,klon
+         zzlev(i,1)=pphis(i)/RG
+      ENDDO
+      DO k=2,klev
+         DO i=1,klon
+            z1=(pplay(i,k-1)+paprs(i,k))/(pplay(i,k-1)-paprs(i,k))
+            z2=(paprs(i,k)  +pplay(i,k))/(paprs(i,k)  -pplay(i,k))
+            zzlev(i,k)=(z1*zzlay(i,k-1)+z2*zzlay(i,k))/(z1+z2)
+         ENDDO
+      ENDDO
+      DO i=1,klon
+         zzlev(i,klev+1)=zzlay(i,klev)+(zzlay(i,klev)-zzlev(i,klev))
+      ENDDO
+
 c====================================================================
 c
@@ -729 +761 @@
       fder = dlw
 
-c     print*,"radsol avant clmain=",radsol(klon/2)
-c     print*,"solsw avant clmain=",solsw(klon/2)
-c     print*,"sollw avant clmain=",sollw(klon/2)
-
 ! ADAPTATION GCM POUR CP(T)
 
@@ -749 +777 @@
      s            ycoefh,yu1,yv1) 
 
-c     print*,"radsol apres clmain=",radsol(klon/2)
-c     print*,"solsw apres clmain=",solsw(klon/2)
-c     print*,"sollw apres clmain=",sollw(klon/2)
-
 CXXX Incrementation des flux
       DO i = 1, klon
@@ -770 +794 @@
       ENDDO
       ENDDO
-c
-c        print*,"d_t_vdf1=",d_t_vdf(1,:)*dtime
-c        print*,"d_t_vdf2=",d_t_vdf(klon/2,:)*dtime
-c        print*,"d_t_vdf3=",d_t_vdf(klon,:)*dtime
-c        print*,"d_u_vdf=",d_u_vdf(klon/2,:)*dtime
-c        print*,"d_v_vdf=",d_v_vdf(klon/2,:)*dtime
 
 C TRACEURS
@@ -810 +828 @@
 c Incrementer la temperature du sol
 c
-c     print*,'Tsol avant clmain:',ftsol(1)
       DO i = 1, klon
          ftsol(i) = ftsol(i) + d_ts(i)
       ENDDO
-c     print*,'Tsol apres clmain:',ftsol(1)
 
 c Calculer la derive du flux infrarouge
@@ -878 +894 @@
          v_seri(:,:) = v_seri(:,:) + d_v_ajs(:,:)
          d_v_ajs(:,:)= d_v_ajs(:,:)/dtime          ! (m/s)/s
-      if (iflag_trac.eq.1) then
+
+         if (iflag_trac.eq.1) then
            tr_seri(:,:,:) = tr_seri(:,:,:) + d_tr_ajs(:,:,:)
            d_tr_ajs(:,:,:)= d_tr_ajs(:,:,:)/dtime  ! /s
-      endif
-
-c        print*,"d_t_ajs1=",d_t_ajs(1,:)*dtime
-c        print*,"d_t_ajs2=",d_t_ajs(klon/2,:)*dtime
-c        print*,"d_t_ajs3=",d_t_ajs(klon,:)*dtime
-c        print*,"d_u_ajs=",d_u_ajs(klon/2,:)*dtime
-c        print*,"d_v_ajs=",d_v_ajs(klon/2,:)*dtime
+         endif
 
       endif
@@ -919 +930 @@
 c     PRINT*,'dtimerad,dtime,radpas',dtimerad,dtime,radpas
             
-c     print*,"radsol avant radlwsw=",radsol(klon/2)
-c     print*,"solsw avant radlwsw=",solsw(klon/2)
-c     print*,"sollw avant radlwsw=",sollw(klon/2)
-c     print*,"avant radlwsw"
-
       CALL radlwsw
-     e            (dist, rmu0, fract, 
+     e            (dist, rmu0, fract, zzlev,
      e             paprs, pplay,ftsol, t_seri,
      s             heat,cool,radsol,
@@ -932 +938 @@
      s             lwnet, swnet)
 
-c     print*,"apres radlwsw"
-c     print*,"radsol apres radlwsw=",radsol(klon/2)
-c     print*,"solsw apres radlwsw=",solsw(klon/2)
-c     print*,"sollw apres radlwsw=",sollw(klon/2)
       itaprad = 0
       DO k = 1, klev
       DO i = 1, klon
-         dtrad(i,k) = heat(i,k)-cool(i,k)
-         dtrad(i,k) = dtrad(i,k)/RDAY  !K/s
-      ENDDO
-      ENDDO
-c        print*,"dtrad1=",dtrad(1,:)*dtime
-c        print*,"dtrad2=",dtrad(klon/2,:)*dtime
-c        print*,"dtrad3=",dtrad(klon,:)*dtime
-      
+         dtrad(i,k) = heat(i,k)-cool(i,k)  !K/s
+      ENDDO
+      ENDDO
+
       ENDIF
       itaprad = itaprad + 1