Head

GCC Code Coverage Report

Directory:	./		Exec	Total	Coverage
File:	dynphy_lonlat/calfis.F	Lines:	198	198	100.0 %
Date:	2023-06-30 12:56:34	Branches:	190	196	96.9 %


!
! $Id: calfis.F 4464 2023-03-09 17:03:51Z lguez $
!
C
C
      SUBROUTINE calfis(lafin,
     $                  jD_cur, jH_cur,
     $                  pucov,
     $                  pvcov,
     $                  pteta,
     $                  pq,
     $                  pmasse,
     $                  pps,
     $                  pp,
     $                  ppk,
     $                  pphis,
     $                  pphi,
     $                  pducov,
     $                  pdvcov,
     $                  pdteta,
     $                  pdq,
     $                  flxw,
     $                  pdufi,
     $                  pdvfi,
     $                  pdhfi,
     $                  pdqfi,
     $                  pdpsfi)
c
c    Auteur :  P. Le Van, F. Hourdin
c   .........
      USE infotrac, ONLY: nqtot, tracers
      USE control_mod, ONLY: planet_type, nsplit_phys
#ifdef CPP_PHYS
      USE callphysiq_mod, ONLY: call_physiq
#endif
      USE comconst_mod, ONLY: cpp, daysec, dtphys, dtvr, kappa, pi
      USE comvert_mod, ONLY: preff, presnivs

      IMPLICIT NONE
c=======================================================================
c
c   1. rearrangement des tableaux et transformation
c      variables dynamiques  >  variables physiques
c   2. calcul des termes physiques
c   3. retransformation des tendances physiques en tendances dynamiques
c
c   remarques:
c   ----------
c
c    - les vents sont donnes dans la physique par leurs composantes
c      naturelles.
c    - la variable thermodynamique de la physique est une variable
c      intensive :   T
c      pour la dynamique on prend    T * ( preff / p(l) ) **kappa
c    - les deux seules variables dependant de la geometrie necessaires
c      pour la physique sont la latitude pour le rayonnement et
c      l'aire de la maille quand on veut integrer une grandeur
c      horizontalement.
c    - les points de la physique sont les points scalaires de la
c      la dynamique; numerotation:
c          1 pour le pole nord
c          (jjm-1)*iim pour l'interieur du domaine
c          ngridmx pour le pole sud
c      ---> ngridmx=2+(jjm-1)*iim
c
c     Input :
c     -------
c       pucov           covariant zonal velocity
c       pvcov           covariant meridional velocity
c       pteta           potential temperature
c       pps             surface pressure
c       pmasse          masse d'air dans chaque maille
c       pts             surface temperature  (K)
c       callrad         clef d'appel au rayonnement
c
c    Output :
c    --------
c        pdufi          tendency for the natural zonal velocity (ms-1)
c        pdvfi          tendency for the natural meridional velocity
c        pdhfi          tendency for the potential temperature
c        pdtsfi         tendency for the surface temperature
c
c        pdtrad         radiative tendencies  \  both input
c        pfluxrad       radiative fluxes      /  and output
c
c=======================================================================
c
c-----------------------------------------------------------------------
c
c    0.  Declarations :
c    ------------------

      include "dimensions.h"
      include "paramet.h"

      INTEGER ngridmx
      PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )

      include "comgeom2.h"
      include "iniprint.h"

