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

Last change on this file since 2526 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
Line 
1!
2! $Id: top_bound.F 1907 2013-11-26 13:10:46Z oboucher $
3!
4      SUBROUTINE top_bound(vcov,ucov,teta,masse,dt)
5      IMPLICIT NONE
6c
7#include "dimensions.h"
8#include "paramet.h"
9#include "comconst.h"
10#include "comvert.h"
11#include "comgeom2.h"
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
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
53#include "comdissipn.h"
54#include "iniprint.h"
55
56c   Arguments:
57c   ----------
58
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
64
65c   Local:
66c   ------
67
68      REAL massebx(iip1,jjp1,llm),masseby(iip1,jjm,llm),zm
69      REAL uzon(jjp1,llm),vzon(jjm,llm),tzon(jjp1,llm)
70     
71      integer i
72      REAL,SAVE :: rdamp(llm) ! quenching coefficient
73      real,save :: lambda(llm) ! inverse or quenching time scale (Hz)
74
75      LOGICAL,SAVE :: first=.true.
76
77      INTEGER j,l
78     
79      if (iflag_top_bound.eq.0) return
80
81      if (first) then
82         if (iflag_top_bound.eq.1) then
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.
89         else if (iflag_top_bound.eq.2) then
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
93     s       *max(presnivs(llm)/presnivs(:)-0.01,0.)
94         endif
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
110         first=.false.
111      endif ! of if (first)
112
113      CALL massbar(masse,massebx,masseby)
114
115      ! compute zonal average of vcov and u
116      if (mode_top_bound.ge.2) then
117       do l=1,llm
118        do j=1,jjm
119          vzon(j,l)=0.
120          zm=0.
121          do i=1,iim
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
124            vzon(j,l)=vzon(j,l)+vcov(i,j,l)*masseby(i,j,l)
125            zm=zm+masseby(i,j,l)
126          enddo
127          vzon(j,l)=vzon(j,l)/zm
128        enddo
129       enddo
130
131       do l=1,llm
132        do j=2,jjm ! excluding poles
133          uzon(j,l)=0.
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
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)
146
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
151          zm=0.
152          tzon(j,l)=0.
153          do i=1,iim
154            tzon(j,l)=tzon(j,l)+teta(i,j,l)*masse(i,j,l)
155            zm=zm+masse(i,j,l)
156          enddo
157          tzon(j,l)=tzon(j,l)/zm
158        enddo
159       enddo
160      endif ! of if (mode_top_bound.ge.3)
161
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
172
173       ! Apply sponge quenching on ucov:
174       do l=1,llm
175        do i=1,iip1
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))
179          enddo
180        enddo
181       enddo
182      endif ! of if (mode_top_bound.ge.1)
183
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   
196      END
Note: See TracBrowser for help on using the repository browser.