source: LMDZ5/branches/testing/libf/dyn3d_common/sortvarc.F

Last change on this file was 2641, checked in by Laurent Fairhead, 8 years ago

Merged trunk changes r2593:2640 into testing branch

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 5.6 KB
RevLine 
[524]1!
[1403]2! $Id: sortvarc.F 2641 2016-09-29 21:26:46Z fairhead $
[524]3!
4      SUBROUTINE sortvarc
5     $(itau,ucov,teta,ps,masse,pk,phis,vorpot,phi,bern,dp,time ,
6     $ vcov )
[2160]7
8      USE control_mod, ONLY: resetvarc
[2641]9      USE comconst_mod, ONLY: dtvr, daysec, g, rad, omeg
10      USE logic_mod, ONLY: read_start
11      USE ener_mod, ONLY: etot,ptot,ztot,stot,ang,
12     &                    etot0,ptot0,ztot0,stot0,ang0,
13     &                    rmsdpdt,rmsv
[524]14      IMPLICIT NONE
15
[2160]16
[524]17c=======================================================================
18c
19c   Auteur:    P. Le Van
20c   -------
21c
22c   Objet:
23c   ------
24c
25c   sortie des variables de controle
26c
27c=======================================================================
28c-----------------------------------------------------------------------
29c   Declarations:
30c   -------------
31
[2160]32      INCLUDE "dimensions.h"
33      INCLUDE "paramet.h"
34      INCLUDE "comgeom.h"
35      INCLUDE "iniprint.h"
[524]36
37c   Arguments:
38c   ----------
39
[2160]40      INTEGER,INTENT(IN) :: itau
41      REAL,INTENT(IN) :: ucov(ip1jmp1,llm)
42      REAL,INTENT(IN) :: teta(ip1jmp1,llm)
43      REAL,INTENT(IN) :: masse(ip1jmp1,llm)
44      REAL,INTENT(IN) :: vcov(ip1jm,llm)
45      REAL,INTENT(IN) :: ps(ip1jmp1)
46      REAL,INTENT(IN) :: phis(ip1jmp1)
47      REAL,INTENT(IN) :: vorpot(ip1jm,llm)
48      REAL,INTENT(IN) :: phi(ip1jmp1,llm)
49      REAL,INTENT(IN) :: bern(ip1jmp1,llm)
50      REAL,INTENT(IN) :: dp(ip1jmp1)
51      REAL,INTENT(IN) :: time
52      REAL,INTENT(IN) :: pk(ip1jmp1,llm)
[524]53
54c   Local:
55c   ------
56
57      REAL vor(ip1jm),bernf(ip1jmp1,llm),ztotl(llm)
58      REAL etotl(llm),stotl(llm),rmsvl(llm),angl(llm),ge(ip1jmp1)
59      REAL cosphi(ip1jm),omegcosp(ip1jm)
60      REAL dtvrs1j,rjour,heure,radsg,radomeg
[2160]61      REAL massebxy(ip1jm,llm)
[524]62      INTEGER  l, ij, imjmp1
63
64      REAL       SSUM
[2160]65      LOGICAL,SAVE :: firstcal=.true.
66      CHARACTER(LEN=*),PARAMETER :: modname="sortvarc"
[524]67
68c-----------------------------------------------------------------------
[2160]69! Ehouarn: when no initialization fields from file, resetvarc should be
70!          set to false
71       if (firstcal) then
72         if (.not.read_start) then
73           resetvarc=.true.
74         endif
75       endif
[524]76
77       dtvrs1j   = dtvr/daysec
[1403]78       rjour     = REAL( INT( itau * dtvrs1j ))
[524]79       heure     = ( itau*dtvrs1j-rjour ) * 24.
80       imjmp1    = iim * jjp1
81       IF(ABS(heure - 24.).LE.0.0001 ) heure = 0.
82c
83       CALL massbarxy ( masse, massebxy )
84
85c   .....  Calcul  de  rmsdpdt  .....
86
87       ge(:)=dp(:)*dp(:)
88
89       rmsdpdt = SSUM(ip1jmp1,ge,1) - SSUM(jjp1,ge,iip1)
90c
91       rmsdpdt = daysec* 1.e-2 * SQRT(rmsdpdt/imjmp1)
92
93       CALL SCOPY( ijp1llm,bern,1,bernf,1 )
94       CALL filtreg(bernf,jjp1,llm,-2,2,.TRUE.,1)
95
96c   .....  Calcul du moment  angulaire   .....
97
98       radsg    = rad /g
99       radomeg  = rad * omeg
100c
101       DO ij=iip2,ip1jm
102          cosphi( ij ) = COS(rlatu((ij-1)/iip1+1))
103          omegcosp(ij) = radomeg   * cosphi(ij)
104       ENDDO
105
106c  ...  Calcul  de l'energie,de l'enstrophie,de l'entropie et de rmsv  .
107
108       DO l=1,llm
109          DO ij = 1,ip1jm
110             vor(ij)=vorpot(ij,l)*vorpot(ij,l)*massebxy(ij,l)
111          ENDDO
112          ztotl(l)=(SSUM(ip1jm,vor,1)-SSUM(jjm,vor,iip1))
113
114          DO ij = 1,ip1jmp1
115             ge(ij)= masse(ij,l)*(phis(ij)+teta(ij,l)*pk(ij,l)  +
116     s        bernf(ij,l)-phi(ij,l))
117          ENDDO
118          etotl(l) = SSUM(ip1jmp1,ge,1) - SSUM(jjp1,ge,iip1)
119
120          DO   ij   = 1, ip1jmp1
121             ge(ij) = masse(ij,l)*teta(ij,l)
122          ENDDO
123          stotl(l)= SSUM(ip1jmp1,ge,1) - SSUM(jjp1,ge,iip1)
124
125          DO ij=1,ip1jmp1
126             ge(ij)=masse(ij,l)*AMAX1(bernf(ij,l)-phi(ij,l),0.)
127          ENDDO
128          rmsvl(l)=2.*(SSUM(ip1jmp1,ge,1)-SSUM(jjp1,ge,iip1))
129
130          DO ij =iip2,ip1jm
131             ge(ij)=(ucov(ij,l)/cu(ij)+omegcosp(ij))*masse(ij,l) *
132     *               cosphi(ij)
133          ENDDO
[2160]134          angl(l) = rad *
[524]135     s    (SSUM(ip1jm-iip1,ge(iip2),1)-SSUM(jjm-1,ge(iip2),iip1))
136      ENDDO
137
138          DO ij=1,ip1jmp1
139            ge(ij)= ps(ij)*aire(ij)
140          ENDDO
141      ptot  = SSUM(ip1jmp1,ge,1)-SSUM(jjp1,ge,iip1)
142      etot  = SSUM(     llm, etotl, 1 )
143      ztot  = SSUM(     llm, ztotl, 1 )
144      stot  = SSUM(     llm, stotl, 1 )
145      rmsv  = SSUM(     llm, rmsvl, 1 )
146      ang   = SSUM(     llm,  angl, 1 )
147
[2160]148      IF (firstcal.and.resetvarc) then
149         WRITE(lunout,3500) itau, rjour, heure, time
150         WRITE(lunout,*) trim(modname),
151     &     ' WARNING!!! Recomputing initial values of : '
152         WRITE(lunout,*) 'ptot,rmsdpdt,etot,ztot,stot,rmsv,ang'
153         WRITE(lunout,*) ptot,rmsdpdt,etot,ztot,stot,rmsv,ang
[524]154         etot0 = etot
155         ptot0 = ptot
156         ztot0 = ztot
157         stot0 = stot
158         ang0  = ang
159      END IF
160
[2160]161      ! compute relative changes in etot,... (except if 'reference' values
162      ! are zero, which can happen when using iniacademic)
163      if (etot0.ne.0) then
164        etot= etot/etot0
165      else
166        etot=1.
167      endif
[524]168      rmsv= SQRT(rmsv/ptot)
[2160]169      if (ptot0.ne.0) then
170        ptot= ptot/ptot0
171      else
172        ptot=1.
173      endif
174      if (ztot0.ne.0) then
175        ztot= ztot/ztot0
176      else
177        ztot=1.
178      endif
179      if (stot0.ne.0) then
180        stot= stot/stot0
181      else
182        stot=1.
183      endif
184      if (ang0.ne.0) then
185        ang = ang /ang0
186      else
187        ang=1.
188      endif
[524]189
[2160]190      firstcal = .false.
[524]191
[2160]192      WRITE(lunout,3500) itau, rjour, heure, time
193      WRITE(lunout,4000) ptot,rmsdpdt,etot,ztot,stot,rmsv,ang
[524]194
[1279]1953500   FORMAT(10("*"),4x,'pas',i7,5x,'jour',f9.0,'heure',f5.1,4x
196     *   ,'date',f14.4,4x,10("*"))
[524]1974000   FORMAT(10x,'masse',4x,'rmsdpdt',7x,'energie',2x,'enstrophie'
198     * ,2x,'entropie',3x,'rmsv',4x,'mt.ang',/,'GLOB  '
199     .  ,f10.6,e13.6,5f10.3/
200     * )
201      END
202
Note: See TracBrowser for help on using the repository browser.