c    Arguments :
c    -----------
      LOGICAL,INTENT(IN) ::  lafin ! .true. for the very last call to physics
      REAL,INTENT(IN):: jD_cur, jH_cur
      REAL,INTENT(IN) :: pvcov(iip1,jjm,llm) ! covariant meridional velocity
      REAL,INTENT(IN) :: pucov(iip1,jjp1,llm) ! covariant zonal velocity
      REAL,INTENT(IN) :: pteta(iip1,jjp1,llm) ! potential temperature
      REAL,INTENT(IN) :: pmasse(iip1,jjp1,llm) ! mass in each cell ! not used
      REAL,INTENT(IN) :: pq(iip1,jjp1,llm,nqtot) ! tracers
      REAL,INTENT(IN) :: pphis(iip1,jjp1) ! surface geopotential
      REAL,INTENT(IN) :: pphi(iip1,jjp1,llm) ! geopotential

      REAL,INTENT(IN) :: pdvcov(iip1,jjm,llm) ! dynamical tendency on vcov
      REAL,INTENT(IN) :: pducov(iip1,jjp1,llm) ! dynamical tendency on ucov
      REAL,INTENT(IN) :: pdteta(iip1,jjp1,llm) ! dynamical tendency on teta
      ! NB: pdteta is used only to compute pcvgt which is in fact not used...
      REAL,INTENT(IN) :: pdq(iip1,jjp1,llm,nqtot) ! dynamical tendency on tracers
      ! NB: pdq is only used to compute pcvgq which is in fact not used...

      REAL,INTENT(IN) :: pps(iip1,jjp1) ! surface pressure (Pa)
      REAL,INTENT(IN) :: pp(iip1,jjp1,llmp1) ! pressure at mesh interfaces (Pa)
      REAL,INTENT(IN) :: ppk(iip1,jjp1,llm) ! Exner at mid-layer
      REAL,INTENT(IN) :: flxw(iip1,jjp1,llm) ! Vertical mass flux on lower mesh interfaces (kg/s) (on llm because flxw(:,:,llm+1)=0)

      ! tendencies (in */s) from the physics
      REAL,INTENT(OUT) :: pdvfi(iip1,jjm,llm) ! tendency on covariant meridional wind
      REAL,INTENT(OUT) :: pdufi(iip1,jjp1,llm) ! tendency on covariant zonal wind
      REAL,INTENT(OUT) :: pdhfi(iip1,jjp1,llm) ! tendency on potential temperature (K/s)
      REAL,INTENT(OUT) :: pdqfi(iip1,jjp1,llm,nqtot) ! tendency on tracers
      REAL,INTENT(OUT) :: pdpsfi(iip1,jjp1) ! tendency on surface pressure (Pa/s)


c    Local variables :
c    -----------------

      INTEGER i,j,l,ig0,ig,iq,itr
      REAL zpsrf(ngridmx)
      REAL zplev(ngridmx,llm+1),zplay(ngridmx,llm)
      REAL zphi(ngridmx,llm),zphis(ngridmx)
c
      REAL zrot(iip1,jjm,llm) ! AdlC May 2014
      REAL zufi(ngridmx,llm), zvfi(ngridmx,llm)
      REAL zrfi(ngridmx,llm) ! relative wind vorticity
      REAL ztfi(ngridmx,llm),zqfi(ngridmx,llm,nqtot)
      REAL zpk(ngridmx,llm)
c
      REAL pcvgu(ngridmx,llm), pcvgv(ngridmx,llm)
      REAL pcvgt(ngridmx,llm), pcvgq(ngridmx,llm,2)
c
      REAL zdufi(ngridmx,llm),zdvfi(ngridmx,llm)
      REAL zdtfi(ngridmx,llm),zdqfi(ngridmx,llm,nqtot)
      REAL zdpsrf(ngridmx)
c
      REAL zdufic(ngridmx,llm),zdvfic(ngridmx,llm)
      REAL zdtfic(ngridmx,llm),zdqfic(ngridmx,llm,nqtot)
      REAL jH_cur_split,zdt_split
      LOGICAL debut_split,lafin_split
      INTEGER isplit

      REAL zsin(iim),zcos(iim),z1(iim)
      REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
      REAL unskap, pksurcp
c
      REAL flxwfi(ngridmx,llm)  ! Flux de masse verticale sur la grille physiq
c

      REAL SSUM

      LOGICAL,SAVE :: firstcal=.true., debut=.true.
!      REAL rdayvrai

c
c-----------------------------------------------------------------------
c
c    1. Initialisations :
c    --------------------
c
c
      IF ( firstcal )  THEN
        debut = .TRUE.
        IF (ngridmx.NE.2+(jjm-1)*iim) THEN
         write(lunout,*) 'STOP dans calfis'
         write(lunout,*)
     &   'La dimension ngridmx doit etre egale a 2 + (jjm-1)*iim'
         write(lunout,*) '  ngridmx  jjm   iim   '
         write(lunout,*) ngridmx,jjm,iim
         call abort_gcm("calfis", "", 1)
        ENDIF
      ELSE
        debut = .FALSE.
      ENDIF ! of IF (firstcal)

c
c
c-----------------------------------------------------------------------
c   40. transformation des variables dynamiques en variables physiques:
c   ---------------------------------------------------------------

c   41. pressions au sol (en Pascals)
c   ----------------------------------


      zpsrf(1) = pps(1,1)

      ig0  = 2
      DO j = 2,jjm
         CALL SCOPY( iim,pps(1,j),1,zpsrf(ig0), 1 )
         ig0 = ig0+iim
      ENDDO

      zpsrf(ngridmx) = pps(1,jjp1)


