Ignore:
Timestamp:
Jul 18, 2013, 10:20:28 AM (11 years ago)
Author:
Ehouarn Millour
Message:

Version testing basee sur la r1794

Testing release based on r1794

Location:
LMDZ5/branches/testing
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/branches/testing

  • LMDZ5/branches/testing/libf/dyn3dmem/top_bound_loc.F

    r1669 r1795  
    1       SUBROUTINE top_bound_loc( vcov,ucov,teta,masse, du,dv,dh )
     1!
     2! $Id: $
     3!
     4      SUBROUTINE top_bound_loc(vcov,ucov,teta,masse,dt)
    25      USE parallel
    36      IMPLICIT NONE
    … …  
    2528c
    2629c=======================================================================
    27 c-----------------------------------------------------------------------
    28 c   Declarations:
    29 c   -------------
     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
    3053
    3154#include "comdissipn.h"
     55#include "iniprint.h"
    3256
    3357c   Arguments:
    3458c   ----------
    3559
    36       REAL ucov(iip1,jjb_u:jje_u,llm),vcov(iip1,jjb_v:jje_v,llm)
    37       REAL teta(iip1,jjb_u:jje_u,llm)
    38       REAL masse(iip1,jjb_u:jje_u,llm)
    39       REAL dv(iip1,jjb_v:jje_v,llm),du(iip1,jjb_u:jje_u,llm)
    40       REAL dh(iip1,jjb_u:jje_u,llm)
     60      real,intent(inout) :: ucov(iip1,jjb_u:jje_u,llm) ! covariant zonal wind
     61      real,intent(inout) :: vcov(iip1,jjb_v:jje_v,llm) ! covariant meridional wind
     62      real,intent(inout) :: teta(iip1,jjb_u:jje_u,llm) ! potential temperature
     63      real,intent(in) :: masse(iip1,jjb_u:jje_u,llm) ! mass of atmosphere
     64      real,intent(in) :: dt ! time step (s) of sponge model
     65
     66!      REAL dv(iip1,jjb_v:jje_v,llm),du(iip1,jjb_u:jje_u,llm)
     67!      REAL dh(iip1,jjb_u:jje_u,llm)
    4168
    4269c   Local:
    … …  
    4774      REAL tzon(jjb_u:jje_u,llm)
    4875     
    49       INTEGER NDAMP
    50       PARAMETER (NDAMP=4)
    5176      integer i
    5277      REAL,SAVE :: rdamp(llm)
    53 !     &   (/(0., i =1,llm-NDAMP),0.125E-5,.25E-5,.5E-5,1.E-5/)
     78      real,save :: lambda(llm) ! inverse or quenching time scale (Hz)
    5479      LOGICAL,SAVE :: first=.true.
    5580      INTEGER j,l,jjb,jje
    … …  
    5782
    5883      if (iflag_top_bound == 0) return
     84
    5985      if (first) then
    6086c$OMP BARRIER
    6187c$OMP MASTER
    6288         if (iflag_top_bound == 1) then
    63 ! couche eponge dans les 4 dernieres couches du modele
    64              rdamp(:)=0.
    65              rdamp(llm)=tau_top_bound
    66              rdamp(llm-1)=tau_top_bound/2.
    67              rdamp(llm-2)=tau_top_bound/4.
    68              rdamp(llm-3)=tau_top_bound/8.
     89! sponge quenching over the topmost 4 atmospheric layers
     90             lambda(:)=0.
     91             lambda(llm)=tau_top_bound
     92             lambda(llm-1)=tau_top_bound/2.
     93             lambda(llm-2)=tau_top_bound/4.
     94             lambda(llm-3)=tau_top_bound/8.
    6995         else if (iflag_top_bound == 2) then
    70 ! couce eponge dans toutes les couches de pression plus faible que
    71 ! 100 fois la pression de la derniere couche
    72              rdamp(:)=tau_top_bound
     96! sponge quenching over topmost layers down to pressures which are
     97! higher than 100 times the topmost layer pressure
     98             lambda(:)=tau_top_bound
    7399     s       *max(presnivs(llm)/presnivs(:)-0.01,0.)
    74100         endif
     101
     102! quenching coefficient rdamp(:)
     103!         rdamp(:)=dt*lambda(:) ! Explicit Euler approx.
     104         rdamp(:)=1.-exp(-lambda(:)*dt)
     105
     106         write(lunout,*)'TOP_BOUND mode',mode_top_bound
     107         write(lunout,*)'Sponge layer coefficients'
     108         write(lunout,*)'p (Pa)  z(km)  tau(s)   1./tau (Hz)'
     109         do l=1,llm
     110           if (rdamp(l).ne.0.) then
     111             write(lunout,'(6(1pe12.4,1x))')
     112     &        presnivs(l),log(preff/presnivs(l))*scaleheight,
     113     &           1./lambda(l),lambda(l)
     114           endif
     115         enddo
    75116         first=.false.
    76          print*,'TOP_BOUND rdamp=',rdamp
    77117c$OMP END MASTER
    78118c$OMP BARRIER
    79       endif
     119      endif ! of if (first)
    80120
    81121
    82122      CALL massbar_loc(masse,massebx,masseby)
    83 C  CALCUL DES CHAMPS EN MOYENNE ZONALE:
    84 
    85       jjb=jj_begin
    86       jje=jj_end
    87       IF (pole_sud) jje=jj_end-1
    88 
    89 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
    90       do l=1,llm
     123
     124      ! compute zonal average of vcov (or set it to zero)
     125      if (mode_top_bound.ge.2) then
     126       jjb=jj_begin
     127       jje=jj_end
     128       IF (pole_sud) jje=jj_end-1
     129c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
     130       do l=1,llm
    91131        do j=jjb,jje
    92132          zm=0.
    93133          vzon(j,l)=0
    94134          do i=1,iim
    95 ! Rm: on peut travailler directement avec la moyenne zonale de vcov
    96 ! plutot qu'avec celle de v car le coefficient cv qui relie les deux
    97 ! ne varie qu'en latitude
     135! NB: we can work using vcov zonal mean rather than v since the
     136! cv coefficient (which relates the two) only varies with latitudes
    98137            vzon(j,l)=vzon(j,l)+vcov(i,j,l)*masseby(i,j,l)
    99138            zm=zm+masseby(i,j,l)
    … …  
    101140          vzon(j,l)=vzon(j,l)/zm
    102141        enddo
    103       enddo
     142       enddo
    104143c$OMP END DO NOWAIT   
    105 
    106 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
    107       do l=1,llm
    108         do j=jjb,jje
    109           do i=1,iip1
    110             dv(i,j,l)=dv(i,j,l)-rdamp(l)*(vcov(i,j,l)-vzon(j,l))
    111           enddo
    112         enddo
    113       enddo
    114 c$OMP END DO NOWAIT
    115 
    116       jjb=jj_begin
    117       jje=jj_end
    118       IF (pole_nord) jjb=jj_begin+1
    119       IF (pole_sud)  jje=jj_end-1
    120 
    121 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
    122       do l=1,llm
     144      else
     145c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
     146       do l=1,llm
     147         vzon(:,l)=0.
     148       enddo
     149c$OMP END DO NOWAIT
     150      endif ! of if (mode_top_bound.ge.2)
     151
     152      ! compute zonal average of u (or set it to zero)
     153      if (mode_top_bound.ge.2) then
     154       jjb=jj_begin
     155       jje=jj_end
     156       IF (pole_nord) jjb=jj_begin+1
     157       IF (pole_sud)  jje=jj_end-1
     158c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
     159       do l=1,llm
    123160        do j=jjb,jje
    124161          uzon(j,l)=0.
    … …  
    130167          uzon(j,l)=uzon(j,l)/zm
    131168        enddo
    132       enddo
    133 c$OMP END DO NOWAIT
    134 
     169       enddo
     170c$OMP END DO NOWAIT
     171      else
     172c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
     173       do l=1,llm
     174         uzon(:,l)=0.
     175       enddo
     176c$OMP END DO NOWAIT
     177      endif ! of if (mode_top_bound.ge.2)
     178
     179      ! compute zonal average of potential temperature, if necessary
     180      if (mode_top_bound.ge.3) then
     181       jjb=jj_begin
     182       jje=jj_end
     183       IF (pole_nord) jjb=jj_begin+1
     184       IF (pole_sud)  jje=jj_end-1
    135185c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
    136       do l=1,llm
     186       do l=1,llm
    137187        do j=jjb,jje
    138188          zm=0.
    … …  
    144194          tzon(j,l)=tzon(j,l)/zm
    145195        enddo
    146       enddo
    147 c$OMP END DO NOWAIT
    148 
    149 C   AMORTISSEMENTS LINEAIRES:
    150 
    151 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
    152       do l=1,llm
     196       enddo
     197c$OMP END DO NOWAIT
     198      endif ! of if (mode_top_bound.ge.3)
     199
     200      if (mode_top_bound.ge.1) then
     201       ! Apply sponge quenching on vcov:
     202       jjb=jj_begin
     203       jje=jj_end
     204       IF (pole_sud) jje=jj_end-1
     205
     206c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
     207       do l=1,llm
    153208        do j=jjb,jje
    154209          do i=1,iip1
    155             du(i,j,l)=du(i,j,l)
    156      s               -rdamp(l)*(ucov(i,j,l)-cu(i,j)*uzon(j,l))
    157             dh(i,j,l)=dh(i,j,l)-rdamp(l)*(teta(i,j,l)-tzon(j,l))
    158           enddo
    159        enddo
    160       enddo
    161 c$OMP END DO NOWAIT
    162      
    163 
    164       RETURN
     210            vcov(i,j,l)=vcov(i,j,l)
     211     &                  -rdamp(l)*(vcov(i,j,l)-vzon(j,l))
     212          enddo
     213        enddo
     214       enddo
     215c$OMP END DO NOWAIT
     216
     217       ! Apply sponge quenching on ucov:
     218       jjb=jj_begin
     219       jje=jj_end
     220       IF (pole_nord) jjb=jj_begin+1
     221       IF (pole_sud)  jje=jj_end-1
     222
     223c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
     224       do l=1,llm
     225        do j=jjb,jje
     226          do i=1,iip1
     227            ucov(i,j,l)=ucov(i,j,l)
     228     &                  -rdamp(l)*(ucov(i,j,l)-cu(i,j)*uzon(j,l))
     229          enddo
     230       enddo
     231       enddo
     232c$OMP END DO NOWAIT
     233      endif ! of if (mode_top_bound.ge.1)
     234
     235      if (mode_top_bound.ge.3) then   
     236       ! Apply sponge quenching on teta:
     237       jjb=jj_begin
     238       jje=jj_end
     239       IF (pole_nord) jjb=jj_begin+1
     240       IF (pole_sud)  jje=jj_end-1
     241
     242c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
     243       do l=1,llm
     244        do j=jjb,jje
     245          do i=1,iip1
     246            teta(i,j,l)=teta(i,j,l)
     247     &                  -rdamp(l)*(teta(i,j,l)-tzon(j,l))
     248          enddo
     249       enddo
     250       enddo
     251c$OMP END DO NOWAIT
     252      endif ! of if (mode_top_bond.ge.3)
     253
    165254      END
Note: See TracChangeset for help on using the changeset viewer.