trunk/dyn3d/vlspltqs.f

SUBROUTINE vlspltqs(q, pente_max, masse, w, pbaru, pbarv, pdt, p, pk, teta)

  ! From LMDZ4/libf/dyn3d/vlspltqs.F, version 1.2 2005/02/24 12:16:57 fairhead

  ! Authors: P. Le Van, F. Hourdin, F. Forget, F. Codron

  ! Schéma d'advection "pseudo amont"
  ! + test sur humidité spécifique : Q advecté < Qsat aval
  ! (F. Codron, 10/99)

  ! q, pbaru, pbarv, w sont des arguments d'entree pour le sous-programme

  ! pente_max facteur de limitation des pentes: 2 en général
  ! 0 pour un schéma amont
  ! pbaru, pbarv, w flux de masse en u , v , w
  ! pdt pas de temps

  ! teta température potentielle, p pression aux interfaces,
  ! pk exner au milieu des couches nécessaire pour calculer Qsat

  USE dimens_m, ONLY : iim, llm
  USE paramet_m, ONLY : iip1, iip2, ijp1llm, ip1jm, ip1jmp1, llmp1
  USE comconst, ONLY : cpp

  IMPLICIT NONE

  ! Arguments:

  REAL masse(ip1jmp1, llm), pente_max
  REAL, intent(in):: pbaru(ip1jmp1, llm), pbarv(ip1jm, llm)
  REAL q(ip1jmp1, llm)
  REAL w(ip1jmp1, llm)
  real, intent(in):: pdt
  REAL, intent(in):: p(ip1jmp1, llmp1)
  real, intent(in):: teta(ip1jmp1, llm)
  real, intent(in):: pk(ip1jmp1, llm)

  ! Local

  INTEGER i, ij, l, j, ii

  REAL qsat(ip1jmp1, llm)
  REAL zm(ip1jmp1, llm)
  REAL mu(ip1jmp1, llm)
  REAL mv(ip1jm, llm)
  REAL mw(ip1jmp1, llm+1)
  REAL zq(ip1jmp1, llm)
  REAL temps1, temps2, temps3
  REAL zzpbar, zzw
  LOGICAL testcpu
  SAVE testcpu
  SAVE temps1, temps2, temps3

  REAL qmin, qmax
  DATA qmin, qmax/0., 1.e33/
  DATA testcpu/.false./
  DATA temps1, temps2, temps3/0., 0., 0./

  !--pour rapport de melange saturant--

  REAL rtt, retv, r2es, r3les, r3ies, r4les, r4ies, play
  REAL ptarg, pdelarg, foeew, zdelta
  REAL tempe(ip1jmp1)

  ! fonction psat(T)
  FOEEW(PTARG, PDELARG) = EXP ((R3LES*(1.-PDELARG)+R3IES*PDELARG) &
       * (PTARG-RTT) / (PTARG-(R4LES*(1.-PDELARG)+R4IES*PDELARG)))

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

  r2es = 380.11733
  r3les = 17.269
  r3ies = 21.875
  r4les = 35.86
  r4ies = 7.66
  retv = 0.6077667
  rtt = 273.16

  !-- Calcul de Qsat en chaque point
  !-- approximation: au milieu des couches play(l)=(p(l)+p(l+1))/2
  ! pour eviter une exponentielle.
  DO l = 1, llm
     DO ij = 1, ip1jmp1
        tempe(ij) = teta(ij, l) * pk(ij, l) /cpp
     ENDDO
     DO ij = 1, ip1jmp1
        zdelta = MAX(0., SIGN(1., rtt - tempe(ij)))
        play = 0.5*(p(ij, l)+p(ij, l+1))
        qsat(ij, l) = MIN(0.5, r2es* FOEEW(tempe(ij), zdelta) / play)
        qsat(ij, l) = qsat(ij, l) / (1. - retv * qsat(ij, l))
     ENDDO
  ENDDO

  zzpbar = 0.5 * pdt
  zzw = pdt
  DO l=1, llm
     DO ij = iip2, ip1jm
        mu(ij, l)=pbaru(ij, l) * zzpbar
     ENDDO
     DO ij=1, ip1jm
        mv(ij, l)=pbarv(ij, l) * zzpbar
     ENDDO
     DO ij=1, ip1jmp1
        mw(ij, l)=w(ij, l) * zzw
     ENDDO
  ENDDO

  DO ij=1, ip1jmp1
     mw(ij, llm+1)=0.
  ENDDO

  CALL SCOPY(ijp1llm, q, 1, zq, 1)
  CALL SCOPY(ijp1llm, masse, 1, zm, 1)

  ! call minmaxq(zq, qmin, qmax, 'avant vlxqs ')
  call vlxqs(zq, pente_max, zm, mu, qsat)

  ! call minmaxq(zq, qmin, qmax, 'avant vlyqs ')

  call vlyqs(zq, pente_max, zm, mv, qsat)

  ! call minmaxq(zq, qmin, qmax, 'avant vlz ')

  call vlz(zq, pente_max, zm, mw)

  ! call minmaxq(zq, qmin, qmax, 'avant vlyqs ')
  ! call minmaxq(zm, qmin, qmax, 'M avant vlyqs ')

  call vlyqs(zq, pente_max, zm, mv, qsat)

  ! call minmaxq(zq, qmin, qmax, 'avant vlxqs ')
  ! call minmaxq(zm, qmin, qmax, 'M avant vlxqs ')

  call vlxqs(zq, pente_max, zm, mu, qsat)

  ! call minmaxq(zq, qmin, qmax, 'apres vlxqs ')
  ! call minmaxq(zm, qmin, qmax, 'M apres vlxqs ')

  DO l=1, llm
     DO ij=1, ip1jmp1
        q(ij, l)=zq(ij, l)
     ENDDO
     DO ij=1, ip1jm+1, iip1
        q(ij+iim, l)=q(ij, l)
     ENDDO
  ENDDO