c   42. pression intercouches et fonction d'Exner:
c
c   -----------------------------------------------------------------
c     .... zplev  definis aux (llm +1) interfaces des couches  ....
c     .... zplay  definis aux (  llm )    milieux des couches  ....
c   -----------------------------------------------------------------

c    ...    Exner = cp * ( p(l) / preff ) ** kappa     ....
c
       unskap   = 1./ kappa
c
      DO l = 1, llm
        zpk(   1,l ) = ppk(1,1,l)
        zplev( 1,l ) = pp(1,1,l)
        ig0 = 2
          DO j = 2, jjm
             DO i =1, iim
              zpk(   ig0,l ) = ppk(i,j,l)
              zplev( ig0,l ) = pp(i,j,l)
              ig0 = ig0 +1
             ENDDO
          ENDDO
        zpk(   ngridmx,l ) = ppk(1,jjp1,l)
        zplev( ngridmx,l ) = pp(1,jjp1,l)
      ENDDO
        zplev( 1,llmp1 ) = pp(1,1,llmp1)
        ig0 = 2
          DO j = 2, jjm
             DO i =1, iim
              zplev( ig0,llmp1 ) = pp(i,j,llmp1)
              ig0 = ig0 +1
             ENDDO
          ENDDO
        zplev( ngridmx,llmp1 ) = pp(1,jjp1,llmp1)
c
c

c   43. temperature naturelle (en K) et pressions milieux couches .
c   ---------------------------------------------------------------

      DO l=1,llm

         pksurcp     =  ppk(1,1,l) / cpp
         zplay(1,l)  =  preff * pksurcp ** unskap
         ztfi(1,l)   =  pteta(1,1,l) *  pksurcp
         pcvgt(1,l)  =  pdteta(1,1,l) * pksurcp / pmasse(1,1,l)
         ig0         = 2

         DO j = 2, jjm
            DO i = 1, iim
              pksurcp        = ppk(i,j,l) / cpp
              zplay(ig0,l)   = preff * pksurcp ** unskap
              ztfi(ig0,l)    = pteta(i,j,l)  * pksurcp
              pcvgt(ig0,l)   = pdteta(i,j,l) * pksurcp / pmasse(i,j,l)
              ig0            = ig0 + 1
            ENDDO
         ENDDO

         pksurcp       = ppk(1,jjp1,l) / cpp
         zplay(ig0,l)  = preff * pksurcp ** unskap
         ztfi (ig0,l)  = pteta(1,jjp1,l)  * pksurcp
         pcvgt(ig0,l)  = pdteta(1,jjp1,l) * pksurcp/ pmasse(1,jjp1,l)

      ENDDO

c   43.bis traceurs
c   ---------------
c
      itr=0
      DO iq=1,nqtot
         IF(.NOT.tracers(iq)%isAdvected) CYCLE
         itr = itr + 1
         DO l=1,llm
            zqfi(1,l,itr) = pq(1,1,l,iq)
            ig0           = 2
            DO j=2,jjm
               DO i = 1, iim
                  zqfi(ig0,l,itr) = pq(i,j,l,iq)
                  ig0             = ig0 + 1
               ENDDO
            ENDDO
            zqfi(ig0,l,itr) = pq(1,jjp1,l,iq)
         ENDDO
      ENDDO

c   convergence dynamique pour les traceurs "EAU"
! Earth-specific treatment of first 2 tracers (water)
       if (planet_type=="earth") then
        DO iq=1,2
         DO l=1,llm
            pcvgq(1,l,iq)= pdq(1,1,l,iq) / pmasse(1,1,l)
            ig0          = 2
            DO j=2,jjm
               DO i = 1, iim
                  pcvgq(ig0,l,iq) = pdq(i,j,l,iq) / pmasse(i,j,l)
                  ig0             = ig0 + 1
               ENDDO
            ENDDO
            pcvgq(ig0,l,iq)= pdq(1,jjp1,l,iq) / pmasse(1,jjp1,l)
         ENDDO
        ENDDO
       endif ! of if (planet_type=="earth")


c   Geopotentiel calcule par rapport a la surface locale:
c   -----------------------------------------------------

      CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,pphi,zphi)
      CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,pphis,zphis)
      DO l=1,llm
         DO ig=1,ngridmx
           zphi(ig,l)=zphi(ig,l)-zphis(ig)
         ENDDO
      ENDDO

