libf/dyn3d/exner_hyb.f90

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	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	! Authors : P. Le Van, F. Hourdin
12
13	! Calcule la fonction d'Exner :
14	! pk = Cp * p ** kappa
15	! 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
19	! Au sommet de l'atmosphère :
20	! p(llm+1) = 0.
21	! "ps" et "pks" sont la pression et la fonction d'Exner au sol.
22
23	! À partir des relations :
24
25	! -------- z
26	!(1) pdz(pk) = kappa pk * dz(p)
27
28	!(2) pk(l) = alpha(l)+ beta(l) * pk(l-1)
29
30	! (voir note de F. Hourdin), on détermine successivement, du haut
31	! vers le bas des couches, les coefficients :
32	! alpha(llm), beta(llm)..., alpha(l), beta(l)..., alpha(2), beta(2)
33	! puis "pk(:, 1)".
34	! Ensuite, on calcule, du bas vers le haut des couches, "pk(:, l)"
35	! 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	INTEGER l
53	REAL unpl2k, dellta((iim + 1) * (jjm + 1))
54
55	REAL ppn(iim), pps(iim)
56
57	!-------------------------------------
58
59	pks = cpp * (ps / preff)**kappa
60	ppn = aire_2d(:iim, 1) * pks(:iim)
61	pps = aire_2d(:iim, jjm + 1) &
62	* pks(1 + (iim + 1) * jjm: iim + (iim + 1) * jjm)
63	pks(:iim + 1) = SUM(ppn) /apoln
64	pks(1+(iim + 1) * jjm:) = SUM(pps) /apols
65
66	unpl2k = 1. + 2 * kappa
67
68	! Calcul des coefficients alpha et beta pour la couche l = llm :
69	alpha(:, llm) = 0.
70	beta(:, llm) = 1./ unpl2k
71
72	! 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	ENDDO
78
79	! Calcul de pk pour la couche 1, près du sol :
80	pk(:, 1) = (p(:, 1) * pks - 0.5 * alpha(:, 2) * p(:, 2)) &
81	/ (p(:, 1) * (1. + kappa) + 0.5 * (beta(:, 2) - unpl2k) * p(:, 2))
82
83	! Calcul de pk(:, l) pour l = 2 à l = llm :
84	DO l = 2, llm
85	pk(:, l) = alpha(:, l) + beta(:, l) * pk(:, l-1)
86	ENDDO
87
88	if (present(pkf)) then
89	pkf = pk
90	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