source: LMDZ5/trunk/libf/dyn3d/top_bound.F @ 2039

Last change on this file since 2039 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 5.9 KB
RevLine 
[1793]1!
2! $Id: top_bound.F 1907 2013-11-26 13:10:46Z fhourdin $
3!
4      SUBROUTINE top_bound(vcov,ucov,teta,masse,dt)
[999]5      IMPLICIT NONE
6c
7#include "dimensions.h"
8#include "paramet.h"
9#include "comconst.h"
[1279]10#include "comvert.h"
11#include "comgeom2.h"
[999]12
13
14c ..  DISSIPATION LINEAIRE A HAUT NIVEAU, RUN MESO,
15C     F. LOTT DEC. 2006
16c                                 (  10/12/06  )
17
18c=======================================================================
19c
20c   Auteur:  F. LOTT 
21c   -------
22c
23c   Objet:
24c   ------
25c
26c   Dissipation linéaire (ex top_bound de la physique)
27c
28c=======================================================================
29
[1793]30! top_bound sponge layer model:
31! Quenching is modeled as: A(t)=Am+A0*exp(-lambda*t)
32! where Am is the zonal average of the field (or zero), and lambda the inverse
33! of the characteristic quenching/relaxation time scale
34! Thus, assuming Am to be time-independent, field at time t+dt is given by:
35! A(t+dt)=A(t)-(A(t)-Am)*(1-exp(-lambda*t))
36! Moreover lambda can be a function of model level (see below), and relaxation
37! can be toward the average zonal field or just zero (see below).
38
39! NB: top_bound sponge is only called from leapfrog if ok_strato=.true.
40
41! sponge parameters: (loaded/set in conf_gcm.F ; stored in comconst.h)
42!    iflag_top_bound=0 for no sponge
43!    iflag_top_bound=1 for sponge over 4 topmost layers
44!    iflag_top_bound=2 for sponge from top to ~1% of top layer pressure
45!    mode_top_bound=0: no relaxation
46!    mode_top_bound=1: u and v relax towards 0
47!    mode_top_bound=2: u and v relax towards their zonal mean
48!    mode_top_bound=3: u,v and pot. temp. relax towards their zonal mean
49!    tau_top_bound : inverse of charactericstic relaxation time scale at
50!                       the topmost layer (Hz)
51
52
[999]53#include "comdissipn.h"
[1793]54#include "iniprint.h"
[999]55
56c   Arguments:
57c   ----------
58
[1793]59      real,intent(inout) :: ucov(iip1,jjp1,llm) ! covariant zonal wind
60      real,intent(inout) :: vcov(iip1,jjm,llm) ! covariant meridional wind
61      real,intent(inout) :: teta(iip1,jjp1,llm) ! potential temperature
62      real,intent(in) :: masse(iip1,jjp1,llm) ! mass of atmosphere
63      real,intent(in) :: dt ! time step (s) of sponge model
[999]64
65c   Local:
66c   ------
67
[1279]68      REAL massebx(iip1,jjp1,llm),masseby(iip1,jjm,llm),zm
[999]69      REAL uzon(jjp1,llm),vzon(jjm,llm),tzon(jjp1,llm)
70     
[1279]71      integer i
[1793]72      REAL,SAVE :: rdamp(llm) ! quenching coefficient
73      real,save :: lambda(llm) ! inverse or quenching time scale (Hz)
[999]74
[1279]75      LOGICAL,SAVE :: first=.true.
76
[999]77      INTEGER j,l
78     
[1279]79      if (iflag_top_bound.eq.0) return
80
81      if (first) then
82         if (iflag_top_bound.eq.1) then
[1793]83! sponge quenching over the topmost 4 atmospheric layers
84             lambda(:)=0.
85             lambda(llm)=tau_top_bound
86             lambda(llm-1)=tau_top_bound/2.
87             lambda(llm-2)=tau_top_bound/4.
88             lambda(llm-3)=tau_top_bound/8.
[1279]89         else if (iflag_top_bound.eq.2) then
[1793]90! sponge quenching over topmost layers down to pressures which are
91! higher than 100 times the topmost layer pressure
92             lambda(:)=tau_top_bound
[1279]93     s       *max(presnivs(llm)/presnivs(:)-0.01,0.)
94         endif
[1793]95
96! quenching coefficient rdamp(:)
97!         rdamp(:)=dt*lambda(:) ! Explicit Euler approx.
98         rdamp(:)=1.-exp(-lambda(:)*dt)
99
100         write(lunout,*)'TOP_BOUND mode',mode_top_bound
101         write(lunout,*)'Sponge layer coefficients'
102         write(lunout,*)'p (Pa)  z(km)  tau(s)   1./tau (Hz)'
103         do l=1,llm
104           if (rdamp(l).ne.0.) then
105             write(lunout,'(6(1pe12.4,1x))')
106     &        presnivs(l),log(preff/presnivs(l))*scaleheight,
107     &           1./lambda(l),lambda(l)
108           endif
109         enddo
[1279]110         first=.false.
[1793]111      endif ! of if (first)
[1279]112
113      CALL massbar(masse,massebx,masseby)
114
[1793]115      ! compute zonal average of vcov and u
116      if (mode_top_bound.ge.2) then
117       do l=1,llm
[999]118        do j=1,jjm
119          vzon(j,l)=0.
[1279]120          zm=0.
[999]121          do i=1,iim
[1793]122! NB: we can work using vcov zonal mean rather than v since the
123! cv coefficient (which relates the two) only varies with latitudes
[1279]124            vzon(j,l)=vzon(j,l)+vcov(i,j,l)*masseby(i,j,l)
125            zm=zm+masseby(i,j,l)
[999]126          enddo
[1279]127          vzon(j,l)=vzon(j,l)/zm
[999]128        enddo
[1793]129       enddo
[999]130
[1793]131       do l=1,llm
132        do j=2,jjm ! excluding poles
[999]133          uzon(j,l)=0.
[1279]134          zm=0.
135          do i=1,iim
136            uzon(j,l)=uzon(j,l)+massebx(i,j,l)*ucov(i,j,l)/cu(i,j)
137            zm=zm+massebx(i,j,l)
138          enddo
139          uzon(j,l)=uzon(j,l)/zm
140        enddo
[1793]141       enddo
142      else ! ucov and vcov will relax towards 0
143        vzon(:,:)=0.
144        uzon(:,:)=0.
145      endif ! of if (mode_top_bound.ge.2)
[1279]146
[1793]147      ! compute zonal average of potential temperature, if necessary
148      if (mode_top_bound.ge.3) then
149       do l=1,llm
150        do j=2,jjm ! excluding poles
[1279]151          zm=0.
[999]152          tzon(j,l)=0.
153          do i=1,iim
[1279]154            tzon(j,l)=tzon(j,l)+teta(i,j,l)*masse(i,j,l)
155            zm=zm+masse(i,j,l)
[999]156          enddo
[1279]157          tzon(j,l)=tzon(j,l)/zm
[999]158        enddo
[1793]159       enddo
160      endif ! of if (mode_top_bound.ge.3)
[999]161
[1793]162      if (mode_top_bound.ge.1) then
163       ! Apply sponge quenching on vcov:
164       do l=1,llm
165        do i=1,iip1
166          do j=1,jjm
167            vcov(i,j,l)=vcov(i,j,l)
168     &                  -rdamp(l)*(vcov(i,j,l)-vzon(j,l))
169          enddo
170        enddo
171       enddo
[999]172
[1793]173       ! Apply sponge quenching on ucov:
174       do l=1,llm
[999]175        do i=1,iip1
[1793]176          do j=2,jjm ! excluding poles
177            ucov(i,j,l)=ucov(i,j,l)
178     &                  -rdamp(l)*(ucov(i,j,l)-cu(i,j)*uzon(j,l))
[999]179          enddo
180        enddo
[1793]181       enddo
182      endif ! of if (mode_top_bound.ge.1)
[999]183
[1793]184      if (mode_top_bound.ge.3) then
185       ! Apply sponge quenching on teta:
186       do l=1,llm
187        do i=1,iip1
188          do j=2,jjm ! excluding poles
189            teta(i,j,l)=teta(i,j,l)
190     &                  -rdamp(l)*(teta(i,j,l)-tzon(j,l))
191          enddo
192        enddo
193       enddo
194      endif ! of if (mode_top_bound.ge.3)
195   
[999]196      END
Note: See TracBrowser for help on using the repository browser.