Ignore:
Timestamp:
Feb 24, 2014, 4:05:47 PM (11 years ago)
Author:
Ehouarn Millour
Message:

Add updating pressure, mass and Exner function (ie: all variables which depend on surface pressure) after adding physics tendencies (which include a surface pressure tendency).
Note that this change induces slight changes in GCM results with respect to previous svn version of the code, even if surface pressure tendency is zero (because of recomputation of polar values as an average over polar points on the dynamics grid).
EM

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/dyn3dpar/caldyn_p.F

    r1907 r1987  
    11!
    2 ! $Header$
    3 !
    4 c
    5 c
     2! $Id$
     3!
    64#undef DEBUG_IO
    7 c#define DEBUG_IO
     5!#define DEBUG_IO
    86
    97      SUBROUTINE caldyn_p
    … …  
    1513      IMPLICIT NONE
    1614
    17 c=======================================================================
    18 c
    19 c  Auteur :  P. Le Van
    20 c
    21 c   Objet:
    22 c   ------
    23 c
    24 c   Calcul des tendances dynamiques.
    25 c
    26 c Modif 04/93 F.Forget
    27 c=======================================================================
    28 
    29 c-----------------------------------------------------------------------
    30 c   0. Declarations:
    31 c   ----------------
     15!=======================================================================
     16!
     17!  Auteur :  P. Le Van
     18!
     19!   Objet:
     20!   ------
     21!
     22!   Calcul des tendances dynamiques.
     23!
     24! Modif 04/93 F.Forget
     25!=======================================================================
     26
     27!-----------------------------------------------------------------------
     28!   0. Declarations:
     29!   ----------------
    3230
    3331#include "dimensions.h"
    … …  
    3735#include "comgeom.h"
    3836
    39 c   Arguments:
    40 c   ----------
    41 
    42       LOGICAL conser
    43 
    44       INTEGER itau
    45       REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm),teta(ip1jmp1,llm)
    46       REAL ps(ip1jmp1),phis(ip1jmp1)
    47       REAL pk(iip1,jjp1,llm),pkf(ip1jmp1,llm)
     37!   Arguments:
     38!   ----------
     39
     40      LOGICAL,INTENT(IN) :: conser ! triggers printing some diagnostics
     41      INTEGER,INTENT(IN) :: itau ! time step index
     42      REAL,INTENT(IN) :: vcov(ip1jm,llm) ! covariant meridional wind
     43      REAL,INTENT(IN) :: ucov(ip1jmp1,llm) ! covariant zonal wind
     44      REAL,INTENT(IN) :: teta(ip1jmp1,llm) ! potential temperature
     45      REAL,INTENT(IN) :: ps(ip1jmp1) ! surface pressure
     46      REAL,INTENT(IN) :: phis(ip1jmp1) ! geopotential at the surface
     47      REAL,INTENT(IN) :: pk(ip1jmp1,llm) ! Exner at mid-layer
     48      REAL,INTENT(IN) :: pkf(ip1jmp1,llm) ! filtered Exner
     49      REAL,INTENT(IN) :: phi(ip1jmp1,llm) ! geopotential
     50      REAL,INTENT(OUT) :: masse(ip1jmp1,llm) ! air mass
     51      REAL,INTENT(OUT) :: dv(ip1jm,llm) ! tendency on vcov
     52      REAL,INTENT(OUT) :: du(ip1jmp1,llm) ! tendency on ucov
     53      REAL,INTENT(OUT) :: dteta(ip1jmp1,llm) ! tenddency on teta
     54      REAL,INTENT(OUT) :: dp(ip1jmp1) ! tendency on ps
     55      REAL,INTENT(OUT) :: w(ip1jmp1,llm) ! vertical velocity
     56      REAL,INTENT(OUT) :: pbaru(ip1jmp1,llm) ! mass flux in the zonal direction
     57      REAL,INTENT(OUT) :: pbarv(ip1jm,llm) ! mass flux in the meridional direction
     58      REAL,INTENT(IN) :: time ! current time
     59
     60!   Local:
     61!   ------
     62
    4863      REAL,SAVE :: vcont(ip1jm,llm),ucont(ip1jmp1,llm)
    49       REAL phi(ip1jmp1,llm),masse(ip1jmp1,llm)
    50       REAL dv(ip1jm,llm),du(ip1jmp1,llm)
    51       REAL dteta(ip1jmp1,llm),dp(ip1jmp1)
    52       REAL w(ip1jmp1,llm)
    53       REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
    54       REAL time
    55 
    56 c   Local:
    57 c   ------
    58 
    5964      REAL,SAVE :: ang(ip1jmp1,llm)
    6065      REAL,SAVE :: p(ip1jmp1,llmp1)
    … …  
    6873      INTEGER   ij,l,ijb,ije,ierr
    6974
    70 c-----------------------------------------------------------------------
    71 c   Calcul des tendances dynamiques:
    72 c   --------------------------------
     75!-----------------------------------------------------------------------
     76!   Compute dynamical tendencies:
     77!--------------------------------
     78
     79      ! compute contravariant winds ucont() and vcont
    7380      CALL covcont_p  ( llm    , ucov    , vcov , ucont, vcont        )
     81      ! compute pressure p()
    7482      CALL pression_p ( ip1jmp1, ap      , bp   ,  ps  , p            )
    75 cym      CALL psextbar (   ps   , psexbarxy                          )
    76 c$OMP BARRIER
     83!ym      CALL psextbar (   ps   , psexbarxy                          )
     84!$OMP BARRIER
     85      ! compute mass in each atmospheric mesh: masse()
    7786      CALL massdair_p (    p   , masse                                )
     87      ! compute X and Y-averages of mass, massebx() and masseby()
    7888      CALL massbar_p  (   masse, massebx , masseby                    )
     89      ! compute XY-average of mass, massebxy()
    7990      call massbarxy_p(   masse, massebxy                             )
     91      ! compute mass fluxes pbaru() and pbarv()
    8092      CALL flumass_p  ( massebx, masseby , vcont, ucont ,pbaru, pbarv )
     93      ! compute dteta() , horizontal converging flux of theta
    8194      CALL dteta1_p   (   teta , pbaru   , pbarv, dteta               )
     95      ! compute convm(), horizontal converging flux of mass
    8296      CALL convmas1_p  (   pbaru, pbarv   , convm                      )
    83 c$OMP BARRIER     
     97!$OMP BARRIER     
    8498      CALL convmas2_p  (   convm                      )
    85 c$OMP BARRIER
     99!$OMP BARRIER
    86100#ifdef DEBUG_IO
    87 c$OMP BARRIER
    88 c$OMP MASTER
     101!$OMP BARRIER
     102!$OMP MASTER
    89103      call WriteField_p('ucont',reshape(ucont,(/iip1,jmp1,llm/)))
    90104      call WriteField_p('vcont',reshape(vcont,(/iip1,jjm,llm/)))
    … …  
    98112      call WriteField_p('dteta',reshape(dteta,(/iip1,jmp1,llm/)))
    99113      call WriteField_p('convm',reshape(convm,(/iip1,jmp1,llm/)))
    100 c$OMP END MASTER
    101 c$OMP BARRIER
     114!$OMP END MASTER
     115!$OMP BARRIER
    102116#endif     
    103117
    104 c$OMP BARRIER
    105 c$OMP MASTER
     118!$OMP BARRIER
     119!$OMP MASTER
    106120      ijb=ij_begin
    107121      ije=ij_end
    108            
     122      ! compute pressure variation due to mass convergence
    109123      DO ij =ijb, ije
    110124         dp( ij ) = convm( ij,1 ) / airesurg( ij )
    111125      ENDDO
    112 c$OMP END MASTER
    113 c$OMP BARRIER
    114 c$OMP FLUSH
     126!$OMP END MASTER
     127!$OMP BARRIER
     128!$OMP FLUSH
     129     
     130      ! compute vertical velocity w()
    115131      CALL vitvert_p ( convm  , w                                  )
     132      ! compute potential vorticity vorpot()
    116133      CALL tourpot_p ( vcov   , ucov  , massebxy  , vorpot         )
     134      ! compute rotation induced du() and dv()
    117135      CALL dudv1_p   ( vorpot , pbaru , pbarv     , du     , dv    )
    118136
    119137#ifdef DEBUG_IO     
    120 c$OMP BARRIER
    121 c$OMP MASTER
     138!$OMP BARRIER
     139!$OMP MASTER
    122140      call WriteField_p('w',reshape(w,(/iip1,jmp1,llm/)))
    123141      call WriteField_p('vorpot',reshape(vorpot,(/iip1,jjm,llm/)))
    124142      call WriteField_p('du',reshape(du,(/iip1,jmp1,llm/)))
    125143      call WriteField_p('dv',reshape(dv,(/iip1,jjm,llm/)))
    126 c$OMP END MASTER
    127 c$OMP BARRIER
     144!$OMP END MASTER
     145!$OMP BARRIER
    128146#endif     
     147     
     148      ! compute kinetic energy ecin()
    129149      CALL enercin_p ( vcov   , ucov  , vcont     , ucont  , ecin  )
     150      ! compute Bernouilli function bern()
    130151      CALL bernoui_p ( ip1jmp1, llm   , phi       , ecin   , bern  )
     152      ! compute and add du() and dv() contributions from Bernouilli and pressure
    131153      CALL dudv2_p   ( teta   , pkf   , bern      , du     , dv    )
    132154
    133155#ifdef DEBUG_IO
    134 c$OMP BARRIER
    135 c$OMP MASTER
     156!$OMP BARRIER
     157!$OMP MASTER
    136158      call WriteField_p('ecin',reshape(ecin,(/iip1,jmp1,llm/)))
    137159      call WriteField_p('bern',reshape(bern,(/iip1,jmp1,llm/)))
    … …  
    139161      call WriteField_p('dv',reshape(dv,(/iip1,jjm,llm/)))
    140162      call WriteField_p('pkf',reshape(pkf,(/iip1,jmp1,llm/)))
    141 c$OMP END MASTER
    142 c$OMP BARRIER
     163!$OMP END MASTER
     164!$OMP BARRIER
    143165#endif
    144166     
    … …  
    149171      if (pole_sud) ije=ij_end
    150172
    151 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)       
     173!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)       
    152174      DO l=1,llm
    153175         DO ij=ijb,ije
    … …  
    155177        ENDDO
    156178      ENDDO
    157 c$OMP END DO
    158 
     179!$OMP END DO
     180
     181      ! compute vertical advection contributions to du(), dv() and dteta()
    159182      CALL advect_new_p(ang,vcov,teta,w,massebx,masseby,du,dv,dteta)
    160183
    161 C  WARNING probleme de peridocite de dv sur les PC/linux. Pb d'arrondi
    162 C          probablement. Observe sur le code compile avec pgf90 3.0-1
     184!  WARNING probleme de peridocite de dv sur les PC/linux. Pb d'arrondi
     185!          probablement. Observe sur le code compile avec pgf90 3.0-1
    163186      ijb=ij_begin
    164187      ije=ij_end
    165188      if (pole_sud) ije=ij_end-iip1
    166189
    167 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     190!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    168191      DO l = 1, llm
    169192         DO ij = ijb, ije, iip1
    170193           IF( dv(ij,l).NE.dv(ij+iim,l) )  THEN
    171 c         PRINT *,'!!!ATTENTION!!! probleme de periodicite sur vcov', 
    172 c    ,   ' dans caldyn'
    173 c         PRINT *,' l,  ij = ', l, ij, ij+iim,dv(ij+iim,l),dv(ij,l)
     194!         PRINT *,'!!!ATTENTION!!! probleme de periodicite sur vcov', 
     195!    ,   ' dans caldyn'
     196!         PRINT *,' l,  ij = ', l, ij, ij+iim,dv(ij+iim,l),dv(ij,l)
    174197          dv(ij+iim,l) = dv(ij,l)
    175198          endif
    176199         enddo
    177200      enddo
    178 c$OMP END DO NOWAIT     
    179 c-----------------------------------------------------------------------
    180 c   Sorties eventuelles des variables de controle:
    181 c   ----------------------------------------------
     201!$OMP END DO NOWAIT     
     202!-----------------------------------------------------------------------
     203!   Output some control variables:
     204!---------------------------------
    182205
    183206      IF( conser )  THEN
    184 c ym ---> exige communication collective ( aussi dans advect)
     207! ym ---> exige communication collective ( aussi dans advect)
    185208        CALL sortvarc
    186      $ ( itau,ucov,teta,ps,masse,pk,phis,vorpot,phi,bern,dp,time,vcov )
     209     & ( itau,ucov,teta,ps,masse,pk,phis,vorpot,phi,bern,dp,time,vcov )
    187210
    188211      ENDIF
    189212
    190       RETURN
    191213      END
Note: See TracChangeset for help on using the changeset viewer.