Ignore:
Timestamp:
Jul 31, 2015, 7:22:21 PM (9 years ago)
Author:
dcugnet
Message:
  • Add parallel capability for ce0l.
  • Small bug in grid_noro fixed (smoothed topography was used instead of unsmoothed one for geopotential computation at north pole).
  • Removed average of mass at poles in etat0dyn_netcdf after start_init_dyn => different results in the zoomed grid case.
  • ok_etat0=n and ok_limit=y combination now works fine (if no initial state is needed, but only limit.nc file). This required:
    • to move grid_noro0 and start_init_noro0 subroutines from etat0dyn_netcdf.F90 to limit_netcdf.F90
    • to create init_ssrf_m.F90 file, so that sub-surfaces can be initialized from limit_netcdf.F90 without any etat0*_netcdf routines call).
  • Simplified somehow the corresponding code, in particular: 1) removed obsolete flags "oldice". 2) removed flag "ibar": barycentric interpolation is used everywhere (except in start_init_subsurf, still calling grille_m - to be changed soon). 3) removed useless CPP_PHY precompilation directives, considering the possibility to run ce0l without physics is useless (ce0l is dedicated to Earth physics).
Location:
LMDZ5/trunk/libf/dyn3d_common
Files:
1 edited
7 moved

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/dyn3d_common/caldyn0.F90

    r2334 r2336  
    1 SUBROUTINE caldyn0(itau,ucov,vcov,teta,ps,pk,phis,phi,w,time)
     1SUBROUTINE caldyn0(itau,ucov,vcov,teta,ps,masse,pk,phis,phi,w,pbaru,pbarv,time)
    22!
    33!-------------------------------------------------------------------------------
    … …  
    1515!===============================================================================
    1616! Arguments:
    17   INTEGER, INTENT(IN)  :: itau                      !---
     17  INTEGER, INTENT(IN)  :: itau                      !--- TIME STEP INDEX
    1818  REAL,    INTENT(IN)  :: vcov (ip1jm    ,llm)      !--- V COVARIANT WIND
    1919  REAL,    INTENT(IN)  :: ucov (ip1jmp1  ,llm)      !--- U COVARIANT WIND
    2020  REAL,    INTENT(IN)  :: teta (ip1jmp1  ,llm)      !--- POTENTIAL TEMPERATURE
    2121  REAL,    INTENT(IN)  :: ps   (ip1jmp1)            !--- GROUND PRESSURE
     22  REAL,    INTENT(OUT) :: masse(ip1jmp1  ,llm)      !--- MASS IN EACH CELL
    2223  REAL,    INTENT(IN)  :: pk   (iip1,jjp1,llm)      !--- PRESSURE
    2324  REAL,    INTENT(IN)  :: phis (ip1jmp1)            !--- GROUND GEOPOTENTIAL
    2425  REAL,    INTENT(IN)  :: phi  (ip1jmp1  ,llm)      !--- 3D GEOPOTENTIAL
    2526  REAL,    INTENT(OUT) :: w    (ip1jmp1  ,llm)      !--- VERTICAL WIND
     27  REAL,    INTENT(OUT) :: pbaru(ip1jmp1  ,llm)      !--- U MASS FLUX
     28  REAL,    INTENT(OUT) :: pbarv(ip1jm    ,llm)      !--- V MASS FLUX
    2629  REAL,    INTENT(IN)  :: time                      !--- TIME
    2730!===============================================================================
    2831! Local variables:
    29   REAL masse(ip1jmp1  ,llm)      !--- MASS IN EACH CELL
    30   REAL pbaru(ip1jmp1  ,llm)      !--- U MASS FLUX
    31   REAL pbarv(ip1jm    ,llm)      !--- V MASS FLUX
    3232  REAL, DIMENSION(ip1jmp1,llmp1) :: p
    3333  REAL, DIMENSION(ip1jmp1,llm)   :: ucont, massebx, ang, ecin, convm, bern

  • LMDZ5/trunk/libf/dyn3d_common/convmas.F90

    r2335 r2336  
     1SUBROUTINE convmas (pbaru, pbarv, convm)
    12!
    2 ! $Header$
    3 !
    4       SUBROUTINE convmas (pbaru, pbarv, convm )
    5 c
    6       IMPLICIT NONE
     3!-------------------------------------------------------------------------------
     4! Authors: P. Le Van , Fr. Hourdin.
     5!-------------------------------------------------------------------------------
     6! Purpose: Compute mass flux convergence at p levels.
     7  IMPLICIT NONE
     8  include "dimensions.h"
     9  include "paramet.h"
     10  include "comgeom.h"
     11  include "logic.h"
     12!===============================================================================
     13! Arguments:
     14  REAL, INTENT(IN)  :: pbaru(ip1jmp1,llm)
     15  REAL, INTENT(IN)  :: pbarv(ip1jm  ,llm)
     16  REAL, INTENT(OUT) :: convm(ip1jmp1,llm)
     17!===============================================================================
     18! Method used:   Computation from top to bottom.
     19!   Mass convergence at level llm is equal to zero and is not stored in convm.
     20!===============================================================================
     21! Local variables:
     22  INTEGER :: l
     23!===============================================================================
    724
    8 c=======================================================================
    9 c
    10 c   Auteurs:  P. Le Van , F. Hourdin  .
    11 c   -------
    12 c
    13 c   Objet:
    14 c   ------
    15 c
    16 c   ********************************************************************
    17 c   .... calcul de la convergence du flux de masse aux niveaux p ...
    18 c   ********************************************************************
    19 c
    20 c
    21 c     pbaru  et  pbarv  sont des arguments d'entree pour le s-pg  ....
    22 c      .....  convm      est  un argument de sortie pour le s-pg  ....
    23 c
    24 c    le calcul se fait de haut en bas,
    25 c    la convergence de masse au niveau p(llm+1) est egale a 0. et
    26 c    n'est pas stockee dans le tableau convm .
    27 c
    28 c
    29 c=======================================================================
    30 c
    31 c   Declarations:
    32 c   -------------
     25!--- Computation of - (d(pbaru)/dx + d(pbarv)/dy )
     26  CALL convflu( pbaru, pbarv, llm, convm )
    3327
    34 #include "dimensions.h"
    35 #include "paramet.h"
    36 #include "comvert.h"
    37 #include "logic.h"
     28!--- Filter
     29  CALL filtreg( convm, jjp1, llm, 2, 2, .TRUE., 1 )
    3830
    39       REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm ),convm(  ip1jmp1,llm )
    40       INTEGER   l,ij
     31!--- Mass convergence is integrated from top to bottom
     32  DO l=llmm1,1,-1
     33    convm(:,l) = convm(:,l) + convm(:,l+1)
     34  END DO
    4135
    42 
    43 c-----------------------------------------------------------------------
    44 c    ....  calcul de - (d(pbaru)/dx + d(pbarv)/dy ) ......
    45 
    46       CALL  convflu( pbaru, pbarv, llm, convm )
    47 
    48 c-----------------------------------------------------------------------
    49 c   filtrage:
    50 c   ---------
    51 
    52        CALL filtreg( convm, jjp1, llm, 2, 2, .true., 1 )
    53 
    54 c    integration de la convergence de masse de haut  en bas ......
    55 
    56       DO      l      = llmm1, 1, -1
    57         DO    ij     = 1, ip1jmp1
    58          convm(ij,l) = convm(ij,l) + convm(ij,l+1)
    59         ENDDO
    60       ENDDO
    61 c
    62       RETURN
    63       END
     36END SUBROUTINE convmas

  • LMDZ5/trunk/libf/dyn3d_common/enercin.F90

    r2335 r2336  
     1SUBROUTINE enercin ( vcov, ucov, vcont, ucont, ecin )
    12!
    2 ! $Header$
     3!-------------------------------------------------------------------------------
     4! Authors: P. Le Van.
     5!-------------------------------------------------------------------------------
     6! Purpose: Compute kinetic energy at sigma levels.
     7  IMPLICIT NONE
     8  include "dimensions.h"
     9  include "paramet.h"
     10  include "comgeom.h"
     11!===============================================================================
     12! Arguments:
     13  REAL, INTENT(IN)  :: vcov    (ip1jm,  llm)
     14  REAL, INTENT(IN)  :: ucov    (ip1jmp1,llm)
     15  REAL, INTENT(IN)  :: vcont   (ip1jm,  llm)
     16  REAL, INTENT(IN)  :: ucont   (ip1jmp1,llm)
     17  REAL, INTENT(OUT) :: ecin    (ip1jmp1,llm)
     18!===============================================================================
     19! Notes:
     20!                 . V
     21!                i,j-1
    322!
    4       SUBROUTINE enercin ( vcov, ucov, vcont, ucont, ecin )
    5       IMPLICIT NONE
     23!      alpha4 .       . alpha1
     24!
     25!
     26!        U .      . P     . U
     27!       i-1,j    i,j      i,j
     28!
     29!      alpha3 .       . alpha2
     30!
     31!
     32!                 . V
     33!                i,j
     34!
     35! Kinetic energy at scalar point P(i,j) (excluding poles) is:
     36!       Ecin = 0.5 * U(i-1,j)**2 *( alpha3 + alpha4 )  +
     37!              0.5 * U(i  ,j)**2 *( alpha1 + alpha2 )  +
     38!              0.5 * V(i,j-1)**2 *( alpha1 + alpha4 )  +
     39!              0.5 * V(i,  j)**2 *( alpha2 + alpha3 )
     40!===============================================================================
     41! Local variables:
     42  INTEGER :: l, ij, i
     43  REAL    :: ecinni(iip1), ecinsi(iip1), ecinpn, ecinps
     44!===============================================================================
     45  DO l=1,llm
     46    DO ij = iip2, ip1jm -1
     47      ecin(ij+1,l)=0.5*(ucov(ij    ,l)*ucont(ij    ,l)*alpha3p4(ij +1)          &
     48                      + ucov(ij+1  ,l)*ucont(ij+1  ,l)*alpha1p2(ij +1)          &
     49                      + vcov(ij-iim,l)*vcont(ij-iim,l)*alpha1p4(ij +1)          &
     50                      + vcov(ij+1  ,l)*vcont(ij+1  ,l)*alpha2p3(ij +1) )
     51    END DO
     52    !--- Correction: ecin(1,j,l)= ecin(iip1,j,l)
     53    DO ij=iip2,ip1jm,iip1; ecin(ij,l) = ecin(ij+iim,l); END DO
    654
    7 c=======================================================================
    8 c
    9 c   Auteur: P. Le Van
    10 c   -------
    11 c
    12 c   Objet:
    13 c   ------
    14 c
    15 c *********************************************************************
    16 c .. calcul de l'energie cinetique aux niveaux s  ......
    17 c *********************************************************************
    18 c  vcov, vcont, ucov et ucont sont des arguments d'entree pour le s-pg .
    19 c  ecin         est  un  argument de sortie pour le s-pg
    20 c
    21 c=======================================================================
     55    !--- North pole
     56    DO i=1,iim
     57      ecinni(i) = vcov(i,l)*vcont(i,l)*aire(i)
     58    END DO
     59    ecinpn = 0.5*SUM(ecinni(1:iim))/apoln
     60    DO ij=1,iip1; ecin(ij,l)=ecinpn; END DO
    2261
    23 #include "dimensions.h"
    24 #include "paramet.h"
    25 #include "comgeom.h"
     62    !--- South pole
     63    DO i=1,iim
     64      ecinsi(i) = vcov(i+ip1jmi1,l)*vcont(i+ip1jmi1,l)*aire(i+ip1jm)
     65    END DO
     66    ecinps = 0.5*SUM(ecinsi(1:iim))/apols
     67    DO ij=1,iip1; ecin(ij+ip1jm,l)=ecinps; END DO
     68  END DO
    2669
    27       REAL vcov( ip1jm,llm ),vcont( ip1jm,llm ),
    28      * ucov( ip1jmp1,llm ),ucont( ip1jmp1,llm ),ecin( ip1jmp1,llm )
     70END SUBROUTINE enercin
    2971
    30       REAL ecinni( iip1 ),ecinsi( iip1 )
    31 
    32       REAL ecinpn, ecinps
    33       INTEGER     l,ij,i
    34 
    35       REAL        SSUM
    36 
    37 
    38 
    39 c                 . V
    40 c                i,j-1
    41 
    42 c      alpha4 .       . alpha1
    43 
    44 
    45 c        U .      . P     . U
    46 c       i-1,j    i,j      i,j
    47 
    48 c      alpha3 .       . alpha2
    49 
    50 
    51 c                 . V
    52 c                i,j
    53 
    54 c   
    55 c  L'energie cinetique au point scalaire P(i,j) ,autre que les poles, est :
    56 c       Ecin = 0.5 * U(i-1,j)**2 *( alpha3 + alpha4 )  +
    57 c              0.5 * U(i  ,j)**2 *( alpha1 + alpha2 )  +
    58 c              0.5 * V(i,j-1)**2 *( alpha1 + alpha4 )  +
    59 c              0.5 * V(i,  j)**2 *( alpha2 + alpha3 )
    60 
    61 
    62       DO 5 l = 1,llm
    63 
    64       DO 1  ij = iip2, ip1jm -1
    65       ecin( ij+1, l )  =    0.5  *
    66      * (   ucov( ij   ,l ) * ucont( ij   ,l ) * alpha3p4( ij +1 )   +
    67      *     ucov( ij+1 ,l ) * ucont( ij+1 ,l ) * alpha1p2( ij +1 )   +
    68      *     vcov(ij-iim,l ) * vcont(ij-iim,l ) * alpha1p4( ij +1 )   +
    69      *     vcov( ij+ 1,l ) * vcont( ij+ 1,l ) * alpha2p3( ij +1 )   )
    70    1  CONTINUE
    71 
    72 c    ... correction pour  ecin(1,j,l)  ....
    73 c    ...   ecin(1,j,l)= ecin(iip1,j,l) ...
    74 
    75 CDIR$ IVDEP
    76       DO 2 ij = iip2, ip1jm, iip1
    77       ecin( ij,l ) = ecin( ij + iim, l )
    78    2  CONTINUE
    79 
    80 c     calcul aux poles  .......
    81 
    82 
    83       DO 3 i = 1, iim
    84       ecinni(i) = vcov(    i  ,  l) * vcont(    i    ,l) * aire(   i   )
    85       ecinsi(i) = vcov(i+ip1jmi1,l) * vcont(i+ip1jmi1,l) * aire(i+ip1jm)
    86    3  CONTINUE
    87 
    88       ecinpn = 0.5 * SSUM( iim,ecinni,1 ) / apoln
    89       ecinps = 0.5 * SSUM( iim,ecinsi,1 ) / apols
    90 
    91       DO 4 ij = 1,iip1
    92       ecin(   ij     , l ) = ecinpn
    93       ecin( ij+ ip1jm, l ) = ecinps
    94    4  CONTINUE
    95 
    96    5  CONTINUE
    97       RETURN
    98       END

  • LMDZ5/trunk/libf/dyn3d_common/flumass.F90

    r2335 r2336  
     1SUBROUTINE flumass (massebx,masseby, vcont, ucont, pbaru, pbarv )
    12!
    2 ! $Header$
    3 !
    4       SUBROUTINE flumass (massebx,masseby, vcont, ucont, pbaru, pbarv )
     3!-------------------------------------------------------------------------------
     4! Authors: P. Le Van , Fr. Hourdin.
     5!-------------------------------------------------------------------------------
     6! Purpose: Compute mass flux at s levels.
     7  IMPLICIT NONE
     8  include "dimensions.h"
     9  include "paramet.h"
     10  include "comgeom.h"
     11!===============================================================================
     12! Arguments:
     13  REAL, INTENT(IN)  :: massebx(ip1jmp1,llm)
     14  REAL, INTENT(IN)  :: masseby(ip1jm  ,llm)
     15  REAL, INTENT(IN)  :: vcont  (ip1jm  ,llm)
     16  REAL, INTENT(IN)  :: ucont  (ip1jmp1,llm)
     17  REAL, INTENT(OUT) :: pbaru  (ip1jmp1,llm)
     18  REAL, INTENT(OUT) :: pbarv  (ip1jm  ,llm)
     19!===============================================================================
     20! Method used:   A 2 equations system is solved.
     21!   * 1st one describes divergence computation at pole point nr. i (i=1 to im):
     22!     (0.5*(pbaru(i)-pbaru(i-1))-pbarv(i))/aire(i) = - SUM(pbarv(n))/aire pole
     23!   * 2nd one specifies that mean mass flux at pole is equal to 0:
     24!     SUM(pbaru(n)*local_area(n))=0
     25! This way, we determine additive constant common to pbary elements representing
     26!   pbaru(0,j,l) in divergence computation equation for point i=1. (i=1 to im)
     27!===============================================================================
     28! Local variables:
     29  REAL    :: sairen, saireun, ctn, ctn0, apbarun(iip1)
     30  REAL    :: saires, saireus, cts, cts0, apbarus(iip1)
     31  INTEGER :: l, i
     32!===============================================================================
     33  DO l=1,llm
     34    pbaru(iip2:ip1jm,l)=massebx(iip2:ip1jm,l)*ucont(iip2:ip1jm,l)
     35    pbarv(   1:ip1jm,l)=masseby(   1:ip1jm,l)*vcont(   1:ip1jm,l)
     36  END DO
    537
    6       IMPLICIT NONE
     38  !--- NORTH POLE
     39  sairen =SUM(aire (1:iim))
     40  saireun=SUM(aireu(1:iim))
     41  DO l = 1,llm
     42    ctn=SUM(pbarv(1:iim,l))/sairen
     43    pbaru(1,l)= pbarv(1,l)-ctn*aire(1)
     44    DO i=2,iim
     45      pbaru(i,l)=pbaru(i-1,l)+pbarv(i,l)-ctn*aire(i)
     46    END DO
     47    DO i=1,iim
     48      apbarun(i)=aireu(i)*pbaru(i,l)
     49    END DO
     50    ctn0 = -SUM(apbarun(1:iim))/saireun
     51    DO i = 1,iim
     52      pbaru(i,l)=2.*(pbaru(i,l)+ctn0)
     53    END DO
     54    pbaru(iip1,l)=pbaru(1,l)
     55  END DO
    756
    8 c=======================================================================
    9 c
    10 c   Auteurs:  P. Le Van, F. Hourdin  .
    11 c   -------
    12 c
    13 c   Objet:
    14 c   ------
    15 c
    16 c *********************************************************************
    17 c     .... calcul du flux de masse  aux niveaux s ......
    18 c *********************************************************************
    19 c   massebx,masseby,vcont et ucont sont des argum. d'entree pour le s-pg .
    20 c       pbaru  et pbarv            sont des argum.de sortie pour le s-pg .
    21 c
    22 c=======================================================================
     57  !--- SOUTH POLE
     58  saires =SUM(aire (ip1jm+1:ip1jmp1-1))
     59  saireus=SUM(aireu(ip1jm+1:ip1jmp1-1))
     60  DO l = 1,llm
     61    cts=SUM(pbarv(ip1jmi1+1:ip1jm-1,l))/saires
     62    pbaru(1+ip1jm,l)=-pbarv(1+ip1jmi1,l)+cts*aire(1+ip1jm)
     63    DO i=2,iim
     64      pbaru(i+ip1jm,l)=pbaru(i-1+ip1jm,l)-pbarv(i+ip1jmi1,l)+cts*aire(i+ip1jm)
     65    END DO
     66    DO i=1,iim
     67      apbarus(i)=aireu(i+ip1jm)*pbaru(i+ip1jm,l)
     68    END DO
     69    cts0 = -SUM(apbarus(1:iim))/saireus
     70    DO i = 1,iim
     71      pbaru(i+ip1jm,l)=2.*(pbaru(i+ip1jm,l)+cts0)
     72    END DO
     73    pbaru(ip1jmp1,l)=pbaru(1+ip1jm,l)
     74  END DO
    2375
    24 
    25 #include "dimensions.h"
    26 #include "paramet.h"
    27 #include "comgeom.h"
    28 
    29       REAL massebx( ip1jmp1,llm ),masseby( ip1jm,llm ) ,
    30      * vcont( ip1jm,llm ),ucont( ip1jmp1,llm ),pbaru( ip1jmp1,llm ),
    31      * pbarv( ip1jm,llm )
    32 
    33       REAL apbarun( iip1 ),apbarus( iip1 )
    34 
    35       REAL sairen,saireun,saires,saireus,ctn,cts,ctn0,cts0
    36       INTEGER  l,ij,i
    37 
    38       REAL       SSUM
    39 
    40 
    41       DO  5 l = 1,llm
    42 
    43       DO  1 ij = iip2,ip1jm
    44       pbaru( ij,l ) = massebx( ij,l ) * ucont( ij,l )
    45    1  CONTINUE
    46 
    47       DO 3 ij = 1,ip1jm
    48       pbarv( ij,l ) = masseby( ij,l ) * vcont( ij,l )
    49    3  CONTINUE
    50 
    51    5  CONTINUE
    52 
    53 c    ................................................................
    54 c     calcul de la composante du flux de masse en x aux poles .......
    55 c    ................................................................
    56 c     par la resolution d'1 systeme de 2 equations .
    57 
    58 c     la premiere equat.decrivant le calcul de la divergence en 1 point i
    59 c     du pole,ce calcul etant itere de i=1 a i=im .
    60 c                 c.a.d   ,
    61 c     ( ( 0.5*pbaru(i)-0.5*pbaru(i-1) - pbarv(i))/aire(i)   =
    62 c                                           - somme de ( pbarv(n) )/aire pole
    63 
    64 c     l'autre equat.specifiant que la moyenne du flux de masse au pole est =0.
    65 c     c.a.d    somme de pbaru(n)*aire locale(n) = 0.
    66 
    67 c     on en revient ainsi a determiner la constante additive commune aux pbaru
    68 c     qui representait pbaru(0,j,l) dans l'equat.du calcul de la diverg.au pt
    69 c     i=1 .
    70 c     i variant de 1 a im
    71 c     n variant de 1 a im
    72 
    73       sairen = SSUM( iim,  aire(   1     ), 1 )
    74       saireun= SSUM( iim, aireu(   1     ), 1 )
    75       saires = SSUM( iim,  aire( ip1jm+1 ), 1 )
    76       saireus= SSUM( iim, aireu( ip1jm+1 ), 1 )
    77 
    78       DO 20 l = 1,llm
    79 
    80       ctn =  SSUM( iim, pbarv(    1     ,l),  1 )/ sairen
    81       cts =  SSUM( iim, pbarv(ip1jmi1+ 1,l),  1 )/ saires
    82 
    83       pbaru(    1   ,l )=   pbarv(    1     ,l ) - ctn * aire(    1    )
    84       pbaru( ip1jm+1,l )= - pbarv( ip1jmi1+1,l ) + cts * aire( ip1jm+1 )
    85 
    86       DO 11 i = 2,iim
    87       pbaru(    i    ,l ) = pbaru(   i - 1   ,l )    +
    88      *                      pbarv(    i      ,l ) - ctn * aire(   i    )
    89 
    90       pbaru( i+ ip1jm,l ) = pbaru( i+ ip1jm-1,l )    -
    91      *                      pbarv( i+ ip1jmi1,l ) + cts * aire(i+ ip1jm)
    92   11  CONTINUE
    93       DO 12 i = 1,iim
    94       apbarun(i) = aireu(    i   ) * pbaru(   i    , l)
    95       apbarus(i) = aireu(i +ip1jm) * pbaru(i +ip1jm, l)
    96   12  CONTINUE
    97       ctn0 = -SSUM( iim,apbarun,1 )/saireun
    98       cts0 = -SSUM( iim,apbarus,1 )/saireus
    99       DO 14 i = 1,iim
    100       pbaru(   i    , l) = 2. * ( pbaru(   i    , l) + ctn0 )
    101       pbaru(i+ ip1jm, l) = 2. * ( pbaru(i +ip1jm, l) + cts0 )
    102   14  CONTINUE
    103 
    104       pbaru(   iip1 ,l ) = pbaru(    1    ,l )
    105       pbaru( ip1jmp1,l ) = pbaru( ip1jm +1,l )
    106   20  CONTINUE
    107 
    108       RETURN
    109       END
     76END SUBROUTINE flumass

  • LMDZ5/trunk/libf/dyn3d_common/massbar.F90

    r2335 r2336  
     1SUBROUTINE massbar(masse,massebx,masseby)
    12!
    2 ! $Header$
     3!-------------------------------------------------------------------------------
     4! Authors: P. Le Van , Fr. Hourdin.
     5!-------------------------------------------------------------------------------
     6! Purpose: Compute air mass mean along X and Y in each cell.
     7! See iniconst for more details.
     8  IMPLICIT NONE
     9  include "dimensions.h"
     10  include "paramet.h"
     11  include "comgeom.h"
     12!===============================================================================
     13! Arguments:
     14  REAL, INTENT(IN)  :: masse  (ip1jmp1,llm)
     15  REAL, INTENT(OUT) :: massebx(ip1jmp1,llm)
     16  REAL, INTENT(OUT) :: masseby(ip1jm  ,llm)
     17!-------------------------------------------------------------------------------
     18! Method used. Each scalar point is associated to 4 area coefficients:
     19!    * alpha1(i,j) at point ( i+1/4,j-1/4 )
     20!    * alpha2(i,j) at point ( i+1/4,j+1/4 )
     21!    * alpha3(i,j) at point ( i-1/4,j+1/4 )
     22!    * alpha4(i,j) at point ( i-1/4,j-1/4 )
     23! where alpha1(i,j) = aire(i+1/4,j-1/4)/ aire(i,j)
    324!
    4       SUBROUTINE massbar(  masse, massebx, masseby )
    5       IMPLICIT NONE
    6 c
    7 c **********************************************************************
    8 c
    9 c  Calcule les moyennes en x et  y de la masse d'air dans chaque maille.
    10 c **********************************************************************
    11 c    Auteurs : P. Le Van , Fr. Hourdin  .
    12 c   ..........
    13 c
    14 c  ..  masse                 est  un argum. d'entree  pour le s-pg ...
    15 c  ..  massebx,masseby      sont des argum. de sortie pour le s-pg ...
    16 c     
    17 c
    18 c     IMPLICIT NONE
    19 c
    20 #include "dimensions.h"
    21 #include "paramet.h"
    22 #include "comconst.h"
    23 #include "comgeom.h"
    24 c
    25       REAL    masse( ip1jmp1,llm ), massebx( ip1jmp1,llm )  ,
    26      *      masseby(   ip1jm,llm )
    27       INTEGER ij,l
    28 c
    29 c
    30 c   Methode pour calculer massebx et masseby .
    31 c   ----------------------------------------
    32 c
    33 c    A chaque point scalaire P (i,j) est affecte 4 coefficients d'aires
    34 c       alpha1(i,j)  calcule  au point ( i+1/4,j-1/4 )
    35 c       alpha2(i,j)  calcule  au point ( i+1/4,j+1/4 )
    36 c       alpha3(i,j)  calcule  au point ( i-1/4,j+1/4 )
    37 c       alpha4(i,j)  calcule  au point ( i-1/4,j-1/4 )
    38 c
    39 c    Avec  alpha1(i,j) = aire(i+1/4,j-1/4)/ aire(i,j)       
    40 c
    41 c    N.B .  Pour plus de details, voir s-pg  ...  iniconst ...
    42 c
    43 c
    44 c
    45 c   alpha4 .         . alpha1    . alpha4
    46 c    (i,j)             (i,j)       (i+1,j)
    47 c
    48 c             P .        U .          . P
    49 c           (i,j)       (i,j)         (i+1,j)
    50 c
    51 c   alpha3 .         . alpha2    .alpha3
    52 c    (i,j)              (i,j)     (i+1,j)
    53 c
    54 c             V .        Z .          . V
    55 c           (i,j)
    56 c
    57 c   alpha4 .         . alpha1    .alpha4
    58 c   (i,j+1)            (i,j+1)   (i+1,j+1)
    59 c
    60 c             P .        U .          . P
    61 c          (i,j+1)                    (i+1,j+1)
    62 c
    63 c
    64 c
    65 c                       On  a :
    66 c
    67 c    massebx(i,j) = masse(i  ,j) * ( alpha1(i  ,j) + alpha2(i,j))   +
    68 c                   masse(i+1,j) * ( alpha3(i+1,j) + alpha4(i+1,j) )
    69 c     localise  au point  ... U (i,j) ...
    70 c
    71 c    masseby(i,j) = masse(i,j  ) * ( alpha2(i,j  ) + alpha3(i,j  )  +
    72 c                   masse(i,j+1) * ( alpha1(i,j+1) + alpha4(i,j+1) 
    73 c     localise  au point  ... V (i,j) ...
    74 c
    75 c
    76 c=======================================================================
     25!   alpha4 .         . alpha1    . alpha4
     26!    (i,j)             (i,j)       (i+1,j)
     27!
     28!             P .        U .          . P
     29!           (i,j)       (i,j)         (i+1,j)
     30!
     31!   alpha3 .         . alpha2    .alpha3
     32!    (i,j)              (i,j)     (i+1,j)
     33!
     34!             V .        Z .          . V
     35!           (i,j)
     36!
     37!   alpha4 .         . alpha1    .alpha4
     38!   (i,j+1)            (i,j+1)   (i+1,j+1)
     39!
     40!             P .        U .          . P
     41!          (i,j+1)                    (i+1,j+1)
     42!
     43!
     44!    massebx(i,j) = masse(i  ,j) * ( alpha1(i  ,j) + alpha2(i,j))   +
     45!                   masse(i+1,j) * ( alpha3(i+1,j) + alpha4(i+1,j) )
     46!     localized at point  ... U (i,j) ...
     47!
     48!    masseby(i,j) = masse(i,j  ) * ( alpha2(i,j  ) + alpha3(i,j  )  +
     49!                   masse(i,j+1) * ( alpha1(i,j+1) + alpha4(i,j+1) 
     50!     localized at point  ... V (i,j) ...
     51!===============================================================================
     52! Local variables:
     53  INTEGER :: ij, l
     54!===============================================================================
     55  DO l=1,llm
     56    DO ij=1,ip1jmp1-1
     57      massebx(ij,l)=masse(ij,l)*alpha1p2(ij)+masse(ij+1   ,l)*alpha3p4(ij+1)
     58    END DO
     59    DO ij=iip1,ip1jmp1,iip1; massebx(ij,l)=massebx(ij-iim,l); END DO
     60    DO ij=1,ip1jm
     61      masseby(ij,l)=masse(ij,l)*alpha2p3(ij)+masse(ij+iip1,l)*alpha1p4(ij+iip1)
     62    END DO
     63  END DO
    7764
    78       DO   100    l = 1 , llm
    79 c
    80         DO  ij = 1, ip1jmp1 - 1
    81          massebx(ij,l) =  masse( ij, l) * alpha1p2( ij  )     +
    82      *                   masse(ij+1, l) * alpha3p4(ij+1 )
    83         ENDDO
     65END SUBROUTINE massbar
    8466
    85 c    .... correction pour massebx( iip1,j) .....
    86 c    ...    massebx(iip1,j)= massebx(1,j) ...
    87 c
    88 CDIR$ IVDEP
    89         DO  ij = iip1, ip1jmp1, iip1
    90          massebx( ij,l ) = massebx( ij - iim,l )
    91         ENDDO
    92 
    93 
    94          DO  ij = 1,ip1jm
    95          masseby( ij,l ) = masse(  ij   , l ) * alpha2p3(   ij    )  +
    96      *                     masse(ij+iip1, l ) * alpha1p4( ij+iip1 )
    97          ENDDO
    98 
    99 100   CONTINUE
    100 c
    101       RETURN
    102       END

  • LMDZ5/trunk/libf/dyn3d_common/massbarxy.F90

    r2335 r2336  
     1SUBROUTINE massbarxy(masse,massebxy)
    12!
    2 ! $Header$
    3 !
    4       SUBROUTINE massbarxy(  masse, massebxy )
    5       IMPLICIT NONE
    6 c
    7 c **********************************************************************
    8 c
    9 c  Calcule les moyennes en x et  y de la masse d'air dans chaque maille.
    10 c **********************************************************************
    11 c    Auteurs : P. Le Van , Fr. Hourdin  .
    12 c   ..........
    13 c
    14 c  ..  masse          est  un  argum. d'entree  pour le s-pg ...
    15 c  ..  massebxy       est  un  argum. de sortie pour le s-pg ...
    16 c     
    17 c
    18 c     IMPLICIT NONE
    19 c
    20 #include "dimensions.h"
    21 #include "paramet.h"
    22 #include "comconst.h"
    23 #include "comgeom.h"
    24 c
    25        REAL  masse( ip1jmp1,llm ), massebxy( ip1jm,llm )
    26        INTEGER ij,l
    27 c
     3!-------------------------------------------------------------------------------
     4! Authors: P. Le Van , Fr. Hourdin.
     5!-------------------------------------------------------------------------------
     6! Purpose: Compute air mass mean along X and Y in each cell.
     7! See iniconst for more details.
     8  IMPLICIT NONE
     9  include "dimensions.h"
     10  include "paramet.h"
     11  include "comconst.h"
     12  include "comgeom.h"
     13!===============================================================================
     14! Arguments:
     15  REAL, INTENT(IN)  :: masse   (ip1jmp1,llm)
     16  REAL, INTENT(OUT) :: massebxy(ip1jm  ,llm)
     17!===============================================================================
     18! Local variables:
     19  INTEGER :: ij, l
     20!===============================================================================
     21  DO l=1,llm
     22    DO ij=1,ip1jm-1
     23      massebxy(ij,l)=masse(ij     ,l)*alpha2(ij     ) + &
     24     +               masse(ij+1   ,l)*alpha3(ij+1   ) + &
     25     +               masse(ij+iip1,l)*alpha1(ij+iip1) + &
     26     +               masse(ij+iip2,l)*alpha4(ij+iip2)
     27    END DO
     28    DO ij=iip1,ip1jm,iip1; massebxy(ij,l)=massebxy(ij-iim,l); END DO
     29  END DO
    2830
    29       DO   100    l = 1 , llm
    30 c
    31       DO 5 ij = 1, ip1jm - 1
    32       massebxy( ij,l ) = masse(    ij  ,l ) * alpha2(   ij    )   +
    33      +                   masse(   ij+1 ,l ) * alpha3(  ij+1   )   +
    34      +                   masse( ij+iip1,l ) * alpha1( ij+iip1 )   +
    35      +                   masse( ij+iip2,l ) * alpha4( ij+iip2 )
    36    5  CONTINUE
    37 
    38 c    ....  correction pour     massebxy( iip1,j )  ........
    39 
    40 CDIR$ IVDEP
    41 
    42       DO 7 ij = iip1, ip1jm, iip1
    43       massebxy( ij,l ) = massebxy( ij - iim,l )
    44    7  CONTINUE
    45 
    46 100   CONTINUE
    47 c
    48       RETURN
    49       END
     31END SUBROUTINE massbarxy

  • LMDZ5/trunk/libf/dyn3d_common/tourpot.F90

    r2335 r2336  
     1SUBROUTINE tourpot ( vcov, ucov, massebxy, vorpot )
    12!
    2 ! $Header$
    3 !
    4       SUBROUTINE tourpot ( vcov, ucov, massebxy, vorpot )
    5       IMPLICIT NONE
     3!-------------------------------------------------------------------------------
     4! Authors: P. Le Van.
     5!-------------------------------------------------------------------------------
     6! Purpose: Compute potential vorticity.
     7  IMPLICIT NONE
     8  include "dimensions.h"
     9  include "paramet.h"
     10  include "comgeom.h"
     11  include "logic.h"
     12!===============================================================================
     13! Arguments:
     14  REAL, INTENT(IN)  :: vcov    (ip1jm,  llm)
     15  REAL, INTENT(IN)  :: ucov    (ip1jmp1,llm)
     16  REAL, INTENT(IN)  :: massebxy(ip1jm,  llm)
     17  REAL, INTENT(OUT) :: vorpot  (ip1jm,  llm)
     18!===============================================================================
     19! Method used:
     20!   vorpot = ( Filtre( d(vcov)/dx - d(ucov)/dy ) + fext ) /psbarxy
     21!===============================================================================
     22! Local variables:
     23  INTEGER :: l, ij
     24  REAL    :: rot(ip1jm,llm)
     25!===============================================================================
    626
    7 c=======================================================================
    8 c
    9 c   Auteur:  P. Le Van
    10 c   -------
    11 c
    12 c   Objet:
    13 c   ------
    14 c
    15 c    *******************************************************************
    16 c    .........      calcul du tourbillon potentiel             .........
    17 c    *******************************************************************
    18 c
    19 c     vcov,ucov,fext et pbarxyfl sont des argum. d'entree pour le s-pg .
    20 c             vorpot            est  un argum.de sortie pour le s-pg .
    21 c
    22 c=======================================================================
     27!--- Wind vorticity ; correction: rot(iip1,j,l) = rot(1,j,l)
     28  DO l=1,llm
     29    DO ij=1,ip1jm-1
     30      rot(ij,l)=vcov(ij+1,l)-vcov(ij,l)+ucov(ij+iip1,l)-ucov(ij,l)
     31    END DO
     32    DO ij=iip1,ip1jm,iip1; rot(ij,l)=rot(ij-iim,l); END DO
     33  END DO
    2334
    24 #include "dimensions.h"
    25 #include "paramet.h"
    26 #include "comgeom.h"
    27 #include "logic.h"
     35!--- Filter
     36  CALL  filtreg(rot,jjm,llm,2,1,.FALSE.,1)
    2837
    29       REAL  rot( ip1jm,llm )
    30       REAL vcov( ip1jm,llm ),ucov( ip1jmp1,llm )
    31       REAL massebxy( ip1jm,llm ),vorpot( ip1jm,llm )
     38!--- Potential vorticity ; correction: rot(iip1,j,l) = rot(1,j,l)
     39  DO l=1,llm
     40    DO ij=1,ip1jm-1
     41      vorpot(ij,l)=(rot(ij,l)+fext(ij))/massebxy(ij,l)
     42    END DO
     43    DO ij=iip1,ip1jm,iip1; vorpot(ij,l)=vorpot(ij-iim,l); END DO
     44  END DO
    3245
    33       INTEGER l, ij
    34 
    35 
    36 
    37 
    38 c  ... vorpot = ( Filtre( d(vcov)/dx - d(ucov)/dy ) + fext ) /psbarxy ..
    39 
    40 
    41 
    42 c    ........  Calcul du rotationnel du vent V  puis filtrage  ........
    43 
    44       DO 5 l = 1,llm
    45 
    46       DO 2 ij = 1, ip1jm - 1
    47       rot( ij,l ) = vcov(ij+1,l)-vcov(ij,l)+ucov(ij+iip1,l)-ucov(ij,l)
    48    2  CONTINUE
    49 
    50 c    ....  correction pour  rot( iip1,j,l )  .....
    51 c    ....     rot(iip1,j,l) = rot(1,j,l)    .....
    52 
    53 CDIR$ IVDEP
    54 
    55       DO 3 ij = iip1, ip1jm, iip1
    56       rot( ij,l ) = rot( ij -iim, l )
    57    3  CONTINUE
    58 
    59    5  CONTINUE
    60 
    61 
    62       CALL  filtreg( rot, jjm, llm, 2, 1, .FALSE., 1 )
    63 
    64 
    65       DO 10 l = 1, llm
    66 
    67       DO 6 ij = 1, ip1jm - 1
    68       vorpot( ij,l ) = ( rot(ij,l) + fext(ij) ) / massebxy(ij,l)
    69    6  CONTINUE
    70 
    71 c    ..... correction pour  vorpot( iip1,j,l)  .....
    72 c    ....   vorpot(iip1,j,l)= vorpot(1,j,l) ....
    73 CDIR$ IVDEP
    74       DO 8 ij = iip1, ip1jm, iip1
    75       vorpot( ij,l ) = vorpot( ij -iim,l )
    76    8  CONTINUE
    77 
    78   10  CONTINUE
    79 
    80       RETURN
    81       END
     46END SUBROUTINE tourpot

  • LMDZ5/trunk/libf/dyn3d_common/vitvert.F90

    r2335 r2336  
     1SUBROUTINE vitvert (convm, w)
    12!
    2 ! $Header$
    3 !
    4       SUBROUTINE vitvert ( convm , w )
    5 c
    6       IMPLICIT NONE
     3!-------------------------------------------------------------------------------
     4! Authors: P. Le Van , Fr. Hourdin.
     5!-------------------------------------------------------------------------------
     6! Purpose: Compute vertical speed at sigma levels.
     7  IMPLICIT NONE
     8  include "dimensions.h"
     9  include "paramet.h"
     10  include "comvert.h"
     11!===============================================================================
     12! Arguments:
     13  REAL, INTENT(IN)  :: convm(ip1jmp1,llm)
     14  REAL, INTENT(OUT) :: w    (ip1jmp1,llm)
     15!===============================================================================
     16! Notes: Vertical speed is oriented from bottom to top.
     17!   * At ground - level sigma(1):     w(i,j,1) = 0.
     18!   * At top    - level sigma(llm+1): w(i,j,l) = 0. (not stored in w)
     19!===============================================================================
     20! Local variables:
     21  INTEGER :: l
     22!===============================================================================
     23  DO l=1,llmm1; w(:,l+1)=convm(:,l+1)-bp(l+1)*convm(:,1); END DO
     24  w(:,1)=0.
    725
    8 c=======================================================================
    9 c
    10 c   Auteurs:  P. Le Van , F. Hourdin .
    11 c   -------
    12 c
    13 c   Objet:
    14 c   ------
    15 c
    16 c    *******************************************************************
    17 c  .... calcul de la vitesse verticale aux niveaux sigma  ....
    18 c    *******************************************************************
    19 c     convm   est un argument  d'entree pour le s-pg  ......
    20 c       w     est un argument de sortie pour le s-pg  ......
    21 c
    22 c    la vitesse verticale est orientee de  haut en bas .
    23 c    au sol, au niveau sigma(1),   w(i,j,1) = 0.
    24 c    au sommet, au niveau sigma(llm+1) , la vit.verticale est aussi
    25 c    egale a 0. et n'est pas stockee dans le tableau w  .
    26 c
    27 c
    28 c=======================================================================
     26END SUBROUTINE vitvert
    2927
    30 #include "dimensions.h"
    31 #include "paramet.h"
    32 #include "comvert.h"
    33 
    34       REAL w(ip1jmp1,llm),convm(ip1jmp1,llm)
    35       INTEGER   l, ij
    36 
    37 
    38 
    39       DO 2  l = 1,llmm1
    40 
    41       DO 1 ij = 1,ip1jmp1
    42       w( ij, l+1 ) = convm( ij, l+1 ) - bp(l+1) * convm( ij, 1 )
    43    1  CONTINUE
    44 
    45    2  CONTINUE
    46 
    47       DO 5 ij  = 1,ip1jmp1
    48       w(ij,1)  = 0.
    49 5     CONTINUE
    50 
    51       RETURN
    52       END
Note: See TracChangeset for help on using the changeset viewer.