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, pp, 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

    !    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 pp(iim + 1,jjm + 1,llm + 1)
    REAL 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 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 unskap, pksurcp

    ! 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     ....

    unskap   = 1./ kappa

    DO l = 1, llm + 1
       zplev(1,l) = pp(1,1,l)
       ig0 = 2
       DO j = 2, jjm
          DO i =1, iim
             zplev(ig0,l) = pp(i,j,l)
             ig0 = ig0 +1
          ENDDO
       ENDDO
       zplev(klon,l) = pp(1,jjm + 1,l)
    ENDDO

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

    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,jjm + 1,l) / cpp
       zplay(ig0,l)  = preff * pksurcp ** unskap
       ztfi (ig0,l)  = pteta(1,jjm + 1,l)  * pksurcp
       pcvgt(ig0,l)  = pdteta(1,jjm + 1,l) * pksurcp/ pmasse(1,jjm + 1,l)

    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			pmasse, pps, pp, ppk, pphis, pphi, pducov, pdvcov, pdteta, pdq, pw, &
11			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
64			! 0. Declarations :
65
66			INTEGER nq
67
68			! Arguments :
69
70			LOGICAL, intent(in):: lafin
71			REAL, intent(in):: heure ! heure de la journée en fraction de jour
72
73			REAL pvcov(iim + 1,jjm,llm)
74			REAL pucov(iim + 1,jjm + 1,llm)
75			REAL pteta(iim + 1,jjm + 1,llm)
76			REAL pmasse(iim + 1,jjm + 1,llm)
77
78			REAL, intent(in):: pq(iim + 1,jjm + 1,llm,nqmx)
79			! (mass fractions of advected fields)
80
81			REAL pphis(iim + 1,jjm + 1)
82			REAL pphi(iim + 1,jjm + 1,llm)
83
84			REAL pdvcov(iim + 1,jjm,llm)
85			REAL pducov(iim + 1,jjm + 1,llm)
86			REAL pdteta(iim + 1,jjm + 1,llm)
87			REAL pdq(iim + 1,jjm + 1,llm,nqmx)
88
89			REAL pw(iim + 1,jjm + 1,llm)
90
91			REAL pps(iim + 1,jjm + 1)
92			REAL pp(iim + 1,jjm + 1,llm + 1)
93			REAL ppk(iim + 1,jjm + 1,llm)
94
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			REAL clesphy0(longcles)
103
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			REAL unskap, pksurcp
127
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			unskap = 1./ kappa
166
167			DO l = 1, llm + 1
168			zplev(1,l) = pp(1,1,l)
169			ig0 = 2
170			DO j = 2, jjm
171			DO i =1, iim
172			zplev(ig0,l) = pp(i,j,l)
173			ig0 = ig0 +1
174			ENDDO
175			ENDDO
176			zplev(klon,l) = pp(1,jjm + 1,l)
177			ENDDO
178
179			! 43. temperature naturelle (en K) et pressions milieux couches .
180
181			DO l=1,llm
182
183			pksurcp = ppk(1,1,l) / cpp
184			zplay(1,l) = preff * pksurcp ** unskap
185			ztfi(1,l) = pteta(1,1,l) * pksurcp
186			pcvgt(1,l) = pdteta(1,1,l) * pksurcp / pmasse(1,1,l)
187			ig0 = 2
188
189			DO j = 2, jjm
190			DO i = 1, iim
191			pksurcp = ppk(i,j,l) / cpp
192			zplay(ig0,l) = preff * pksurcp ** unskap
193			ztfi(ig0,l) = pteta(i,j,l) * pksurcp
194			pcvgt(ig0,l) = pdteta(i,j,l) * pksurcp / pmasse(i,j,l)
195			ig0 = ig0 + 1
196			ENDDO
197			ENDDO
198
199			pksurcp = ppk(1,jjm + 1,l) / cpp
200			zplay(ig0,l) = preff * pksurcp ** unskap
201			ztfi (ig0,l) = pteta(1,jjm + 1,l) * pksurcp
202			pcvgt(ig0,l) = pdteta(1,jjm + 1,l) * pksurcp/ pmasse(1,jjm + 1,l)
203
204			ENDDO
205
206			! 43.bis traceurs
207
208			DO iq=1,nq
209			iiq=niadv(iq)
210			DO l=1,llm
211			zqfi(1,l,iq) = pq(1,1,l,iiq)
212			ig0 = 2
213			DO j=2,jjm
214			DO i = 1, iim
215			zqfi(ig0,l,iq) = pq(i,j,l,iiq)
216			ig0 = ig0 + 1
217			ENDDO
218			ENDDO
219			zqfi(ig0,l,iq) = pq(1,jjm + 1,l,iiq)
220			ENDDO
221			ENDDO
222
223			! convergence dynamique pour les traceurs "EAU"
224
225			DO iq=1,2
226			DO l=1,llm
227			pcvgq(1,l,iq)= pdq(1,1,l,iq) / pmasse(1,1,l)
228			ig0 = 2
229			DO j=2,jjm
230			DO i = 1, iim
231			pcvgq(ig0,l,iq) = pdq(i,j,l,iq) / pmasse(i,j,l)
232			ig0 = ig0 + 1
233			ENDDO
234			ENDDO
235			pcvgq(ig0,l,iq)= pdq(1,jjm + 1,l,iq) / pmasse(1,jjm + 1,l)
236			ENDDO
237			ENDDO
238
239			! Geopotentiel calcule par rapport a la surface locale:
240
241			forall (l = 1:llm) zphi(:, l) = pack(pphi(:, :, l), dyn_phy)
242			zphis = pack(pphis, dyn_phy)
243			DO l=1,llm
244			DO ig=1,klon
245			zphi(ig,l)=zphi(ig,l)-zphis(ig)
246			ENDDO
247			ENDDO
248
249			! .... Calcul de la vitesse verticale (en Pams ou Kg/s) ....
250
251			DO l=1,llm
252			pvervel(1,l)=pw(1,1,l) * g /apoln
253			ig0=2
254			DO j=2,jjm
255			DO i = 1, iim
256			pvervel(ig0,l) = pw(i,j,l) * g * unsaire_2d(i,j)
257			ig0 = ig0 + 1
258			ENDDO
259			ENDDO
260			pvervel(ig0,l)=pw(1,jjm + 1,l) * g /apols
261			ENDDO
262
263			! 45. champ u:
264
265			DO l=1,llm
266
267			DO j=2,jjm
268			ig0 = 1+(j-2)*iim
269			zufi(ig0+1,l)= 0.5 * &
270			(pucov(iim,j,l)/cu_2d(iim,j) + pucov(1,j,l)/cu_2d(1,j))
271			pcvgu(ig0+1,l)= 0.5 * &
272			(pducov(iim,j,l)/cu_2d(iim,j) + pducov(1,j,l)/cu_2d(1,j))
273			DO i=2,iim
274			zufi(ig0+i,l)= 0.5 * &
275			(pucov(i-1,j,l)/cu_2d(i-1,j) &
276			+ pucov(i,j,l)/cu_2d(i,j))
277			pcvgu(ig0+i,l)= 0.5 * &
278			(pducov(i-1,j,l)/cu_2d(i-1,j) &
279			+ pducov(i,j,l)/cu_2d(i,j))
280			end DO
281			end DO
282
283			end DO
284
285			! 46.champ v:
286
287			DO l=1,llm
288			DO j=2,jjm
289			ig0=1+(j-2)*iim
290			DO i=1,iim
291			zvfi(ig0+i,l)= 0.5 * &
292			(pvcov(i,j-1,l)/cv_2d(i,j-1) &
293			+ pvcov(i,j,l)/cv_2d(i,j))
294			pcvgv(ig0+i,l)= 0.5 * &
295			(pdvcov(i,j-1,l)/cv_2d(i,j-1) &
296			+ pdvcov(i,j,l)/cv_2d(i,j))
297			ENDDO
298			ENDDO
299			ENDDO
300
301			! 47. champs de vents aux pole nord
302			! U = 1 / pi * integrale [ v * cos(long) * d long ]
303			! V = 1 / pi * integrale [ v * sin(long) * d long ]
304
305			DO l=1,llm
306
307			z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,1,l)/cv_2d(1,1)
308			z1bis(1)=(rlonu(1)-rlonu(iim)+2.pi)pdvcov(1,1,l)/cv_2d(1,1)
309			DO i=2,iim
310			z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv_2d(i,1)
311			z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,1,l)/cv_2d(i,1)
312			ENDDO
313
314			DO i=1,iim
315			zcos(i) = COS(rlonv(i))*z1(i)
316			zcosbis(i)= COS(rlonv(i))*z1bis(i)
317			zsin(i) = SIN(rlonv(i))*z1(i)
318			zsinbis(i)= SIN(rlonv(i))*z1bis(i)
319			ENDDO
320
321			zufi(1,l) = SSUM(iim,zcos,1)/pi
322			pcvgu(1,l) = SSUM(iim,zcosbis,1)/pi
323			zvfi(1,l) = SSUM(iim,zsin,1)/pi
324			pcvgv(1,l) = SSUM(iim,zsinbis,1)/pi
325
326			ENDDO
327
328			! 48. champs de vents aux pole sud:
329			! U = 1 / pi * integrale [ v * cos(long) * d long ]
330			! V = 1 / pi * integrale [ v * sin(long) * d long ]
331
332			DO l=1,llm
333
334			z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,jjm,l) &
335			/cv_2d(1,jjm)
336			z1bis(1)=(rlonu(1)-rlonu(iim)+2.pi)pdvcov(1,jjm,l) &
337			/cv_2d(1,jjm)
338			DO i=2,iim
339			z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv_2d(i,jjm)
340			z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,jjm,l)/cv_2d(i,jjm)
341			ENDDO
342
343			DO i=1,iim
344			zcos(i) = COS(rlonv(i))*z1(i)
345			zcosbis(i) = COS(rlonv(i))*z1bis(i)
346			zsin(i) = SIN(rlonv(i))*z1(i)
347			zsinbis(i) = SIN(rlonv(i))*z1bis(i)
348			ENDDO
349
350			zufi(klon,l) = SSUM(iim,zcos,1)/pi
351			pcvgu(klon,l) = SSUM(iim,zcosbis,1)/pi
352			zvfi(klon,l) = SSUM(iim,zsin,1)/pi
353			pcvgv(klon,l) = SSUM(iim,zsinbis,1)/pi
354
355			ENDDO
356
357			!IM calcul PV a teta=350, 380, 405K
358			CALL PVtheta(klon,llm,pucov,pvcov,pteta, &
359			ztfi,zplay,zplev, &
360			ntetaSTD,rtetaSTD,PVteta)
361
362			! Appel de la physique:
363
364			CALL physiq(nq, firstcal, lafin, rdayvrai, heure, dtphys, &
365			zplev, zplay, zphi, zphis, presnivs, clesphy0, zufi, zvfi, &
366			ztfi, zqfi, pvervel, zdufi, zdvfi, zdtfi, zdqfi, zdpsrf, pducov, &
367			PVteta) ! IM diagnostique PVteta, Amip2
368
369			! transformation des tendances physiques en tendances dynamiques:
370
371			! tendance sur la pression :
372
373			pdpsfi = gr_fi_dyn(zdpsrf)
374
375			! 62. enthalpie potentielle
376
377			DO l=1,llm
378
379			DO i=1,iim + 1
380			pdhfi(i,1,l) = cpp * zdtfi(1,l) / ppk(i, 1 ,l)
381			pdhfi(i,jjm + 1,l) = cpp * zdtfi(klon,l)/ ppk(i,jjm + 1,l)
382			ENDDO
383
384			DO j=2,jjm
385			ig0=1+(j-2)*iim
386			DO i=1,iim
387			pdhfi(i,j,l) = cpp * zdtfi(ig0+i,l) / ppk(i,j,l)
388			ENDDO
389			pdhfi(iim + 1,j,l) = pdhfi(1,j,l)
390			ENDDO
391
392			ENDDO
393
394			! 62. humidite specifique
395
396			DO iq=1,nqmx
397			DO l=1,llm
398			DO i=1,iim + 1
399			pdqfi(i,1,l,iq) = zdqfi(1,l,iq)
400			pdqfi(i,jjm + 1,l,iq) = zdqfi(klon,l,iq)
401			ENDDO
402			DO j=2,jjm
403			ig0=1+(j-2)*iim
404			DO i=1,iim
405			pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,iq)
406			ENDDO
407			pdqfi(iim + 1,j,l,iq) = pdqfi(1,j,l,iq)
408			ENDDO
409			ENDDO
410			ENDDO
411
412			! 63. traceurs
413
414			! initialisation des tendances
415			pdqfi=0.
416
417			DO iq=1,nq
418			iiq=niadv(iq)
419			DO l=1,llm
420			DO i=1,iim + 1
421			pdqfi(i,1,l,iiq) = zdqfi(1,l,iq)
422			pdqfi(i,jjm + 1,l,iiq) = zdqfi(klon,l,iq)
423			ENDDO
424			DO j=2,jjm
425			ig0=1+(j-2)*iim
426			DO i=1,iim
427			pdqfi(i,j,l,iiq) = zdqfi(ig0+i,l,iq)
428			ENDDO
429			pdqfi(iim + 1,j,l,iiq) = pdqfi(1,j,l,iq)
430			ENDDO
431			ENDDO
432			ENDDO
433
434			! 65. champ u:
435
436			DO l=1,llm
437
438			DO i=1,iim + 1
439			pdufi(i,1,l) = 0.
440			pdufi(i,jjm + 1,l) = 0.
441			ENDDO
442
443			DO j=2,jjm
444			ig0=1+(j-2)*iim
445			DO i=1,iim-1
446			pdufi(i,j,l)= &
447			0.5(zdufi(ig0+i,l)+zdufi(ig0+i+1,l))cu_2d(i,j)
448			ENDDO
449			pdufi(iim,j,l)= &
450			0.5(zdufi(ig0+1,l)+zdufi(ig0+iim,l))cu_2d(iim,j)
451			pdufi(iim + 1,j,l)=pdufi(1,j,l)
452			ENDDO
453
454			ENDDO
455
456			! 67. champ v:
457
458			DO l=1,llm
459
460			DO j=2,jjm-1
461			ig0=1+(j-2)*iim
462			DO i=1,iim
463			pdvfi(i,j,l)= &
464			0.5(zdvfi(ig0+i,l)+zdvfi(ig0+i+iim,l))cv_2d(i,j)
465			ENDDO
466			pdvfi(iim + 1,j,l) = pdvfi(1,j,l)
467			ENDDO
468			ENDDO
469
470			! 68. champ v pres des poles:
471			! v = U * cos(long) + V * SIN(long)
472
473			DO l=1,llm
474
475			DO i=1,iim
476			pdvfi(i,1,l)= &
477			zdufi(1,l)COS(rlonv(i))+zdvfi(1,l)SIN(rlonv(i))
478			pdvfi(i,jjm,l)=zdufi(klon,l)*COS(rlonv(i)) &
479			+zdvfi(klon,l)*SIN(rlonv(i))
480			pdvfi(i,1,l)= &
481			0.5(pdvfi(i,1,l)+zdvfi(i+1,l))cv_2d(i,1)
482			pdvfi(i,jjm,l)= &
483			0.5(pdvfi(i,jjm,l)+zdvfi(klon-iim-1+i,l))cv_2d(i,jjm)
484			ENDDO
485
486			pdvfi(iim + 1,1,l) = pdvfi(1,1,l)
487			pdvfi(iim + 1,jjm,l)= pdvfi(1,jjm,l)
488
489			ENDDO
490
491			firstcal = .FALSE.
492
493			END SUBROUTINE calfis
494
495			end module calfis_m