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).
File:
1 moved

Legend:

Unmodified
Added
Removed

  • 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
Note: See TracChangeset for help on using the changeset viewer.