source: LMDZ5/trunk/libf/dyn3dpar/top_bound_p.F @ 3012

Last change on this file since 3012 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: 7.3 KB
RevLine 
[1793]1!
2! $Id: top_bound_p.F 2600 2016-07-23 05:45:38Z fairhead $
3!
4      SUBROUTINE top_bound_p(vcov,ucov,teta,masse,dt)
[1823]5      USE parallel_lmdz
[2597]6      USE comconst_mod, ONLY: iflag_top_bound, mode_top_bound,
7     &                        tau_top_bound
[2600]8      USE comvert_mod, ONLY: presnivs, preff, scaleheight
9     
[1000]10      IMPLICIT NONE
11c
12#include "dimensions.h"
13#include "paramet.h"
[1279]14#include "comgeom2.h"
[1000]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
[1793]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
[2597]44! sponge parameters: (loaded/set in conf_gcm.F ; stored in comconst_mod)
[1793]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
[1000]56#include "comdissipn.h"
[1793]57#include "iniprint.h"
[1000]58
59c   Arguments:
60c   ----------
61
[1793]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
[1000]67
68c   Local:
69c   ------
[1279]70      REAL massebx(iip1,jjp1,llm),masseby(iip1,jjm,llm),zm
[1000]71      REAL uzon(jjp1,llm),vzon(jjm,llm),tzon(jjp1,llm)
72     
[2600]73      integer i
[1793]74      REAL,SAVE :: rdamp(llm) ! quenching coefficient
75      real,save :: lambda(llm) ! inverse or quenching time scale (Hz)
[1279]76      LOGICAL,SAVE :: first=.true.
[1000]77      INTEGER j,l,jjb,jje
78
79
[1279]80      if (iflag_top_bound == 0) return
[1793]81
[1279]82      if (first) then
83c$OMP BARRIER
84c$OMP MASTER
85         if (iflag_top_bound == 1) then
[1793]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.
[1279]92         else if (iflag_top_bound == 2) then
[1793]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
[1279]96     s       *max(presnivs(llm)/presnivs(:)-0.01,0.)
97         endif
[1793]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
[1279]113         first=.false.
114c$OMP END MASTER
115c$OMP BARRIER
[1793]116      endif ! of if (first)
[1279]117
118
119      CALL massbar_p(masse,massebx,masseby)
[1000]120
[1793]121      ! compute zonal average of vcov (or set it to zero)
122      if (mode_top_bound.ge.2) then
123       jjb=jj_begin
124       jje=jj_end
125       IF (pole_sud) jje=jj_end-1
[1000]126c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
[1793]127       do l=1,llm
[1000]128        do j=jjb,jje
[1279]129          zm=0.
130          vzon(j,l)=0
[1000]131          do i=1,iim
[1793]132! NB: we can work using vcov zonal mean rather than v since the
133! cv coefficient (which relates the two) only varies with latitudes
[1279]134            vzon(j,l)=vzon(j,l)+vcov(i,j,l)*masseby(i,j,l)
135            zm=zm+masseby(i,j,l)
[1000]136          enddo
[1279]137          vzon(j,l)=vzon(j,l)/zm
[1000]138        enddo
[1793]139       enddo
[1000]140c$OMP END DO NOWAIT   
[1793]141      else
[1000]142c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
[1793]143       do l=1,llm
144         vzon(:,l)=0.
145       enddo
[1000]146c$OMP END DO NOWAIT
[1793]147      endif ! of if (mode_top_bound.ge.2)
[1000]148
[1793]149      ! compute zonal average of u (or set it to zero)
150      if (mode_top_bound.ge.2) then
151       jjb=jj_begin
152       jje=jj_end
153       IF (pole_nord) jjb=jj_begin+1
154       IF (pole_sud)  jje=jj_end-1
[1000]155c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
[1793]156       do l=1,llm
[1000]157        do j=jjb,jje
158          uzon(j,l)=0.
[1279]159          zm=0.
160          do i=1,iim
161            uzon(j,l)=uzon(j,l)+massebx(i,j,l)*ucov(i,j,l)/cu(i,j)
162            zm=zm+massebx(i,j,l)
163          enddo
164          uzon(j,l)=uzon(j,l)/zm
165        enddo
[1793]166       enddo
[1279]167c$OMP END DO NOWAIT
[1793]168      else
169c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
170       do l=1,llm
171         uzon(:,l)=0.
172       enddo
173c$OMP END DO NOWAIT
174      endif ! of if (mode_top_bound.ge.2)
[1279]175
[1793]176      ! compute zonal average of potential temperature, if necessary
177      if (mode_top_bound.ge.3) then
178       jjb=jj_begin
179       jje=jj_end
180       IF (pole_nord) jjb=jj_begin+1
181       IF (pole_sud)  jje=jj_end-1
[1279]182c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
[1793]183       do l=1,llm
[1279]184        do j=jjb,jje
185          zm=0.
[1000]186          tzon(j,l)=0.
187          do i=1,iim
[1279]188            tzon(j,l)=tzon(j,l)+teta(i,j,l)*masse(i,j,l)
189            zm=zm+masse(i,j,l)
[1000]190          enddo
[1279]191          tzon(j,l)=tzon(j,l)/zm
[1000]192        enddo
[1793]193       enddo
[1000]194c$OMP END DO NOWAIT
[1793]195      endif ! of if (mode_top_bound.ge.3)
[1000]196
[1793]197      if (mode_top_bound.ge.1) then
198       ! Apply sponge quenching on vcov:
199       jjb=jj_begin
200       jje=jj_end
201       IF (pole_sud) jje=jj_end-1
[1000]202
203c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
[1793]204       do l=1,llm
[1000]205        do j=jjb,jje
206          do i=1,iip1
[1793]207            vcov(i,j,l)=vcov(i,j,l)
208     &                  -rdamp(l)*(vcov(i,j,l)-vzon(j,l))
[1000]209          enddo
[1793]210        enddo
[1000]211       enddo
212c$OMP END DO NOWAIT
213
[1793]214       ! Apply sponge quenching on ucov:
215       jjb=jj_begin
216       jje=jj_end
217       IF (pole_nord) jjb=jj_begin+1
218       IF (pole_sud)  jje=jj_end-1
219
220c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
221       do l=1,llm
222        do j=jjb,jje
223          do i=1,iip1
224            ucov(i,j,l)=ucov(i,j,l)
225     &                  -rdamp(l)*(ucov(i,j,l)-cu(i,j)*uzon(j,l))
226          enddo
227       enddo
228       enddo
229c$OMP END DO NOWAIT
230      endif ! of if (mode_top_bound.ge.1)
231
232      if (mode_top_bound.ge.3) then   
233       ! Apply sponge quenching on teta:
234       jjb=jj_begin
235       jje=jj_end
236       IF (pole_nord) jjb=jj_begin+1
237       IF (pole_sud)  jje=jj_end-1
238
239c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
240       do l=1,llm
241        do j=jjb,jje
242          do i=1,iip1
243            teta(i,j,l)=teta(i,j,l)
244     &                  -rdamp(l)*(teta(i,j,l)-tzon(j,l))
245          enddo
246       enddo
247       enddo
248c$OMP END DO NOWAIT
249      endif ! of if (mode_top_bond.ge.3)
250
[1000]251      END
Note: See TracBrowser for help on using the repository browser.