trunk/dyn3d/exner_hyb.f

module exner_hyb_m

  IMPLICIT NONE

contains

  SUBROUTINE exner_hyb(ps, p, pks, pk, pkf)

    ! From dyn3d/exner_hyb.F, v 1.1.1.1 2004/05/19 12:53:07

    ! Authors : P. Le Van, F. Hourdin

    ! Calcule la fonction d'Exner :
    ! pk = Cp * p ** kappa
    ! aux milieux des "llm" couches.
    ! "Pk(l)" est calculé au milieu de la couche "l", entre les pressions
    ! "p(l)" et "p(l+1)", définies aux interfaces des couches.

    ! Au sommet de l'atmosphère :
    ! p(llm+1) = 0.
    ! "ps" et "pks" sont la pression et la fonction d'Exner au sol.

    ! À partir des relations :

    !   -------- z
    !(1) p*dz(pk) = kappa * pk * dz(p)

    !(2) pk(l) = alpha(l)+ beta(l) * pk(l-1)

    ! (voir note de F. Hourdin), on détermine successivement, du haut
    ! vers le bas des couches, les coefficients :
    ! alpha(llm), beta(llm)..., alpha(l), beta(l)..., alpha(2), beta(2)
    ! puis "pk(:, 1)".
    ! Ensuite, on calcule, du bas vers le haut des couches, "pk(:, l)"
    ! donné par la relation (2), pour l = 2 à l = llm.

    use dimens_m, only: iim, jjm, llm
    use comconst, only: kappa, cpp
    use comvert, only: preff
    use comgeom, only: aire_2d, apoln, apols

    REAL, intent(in):: ps((iim + 1) * (jjm + 1))
    REAL, intent(in):: p((iim + 1) * (jjm + 1), llm + 1)

    real, intent(out):: pks((iim + 1) * (jjm + 1))
    real, intent(out):: pk((iim + 1) * (jjm + 1), llm)
    real, intent(out), optional:: pkf((iim + 1) * (jjm + 1), llm)

    ! Variables locales

    real alpha((iim + 1) * (jjm + 1), llm), beta((iim + 1) * (jjm + 1), llm)
    INTEGER l
    REAL unpl2k, dellta((iim + 1) * (jjm + 1))

    REAL ppn(iim), pps(iim)

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

    pks = cpp * (ps / preff)**kappa
    ppn = aire_2d(:iim, 1) * pks(:iim)
    pps = aire_2d(:iim, jjm + 1) &
         * pks(1 + (iim + 1) * jjm: iim + (iim + 1) * jjm)
    pks(:iim + 1) = SUM(ppn) /apoln
    pks(1+(iim + 1) * jjm:) = SUM(pps) /apols

    unpl2k = 1. + 2 * kappa

    ! Calcul des coefficients alpha et beta pour la couche l = llm :
    alpha(:, llm) = 0.
    beta(:, llm) = 1./ unpl2k

    ! Calcul des coefficients alpha et beta pour l = llm-1 à l = 2 :
    DO l = llm - 1, 2, -1
       dellta = p(:, l) * unpl2k + p(:, l+1) * (beta(:, l+1) - unpl2k)
       alpha(:, l) = - p(:, l+1) / dellta * alpha(:, l+1)
       beta(:, l) = p(:, l) / dellta   
    ENDDO

    ! Calcul de pk pour la couche 1, près du sol :
    pk(:, 1) = (p(:, 1) * pks - 0.5 * alpha(:, 2) * p(:, 2))  &
         / (p(:, 1) * (1. + kappa) + 0.5 * (beta(:, 2) - unpl2k) * p(:, 2))

    ! Calcul de pk(:, l) pour l = 2 à l = llm :
    DO l = 2, llm
       pk(:, l) = alpha(:, l) + beta(:, l) * pk(:, l-1)
    ENDDO

    if (present(pkf)) then
       pkf = pk
       CALL filtreg(pkf, jjm + 1, llm, 2, 1, .TRUE., 1)
    end if

  END SUBROUTINE exner_hyb

