trunk/dyn3d/calfis.f

module calfis_m

  ! Clean: no C preprocessor directive, no include line

  IMPLICIT NONE

contains

  SUBROUTINE calfis(lafin, rdayvrai, heure, pucov, pvcov, pteta, q, &
       pmasse, pps, ppk, pphis, pphi, pducov, pdvcov, pdteta, pdq, pw, &
       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
    use comgeom, only: apoln, cu_2d, cv_2d, unsaire_2d, apols, rlonu, rlonv
    use iniadvtrac_m, only: niadv
    use grid_change, only: dyn_phy, gr_fi_dyn
    use physiq_m, only: physiq
    use pressure_var, only: p3d, pls

    !    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):: q(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

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