source: LMDZ5/trunk/libf/dyn3dmem/top_bound_loc.F

Last change on this file 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
File size: 7.4 KB
Line 
1!
2! $Id: $
3!
4      SUBROUTINE top_bound_loc(vcov,ucov,teta,masse,dt)
5      USE parallel_lmdz
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,jjb_u:jje_u,llm) ! covariant zonal wind
63      real,intent(inout) :: vcov(iip1,jjb_v:jje_v,llm) ! covariant meridional wind
64      real,intent(inout) :: teta(iip1,jjb_u:jje_u,llm) ! potential temperature
65      real,intent(in) :: masse(iip1,jjb_u:jje_u,llm) ! mass of atmosphere
66      real,intent(in) :: dt ! time step (s) of sponge model
67
68!      REAL dv(iip1,jjb_v:jje_v,llm),du(iip1,jjb_u:jje_u,llm)
69!      REAL dh(iip1,jjb_u:jje_u,llm)
70
71c   Local:
72c   ------
73      REAL massebx(iip1,jjb_u:jje_u,llm),masseby(iip1,jjb_v:jje_v,llm)
74      REAL zm
75      REAL uzon(jjb_u:jje_u,llm),vzon(jjb_v:jje_v,llm)
76      REAL tzon(jjb_u:jje_u,llm)
77     
78      integer i
79      REAL,SAVE :: rdamp(llm)
80      real,save :: lambda(llm) ! inverse or quenching time scale (Hz)
81      LOGICAL,SAVE :: first=.true.
82      INTEGER j,l,jjb,jje
83
84
85      if (iflag_top_bound == 0) return
86
87      if (first) then
88c$OMP BARRIER
89c$OMP MASTER
90         if (iflag_top_bound == 1) then
91! sponge quenching over the topmost 4 atmospheric layers
92             lambda(:)=0.
93             lambda(llm)=tau_top_bound
94             lambda(llm-1)=tau_top_bound/2.
95             lambda(llm-2)=tau_top_bound/4.
96             lambda(llm-3)=tau_top_bound/8.
97         else if (iflag_top_bound == 2) then
98! sponge quenching over topmost layers down to pressures which are
99! higher than 100 times the topmost layer pressure
100             lambda(:)=tau_top_bound
101     s       *max(presnivs(llm)/presnivs(:)-0.01,0.)
102         endif
103
104! quenching coefficient rdamp(:)
105!         rdamp(:)=dt*lambda(:) ! Explicit Euler approx.
106         rdamp(:)=1.-exp(-lambda(:)*dt)
107
108         write(lunout,*)'TOP_BOUND mode',mode_top_bound
109         write(lunout,*)'Sponge layer coefficients'
110         write(lunout,*)'p (Pa)  z(km)  tau(s)   1./tau (Hz)'
111         do l=1,llm
112           if (rdamp(l).ne.0.) then
113             write(lunout,'(6(1pe12.4,1x))')
114     &        presnivs(l),log(preff/presnivs(l))*scaleheight,
115     &           1./lambda(l),lambda(l)
116           endif
117         enddo
118         first=.false.
119c$OMP END MASTER
120c$OMP BARRIER
121      endif ! of if (first)
122
123
124      CALL massbar_loc(masse,massebx,masseby)
125
126      ! compute zonal average of vcov (or set it to zero)
127      if (mode_top_bound.ge.2) then
128       jjb=jj_begin
129       jje=jj_end
130       IF (pole_sud) jje=jj_end-1
131c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
132       do l=1,llm
133        do j=jjb,jje
134          zm=0.
135          vzon(j,l)=0
136          do i=1,iim
137! NB: we can work using vcov zonal mean rather than v since the
138! cv coefficient (which relates the two) only varies with latitudes
139            vzon(j,l)=vzon(j,l)+vcov(i,j,l)*masseby(i,j,l)
140            zm=zm+masseby(i,j,l)
141          enddo
142          vzon(j,l)=vzon(j,l)/zm
143        enddo
144       enddo
145c$OMP END DO NOWAIT   
146      else
147c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
148       do l=1,llm
149         vzon(:,l)=0.
150       enddo
151c$OMP END DO NOWAIT
152      endif ! of if (mode_top_bound.ge.2)
153
154      ! compute zonal average of u (or set it to zero)
155      if (mode_top_bound.ge.2) then
156       jjb=jj_begin
157       jje=jj_end
158       IF (pole_nord) jjb=jj_begin+1
159       IF (pole_sud)  jje=jj_end-1
160c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
161       do l=1,llm
162        do j=jjb,jje
163          uzon(j,l)=0.
164          zm=0.
165          do i=1,iim
166            uzon(j,l)=uzon(j,l)+massebx(i,j,l)*ucov(i,j,l)/cu(i,j)
167            zm=zm+massebx(i,j,l)
168          enddo
169          uzon(j,l)=uzon(j,l)/zm
170        enddo
171       enddo
172c$OMP END DO NOWAIT
173      else
174c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
175       do l=1,llm
176         uzon(:,l)=0.
177       enddo
178c$OMP END DO NOWAIT
179      endif ! of if (mode_top_bound.ge.2)
180
181      ! compute zonal average of potential temperature, if necessary
182      if (mode_top_bound.ge.3) then
183       jjb=jj_begin
184       jje=jj_end
185       IF (pole_nord) jjb=jj_begin+1
186       IF (pole_sud)  jje=jj_end-1
187c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
188       do l=1,llm
189        do j=jjb,jje
190          zm=0.
191          tzon(j,l)=0.
192          do i=1,iim
193            tzon(j,l)=tzon(j,l)+teta(i,j,l)*masse(i,j,l)
194            zm=zm+masse(i,j,l)
195          enddo
196          tzon(j,l)=tzon(j,l)/zm
197        enddo
198       enddo
199c$OMP END DO NOWAIT
200      endif ! of if (mode_top_bound.ge.3)
201
202      if (mode_top_bound.ge.1) then
203       ! Apply sponge quenching on vcov:
204       jjb=jj_begin
205       jje=jj_end
206       IF (pole_sud) jje=jj_end-1
207
208c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
209       do l=1,llm
210        do j=jjb,jje
211          do i=1,iip1
212            vcov(i,j,l)=vcov(i,j,l)
213     &                  -rdamp(l)*(vcov(i,j,l)-vzon(j,l))
214          enddo
215        enddo
216       enddo
217c$OMP END DO NOWAIT
218
219       ! Apply sponge quenching on ucov:
220       jjb=jj_begin
221       jje=jj_end
222       IF (pole_nord) jjb=jj_begin+1
223       IF (pole_sud)  jje=jj_end-1
224
225c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
226       do l=1,llm
227        do j=jjb,jje
228          do i=1,iip1
229            ucov(i,j,l)=ucov(i,j,l)
230     &                  -rdamp(l)*(ucov(i,j,l)-cu(i,j)*uzon(j,l))
231          enddo
232       enddo
233       enddo
234c$OMP END DO NOWAIT
235      endif ! of if (mode_top_bound.ge.1)
236
237      if (mode_top_bound.ge.3) then   
238       ! Apply sponge quenching on teta:
239       jjb=jj_begin
240       jje=jj_end
241       IF (pole_nord) jjb=jj_begin+1
242       IF (pole_sud)  jje=jj_end-1
243
244c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
245       do l=1,llm
246        do j=jjb,jje
247          do i=1,iip1
248            teta(i,j,l)=teta(i,j,l)
249     &                  -rdamp(l)*(teta(i,j,l)-tzon(j,l))
250          enddo
251       enddo
252       enddo
253c$OMP END DO NOWAIT
254      endif ! of if (mode_top_bond.ge.3)
255
256      END
Note: See TracBrowser for help on using the repository browser.