end module exner_hyb_m
1	guez	3	module exner_hyb_m
2
3			IMPLICIT NONE
4
5			contains
6
7			SUBROUTINE exner_hyb(ps, p, pks, pk, pkf)
8
9			! From dyn3d/exner_hyb.F, v 1.1.1.1 2004/05/19 12:53:07
10
11	guez	10	! Authors : P. Le Van, F. Hourdin
12	guez	3
13			! Calcule la fonction d'Exner :
14			! pk = Cp * p ** kappa
15	guez	10	! aux milieux des "llm" couches.
16			! "Pk(l)" est calculé au milieu de la couche "l", entre les pressions
17			! "p(l)" et "p(l+1)", définies aux interfaces des couches.
18	guez	3
19			! Au sommet de l'atmosphère :
20			! p(llm+1) = 0.
21	guez	10	! "ps" et "pks" sont la pression et la fonction d'Exner au sol.
22	guez	3
23			! À partir des relations :
24
25			! -------- z
26	guez	10	!(1) pdz(pk) = kappa pk * dz(p)
27	guez	3
28	guez	10	!(2) pk(l) = alpha(l)+ beta(l) * pk(l-1)
29	guez	3
30	guez	10	! (voir note de F. Hourdin), on détermine successivement, du haut
31	guez	3	! vers le bas des couches, les coefficients :
32			! alpha(llm), beta(llm)..., alpha(l), beta(l)..., alpha(2), beta(2)
33	guez	10	! puis "pk(:, 1)".
34			! Ensuite, on calcule, du bas vers le haut des couches, "pk(:, l)"
35	guez	3	! donné par la relation (2), pour l = 2 à l = llm.
36
37			use dimens_m, only: iim, jjm, llm
38			use comconst, only: kappa, cpp
39			use comvert, only: preff
40			use comgeom, only: aire_2d, apoln, apols
41
42			REAL, intent(in):: ps((iim + 1) * (jjm + 1))
43			REAL, intent(in):: p((iim + 1) * (jjm + 1), llm + 1)
44
45			real, intent(out):: pks((iim + 1) * (jjm + 1))
46			real, intent(out):: pk((iim + 1) * (jjm + 1), llm)
47			real, intent(out), optional:: pkf((iim + 1) * (jjm + 1), llm)
48
49			! Variables locales
50
51			real alpha((iim + 1) * (jjm + 1), llm), beta((iim + 1) * (jjm + 1), llm)
52	guez	10	INTEGER l
53			REAL unpl2k, dellta((iim + 1) * (jjm + 1))
54	guez	3
55			REAL ppn(iim), pps(iim)
56
57			!-------------------------------------
58
59	guez	10	pks = cpp * (ps / preff)**kappa
60			ppn = aire_2d(:iim, 1) * pks(:iim)
61			pps = aire_2d(:iim, jjm + 1) &
62	guez	3	* pks(1 + (iim + 1) * jjm: iim + (iim + 1) * jjm)
63	guez	10	pks(:iim + 1) = SUM(ppn) /apoln
64			pks(1+(iim + 1) * jjm:) = SUM(pps) /apols
65	guez	3
66			unpl2k = 1. + 2 * kappa
67
68	guez	10	! Calcul des coefficients alpha et beta pour la couche l = llm :
69	guez	3	alpha(:, llm) = 0.
70	guez	10	beta(:, llm) = 1./ unpl2k
71	guez	3
72	guez	10	! Calcul des coefficients alpha et beta pour l = llm-1 à l = 2 :
73			DO l = llm - 1, 2, -1
74			dellta = p(:, l) * unpl2k + p(:, l+1) * (beta(:, l+1) - unpl2k)
75			alpha(:, l) = - p(:, l+1) / dellta * alpha(:, l+1)
76			beta(:, l) = p(:, l) / dellta
77	guez	3	ENDDO
78
79			! Calcul de pk pour la couche 1, près du sol :
80	guez	10	pk(:, 1) = (p(:, 1) * pks - 0.5 * alpha(:, 2) * p(:, 2)) &
81	guez	3	/ (p(:, 1) * (1. + kappa) + 0.5 * (beta(:, 2) - unpl2k) * p(:, 2))
82
83	guez	10	! Calcul de pk(:, l) pour l = 2 à l = llm :
84	guez	3	DO l = 2, llm
85	guez	10	pk(:, l) = alpha(:, l) + beta(:, l) * pk(:, l-1)
86	guez	3	ENDDO
87
88			if (present(pkf)) then
89	guez	10	pkf = pk
90	guez	3	CALL filtreg(pkf, jjm + 1, llm, 2, 1, .TRUE., 1)
91			end if
92
93			END SUBROUTINE exner_hyb
94
95			end module exner_hyb_m