source: LMDZ6/branches/IPSLCM6.0.13/libf/dyn3d/leapfrog.F @ 3012

Last change on this file since 3012 was 2622, checked in by Ehouarn Millour, 8 years ago

Some code tidying: turn ener.h into ener_mod.F90
EM

  • 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: 28.6 KB
Line 
1!
2! $Id: leapfrog.F 2622 2016-09-04 06:12:02Z fairhead $
3!
4c
5c
6      SUBROUTINE leapfrog(ucov,vcov,teta,ps,masse,phis,q,time_0)
7
8
9cIM : pour sortir les param. du modele dans un fis. netcdf 110106
10#ifdef CPP_IOIPSL
11      use IOIPSL
12#endif
13      USE infotrac, ONLY: nqtot,ok_iso_verif
14      USE guide_mod, ONLY : guide_main
15      USE write_field, ONLY: writefield
16      USE control_mod, ONLY: nday, day_step, planet_type, offline,
17     &                       iconser, iphysiq, iperiod, dissip_period,
18     &                       iecri, ip_ebil_dyn, ok_dynzon, ok_dyn_ins,
19     &                       periodav, ok_dyn_ave, output_grads_dyn
20      use exner_hyb_m, only: exner_hyb
21      use exner_milieu_m, only: exner_milieu
22      USE comvert_mod, ONLY: ap,bp,pressure_exner,presnivs
23      USE comconst_mod, ONLY: cpp, dtphys, dtvr, pi, ihf
24      USE logic_mod, ONLY: iflag_phys,ok_guide,forward,leapf,apphys,
25     &                     statcl,conser,apdiss,purmats,ok_strato
26      USE temps_mod, ONLY: jD_ref,jH_ref,itaufin,day_ini,day_ref,
27     &                        start_time,dt
28
29      IMPLICIT NONE
30
31c      ......   Version  du 10/01/98    ..........
32
33c             avec  coordonnees  verticales hybrides
34c   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
35
36c=======================================================================
37c
38c   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
39c   -------
40c
41c   Objet:
42c   ------
43c
44c   GCM LMD nouvelle grille
45c
46c=======================================================================
47c
48c  ... Dans inigeom , nouveaux calculs pour les elongations  cu , cv
49c      et possibilite d'appeler une fonction f(y)  a derivee tangente
50c      hyperbolique a la  place de la fonction a derivee sinusoidale.
51
52c  ... Possibilite de choisir le shema pour l'advection de
53c        q  , en modifiant iadv dans traceur.def  (10/02) .
54c
55c      Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
56c      Pour Van-Leer iadv=10
57c
58c-----------------------------------------------------------------------
59c   Declarations:
60c   -------------
61
62      include "dimensions.h"
63      include "paramet.h"
64      include "comdissnew.h"
65      include "comgeom.h"
66      include "description.h"
67      include "iniprint.h"
68      include "academic.h"
69
70      REAL,INTENT(IN) :: time_0 ! not used
71
72c   dynamical variables:
73      REAL,INTENT(INOUT) :: ucov(ip1jmp1,llm)    ! zonal covariant wind
74      REAL,INTENT(INOUT) :: vcov(ip1jm,llm)      ! meridional covariant wind
75      REAL,INTENT(INOUT) :: teta(ip1jmp1,llm)    ! potential temperature
76      REAL,INTENT(INOUT) :: ps(ip1jmp1)          ! surface pressure (Pa)
77      REAL,INTENT(INOUT) :: masse(ip1jmp1,llm)   ! air mass
78      REAL,INTENT(INOUT) :: phis(ip1jmp1)        ! geopotentiat at the surface
79      REAL,INTENT(INOUT) :: q(ip1jmp1,llm,nqtot) ! advected tracers
80
81      REAL p (ip1jmp1,llmp1  )               ! interlayer pressure
82      REAL pks(ip1jmp1)                      ! exner at the surface
83      REAL pk(ip1jmp1,llm)                   ! exner at mid-layer
84      REAL pkf(ip1jmp1,llm)                  ! filtered exner at mid-layer
85      REAL phi(ip1jmp1,llm)                  ! geopotential
86      REAL w(ip1jmp1,llm)                    ! vertical velocity
87
88      real zqmin,zqmax
89
90c variables dynamiques intermediaire pour le transport
91      REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse
92
93c   variables dynamiques au pas -1
94      REAL vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
95      REAL tetam1(ip1jmp1,llm),psm1(ip1jmp1)
96      REAL massem1(ip1jmp1,llm)
97
98c   tendances dynamiques
99      REAL dv(ip1jm,llm),du(ip1jmp1,llm)
100      REAL dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqtot),dp(ip1jmp1)
101
102c   tendances de la dissipation
103      REAL dvdis(ip1jm,llm),dudis(ip1jmp1,llm)
104      REAL dtetadis(ip1jmp1,llm)
105
106c   tendances physiques
107      REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
108      REAL dtetafi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqtot),dpfi(ip1jmp1)
109
110c   variables pour le fichier histoire
111      REAL dtav      ! intervalle de temps elementaire
112
113      REAL tppn(iim),tpps(iim),tpn,tps
114c
115      INTEGER itau,itaufinp1,iav
116!      INTEGER  iday ! jour julien
117      REAL       time
118
119      REAL  SSUM
120!     REAL finvmaold(ip1jmp1,llm)
121
122cym      LOGICAL  lafin
123      LOGICAL :: lafin=.false.
124      INTEGER ij,iq,l
125      INTEGER ik
126
127      real time_step, t_wrt, t_ops
128
129!      REAL rdayvrai,rdaym_ini
130! jD_cur: jour julien courant
131! jH_cur: heure julienne courante
132      REAL :: jD_cur, jH_cur
133      INTEGER :: an, mois, jour
134      REAL :: secondes
135
136      LOGICAL first,callinigrads
137cIM : pour sortir les param. du modele dans un fis. netcdf 110106
138      save first
139      data first/.true./
140      real dt_cum
141      character*10 infile
142      integer zan, tau0, thoriid
143      integer nid_ctesGCM
144      save nid_ctesGCM
145      real degres
146      real rlong(iip1), rlatg(jjp1)
147      real zx_tmp_2d(iip1,jjp1)
148      integer ndex2d(iip1*jjp1)
149      logical ok_sync
150      parameter (ok_sync = .true.)
151      logical physic
152
153      data callinigrads/.true./
154      character*10 string10
155
156      REAL :: flxw(ip1jmp1,llm)  ! flux de masse verticale
157
158c+jld variables test conservation energie
159      REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
160C     Tendance de la temp. potentiel d (theta)/ d t due a la
161C     tansformation d'energie cinetique en energie thermique
162C     cree par la dissipation
163      REAL dtetaecdt(ip1jmp1,llm)
164      REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
165      REAL vnat(ip1jm,llm),unat(ip1jmp1,llm)
166      REAL      d_h_vcol, d_qt, d_qw, d_ql, d_ec
167      CHARACTER*15 ztit
168!IM   INTEGER   ip_ebil_dyn  ! PRINT level for energy conserv. diag.
169!IM   SAVE      ip_ebil_dyn
170!IM   DATA      ip_ebil_dyn/0/
171c-jld
172
173      character*80 dynhist_file, dynhistave_file
174      character(len=*),parameter :: modname="leapfrog"
175      character*80 abort_message
176
177      logical dissip_conservative
178      save dissip_conservative
179      data dissip_conservative/.true./
180
181      LOGICAL prem
182      save prem
183      DATA prem/.true./
184      INTEGER testita
185      PARAMETER (testita = 9)
186
187      logical , parameter :: flag_verif = .false.
188     
189
190      integer itau_w   ! pas de temps ecriture = itap + itau_phy
191
192
193      if (nday>=0) then
194         itaufin   = nday*day_step
195      else
196         itaufin   = -nday
197      endif
198      itaufinp1 = itaufin +1
199      itau = 0
200      physic=.true.
201      if (iflag_phys==0.or.iflag_phys==2) physic=.false.
202
203c      iday = day_ini+itau/day_step
204c      time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
205c         IF(time.GT.1.) THEN
206c          time = time-1.
207c          iday = iday+1
208c         ENDIF
209
210
211c-----------------------------------------------------------------------
212c   On initialise la pression et la fonction d'Exner :
213c   --------------------------------------------------
214
215      dq(:,:,:)=0.
216      CALL pression ( ip1jmp1, ap, bp, ps, p       )
217      if (pressure_exner) then
218        CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
219      else
220        CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
221      endif
222
223c-----------------------------------------------------------------------
224c   Debut de l'integration temporelle:
225c   ----------------------------------
226
227   1  CONTINUE ! Matsuno Forward step begins here
228
229c   date: (NB: date remains unchanged for Backward step)
230c   -----
231
232      jD_cur = jD_ref + day_ini - day_ref +                             &
233     &          (itau+1)/day_step
234      jH_cur = jH_ref + start_time +                                    &
235     &          mod(itau+1,day_step)/float(day_step)
236      jD_cur = jD_cur + int(jH_cur)
237      jH_cur = jH_cur - int(jH_cur)
238
239        if (ok_iso_verif) then
240           call check_isotopes_seq(q,ip1jmp1,'leapfrog 321')
241        endif !if (ok_iso_verif) then
242
243#ifdef CPP_IOIPSL
244      if (ok_guide) then
245        call guide_main(itau,ucov,vcov,teta,q,masse,ps)
246      endif
247#endif
248
249
250c
251c     IF( MOD( itau, 10* day_step ).EQ.0 )  THEN
252c       CALL  test_period ( ucov,vcov,teta,q,p,phis )
253c       PRINT *,' ----   Test_period apres continue   OK ! -----', itau
254c     ENDIF
255c
256
257! Save fields obtained at previous time step as '...m1'
258      CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
259      CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
260      CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
261      CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
262      CALL SCOPY( ip1jmp1, ps  , 1,   psm1 , 1 )
263
264      forward = .TRUE.
265      leapf   = .FALSE.
266      dt      =  dtvr
267
268c   ...    P.Le Van .26/04/94  ....
269! Ehouarn: finvmaold is actually not used
270!      CALL SCOPY   ( ijp1llm,   masse, 1, finvmaold,     1 )
271!      CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
272
273        if (ok_iso_verif) then
274           call check_isotopes_seq(q,ip1jmp1,'leapfrog 400')
275        endif !if (ok_iso_verif) then
276
277   2  CONTINUE ! Matsuno backward or leapfrog step begins here
278
279c-----------------------------------------------------------------------
280
281c   date: (NB: only leapfrog step requires recomputing date)
282c   -----
283
284      IF (leapf) THEN
285        jD_cur = jD_ref + day_ini - day_ref +
286     &            (itau+1)/day_step
287        jH_cur = jH_ref + start_time +
288     &            mod(itau+1,day_step)/float(day_step)
289        jD_cur = jD_cur + int(jH_cur)
290        jH_cur = jH_cur - int(jH_cur)
291      ENDIF
292
293
294c   gestion des appels de la physique et des dissipations:
295c   ------------------------------------------------------
296c
297c   ...    P.Le Van  ( 6/02/95 )  ....
298
299      apphys = .FALSE.
300      statcl = .FALSE.
301      conser = .FALSE.
302      apdiss = .FALSE.
303
304      IF( purmats ) THEN
305      ! Purely Matsuno time stepping
306         IF( MOD(itau,iconser) .EQ.0.AND.  forward    ) conser = .TRUE.
307         IF( MOD(itau,dissip_period ).EQ.0.AND..NOT.forward )
308     s        apdiss = .TRUE.
309         IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
310     s          .and. physic                        ) apphys = .TRUE.
311      ELSE
312      ! Leapfrog/Matsuno time stepping
313         IF( MOD(itau   ,iconser) .EQ. 0              ) conser = .TRUE.
314         IF( MOD(itau+1,dissip_period).EQ.0 .AND. .NOT. forward )
315     s        apdiss = .TRUE.
316         IF( MOD(itau+1,iphysiq).EQ.0.AND.physic       ) apphys=.TRUE.
317      END IF
318
319! Ehouarn: for Shallow Water case (ie: 1 vertical layer),
320!          supress dissipation step
321      if (llm.eq.1) then
322        apdiss=.false.
323      endif
324
325
326        if (ok_iso_verif) then
327           call check_isotopes_seq(q,ip1jmp1,'leapfrog 589')
328        endif !if (ok_iso_verif) then
329
330c-----------------------------------------------------------------------
331c   calcul des tendances dynamiques:
332c   --------------------------------
333
334      ! compute geopotential phi()
335      CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi   )
336
337      time = jD_cur + jH_cur
338      CALL caldyn
339     $  ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
340     $    phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time )
341
342
343c-----------------------------------------------------------------------
344c   calcul des tendances advection des traceurs (dont l'humidite)
345c   -------------------------------------------------------------
346
347        if (ok_iso_verif) then
348           call check_isotopes_seq(q,ip1jmp1,
349     &           'leapfrog 686: avant caladvtrac')
350        endif !if (ok_iso_verif) then
351
352      IF( forward. OR . leapf )  THEN
353! Ehouarn: NB: fields sent to advtrac are those at the beginning of the time step
354         CALL caladvtrac(q,pbaru,pbarv,
355     *        p, masse, dq,  teta,
356     .        flxw, pk)
357          !write(*,*) 'caladvtrac 346'
358
359         
360         IF (offline) THEN
361Cmaf stokage du flux de masse pour traceurs OFF-LINE
362
363#ifdef CPP_IOIPSL
364           CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
365     .   dtvr, itau)
366#endif
367
368
369         ENDIF ! of IF (offline)
370c
371      ENDIF ! of IF( forward. OR . leapf )
372
373
374c-----------------------------------------------------------------------
375c   integrations dynamique et traceurs:
376c   ----------------------------------
377
378        if (ok_iso_verif) then
379           write(*,*) 'leapfrog 720'
380           call check_isotopes_seq(q,ip1jmp1,'leapfrog 756')
381        endif !if (ok_iso_verif) then
382       
383       CALL integrd ( nqtot,vcovm1,ucovm1,tetam1,psm1,massem1 ,
384     $         dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis )
385!     $              finvmaold                                    )
386
387       if (ok_iso_verif) then
388          write(*,*) 'leapfrog 724'
389           call check_isotopes_seq(q,ip1jmp1,'leapfrog 762')
390        endif !if (ok_iso_verif) then
391
392c .P.Le Van (26/04/94  ajout de  finvpold dans l'appel d'integrd)
393c
394c-----------------------------------------------------------------------
395c   calcul des tendances physiques:
396c   -------------------------------
397c    ########   P.Le Van ( Modif le  6/02/95 )   ###########
398c
399       IF( purmats )  THEN
400          IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
401       ELSE
402          IF( itau+1. EQ. itaufin )              lafin = .TRUE.
403       ENDIF
404c
405c
406       IF( apphys )  THEN
407c
408c     .......   Ajout   P.Le Van ( 17/04/96 )   ...........
409c
410
411         CALL pression (  ip1jmp1, ap, bp, ps,  p      )
412         if (pressure_exner) then
413           CALL exner_hyb(  ip1jmp1, ps, p,pks, pk, pkf )
414         else
415           CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
416         endif
417
418! Appel a geopot ajoute le 2014/05/08 pour garantir la convergence numerique
419! avec dyn3dmem
420         CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi   )
421
422!           rdaym_ini  = itau * dtvr / daysec
423!           rdayvrai   = rdaym_ini  + day_ini
424!           jD_cur = jD_ref + day_ini - day_ref
425!     $        + int (itau * dtvr / daysec)
426!           jH_cur = jH_ref +                                            &
427!     &              (itau * dtvr / daysec - int(itau * dtvr / daysec))
428           jD_cur = jD_ref + day_ini - day_ref +                        &
429     &          (itau+1)/day_step
430
431           IF (planet_type .eq."generic") THEN
432              ! AS: we make jD_cur to be pday
433              jD_cur = int(day_ini + itau/day_step)
434           ENDIF
435
436           jH_cur = jH_ref + start_time +                               &
437     &              mod(itau+1,day_step)/float(day_step)
438           jD_cur = jD_cur + int(jH_cur)
439           jH_cur = jH_cur - int(jH_cur)
440!         write(lunout,*)'itau, jD_cur = ', itau, jD_cur, jH_cur
441!         call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
442!         write(lunout,*)'current date = ',an, mois, jour, secondes
443
444c rajout debug
445c       lafin = .true.
446
447
448c   Inbterface avec les routines de phylmd (phymars ... )
449c   -----------------------------------------------------
450
451c+jld
452
453c  Diagnostique de conservation de l'énergie : initialisation
454         IF (ip_ebil_dyn.ge.1 ) THEN
455          ztit='bil dyn'
456! Ehouarn: be careful, diagedyn is Earth-specific!
457           IF (planet_type.eq."earth") THEN
458            CALL diagedyn(ztit,2,1,1,dtphys
459     &    , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
460           ENDIF
461         ENDIF ! of IF (ip_ebil_dyn.ge.1 )
462c-jld
463#ifdef CPP_IOIPSL
464cIM decommenter les 6 lignes suivantes pour sortir quelques parametres dynamiques de LMDZ
465cIM uncomment next 6 lines to get some parameters for LMDZ dynamics
466c        IF (first) THEN
467c         first=.false.
468c#include "ini_paramLMDZ_dyn.h"
469c        ENDIF
470c
471c#include "write_paramLMDZ_dyn.h"
472c
473#endif
474! #endif of #ifdef CPP_IOIPSL
475#ifdef CPP_PHYS
476         CALL calfis( lafin , jD_cur, jH_cur,
477     $               ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
478     $               du,dv,dteta,dq,
479     $               flxw,dufi,dvfi,dtetafi,dqfi,dpfi  )
480#endif
481c      ajout des tendances physiques:
482c      ------------------------------
483          CALL addfi( dtphys, leapf, forward   ,
484     $                  ucov, vcov, teta , q   ,ps ,
485     $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
486          ! since addfi updates ps(), also update p(), masse() and pk()
487          CALL pression (ip1jmp1,ap,bp,ps,p)
488          CALL massdair(p,masse)
489          if (pressure_exner) then
490            CALL exner_hyb(ip1jmp1,ps,p,pks,pk,pkf)
491          else
492            CALL exner_milieu(ip1jmp1,ps,p,pks,pk,pkf)
493          endif
494
495         IF (ok_strato) THEN
496           CALL top_bound( vcov,ucov,teta,masse,dtphys)
497         ENDIF
498       
499c
500c  Diagnostique de conservation de l'énergie : difference
501         IF (ip_ebil_dyn.ge.1 ) THEN
502          ztit='bil phys'
503          IF (planet_type.eq."earth") THEN
504           CALL diagedyn(ztit,2,1,1,dtphys
505     &     , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
506          ENDIF
507         ENDIF ! of IF (ip_ebil_dyn.ge.1 )
508
509       ENDIF ! of IF( apphys )
510
511      IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
512!   Academic case : Simple friction and Newtonan relaxation
513!   -------------------------------------------------------
514        DO l=1,llm   
515          DO ij=1,ip1jmp1
516           teta(ij,l)=teta(ij,l)-dtvr*
517     &      (teta(ij,l)-tetarappel(ij,l))*(knewt_g+knewt_t(l)*clat4(ij))
518          ENDDO
519        ENDDO ! of DO l=1,llm
520       
521        if (planet_type.eq."giant") then
522          ! add an intrinsic heat flux at the base of the atmosphere
523          teta(:,1)=teta(:,1)+dtvr*aire(:)*ihf/cpp/masse(:,1)
524        endif
525
526        call friction(ucov,vcov,dtvr)
527       
528        ! Sponge layer (if any)
529        IF (ok_strato) THEN
530!          dufi(:,:)=0.
531!          dvfi(:,:)=0.
532!          dtetafi(:,:)=0.
533!          dqfi(:,:,:)=0.
534!          dpfi(:)=0.
535!          CALL top_bound(vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
536           CALL top_bound( vcov,ucov,teta,masse,dtvr)
537!          CALL addfi( dtvr, leapf, forward   ,
538!     $                  ucov, vcov, teta , q   ,ps ,
539!     $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
540        ENDIF ! of IF (ok_strato)
541      ENDIF ! of IF (iflag_phys.EQ.2)
542
543
544c-jld
545
546        CALL pression ( ip1jmp1, ap, bp, ps, p                  )
547        if (pressure_exner) then
548          CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
549        else
550          CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
551        endif
552        CALL massdair(p,masse)
553
554        if (ok_iso_verif) then
555           call check_isotopes_seq(q,ip1jmp1,'leapfrog 1196')
556        endif !if (ok_iso_verif) then
557
558c-----------------------------------------------------------------------
559c   dissipation horizontale et verticale  des petites echelles:
560c   ----------------------------------------------------------
561
562      IF(apdiss) THEN
563
564
565c   calcul de l'energie cinetique avant dissipation
566        call covcont(llm,ucov,vcov,ucont,vcont)
567        call enercin(vcov,ucov,vcont,ucont,ecin0)
568
569c   dissipation
570        CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
571        ucov=ucov+dudis
572        vcov=vcov+dvdis
573c       teta=teta+dtetadis
574
575
576c------------------------------------------------------------------------
577        if (dissip_conservative) then
578C       On rajoute la tendance due a la transform. Ec -> E therm. cree
579C       lors de la dissipation
580            call covcont(llm,ucov,vcov,ucont,vcont)
581            call enercin(vcov,ucov,vcont,ucont,ecin)
582            dtetaecdt= (ecin0-ecin)/ pk
583c           teta=teta+dtetaecdt
584            dtetadis=dtetadis+dtetaecdt
585        endif
586        teta=teta+dtetadis
587c------------------------------------------------------------------------
588
589
590c    .......        P. Le Van (  ajout  le 17/04/96  )   ...........
591c   ...      Calcul de la valeur moyenne, unique de h aux poles  .....
592c
593
594        DO l  =  1, llm
595          DO ij =  1,iim
596           tppn(ij)  = aire(  ij    ) * teta(  ij    ,l)
597           tpps(ij)  = aire(ij+ip1jm) * teta(ij+ip1jm,l)
598          ENDDO
599           tpn  = SSUM(iim,tppn,1)/apoln
600           tps  = SSUM(iim,tpps,1)/apols
601
602          DO ij = 1, iip1
603           teta(  ij    ,l) = tpn
604           teta(ij+ip1jm,l) = tps
605          ENDDO
606        ENDDO
607
608        if (1 == 0) then
609!!! Ehouarn: lines here 1) kill 1+1=2 in the dynamics
610!!!                     2) should probably not be here anyway
611!!! but are kept for those who would want to revert to previous behaviour
612           DO ij =  1,iim
613             tppn(ij)  = aire(  ij    ) * ps (  ij    )
614             tpps(ij)  = aire(ij+ip1jm) * ps (ij+ip1jm)
615           ENDDO
616             tpn  = SSUM(iim,tppn,1)/apoln
617             tps  = SSUM(iim,tpps,1)/apols
618
619           DO ij = 1, iip1
620             ps(  ij    ) = tpn
621             ps(ij+ip1jm) = tps
622           ENDDO
623        endif ! of if (1 == 0)
624
625      END IF ! of IF(apdiss)
626
627c ajout debug
628c              IF( lafin ) then 
629c                abort_message = 'Simulation finished'
630c                call abort_gcm(modname,abort_message,0)
631c              ENDIF
632       
633c   ********************************************************************
634c   ********************************************************************
635c   .... fin de l'integration dynamique  et physique pour le pas itau ..
636c   ********************************************************************
637c   ********************************************************************
638
639c   preparation du pas d'integration suivant  ......
640
641        if (ok_iso_verif) then
642           call check_isotopes_seq(q,ip1jmp1,'leapfrog 1509')
643        endif !if (ok_iso_verif) then
644
645      IF ( .NOT.purmats ) THEN
646c       ........................................................
647c       ..............  schema matsuno + leapfrog  ..............
648c       ........................................................
649
650            IF(forward. OR. leapf) THEN
651              itau= itau + 1
652c              iday= day_ini+itau/day_step
653c              time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
654c                IF(time.GT.1.) THEN
655c                  time = time-1.
656c                  iday = iday+1
657c                ENDIF
658            ENDIF
659
660
661            IF( itau. EQ. itaufinp1 ) then 
662              if (flag_verif) then
663                write(79,*) 'ucov',ucov
664                write(80,*) 'vcov',vcov
665                write(81,*) 'teta',teta
666                write(82,*) 'ps',ps
667                write(83,*) 'q',q
668                WRITE(85,*) 'q1 = ',q(:,:,1)
669                WRITE(86,*) 'q3 = ',q(:,:,3)
670              endif
671
672              abort_message = 'Simulation finished'
673
674              call abort_gcm(modname,abort_message,0)
675            ENDIF
676c-----------------------------------------------------------------------
677c   ecriture du fichier histoire moyenne:
678c   -------------------------------------
679
680            IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
681               IF(itau.EQ.itaufin) THEN
682                  iav=1
683               ELSE
684                  iav=0
685               ENDIF
686               
687!              ! Ehouarn: re-compute geopotential for outputs
688               CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
689
690               IF (ok_dynzon) THEN
691#ifdef CPP_IOIPSL
692                 CALL bilan_dyn(2,dtvr*iperiod,dtvr*day_step*periodav,
693     &                 ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
694#endif
695               END IF
696               IF (ok_dyn_ave) THEN
697#ifdef CPP_IOIPSL
698                 CALL writedynav(itau,vcov,
699     &                 ucov,teta,pk,phi,q,masse,ps,phis)
700#endif
701               ENDIF
702
703            ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
704
705        if (ok_iso_verif) then
706           call check_isotopes_seq(q,ip1jmp1,'leapfrog 1584')
707        endif !if (ok_iso_verif) then
708
709c-----------------------------------------------------------------------
710c   ecriture de la bande histoire:
711c   ------------------------------
712
713            IF( MOD(itau,iecri).EQ.0) THEN
714             ! Ehouarn: output only during LF or Backward Matsuno
715             if (leapf.or.(.not.leapf.and.(.not.forward))) then
716              CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
717              unat=0.
718              do l=1,llm
719                unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
720                vnat(:,l)=vcov(:,l)/cv(:)
721              enddo
722#ifdef CPP_IOIPSL
723              if (ok_dyn_ins) then
724!               write(lunout,*) "leapfrog: call writehist, itau=",itau
725               CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
726!               call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
727!               call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
728!              call WriteField('teta',reshape(teta,(/iip1,jmp1,llm/)))
729!               call WriteField('ps',reshape(ps,(/iip1,jmp1/)))
730!               call WriteField('masse',reshape(masse,(/iip1,jmp1,llm/)))
731              endif ! of if (ok_dyn_ins)
732#endif
733! For some Grads outputs of fields
734              if (output_grads_dyn) then
735#include "write_grads_dyn.h"
736              endif
737             endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
738            ENDIF ! of IF(MOD(itau,iecri).EQ.0)
739
740            IF(itau.EQ.itaufin) THEN
741
742
743!              if (planet_type.eq."earth") then
744! Write an Earth-format restart file
745                CALL dynredem1("restart.nc",start_time,
746     &                         vcov,ucov,teta,q,masse,ps)
747!              endif ! of if (planet_type.eq."earth")
748
749              CLOSE(99)
750              !!! Ehouarn: Why not stop here and now?
751            ENDIF ! of IF (itau.EQ.itaufin)
752
753c-----------------------------------------------------------------------
754c   gestion de l'integration temporelle:
755c   ------------------------------------
756
757            IF( MOD(itau,iperiod).EQ.0 )    THEN
758                    GO TO 1
759            ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
760
761                   IF( forward )  THEN
762c      fin du pas forward et debut du pas backward
763
764                      forward = .FALSE.
765                        leapf = .FALSE.
766                           GO TO 2
767
768                   ELSE
769c      fin du pas backward et debut du premier pas leapfrog
770
771                        leapf =  .TRUE.
772                        dt  =  2.*dtvr
773                        GO TO 2
774                   END IF ! of IF (forward)
775            ELSE
776
777c      ......   pas leapfrog  .....
778
779                 leapf = .TRUE.
780                 dt  = 2.*dtvr
781                 GO TO 2
782            END IF ! of IF (MOD(itau,iperiod).EQ.0)
783                   !    ELSEIF (MOD(itau-1,iperiod).EQ.0)
784
785      ELSE ! of IF (.not.purmats)
786
787        if (ok_iso_verif) then
788           call check_isotopes_seq(q,ip1jmp1,'leapfrog 1664')
789        endif !if (ok_iso_verif) then
790
791c       ........................................................
792c       ..............       schema  matsuno        ...............
793c       ........................................................
794            IF( forward )  THEN
795
796             itau =  itau + 1
797c             iday = day_ini+itau/day_step
798c             time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
799c
800c                  IF(time.GT.1.) THEN
801c                   time = time-1.
802c                   iday = iday+1
803c                  ENDIF
804
805               forward =  .FALSE.
806               IF( itau. EQ. itaufinp1 ) then 
807                 abort_message = 'Simulation finished'
808                 call abort_gcm(modname,abort_message,0)
809               ENDIF
810               GO TO 2
811
812            ELSE ! of IF(forward) i.e. backward step
813 
814        if (ok_iso_verif) then
815           call check_isotopes_seq(q,ip1jmp1,'leapfrog 1698')
816        endif !if (ok_iso_verif) then 
817
818              IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
819               IF(itau.EQ.itaufin) THEN
820                  iav=1
821               ELSE
822                  iav=0
823               ENDIF
824
825!              ! Ehouarn: re-compute geopotential for outputs
826               CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
827
828               IF (ok_dynzon) THEN
829#ifdef CPP_IOIPSL
830                 CALL bilan_dyn(2,dtvr*iperiod,dtvr*day_step*periodav,
831     &                 ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
832#endif
833               ENDIF
834               IF (ok_dyn_ave) THEN
835#ifdef CPP_IOIPSL
836                 CALL writedynav(itau,vcov,
837     &                 ucov,teta,pk,phi,q,masse,ps,phis)
838#endif
839               ENDIF
840
841              ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
842
843              IF(MOD(itau,iecri         ).EQ.0) THEN
844c              IF(MOD(itau,iecri*day_step).EQ.0) THEN
845                CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
846                unat=0.
847                do l=1,llm
848                  unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
849                  vnat(:,l)=vcov(:,l)/cv(:)
850                enddo
851#ifdef CPP_IOIPSL
852              if (ok_dyn_ins) then
853!                write(lunout,*) "leapfrog: call writehist (b)",
854!     &                        itau,iecri
855                CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
856              endif ! of if (ok_dyn_ins)
857#endif
858! For some Grads outputs
859                if (output_grads_dyn) then
860#include "write_grads_dyn.h"
861                endif
862
863              ENDIF ! of IF(MOD(itau,iecri         ).EQ.0)
864
865              IF(itau.EQ.itaufin) THEN
866!                if (planet_type.eq."earth") then
867                  CALL dynredem1("restart.nc",start_time,
868     &                           vcov,ucov,teta,q,masse,ps)
869!                endif ! of if (planet_type.eq."earth")
870              ENDIF ! of IF(itau.EQ.itaufin)
871
872              forward = .TRUE.
873              GO TO  1
874
875            ENDIF ! of IF (forward)
876
877      END IF ! of IF(.not.purmats)
878
879      STOP
880      END
Note: See TracBrowser for help on using the repository browser.