source: LMDZ6/branches/LMDZ-INCA-Dyn/libf/dyn3d/top_bound.F @ 4976

Last change on this file since 4976 was 2600, checked in by Ehouarn Millour, 8 years ago

Cleanup in the dynamics: turn comvert.h into module comvert_mod.F90
EM

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