c   ....  Calcul de la vitesse  verticale  ( en Pa*m*s  ou Kg/s )  ....
c JG : ancien calcule de omega utilise dans physiq.F. Maintenant le flux
c    de masse est calclue dans advtrac.F
c      DO l=1,llm
c        pvervel(1,l)=pw(1,1,l) * g /apoln
c        ig0=2
c       DO j=2,jjm
c           DO i = 1, iim
c              pvervel(ig0,l) = pw(i,j,l) * g * unsaire(i,j)
c              ig0 = ig0 + 1
c           ENDDO
c       ENDDO
c        pvervel(ig0,l)=pw(1,jjp1,l) * g /apols
c      ENDDO

c
c   45. champ u:
c   ------------

      DO 50 l=1,llm

         DO 25 j=2,jjm
            ig0 = 1+(j-2)*iim
            zufi(ig0+1,l)= 0.5 *
     $      ( pucov(iim,j,l)/cu(iim,j) + pucov(1,j,l)/cu(1,j) )
            pcvgu(ig0+1,l)= 0.5 *
     $      ( pducov(iim,j,l)/cu(iim,j) + pducov(1,j,l)/cu(1,j) )
            DO 10 i=2,iim
               zufi(ig0+i,l)= 0.5 *
     $         ( pucov(i-1,j,l)/cu(i-1,j) + pucov(i,j,l)/cu(i,j) )
               pcvgu(ig0+i,l)= 0.5 *
     $         ( pducov(i-1,j,l)/cu(i-1,j) + pducov(i,j,l)/cu(i,j) )
10         CONTINUE
25      CONTINUE

50    CONTINUE


C  Alvaro de la Camara (May 2014)
C  46.1 Calcul de la vorticite et passage sur la grille physique
C  --------------------------------------------------------------
      DO l=1,llm
        do i=1,iim
          do j=1,jjm
            zrot(i,j,l) = (pvcov(i+1,j,l) - pvcov(i,j,l)
     $                   + pucov(i,j+1,l) - pucov(i,j,l))
     $                   / (cu(i,j)+cu(i,j+1))
     $                   / (cv(i+1,j)+cv(i,j)) *4
          enddo
        enddo
      ENDDO

c   46.champ v:
c   -----------

      DO l=1,llm
         DO j=2,jjm
            ig0=1+(j-2)*iim
            DO i=1,iim
               zvfi(ig0+i,l)= 0.5 *
     $         ( pvcov(i,j-1,l)/cv(i,j-1) + pvcov(i,j,l)/cv(i,j) )
               pcvgv(ig0+i,l)= 0.5 *
     $         ( pdvcov(i,j-1,l)/cv(i,j-1) + pdvcov(i,j,l)/cv(i,j) )
            ENDDO
               zrfi(ig0 + 1,l)= 0.25 *(zrot(iim,j-1,l)+zrot(iim,j,l)
     &                                +zrot(1,j-1,l)+zrot(1,j,l))
            DO i=2,iim
               zrfi(ig0 + i,l)= 0.25 *(zrot(i-1,j-1,l)+zrot(i-1,j,l)
     $                   +zrot(i,j-1,l)+zrot(i,j,l))   !  AdlC MAY 2014
            ENDDO
         ENDDO
      ENDDO


c   47. champs de vents aux pole nord
c   ------------------------------
c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]

      DO l=1,llm

         z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,1,l)/cv(1,1)
         z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,1,l)/cv(1,1)
         DO i=2,iim
            z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
            z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,1,l)/cv(i,1)
         ENDDO

         DO i=1,iim
            zcos(i)   = COS(rlonv(i))*z1(i)
            zcosbis(i)= COS(rlonv(i))*z1bis(i)
            zsin(i)   = SIN(rlonv(i))*z1(i)
            zsinbis(i)= SIN(rlonv(i))*z1bis(i)
         ENDDO

         zufi(1,l)  = SSUM(iim,zcos,1)/pi
         pcvgu(1,l) = SSUM(iim,zcosbis,1)/pi
         zvfi(1,l)  = SSUM(iim,zsin,1)/pi
         pcvgv(1,l) = SSUM(iim,zsinbis,1)/pi
         zrfi(1, l) = 0.
      ENDDO


