trunk/dyn3d/sortvarc.f

SUBROUTINE sortvarc(itau,ucov,teta,ps,masse,pk,phis,vorpot,phi,bern,dp, &
     time)


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

  !   Auteur:    P. Le Van
  !   Objet:
  !   sortie des variables de controle

  USE dimens_m, ONLY : iim, jjm, llm
  USE paramet_m, ONLY : iip1, iip2, ijp1llm, ip1jm, ip1jmp1, jjp1
  USE comconst, ONLY : daysec, dtvr, g, omeg, rad
  USE comgeom, ONLY : aire, cu, rlatu
  USE temps, ONLY : day_ini
  USE ener, ONLY : ang, ang0, etot, etot0, ptot, ptot0, rmsdpdt, rmsv, &
       stot, stot0, ztot, ztot0

  IMPLICIT NONE

  !   Arguments:

  INTEGER, INTENT (IN) :: itau
  REAL :: ucov(ip1jmp1,llm), teta(ip1jmp1,llm), masse(ip1jmp1,llm)
  REAL :: ps(ip1jmp1), phis(ip1jmp1)
  REAL :: vorpot(ip1jm,llm)
  REAL :: phi(ip1jmp1,llm), bern(ip1jmp1,llm)
  REAL :: dp(ip1jmp1)
  REAL :: time
  REAL, INTENT (IN) :: pk(ip1jmp1,llm)

  !   Local:

  REAL :: vor(ip1jm), bernf(ip1jmp1,llm), ztotl(llm)
  REAL :: etotl(llm), stotl(llm), rmsvl(llm), angl(llm), ge(ip1jmp1)
  REAL :: cosphi(ip1jm), omegcosp(ip1jm)
  REAL :: dtvrs1j, rjour, heure, radsg, radomeg
  REAL :: rday, massebxy(ip1jm,llm)
  INTEGER :: l, ij, imjmp1

  REAL :: ssum

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

  dtvrs1j = dtvr/daysec
  rjour = float(int(itau*dtvrs1j))
  heure = (itau*dtvrs1j-rjour)*24.
  imjmp1 = iim*jjp1
  IF (abs(heure-24.)<=0.0001) heure = 0.

  CALL massbarxy(masse,massebxy)

  !   .....  Calcul  de  rmsdpdt  .....

  ge(:) = dp(:)*dp(:)

  rmsdpdt = ssum(ip1jmp1,ge,1) - ssum(jjp1,ge,iip1)

  rmsdpdt = daysec*1.E-2*sqrt(rmsdpdt/imjmp1)

  CALL scopy(ijp1llm,bern,1,bernf,1)
  CALL filtreg(bernf,jjp1,llm,-2,2,.TRUE.,1)

  !   .....  Calcul du moment  angulaire   .....

  radsg = rad/g
  radomeg = rad*omeg

  DO ij = iip2, ip1jm
     cosphi(ij) = cos(rlatu((ij-1)/iip1+1))
     omegcosp(ij) = radomeg*cosphi(ij)
  END DO

  !  ...  Calcul  de l'energie,de l'enstrophie,de l'entropie et de rmsv  .

  DO l = 1, llm
     DO ij = 1, ip1jm
        vor(ij) = vorpot(ij,l)*vorpot(ij,l)*massebxy(ij,l)
     END DO
     ztotl(l) = (ssum(ip1jm,vor,1)-ssum(jjm,vor,iip1))

     DO ij = 1, ip1jmp1
        ge(ij) = masse(ij,l)*(phis(ij)+teta(ij,l)*pk(ij,l)+bernf(ij,l)-phi( &
             ij,l))
     END DO
     etotl(l) = ssum(ip1jmp1,ge,1) - ssum(jjp1,ge,iip1)

     DO ij = 1, ip1jmp1
        ge(ij) = masse(ij,l)*teta(ij,l)
     END DO
     stotl(l) = ssum(ip1jmp1,ge,1) - ssum(jjp1,ge,iip1)

     DO ij = 1, ip1jmp1
        ge(ij) = masse(ij,l)*amax1(bernf(ij,l)-phi(ij,l),0.)
     END DO
     rmsvl(l) = 2.*(ssum(ip1jmp1,ge,1)-ssum(jjp1,ge,iip1))

     DO ij = iip2, ip1jm
        ge(ij) = (ucov(ij,l)/cu(ij)+omegcosp(ij))*masse(ij,l)*cosphi(ij)
     END DO
     angl(l) = radsg*(ssum(ip1jm-iip1,ge(iip2),1)-ssum(jjm-1,ge(iip2),iip1) &
          )
  END DO

  DO ij = 1, ip1jmp1
     ge(ij) = ps(ij)*aire(ij)
  END DO
  ptot = ssum(ip1jmp1,ge,1) - ssum(jjp1,ge,iip1)
  etot = ssum(llm,etotl,1)
  ztot = ssum(llm,ztotl,1)
  stot = ssum(llm,stotl,1)
  rmsv = ssum(llm,rmsvl,1)
  ang = ssum(llm,angl,1)

  rday = float(int(day_ini+time))

  IF (ptot0==0.) THEN
     PRINT 3500, itau, rday, heure, time
     PRINT *, 'WARNING!!! On recalcule les valeurs initiales de :'
     PRINT *, 'ptot,rmsdpdt,etot,ztot,stot,rmsv,ang'
     PRINT *, ptot, rmsdpdt, etot, ztot, stot, rmsv, ang
     etot0 = etot
     ptot0 = ptot
     ztot0 = ztot
     stot0 = stot
     ang0 = ang
  END IF

  etot = etot/etot0
  rmsv = sqrt(rmsv/ptot)
  ptot = ptot/ptot0
  ztot = ztot/ztot0
  stot = stot/stot0
  ang = ang/ang0


  PRINT 3500, itau, rday, heure, time
  PRINT 4000, ptot, rmsdpdt, etot, ztot, stot, rmsv, ang

  RETURN

