Ignore:
Timestamp:
Mar 28, 2016, 5:27:51 PM (9 years ago)
Author:
emillour
Message:

All models: More updates to make planetary codes (+Earth) setups converge.

  • in dyn3d_common:
  • convmas.F => convmas.F90
  • enercin.F => enercin.F90
  • flumass.F => flumass.F90
  • massbar.F => massbar.F90
  • tourpot.F => tourpot.F90
  • vitvert.F => vitvert.F90
  • in misc:
  • move "q_sat" from "dyn3d_common" to "misc" (in Earth model, it is also called by the physics)
  • move "write_field" from "dyn3d_common" to "misc"(may be called from physics or dynamics and depends on neither).
  • in phy_common:
  • move "write_field_phy" here since it may be called from any physics package)
  • add module "regular_lonlat_mod" to store global information on lon-lat grid
  • in dynlonlat_phylonlat/phy*:
  • turn "iniphysiq.F90" into module "iniphysiq_mod.F90" (and of course adapt gcm.F[90] and 1D models accordingly)

EM

Location:
trunk/LMDZ.COMMON/libf/dyn3d_common
Files:
7 moved

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.COMMON/libf/dyn3d_common/convmas.F90

    r1520 r1523  
     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!===============================================================================
     12! Arguments:
     13  REAL, INTENT(IN)  :: pbaru(ip1jmp1,llm)
     14  REAL, INTENT(IN)  :: pbarv(ip1jm  ,llm)
     15  REAL, INTENT(OUT) :: convm(ip1jmp1,llm)
     16!===============================================================================
     17! Method used:   Computation from top to bottom.
     18!   Mass convergence at level llm is equal to zero and is not stored in convm.
     19!===============================================================================
     20! Local variables:
     21  INTEGER :: l
     22!===============================================================================
    723
    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   -------------
     24!--- Computation of - (d(pbaru)/dx + d(pbarv)/dy )
     25  CALL convflu( pbaru, pbarv, llm, convm )
    3326
    34 #include "dimensions.h"
    35 #include "paramet.h"
     27!--- Filter
     28  CALL filtreg( convm, jjp1, llm, 2, 2, .TRUE., 1 )
    3629
    37       REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm ),convm(  ip1jmp1,llm )
    38       INTEGER   l,ij
     30!--- Mass convergence is integrated from top to bottom
     31  DO l=llmm1,1,-1
     32    convm(:,l) = convm(:,l) + convm(:,l+1)
     33  END DO
    3934
    40 
    41 c-----------------------------------------------------------------------
    42 c    ....  calcul de - (d(pbaru)/dx + d(pbarv)/dy ) ......
    43 
    44       CALL  convflu( pbaru, pbarv, llm, convm )
    45 
    46 c-----------------------------------------------------------------------
    47 c   filtrage:
    48 c   ---------
    49 
    50        CALL filtreg( convm, jjp1, llm, 2, 2, .true., 1 )
    51 
    52 c    integration de la convergence de masse de haut  en bas ......
    53 
    54       DO      l      = llmm1, 1, -1
    55         DO    ij     = 1, ip1jmp1
    56          convm(ij,l) = convm(ij,l) + convm(ij,l+1)
    57         ENDDO
    58       ENDDO
    59 c
    60       RETURN
    61       END
     35END SUBROUTINE convmas

  • trunk/LMDZ.COMMON/libf/dyn3d_common/enercin.F90

    r1520 r1523  
     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

  • trunk/LMDZ.COMMON/libf/dyn3d_common/flumass.F90

    r1520 r1523  
     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

  • trunk/LMDZ.COMMON/libf/dyn3d_common/massbar.F90

    r1520 r1523  
     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 "comgeom.h"
    23 c
    24       REAL    masse( ip1jmp1,llm ), massebx( ip1jmp1,llm )  ,
    25      *      masseby(   ip1jm,llm )
    26       INTEGER ij,l
    27 c
    28 c
    29 c   Methode pour calculer massebx et masseby .
    30 c   ----------------------------------------
    31 c
    32 c    A chaque point scalaire P (i,j) est affecte 4 coefficients d'aires
    33 c       alpha1(i,j)  calcule  au point ( i+1/4,j-1/4 )
    34 c       alpha2(i,j)  calcule  au point ( i+1/4,j+1/4 )
    35 c       alpha3(i,j)  calcule  au point ( i-1/4,j+1/4 )
    36 c       alpha4(i,j)  calcule  au point ( i-1/4,j-1/4 )
    37 c
    38 c    Avec  alpha1(i,j) = aire(i+1/4,j-1/4)/ aire(i,j)       
    39 c
    40 c    N.B .  Pour plus de details, voir s-pg  ...  iniconst ...
    41 c
    42 c
    43 c
    44 c   alpha4 .         . alpha1    . alpha4
    45 c    (i,j)             (i,j)       (i+1,j)
    46 c
    47 c             P .        U .          . P
    48 c           (i,j)       (i,j)         (i+1,j)
    49 c
    50 c   alpha3 .         . alpha2    .alpha3
    51 c    (i,j)              (i,j)     (i+1,j)
    52 c
    53 c             V .        Z .          . V
    54 c           (i,j)
    55 c
    56 c   alpha4 .         . alpha1    .alpha4
    57 c   (i,j+1)            (i,j+1)   (i+1,j+1)
    58 c
    59 c             P .        U .          . P
    60 c          (i,j+1)                    (i+1,j+1)
    61 c
    62 c
    63 c
    64 c                       On  a :
    65 c
    66 c    massebx(i,j) = masse(i  ,j) * ( alpha1(i  ,j) + alpha2(i,j))   +
    67 c                   masse(i+1,j) * ( alpha3(i+1,j) + alpha4(i+1,j) )
    68 c     localise  au point  ... U (i,j) ...
    69 c
    70 c    masseby(i,j) = masse(i,j  ) * ( alpha2(i,j  ) + alpha3(i,j  )  +
    71 c                   masse(i,j+1) * ( alpha1(i,j+1) + alpha4(i,j+1) 
    72 c     localise  au point  ... V (i,j) ...
    73 c
    74 c
    75 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
    7664
    77       DO   100    l = 1 , llm
    78 c
    79         DO  ij = 1, ip1jmp1 - 1
    80          massebx(ij,l) =  masse( ij, l) * alpha1p2( ij  )     +
    81      *                   masse(ij+1, l) * alpha3p4(ij+1 )
    82         ENDDO
     65END SUBROUTINE massbar
    8366
    84 c    .... correction pour massebx( iip1,j) .....
    85 c    ...    massebx(iip1,j)= massebx(1,j) ...
    86 c
    87 CDIR$ IVDEP
    88         DO  ij = iip1, ip1jmp1, iip1
    89          massebx( ij,l ) = massebx( ij - iim,l )
    90         ENDDO
    91 
    92 
    93          DO  ij = 1,ip1jm
    94          masseby( ij,l ) = masse(  ij   , l ) * alpha2p3(   ij    )  +
    95      *                     masse(ij+iip1, l ) * alpha1p4( ij+iip1 )
    96          ENDDO
    97 
    98 100   CONTINUE
    99 c
    100       RETURN
    101       END

  • trunk/LMDZ.COMMON/libf/dyn3d_common/massbarxy.F90

    r1520 r1523  
     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 "comgeom.h"
    23 c
    24        REAL  masse( ip1jmp1,llm ), massebxy( ip1jm,llm )
    25        INTEGER ij,l
    26 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 "comgeom.h"
     12!===============================================================================
     13! Arguments:
     14  REAL, INTENT(IN)  :: masse   (ip1jmp1,llm)
     15  REAL, INTENT(OUT) :: massebxy(ip1jm  ,llm)
     16!===============================================================================
     17! Local variables:
     18  INTEGER :: ij, l
     19!===============================================================================
     20  DO l=1,llm
     21    DO ij=1,ip1jm-1
     22      massebxy(ij,l)=masse(ij     ,l)*alpha2(ij     ) + &
     23                     masse(ij+1   ,l)*alpha3(ij+1   ) + &
     24                     masse(ij+iip1,l)*alpha1(ij+iip1) + &
     25                     masse(ij+iip2,l)*alpha4(ij+iip2)
     26    END DO
     27    DO ij=iip1,ip1jm,iip1; massebxy(ij,l)=massebxy(ij-iim,l); END DO
     28  END DO
    2729
    28       DO   100    l = 1 , llm
    29 c
    30       DO 5 ij = 1, ip1jm - 1
    31       massebxy( ij,l ) = masse(    ij  ,l ) * alpha2(   ij    )   +
    32      +                   masse(   ij+1 ,l ) * alpha3(  ij+1   )   +
    33      +                   masse( ij+iip1,l ) * alpha1( ij+iip1 )   +
    34      +                   masse( ij+iip2,l ) * alpha4( ij+iip2 )
    35    5  CONTINUE
    36 
    37 c    ....  correction pour     massebxy( iip1,j )  ........
    38 
    39 CDIR$ IVDEP
    40 
    41       DO 7 ij = iip1, ip1jm, iip1
    42       massebxy( ij,l ) = massebxy( ij - iim,l )
    43    7  CONTINUE
    44 
    45 100   CONTINUE
    46 c
    47       RETURN
    48       END
     30END SUBROUTINE massbarxy

  • trunk/LMDZ.COMMON/libf/dyn3d_common/tourpot.F90

    r1520 r1523  
     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!===============================================================================
     12! Arguments:
     13  REAL, INTENT(IN)  :: vcov    (ip1jm,  llm)
     14  REAL, INTENT(IN)  :: ucov    (ip1jmp1,llm)
     15  REAL, INTENT(IN)  :: massebxy(ip1jm,  llm)
     16  REAL, INTENT(OUT) :: vorpot  (ip1jm,  llm)
     17!===============================================================================
     18! Method used:
     19!   vorpot = ( Filtre( d(vcov)/dx - d(ucov)/dy ) + fext ) /psbarxy
     20!===============================================================================
     21! Local variables:
     22  INTEGER :: l, ij
     23  REAL    :: rot(ip1jm,llm)
     24!===============================================================================
    625
    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=======================================================================
     26!--- Wind vorticity ; correction: rot(iip1,j,l) = rot(1,j,l)
     27  DO l=1,llm
     28    DO ij=1,ip1jm-1
     29      rot(ij,l)=vcov(ij+1,l)-vcov(ij,l)+ucov(ij+iip1,l)-ucov(ij,l)
     30    END DO
     31    DO ij=iip1,ip1jm,iip1; rot(ij,l)=rot(ij-iim,l); END DO
     32  END DO
    2333
    24 #include "dimensions.h"
    25 #include "paramet.h"
    26 #include "comgeom.h"
     34!--- Filter
     35  CALL  filtreg(rot,jjm,llm,2,1,.FALSE.,1)
    2736
    28       REAL  rot( ip1jm,llm )
    29       REAL vcov( ip1jm,llm ),ucov( ip1jmp1,llm )
    30       REAL massebxy( ip1jm,llm ),vorpot( ip1jm,llm )
     37!--- Potential vorticity ; correction: rot(iip1,j,l) = rot(1,j,l)
     38  DO l=1,llm
     39    DO ij=1,ip1jm-1
     40      vorpot(ij,l)=(rot(ij,l)+fext(ij))/massebxy(ij,l)
     41    END DO
     42    DO ij=iip1,ip1jm,iip1; vorpot(ij,l)=vorpot(ij-iim,l); END DO
     43  END DO
    3144
    32       INTEGER l, ij
    33 
    34 
    35 
    36 
    37 c  ... vorpot = ( Filtre( d(vcov)/dx - d(ucov)/dy ) + fext ) /psbarxy ..
    38 
    39 
    40 
    41 c    ........  Calcul du rotationnel du vent V  puis filtrage  ........
    42 
    43       DO 5 l = 1,llm
    44 
    45       DO 2 ij = 1, ip1jm - 1
    46       rot( ij,l ) = vcov(ij+1,l)-vcov(ij,l)+ucov(ij+iip1,l)-ucov(ij,l)
    47    2  CONTINUE
    48 
    49 c    ....  correction pour  rot( iip1,j,l )  .....
    50 c    ....     rot(iip1,j,l) = rot(1,j,l)    .....
    51 
    52 CDIR$ IVDEP
    53 
    54       DO 3 ij = iip1, ip1jm, iip1
    55       rot( ij,l ) = rot( ij -iim, l )
    56    3  CONTINUE
    57 
    58    5  CONTINUE
    59 
    60 
    61       CALL  filtreg( rot, jjm, llm, 2, 1, .FALSE., 1 )
    62 
    63 
    64       DO 10 l = 1, llm
    65 
    66       DO 6 ij = 1, ip1jm - 1
    67       vorpot( ij,l ) = ( rot(ij,l) + fext(ij) ) / massebxy(ij,l)
    68    6  CONTINUE
    69 
    70 c    ..... correction pour  vorpot( iip1,j,l)  .....
    71 c    ....   vorpot(iip1,j,l)= vorpot(1,j,l) ....
    72 CDIR$ IVDEP
    73       DO 8 ij = iip1, ip1jm, iip1
    74       vorpot( ij,l ) = vorpot( ij -iim,l )
    75    8  CONTINUE
    76 
    77   10  CONTINUE
    78 
    79       RETURN
    80       END
     45END SUBROUTINE tourpot

  • trunk/LMDZ.COMMON/libf/dyn3d_common/vitvert.F90

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