Changeset 1795 for LMDZ5/branches/testing/libf/dyn3d/top_bound.F

Timestamp:

Jul 18, 2013, 10:20:28 AM (11 years ago)

Author:

Ehouarn Millour

Message:

Version testing basee sur la r1794

---

Testing release based on r1794

Location:

LMDZ5/branches/testing

Files:

• . (modified) (1 prop)
• libf/dyn3d/top_bound.F (modified) (7 diffs)

LMDZ5/branches/testing/libf/dyn3d/top_bound.F

@@ -1 @@
-      SUBROUTINE top_bound( vcov,ucov,teta,masse, du,dv,dh )
+!
+! $Id$
+!
+      SUBROUTINE top_bound(vcov,ucov,teta,masse,dt)
       IMPLICIT NONE
 c
@@ -24 +27 @@
 c
 c=======================================================================
-c-----------------------------------------------------------------------
-c   Declarations:
-c   -------------
 
-! #include "comgeom.h"
+! top_bound sponge layer model:
+! Quenching is modeled as: A(t)=Am+A0*exp(-lambda*t)
+! where Am is the zonal average of the field (or zero), and lambda the inverse
+! of the characteristic quenching/relaxation time scale
+! Thus, assuming Am to be time-independent, field at time t+dt is given by:
+! A(t+dt)=A(t)-(A(t)-Am)*(1-exp(-lambda*t))
+! Moreover lambda can be a function of model level (see below), and relaxation
+! can be toward the average zonal field or just zero (see below).
+
+! NB: top_bound sponge is only called from leapfrog if ok_strato=.true.
+
+! sponge parameters: (loaded/set in conf_gcm.F ; stored in comconst.h)
+!    iflag_top_bound=0 for no sponge
+!    iflag_top_bound=1 for sponge over 4 topmost layers
+!    iflag_top_bound=2 for sponge from top to ~1% of top layer pressure
+!    mode_top_bound=0: no relaxation
+!    mode_top_bound=1: u and v relax towards 0
+!    mode_top_bound=2: u and v relax towards their zonal mean
+!    mode_top_bound=3: u,v and pot. temp. relax towards their zonal mean
+!    tau_top_bound : inverse of charactericstic relaxation time scale at
+!                       the topmost layer (Hz)
+
+
 #include "comdissipn.h"
+#include "iniprint.h"
 
 c   Arguments:
 c   ----------
 
-      REAL ucov(iip1,jjp1,llm),vcov(iip1,jjm,llm),teta(iip1,jjp1,llm)
-      REAL masse(iip1,jjp1,llm)
-      REAL dv(iip1,jjm,llm),du(iip1,jjp1,llm),dh(iip1,jjp1,llm)
+      real,intent(inout) :: ucov(iip1,jjp1,llm) ! covariant zonal wind
+      real,intent(inout) :: vcov(iip1,jjm,llm) ! covariant meridional wind
+      real,intent(inout) :: teta(iip1,jjp1,llm) ! potential temperature
+      real,intent(in) :: masse(iip1,jjp1,llm) ! mass of atmosphere 
+      real,intent(in) :: dt ! time step (s) of sponge model
 
 c   Local:
@@ -44 +69 @@
       REAL uzon(jjp1,llm),vzon(jjm,llm),tzon(jjp1,llm)
       
-      INTEGER NDAMP
-      PARAMETER (NDAMP=4)
       integer i
-      REAL,SAVE :: rdamp(llm)
-!     &   (/(0., i =1,llm-NDAMP),0.125E-5,.25E-5,.5E-5,1.E-5/) 
+      REAL,SAVE :: rdamp(llm) ! quenching coefficient
+      real,save :: lambda(llm) ! inverse or quenching time scale (Hz)
 
       LOGICAL,SAVE :: first=.true.
 
       INTEGER j,l
-
-
-C  CALCUL DES CHAMPS EN MOYENNE ZONALE:
       
       if (iflag_top_bound.eq.0) return
@@ -61 +81 @@
       if (first) then
          if (iflag_top_bound.eq.1) then
-! couche eponge dans les 4 dernieres couches du modele
-             rdamp(:)=0.
-             rdamp(llm)=tau_top_bound
-             rdamp(llm-1)=tau_top_bound/2.
-             rdamp(llm-2)=tau_top_bound/4.
-             rdamp(llm-3)=tau_top_bound/8.
+! sponge quenching over the topmost 4 atmospheric layers
+             lambda(:)=0.
+             lambda(llm)=tau_top_bound
+             lambda(llm-1)=tau_top_bound/2.
+             lambda(llm-2)=tau_top_bound/4.
+             lambda(llm-3)=tau_top_bound/8.
          else if (iflag_top_bound.eq.2) then