3500 FORMAT (4X,'pas',I7,5X,'jour',F5.0,'heure',F5.1,4X,'date',F10.5)
4000 FORMAT (10X,'masse',4X,'rmsdpdt',7X,'energie',2X,'enstrophie',2X, &
       'entropie',3X,'rmsv',4X,'mt.ang',/,'GLOB  ',F10.6,E13.6,5F10.3/)

END SUBROUTINE sortvarc
1	SUBROUTINE sortvarc(itau,ucov,teta,ps,masse,pk,phis,vorpot,phi,bern,dp, &
2	time)
3
4
5	! From dyn3d/sortvarc.F,v 1.1.1.1 2004/05/19 12:53:07
6
7	! Auteur: P. Le Van
8	! Objet:
9	! sortie des variables de controle
10
11	USE dimens_m, ONLY : iim, jjm, llm
12	USE paramet_m, ONLY : iip1, iip2, ijp1llm, ip1jm, ip1jmp1, jjp1
13	USE comconst, ONLY : daysec, dtvr, g, omeg, rad
14	USE comgeom, ONLY : aire, cu, rlatu
15	USE temps, ONLY : day_ini
16	USE ener, ONLY : ang, ang0, etot, etot0, ptot, ptot0, rmsdpdt, rmsv, &
17	stot, stot0, ztot, ztot0
18
19	IMPLICIT NONE
20
21	! Arguments:
22
23	INTEGER, INTENT (IN) :: itau
24	REAL :: ucov(ip1jmp1,llm), teta(ip1jmp1,llm), masse(ip1jmp1,llm)
25	REAL :: ps(ip1jmp1), phis(ip1jmp1)
26	REAL :: vorpot(ip1jm,llm)
27	REAL :: phi(ip1jmp1,llm), bern(ip1jmp1,llm)
28	REAL :: dp(ip1jmp1)
29	REAL :: time
30	REAL, INTENT (IN) :: pk(ip1jmp1,llm)
31
32	! Local:
33
34	REAL :: vor(ip1jm), bernf(ip1jmp1,llm), ztotl(llm)
35	REAL :: etotl(llm), stotl(llm), rmsvl(llm), angl(llm), ge(ip1jmp1)
36	REAL :: cosphi(ip1jm), omegcosp(ip1jm)
37	REAL :: dtvrs1j, rjour, heure, radsg, radomeg
38	REAL :: rday, massebxy(ip1jm,llm)
39	INTEGER :: l, ij, imjmp1
40
41	REAL :: ssum
42
43	!-----------------------------------------------------------------------
44
45	dtvrs1j = dtvr/daysec
46	rjour = float(int(itau*dtvrs1j))
47	heure = (itaudtvrs1j-rjour)24.
48	imjmp1 = iim*jjp1
49	IF (abs(heure-24.)<=0.0001) heure = 0.
50
51	CALL massbarxy(masse,massebxy)
52
53	! ..... Calcul de rmsdpdt .....
54
55	ge(:) = dp(:)*dp(:)
56
57	rmsdpdt = ssum(ip1jmp1,ge,1) - ssum(jjp1,ge,iip1)
58
59	rmsdpdt = daysec1.E-2sqrt(rmsdpdt/imjmp1)
60
61	CALL scopy(ijp1llm,bern,1,bernf,1)
62	CALL filtreg(bernf,jjp1,llm,-2,2,.TRUE.,1)
63
64	! ..... Calcul du moment angulaire .....
65
66	radsg = rad/g
67	radomeg = rad*omeg
68
69	DO ij = iip2, ip1jm
70	cosphi(ij) = cos(rlatu((ij-1)/iip1+1))
71	omegcosp(ij) = radomeg*cosphi(ij)
72	END DO
73
74	! ... Calcul de l'energie,de l'enstrophie,de l'entropie et de rmsv .
75
76	DO l = 1, llm
77	DO ij = 1, ip1jm
78	vor(ij) = vorpot(ij,l)vorpot(ij,l)massebxy(ij,l)
79	END DO
80	ztotl(l) = (ssum(ip1jm,vor,1)-ssum(jjm,vor,iip1))
81
82	DO ij = 1, ip1jmp1
83	ge(ij) = masse(ij,l)(phis(ij)+teta(ij,l)pk(ij,l)+bernf(ij,l)-phi( &
84	ij,l))
85	END DO
86	etotl(l) = ssum(ip1jmp1,ge,1) - ssum(jjp1,ge,iip1)
87
88	DO ij = 1, ip1jmp1
89	ge(ij) = masse(ij,l)*teta(ij,l)
90	END DO
91	stotl(l) = ssum(ip1jmp1,ge,1) - ssum(jjp1,ge,iip1)
92
93	DO ij = 1, ip1jmp1
94	ge(ij) = masse(ij,l)*amax1(bernf(ij,l)-phi(ij,l),0.)
95	END DO
96	rmsvl(l) = 2.*(ssum(ip1jmp1,ge,1)-ssum(jjp1,ge,iip1))
97
98	DO ij = iip2, ip1jm
99	ge(ij) = (ucov(ij,l)/cu(ij)+omegcosp(ij))masse(ij,l)cosphi(ij)
100	END DO
101	angl(l) = radsg*(ssum(ip1jm-iip1,ge(iip2),1)-ssum(jjm-1,ge(iip2),iip1) &
102	)
103	END DO
104
105	DO ij = 1, ip1jmp1
106	ge(ij) = ps(ij)*aire(ij)
107	END DO
108	ptot = ssum(ip1jmp1,ge,1) - ssum(jjp1,ge,iip1)
109	etot = ssum(llm,etotl,1)
110	ztot = ssum(llm,ztotl,1)
111	stot = ssum(llm,stotl,1)
112	rmsv = ssum(llm,rmsvl,1)
113	ang = ssum(llm,angl,1)
114
115	rday = float(int(day_ini+time))
116
117	IF (ptot0==0.) THEN
118	PRINT 3500, itau, rday, heure, time
119	PRINT *, 'WARNING!!! On recalcule les valeurs initiales de :'
120	PRINT *, 'ptot,rmsdpdt,etot,ztot,stot,rmsv,ang'
121	PRINT *, ptot, rmsdpdt, etot, ztot, stot, rmsv, ang
122	etot0 = etot
123	ptot0 = ptot
124	ztot0 = ztot
125	stot0 = stot
126	ang0 = ang
127	END IF
128
129	etot = etot/etot0
130	rmsv = sqrt(rmsv/ptot)
131	ptot = ptot/ptot0
132	ztot = ztot/ztot0
133	stot = stot/stot0
134	ang = ang/ang0
135
136
137	PRINT 3500, itau, rday, heure, time
138	PRINT 4000, ptot, rmsdpdt, etot, ztot, stot, rmsv, ang
139
140	RETURN
141
142	3500 FORMAT (4X,'pas',I7,5X,'jour',F5.0,'heure',F5.1,4X,'date',F10.5)
143	4000 FORMAT (10X,'masse',4X,'rmsdpdt',7X,'energie',2X,'enstrophie',2X, &
144	'entropie',3X,'rmsv',4X,'mt.ang',/,'GLOB ',F10.6,E13.6,5F10.3/)
145
146	END SUBROUTINE sortvarc