c   48. champs de vents aux pole sud:
c   ---------------------------------
c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]

      DO l=1,llm

         z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,jjm,l)/cv(1,jjm)
         z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,jjm,l)/cv(1,jjm)
         DO i=2,iim
            z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
            z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,jjm,l)/cv(i,jjm)
         ENDDO

         DO i=1,iim
            zcos(i)    = COS(rlonv(i))*z1(i)
            zcosbis(i) = COS(rlonv(i))*z1bis(i)
            zsin(i)    = SIN(rlonv(i))*z1(i)
            zsinbis(i) = SIN(rlonv(i))*z1bis(i)
         ENDDO

         zufi(ngridmx,l)  = SSUM(iim,zcos,1)/pi
         pcvgu(ngridmx,l) = SSUM(iim,zcosbis,1)/pi
         zvfi(ngridmx,l)  = SSUM(iim,zsin,1)/pi
         pcvgv(ngridmx,l) = SSUM(iim,zsinbis,1)/pi
         zrfi(ngridmx, l) = 0.
      ENDDO
c
c On change de grille, dynamique vers physiq, pour le flux de masse verticale
      CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,flxw,flxwfi)

c-----------------------------------------------------------------------
c   Appel de la physique:
c   ---------------------



!      write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
      zdt_split=dtphys/nsplit_phys
      zdufic(:,:)=0.
      zdvfic(:,:)=0.
      zdtfic(:,:)=0.
      zdqfic(:,:,:)=0.

#ifdef CPP_PHYS

       do isplit=1,nsplit_phys

         jH_cur_split=jH_cur+(isplit-1) * dtvr / (daysec *nsplit_phys)
         debut_split=debut.and.isplit==1
         lafin_split=lafin.and.isplit==nsplit_phys

!      if (planet_type=="earth") then
        CALL call_physiq(ngridmx,llm,nqtot,tracers(:)%name,
     &                   debut_split,lafin_split,
     &                   jD_cur,jH_cur_split,zdt_split,
     &                   zplev,zplay,
     &                   zpk,zphi,zphis,
     &                   presnivs,
     &                   zufi,zvfi,zrfi,ztfi,zqfi,
     &                   flxwfi,pducov,
     &                   zdufi,zdvfi,zdtfi,zdqfi,zdpsrf)
!
!      else if ( planet_type=="generic" ) then
!
!         CALL physiq (ngridmx,     !! ngrid
!     .             llm,            !! nlayer
!     .             nqtot,          !! nq
!     .             tracers(:)%name,!! tracer names from dynamical core (given in infotrac)
!     .             debut_split,    !! firstcall
!     .             lafin_split,    !! lastcall
!     .             jD_cur,         !! pday. see leapfrog
!     .             jH_cur_split,   !! ptime "fraction of day"
!     .             zdt_split,      !! ptimestep
!     .             zplev,          !! pplev
!     .             zplay,          !! pplay
!     .             zphi,           !! pphi
!     .             zufi,           !! pu
!     .             zvfi,           !! pv
!     .             ztfi,           !! pt
!     .             zqfi,           !! pq
!     .             flxwfi,         !! pw !! or 0. anyway this is for diagnostic. not used in physiq.
!     .             zdufi,          !! pdu
!     .             zdvfi,          !! pdv
!     .             zdtfi,          !! pdt
!     .             zdqfi,          !! pdq
!     .             zdpsrf,         !! pdpsrf
!     .             tracerdyn)      !! tracerdyn <-- utilite ???
!
!      endif ! of if (planet_type=="earth")

         zufi(:,:)=zufi(:,:)+zdufi(:,:)*zdt_split
         zvfi(:,:)=zvfi(:,:)+zdvfi(:,:)*zdt_split
         ztfi(:,:)=ztfi(:,:)+zdtfi(:,:)*zdt_split
         zqfi(:,:,:)=zqfi(:,:,:)+zdqfi(:,:,:)*zdt_split

         zdufic(:,:)=zdufic(:,:)+zdufi(:,:)
         zdvfic(:,:)=zdvfic(:,:)+zdvfi(:,:)
         zdtfic(:,:)=zdtfic(:,:)+zdtfi(:,:)
         zdqfic(:,:,:)=zdqfic(:,:,:)+zdqfi(:,:,:)

       enddo ! of do isplit=1,nsplit_phys

#endif
! of #ifdef CPP_PHYS

      zdufi(:,:)=zdufic(:,:)/nsplit_phys
      zdvfi(:,:)=zdvfic(:,:)/nsplit_phys
      zdtfi(:,:)=zdtfic(:,:)/nsplit_phys
      zdqfi(:,:,:)=zdqfic(:,:,:)/nsplit_phys


