source: trunk/LMDZ.PLUTO.old/libf/dyn3d/nogcm.F_saved

Last change on this file was 3175, checked in by emillour, 23 months ago

Pluto PCM:
Add the old Pluto LMDZ for reference (required prior step to making
an LMDZ.PLUTO using the same framework as the other physics packages).
TB+EM
File size: 14.8 KB
Line 
1      PROGRAM nogcm
2
3      IMPLICIT NONE
4
5c      ......   Version  du 01/09/15    ..........
6
7c             avec  coordonnees  verticales hybrides
8
9c=======================================================================
10c
11c   Auteur:  F. Forget, T. Bertrand
12c   -------
13c
14c   Objet:
15c   ------
16c
17c   GCM LMD sans dynamique, pour modele saisonnier de surface
18c   
19c=======================================================================
20c
21c-----------------------------------------------------------------------
22c   Declarations:
23c   -------------
24
25#include "dimensions.h"
26#include "paramet.h"
27#include "comconst.h"
28#include "comdissnew.h"
29#include "comvert.h"
30#include "comgeom.h"
31#include "logic.h"
32#include "temps.h"
33#include "control.h"
34#include "ener.h"
35#include "netcdf.inc"
36#include "description.h"
37#include "serre.h"
38#include "tracstoke.h"
39#include "sponge.h"
40#include "advtrac.h"
41#include "callkeys.h"
42#include "tracer.h"
43
44      INTEGER*4  iday ! jour julien
45      REAL       time ! Heure de la journee en fraction d''1 jour
46      REAL zdtvr
47
48
49c   variables dynamiques
50      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
51      real, dimension(ip1jmp1,llm) :: teta   ! temperature potentielle
52      REAL q(ip1jmp1,llm,nqmx)               ! champs advectes
53      REAL ps(ip1jmp1)                       ! pression  au sol
54      REAL kp(ip1jmp1)                       ! TB15 = exp (-z/H)
55      REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
56      REAL pks(ip1jmp1)                      ! exner au  sol
57      REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
58      REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
59      REAL masse(ip1jmp1,llm)                ! masse d''air
60      REAL phis(ip1jmp1)                     ! geopotentiel au sol
61      REAL phi(ip1jmp1,llm)                  ! geopotentiel
62      REAL w(ip1jmp1,llm)                    ! vitesse verticale
63
64c   tendances dynamiques
65      REAL dv(ip1jm,llm),du(ip1jmp1,llm)
66      REAL dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqmx)
67
68c   tendances physiques
69      REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
70      REAL dhfi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqmx),dpfi(ip1jmp1)
71
72c   variables pour le fichier histoire
73      REAL dtav      ! intervalle de temps elementaire
74
75      INTEGER itau,itaufinp1
76
77      EXTERNAL traceur
78      EXTERNAL iniconst
79      EXTERNAL SCOPY
80      EXTERNAL inigeom
81      EXTERNAL exner_hyb,addit
82      EXTERNAL defrun_new, test_period
83      REAL time_0
84
85      LOGICAL lafin
86      INTEGER ij
87
88      REAL rdayvrai,rdaym_ini,rday_ecri
89      REAL beta(ip1jmp1,llm)
90
91      character*20 modname
92      character*80 abort_message
93
94      INTEGER nq, numvanle
95
96      REAL psmean                    ! pression moyenne
97      REAL p0                        ! pression de reference
98      REAL p00                        ! globalaverage(kp)
99      REAL qmean_ch4                    ! mass mean mixing ratio vap ch4
100      REAL qmean_co                    ! mass mean mixing ratio vap co
101      REAL pqch4(ip1jmp1)           ! average methane mass index : ps*qch4
102      REAL pqco(ip1jmp1)           ! average methane mass index : ps*q_co
103      REAL oldps(ip1jmp1)           ! saving old pressure ps to calculate qch4
104
105! flag to set/remove calls to groupeun
106      real globaverage2d    !tau_n2, tau_ch4,tau_co
107      logical firstphys
108      data firstphys/.true./
109      save firstphys
110c-----------------------------------------------------------------------
111c   Initialisations:
112c   ----------------
113
114      modname = 'gcm'
115      descript = 'Run GCM LMDZ'
116      lafin    = .FALSE.
117
118c-----------------------------------------------------------------------
119c  Initialize tracers using iniadvtrac (Ehouarn, oct 2008)
120      CALL iniadvtrac(nq,numvanle)
121
122      CALL dynetat0("start.nc",nqmx,vcov,ucov,
123     .              teta,q,masse,ps,phis, time_0)
124      CALL defrun_new( 99, .TRUE. )
125      write(*,*) 'TB15: test0 r, cpp=',r,cpp
126
127c  on recalcule eventuellement le pas de temps
128
129      IF(MOD(day_step,iperiod).NE.0)
130     * STOP'Il faut choisir un nb de pas par jour multiple de iperiod'
131
132c      IF(MOD(day_step,iphysiq).NE.0)
133c     * STOP'Il faut choisir un nb de pas par jour multiple de iphysiq'
134
135      zdtvr    = daysec/FLOAT(day_step)
136
137        IF(dtvr.NE.zdtvr) THEN
138         PRINT*,'WARNING!!! changement de pas de temps',dtvr,'>',zdtvr
139        ENDIF
140
141c  nombre d etats dans les fichiers demarrage et histoire
142
143      dtvr = zdtvr
144      CALL iniconst
145      CALL inigeom
146c      CALL inifilr
147
148c   ......    P.Le Van    ( modif  le 29/04/97 )   ......... 
149
150c      CALL inidissip ( lstardis, nitergdiv, nitergrot, niterh   ,
151c     *                tetagdiv, tetagrot , tetatemp              )
152
153      CALL pression ( ip1jmp1, ap, bp, ps, p       )
154      call dump2d(iip1,jjp1,ps,'PRESSION SURFACE')
155      CALL exner_hyb( ip1jmp1, ps, p,beta, pks, pk, pkf )
156
157c
158c  numero de stockage pour les fichiers de redemarrage:
159
160c-----------------------------------------------------------------------
161c   temps de depart et de fin:
162c   --------------------------
163
164      itau = 0
165      iday = day_ini+itau/day_step
166      time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
167
168
169         IF(time.GT.1.) THEN
170          time = time-1.
171          iday = iday+1
172         ENDIF
173      itaufin   = nday*day_step
174      itaufinp1 = itaufin +1
175
176      day_end = day_ini + nday
177      write(*,*) "TB16:",day_ini,nday,day_end
178      PRINT 300, itau,itaufin,day_ini,day_end
179 300  FORMAT('1'/,15x,'run du pas',i7,2x,'au pas',i7,2x, 
180     . 'c''est a dire du jour',i7,3x,'au jour',i7//)
181
182      CALL dynredem0("restart.nc",day_end,anne_ini,phis,nqmx)
183 
184      ecripar = .TRUE.
185
186      dtav = iperiod*dtvr/daysec
187
188
189c   Quelques initialisations pour les traceurs
190      call initial0(ijp1llm*nqmx,dq)
191c     istdyn=day_step/4     ! stockage toutes les 6h=1jour/4
192c     istphy=istdyn/iphysiq
193
194c      CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi   )
195
196c-----------------------------------------------------------------------
197c   Debut de l integration temporelle:
198c   ----------------------------------
199
200c     kp=exp(-z/H)  Needed to define a reference pressure :
201c     P0=pmean/globalaverage(kp) 
202c     thus P(i) = P0*exp(-z(i)/H) = P0*kp(i)
203c     It is checked that globalaverage2d(Pi)=pmean
204      DO ij=1,ip1jmp1
205             kp(ij) = exp(-phis(ij)/(r*38.))
206      ENDDO
207      p00=globaverage2d(kp) ! mean pres at ref level
208
209c      tau_n2 = 1. ! constante de rappel for pressure n2 (s)
210c      tau_ch4 = 1.E7 ! constante de rappel for mix ratio qch4 (s)
211c      tau_co = 1. !E5 ! constante de rappel for mix ratio qco (s)
212
213      dt      =  dtvr
214      dtphys  = iphysiq * dtvr
215      DO itau = 1,itaufin
216
217        apphys = ( MOD(itau+1,iphysiq ).EQ.0.AND. physic)
218c-----------------------------------------------------------------------
219c   calcul des tendances physiques:
220c -------------------------------
221        IF( itau+1.EQ.itaufin) lafin = .TRUE.
222c
223        IF( apphys )  THEN
224           rdaym_ini  = (itau) * dtvr / daysec
225           rdayvrai   = rdaym_ini  + day_ini
226           IF ( ecritphy.LT.1. )  THEN
227             rday_ecri = rdaym_ini
228           ELSE
229             rday_ecri = INT(rdaym_ini)+INT(day_ini)
230           ENDIF
231c           write(20,*) 'TB15 : ps1 nogcm= ',globaverage2d(ps)  ! mean pressure
232
233           CALL calfis( nqmx, lafin ,rdayvrai,rday_ecri,time  ,
234     $                 ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
235     $     du,dv,dteta,dq,w,dufi,dvfi,dhfi,dqfi,dpfi,tracer)
236c           write(21,*) 'TB15 : ps2 nogcm= ',globaverage2d(ps)  ! mean pressure
237
238c          saving pressure pplev :
239           DO ij=1,ip1jmp1
240               oldps(ij)=ps(ij)             
241           ENDDO
242
243c          increment q_ch4 with physical tendancy
244           if (methane) then
245             DO ij=1,ip1jmp1
246                 q(ij,1,igcm_ch4_gas)=q(ij,1,igcm_ch4_gas)+
247     &                    dqfi(ij,1,igcm_ch4_gas)*dtphys
248             ENDDO
249c             write(31,*) 'TB15 : qch4 t1 = ',
250c     &                     globaverage2d(q(:,1,igcm_ch4_gas)*ps(:)/g)
251           endif
252
253c          increment q_co with physical tendancy
254           if (carbox) then
255             DO ij=1,ip1jmp1
256                 q(ij,1,igcm_co_gas)=q(ij,1,igcm_co_gas)+
257     &                    dqfi(ij,1,igcm_co_gas)*dtphys
258             ENDDO
259!             write(41,*) 'TB15 : qco t1 = ',
260!     &                     globaverage2d(q(:,1,igcm_co_gas)*ps(:)/g)
261           ENDIF
262
263c          update pressure
264           DO ij=1,ip1jmp1
265             ps(ij) = ps(ij) + dpfi(ij)*dtphys
266           ENDDO
267
268           psmean= globaverage2d(ps)  ! mean pressure
269           p0=psmean/p00
270           DO ij=1,ip1jmp1
271             ! Rappel newtonien vers psmean
272             ps(ij)=ps(ij) +(p0*kp(ij)
273     &                  -ps(ij))*(1-exp(-dtphys/tau_n2))
274           ENDDO
275
276           write(72,*) 'TB15 : ps redistributed = ',psmean  ! mean pressure
277           write(72,*) ' '
278
279c          TB15 : test pressure negative
280           DO ij=1,ip1jmp1
281             IF (ps(ij).le.0.) then
282                 write(*,*) 'Negative pressure :'
283                 write(*,*) 'ps = ',ps(ij),' ig = ',ij
284                 write(*,*) 'pmean = ',p0*kp(ij)
285                 write(*,*) 'tau_n2 = ',tau_n2
286                 STOP
287             ENDIF
288           ENDDO
289
290           if (methane) then
291             ! Simulate redistribution by dynamics for qch4
292             DO ij=1,ip1jmp1
293c              update mmr q for new pressure
294c              (both q and dqfi were calculated with old pressure in physiq)
295               q(ij,1,igcm_ch4_gas)=q(ij,1,igcm_ch4_gas)*oldps(ij)/
296     &                                                           ps(ij)
297               ! average methane mass index : ps * qch4
298               pqch4(ij)= ps(ij) * q(ij,1,igcm_ch4_gas)
299             end do
300c             write(32,*) 'TB15 : qch4 t2 = ',
301c     &                     globaverage2d(q(:,1,igcm_ch4_gas)*ps(:)/g)
302             ! weighted averaged CH4 mass mixing ratio:
303             qmean_ch4= globaverage2d(pqch4) / globaverage2d(ps)
304             DO ij=1,ip1jmp1
305               ! Rappel newtonien vers qmean
306               q(ij,1,igcm_ch4_gas)=q(ij,1,igcm_ch4_gas) +(qmean_ch4
307     &                  -q(ij,1,igcm_ch4_gas))*(1.-exp(-dtphys/tau_ch4))
308             ENDDO
309c             write(33,*) 'TB15 : qch4 t3 = ',
310c     &                     globaverage2d(q(:,1,igcm_ch4_gas)*ps(:)/g)
311           endif
312
313           if (carbox) then
314             ! Simulate redistribution by dynamics for q_co
315             DO ij=1,ip1jmp1
316c              update mmr q for new pressure
317c              (both q and dqfi were calculated with old pressure in physiq)
318               q(ij,1,igcm_co_gas)=q(ij,1,igcm_co_gas)*oldps(ij)/
319     &                                                           ps(ij)
320               ! average methane mass index : ps * q_co
321               pqco(ij)= ps(ij) * q(ij,1,igcm_co_gas)
322             end do
323!             write(42,*) 'TB15 : qco t2 = ',
324!     &                     globaverage2d(q(:,1,igcm_co_gas)*ps(:)/g)
325             ! weighted averaged CO mass mixing ratio:
326             qmean_co= globaverage2d(pqco) / globaverage2d(ps)
327             DO ij=1,ip1jmp1
328               ! Rappel newtonien vers qmean
329               q(ij,1,igcm_co_gas)=q(ij,1,igcm_co_gas) +(qmean_co
330     &                  -q(ij,1,igcm_co_gas))*(1.-exp(-dtphys/tau_co))
331             ENDDO
332!             write(43,*) 'TB15 : qco t3 = ',
333!     &                     globaverage2d(q(:,1,igcm_co_gas)*ps(:)/g)
334           endif
335
336ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
337
338
339           CALL pression ( ip1jmp1, ap, bp, ps, p                  )
340           CALL massdair (     p  , masse         )
341           CALL exner_hyb( ip1jmp1, ps, p,beta, pks, pk, pkf )
342
343          if(firstphys) then
344
345c            Option Special "nogcm"
346c            ---------------
347             corrk      = .false.
348             calldifv = .false.
349             calladj  = .false.
350             callconduct=.false.
351             IF (paleo) then
352              day_end=int(day_ini+paleoyears*365.25/6.3872)+1
353              CALL dynredem0("restart.nc",day_end,anne_ini,phis,nqmx)
354             ENDIF
355
356             write(*,*) '**********************************************'
357             write(*,*) 'WARNING: For nogcm, these options are forced:'
358             write(*,*) '**********************************************'
359             write(*,*) 'Option corrk = ', corrk
360             write(*,*) 'Option  calldifv = ',  calldifv
361             write(*,*) 'Option  calladj = ', calladj
362             write(*,*) 'Option  callconduct = ', callconduct
363             write(*,*) 'tau N2 - CH4 - CO = ', tau_n2,tau_ch4,tau_co
364             write(*,*) '**********************************************'
365             firstphys=.false.
366          end if
367        ENDIF
368
369c   ********************************************************************
370c   .... fin de l'integration dynamique  et physique pour le pas itau ..
371c   ********************************************************************
372
373c       preparation du pas d'integration suivant  ......
374
375        iday= day_ini+itau/day_step
376        time= FLOAT(itau+1-(iday-day_ini)*day_step)/day_step+time_0
377        IF(time.GT.1.) THEN
378             time = time-1.
379             iday = iday+1
380        ENDIF
381       
382c-----------------------------------------------------------------------
383
384        IF(itau.EQ.itaufin) THEN
385           PRINT *,' Appel test_period avant redem ', itau
386
387c           TBch4 : coupe test period car nous embete pour des pouilleme
388c           sur le ch4_gas
389           CALL test_period ( ucov,vcov,teta,q,p,phis )
390           CALL dynredem1("restart.nc",0.0,
391     .                     vcov,ucov,teta,q,nqmx,masse,ps)
392              CLOSE(99)
393        ENDIF
394
395      END DO   ! fin de la boucle temporelle
396      abort_message = 'Simulation finished'
397      call abort_gcm(modname,abort_message,0)
398
399      STOP
400      END
401
402!*********************************************************************************
403      FUNCTION globaverage2d(var)
404      IMPLICIT NONE
405#include "dimensions.h"
406#include "paramet.h"
407#include "comgeom2.h"
408
409       REAL var(iip1,jjp1), globaverage2d
410       integer i,j
411       REAL airetot
412       save airetot
413       logical firstcall
414       data firstcall/.true./
415       save firstcall
416
417       if (firstcall) then
418         firstcall=.false.
419c         write(30,*) 'TB15 : aire = ',aire
420         airetot =0.
421         DO j=2,jjp1-1      ! lat
422           DO i = 1, iim    ! lon
423             airetot = airetot + aire(i,j)
424           END DO
425         END DO
426         airetot=airetot + aire(1,1)+aire(1,jjp1)
427c         write(*,*) 'TB15 : nogcm totarea= ',airetot
428       end if
429
430       globaverage2d = 0.
431       DO j=2,jjp1-1
432         DO i = 1, iim
433           globaverage2d = globaverage2d + var(i,j)*aire(i,j)
434         END DO
435       END DO
436       globaverage2d = globaverage2d + var(1,1)*aire(1,1)
437       globaverage2d = globaverage2d + var(1,jjp1)*aire(1,jjp1)
438       globaverage2d = globaverage2d/airetot
439      return
440      end 
Note: See TracBrowser for help on using the repository browser.