/[lmdze]/trunk/dyn3d/sortvarc.f
ViewVC logotype

Diff of /trunk/dyn3d/sortvarc.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/libf/dyn3d/sortvarc.f revision 3 by guez, Wed Feb 27 13:16:39 2008 UTC trunk/libf/dyn3d/sortvarc.f90 revision 23 by guez, Mon Dec 14 15:25:16 2009 UTC
# Line 1  Line 1 
1  !  SUBROUTINE sortvarc(itau,ucov,teta,ps,masse,pk,phis,vorpot,phi,bern,dp, &
2  ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/sortvarc.F,v 1.1.1.1 2004/05/19 12:53:07 lmdzadmin Exp $       time)
 !  
       SUBROUTINE sortvarc  
      $(itau,ucov,teta,ps,masse,pk,phis,vorpot,phi,bern,dp,time ,  
      $ vcov )  
   
 c=======================================================================  
 c  
 c   Auteur:    P. Le Van  
 c   -------  
 c  
 c   Objet:  
 c   ------  
 c  
 c   sortie des variables de controle  
 c  
 c=======================================================================  
 c-----------------------------------------------------------------------  
 c   Declarations:  
 c   -------------  
   
       use dimens_m  
       use paramet_m  
       use comconst  
       use comvert  
       use logic  
       use comgeom  
       use temps  
       use ener  
       IMPLICIT NONE  
   
 c   Arguments:  
 c   ----------  
   
       INTEGER, intent(in):: itau  
       REAL ucov(ip1jmp1,llm),teta(ip1jmp1,llm),masse(ip1jmp1,llm)  
       REAL vcov(ip1jm,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)  
   
 c   Local:  
 c   ------  
   
       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  
   
 c-----------------------------------------------------------------------  
   
        dtvrs1j   = dtvr/daysec  
        rjour     = FLOAT( INT( itau * dtvrs1j ))  
        heure     = ( itau*dtvrs1j-rjour ) * 24.  
        imjmp1    = iim * jjp1  
        IF(ABS(heure - 24.).LE.0.0001 ) heure = 0.  
 c  
        CALL massbarxy ( masse, massebxy )  
   
 c   .....  Calcul  de  rmsdpdt  .....  
   
        ge(:)=dp(:)*dp(:)  
   
        rmsdpdt = SSUM(ip1jmp1,ge,1) - SSUM(jjp1,ge,iip1)  
 c  
        rmsdpdt = daysec* 1.e-2 * SQRT(rmsdpdt/imjmp1)  
   
        CALL SCOPY( ijp1llm,bern,1,bernf,1 )  
        CALL filtreg(bernf,jjp1,llm,-2,2,.TRUE.,1)  
   
 c   .....  Calcul du moment  angulaire   .....  
   
        radsg    = rad /g  
        radomeg  = rad * omeg  
 c  
        DO ij=iip2,ip1jm  
           cosphi( ij ) = COS(rlatu((ij-1)/iip1+1))  
           omegcosp(ij) = radomeg   * cosphi(ij)  
        ENDDO  
   
 c  ...  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)  
           ENDDO  
           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)  +  
      s        bernf(ij,l)-phi(ij,l))  
           ENDDO  
           etotl(l) = SSUM(ip1jmp1,ge,1) - SSUM(jjp1,ge,iip1)  
   
           DO   ij   = 1, ip1jmp1  
              ge(ij) = masse(ij,l)*teta(ij,l)  
           ENDDO  
           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.)  
           ENDDO  
           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)  
           ENDDO  
           angl(l) = radsg *  
      s    (SSUM(ip1jm-iip1,ge(iip2),1)-SSUM(jjm-1,ge(iip2),iip1))  
       ENDDO  
   
           DO ij=1,ip1jmp1  
             ge(ij)= ps(ij)*aire(ij)  
           ENDDO  
       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 ))  
 c  
       IF(ptot0.eq.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  
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 = (itau*dtvrs1j-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 = daysec*1.E-2*sqrt(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

Legend:
Removed from v.3  
changed lines
  Added in v.23

  ViewVC Help
Powered by ViewVC 1.1.21