source: LMDZ4/branches/V3_test/libf/dyn3d/sortvarc.F @ 4360

Last change on this file since 4360 was 524, checked in by lmdzadmin, 20 years ago

Initial revision

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 4.4 KB
Line 
1!
2! $Header$
3!
4      SUBROUTINE sortvarc
5     $(itau,ucov,teta,ps,masse,pk,phis,vorpot,phi,bern,dp,time ,
6     $ vcov )
7      IMPLICIT NONE
8
9c=======================================================================
10c
11c   Auteur:    P. Le Van
12c   -------
13c
14c   Objet:
15c   ------
16c
17c   sortie des variables de controle
18c
19c=======================================================================
20c-----------------------------------------------------------------------
21c   Declarations:
22c   -------------
23
24#include "dimensions.h"
25#include "paramet.h"
26#include "comconst.h"
27#include "comvert.h"
28#include "comgeom.h"
29#include "ener.h"
30#include "logic.h"
31#include "temps.h"
32
33c   Arguments:
34c   ----------
35
36      INTEGER itau
37      REAL ucov(ip1jmp1,llm),teta(ip1jmp1,llm),masse(ip1jmp1,llm)
38      REAL vcov(ip1jm,llm)
39      REAL ps(ip1jmp1),phis(ip1jmp1)
40      REAL vorpot(ip1jm,llm)
41      REAL phi(ip1jmp1,llm),bern(ip1jmp1,llm)
42      REAL dp(ip1jmp1)
43      REAL time
44      REAL pk(ip1jmp1,llm)
45
46c   Local:
47c   ------
48
49      REAL vor(ip1jm),bernf(ip1jmp1,llm),ztotl(llm)
50      REAL etotl(llm),stotl(llm),rmsvl(llm),angl(llm),ge(ip1jmp1)
51      REAL cosphi(ip1jm),omegcosp(ip1jm)
52      REAL dtvrs1j,rjour,heure,radsg,radomeg
53      REAL rday, massebxy(ip1jm,llm)
54      INTEGER  l, ij, imjmp1
55
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       ge(:)=dp(:)*dp(:)
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',/,'GLOB  '
162     .  ,f10.6,e13.6,5f10.3/
163     * )
164      END
165
Note: See TracBrowser for help on using the repository browser.