500   CONTINUE

c-----------------------------------------------------------------------
c   transformation des tendances physiques en tendances dynamiques:
c   ---------------------------------------------------------------

c  tendance sur la pression :
c  -----------------------------------

      CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,zdpsrf,pdpsfi)
c
c   62. enthalpie potentielle
c   ---------------------

      DO l=1,llm

         DO i=1,iip1
          pdhfi(i,1,l)    = cpp *  zdtfi(1,l)      / ppk(i, 1  ,l)
          pdhfi(i,jjp1,l) = cpp *  zdtfi(ngridmx,l)/ ppk(i,jjp1,l)
         ENDDO

         DO j=2,jjm
            ig0=1+(j-2)*iim
            DO i=1,iim
               pdhfi(i,j,l) = cpp * zdtfi(ig0+i,l) / ppk(i,j,l)
            ENDDO
               pdhfi(iip1,j,l) =  pdhfi(1,j,l)
         ENDDO

      ENDDO


c   62. humidite specifique
c   ---------------------
! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
!      DO iq=1,nqtot
!         DO l=1,llm
!            DO i=1,iip1
!               pdqfi(i,1,l,iq)    = zdqfi(1,l,iq)
!               pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,iq)
!            ENDDO
!            DO j=2,jjm
!               ig0=1+(j-2)*iim
!               DO i=1,iim
!                  pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,iq)
!               ENDDO
!               pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,iq)
!            ENDDO
!         ENDDO
!      ENDDO

c   63. traceurs
c   ------------
C     initialisation des tendances
      pdqfi(:,:,:,:)=0.
C
      itr = 0
      DO iq=1,nqtot
         IF(.NOT.tracers(iq)%isAdvected) CYCLE
         itr = itr + 1
         DO l=1,llm
            DO i=1,iip1
               pdqfi(i,1,l,iq)    = zdqfi(1,l,itr)
               pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,itr)
            ENDDO
            DO j=2,jjm
               ig0=1+(j-2)*iim
               DO i=1,iim
                  pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,itr)
               ENDDO
               pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,itr)
            ENDDO
         ENDDO
      ENDDO

c   65. champ u:
c   ------------

      DO l=1,llm

         DO i=1,iip1
            pdufi(i,1,l)    = 0.
            pdufi(i,jjp1,l) = 0.
         ENDDO

         DO j=2,jjm
            ig0=1+(j-2)*iim
            DO i=1,iim-1
               pdufi(i,j,l)=
     $         0.5*(zdufi(ig0+i,l)+zdufi(ig0+i+1,l))*cu(i,j)
            ENDDO
            pdufi(iim,j,l)=
     $      0.5*(zdufi(ig0+1,l)+zdufi(ig0+iim,l))*cu(iim,j)
            pdufi(iip1,j,l)=pdufi(1,j,l)
         ENDDO

      ENDDO


c   67. champ v:
c   ------------

      DO l=1,llm

         DO j=2,jjm-1
            ig0=1+(j-2)*iim
            DO i=1,iim
               pdvfi(i,j,l)=
     $         0.5*(zdvfi(ig0+i,l)+zdvfi(ig0+i+iim,l))*cv(i,j)
            ENDDO
            pdvfi(iip1,j,l) = pdvfi(1,j,l)
         ENDDO
      ENDDO


c   68. champ v pres des poles:
c   ---------------------------
c      v = U * cos(long) + V * SIN(long)

      DO l=1,llm

         DO i=1,iim
            pdvfi(i,1,l)=
     $      zdufi(1,l)*COS(rlonv(i))+zdvfi(1,l)*SIN(rlonv(i))
            pdvfi(i,jjm,l)=zdufi(ngridmx,l)*COS(rlonv(i))
     $      +zdvfi(ngridmx,l)*SIN(rlonv(i))
            pdvfi(i,1,l)=
     $      0.5*(pdvfi(i,1,l)+zdvfi(i+1,l))*cv(i,1)
            pdvfi(i,jjm,l)=
     $      0.5*(pdvfi(i,jjm,l)+zdvfi(ngridmx-iip1+i,l))*cv(i,jjm)
          ENDDO

         pdvfi(iip1,1,l)  = pdvfi(1,1,l)
         pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)

      ENDDO

c-----------------------------------------------------------------------

700   CONTINUE

      firstcal = .FALSE.

      RETURN
      END


Generated by: GCOVR (Version 4.2)

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