source: LMDZ5/branches/AI-cosp/libf/dyn3dpar/top_bound_p.F @ 3717

Last change on this file since 3717 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: 7.1 KB
Line 
1!
2! $Id: top_bound_p.F 1907 2013-11-26 13:10:46Z adurocher $
3!
4      SUBROUTINE top_bound_p(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,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
66c   Local:
67c   ------
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      LOGICAL,SAVE :: first=.true.
75      INTEGER j,l,jjb,jje
76
77
78      if (iflag_top_bound == 0) return
79
80      if (first) then
81c$OMP BARRIER
82c$OMP MASTER
83         if (iflag_top_bound == 1) then
84! sponge quenching over the topmost 4 atmospheric layers
85             lambda(:)=0.
86             lambda(llm)=tau_top_bound
87             lambda(llm-1)=tau_top_bound/2.
88             lambda(llm-2)=tau_top_bound/4.
89             lambda(llm-3)=tau_top_bound/8.
90         else if (iflag_top_bound == 2) then
91! sponge quenching over topmost layers down to pressures which are
92! higher than 100 times the topmost layer pressure
93             lambda(:)=tau_top_bound
94     s       *max(presnivs(llm)/presnivs(:)-0.01,0.)
95         endif
96
97! quenching coefficient rdamp(:)
98!         rdamp(:)=dt*lambda(:) ! Explicit Euler approx.
99         rdamp(:)=1.-exp(-lambda(:)*dt)
100
101         write(lunout,*)'TOP_BOUND mode',mode_top_bound
102         write(lunout,*)'Sponge layer coefficients'
103         write(lunout,*)'p (Pa)  z(km)  tau(s)   1./tau (Hz)'
104         do l=1,llm
105           if (rdamp(l).ne.0.) then
106             write(lunout,'(6(1pe12.4,1x))')
107     &        presnivs(l),log(preff/presnivs(l))*scaleheight,
108     &           1./lambda(l),lambda(l)
109           endif
110         enddo
111         first=.false.
112c$OMP END MASTER
113c$OMP BARRIER
114      endif ! of if (first)
115
116
117      CALL massbar_p(masse,massebx,masseby)
118
119      ! compute zonal average of vcov (or set it to zero)
120      if (mode_top_bound.ge.2) then
121       jjb=jj_begin
122       jje=jj_end
123       IF (pole_sud) jje=jj_end-1
124c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
125       do l=1,llm
126        do j=jjb,jje
127          zm=0.
128          vzon(j,l)=0
129          do i=1,iim
130! NB: we can work using vcov zonal mean rather than v since the
131! cv coefficient (which relates the two) only varies with latitudes
132            vzon(j,l)=vzon(j,l)+vcov(i,j,l)*masseby(i,j,l)
133            zm=zm+masseby(i,j,l)
134          enddo
135          vzon(j,l)=vzon(j,l)/zm
136        enddo
137       enddo
138c$OMP END DO NOWAIT   
139      else
140c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
141       do l=1,llm
142         vzon(:,l)=0.
143       enddo
144c$OMP END DO NOWAIT
145      endif ! of if (mode_top_bound.ge.2)
146
147      ! compute zonal average of u (or set it to zero)
148      if (mode_top_bound.ge.2) then
149       jjb=jj_begin
150       jje=jj_end
151       IF (pole_nord) jjb=jj_begin+1
152       IF (pole_sud)  jje=jj_end-1
153c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
154       do l=1,llm
155        do j=jjb,jje
156          uzon(j,l)=0.
157          zm=0.
158          do i=1,iim
159            uzon(j,l)=uzon(j,l)+massebx(i,j,l)*ucov(i,j,l)/cu(i,j)
160            zm=zm+massebx(i,j,l)
161          enddo
162          uzon(j,l)=uzon(j,l)/zm
163        enddo
164       enddo
165c$OMP END DO NOWAIT
166      else
167c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
168       do l=1,llm
169         uzon(:,l)=0.
170       enddo
171c$OMP END DO NOWAIT
172      endif ! of if (mode_top_bound.ge.2)
173
174      ! compute zonal average of potential temperature, if necessary
175      if (mode_top_bound.ge.3) then
176       jjb=jj_begin
177       jje=jj_end
178       IF (pole_nord) jjb=jj_begin+1
179       IF (pole_sud)  jje=jj_end-1
180c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
181       do l=1,llm
182        do j=jjb,jje
183          zm=0.
184          tzon(j,l)=0.
185          do i=1,iim
186            tzon(j,l)=tzon(j,l)+teta(i,j,l)*masse(i,j,l)
187            zm=zm+masse(i,j,l)
188          enddo
189          tzon(j,l)=tzon(j,l)/zm
190        enddo
191       enddo
192c$OMP END DO NOWAIT
193      endif ! of if (mode_top_bound.ge.3)
194
195      if (mode_top_bound.ge.1) then
196       ! Apply sponge quenching on vcov:
197       jjb=jj_begin
198       jje=jj_end
199       IF (pole_sud) jje=jj_end-1
200
201c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
202       do l=1,llm
203        do j=jjb,jje
204          do i=1,iip1
205            vcov(i,j,l)=vcov(i,j,l)
206     &                  -rdamp(l)*(vcov(i,j,l)-vzon(j,l))
207          enddo
208        enddo
209       enddo
210c$OMP END DO NOWAIT
211
212       ! Apply sponge quenching on ucov:
213       jjb=jj_begin
214       jje=jj_end
215       IF (pole_nord) jjb=jj_begin+1
216       IF (pole_sud)  jje=jj_end-1
217
218c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
219       do l=1,llm
220        do j=jjb,jje
221          do i=1,iip1
222            ucov(i,j,l)=ucov(i,j,l)
223     &                  -rdamp(l)*(ucov(i,j,l)-cu(i,j)*uzon(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.