END SUBROUTINE vlspltqs
1	guez	44	SUBROUTINE vlspltqs(q, pente_max, masse, w, pbaru, pbarv, pdt, p, pk, teta)
2	guez	3
3	guez	44	! From LMDZ4/libf/dyn3d/vlspltqs.F, version 1.2 2005/02/24 12:16:57 fairhead
4	guez	3
5	guez	44	! Authors: P. Le Van, F. Hourdin, F. Forget, F. Codron
6	guez	3
7	guez	44	! Schéma d'advection "pseudo amont"
8			! + test sur humidité spécifique : Q advecté < Qsat aval
9			! (F. Codron, 10/99)
10	guez	3
11	guez	44	! q, pbaru, pbarv, w sont des arguments d'entree pour le sous-programme
12	guez	3
13	guez	44	! pente_max facteur de limitation des pentes: 2 en général
14			! 0 pour un schéma amont
15			! pbaru, pbarv, w flux de masse en u , v , w
16			! pdt pas de temps
17	guez	3
18	guez	44	! teta température potentielle, p pression aux interfaces,
19			! pk exner au milieu des couches nécessaire pour calculer Qsat
20	guez	3
21	guez	44	USE dimens_m, ONLY : iim, llm
22			USE paramet_m, ONLY : iip1, iip2, ijp1llm, ip1jm, ip1jmp1, llmp1
23			USE comconst, ONLY : cpp
24	guez	3
25	guez	44	IMPLICIT NONE
26	guez	3
27	guez	44	! Arguments:
28	guez	3
29	guez	44	REAL masse(ip1jmp1, llm), pente_max
30			REAL, intent(in):: pbaru(ip1jmp1, llm), pbarv(ip1jm, llm)
31			REAL q(ip1jmp1, llm)
32			REAL w(ip1jmp1, llm)
33			real, intent(in):: pdt
34			REAL, intent(in):: p(ip1jmp1, llmp1)
35			real, intent(in):: teta(ip1jmp1, llm)
36	guez	71	real, intent(in):: pk(ip1jmp1, llm)
37	guez	3
38	guez	44	! Local
39	guez	3
40	guez	44	INTEGER i, ij, l, j, ii
41	guez	3
42	guez	44	REAL qsat(ip1jmp1, llm)
43			REAL zm(ip1jmp1, llm)
44			REAL mu(ip1jmp1, llm)
45			REAL mv(ip1jm, llm)
46			REAL mw(ip1jmp1, llm+1)
47			REAL zq(ip1jmp1, llm)
48			REAL temps1, temps2, temps3
49			REAL zzpbar, zzw
50			LOGICAL testcpu
51			SAVE testcpu
52			SAVE temps1, temps2, temps3
53	guez	3
54	guez	44	REAL qmin, qmax
55			DATA qmin, qmax/0., 1.e33/
56			DATA testcpu/.false./
57			DATA temps1, temps2, temps3/0., 0., 0./
58	guez	3
59	guez	44	!--pour rapport de melange saturant--
60	guez	3
61	guez	44	REAL rtt, retv, r2es, r3les, r3ies, r4les, r4ies, play
62			REAL ptarg, pdelarg, foeew, zdelta
63			REAL tempe(ip1jmp1)
64	guez	3
65	guez	44	! fonction psat(T)
66			FOEEW(PTARG, PDELARG) = EXP ((R3LES(1.-PDELARG)+R3IESPDELARG) &
67			* (PTARG-RTT) / (PTARG-(R4LES(1.-PDELARG)+R4IESPDELARG)))
68	guez	3
69	guez	44	!------------------------------------------------------------------
70	guez	3
71	guez	44	r2es = 380.11733
72			r3les = 17.269
73			r3ies = 21.875
74			r4les = 35.86
75			r4ies = 7.66
76			retv = 0.6077667
77			rtt = 273.16
78	guez	3
79	guez	44	!-- Calcul de Qsat en chaque point
80			!-- approximation: au milieu des couches play(l)=(p(l)+p(l+1))/2
81			! pour eviter une exponentielle.
82			DO l = 1, llm
83			DO ij = 1, ip1jmp1
84			tempe(ij) = teta(ij, l) * pk(ij, l) /cpp
85			ENDDO
86			DO ij = 1, ip1jmp1
87			zdelta = MAX(0., SIGN(1., rtt - tempe(ij)))
88			play = 0.5*(p(ij, l)+p(ij, l+1))
89			qsat(ij, l) = MIN(0.5, r2es* FOEEW(tempe(ij), zdelta) / play)
90			qsat(ij, l) = qsat(ij, l) / (1. - retv * qsat(ij, l))
91			ENDDO
92			ENDDO
93	guez	3
94	guez	44	zzpbar = 0.5 * pdt
95			zzw = pdt
96			DO l=1, llm
97			DO ij = iip2, ip1jm
98			mu(ij, l)=pbaru(ij, l) * zzpbar
99			ENDDO
100			DO ij=1, ip1jm
101			mv(ij, l)=pbarv(ij, l) * zzpbar
102			ENDDO
103			DO ij=1, ip1jmp1
104			mw(ij, l)=w(ij, l) * zzw
105			ENDDO
106			ENDDO
107	guez	3
108	guez	44	DO ij=1, ip1jmp1
109			mw(ij, llm+1)=0.
110			ENDDO
111	guez	3
112	guez	44	CALL SCOPY(ijp1llm, q, 1, zq, 1)
113			CALL SCOPY(ijp1llm, masse, 1, zm, 1)
114	guez	3
115	guez	44	! call minmaxq(zq, qmin, qmax, 'avant vlxqs ')
116			call vlxqs(zq, pente_max, zm, mu, qsat)
117	guez	3
118	guez	44	! call minmaxq(zq, qmin, qmax, 'avant vlyqs ')
119	guez	3
120	guez	44	call vlyqs(zq, pente_max, zm, mv, qsat)
121	guez	3
122	guez	44	! call minmaxq(zq, qmin, qmax, 'avant vlz ')
123	guez	3
124	guez	44	call vlz(zq, pente_max, zm, mw)
125
126			! call minmaxq(zq, qmin, qmax, 'avant vlyqs ')
127			! call minmaxq(zm, qmin, qmax, 'M avant vlyqs ')
128
129			call vlyqs(zq, pente_max, zm, mv, qsat)
130
131			! call minmaxq(zq, qmin, qmax, 'avant vlxqs ')
132			! call minmaxq(zm, qmin, qmax, 'M avant vlxqs ')
133
134			call vlxqs(zq, pente_max, zm, mu, qsat)
135
136			! call minmaxq(zq, qmin, qmax, 'apres vlxqs ')
137			! call minmaxq(zm, qmin, qmax, 'M apres vlxqs ')
138
139			DO l=1, llm
140			DO ij=1, ip1jmp1
141			q(ij, l)=zq(ij, l)
142			ENDDO
143			DO ij=1, ip1jm+1, iip1
144			q(ij+iim, l)=q(ij, l)
145			ENDDO
146			ENDDO
147
148			END SUBROUTINE vlspltqs