source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/dyn3dpar/top_bound_p.F @ 245

Last change on this file since 245 was 222, checked in by ymipsl, 10 years ago

Creating temporary dynamico/lmdz/saturn branche

YM
File size: 7.3 KB
Line 
1!
2! $Id: top_bound_p.F 1793 2013-07-18 07:13:18Z emillour $
3!
4      SUBROUTINE top_bound_p(vcov,ucov,teta,masse,dt,ducov)
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*dt)
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,jjp1,llm) ! covariant zonal wind
61      real,intent(inout) :: vcov(iip1,jjm,llm) ! covariant meridional wind
62      real,intent(inout) :: teta(iip1,jjp1,llm) ! potential temperature
63      real,intent(in) :: masse(iip1,jjp1,llm) ! mass of atmosphere
64      real,intent(in) :: dt ! time step (s) of sponge model
65      real,intent(out) :: ducov(iip1,jjp1,llm) ! increment on ucov due to sponge
66
67c   Local:
68c   ------
69      REAL massebx(iip1,jjp1,llm),masseby(iip1,jjm,llm),zm
70      REAL uzon(jjp1,llm),vzon(jjm,llm),tzon(jjp1,llm)
71     
72      integer i 
73      REAL,SAVE :: rdamp(llm) ! quenching coefficient
74      real,save :: lambda(llm) ! inverse or quenching time scale (Hz)
75      LOGICAL,SAVE :: first=.true.
76      INTEGER j,l,jjb,jje
77
78
79      if (first) then
80c$OMP BARRIER
81c$OMP MASTER
82         if (iflag_top_bound == 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 == 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.
111c$OMP END MASTER
112c$OMP BARRIER
113      endif ! of if (first)
114
115
116      CALL massbar_p(masse,massebx,masseby)
117
118      ! compute zonal average of vcov (or set it to zero)
119      if (mode_top_bound.ge.2) then
120       jjb=jj_begin
121       jje=jj_end
122       IF (pole_sud) jje=jj_end-1
123c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
124       do l=1,llm
125        do j=jjb,jje
126          zm=0.
127          vzon(j,l)=0
128          do i=1,iim
129! NB: we can work using vcov zonal mean rather than v since the
130! cv coefficient (which relates the two) only varies with latitudes
131            vzon(j,l)=vzon(j,l)+vcov(i,j,l)*masseby(i,j,l)
132            zm=zm+masseby(i,j,l)
133          enddo
134          vzon(j,l)=vzon(j,l)/zm
135        enddo
136       enddo
137c$OMP END DO NOWAIT   
138      else
139c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
140       do l=1,llm
141         vzon(:,l)=0.
142       enddo
143c$OMP END DO NOWAIT
144      endif ! of if (mode_top_bound.ge.2)
145
146      ! compute zonal average of u (or set it to zero)
147      if (mode_top_bound.ge.2) then
148       jjb=jj_begin
149       jje=jj_end
150       IF (pole_nord) jjb=jj_begin+1
151       IF (pole_sud)  jje=jj_end-1
152c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
153       do l=1,llm
154        do j=jjb,jje
155          uzon(j,l)=0.
156          zm=0.
157          do i=1,iim
158            uzon(j,l)=uzon(j,l)+massebx(i,j,l)*ucov(i,j,l)/cu(i,j)
159            zm=zm+massebx(i,j,l)
160          enddo
161          uzon(j,l)=uzon(j,l)/zm
162        enddo
163       enddo
164c$OMP END DO NOWAIT
165      else
166c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
167       do l=1,llm
168         uzon(:,l)=0.
169       enddo
170c$OMP END DO NOWAIT
171      endif ! of if (mode_top_bound.ge.2)
172
173      ! compute zonal average of potential temperature, if necessary
174      if (mode_top_bound.ge.3) then
175       jjb=jj_begin
176       jje=jj_end
177       IF (pole_nord) jjb=jj_begin+1
178       IF (pole_sud)  jje=jj_end-1
179c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
180       do l=1,llm
181        do j=jjb,jje
182          zm=0.
183          tzon(j,l)=0.
184          do i=1,iim
185            tzon(j,l)=tzon(j,l)+teta(i,j,l)*masse(i,j,l)
186            zm=zm+masse(i,j,l)
187          enddo
188          tzon(j,l)=tzon(j,l)/zm
189        enddo
190       enddo
191c$OMP END DO NOWAIT
192      endif ! of if (mode_top_bound.ge.3)
193
194      if (mode_top_bound.ge.1) then
195       ! Apply sponge quenching on vcov:
196       jjb=jj_begin
197       jje=jj_end
198       IF (pole_sud) jje=jj_end-1
199
200c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
201       do l=1,llm
202        do j=jjb,jje
203          do i=1,iip1
204            vcov(i,j,l)=vcov(i,j,l)
205     &                  -rdamp(l)*(vcov(i,j,l)-vzon(j,l))
206          enddo
207        enddo
208       enddo
209c$OMP END DO NOWAIT
210
211       ! Apply sponge quenching on ucov:
212       jjb=jj_begin
213       jje=jj_end
214       IF (pole_nord) jjb=jj_begin+1
215       IF (pole_sud)  jje=jj_end-1
216
217c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
218       do l=1,llm
219        do j=jjb,jje
220          do i=1,iip1
221            ducov(i,j,l)=-rdamp(l)*(ucov(i,j,l)-cu(i,j)*uzon(j,l))
222            ucov(i,j,l)=ucov(i,j,l)
223     &                  +ducov(i,j,l)
224          enddo
225       enddo
226       enddo
227c$OMP END DO NOWAIT
228      endif ! of if (mode_top_bound.ge.1)
229
230      if (mode_top_bound.ge.3) then   
231       ! Apply sponge quenching on teta:
232       jjb=jj_begin
233       jje=jj_end
234       IF (pole_nord) jjb=jj_begin+1
235       IF (pole_sud)  jje=jj_end-1
236
237c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
238       do l=1,llm
239        do j=jjb,jje
240          do i=1,iip1
241            teta(i,j,l)=teta(i,j,l)
242     &                  -rdamp(l)*(teta(i,j,l)-tzon(j,l))
243          enddo
244       enddo
245       enddo
246c$OMP END DO NOWAIT
247      endif ! of if (mode_top_bond.ge.3)
248
249      END
Note: See TracBrowser for help on using the repository browser.