source: LMDZ5/trunk/libf/dyn3d/top_bound.F @ 5179

Last change on this file since 5179 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: 6.0 KB
RevLine 
[1793]1!
2! $Id: top_bound.F 2600 2016-07-23 05:45:38Z abarral $
3!
4      SUBROUTINE top_bound(vcov,ucov,teta,masse,dt)
[2597]5     
6      USE comconst_mod, ONLY: iflag_top_bound, mode_top_bound,
7     &                        tau_top_bound
[2600]8      USE comvert_mod, ONLY: presnivs, preff, scaleheight
[2597]9     
[999]10      IMPLICIT NONE
11c
[2597]12      include "dimensions.h"
13      include "paramet.h"
14      include "comgeom2.h"
[999]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
[999]56#include "comdissipn.h"
[1793]57#include "iniprint.h"
[999]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
[999]67
68c   Local:
69c   ------
70
[1279]71      REAL massebx(iip1,jjp1,llm),masseby(iip1,jjm,llm),zm
[999]72      REAL uzon(jjp1,llm),vzon(jjm,llm),tzon(jjp1,llm)
73     
[1279]74      integer i
[1793]75      REAL,SAVE :: rdamp(llm) ! quenching coefficient
76      real,save :: lambda(llm) ! inverse or quenching time scale (Hz)
[999]77
[1279]78      LOGICAL,SAVE :: first=.true.
79
[999]80      INTEGER j,l
81     
[1279]82      if (iflag_top_bound.eq.0) return
83
84      if (first) then
85         if (iflag_top_bound.eq.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.eq.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.
[1793]114      endif ! of if (first)
[1279]115
116      CALL massbar(masse,massebx,masseby)
117
[1793]118      ! compute zonal average of vcov and u
119      if (mode_top_bound.ge.2) then
120       do l=1,llm
[999]121        do j=1,jjm
122          vzon(j,l)=0.
[1279]123          zm=0.
[999]124          do i=1,iim
[1793]125! NB: we can work using vcov zonal mean rather than v since the
126! cv coefficient (which relates the two) only varies with latitudes
[1279]127            vzon(j,l)=vzon(j,l)+vcov(i,j,l)*masseby(i,j,l)
128            zm=zm+masseby(i,j,l)
[999]129          enddo
[1279]130          vzon(j,l)=vzon(j,l)/zm
[999]131        enddo
[1793]132       enddo
[999]133
[1793]134       do l=1,llm
135        do j=2,jjm ! excluding poles
[999]136          uzon(j,l)=0.
[1279]137          zm=0.
138          do i=1,iim
139            uzon(j,l)=uzon(j,l)+massebx(i,j,l)*ucov(i,j,l)/cu(i,j)
140            zm=zm+massebx(i,j,l)
141          enddo
142          uzon(j,l)=uzon(j,l)/zm
143        enddo
[1793]144       enddo
145      else ! ucov and vcov will relax towards 0
146        vzon(:,:)=0.
147        uzon(:,:)=0.
148      endif ! of if (mode_top_bound.ge.2)
[1279]149
[1793]150      ! compute zonal average of potential temperature, if necessary
151      if (mode_top_bound.ge.3) then
152       do l=1,llm
153        do j=2,jjm ! excluding poles
[1279]154          zm=0.
[999]155          tzon(j,l)=0.
156          do i=1,iim
[1279]157            tzon(j,l)=tzon(j,l)+teta(i,j,l)*masse(i,j,l)
158            zm=zm+masse(i,j,l)
[999]159          enddo
[1279]160          tzon(j,l)=tzon(j,l)/zm
[999]161        enddo
[1793]162       enddo
163      endif ! of if (mode_top_bound.ge.3)
[999]164
[1793]165      if (mode_top_bound.ge.1) then
166       ! Apply sponge quenching on vcov:
167       do l=1,llm
168        do i=1,iip1
169          do j=1,jjm
170            vcov(i,j,l)=vcov(i,j,l)
171     &                  -rdamp(l)*(vcov(i,j,l)-vzon(j,l))
172          enddo
173        enddo
174       enddo
[999]175
[1793]176       ! Apply sponge quenching on ucov:
177       do l=1,llm
[999]178        do i=1,iip1
[1793]179          do j=2,jjm ! excluding poles
180            ucov(i,j,l)=ucov(i,j,l)
181     &                  -rdamp(l)*(ucov(i,j,l)-cu(i,j)*uzon(j,l))
[999]182          enddo
183        enddo
[1793]184       enddo
185      endif ! of if (mode_top_bound.ge.1)
[999]186
[1793]187      if (mode_top_bound.ge.3) then
188       ! Apply sponge quenching on teta:
189       do l=1,llm
190        do i=1,iip1
191          do j=2,jjm ! excluding poles
192            teta(i,j,l)=teta(i,j,l)
193     &                  -rdamp(l)*(teta(i,j,l)-tzon(j,l))
194          enddo
195        enddo
196       enddo
197      endif ! of if (mode_top_bound.ge.3)
198   
[999]199      END
Note: See TracBrowser for help on using the repository browser.