source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/config/ppsrc/dyn3d_common/sortvarc.f @ 224

Last change on this file since 224 was 224, checked in by ymipsl, 10 years ago
File size: 15.2 KB
Line 
1!
2! $Id: sortvarc.F 1403 2010-07-01 09:02:53Z fairhead $
3!
4      SUBROUTINE sortvarc
5     $(itau,ucov,teta,ps,masse,pk,phis,vorpot,phi,bern,dp,time ,
6     $ vcov )
7
8      use control_mod,only:resetvarc 
9      IMPLICIT NONE
10
11c=======================================================================
12c
13c   Auteur:    P. Le Van
14c   -------
15c
16c   Objet:
17c   ------
18c
19c   sortie des variables de controle
20c
21c=======================================================================
22c-----------------------------------------------------------------------
23c   Declarations:
24c   -------------
25
26!-----------------------------------------------------------------------
27!   INCLUDE 'dimensions.h'
28!
29!   dimensions.h contient les dimensions du modele
30!   ndm est tel que iim=2**ndm
31!-----------------------------------------------------------------------
32
33      INTEGER iim,jjm,llm,ndm
34
35      PARAMETER (iim= 128,jjm=96,llm=64,ndm=1)
36
37!-----------------------------------------------------------------------
38!
39! $Header$
40!
41!
42!  ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
43!                 veillez  n'utiliser que des ! pour les commentaires
44!                 et  bien positionner les & des lignes de continuation
45!                 (les placer en colonne 6 et en colonne 73)
46!
47!
48!-----------------------------------------------------------------------
49!   INCLUDE 'paramet.h'
50
51      INTEGER  iip1,iip2,iip3,jjp1,llmp1,llmp2,llmm1
52      INTEGER  kftd,ip1jm,ip1jmp1,ip1jmi1,ijp1llm
53      INTEGER  ijmllm,mvar
54      INTEGER jcfil,jcfllm
55
56      PARAMETER( iip1= iim+1,iip2=iim+2,iip3=iim+3                       &
57     &    ,jjp1=jjm+1-1/jjm)
58      PARAMETER( llmp1 = llm+1,  llmp2 = llm+2, llmm1 = llm-1 )
59      PARAMETER( kftd  = iim/2 -ndm )
60      PARAMETER( ip1jm  = iip1*jjm,  ip1jmp1= iip1*jjp1 )
61      PARAMETER( ip1jmi1= ip1jm - iip1 )
62      PARAMETER( ijp1llm= ip1jmp1 * llm, ijmllm= ip1jm * llm )
63      PARAMETER( mvar= ip1jmp1*( 2*llm+1) + ijmllm )
64      PARAMETER( jcfil=jjm/2+5, jcfllm=jcfil*llm )
65
66!-----------------------------------------------------------------------
67!
68! $Id: comconst.h 1437 2010-09-30 08:29:10Z emillour $
69!
70!-----------------------------------------------------------------------
71! INCLUDE comconst.h
72
73      COMMON/comconsti/im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl,          &
74     &                 iflag_top_bound,mode_top_bound
75      COMMON/comconstr/dtvr,daysec,                                     &
76     & pi,dtphys,dtdiss,rad,r,kappa,cotot,unsim,g,omeg                  &
77     & ,dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta        &
78     & ,dissip_pupstart  ,tau_top_bound,                                &
79     & daylen,molmass, ihf
80      COMMON/cpdetvenus/cpp,nu_venus,t0_venus
81
82      INTEGER im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl
83      REAL dtvr ! dynamical time step (in s)
84      REAL daysec !length (in s) of a standard day
85      REAL pi    ! something like 3.14159....
86      REAL dtphys ! (s) time step for the physics
87      REAL dtdiss ! (s) time step for the dissipation
88      REAL rad ! (m) radius of the planet
89      REAL r ! Reduced Gas constant r=R/mu
90             ! with R=8.31.. J.K-1.mol-1, mu: mol mass of atmosphere (kg/mol)
91      REAL cpp   ! Cp
92      REAL kappa ! kappa=R/Cp
93      REAL cotot
94      REAL unsim ! = 1./iim
95      REAL g ! (m/s2) gravity
96      REAL omeg ! (rad/s) rotation rate of the planet
97! Dissipation factors, for Earth model:
98      REAL dissip_factz,dissip_zref !dissip_deltaz
99! Dissipation factors, for other planets:
100      REAL dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta
101      REAL dissip_pupstart
102      INTEGER iflag_top_bound,mode_top_bound
103      REAL tau_top_bound
104      REAL daylen ! length of solar day, in 'standard' day length
105      REAL molmass ! (g/mol) molar mass of the atmosphere
106
107      REAL nu_venus,t0_venus ! coeffs needed for Cp(T), Venus atmosphere
108      REAL ihf  ! (W/m2) intrinsic heat flux for giant planets
109
110
111!-----------------------------------------------------------------------
112!
113! $Id: comvert.h 1654 2012-09-24 15:07:18Z aslmd $
114!
115!-----------------------------------------------------------------------
116!   INCLUDE 'comvert.h'
117
118      COMMON/comvertr/ap(llm+1),bp(llm+1),presnivs(llm),dpres(llm),     &
119     &               pa,preff,nivsigs(llm),nivsig(llm+1),               &
120     &               aps(llm),bps(llm),scaleheight,pseudoalt(llm)
121
122      common/comverti/disvert_type, pressure_exner
123
124      real ap     ! hybrid pressure contribution at interlayers
125      real bp     ! hybrid sigma contribution at interlayer
126      real presnivs ! (reference) pressure at mid-layers
127      real dpres
128      real pa     ! reference pressure (Pa) at which hybrid coordinates
129                  ! become purely pressure
130      real preff  ! reference surface pressure (Pa)
131      real nivsigs
132      real nivsig
133      real aps    ! hybrid pressure contribution at mid-layers
134      real bps    ! hybrid sigma contribution at mid-layers
135      real scaleheight ! atmospheric (reference) scale height (km)
136      real pseudoalt ! pseudo-altitude of model levels (km), based on presnivs(),
137                     ! preff and scaleheight
138
139      integer disvert_type ! type of vertical discretization:
140                           ! 1: Earth (default for planet_type==earth),
141                           !     automatic generation
142                           ! 2: Planets (default for planet_type!=earth),
143                           !     using 'z2sig.def' (or 'esasig.def) file
144
145      logical pressure_exner
146!     compute pressure inside layers using Exner function, else use mean
147!     of pressure values at interfaces
148
149 !-----------------------------------------------------------------------
150!
151! $Header$
152!
153!CDK comgeom
154      COMMON/comgeom/                                                   &
155     & cu(ip1jmp1),cv(ip1jm),unscu2(ip1jmp1),unscv2(ip1jm),             &
156     & aire(ip1jmp1),airesurg(ip1jmp1),aireu(ip1jmp1),                  &
157     & airev(ip1jm),unsaire(ip1jmp1),apoln,apols,                       &
158     & unsairez(ip1jm),airuscv2(ip1jm),airvscu2(ip1jm),                 &
159     & aireij1(ip1jmp1),aireij2(ip1jmp1),aireij3(ip1jmp1),              &
160     & aireij4(ip1jmp1),alpha1(ip1jmp1),alpha2(ip1jmp1),                &
161     & alpha3(ip1jmp1),alpha4(ip1jmp1),alpha1p2(ip1jmp1),               &
162     & alpha1p4(ip1jmp1),alpha2p3(ip1jmp1),alpha3p4(ip1jmp1),           &
163     & fext(ip1jm),constang(ip1jmp1),rlatu(jjp1),rlatv(jjm),            &
164     & rlonu(iip1),rlonv(iip1),cuvsurcv(ip1jm),cvsurcuv(ip1jm),         &
165     & cvusurcu(ip1jmp1),cusurcvu(ip1jmp1),cuvscvgam1(ip1jm),           &
166     & cuvscvgam2(ip1jm),cvuscugam1(ip1jmp1),                           &
167     & cvuscugam2(ip1jmp1),cvscuvgam(ip1jm),cuscvugam(ip1jmp1),         &
168     & unsapolnga1,unsapolnga2,unsapolsga1,unsapolsga2,                 &
169     & unsair_gam1(ip1jmp1),unsair_gam2(ip1jmp1),unsairz_gam(ip1jm),    &
170     & aivscu2gam(ip1jm),aiuscv2gam(ip1jm),xprimu(iip1),xprimv(iip1)
171
172!
173        REAL                                                            &
174     & cu,cv,unscu2,unscv2,aire,airesurg,aireu,airev,unsaire,apoln     ,&
175     & apols,unsairez,airuscv2,airvscu2,aireij1,aireij2,aireij3,aireij4,&
176     & alpha1,alpha2,alpha3,alpha4,alpha1p2,alpha1p4,alpha2p3,alpha3p4 ,&
177     & fext,constang,rlatu,rlatv,rlonu,rlonv,cuvscvgam1,cuvscvgam2     ,&
178     & cvuscugam1,cvuscugam2,cvscuvgam,cuscvugam,unsapolnga1,unsapolnga2&
179     & ,unsapolsga1,unsapolsga2,unsair_gam1,unsair_gam2,unsairz_gam    ,&
180     & aivscu2gam ,aiuscv2gam,cuvsurcv,cvsurcuv,cvusurcu,cusurcvu,xprimu&
181     & , xprimv
182!
183!
184! $Id: ener.h 1447 2010-10-22 16:18:27Z jghattas $
185!
186!  ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
187!                 veillez à n'utiliser que des ! pour les commentaires
188!                 et à bien positionner les & des lignes de continuation
189!                 (les placer en colonne 6 et en colonne 73)
190!
191! INCLUDE 'ener.h'
192
193      COMMON/ener/ang0,etot0,ptot0,ztot0,stot0,                         &
194     &            ang,etot,ptot,ztot,stot,rmsdpdt ,                     &
195     &            rmsv,gtot(llmm1)
196
197      REAL ang0,etot0,ptot0,ztot0,stot0,                                &
198     &     ang,etot,ptot,ztot,stot,rmsdpdt,rmsv,gtot
199
200!-----------------------------------------------------------------------
201!
202! $Id: logic.h 1520 2011-05-23 11:37:09Z emillour $
203!
204!
205! NB: keep items of different kinds in seperate common blocs to avoid
206!     "misaligned commons" issues
207!-----------------------------------------------------------------------
208! INCLUDE 'logic.h'
209
210      COMMON/logicl/ purmats,forward,leapf,apphys,                      &
211     &  statcl,conser,apdiss,apdelq,saison,ecripar,fxyhypb,ysinus       &
212     &  ,read_start,ok_guide,ok_strato,tidal,ok_gradsfile               &
213     &  ,ok_limit,ok_etat0,hybrid                                       &
214     &  ,moyzon_mu,moyzon_ch
215
216      COMMON/logici/ iflag_phys,iflag_trac
217     
218      LOGICAL purmats,forward,leapf,apphys,statcl,conser,               &
219     & apdiss,apdelq,saison,ecripar,fxyhypb,ysinus                      &
220     &  ,read_start,ok_guide,ok_strato,tidal,ok_gradsfile               &
221     &  ,ok_limit,ok_etat0
222      logical hybrid ! vertical coordinate is hybrid if true (sigma otherwise)
223                     ! (only used if disvert_type==2)
224      logical moyzon_mu,moyzon_ch ! used for zonal averages in Titan
225
226      integer iflag_phys,iflag_trac
227!$OMP THREADPRIVATE(/logicl/)
228!$OMP THREADPRIVATE(/logici/)
229!-----------------------------------------------------------------------
230!
231! $Id: temps.h 1577 2011-10-20 15:06:47Z fairhead $
232!
233!  ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
234!                 veillez  n'utiliser que des ! pour les commentaires
235!                 et  bien positionner les & des lignes de continuation
236!                 (les placer en colonne 6 et en colonne 73)
237!
238!
239! jD_ref = jour julien de la date de reference (lancement de l'experience)
240! hD_ref = "heure" julienne de la date de reference
241!-----------------------------------------------------------------------
242! INCLUDE 'temps.h'
243
244      COMMON/temps_r/dt,jD_ref,jH_ref,start_time,hour_ini
245      COMMON/temps_i/day_ini,day_end,annee_ref,day_ref,                 &
246     &             itau_dyn,itau_phy,itaufin
247      COMMON/temps_c/calend
248
249
250      INTEGER   itaufin ! total number of dynamical steps for the run
251      INTEGER   itau_dyn, itau_phy
252      INTEGER   day_ini ! initial day # of simulation sequence
253      INTEGER   day_end ! final day # ; i.e. day # when this simulation ends
254      INTEGER   annee_ref
255      INTEGER   day_ref
256      REAL      dt ! (dynamics) time step (changes if doing Matsuno or LF step)
257      REAL      jD_ref, jH_ref, start_time
258      CHARACTER (len=10) :: calend
259
260      ! Additionnal Mars stuff:
261      real hour_ini ! initial fraction of day of simulation sequence (0=<hour_ini<1)
262
263!-----------------------------------------------------------------------
264
265c   Arguments:
266c   ----------
267
268      INTEGER itau
269      REAL ucov(ip1jmp1,llm),teta(ip1jmp1,llm),masse(ip1jmp1,llm)
270      REAL vcov(ip1jm,llm)
271      REAL ps(ip1jmp1),phis(ip1jmp1)
272      REAL vorpot(ip1jm,llm)
273      REAL phi(ip1jmp1,llm),bern(ip1jmp1,llm)
274      REAL dp(ip1jmp1)
275      REAL time
276      REAL pk(ip1jmp1,llm)
277
278c   Local:
279c   ------
280
281      REAL vor(ip1jm),bernf(ip1jmp1,llm),ztotl(llm)
282      REAL etotl(llm),stotl(llm),rmsvl(llm),angl(llm),ge(ip1jmp1)
283      REAL cosphi(ip1jm),omegcosp(ip1jm)
284      REAL dtvrs1j,rjour,heure,radsg,radomeg
285      REAL massebxy(ip1jm,llm)
286      INTEGER  l, ij, imjmp1
287
288      REAL       SSUM
289
290      logical  firstcal
291      data     firstcal/.true./
292      save     firstcal
293
294c-----------------------------------------------------------------------
295! Ehouarn: when no initialization fields from file, resetvarc should be
296!          set to false
297       if (firstcal) then
298         if (.not.read_start) then
299           resetvarc=.true.
300         endif
301       endif
302
303       dtvrs1j   = dtvr/daysec
304       rjour     = REAL( INT( itau * dtvrs1j ))
305       heure     = ( itau*dtvrs1j-rjour ) * 24.
306       imjmp1    = iim * jjp1
307       IF(ABS(heure - 24.).LE.0.0001 ) heure = 0.
308c
309       CALL massbarxy ( masse, massebxy )
310
311c   .....  Calcul  de  rmsdpdt  .....
312
313       ge(:)=dp(:)*dp(:)
314
315       rmsdpdt = SSUM(ip1jmp1,ge,1) - SSUM(jjp1,ge,iip1)
316c
317       rmsdpdt = daysec* 1.e-2 * SQRT(rmsdpdt/imjmp1) 
318
319       CALL SCOPY( ijp1llm,bern,1,bernf,1 )
320       CALL filtreg(bernf,jjp1,llm,-2,2,.TRUE.,1)
321
322c   .....  Calcul du moment  angulaire   .....
323
324       radsg    = rad /g
325       radomeg  = rad * omeg
326c
327       DO ij=iip2,ip1jm
328          cosphi( ij ) = COS(rlatu((ij-1)/iip1+1))
329          omegcosp(ij) = radomeg   * cosphi(ij)
330       ENDDO
331
332c  ...  Calcul  de l'energie,de l'enstrophie,de l'entropie et de rmsv  .
333
334       DO l=1,llm
335          DO ij = 1,ip1jm
336             vor(ij)=vorpot(ij,l)*vorpot(ij,l)*massebxy(ij,l)
337          ENDDO
338          ztotl(l)=(SSUM(ip1jm,vor,1)-SSUM(jjm,vor,iip1))
339
340          DO ij = 1,ip1jmp1
341             ge(ij)= masse(ij,l)*(phis(ij)+teta(ij,l)*pk(ij,l)  +
342     s        bernf(ij,l)-phi(ij,l))
343          ENDDO
344          etotl(l) = SSUM(ip1jmp1,ge,1) - SSUM(jjp1,ge,iip1)
345
346          DO   ij   = 1, ip1jmp1
347             ge(ij) = masse(ij,l)*teta(ij,l)
348          ENDDO
349          stotl(l)= SSUM(ip1jmp1,ge,1) - SSUM(jjp1,ge,iip1)
350
351          DO ij=1,ip1jmp1
352             ge(ij)=masse(ij,l)*AMAX1(bernf(ij,l)-phi(ij,l),0.)
353          ENDDO
354          rmsvl(l)=2.*(SSUM(ip1jmp1,ge,1)-SSUM(jjp1,ge,iip1))
355
356          DO ij =iip2,ip1jm
357             ge(ij)=(ucov(ij,l)/cu(ij)+omegcosp(ij))*masse(ij,l) *
358     *               cosphi(ij)
359          ENDDO
360          angl(l) = rad *
361     s    (SSUM(ip1jm-iip1,ge(iip2),1)-SSUM(jjm-1,ge(iip2),iip1))
362      ENDDO
363
364          DO ij=1,ip1jmp1
365            ge(ij)= ps(ij)*aire(ij)
366          ENDDO
367      ptot  = SSUM(ip1jmp1,ge,1)-SSUM(jjp1,ge,iip1)
368      etot  = SSUM(     llm, etotl, 1 )
369      ztot  = SSUM(     llm, ztotl, 1 )
370      stot  = SSUM(     llm, stotl, 1 )
371      rmsv  = SSUM(     llm, rmsvl, 1 )
372      ang   = SSUM(     llm,  angl, 1 )
373
374      IF (firstcal.and.resetvarc) then
375         PRINT 3500, itau, rjour, heure,time
376         PRINT*,'WARNING!!! On recalcule les valeurs initiales de :'
377         PRINT*,'ptot,rmsdpdt,etot,ztot,stot,rmsv,ang'
378         PRINT *, ptot,rmsdpdt,etot,ztot,stot,rmsv,ang
379         etot0 = etot
380         ptot0 = ptot
381         ztot0 = ztot
382         stot0 = stot
383         ang0  = ang
384      END IF
385
386      ! compute relative changes in etot,... (except if 'reference' values
387      ! are zero, which can happen when using iniacademic)
388      if (etot0.ne.0) then
389        etot= etot/etot0
390      else
391        etot=1.
392      endif
393      rmsv= SQRT(rmsv/ptot)
394      if (ptot0.ne.0) then
395        ptot= ptot/ptot0
396      else
397        ptot=1.
398      endif
399      if (ztot0.ne.0) then
400        ztot= ztot/ztot0
401      else
402        ztot=1.
403      endif
404      if (stot0.ne.0) then
405        stot= stot/stot0
406      else
407        stot=1.
408      endif
409      if (ang0.ne.0) then
410        ang = ang /ang0
411      else
412        ang=1.
413      endif
414
415      firstcal = .false.
416
417      PRINT 3500, itau, rjour, heure, time
418      PRINT 4000, ptot,rmsdpdt,etot,ztot,stot,rmsv,ang
419
420      RETURN
421
4223500   FORMAT(10("*"),4x,'pas',i7,5x,'jour',f9.0,'heure',f5.1,4x 
423     *   ,'date',f14.4,4x,10("*"))
4244000   FORMAT(10x,'masse',4x,'rmsdpdt',7x,'energie',2x,'enstrophie'
425     * ,2x,'entropie',3x,'rmsv',4x,'mt.ang',/,'GLOB  '
426     .  ,f10.6,e13.6,5f10.3/
427     * )
428      END
429
Note: See TracBrowser for help on using the repository browser.