Ignore:
Timestamp:
Jul 23, 2024, 7:14:34 PM (8 weeks ago)
Author:
abarral
Message:

Replace 1DUTILS.h by module lmdz_1dutils.f90
Replace 1DConv.h by module lmdz_old_1dconv.f90 (it's only used by old_* files)
Convert *.F to *.f90
Fix gradsdef.h formatting
Remove unnecessary "RETURN" at the end of functions/subroutines

File:
1 moved

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/Amaury_dev/libf/dynphy_lonlat/calfis.f90

    r5104 r5105  
    22! $Id$
    33
    4 C
    5 C
    6       SUBROUTINE calfis(lafin,
    7      $                  jD_cur, jH_cur,
    8      $                  pucov,
    9      $                  pvcov,
    10      $                  pteta,
    11      $                  pq,
    12      $                  pmasse,
    13      $                  pps,
    14      $                  pp,
    15      $                  ppk,
    16      $                  pphis,
    17      $                  pphi,
    18      $                  pducov,
    19      $                  pdvcov,
    20      $                  pdteta,
    21      $                  pdq,
    22      $                  flxw,
    23      $                  pdufi,
    24      $                  pdvfi,
    25      $                  pdhfi,
    26      $                  pdqfi,
    27      $                  pdpsfi)
    28 c
    29 c    Auteur :  P. Le Van, F. Hourdin
    30 c   .........
    31       USE infotrac, ONLY: nqtot, tracers
    32       USE control_mod, ONLY: planet_type, nsplit_phys
    33       USE callphysiq_mod, ONLY: call_physiq
    34       USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_PHYS
    35       USE comconst_mod, ONLY: cpp, daysec, dtphys, dtvr, kappa, pi
    36       USE comvert_mod, ONLY: preff, presnivs
    37 
    38       IMPLICIT NONE
    39 c=======================================================================
    40 c
    41 c   1. rearrangement des tableaux et transformation
    42 c      variables dynamiques  >  variables physiques
    43 c   2. calcul des termes physiques
    44 c   3. retransformation des tendances physiques en tendances dynamiques
    45 c
    46 c   remarques:
    47 c   ----------
    48 c
    49 c    - les vents sont donnes dans la physique par leurs composantes
    50 c      naturelles.
    51 c    - la variable thermodynamique de la physique est une variable
    52 c      intensive :   T
    53 c      pour la dynamique on prend    T * ( preff / p(l) ) **kappa
    54 c    - les deux seules variables dependant de la geometrie necessaires
    55 c      pour la physique sont la latitude pour le rayonnement et
    56 c      l'aire de la maille quand on veut integrer une grandeur
    57 c      horizontalement.
    58 c    - les points de la physique sont les points scalaires de la
    59 c      la dynamique; numerotation:
    60 c          1 pour le pole nord
    61 c          (jjm-1)*iim pour l'interieur du domaine
    62 c          ngridmx pour le pole sud
    63 c      ---> ngridmx=2+(jjm-1)*iim
    64 c
    65 c     Input :
    66 c     -------
    67 c       pucov           covariant zonal velocity
    68 c       pvcov           covariant meridional velocity
    69 c       pteta           potential temperature
    70 c       pps             surface pressure
    71 c       pmasse          masse d'air dans chaque maille
    72 c       pts             surface temperature  (K)
    73 c       callrad         clef d'appel au rayonnement
    74 c
    75 c    Output :
    76 c    --------
    77 c        pdufi          tendency for the natural zonal velocity (ms-1)
    78 c        pdvfi          tendency for the natural meridional velocity
    79 c        pdhfi          tendency for the potential temperature
    80 c        pdtsfi         tendency for the surface temperature
    81 c
    82 c        pdtrad         radiative tendencies  \  both input
    83 c        pfluxrad       radiative fluxes      /  and output
    84 c
    85 c=======================================================================
    86 c
    87 c-----------------------------------------------------------------------
    88 c
    89 c    0.  Declarations :
    90 c    ------------------
    91 
    92       include "dimensions.h"
    93       include "paramet.h"
    94 
    95       INTEGER ngridmx
    96       PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
    97 
    98       include "comgeom2.h"
    99       include "iniprint.h"
    100 
    101 c    Arguments :
    102 c    -----------
    103       LOGICAL,INTENT(IN) ::  lafin ! .TRUE. for the very last CALL to physics
    104       REAL,INTENT(IN):: jD_cur, jH_cur
    105       REAL,INTENT(IN) :: pvcov(iip1,jjm,llm) ! covariant meridional velocity
    106       REAL,INTENT(IN) :: pucov(iip1,jjp1,llm) ! covariant zonal velocity
    107       REAL,INTENT(IN) :: pteta(iip1,jjp1,llm) ! potential temperature
    108       REAL,INTENT(IN) :: pmasse(iip1,jjp1,llm) ! mass in each cell ! not used
    109       REAL,INTENT(IN) :: pq(iip1,jjp1,llm,nqtot) ! tracers
    110       REAL,INTENT(IN) :: pphis(iip1,jjp1) ! surface geopotential
    111       REAL,INTENT(IN) :: pphi(iip1,jjp1,llm) ! geopotential
    112 
    113       REAL,INTENT(IN) :: pdvcov(iip1,jjm,llm) ! dynamical tendency on vcov
    114       REAL,INTENT(IN) :: pducov(iip1,jjp1,llm) ! dynamical tendency on ucov
    115       REAL,INTENT(IN) :: pdteta(iip1,jjp1,llm) ! dynamical tendency on teta
    116       ! NB: pdteta is used only to compute pcvgt which is in fact not used...
    117       REAL,INTENT(IN) :: pdq(iip1,jjp1,llm,nqtot) ! dynamical tendency on tracers
    118       ! NB: pdq is only used to compute pcvgq which is in fact not used...
    119 
    120       REAL,INTENT(IN) :: pps(iip1,jjp1) ! surface pressure (Pa)
    121       REAL,INTENT(IN) :: pp(iip1,jjp1,llmp1) ! pressure at mesh interfaces (Pa)
    122       REAL,INTENT(IN) :: ppk(iip1,jjp1,llm) ! Exner at mid-layer
    123       REAL,INTENT(IN) :: flxw(iip1,jjp1,llm) ! Vertical mass flux on lower mesh interfaces (kg/s) (on llm because flxw(:,:,llm+1)=0)
    124 
    125       ! tendencies (in */s) from the physics
    126       REAL,INTENT(OUT) :: pdvfi(iip1,jjm,llm) ! tendency on covariant meridional wind
    127       REAL,INTENT(OUT) :: pdufi(iip1,jjp1,llm) ! tendency on covariant zonal wind
    128       REAL,INTENT(OUT) :: pdhfi(iip1,jjp1,llm) ! tendency on potential temperature (K/s)
    129       REAL,INTENT(OUT) :: pdqfi(iip1,jjp1,llm,nqtot) ! tendency on tracers
    130       REAL,INTENT(OUT) :: pdpsfi(iip1,jjp1) ! tendency on surface pressure (Pa/s)
    131 
    132 
    133 c    Local variables :
    134 c    -----------------
    135 
    136       INTEGER i,j,l,ig0,ig,iq,itr
    137       REAL zpsrf(ngridmx)
    138       REAL zplev(ngridmx,llm+1),zplay(ngridmx,llm)
    139       REAL zphi(ngridmx,llm),zphis(ngridmx)
    140 c
    141       REAL zrot(iip1,jjm,llm) ! AdlC May 2014
    142       REAL zufi(ngridmx,llm), zvfi(ngridmx,llm)
    143       REAL zrfi(ngridmx,llm) ! relative wind vorticity
    144       REAL ztfi(ngridmx,llm),zqfi(ngridmx,llm,nqtot)
    145       REAL zpk(ngridmx,llm)
    146 c
    147       REAL pcvgu(ngridmx,llm), pcvgv(ngridmx,llm)
    148       REAL pcvgt(ngridmx,llm), pcvgq(ngridmx,llm,2)
    149 c
    150       REAL zdufi(ngridmx,llm),zdvfi(ngridmx,llm)
    151       REAL zdtfi(ngridmx,llm),zdqfi(ngridmx,llm,nqtot)
    152       REAL zdpsrf(ngridmx)
    153 c
    154       REAL zdufic(ngridmx,llm),zdvfic(ngridmx,llm)
    155       REAL zdtfic(ngridmx,llm),zdqfic(ngridmx,llm,nqtot)
    156       REAL jH_cur_split,zdt_split
    157       LOGICAL debut_split,lafin_split
    158       INTEGER isplit
    159 
    160       REAL zsin(iim),zcos(iim),z1(iim)
    161       REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
    162       REAL unskap, pksurcp
    163 c
    164       REAL flxwfi(ngridmx,llm)  ! Flux de masse verticale sur la grille physiq
    165 c
    166 
    167       REAL SSUM
    168 
    169       LOGICAL,SAVE :: firstcal=.TRUE., debut=.TRUE.
    170 !      REAL rdayvrai
    171 
    172 c
    173 c-----------------------------------------------------------------------
    174 c
    175 c    1. Initialisations :
    176 c    --------------------
    177 c
    178 c
    179       IF ( firstcal )  THEN
    180         debut = .TRUE.
    181         IF (ngridmx/=2+(jjm-1)*iim) THEN
    182          write(lunout,*) 'STOP dans calfis'
    183          write(lunout,*)
    184      &   'La dimension ngridmx doit etre egale a 2 + (jjm-1)*iim'
    185          write(lunout,*) '  ngridmx  jjm   iim   '
    186          write(lunout,*) ngridmx,jjm,iim
    187          CALL abort_gcm("calfis", "", 1)
    188         ENDIF
    189       ELSE
    190         debut = .FALSE.
    191       ENDIF ! of IF (firstcal)
    192 
    193 c
    194 c
    195 c-----------------------------------------------------------------------
    196 c   40. transformation des variables dynamiques en variables physiques:
    197 c   ---------------------------------------------------------------
    198 
    199 c   41. pressions au sol (en Pascals)
    200 c   ----------------------------------
    201 
    202 
    203       zpsrf(1) = pps(1,1)
    204 
    205       ig0  = 2
    206       DO j = 2,jjm
    207          CALL SCOPY( iim,pps(1,j),1,zpsrf(ig0), 1 )
    208          ig0 = ig0+iim
    209       ENDDO
    210 
    211       zpsrf(ngridmx) = pps(1,jjp1)
    212 
    213 
    214 c   42. pression intercouches et fonction d'Exner:
    215 c
    216 c   -----------------------------------------------------------------
    217 c     .... zplev  definis aux (llm +1) interfaces des couches  ....
    218 c     .... zplay  definis aux (  llm )    milieux des couches  ....
    219 c   -----------------------------------------------------------------
    220 
    221 c    ...    Exner = cp * ( p(l) / preff ) ** kappa     ....
    222 c
    223        unskap   = 1./ kappa
    224 c
    225       DO l = 1, llm
    226         zpk(   1,l ) = ppk(1,1,l)
    227         zplev( 1,l ) = pp(1,1,l)
    228         ig0 = 2
    229           DO j = 2, jjm
    230              DO i =1, iim
    231               zpk(   ig0,l ) = ppk(i,j,l)
    232               zplev( ig0,l ) = pp(i,j,l)
    233               ig0 = ig0 +1
    234              ENDDO
    235           ENDDO
    236         zpk(   ngridmx,l ) = ppk(1,jjp1,l)
    237         zplev( ngridmx,l ) = pp(1,jjp1,l)
    238       ENDDO
    239         zplev( 1,llmp1 ) = pp(1,1,llmp1)
    240         ig0 = 2
    241           DO j = 2, jjm
    242              DO i =1, iim
    243               zplev( ig0,llmp1 ) = pp(i,j,llmp1)
    244               ig0 = ig0 +1
    245              ENDDO
    246           ENDDO
    247         zplev( ngridmx,llmp1 ) = pp(1,jjp1,llmp1)
    248 c
    249 c
    250 
    251 c   43. temperature naturelle (en K) et pressions milieux couches .
    252 c   ---------------------------------------------------------------
    253 
    254       DO l=1,llm
    255 
    256          pksurcp     =  ppk(1,1,l) / cpp
    257          zplay(1,l)  =  preff * pksurcp ** unskap
    258          ztfi(1,l)   =  pteta(1,1,l) *  pksurcp
    259          pcvgt(1,l)  =  pdteta(1,1,l) * pksurcp / pmasse(1,1,l)
    260          ig0         = 2
    261 
    262          DO j = 2, jjm
    263             DO i = 1, iim
    264               pksurcp        = ppk(i,j,l) / cpp
    265               zplay(ig0,l)   = preff * pksurcp ** unskap
    266               ztfi(ig0,l)    = pteta(i,j,l)  * pksurcp
    267               pcvgt(ig0,l)   = pdteta(i,j,l) * pksurcp / pmasse(i,j,l)
    268               ig0            = ig0 + 1
    269             ENDDO
    270          ENDDO
    271 
    272          pksurcp       = ppk(1,jjp1,l) / cpp
    273          zplay(ig0,l)  = preff * pksurcp ** unskap
    274          ztfi (ig0,l)  = pteta(1,jjp1,l)  * pksurcp
    275          pcvgt(ig0,l)  = pdteta(1,jjp1,l) * pksurcp/ pmasse(1,jjp1,l)
    276 
    277       ENDDO
    278 
    279 c   43.bis traceurs
    280 c   ---------------
    281 c
    282       itr=0
    283       DO iq=1,nqtot
    284          IF(.NOT.tracers(iq)%isAdvected) CYCLE
    285          itr = itr + 1
    286          DO l=1,llm
    287             zqfi(1,l,itr) = pq(1,1,l,iq)
    288             ig0           = 2
    289             DO j=2,jjm
    290                DO i = 1, iim
    291                   zqfi(ig0,l,itr) = pq(i,j,l,iq)
    292                   ig0             = ig0 + 1
    293                ENDDO
    294             ENDDO
    295             zqfi(ig0,l,itr) = pq(1,jjp1,l,iq)
     4!
     5!
     6SUBROUTINE calfis(lafin, &
     7        jD_cur, jH_cur, &
     8        pucov, &
     9        pvcov, &
     10        pteta, &
     11        pq, &
     12        pmasse, &
     13        pps, &
     14        pp, &
     15        ppk, &
     16        pphis, &
     17        pphi, &
     18        pducov, &
     19        pdvcov, &
     20        pdteta, &
     21        pdq, &
     22        flxw, &
     23        pdufi, &
     24        pdvfi, &
     25        pdhfi, &
     26        pdqfi, &
     27        pdpsfi)
     28  !
     29  !    Auteur :  P. Le Van, F. Hourdin
     30  !   .........
     31  USE infotrac, ONLY: nqtot, tracers
     32  USE control_mod, ONLY: planet_type, nsplit_phys
     33  USE callphysiq_mod, ONLY: call_physiq
     34  USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_PHYS
     35  USE comconst_mod, ONLY: cpp, daysec, dtphys, dtvr, kappa, pi
     36  USE comvert_mod, ONLY: preff, presnivs
     37
     38  IMPLICIT NONE
     39  !=======================================================================
     40  !
     41  !   1. rearrangement des tableaux et transformation
     42  !  variables dynamiques  >  variables physiques
     43  !   2. calcul des termes physiques
     44  !   3. retransformation des tendances physiques en tendances dynamiques
     45  !
     46  !   remarques:
     47  !   ----------
     48  !
     49  !    - les vents sont donnes dans la physique par leurs composantes
     50  !  naturelles.
     51  !    - la variable thermodynamique de la physique est une variable
     52  !  intensive :   T
     53  !  pour la dynamique on prend    T * ( preff / p(l) ) **kappa
     54  !    - les deux seules variables dependant de la geometrie necessaires
     55  !  pour la physique sont la latitude pour le rayonnement et
     56  !  l'aire de la maille quand on veut integrer une grandeur
     57  !  horizontalement.
     58  !    - les points de la physique sont les points scalaires de la
     59  !  la dynamique; numerotation:
     60  !      1 pour le pole nord
     61  !      (jjm-1)*iim pour l'interieur du domaine
     62  !      ngridmx pour le pole sud
     63  !  ---> ngridmx=2+(jjm-1)*iim
     64  !
     65  ! Input :
     66  ! -------
     67  !   pucov           covariant zonal velocity
     68  !   pvcov           covariant meridional velocity
     69  !   pteta           potential temperature
     70  !   pps             surface pressure
     71  !   pmasse          masse d'air dans chaque maille
     72  !   pts             surface temperature  (K)
     73  !   callrad         clef d'appel au rayonnement
     74  !
     75  !    Output :
     76  !    --------
     77  !    pdufi          tendency for the natural zonal velocity (ms-1)
     78  !    pdvfi          tendency for the natural meridional velocity
     79  !    pdhfi          tendency for the potential temperature
     80  !    pdtsfi         tendency for the surface temperature
     81  !
     82  !    pdtrad         radiative tendencies  \  both input
     83  !    pfluxrad       radiative fluxes      /  and output
     84  !
     85  !=======================================================================
     86  !
     87  !-----------------------------------------------------------------------
     88  !
     89  !    0.  Declarations :
     90  !    ------------------
     91
     92  include "dimensions.h"
     93  include "paramet.h"
     94
     95  INTEGER :: ngridmx
     96  PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
     97
     98  include "comgeom2.h"
     99  include "iniprint.h"
     100
     101  !    Arguments :
     102  !    -----------
     103  LOGICAL,INTENT(IN) ::  lafin ! .TRUE. for the very last CALL to physics
     104  REAL,INTENT(IN):: jD_cur, jH_cur
     105  REAL,INTENT(IN) :: pvcov(iip1,jjm,llm) ! covariant meridional velocity
     106  REAL,INTENT(IN) :: pucov(iip1,jjp1,llm) ! covariant zonal velocity
     107  REAL,INTENT(IN) :: pteta(iip1,jjp1,llm) ! potential temperature
     108  REAL,INTENT(IN) :: pmasse(iip1,jjp1,llm) ! mass in each cell ! not used
     109  REAL,INTENT(IN) :: pq(iip1,jjp1,llm,nqtot) ! tracers
     110  REAL,INTENT(IN) :: pphis(iip1,jjp1) ! surface geopotential
     111  REAL,INTENT(IN) :: pphi(iip1,jjp1,llm) ! geopotential
     112
     113  REAL,INTENT(IN) :: pdvcov(iip1,jjm,llm) ! dynamical tendency on vcov
     114  REAL,INTENT(IN) :: pducov(iip1,jjp1,llm) ! dynamical tendency on ucov
     115  REAL,INTENT(IN) :: pdteta(iip1,jjp1,llm) ! dynamical tendency on teta
     116  ! ! NB: pdteta is used only to compute pcvgt which is in fact not used...
     117  REAL,INTENT(IN) :: pdq(iip1,jjp1,llm,nqtot) ! dynamical tendency on tracers
     118  ! ! NB: pdq is only used to compute pcvgq which is in fact not used...
     119
     120  REAL,INTENT(IN) :: pps(iip1,jjp1) ! surface pressure (Pa)
     121  REAL,INTENT(IN) :: pp(iip1,jjp1,llmp1) ! pressure at mesh interfaces (Pa)
     122  REAL,INTENT(IN) :: ppk(iip1,jjp1,llm) ! Exner at mid-layer
     123  REAL,INTENT(IN) :: flxw(iip1,jjp1,llm) ! Vertical mass flux on lower mesh interfaces (kg/s) (on llm because flxw(:,:,llm+1)=0)
     124
     125  ! ! tendencies (in */s) from the physics
     126  REAL,INTENT(OUT) :: pdvfi(iip1,jjm,llm) ! tendency on covariant meridional wind
     127  REAL,INTENT(OUT) :: pdufi(iip1,jjp1,llm) ! tendency on covariant zonal wind
     128  REAL,INTENT(OUT) :: pdhfi(iip1,jjp1,llm) ! tendency on potential temperature (K/s)
     129  REAL,INTENT(OUT) :: pdqfi(iip1,jjp1,llm,nqtot) ! tendency on tracers
     130  REAL,INTENT(OUT) :: pdpsfi(iip1,jjp1) ! tendency on surface pressure (Pa/s)
     131
     132
     133  !    Local variables :
     134  !    -----------------
     135
     136  INTEGER :: i,j,l,ig0,ig,iq,itr
     137  REAL :: zpsrf(ngridmx)
     138  REAL :: zplev(ngridmx,llm+1),zplay(ngridmx,llm)
     139  REAL :: zphi(ngridmx,llm),zphis(ngridmx)
     140  !
     141  REAL :: zrot(iip1,jjm,llm) ! AdlC May 2014
     142  REAL :: zufi(ngridmx,llm), zvfi(ngridmx,llm)
     143  REAL :: zrfi(ngridmx,llm) ! relative wind vorticity
     144  REAL :: ztfi(ngridmx,llm),zqfi(ngridmx,llm,nqtot)
     145  REAL :: zpk(ngridmx,llm)
     146  !
     147  REAL :: pcvgu(ngridmx,llm), pcvgv(ngridmx,llm)
     148  REAL :: pcvgt(ngridmx,llm), pcvgq(ngridmx,llm,2)
     149  !
     150  REAL :: zdufi(ngridmx,llm),zdvfi(ngridmx,llm)
     151  REAL :: zdtfi(ngridmx,llm),zdqfi(ngridmx,llm,nqtot)
     152  REAL :: zdpsrf(ngridmx)
     153  !
     154  REAL :: zdufic(ngridmx,llm),zdvfic(ngridmx,llm)
     155  REAL :: zdtfic(ngridmx,llm),zdqfic(ngridmx,llm,nqtot)
     156  REAL :: jH_cur_split,zdt_split
     157  LOGICAL :: debut_split,lafin_split
     158  INTEGER :: isplit
     159
     160  REAL :: zsin(iim),zcos(iim),z1(iim)
     161  REAL :: zsinbis(iim),zcosbis(iim),z1bis(iim)
     162  REAL :: unskap, pksurcp
     163  !
     164  REAL :: flxwfi(ngridmx,llm)  ! Flux de masse verticale sur la grille physiq
     165  !
     166
     167  REAL :: SSUM
     168
     169  LOGICAL,SAVE :: firstcal=.TRUE., debut=.TRUE.
     170   ! REAL rdayvrai
     171
     172  !
     173  !-----------------------------------------------------------------------
     174  !
     175  !    1. Initialisations :
     176  !    --------------------
     177  !
     178  !
     179  IF ( firstcal )  THEN
     180    debut = .TRUE.
     181    IF (ngridmx/=2+(jjm-1)*iim) THEN
     182     write(lunout,*) 'STOP dans calfis'
     183     write(lunout,*) &
     184           'La dimension ngridmx doit etre egale a 2 + (jjm-1)*iim'
     185     write(lunout,*) '  ngridmx  jjm   iim   '
     186     write(lunout,*) ngridmx,jjm,iim
     187     CALL abort_gcm("calfis", "", 1)
     188    ENDIF
     189  ELSE
     190    debut = .FALSE.
     191  ENDIF ! of IF (firstcal)
     192
     193  !
     194  !
     195  !-----------------------------------------------------------------------
     196  !   40. transformation des variables dynamiques en variables physiques:
     197  !   ---------------------------------------------------------------
     198
     199  !   41. pressions au sol (en Pascals)
     200  !   ----------------------------------
     201
     202
     203  zpsrf(1) = pps(1,1)
     204
     205  ig0  = 2
     206  DO j = 2,jjm
     207     CALL SCOPY( iim,pps(1,j),1,zpsrf(ig0), 1 )
     208     ig0 = ig0+iim
     209  ENDDO
     210
     211  zpsrf(ngridmx) = pps(1,jjp1)
     212
     213
     214  !   42. pression intercouches et fonction d'Exner:
     215  !
     216  !   -----------------------------------------------------------------
     217  ! .... zplev  definis aux (llm +1) interfaces des couches  ....
     218  ! .... zplay  definis aux (  llm )    milieux des couches  ....
     219  !   -----------------------------------------------------------------
     220
     221  !    ...    Exner = cp * ( p(l) / preff ) ** kappa     ....
     222  !
     223   unskap   = 1./ kappa
     224  !
     225  DO l = 1, llm
     226    zpk(   1,l ) = ppk(1,1,l)
     227    zplev( 1,l ) = pp(1,1,l)
     228    ig0 = 2
     229      DO j = 2, jjm
     230         DO i =1, iim
     231          zpk(   ig0,l ) = ppk(i,j,l)
     232          zplev( ig0,l ) = pp(i,j,l)
     233          ig0 = ig0 +1
    296234         ENDDO
    297235      ENDDO
    298 
    299 c   convergence dynamique pour les traceurs "EAU"
    300 ! Earth-specific treatment of first 2 tracers (water)
    301        if (planet_type=="earth") then
    302         DO iq=1,2
    303          DO l=1,llm
    304             pcvgq(1,l,iq)= pdq(1,1,l,iq) / pmasse(1,1,l)
    305             ig0          = 2
    306             DO j=2,jjm
    307                DO i = 1, iim
    308                   pcvgq(ig0,l,iq) = pdq(i,j,l,iq) / pmasse(i,j,l)
    309                   ig0             = ig0 + 1
    310                ENDDO
    311             ENDDO
    312             pcvgq(ig0,l,iq)= pdq(1,jjp1,l,iq) / pmasse(1,jjp1,l)
    313          ENDDO
    314         ENDDO
    315        endif ! of if (planet_type=="earth")
    316 
    317 
    318 c   Geopotentiel calcule par rapport a la surface locale:
    319 c   -----------------------------------------------------
    320 
    321       CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,pphi,zphi)
    322       CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,pphis,zphis)
    323       DO l=1,llm
    324          DO ig=1,ngridmx
    325            zphi(ig,l)=zphi(ig,l)-zphis(ig)
     236    zpk(   ngridmx,l ) = ppk(1,jjp1,l)
     237    zplev( ngridmx,l ) = pp(1,jjp1,l)
     238  ENDDO
     239    zplev( 1,llmp1 ) = pp(1,1,llmp1)
     240    ig0 = 2
     241      DO j = 2, jjm
     242         DO i =1, iim
     243          zplev( ig0,llmp1 ) = pp(i,j,llmp1)
     244          ig0 = ig0 +1
    326245         ENDDO
    327246      ENDDO
    328 
    329 c   ....  Calcul de la vitesse  verticale  ( en Pa*m*s  ou Kg/s )  ....
    330 c JG : ancien calcule de omega utilise dans physiq.F. Maintenant le flux
    331 c    de masse est calclue dans advtrac.F
    332 c      DO l=1,llm
    333 c        pvervel(1,l)=pw(1,1,l) * g /apoln
    334 c        ig0=2
    335 c       DO j=2,jjm
    336 c           DO i = 1, iim
    337 c              pvervel(ig0,l) = pw(i,j,l) * g * unsaire(i,j)
    338 c              ig0 = ig0 + 1
    339 c           ENDDO
    340 c       ENDDO
    341 c        pvervel(ig0,l)=pw(1,jjp1,l) * g /apols
    342 c      ENDDO
    343 
    344 c
    345 c   45. champ u:
    346 c   ------------
    347 
    348       DO l=1,llm
    349 
    350          DO j=2,jjm
    351             ig0 = 1+(j-2)*iim
    352             zufi(ig0+1,l)= 0.5 *
    353      $      ( pucov(iim,j,l)/cu(iim,j) + pucov(1,j,l)/cu(1,j) )
    354             pcvgu(ig0+1,l)= 0.5 *
    355      $      ( pducov(iim,j,l)/cu(iim,j) + pducov(1,j,l)/cu(1,j) )
    356             DO i=2,iim
    357                zufi(ig0+i,l)= 0.5 *
    358      $         ( pucov(i-1,j,l)/cu(i-1,j) + pucov(i,j,l)/cu(i,j) )
    359                pcvgu(ig0+i,l)= 0.5 *
    360      $         ( pducov(i-1,j,l)/cu(i-1,j) + pducov(i,j,l)/cu(i,j) )
    361       END DO
    362       END DO
    363 
    364       END DO
    365 
    366 
    367 C  Alvaro de la Camara (May 2014)
    368 C  46.1 Calcul de la vorticite et passage sur la grille physique
    369 C  --------------------------------------------------------------
    370       DO l=1,llm
    371         do i=1,iim
    372           do j=1,jjm
    373             zrot(i,j,l) = (pvcov(i+1,j,l) - pvcov(i,j,l)
    374      $                   + pucov(i,j+1,l) - pucov(i,j,l))
    375      $                   / (cu(i,j)+cu(i,j+1))
    376      $                   / (cv(i+1,j)+cv(i,j)) *4
    377           enddo
    378         enddo
     247    zplev( ngridmx,llmp1 ) = pp(1,jjp1,llmp1)
     248  !
     249  !
     250
     251  !   43. temperature naturelle (en K) et pressions milieux couches .
     252  !   ---------------------------------------------------------------
     253
     254  DO l=1,llm
     255
     256     pksurcp     =  ppk(1,1,l) / cpp
     257     zplay(1,l)  =  preff * pksurcp ** unskap
     258     ztfi(1,l)   =  pteta(1,1,l) *  pksurcp
     259     pcvgt(1,l)  =  pdteta(1,1,l) * pksurcp / pmasse(1,1,l)
     260     ig0         = 2
     261
     262     DO j = 2, jjm
     263        DO i = 1, iim
     264          pksurcp        = ppk(i,j,l) / cpp
     265          zplay(ig0,l)   = preff * pksurcp ** unskap
     266          ztfi(ig0,l)    = pteta(i,j,l)  * pksurcp
     267          pcvgt(ig0,l)   = pdteta(i,j,l) * pksurcp / pmasse(i,j,l)
     268          ig0            = ig0 + 1
     269        ENDDO
     270     ENDDO
     271
     272     pksurcp       = ppk(1,jjp1,l) / cpp
     273     zplay(ig0,l)  = preff * pksurcp ** unskap
     274     ztfi (ig0,l)  = pteta(1,jjp1,l)  * pksurcp
     275     pcvgt(ig0,l)  = pdteta(1,jjp1,l) * pksurcp/ pmasse(1,jjp1,l)
     276
     277  ENDDO
     278
     279  !   43.bis traceurs
     280  !   ---------------
     281  !
     282  itr=0
     283  DO iq=1,nqtot
     284     IF(.NOT.tracers(iq)%isAdvected) CYCLE
     285     itr = itr + 1
     286     DO l=1,llm
     287        zqfi(1,l,itr) = pq(1,1,l,iq)
     288        ig0           = 2
     289        DO j=2,jjm
     290           DO i = 1, iim
     291              zqfi(ig0,l,itr) = pq(i,j,l,iq)
     292              ig0             = ig0 + 1
     293           ENDDO
     294        ENDDO
     295        zqfi(ig0,l,itr) = pq(1,jjp1,l,iq)
     296     ENDDO
     297  ENDDO
     298
     299  !   convergence dynamique pour les traceurs "EAU"
     300  ! Earth-specific treatment of first 2 tracers (water)
     301   if (planet_type=="earth") then
     302    DO iq=1,2
     303     DO l=1,llm
     304        pcvgq(1,l,iq)= pdq(1,1,l,iq) / pmasse(1,1,l)
     305        ig0          = 2
     306        DO j=2,jjm
     307           DO i = 1, iim
     308              pcvgq(ig0,l,iq) = pdq(i,j,l,iq) / pmasse(i,j,l)
     309              ig0             = ig0 + 1
     310           ENDDO
     311        ENDDO
     312        pcvgq(ig0,l,iq)= pdq(1,jjp1,l,iq) / pmasse(1,jjp1,l)
     313     ENDDO
     314    ENDDO
     315   endif ! of if (planet_type=="earth")
     316
     317
     318  !   Geopotentiel calcule par rapport a la surface locale:
     319  !   -----------------------------------------------------
     320
     321  CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,pphi,zphi)
     322  CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,pphis,zphis)
     323  DO l=1,llm
     324     DO ig=1,ngridmx
     325       zphi(ig,l)=zphi(ig,l)-zphis(ig)
     326     ENDDO
     327  ENDDO
     328
     329  !   ....  Calcul de la vitesse  verticale  ( en Pa*m*s  ou Kg/s )  ....
     330  ! JG : ancien calcule de omega utilise dans physiq.F. Maintenant le flux
     331  !    de masse est calclue dans advtrac.F
     332   ! DO l=1,llm
     333   !   pvervel(1,l)=pw(1,1,l) * g /apoln
     334   !   ig0=2
     335   !  DO j=2,jjm
     336   !      DO i = 1, iim
     337   !         pvervel(ig0,l) = pw(i,j,l) * g * unsaire(i,j)
     338   !         ig0 = ig0 + 1
     339   !      ENDDO
     340  !       ENDDO
     341   !   pvervel(ig0,l)=pw(1,jjp1,l) * g /apols
     342   ! ENDDO
     343
     344  !
     345  !   45. champ u:
     346  !   ------------
     347
     348  DO l=1,llm
     349
     350     DO j=2,jjm
     351        ig0 = 1+(j-2)*iim
     352        zufi(ig0+1,l)= 0.5 * &
     353              ( pucov(iim,j,l)/cu(iim,j) + pucov(1,j,l)/cu(1,j) )
     354        pcvgu(ig0+1,l)= 0.5 * &
     355              ( pducov(iim,j,l)/cu(iim,j) + pducov(1,j,l)/cu(1,j) )
     356        DO i=2,iim
     357           zufi(ig0+i,l)= 0.5 * &
     358                 ( pucov(i-1,j,l)/cu(i-1,j) + pucov(i,j,l)/cu(i,j) )
     359           pcvgu(ig0+i,l)= 0.5 * &
     360                 ( pducov(i-1,j,l)/cu(i-1,j) + pducov(i,j,l)/cu(i,j) )
     361  END DO
     362  END DO
     363
     364  END DO
     365
     366
     367  !  Alvaro de la Camara (May 2014)
     368  !  46.1 Calcul de la vorticite et passage sur la grille physique
     369  !  --------------------------------------------------------------
     370  DO l=1,llm
     371    do i=1,iim
     372      do j=1,jjm
     373        zrot(i,j,l) = (pvcov(i+1,j,l) - pvcov(i,j,l) &
     374              + pucov(i,j+1,l) - pucov(i,j,l)) &
     375              / (cu(i,j)+cu(i,j+1)) &
     376              / (cv(i+1,j)+cv(i,j)) *4
     377      enddo
     378    enddo
     379  ENDDO
     380
     381  !   46.champ v:
     382  !   -----------
     383
     384  DO l=1,llm
     385     DO j=2,jjm
     386        ig0=1+(j-2)*iim
     387        DO i=1,iim
     388           zvfi(ig0+i,l)= 0.5 * &
     389                 ( pvcov(i,j-1,l)/cv(i,j-1) + pvcov(i,j,l)/cv(i,j) )
     390           pcvgv(ig0+i,l)= 0.5 * &
     391                 ( pdvcov(i,j-1,l)/cv(i,j-1) + pdvcov(i,j,l)/cv(i,j) )
     392        ENDDO
     393           zrfi(ig0 + 1,l)= 0.25 *(zrot(iim,j-1,l)+zrot(iim,j,l) &
     394                 +zrot(1,j-1,l)+zrot(1,j,l))
     395        DO i=2,iim
     396           zrfi(ig0 + i,l)= 0.25 *(zrot(i-1,j-1,l)+zrot(i-1,j,l) &
     397                 +zrot(i,j-1,l)+zrot(i,j,l))   !  AdlC MAY 2014
     398        ENDDO
     399     ENDDO
     400  ENDDO
     401
     402
     403  !   47. champs de vents aux pole nord
     404  !   ------------------------------
     405     ! U = 1 / pi  *  integrale [ v * cos(long) * d long ]
     406     ! V = 1 / pi  *  integrale [ v * sin(long) * d long ]
     407
     408  DO l=1,llm
     409
     410     z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,1,l)/cv(1,1)
     411     z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,1,l)/cv(1,1)
     412     DO i=2,iim
     413        z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
     414        z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,1,l)/cv(i,1)
     415     ENDDO
     416
     417     DO i=1,iim
     418        zcos(i)   = COS(rlonv(i))*z1(i)
     419        zcosbis(i)= COS(rlonv(i))*z1bis(i)
     420        zsin(i)   = SIN(rlonv(i))*z1(i)
     421        zsinbis(i)= SIN(rlonv(i))*z1bis(i)
     422     ENDDO
     423
     424     zufi(1,l)  = SSUM(iim,zcos,1)/pi
     425     pcvgu(1,l) = SSUM(iim,zcosbis,1)/pi
     426     zvfi(1,l)  = SSUM(iim,zsin,1)/pi
     427     pcvgv(1,l) = SSUM(iim,zsinbis,1)/pi
     428     zrfi(1, l) = 0.
     429  ENDDO
     430
     431
     432  !   48. champs de vents aux pole sud:
     433  !   ---------------------------------
     434     ! U = 1 / pi  *  integrale [ v * cos(long) * d long ]
     435     ! V = 1 / pi  *  integrale [ v * sin(long) * d long ]
     436
     437  DO l=1,llm
     438
     439     z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,jjm,l)/cv(1,jjm)
     440     z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,jjm,l)/cv(1,jjm)
     441     DO i=2,iim
     442        z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
     443        z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,jjm,l)/cv(i,jjm)
     444     ENDDO
     445
     446     DO i=1,iim
     447        zcos(i)    = COS(rlonv(i))*z1(i)
     448        zcosbis(i) = COS(rlonv(i))*z1bis(i)
     449        zsin(i)    = SIN(rlonv(i))*z1(i)
     450        zsinbis(i) = SIN(rlonv(i))*z1bis(i)
     451     ENDDO
     452
     453     zufi(ngridmx,l)  = SSUM(iim,zcos,1)/pi
     454     pcvgu(ngridmx,l) = SSUM(iim,zcosbis,1)/pi
     455     zvfi(ngridmx,l)  = SSUM(iim,zsin,1)/pi
     456     pcvgv(ngridmx,l) = SSUM(iim,zsinbis,1)/pi
     457     zrfi(ngridmx, l) = 0.
     458  ENDDO
     459  !
     460  ! On change de grille, dynamique vers physiq, pour le flux de masse verticale
     461  CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,flxw,flxwfi)
     462
     463  !-----------------------------------------------------------------------
     464  !   Appel de la physique:
     465  !   ---------------------
     466
     467
     468
     469   ! write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
     470  zdt_split=dtphys/nsplit_phys
     471  zdufic(:,:)=0.
     472  zdvfic(:,:)=0.
     473  zdtfic(:,:)=0.
     474  zdqfic(:,:,:)=0.
     475
     476   IF (CPPKEY_PHYS) THEN
     477
     478   do isplit=1,nsplit_phys
     479
     480     jH_cur_split=jH_cur+(isplit-1) * dtvr / (daysec *nsplit_phys)
     481     debut_split=debut.and.isplit==1
     482     lafin_split=lafin.and.isplit==nsplit_phys
     483
     484   ! if (planet_type=="earth") then
     485    CALL call_physiq(ngridmx,llm,nqtot,tracers(:)%name, &
     486          debut_split,lafin_split, &
     487          jD_cur,jH_cur_split,zdt_split, &
     488          zplev,zplay, &
     489          zpk,zphi,zphis, &
     490          presnivs, &
     491          zufi,zvfi,zrfi,ztfi,zqfi, &
     492          flxwfi,pducov, &
     493          zdufi,zdvfi,zdtfi,zdqfi,zdpsrf)
     494
     495   ! else if ( planet_type=="generic" ) then
     496
     497   !    CALL physiq (ngridmx,     !! ngrid
     498  ! .             llm,            !! nlayer
     499  ! .             nqtot,          !! nq
     500  ! .             tracers(:)%name,!! tracer names from dynamical core (given in infotrac)
     501  ! .             debut_split,    !! firstcall
     502  ! .             lafin_split,    !! lastcall
     503  ! .             jD_cur,         !! pday. see leapfrog
     504  ! .             jH_cur_split,   !! ptime "fraction of day"
     505  ! .             zdt_split,      !! ptimestep
     506  ! .             zplev,          !! pplev
     507  ! .             zplay,          !! pplay
     508  ! .             zphi,           !! pphi
     509  ! .             zufi,           !! pu
     510  ! .             zvfi,           !! pv
     511  ! .             ztfi,           !! pt
     512  ! .             zqfi,           !! pq
     513  ! .             flxwfi,         !! pw !! or 0. anyway this is for diagnostic. not used in physiq.
     514  ! .             zdufi,          !! pdu
     515  ! .             zdvfi,          !! pdv
     516  ! .             zdtfi,          !! pdt
     517  ! .             zdqfi,          !! pdq
     518  ! .             zdpsrf,         !! pdpsrf
     519  ! .             tracerdyn)      !! tracerdyn <-- utilite ???
     520
     521  !  endif ! of if (planet_type=="earth")
     522
     523     zufi(:,:)=zufi(:,:)+zdufi(:,:)*zdt_split
     524     zvfi(:,:)=zvfi(:,:)+zdvfi(:,:)*zdt_split
     525     ztfi(:,:)=ztfi(:,:)+zdtfi(:,:)*zdt_split
     526     zqfi(:,:,:)=zqfi(:,:,:)+zdqfi(:,:,:)*zdt_split
     527
     528     zdufic(:,:)=zdufic(:,:)+zdufi(:,:)
     529     zdvfic(:,:)=zdvfic(:,:)+zdvfi(:,:)
     530     zdtfic(:,:)=zdtfic(:,:)+zdtfi(:,:)
     531     zdqfic(:,:,:)=zdqfic(:,:,:)+zdqfi(:,:,:)
     532
     533   enddo ! of do isplit=1,nsplit_phys
     534
     535   END IF
     536
     537  zdufi(:,:)=zdufic(:,:)/nsplit_phys
     538  zdvfi(:,:)=zdvfic(:,:)/nsplit_phys
     539  zdtfi(:,:)=zdtfic(:,:)/nsplit_phys
     540  zdqfi(:,:,:)=zdqfic(:,:,:)/nsplit_phys
     541
     542  !-----------------------------------------------------------------------
     543  !   transformation des tendances physiques en tendances dynamiques:
     544  !   ---------------------------------------------------------------
     545
     546  !  tendance sur la pression :
     547  !  -----------------------------------
     548
     549  CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,zdpsrf,pdpsfi)
     550  !
     551  !   62. enthalpie potentielle
     552  !   ---------------------
     553
     554  DO l=1,llm
     555
     556     DO i=1,iip1
     557      pdhfi(i,1,l)    = cpp *  zdtfi(1,l)      / ppk(i, 1  ,l)
     558      pdhfi(i,jjp1,l) = cpp *  zdtfi(ngridmx,l)/ ppk(i,jjp1,l)
     559     ENDDO
     560
     561     DO j=2,jjm
     562        ig0=1+(j-2)*iim
     563        DO i=1,iim
     564           pdhfi(i,j,l) = cpp * zdtfi(ig0+i,l) / ppk(i,j,l)
     565        ENDDO
     566           pdhfi(iip1,j,l) =  pdhfi(1,j,l)
     567     ENDDO
     568
     569  ENDDO
     570
     571
     572  !   62. humidite specifique
     573  !   ---------------------
     574  ! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
     575   ! DO iq=1,nqtot
     576   !    DO l=1,llm
     577   !       DO i=1,iip1
     578   !          pdqfi(i,1,l,iq)    = zdqfi(1,l,iq)
     579   !          pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,iq)
     580   !       ENDDO
     581   !       DO j=2,jjm
     582   !          ig0=1+(j-2)*iim
     583   !          DO i=1,iim
     584   !             pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,iq)
     585   !          ENDDO
     586   !          pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,iq)
     587   !       ENDDO
     588   !    ENDDO
     589   ! ENDDO
     590
     591  !   63. traceurs
     592  !   ------------
     593  ! initialisation des tendances
     594  pdqfi(:,:,:,:)=0.
     595  !
     596  itr = 0
     597  DO iq=1,nqtot
     598     IF(.NOT.tracers(iq)%isAdvected) CYCLE
     599     itr = itr + 1
     600     DO l=1,llm
     601        DO i=1,iip1
     602           pdqfi(i,1,l,iq)    = zdqfi(1,l,itr)
     603           pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,itr)
     604        ENDDO
     605        DO j=2,jjm
     606           ig0=1+(j-2)*iim
     607           DO i=1,iim
     608              pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,itr)
     609           ENDDO
     610           pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,itr)
     611        ENDDO
     612     ENDDO
     613  ENDDO
     614
     615  !   65. champ u:
     616  !   ------------
     617
     618  DO l=1,llm
     619
     620     DO i=1,iip1
     621        pdufi(i,1,l)    = 0.
     622        pdufi(i,jjp1,l) = 0.
     623     ENDDO
     624
     625     DO j=2,jjm
     626        ig0=1+(j-2)*iim
     627        DO i=1,iim-1
     628           pdufi(i,j,l)= &
     629                 0.5*(zdufi(ig0+i,l)+zdufi(ig0+i+1,l))*cu(i,j)
     630        ENDDO
     631        pdufi(iim,j,l)= &
     632              0.5*(zdufi(ig0+1,l)+zdufi(ig0+iim,l))*cu(iim,j)
     633        pdufi(iip1,j,l)=pdufi(1,j,l)
     634     ENDDO
     635
     636  ENDDO
     637
     638
     639  !   67. champ v:
     640  !   ------------
     641
     642  DO l=1,llm
     643
     644     DO j=2,jjm-1
     645        ig0=1+(j-2)*iim
     646        DO i=1,iim
     647           pdvfi(i,j,l)= &
     648                 0.5*(zdvfi(ig0+i,l)+zdvfi(ig0+i+iim,l))*cv(i,j)
     649        ENDDO
     650        pdvfi(iip1,j,l) = pdvfi(1,j,l)
     651     ENDDO
     652  ENDDO
     653
     654
     655  !   68. champ v pres des poles:
     656  !   ---------------------------
     657   ! v = U * cos(long) + V * SIN(long)
     658
     659  DO l=1,llm
     660
     661     DO i=1,iim
     662        pdvfi(i,1,l)= &
     663              zdufi(1,l)*COS(rlonv(i))+zdvfi(1,l)*SIN(rlonv(i))
     664        pdvfi(i,jjm,l)=zdufi(ngridmx,l)*COS(rlonv(i)) &
     665              +zdvfi(ngridmx,l)*SIN(rlonv(i))
     666        pdvfi(i,1,l)= &
     667              0.5*(pdvfi(i,1,l)+zdvfi(i+1,l))*cv(i,1)
     668        pdvfi(i,jjm,l)= &
     669              0.5*(pdvfi(i,jjm,l)+zdvfi(ngridmx-iip1+i,l))*cv(i,jjm)
    379670      ENDDO
    380671
    381 c   46.champ v:
    382 c   -----------
    383 
    384       DO l=1,llm
    385          DO j=2,jjm
    386             ig0=1+(j-2)*iim
    387             DO i=1,iim
    388                zvfi(ig0+i,l)= 0.5 *
    389      $         ( pvcov(i,j-1,l)/cv(i,j-1) + pvcov(i,j,l)/cv(i,j) )
    390                pcvgv(ig0+i,l)= 0.5 *
    391      $         ( pdvcov(i,j-1,l)/cv(i,j-1) + pdvcov(i,j,l)/cv(i,j) )
    392             ENDDO
    393                zrfi(ig0 + 1,l)= 0.25 *(zrot(iim,j-1,l)+zrot(iim,j,l)
    394      &                                +zrot(1,j-1,l)+zrot(1,j,l))
    395             DO i=2,iim
    396                zrfi(ig0 + i,l)= 0.25 *(zrot(i-1,j-1,l)+zrot(i-1,j,l)
    397      $                   +zrot(i,j-1,l)+zrot(i,j,l))   !  AdlC MAY 2014
    398             ENDDO
    399          ENDDO
    400       ENDDO
    401 
    402 
    403 c   47. champs de vents aux pole nord
    404 c   ------------------------------
    405 c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
    406 c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]
    407 
    408       DO l=1,llm
    409 
    410          z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,1,l)/cv(1,1)
    411          z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,1,l)/cv(1,1)
    412          DO i=2,iim
    413             z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
    414             z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,1,l)/cv(i,1)
    415          ENDDO
    416 
    417          DO i=1,iim
    418             zcos(i)   = COS(rlonv(i))*z1(i)
    419             zcosbis(i)= COS(rlonv(i))*z1bis(i)
    420             zsin(i)   = SIN(rlonv(i))*z1(i)
    421             zsinbis(i)= SIN(rlonv(i))*z1bis(i)
    422          ENDDO
    423 
    424          zufi(1,l)  = SSUM(iim,zcos,1)/pi
    425          pcvgu(1,l) = SSUM(iim,zcosbis,1)/pi
    426          zvfi(1,l)  = SSUM(iim,zsin,1)/pi
    427          pcvgv(1,l) = SSUM(iim,zsinbis,1)/pi
    428          zrfi(1, l) = 0.
    429       ENDDO
    430 
    431 
    432 c   48. champs de vents aux pole sud:
    433 c   ---------------------------------
    434 c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
    435 c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]
    436 
    437       DO l=1,llm
    438 
    439          z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,jjm,l)/cv(1,jjm)
    440          z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,jjm,l)/cv(1,jjm)
    441          DO i=2,iim
    442             z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
    443             z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,jjm,l)/cv(i,jjm)
    444          ENDDO
    445 
    446          DO i=1,iim
    447             zcos(i)    = COS(rlonv(i))*z1(i)
    448             zcosbis(i) = COS(rlonv(i))*z1bis(i)
    449             zsin(i)    = SIN(rlonv(i))*z1(i)
    450             zsinbis(i) = SIN(rlonv(i))*z1bis(i)
    451          ENDDO
    452 
    453          zufi(ngridmx,l)  = SSUM(iim,zcos,1)/pi
    454          pcvgu(ngridmx,l) = SSUM(iim,zcosbis,1)/pi
    455          zvfi(ngridmx,l)  = SSUM(iim,zsin,1)/pi
    456          pcvgv(ngridmx,l) = SSUM(iim,zsinbis,1)/pi
    457          zrfi(ngridmx, l) = 0.
    458       ENDDO
    459 c
    460 c On change de grille, dynamique vers physiq, pour le flux de masse verticale
    461       CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,flxw,flxwfi)
    462 
    463 c-----------------------------------------------------------------------
    464 c   Appel de la physique:
    465 c   ---------------------
    466 
    467 
    468 
    469 !      write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
    470       zdt_split=dtphys/nsplit_phys
    471       zdufic(:,:)=0.
    472       zdvfic(:,:)=0.
    473       zdtfic(:,:)=0.
    474       zdqfic(:,:,:)=0.
    475 
    476        IF (CPPKEY_PHYS) THEN
    477 
    478        do isplit=1,nsplit_phys
    479 
    480          jH_cur_split=jH_cur+(isplit-1) * dtvr / (daysec *nsplit_phys)
    481          debut_split=debut.and.isplit==1
    482          lafin_split=lafin.and.isplit==nsplit_phys
    483 
    484 !      if (planet_type=="earth") then
    485         CALL call_physiq(ngridmx,llm,nqtot,tracers(:)%name,
    486      &                   debut_split,lafin_split,
    487      &                   jD_cur,jH_cur_split,zdt_split,
    488      &                   zplev,zplay,
    489      &                   zpk,zphi,zphis,
    490      &                   presnivs,
    491      &                   zufi,zvfi,zrfi,ztfi,zqfi,
    492      &                   flxwfi,pducov,
    493      &                   zdufi,zdvfi,zdtfi,zdqfi,zdpsrf)
    494 
    495 !      else if ( planet_type=="generic" ) then
    496 
    497 !         CALL physiq (ngridmx,     !! ngrid
    498 !     .             llm,            !! nlayer
    499 !     .             nqtot,          !! nq
    500 !     .             tracers(:)%name,!! tracer names from dynamical core (given in infotrac)
    501 !     .             debut_split,    !! firstcall
    502 !     .             lafin_split,    !! lastcall
    503 !     .             jD_cur,         !! pday. see leapfrog
    504 !     .             jH_cur_split,   !! ptime "fraction of day"
    505 !     .             zdt_split,      !! ptimestep
    506 !     .             zplev,          !! pplev
    507 !     .             zplay,          !! pplay
    508 !     .             zphi,           !! pphi
    509 !     .             zufi,           !! pu
    510 !     .             zvfi,           !! pv
    511 !     .             ztfi,           !! pt
    512 !     .             zqfi,           !! pq
    513 !     .             flxwfi,         !! pw !! or 0. anyway this is for diagnostic. not used in physiq.
    514 !     .             zdufi,          !! pdu
    515 !     .             zdvfi,          !! pdv
    516 !     .             zdtfi,          !! pdt
    517 !     .             zdqfi,          !! pdq
    518 !     .             zdpsrf,         !! pdpsrf
    519 !     .             tracerdyn)      !! tracerdyn <-- utilite ???
    520 
    521 !      endif ! of if (planet_type=="earth")
    522 
    523          zufi(:,:)=zufi(:,:)+zdufi(:,:)*zdt_split
    524          zvfi(:,:)=zvfi(:,:)+zdvfi(:,:)*zdt_split
    525          ztfi(:,:)=ztfi(:,:)+zdtfi(:,:)*zdt_split
    526          zqfi(:,:,:)=zqfi(:,:,:)+zdqfi(:,:,:)*zdt_split
    527 
    528          zdufic(:,:)=zdufic(:,:)+zdufi(:,:)
    529          zdvfic(:,:)=zdvfic(:,:)+zdvfi(:,:)
    530          zdtfic(:,:)=zdtfic(:,:)+zdtfi(:,:)
    531          zdqfic(:,:,:)=zdqfic(:,:,:)+zdqfi(:,:,:)
    532 
    533        enddo ! of do isplit=1,nsplit_phys
    534 
    535        END IF
    536 
    537       zdufi(:,:)=zdufic(:,:)/nsplit_phys
    538       zdvfi(:,:)=zdvfic(:,:)/nsplit_phys
    539       zdtfi(:,:)=zdtfic(:,:)/nsplit_phys
    540       zdqfi(:,:,:)=zdqfic(:,:,:)/nsplit_phys
    541 
    542 c-----------------------------------------------------------------------
    543 c   transformation des tendances physiques en tendances dynamiques:
    544 c   ---------------------------------------------------------------
    545 
    546 c  tendance sur la pression :
    547 c  -----------------------------------
    548 
    549       CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,zdpsrf,pdpsfi)
    550 c
    551 c   62. enthalpie potentielle
    552 c   ---------------------
    553 
    554       DO l=1,llm
    555 
    556          DO i=1,iip1
    557           pdhfi(i,1,l)    = cpp *  zdtfi(1,l)      / ppk(i, 1  ,l)
    558           pdhfi(i,jjp1,l) = cpp *  zdtfi(ngridmx,l)/ ppk(i,jjp1,l)
    559          ENDDO
    560 
    561          DO j=2,jjm
    562             ig0=1+(j-2)*iim
    563             DO i=1,iim
    564                pdhfi(i,j,l) = cpp * zdtfi(ig0+i,l) / ppk(i,j,l)
    565             ENDDO
    566                pdhfi(iip1,j,l) =  pdhfi(1,j,l)
    567          ENDDO
    568 
    569       ENDDO
    570 
    571 
    572 c   62. humidite specifique
    573 c   ---------------------
    574 ! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
    575 !      DO iq=1,nqtot
    576 !         DO l=1,llm
    577 !            DO i=1,iip1
    578 !               pdqfi(i,1,l,iq)    = zdqfi(1,l,iq)
    579 !               pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,iq)
    580 !            ENDDO
    581 !            DO j=2,jjm
    582 !               ig0=1+(j-2)*iim
    583 !               DO i=1,iim
    584 !                  pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,iq)
    585 !               ENDDO
    586 !               pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,iq)
    587 !            ENDDO
    588 !         ENDDO
    589 !      ENDDO
    590 
    591 c   63. traceurs
    592 c   ------------
    593 C     initialisation des tendances
    594       pdqfi(:,:,:,:)=0.
    595 C
    596       itr = 0
    597       DO iq=1,nqtot
    598          IF(.NOT.tracers(iq)%isAdvected) CYCLE
    599          itr = itr + 1
    600          DO l=1,llm
    601             DO i=1,iip1
    602                pdqfi(i,1,l,iq)    = zdqfi(1,l,itr)
    603                pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,itr)
    604             ENDDO
    605             DO j=2,jjm
    606                ig0=1+(j-2)*iim
    607                DO i=1,iim
    608                   pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,itr)
    609                ENDDO
    610                pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,itr)
    611             ENDDO
    612          ENDDO
    613       ENDDO
    614 
    615 c   65. champ u:
    616 c   ------------
    617 
    618       DO l=1,llm
    619 
    620          DO i=1,iip1
    621             pdufi(i,1,l)    = 0.
    622             pdufi(i,jjp1,l) = 0.
    623          ENDDO
    624 
    625          DO j=2,jjm
    626             ig0=1+(j-2)*iim
    627             DO i=1,iim-1
    628                pdufi(i,j,l)=
    629      $         0.5*(zdufi(ig0+i,l)+zdufi(ig0+i+1,l))*cu(i,j)
    630             ENDDO
    631             pdufi(iim,j,l)=
    632      $      0.5*(zdufi(ig0+1,l)+zdufi(ig0+iim,l))*cu(iim,j)
    633             pdufi(iip1,j,l)=pdufi(1,j,l)
    634          ENDDO
    635 
    636       ENDDO
    637 
    638 
    639 c   67. champ v:
    640 c   ------------
    641 
    642       DO l=1,llm
    643 
    644          DO j=2,jjm-1
    645             ig0=1+(j-2)*iim
    646             DO i=1,iim
    647                pdvfi(i,j,l)=
    648      $         0.5*(zdvfi(ig0+i,l)+zdvfi(ig0+i+iim,l))*cv(i,j)
    649             ENDDO
    650             pdvfi(iip1,j,l) = pdvfi(1,j,l)
    651          ENDDO
    652       ENDDO
    653 
    654 
    655 c   68. champ v pres des poles:
    656 c   ---------------------------
    657 c      v = U * cos(long) + V * SIN(long)
    658 
    659       DO l=1,llm
    660 
    661          DO i=1,iim
    662             pdvfi(i,1,l)=
    663      $      zdufi(1,l)*COS(rlonv(i))+zdvfi(1,l)*SIN(rlonv(i))
    664             pdvfi(i,jjm,l)=zdufi(ngridmx,l)*COS(rlonv(i))
    665      $      +zdvfi(ngridmx,l)*SIN(rlonv(i))
    666             pdvfi(i,1,l)=
    667      $      0.5*(pdvfi(i,1,l)+zdvfi(i+1,l))*cv(i,1)
    668             pdvfi(i,jjm,l)=
    669      $      0.5*(pdvfi(i,jjm,l)+zdvfi(ngridmx-iip1+i,l))*cv(i,jjm)
    670           ENDDO
    671 
    672          pdvfi(iip1,1,l)  = pdvfi(1,1,l)
    673          pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
    674 
    675       ENDDO
    676 
    677 c-----------------------------------------------------------------------
    678       firstcal = .FALSE.
    679 
    680       RETURN
    681       END
     672     pdvfi(iip1,1,l)  = pdvfi(1,1,l)
     673     pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
     674
     675  ENDDO
     676
     677  !-----------------------------------------------------------------------
     678  firstcal = .FALSE.
     679
     680  RETURN
     681END SUBROUTINE calfis
Note: See TracChangeset for help on using the changeset viewer.