trunk/dyn3d/calfis.f90

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

    !    0.  Declarations :

    INTEGER, intent(in):: 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

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