source: trunk/LMDZ.PLUTO.old/libf/dyn3d/nogcm.F @ 3436

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