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/dyn3dpar/top_bound_p.F

    r1279 r1795  
    1       SUBROUTINE top_bound_p( vcov,ucov,teta,masse, du,dv,dh )
     1!
     2! $Id$
     3!
     4      SUBROUTINE top_bound_p(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,jjp1,llm),vcov(iip1,jjm,llm),teta(iip1,jjp1,llm)
    37       REAL masse(iip1,jjp1,llm)
    38       REAL dv(iip1,jjm,llm),du(iip1,jjp1,llm),dh(iip1,jjp1,llm)
     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
    3965
    4066c   Local:
    … …  
    4369      REAL uzon(jjp1,llm),vzon(jjm,llm),tzon(jjp1,llm)
    4470     
    45       INTEGER NDAMP
    46       PARAMETER (NDAMP=4)
    4771      integer i
    48       REAL,SAVE :: rdamp(llm)
    49 !     &   (/(0., i =1,llm-NDAMP),0.125E-5,.25E-5,.5E-5,1.E-5/)
     72      REAL,SAVE :: rdamp(llm) ! quenching coefficient
     73      real,save :: lambda(llm) ! inverse or quenching time scale (Hz)
    5074      LOGICAL,SAVE :: first=.true.
    5175      INTEGER j,l,jjb,jje
    … …  
    5377
    5478      if (iflag_top_bound == 0) return
     79
    5580      if (first) then
    5681c$OMP BARRIER
    5782c$OMP MASTER
    5883         if (iflag_top_bound == 1) then
    59 ! couche eponge dans les 4 dernieres couches du modele
    60              rdamp(:)=0.
    61              rdamp(llm)=tau_top_bound
    62              rdamp(llm-1)=tau_top_bound/2.
    63              rdamp(llm-2)=tau_top_bound/4.
    64              rdamp(llm-3)=tau_top_bound/8.
     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.
    6590         else if (iflag_top_bound == 2) then
    66 ! couce eponge dans toutes les couches de pression plus faible que
    67 ! 100 fois la pression de la derniere couche
    68              rdamp(:)=tau_top_bound
     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
    6994     s       *max(presnivs(llm)/presnivs(:)-0.01,0.)
    7095         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
    71111         first=.false.
    72          print*,'TOP_BOUND rdamp=',rdamp
    73112c$OMP END MASTER
    74113c$OMP BARRIER
    75       endif
     114      endif ! of if (first)
    76115
    77116
    78117      CALL massbar_p(masse,massebx,masseby)
    79 C  CALCUL DES CHAMPS EN MOYENNE ZONALE:
    80 
    81       jjb=jj_begin
    82       jje=jj_end
    83       IF (pole_sud) jje=jj_end-1
    84 
    85 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
    86       do l=1,llm
     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
    87126        do j=jjb,jje
    88127          zm=0.
    89128          vzon(j,l)=0
    90129          do i=1,iim
    91 ! Rm: on peut travailler directement avec la moyenne zonale de vcov
    92 ! plutot qu'avec celle de v car le coefficient cv qui relie les deux
    93 ! ne varie qu'en latitude
     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
    94132            vzon(j,l)=vzon(j,l)+vcov(i,j,l)*masseby(i,j,l)
    95133            zm=zm+masseby(i,j,l)
    … …  
    97135          vzon(j,l)=vzon(j,l)/zm
    98136        enddo
    99       enddo
     137       enddo
    100138c$OMP END DO NOWAIT   
    101 
    102 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
    103       do l=1,llm
    104         do j=jjb,jje
    105           do i=1,iip1
    106             dv(i,j,l)=dv(i,j,l)-rdamp(l)*(vcov(i,j,l)-vzon(j,l))
    107           enddo
    108         enddo
    109       enddo
    110 c$OMP END DO NOWAIT
    111 
    112       jjb=jj_begin
    113       jje=jj_end
    114       IF (pole_nord) jjb=jj_begin+1
    115       IF (pole_sud)  jje=jj_end-1
    116 
    117 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
    118       do l=1,llm
     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
    119155        do j=jjb,jje
    120156          uzon(j,l)=0.
    … …  
    126162          uzon(j,l)=uzon(j,l)/zm
    127163        enddo
    128       enddo
    129 c$OMP END DO NOWAIT
    130 
     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
    131180c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
    132       do l=1,llm
     181       do l=1,llm
    133182        do j=jjb,jje
    134183          zm=0.
    … …  
    140189          tzon(j,l)=tzon(j,l)/zm
    141190        enddo
    142       enddo
    143 c$OMP END DO NOWAIT
    144 
    145 C   AMORTISSEMENTS LINEAIRES:
    146 
    147 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
    148       do l=1,llm
     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
    149203        do j=jjb,jje
    150204          do i=1,iip1
    151             du(i,j,l)=du(i,j,l)
    152      s               -rdamp(l)*(ucov(i,j,l)-cu(i,j)*uzon(j,l))
    153             dh(i,j,l)=dh(i,j,l)-rdamp(l)*(teta(i,j,l)-tzon(j,l))
    154           enddo
    155        enddo
    156       enddo
    157 c$OMP END DO NOWAIT
    158      
    159 
    160       RETURN
     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
    161249      END
Note: See TracChangeset for help on using the changeset viewer.