source: LMDZ5/branches/LMDZ5_SPLA/libf/dyn3dmem/top_bound_loc.F @ 2508

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