-! couce eponge dans toutes les couches de pression plus faible que
-! 100 fois la pression de la derniere couche
-             rdamp(:)=tau_top_bound
+! sponge quenching over topmost layers down to pressures which are
+! higher than 100 times the topmost layer pressure
+             lambda(:)=tau_top_bound
      s       *max(presnivs(llm)/presnivs(:)-0.01,0.)
          endif
+
+! quenching coefficient rdamp(:)
+!         rdamp(:)=dt*lambda(:) ! Explicit Euler approx.
+         rdamp(:)=1.-exp(-lambda(:)*dt)
+
+         write(lunout,*)'TOP_BOUND mode',mode_top_bound
+         write(lunout,*)'Sponge layer coefficients'
+         write(lunout,*)'p (Pa)  z(km)  tau(s)   1./tau (Hz)'
+         do l=1,llm
+           if (rdamp(l).ne.0.) then
+             write(lunout,'(6(1pe12.4,1x))')
+     &        presnivs(l),log(preff/presnivs(l))*scaleheight,
+     &           1./lambda(l),lambda(l)
+           endif
+         enddo
          first=.false.
-         print*,'TOP_BOUND rdamp=',rdamp
-      endif
+      endif ! of if (first)
 
       CALL massbar(masse,massebx,masseby)
 
-      do l=1,llm
+      ! compute zonal average of vcov and u
+      if (mode_top_bound.ge.2) then
+       do l=1,llm
         do j=1,jjm
           vzon(j,l)=0.
           zm=0.
           do i=1,iim
-! Rm: on peut travailler directement avec la moyenne zonale de vcov
-! plutot qu'avec celle de v car le coefficient cv qui relie les deux
-! ne varie qu'en latitude
+! NB: we can work using vcov zonal mean rather than v since the
+! cv coefficient (which relates the two) only varies with latitudes 
             vzon(j,l)=vzon(j,l)+vcov(i,j,l)*masseby(i,j,l)
             zm=zm+masseby(i,j,l)
@@ -92 +127 @@
           vzon(j,l)=vzon(j,l)/zm
         enddo
-      enddo
+       enddo
 
-      do l=1,llm
-        do i=1,iip1
-          do j=1,jjm
-            dv(i,j,l)=dv(i,j,l)-rdamp(l)*(vcov(i,j,l)-vzon(j,l))
-          enddo
-        enddo
-      enddo
-
-      do l=1,llm
-        do j=2,jjm
+       do l=1,llm
+        do j=2,jjm ! excluding poles
           uzon(j,l)=0.
           zm=0.
@@ -112 +139 @@
           uzon(j,l)=uzon(j,l)/zm
         enddo
-      enddo
+       enddo
+      else ! ucov and vcov will relax towards 0
+        vzon(:,:)=0.
+        uzon(:,:)=0.
+      endif ! of if (mode_top_bound.ge.2)
 
-      do l=1,llm
-        do j=2,jjm
+      ! compute zonal average of potential temperature, if necessary
+      if (mode_top_bound.ge.3) then
+       do l=1,llm
+        do j=2,jjm ! excluding poles
           zm=0.
           tzon(j,l)=0.
@@ -124 +157 @@
           tzon(j,l)=tzon(j,l)/zm
         enddo
-      enddo
+       enddo
+      endif ! of if (mode_top_bound.ge.3)
 
-C   AMORTISSEMENTS LINEAIRES:
-
-      do l=1,llm
+      if (mode_top_bound.ge.1) then
+       ! Apply sponge quenching on vcov:
+       do l=1,llm
         do i=1,iip1
-          do j=2,jjm
-            du(i,j,l)=du(i,j,l)
-     s               -rdamp(l)*(ucov(i,j,l)-cu(i,j)*uzon(j,l))
-            dh(i,j,l)=dh(i,j,l)-rdamp(l)*(teta(i,j,l)-tzon(j,l))
+          do j=1,jjm
+            vcov(i,j,l)=vcov(i,j,l)
+     &                  -rdamp(l)*(vcov(i,j,l)-vzon(j,l))
           enddo
         enddo
-      enddo
-      
+       enddo
 
-      RETURN
+       ! Apply sponge quenching on ucov:
+       do l=1,llm
+        do i=1,iip1
+          do j=2,jjm ! excluding poles
+            ucov(i,j,l)=ucov(i,j,l)
+     &                  -rdamp(l)*(ucov(i,j,l)-cu(i,j)*uzon(j,l))
+          enddo
+        enddo
+       enddo
+      endif ! of if (mode_top_bound.ge.1)
+
+      if (mode_top_bound.ge.3) then
+       ! Apply sponge quenching on teta:
+       do l=1,llm
+        do i=1,iip1
+          do j=2,jjm ! excluding poles
+            teta(i,j,l)=teta(i,j,l)
+     &                  -rdamp(l)*(teta(i,j,l)-tzon(j,l))
+          enddo
+        enddo
+       enddo
+      endif ! of if (mode_top_bound.ge.3)
+    
       END