trunk/dyn3d/calfis.f

module calfis_m

  ! Clean: no C preprocessor directive, no include line

  IMPLICIT NONE

contains

  SUBROUTINE calfis(nq, lafin, rdayvrai, heure, pucov, pvcov, pteta, pq, &
       pmasse, pps, ppk, pphis, pphi, pducov, pdvcov, pdteta, pdq, pw, &
       clesphy0, pdufi, pdvfi, pdhfi, pdqfi, pdpsfi)

    ! From dyn3d/calfis.F,v 1.3 2005/05/25 13:10:09

    ! Auteurs : P. Le Van, F. Hourdin

    !   1. rearrangement des tableaux et transformation
    !      variables dynamiques  >  variables physiques
    !   2. calcul des termes physiques
    !   3. retransformation des tendances physiques en tendances dynamiques

    !   remarques:
    !   ----------

    !    - les vents sont donnes dans la physique par leurs composantes 
    !      naturelles.
    !    - la variable thermodynamique de la physique est une variable
    !      intensive :   T 
    !      pour la dynamique on prend    T * (preff / p(l)) **kappa
    !    - les deux seules variables dependant de la geometrie necessaires
    !      pour la physique sont la latitude pour le rayonnement et 
    !      l'aire de la maille quand on veut integrer une grandeur 
    !      horizontalement.

    !     Input :
    !     -------
    !       pucov           covariant zonal velocity
    !       pvcov           covariant meridional velocity 
    !       pteta           potential temperature
    !       pps             surface pressure
    !       pmasse          masse d'air dans chaque maille
    !       pts             surface temperature  (K)
    !       callrad         clef d'appel au rayonnement

    !    Output :
    !    --------
    !        pdufi          tendency for the natural zonal velocity (ms-1)
    !        pdvfi          tendency for the natural meridional velocity 
    !        pdhfi          tendency for the potential temperature
    !        pdtsfi         tendency for the surface temperature

    !        pdtrad         radiative tendencies  \  both input
    !        pfluxrad       radiative fluxes      /  and output

    use dimens_m, only: iim, jjm, llm, nqmx
    use dimphy, only: klon
    use comconst, only: kappa, cpp, dtphys, g, pi
    use comvert, only: preff, presnivs
    use comgeom, only: apoln, cu_2d, cv_2d, unsaire_2d, apols, rlonu, rlonv
    use advtrac_m, only: niadv
    use grid_change, only: dyn_phy, gr_fi_dyn
    use physiq_m, only: physiq
    use pressure_var, only: p3d, pls

    !    0.  Declarations :

    INTEGER nq

    !    Arguments :

    LOGICAL, intent(in):: lafin
    REAL, intent(in):: heure ! heure de la journée en fraction de jour

    REAL pvcov(iim + 1,jjm,llm)
    REAL pucov(iim + 1,jjm + 1,llm)
    REAL pteta(iim + 1,jjm + 1,llm)
    REAL pmasse(iim + 1,jjm + 1,llm)

    REAL, intent(in):: pq(iim + 1,jjm + 1,llm,nqmx)
    ! (mass fractions of advected fields)

    REAL pphis(iim + 1,jjm + 1)
    REAL pphi(iim + 1,jjm + 1,llm)

    REAL pdvcov(iim + 1,jjm,llm)
    REAL pducov(iim + 1,jjm + 1,llm)
    REAL pdteta(iim + 1,jjm + 1,llm)
    REAL pdq(iim + 1,jjm + 1,llm,nqmx)

    REAL pw(iim + 1,jjm + 1,llm)

    REAL pps(iim + 1,jjm + 1)
    REAL, intent(in):: ppk(iim + 1,jjm + 1,llm)

    REAL pdvfi(iim + 1,jjm,llm)
    REAL pdufi(iim + 1,jjm + 1,llm)
    REAL pdhfi(iim + 1,jjm + 1,llm)
    REAL pdqfi(iim + 1,jjm + 1,llm,nqmx)
    REAL pdpsfi(iim + 1,jjm + 1)

    INTEGER, PARAMETER:: longcles = 20
    REAL, intent(in):: clesphy0(longcles)

    !    Local variables :

    INTEGER i,j,l,ig0,ig,iq,iiq
    REAL zpsrf(klon)
    REAL zplev(klon,llm+1),zplay(klon,llm)
    REAL zphi(klon,llm),zphis(klon)

    REAL zufi(klon,llm), zvfi(klon,llm)
    REAL ztfi(klon,llm) ! temperature
    real zqfi(klon,llm,nqmx) ! mass fractions of advected fields

    REAL pcvgu(klon,llm), pcvgv(klon,llm)
    REAL pcvgt(klon,llm), pcvgq(klon,llm,2)

    REAL pvervel(klon,llm)

    REAL zdufi(klon,llm),zdvfi(klon,llm)
    REAL zdtfi(klon,llm),zdqfi(klon,llm,nqmx)
    REAL zdpsrf(klon)

    REAL zsin(iim),zcos(iim),z1(iim)
    REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
    REAL pksurcp(iim + 1,jjm + 1)

    ! I. Musat: diagnostic PVteta, Amip2
    INTEGER, PARAMETER:: ntetaSTD=3
    REAL:: rtetaSTD(ntetaSTD) = (/350., 380., 405./)
    REAL PVteta(klon,ntetaSTD)

    REAL SSUM

    LOGICAL:: firstcal = .true.
    REAL, intent(in):: rdayvrai

    !-----------------------------------------------------------------------

    !!print *, "Call sequence information: calfis"

    !    1. Initialisations :
    !   latitude, longitude et aires des mailles pour la physique:

    !   40. transformation des variables dynamiques en variables physiques:
    !   41. pressions au sol (en Pascals)

    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(klon) = pps(1,jjm + 1)

    !   42. pression intercouches :

    !     .... zplev  definis aux (llm +1) interfaces des couches  ....
    !     .... zplay  definis aux (llm)    milieux des couches  .... 

    !    ...    Exner = cp * (p(l) / preff) ** kappa     ....

    forall (l = 1: llm+1) zplev(:, l) = pack(p3d(:, :, l), dyn_phy)

    !   43. temperature naturelle (en K) et pressions milieux couches .
    DO l=1,llm
       pksurcp     =  ppk(:, :, l) / cpp
       pls(:, :, l) = preff * pksurcp**(1./ kappa)
       zplay(:, l) = pack(pls(:, :, l), dyn_phy)
       ztfi(:, l) = pack(pteta(:, :, l) * pksurcp, dyn_phy)
       pcvgt(:, l) = pack(pdteta(:, :, l) * pksurcp / pmasse(:, :, l), dyn_phy)
    ENDDO

    !   43.bis traceurs

    DO iq=1,nq
       iiq=niadv(iq) 
       DO l=1,llm
          zqfi(1,l,iq) = pq(1,1,l,iiq)
          ig0          = 2
          DO j=2,jjm
             DO i = 1, iim
                zqfi(ig0,l,iq)  = pq(i,j,l,iiq)
                ig0             = ig0 + 1
             ENDDO
          ENDDO
          zqfi(ig0,l,iq) = pq(1,jjm + 1,l,iiq)
       ENDDO
    ENDDO

    !   convergence dynamique pour les traceurs "EAU"

    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,jjm + 1,l,iq) / pmasse(1,jjm + 1,l)
       ENDDO
    ENDDO

    !   Geopotentiel calcule par rapport a la surface locale:

    forall (l = 1:llm) zphi(:, l) = pack(pphi(:, :, l), dyn_phy)
    zphis = pack(pphis, dyn_phy)
    DO l=1,llm
       DO ig=1,klon
          zphi(ig,l)=zphi(ig,l)-zphis(ig)
       ENDDO
    ENDDO

    !   ....  Calcul de la vitesse  verticale  (en Pa*m*s  ou Kg/s)  ....

    DO l=1,llm
       pvervel(1,l)=pw(1,1,l) * g /apoln
       ig0=2
       DO j=2,jjm
          DO i = 1, iim
             pvervel(ig0,l) = pw(i,j,l) * g * unsaire_2d(i,j)
             ig0 = ig0 + 1
          ENDDO
       ENDDO
       pvervel(ig0,l)=pw(1,jjm + 1,l) * g /apols
    ENDDO

    !   45. champ u:

    DO  l=1,llm

       DO  j=2,jjm
          ig0 = 1+(j-2)*iim
          zufi(ig0+1,l)= 0.5 *  &
               (pucov(iim,j,l)/cu_2d(iim,j) + pucov(1,j,l)/cu_2d(1,j))
          pcvgu(ig0+1,l)= 0.5 *  &
               (pducov(iim,j,l)/cu_2d(iim,j) + pducov(1,j,l)/cu_2d(1,j))
          DO i=2,iim
             zufi(ig0+i,l)= 0.5 * &
                  (pucov(i-1,j,l)/cu_2d(i-1,j) &
                  + pucov(i,j,l)/cu_2d(i,j))
             pcvgu(ig0+i,l)= 0.5 * &
                  (pducov(i-1,j,l)/cu_2d(i-1,j) &
                  + pducov(i,j,l)/cu_2d(i,j))
          end DO
       end DO

    end DO

    !   46.champ v:

    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_2d(i,j-1) &
                  + pvcov(i,j,l)/cv_2d(i,j))
             pcvgv(ig0+i,l)= 0.5 * &
                  (pdvcov(i,j-1,l)/cv_2d(i,j-1) &
                  + pdvcov(i,j,l)/cv_2d(i,j))
          ENDDO
       ENDDO
    ENDDO

    !   47. champs de vents aux pole nord   
    !        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
    !        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_2d(1,1)
       z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,1,l)/cv_2d(1,1)
       DO i=2,iim
          z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv_2d(i,1)
          z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,1,l)/cv_2d(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

    ENDDO

    !   48. champs de vents aux pole sud:
    !        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
    !        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_2d(1,jjm)
       z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,jjm,l) &
            /cv_2d(1,jjm)
       DO i=2,iim
          z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv_2d(i,jjm)
          z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,jjm,l)/cv_2d(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(klon,l)  = SSUM(iim,zcos,1)/pi
       pcvgu(klon,l) = SSUM(iim,zcosbis,1)/pi
       zvfi(klon,l)  = SSUM(iim,zsin,1)/pi
       pcvgv(klon,l) = SSUM(iim,zsinbis,1)/pi

    ENDDO

    !IM calcul PV a teta=350, 380, 405K
    CALL PVtheta(klon,llm,pucov,pvcov,pteta, &
         ztfi,zplay,zplev, &
         ntetaSTD,rtetaSTD,PVteta)

    !   Appel de la physique:

    CALL physiq(nq, firstcal, lafin, rdayvrai, heure, dtphys, &
         zplev, zplay, zphi, zphis, presnivs, clesphy0, zufi, zvfi, &
         ztfi, zqfi, pvervel, zdufi, zdvfi, zdtfi, zdqfi, zdpsrf, pducov, &
         PVteta) ! IM diagnostique PVteta, Amip2

    !   transformation des tendances physiques en tendances dynamiques:

    !  tendance sur la pression :

    pdpsfi = gr_fi_dyn(zdpsrf)

    !   62. enthalpie potentielle

    DO l=1,llm

       DO i=1,iim + 1
          pdhfi(i,1,l)    = cpp *  zdtfi(1,l)      / ppk(i, 1  ,l)
          pdhfi(i,jjm + 1,l) = cpp *  zdtfi(klon,l)/ ppk(i,jjm + 1,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(iim + 1,j,l) =  pdhfi(1,j,l)
       ENDDO

    ENDDO

    !   62. humidite specifique

    DO iq=1,nqmx
       DO l=1,llm
          DO i=1,iim + 1
             pdqfi(i,1,l,iq)    = zdqfi(1,l,iq)
             pdqfi(i,jjm + 1,l,iq) = zdqfi(klon,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(iim + 1,j,l,iq) = pdqfi(1,j,l,iq)
          ENDDO
       ENDDO
    ENDDO

    !   63. traceurs

    !     initialisation des tendances
    pdqfi=0.

    DO iq=1,nq
       iiq=niadv(iq)
       DO l=1,llm
          DO i=1,iim + 1
             pdqfi(i,1,l,iiq)    = zdqfi(1,l,iq)
             pdqfi(i,jjm + 1,l,iiq) = zdqfi(klon,l,iq)
          ENDDO
          DO j=2,jjm
             ig0=1+(j-2)*iim
             DO i=1,iim
                pdqfi(i,j,l,iiq) = zdqfi(ig0+i,l,iq)
             ENDDO
             pdqfi(iim + 1,j,l,iiq) = pdqfi(1,j,l,iq)
          ENDDO
       ENDDO
    ENDDO

    !   65. champ u:

    DO l=1,llm

       DO i=1,iim + 1
          pdufi(i,1,l)    = 0.
          pdufi(i,jjm + 1,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_2d(i,j)
          ENDDO
          pdufi(iim,j,l)= &
               0.5*(zdufi(ig0+1,l)+zdufi(ig0+iim,l))*cu_2d(iim,j)
          pdufi(iim + 1,j,l)=pdufi(1,j,l)
       ENDDO

    ENDDO

    !   67. champ v:

    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_2d(i,j)
          ENDDO
          pdvfi(iim + 1,j,l) = pdvfi(1,j,l)
       ENDDO
    ENDDO

    !   68. champ v pres des poles:
    !      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(klon,l)*COS(rlonv(i)) &
               +zdvfi(klon,l)*SIN(rlonv(i))
          pdvfi(i,1,l)= &
               0.5*(pdvfi(i,1,l)+zdvfi(i+1,l))*cv_2d(i,1)
          pdvfi(i,jjm,l)= &
               0.5*(pdvfi(i,jjm,l)+zdvfi(klon-iim-1+i,l))*cv_2d(i,jjm)
       ENDDO

       pdvfi(iim + 1,1,l)  = pdvfi(1,1,l)
       pdvfi(iim + 1,jjm,l)= pdvfi(1,jjm,l)

    ENDDO

    firstcal = .FALSE.

  END SUBROUTINE calfis

end module calfis_m
1	guez	3	module calfis_m
2
3			! Clean: no C preprocessor directive, no include line
4
5			IMPLICIT NONE
6
7			contains
8
9			SUBROUTINE calfis(nq, lafin, rdayvrai, heure, pucov, pvcov, pteta, pq, &
10	guez	10	pmasse, pps, ppk, pphis, pphi, pducov, pdvcov, pdteta, pdq, pw, &
11	guez	3	clesphy0, pdufi, pdvfi, pdhfi, pdqfi, pdpsfi)
12
13			! From dyn3d/calfis.F,v 1.3 2005/05/25 13:10:09
14
15			! Auteurs : P. Le Van, F. Hourdin
16
17			! 1. rearrangement des tableaux et transformation
18			! variables dynamiques > variables physiques
19			! 2. calcul des termes physiques
20			! 3. retransformation des tendances physiques en tendances dynamiques
21
22			! remarques:
23			! ----------
24
25			! - les vents sont donnes dans la physique par leurs composantes
26			! naturelles.
27			! - la variable thermodynamique de la physique est une variable
28			! intensive : T
29			! pour la dynamique on prend T * (preff / p(l)) **kappa
30			! - les deux seules variables dependant de la geometrie necessaires
31			! pour la physique sont la latitude pour le rayonnement et
32			! l'aire de la maille quand on veut integrer une grandeur
33			! horizontalement.
34
35			! Input :
36			! -------
37			! pucov covariant zonal velocity
38			! pvcov covariant meridional velocity
39			! pteta potential temperature
40			! pps surface pressure
41			! pmasse masse d'air dans chaque maille
42			! pts surface temperature (K)
43			! callrad clef d'appel au rayonnement
44
45			! Output :
46			! --------
47			! pdufi tendency for the natural zonal velocity (ms-1)
48			! pdvfi tendency for the natural meridional velocity
49			! pdhfi tendency for the potential temperature
50			! pdtsfi tendency for the surface temperature
51
52			! pdtrad radiative tendencies \ both input
53			! pfluxrad radiative fluxes / and output
54
55			use dimens_m, only: iim, jjm, llm, nqmx
56			use dimphy, only: klon
57			use comconst, only: kappa, cpp, dtphys, g, pi
58			use comvert, only: preff, presnivs
59			use comgeom, only: apoln, cu_2d, cv_2d, unsaire_2d, apols, rlonu, rlonv
60			use advtrac_m, only: niadv
61			use grid_change, only: dyn_phy, gr_fi_dyn
62			use physiq_m, only: physiq
63	guez	10	use pressure_var, only: p3d, pls
64	guez	3
65			! 0. Declarations :
66
67			INTEGER nq
68
69			! Arguments :
70
71			LOGICAL, intent(in):: lafin
72			REAL, intent(in):: heure ! heure de la journée en fraction de jour
73
74			REAL pvcov(iim + 1,jjm,llm)
75			REAL pucov(iim + 1,jjm + 1,llm)
76			REAL pteta(iim + 1,jjm + 1,llm)
77			REAL pmasse(iim + 1,jjm + 1,llm)
78
79			REAL, intent(in):: pq(iim + 1,jjm + 1,llm,nqmx)
80			! (mass fractions of advected fields)
81
82			REAL pphis(iim + 1,jjm + 1)
83			REAL pphi(iim + 1,jjm + 1,llm)
84
85			REAL pdvcov(iim + 1,jjm,llm)
86			REAL pducov(iim + 1,jjm + 1,llm)
87			REAL pdteta(iim + 1,jjm + 1,llm)
88			REAL pdq(iim + 1,jjm + 1,llm,nqmx)
89
90			REAL pw(iim + 1,jjm + 1,llm)
91
92			REAL pps(iim + 1,jjm + 1)
93	guez	10	REAL, intent(in):: ppk(iim + 1,jjm + 1,llm)
94	guez	3
95			REAL pdvfi(iim + 1,jjm,llm)
96			REAL pdufi(iim + 1,jjm + 1,llm)
97			REAL pdhfi(iim + 1,jjm + 1,llm)
98			REAL pdqfi(iim + 1,jjm + 1,llm,nqmx)
99			REAL pdpsfi(iim + 1,jjm + 1)
100
101			INTEGER, PARAMETER:: longcles = 20
102	guez	12	REAL, intent(in):: clesphy0(longcles)
103	guez	3
104			! Local variables :
105
106			INTEGER i,j,l,ig0,ig,iq,iiq
107			REAL zpsrf(klon)
108			REAL zplev(klon,llm+1),zplay(klon,llm)
109			REAL zphi(klon,llm),zphis(klon)
110
111			REAL zufi(klon,llm), zvfi(klon,llm)
112			REAL ztfi(klon,llm) ! temperature
113			real zqfi(klon,llm,nqmx) ! mass fractions of advected fields
114
115			REAL pcvgu(klon,llm), pcvgv(klon,llm)
116			REAL pcvgt(klon,llm), pcvgq(klon,llm,2)
117
118			REAL pvervel(klon,llm)
119
120			REAL zdufi(klon,llm),zdvfi(klon,llm)
121			REAL zdtfi(klon,llm),zdqfi(klon,llm,nqmx)
122			REAL zdpsrf(klon)
123
124			REAL zsin(iim),zcos(iim),z1(iim)
125			REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
126	guez	10	REAL pksurcp(iim + 1,jjm + 1)
127	guez	3
128			! I. Musat: diagnostic PVteta, Amip2
129			INTEGER, PARAMETER:: ntetaSTD=3
130			REAL:: rtetaSTD(ntetaSTD) = (/350., 380., 405./)
131			REAL PVteta(klon,ntetaSTD)
132
133			REAL SSUM
134
135			LOGICAL:: firstcal = .true.
136	guez	7	REAL, intent(in):: rdayvrai
137	guez	3
138			!-----------------------------------------------------------------------
139
140			!!print *, "Call sequence information: calfis"
141
142			! 1. Initialisations :
143			! latitude, longitude et aires des mailles pour la physique:
144
145			! 40. transformation des variables dynamiques en variables physiques:
146			! 41. pressions au sol (en Pascals)
147
148			zpsrf(1) = pps(1,1)
149
150			ig0 = 2
151			DO j = 2,jjm
152			CALL SCOPY(iim,pps(1,j),1,zpsrf(ig0), 1)
153			ig0 = ig0+iim
154			ENDDO
155
156			zpsrf(klon) = pps(1,jjm + 1)
157
158			! 42. pression intercouches :
159
160			! .... zplev definis aux (llm +1) interfaces des couches ....
161			! .... zplay definis aux (llm) milieux des couches ....
162
163			! ... Exner = cp * (p(l) / preff) ** kappa ....
164
165	guez	10	forall (l = 1: llm+1) zplev(:, l) = pack(p3d(:, :, l), dyn_phy)
166	guez	3
167			! 43. temperature naturelle (en K) et pressions milieux couches .
168			DO l=1,llm
169	guez	10	pksurcp = ppk(:, :, l) / cpp
170			pls(:, :, l) = preff * pksurcp**(1./ kappa)
171			zplay(:, l) = pack(pls(:, :, l), dyn_phy)
172			ztfi(:, l) = pack(pteta(:, :, l) * pksurcp, dyn_phy)
173			pcvgt(:, l) = pack(pdteta(:, :, l) * pksurcp / pmasse(:, :, l), dyn_phy)
174	guez	3	ENDDO
175
176			! 43.bis traceurs
177
178			DO iq=1,nq
179			iiq=niadv(iq)
180			DO l=1,llm
181			zqfi(1,l,iq) = pq(1,1,l,iiq)
182			ig0 = 2
183			DO j=2,jjm
184			DO i = 1, iim
185			zqfi(ig0,l,iq) = pq(i,j,l,iiq)
186			ig0 = ig0 + 1
187			ENDDO
188			ENDDO
189			zqfi(ig0,l,iq) = pq(1,jjm + 1,l,iiq)
190			ENDDO
191			ENDDO
192
193			! convergence dynamique pour les traceurs "EAU"
194
195			DO iq=1,2
196			DO l=1,llm
197			pcvgq(1,l,iq)= pdq(1,1,l,iq) / pmasse(1,1,l)
198			ig0 = 2
199			DO j=2,jjm
200			DO i = 1, iim
201			pcvgq(ig0,l,iq) = pdq(i,j,l,iq) / pmasse(i,j,l)
202			ig0 = ig0 + 1
203			ENDDO
204			ENDDO
205			pcvgq(ig0,l,iq)= pdq(1,jjm + 1,l,iq) / pmasse(1,jjm + 1,l)
206			ENDDO
207			ENDDO
208
209			! Geopotentiel calcule par rapport a la surface locale:
210
211			forall (l = 1:llm) zphi(:, l) = pack(pphi(:, :, l), dyn_phy)
212			zphis = pack(pphis, dyn_phy)
213			DO l=1,llm
214			DO ig=1,klon
215			zphi(ig,l)=zphi(ig,l)-zphis(ig)
216			ENDDO
217			ENDDO
218
219			! .... Calcul de la vitesse verticale (en Pams ou Kg/s) ....
220
221			DO l=1,llm
222			pvervel(1,l)=pw(1,1,l) * g /apoln
223			ig0=2
224			DO j=2,jjm
225			DO i = 1, iim
226			pvervel(ig0,l) = pw(i,j,l) * g * unsaire_2d(i,j)
227			ig0 = ig0 + 1
228			ENDDO
229			ENDDO
230			pvervel(ig0,l)=pw(1,jjm + 1,l) * g /apols
231			ENDDO
232
233			! 45. champ u:
234
235			DO l=1,llm
236
237			DO j=2,jjm
238			ig0 = 1+(j-2)*iim
239			zufi(ig0+1,l)= 0.5 * &
240			(pucov(iim,j,l)/cu_2d(iim,j) + pucov(1,j,l)/cu_2d(1,j))
241			pcvgu(ig0+1,l)= 0.5 * &
242			(pducov(iim,j,l)/cu_2d(iim,j) + pducov(1,j,l)/cu_2d(1,j))
243			DO i=2,iim
244			zufi(ig0+i,l)= 0.5 * &
245			(pucov(i-1,j,l)/cu_2d(i-1,j) &
246			+ pucov(i,j,l)/cu_2d(i,j))
247			pcvgu(ig0+i,l)= 0.5 * &
248			(pducov(i-1,j,l)/cu_2d(i-1,j) &
249			+ pducov(i,j,l)/cu_2d(i,j))
250			end DO
251			end DO
252
253			end DO
254
255			! 46.champ v:
256
257			DO l=1,llm
258			DO j=2,jjm
259			ig0=1+(j-2)*iim
260			DO i=1,iim
261			zvfi(ig0+i,l)= 0.5 * &
262			(pvcov(i,j-1,l)/cv_2d(i,j-1) &
263			+ pvcov(i,j,l)/cv_2d(i,j))
264			pcvgv(ig0+i,l)= 0.5 * &
265			(pdvcov(i,j-1,l)/cv_2d(i,j-1) &
266			+ pdvcov(i,j,l)/cv_2d(i,j))
267			ENDDO
268			ENDDO
269			ENDDO
270
271			! 47. champs de vents aux pole nord
272			! U = 1 / pi * integrale [ v * cos(long) * d long ]
273			! V = 1 / pi * integrale [ v * sin(long) * d long ]
274
275			DO l=1,llm
276
277			z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,1,l)/cv_2d(1,1)
278			z1bis(1)=(rlonu(1)-rlonu(iim)+2.pi)pdvcov(1,1,l)/cv_2d(1,1)
279			DO i=2,iim
280			z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv_2d(i,1)
281			z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,1,l)/cv_2d(i,1)
282			ENDDO
283
284			DO i=1,iim
285			zcos(i) = COS(rlonv(i))*z1(i)
286			zcosbis(i)= COS(rlonv(i))*z1bis(i)
287			zsin(i) = SIN(rlonv(i))*z1(i)
288			zsinbis(i)= SIN(rlonv(i))*z1bis(i)
289			ENDDO
290
291			zufi(1,l) = SSUM(iim,zcos,1)/pi
292			pcvgu(1,l) = SSUM(iim,zcosbis,1)/pi
293			zvfi(1,l) = SSUM(iim,zsin,1)/pi
294			pcvgv(1,l) = SSUM(iim,zsinbis,1)/pi
295
296			ENDDO
297
298			! 48. champs de vents aux pole sud:
299			! U = 1 / pi * integrale [ v * cos(long) * d long ]
300			! V = 1 / pi * integrale [ v * sin(long) * d long ]
301
302			DO l=1,llm
303
304			z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,jjm,l) &
305			/cv_2d(1,jjm)
306			z1bis(1)=(rlonu(1)-rlonu(iim)+2.pi)pdvcov(1,jjm,l) &
307			/cv_2d(1,jjm)
308			DO i=2,iim
309			z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv_2d(i,jjm)
310			z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,jjm,l)/cv_2d(i,jjm)
311			ENDDO
312
313			DO i=1,iim
314			zcos(i) = COS(rlonv(i))*z1(i)
315			zcosbis(i) = COS(rlonv(i))*z1bis(i)
316			zsin(i) = SIN(rlonv(i))*z1(i)
317			zsinbis(i) = SIN(rlonv(i))*z1bis(i)
318			ENDDO
319
320			zufi(klon,l) = SSUM(iim,zcos,1)/pi
321			pcvgu(klon,l) = SSUM(iim,zcosbis,1)/pi
322			zvfi(klon,l) = SSUM(iim,zsin,1)/pi
323			pcvgv(klon,l) = SSUM(iim,zsinbis,1)/pi
324
325			ENDDO
326
327			!IM calcul PV a teta=350, 380, 405K
328			CALL PVtheta(klon,llm,pucov,pvcov,pteta, &
329			ztfi,zplay,zplev, &
330			ntetaSTD,rtetaSTD,PVteta)
331
332			! Appel de la physique:
333
334			CALL physiq(nq, firstcal, lafin, rdayvrai, heure, dtphys, &
335			zplev, zplay, zphi, zphis, presnivs, clesphy0, zufi, zvfi, &
336			ztfi, zqfi, pvervel, zdufi, zdvfi, zdtfi, zdqfi, zdpsrf, pducov, &
337			PVteta) ! IM diagnostique PVteta, Amip2
338
339			! transformation des tendances physiques en tendances dynamiques:
340
341			! tendance sur la pression :
342
343			pdpsfi = gr_fi_dyn(zdpsrf)
344
345			! 62. enthalpie potentielle
346
347			DO l=1,llm
348
349			DO i=1,iim + 1
350			pdhfi(i,1,l) = cpp * zdtfi(1,l) / ppk(i, 1 ,l)
351			pdhfi(i,jjm + 1,l) = cpp * zdtfi(klon,l)/ ppk(i,jjm + 1,l)
352			ENDDO
353
354			DO j=2,jjm
355			ig0=1+(j-2)*iim
356			DO i=1,iim
357			pdhfi(i,j,l) = cpp * zdtfi(ig0+i,l) / ppk(i,j,l)
358			ENDDO
359			pdhfi(iim + 1,j,l) = pdhfi(1,j,l)
360			ENDDO
361
362			ENDDO
363
364			! 62. humidite specifique
365
366			DO iq=1,nqmx
367			DO l=1,llm
368			DO i=1,iim + 1
369			pdqfi(i,1,l,iq) = zdqfi(1,l,iq)
370			pdqfi(i,jjm + 1,l,iq) = zdqfi(klon,l,iq)
371			ENDDO
372			DO j=2,jjm
373			ig0=1+(j-2)*iim
374			DO i=1,iim
375			pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,iq)
376			ENDDO
377			pdqfi(iim + 1,j,l,iq) = pdqfi(1,j,l,iq)
378			ENDDO
379			ENDDO
380			ENDDO
381
382			! 63. traceurs
383
384			! initialisation des tendances
385			pdqfi=0.
386
387			DO iq=1,nq
388			iiq=niadv(iq)
389			DO l=1,llm
390			DO i=1,iim + 1
391			pdqfi(i,1,l,iiq) = zdqfi(1,l,iq)
392			pdqfi(i,jjm + 1,l,iiq) = zdqfi(klon,l,iq)
393			ENDDO
394			DO j=2,jjm
395			ig0=1+(j-2)*iim
396			DO i=1,iim
397			pdqfi(i,j,l,iiq) = zdqfi(ig0+i,l,iq)
398			ENDDO
399			pdqfi(iim + 1,j,l,iiq) = pdqfi(1,j,l,iq)
400			ENDDO
401			ENDDO
402			ENDDO
403
404			! 65. champ u:
405
406			DO l=1,llm
407
408			DO i=1,iim + 1
409			pdufi(i,1,l) = 0.
410			pdufi(i,jjm + 1,l) = 0.
411			ENDDO
412
413			DO j=2,jjm
414			ig0=1+(j-2)*iim
415			DO i=1,iim-1
416			pdufi(i,j,l)= &
417			0.5(zdufi(ig0+i,l)+zdufi(ig0+i+1,l))cu_2d(i,j)
418			ENDDO
419			pdufi(iim,j,l)= &
420			0.5(zdufi(ig0+1,l)+zdufi(ig0+iim,l))cu_2d(iim,j)
421			pdufi(iim + 1,j,l)=pdufi(1,j,l)
422			ENDDO
423
424			ENDDO
425
426			! 67. champ v:
427
428			DO l=1,llm
429
430			DO j=2,jjm-1
431			ig0=1+(j-2)*iim
432			DO i=1,iim
433			pdvfi(i,j,l)= &
434			0.5(zdvfi(ig0+i,l)+zdvfi(ig0+i+iim,l))cv_2d(i,j)
435			ENDDO
436			pdvfi(iim + 1,j,l) = pdvfi(1,j,l)
437			ENDDO
438			ENDDO
439
440			! 68. champ v pres des poles:
441			! v = U * cos(long) + V * SIN(long)
442
443			DO l=1,llm
444
445			DO i=1,iim
446			pdvfi(i,1,l)= &
447			zdufi(1,l)COS(rlonv(i))+zdvfi(1,l)SIN(rlonv(i))
448			pdvfi(i,jjm,l)=zdufi(klon,l)*COS(rlonv(i)) &
449			+zdvfi(klon,l)*SIN(rlonv(i))
450			pdvfi(i,1,l)= &
451			0.5(pdvfi(i,1,l)+zdvfi(i+1,l))cv_2d(i,1)
452			pdvfi(i,jjm,l)= &
453			0.5(pdvfi(i,jjm,l)+zdvfi(klon-iim-1+i,l))cv_2d(i,jjm)
454			ENDDO
455
456			pdvfi(iim + 1,1,l) = pdvfi(1,1,l)
457			pdvfi(iim + 1,jjm,l)= pdvfi(1,jjm,l)
458
459			ENDDO
460
461			firstcal = .FALSE.
462
463			END SUBROUTINE calfis
464
465			end module calfis_m