source: trunk/LMDZ.GENERIC/libf/dyn3d/sortvarc.F @ 773

Last change on this file since 773 was 135, checked in by aslmd, 14 years ago

CHANGEMENT ARBORESCENCE ETAPE 2 -- NON COMPLET

File size: 4.4 KB
Line 
1      SUBROUTINE sortvarc
2     $(itau,ucov,teta,ps,masse,pk,phis,vorpot,phi,bern,dp,time ,
3     $ vcov )
4      IMPLICIT NONE
5
6c=======================================================================
7c
8c   Auteur:    P. Le Van
9c   -------
10c
11c   Objet:
12c   ------
13c
14c   sortie des variables de controle
15c
16c=======================================================================
17c-----------------------------------------------------------------------
18c   Declarations:
19c   -------------
20
21#include "dimensions.h"
22#include "paramet.h"
23#include "comconst.h"
24#include "comvert.h"
25#include "comgeom.h"
26#include "ener.h"
27#include "logic.h"
28#include "temps.h"
29
30c   Arguments:
31c   ----------
32
33      INTEGER itau
34      REAL ucov(ip1jmp1,llm),teta(ip1jmp1,llm),masse(ip1jmp1,llm)
35      REAL vcov(ip1jm,llm)
36      REAL ps(ip1jmp1),phis(ip1jmp1)
37      REAL vorpot(ip1jm,llm)
38      REAL phi(ip1jmp1,llm),bern(ip1jmp1,llm)
39      REAL dp(ip1jmp1)
40      REAL time
41      REAL pk(ip1jmp1,llm)
42
43c   Local:
44c   ------
45
46      REAL vor(ip1jm),bernf(ip1jmp1,llm),ztotl(llm)
47      REAL etotl(llm),stotl(llm),rmsvl(llm),angl(llm),ge(ip1jmp1)
48      REAL cosphi(ip1jm),omegcosp(ip1jm)
49      REAL dtvrs1j,rjour,heure,radsg,radomeg
50      REAL rday, massebxy(ip1jm,llm)
51      INTEGER  l, ij, imjmp1
52
53      EXTERNAL  filtreg, massbarxy
54c     EXTERNAL FLUSH
55      EXTERNAL   SSUM, SCOPY
56      REAL       SSUM
57
58c-----------------------------------------------------------------------
59
60       dtvrs1j   = dtvr/daysec
61       rjour     = FLOAT( INT( itau * dtvrs1j ))
62       heure     = ( itau*dtvrs1j-rjour ) * 24.
63       imjmp1    = iim * jjp1
64       IF(ABS(heure - 24.).LE.0.0001 ) heure = 0.
65c
66       CALL massbarxy ( masse, massebxy )
67
68c   .....  Calcul  de  rmsdpdt  .....
69
70       CALL multipl(ip1jmp1,dp,dp,ge)
71
72       rmsdpdt = SSUM(ip1jmp1,ge,1) - SSUM(jjp1,ge,iip1)
73c
74       rmsdpdt = daysec* 1.e-2 * SQRT(rmsdpdt/imjmp1)
75
76       CALL SCOPY( ijp1llm,bern,1,bernf,1 )
77       CALL filtreg(bernf,jjp1,llm,-2,2,.TRUE.,1)
78
79c   .....  Calcul du moment  angulaire   .....
80
81       radsg    = rad /g
82       radomeg  = rad * omeg
83c
84       DO ij=iip2,ip1jm
85          cosphi( ij ) = COS(rlatu((ij-1)/iip1+1))
86          omegcosp(ij) = radomeg   * cosphi(ij)
87       ENDDO
88
89c  ...  Calcul  de l'energie,de l'enstrophie,de l'entropie et de rmsv  .
90
91       DO l=1,llm
92          DO ij = 1,ip1jm
93             vor(ij)=vorpot(ij,l)*vorpot(ij,l)*massebxy(ij,l)
94          ENDDO
95          ztotl(l)=(SSUM(ip1jm,vor,1)-SSUM(jjm,vor,iip1))
96
97          DO ij = 1,ip1jmp1
98             ge(ij)= masse(ij,l)*(phis(ij)+teta(ij,l)*pk(ij,l)  +
99     s        bernf(ij,l)-phi(ij,l))
100          ENDDO
101          etotl(l) = SSUM(ip1jmp1,ge,1) - SSUM(jjp1,ge,iip1)
102
103          DO   ij   = 1, ip1jmp1
104             ge(ij) = masse(ij,l)*teta(ij,l)
105          ENDDO
106          stotl(l)= SSUM(ip1jmp1,ge,1) - SSUM(jjp1,ge,iip1)
107
108          DO ij=1,ip1jmp1
109             ge(ij)=masse(ij,l)*AMAX1(bernf(ij,l)-phi(ij,l),0.)
110          ENDDO
111          rmsvl(l)=2.*(SSUM(ip1jmp1,ge,1)-SSUM(jjp1,ge,iip1))
112
113          DO ij =iip2,ip1jm
114             ge(ij)=(ucov(ij,l)/cu(ij)+omegcosp(ij))*masse(ij,l) *
115     *               cosphi(ij)
116          ENDDO
117          angl(l) = radsg *
118     s    (SSUM(ip1jm-iip1,ge(iip2),1)-SSUM(jjm-1,ge(iip2),iip1))
119      ENDDO
120
121          DO ij=1,ip1jmp1
122            ge(ij)= ps(ij)*aire(ij)
123          ENDDO
124      ptot  = SSUM(ip1jmp1,ge,1)-SSUM(jjp1,ge,iip1)
125      etot  = SSUM(     llm, etotl, 1 )
126      ztot  = SSUM(     llm, ztotl, 1 )
127      stot  = SSUM(     llm, stotl, 1 )
128      rmsv  = SSUM(     llm, rmsvl, 1 )
129      ang   = SSUM(     llm,  angl, 1 )
130
131      rday = FLOAT(INT ( day_ini + time ))
132c
133      IF(ptot0.eq.0.)  THEN
134         PRINT 3500, itau, rday, heure,time
135         PRINT*,'WARNING!!! On recalcule les valeurs initiales de :'
136         PRINT*,'ptot,rmsdpdt,etot,ztot,stot,rmsv,ang'
137         PRINT *, ptot,rmsdpdt,etot,ztot,stot,rmsv,ang
138         etot0 = etot
139         ptot0 = ptot
140         ztot0 = ztot
141         stot0 = stot
142         ang0  = ang
143      END IF
144
145      etot= etot/etot0
146      rmsv= SQRT(rmsv/ptot)
147      ptot= ptot/ptot0
148      ztot= ztot/ztot0
149      stot= stot/stot0
150      ang = ang /ang0
151
152
153      PRINT 3500, itau, rday, heure, time
154      PRINT 4000, ptot,rmsdpdt,etot,ztot,stot,rmsv,ang
155
156      RETURN
157
1583500   FORMAT('0',10(1h*),4x,'pas',i7,5x,'jour',f5.0,'heure',f5.1,4x
159     *   ,'date',f10.5,4x,10(1h*))
1604000   FORMAT(10x,'masse',4x,'rmsdpdt',7x,'energie',2x,'enstrophie'
161     * ,2x,'entropie',3x,'rmsv',4x,'mt.ang'/6x,f10.6,e13.6,5f10.3/
162     * )
163      END
164
Note: See TracBrowser for help on using the repository browser.