source: LMDZ.3.3/branches/LF/libf/dyn3d/offline.F

Last change on this file was 2, checked in by lmdz, 25 years ago

Initial revision

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 19.0 KB
Line 
1      PROGRAM offline
2      USE ioipsl
3      IMPLICIT NONE
4
5c      ......   Version  du 17/04/96    ..........
6c=======================================================================
7c
8c   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
9c   -------
10c   Modif special traceur F.Forget 05/94
11c
12c   Modifs M-A Filiberti 07/99
13C
14c   Objet:
15c   ------
16c
17c   GCM LMD nouvelle grille
18c
19c=======================================================================
20
21c   ...      modification de l'integration de q ( 26/04/94 )      ....
22
23c-----------------------------------------------------------------------
24c   Declarations:
25c   -------------
26
27#include "dimensions.h"
28#include "paramet.h"
29#include "comconst.h"
30#include "comdissnew.h"
31#include "comvert.h"
32#include "comgeom2.h"
33#include "logic.h"
34#include "temps.h"
35#include "control.h"
36#include "ener.h"
37#include "description.h"
38#include "serre.h"
39#include "tracstoke.h"
40      integer ngridmx,nbsrf
41      PARAMETER (ngridmx=iim*(jjm-1)+2,nbsrf=4)
42c-----------------------------------------------------------------------
43
44      INTEGER*4  iday ! jour julien
45      REAL       time ! Heure de la journee en fraction d'1 jour
46      REAL dtdyn,dtav
47      INTEGER iday_step
48
49      INTEGER         longcles
50      PARAMETER     ( longcles = 20 )
51      REAL  clesphy0( longcles )
52      logical clefinit
53
54
55c   variables dynamiques
56      REAL vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants
57      REAL q(iip1,jjp1,llm,nqmx)               ! champs advectes
58      REAL pks(iip1,jjp1)                      ! exner (f pour filtre)
59      REAL phis(iip1,jjp1)                     ! geopotentiel au sol
60      REAL phi(iip1,jjp1,llm)                  ! geopotentiel
61      REAL w(iip1,jjp1,llm)                    ! vitesse verticale
62C
63C   variables physiques
64      real zmfd(ngridmx,llm),zde_d(ngridmx,llm),zen_d(ngridmx,llm)
65      REAL qfi(ngridmx,llm,nqmx)               ! champs advectes
66      real zmfu(ngridmx,llm),zde_u(ngridmx,llm),zen_u(ngridmx,llm)
67      real coefkz(ngridmx,llm)
68      real frac_impa(ngridmx,llm),frac_nucl(ngridmx,llm)
69      real pphis (ngridmx)
70      REAL airefi(ngridmx)
71      REAL latfi(ngridmx),lonfi(ngridmx) ! latitude et long pour chaque point
72      real yu1(ngridmx), yv1(ngridmx)
73      real ftsol(ngridmx,nbsrf),pctsrf(ngridmx,nbsrf)
74
75C
76      character*10 file
77
78c variables dynamiques intermediaire pour le transport
79      REAL pbaru(iip1,jjp1,llm),pbarv(iip1,jjm,llm) !flux de masse
80
81      REAL masse(iip1,jjp1,llm)
82      REAL massebx(iip1,jjp1,llm),masseby(iip1,jjm,llm)
83      REAL teta(iip1,jjp1,llm)
84
85      INTEGER itau,irec,nrec,ldec,itauav
86      integer recmin,recmax,recstkmin,recstkmax
87      character*1 injecte,stoke
88      character*10 nom
89
90      integer jour0,isplit,nsplit_dyn,nsplit_phy,nsplit
91      logical debut,lectstart,rnpb
92
93      EXTERNAL inidissip,iniconst,inifilr
94      EXTERNAL inigeom
95      EXTERNAL exner_hyb,pression
96      EXTERNAL defrun_new
97
98      REAL time_0 ! temps de debut de simulation (% journee)
99
100      character*2 str2
101      INTEGER iq,j,i,l,nq,mode
102
103      real paprs(ngridmx,llm+1),pplay(ngridmx,llm),tempfi(ngridmx,llm)
104      real p(iip1,jjp1,llmp1),pk(iip1,jjp1,llm)
105
106      real lon_stat(nqmx),lat_stat(nqmx),zlon,zlat
107      integer i_stat(nqmx),j_stat(nqmx)
108      real intensite(nqmx)
109
110      integer iqmin,iqmax
111      integer npasinject
112      data npasinject/1/
113
114      integer histid, histvid, histaveid
115      real time_step, t_wrt, t_ops
116      integer itime_step
117      character*80 dynhist_file, dynhistave_file
118
119
120      logical avant
121
122      logical debutphy,redecoupage
123      DATA debutphy/.true./
124     
125      real unskap,pksurcp,smass(iip1,jjp1)
126      integer ig0,ii,jj,iii,ij
127
128     
129      REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
130      REAL ps(iip1,jjp1),pkf(ip1jmp1,llm)
131
132
133c-----------------------------------------------------------------------
134c   Initialisations:
135c   ----------------
136
137      descript = 'Run GCM LMDZ'
138
139c-----------------------------------------------------------------------
140
141
142      rad = 6400000
143      omeg = 7.272205e-05
144      g = 9.8
145      kappa = 0.285716
146      daysec = 86400
147      cpp = 1004.70885
148Cmaf a verifier....
149      preff = 101325.
150      pa= 50000.
151C
152
153      print*,'Nombre effectif de traceurs'
154      read(*,*) nq
155      write(*,*) nq
156      print*,'Splitting du pas de temps dynamique'
157      read(*,*) nsplit_dyn
158      write(*,*) nsplit_dyn
159      print*,'Splitting du pas de temps physique'
160      read(*,*) nsplit_phy
161      write(*,*) nsplit_phy
162      print*,'jour0, par rapport au debut du fichier fluxmass'
163      read(*,*) jour0
164      write(*,*) jour0
165      print*,'redecoupage (T) ou standard (F)'
166      read(*,'(l1)') redecoupage
167      write(*,*) redecoupage
168      print*,'Lecture fichier start (T)?'
169      read(*,'(l1)') lectstart
170      write(*,'(l1)') lectstart
171      print*,'Traceurs 1 et 2 = Rn222 et Pb'
172      read(*,'(l1)') rnpb
173      write(*,'(l1)') rnpb
174      print*,'Avant T ou arriere F?'
175      read(*,'(l1)') avant
176      write(*,'(l1)') avant
177      print*,'Index min et max des traceurs a initaliser'
178      read(*,*) iqmin,iqmax
179      write(*,*) iqmin,iqmax
180      print*,'localisation (lon, lat, intensite)'
181      if(iqmax.le.0) then
182         iqmin=iqmax+1
183      endif
184      do iq=iqmin,iqmax
185         read(*,*) lon_stat(iq),lat_stat(iq),intensite(iq)
186         write(*,*) lon_stat(iq),lat_stat(iq),intensite(iq)
187      enddo
188      print*,'Injection sur combien de pas de temps'
189      read(*,*) npasinject
190      write(*,*) npasinject
191      print*,'Records min et max pour le stokage et couche '
192     s    ,'verticale (specifique reseaux), def 0 0 1'
193      read(*,*) recstkmin,recstkmax,ldec
194      write(*,*) recstkmin,recstkmax,ldec
195
196c   redefinition du splitting
197      if(mod(nsplit_dyn,nsplit_phy).eq.0) then
198         nsplit=nsplit_phy
199      elseif (mod(nsplit_phy,nsplit_dyn).eq.0) then
200         nsplit=nsplit_dyn
201      else
202         stop'prendre un des deux nsplit doit diviser l autre'
203      endif
204      nsplit_dyn=nsplit_dyn/nsplit
205      nsplit_phy=nsplit_phy/nsplit
206      print*,'nsplit,nsplit_dyn,nsplit_phy'
207     s       ,nsplit,nsplit_dyn,nsplit_phy
208c   Lecture eventuelle
209      time_0=0.
210
211      if(nq.gt.nqmx) stop'surdimensioner nqmx'
212
213      print*,'av initial0'
214      call initial0(nq*ijp1llm,q)
215
216      print*,'av lectstart'
217      if(lectstart) then
218         CALL dynetat0("start.nc",nq,vcov,ucov,
219     .              teta,q,masse,ps,phis, time_0)
220         print*,'Lecture du start'
221      else
222         day_ini=0
223         anne_ini=0
224      endif
225      day_end=day_ini+nday
226
227
228      print*,'av defrun'
229c     open (99,file='run.def',status='old',form='formatted')
230      clefinit=.true.
231      CALL defrun_new( 99, clefinit, clesphy0 )
232c     close(99)
233
234      print*,'av iniconst'
235c   lecture du jour de demarrage
236c    premiere initialisation, eventuellement bidon
237      CALL iniconst
238      CALL inigeom
239
240C
241c  Lecture des entetes des fichiers de flux
242      if(redecoupage) then
243         call redecoupe(0,masse,pbaru,pbarv,w,teta,phi,
244     s     nrec,avant,airefi,
245     s     zmfu, zmfd, zen_u, zde_u,zen_d, zde_d, coefkz,
246     s     yu1,yv1,ftsol,pctsrf,
247     s     frac_impa,frac_nucl)
248      else
249         call lectflux(0,masse,pbaru,pbarv,w,teta,phi,
250     s     nrec,avant,airefi,
251     s     zmfu, zmfd, zen_u, zde_u,zen_d, zde_d, coefkz,
252     s     yu1,yv1,ftsol,pctsrf,
253     s     frac_impa,frac_nucl,phis)
254         call inigeom
255      endif
256
257c   controle sur les records
258      if(avant) then
259         recmin=2+(jour0-1)*iday_step
260         recmax=1+(jour0+nday-1)*iday_step
261      else
262         recmin=(jour0-nday)*iday_step+2
263         recmax=jour0*iday_step+1
264      endif
265      if(recmin.le.1.or.recmax.ge.nrec) then
266         print*,'Probleme de lecture '
267         print*,'recmin,recmax,nrec',recmin,recmax,nrec
268         stop
269      endif
270
271C on recalcule le nombre de pas de temps dans une journee
272c   istdyn est la frequence de stokage en "pas de temps dynamique" dans
273c   le modele on-line. dtvr est le pas de temps du meme modele en s.
274c iday-step est le nombre de lectures par jour (grand pas de temps).
275      iday_step=daysec/(istdyn*dtvr)
276C
277      CALL iniconst
278c  on calcule le pas de temps
279      dtvr    = daysec/FLOAT(iday_step)
280      dtphys=dtvr/(nsplit*nsplit_phy) ! a mettre apres iniconst!
281      print*,'Pas de temps physique ',dtphys
282
283
284C
285
286      call gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
287
288      print *,'Pas de temps dynamique', dtvr
289      print*,'Pas de temps physique ',dtphys
290
291c-----------------------------------------------------------------------
292C maf est ce utile ?
293
294      call suphec
295
296      CALL inifilr
297      print*,'Pas de temps physique ',dtphys
298      CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   ,
299     *                tetagdiv, tetagrot , tetatemp              )
300      print*,'Pas de temps physique ',dtphys
301
302c-----------------------------------------------------------------------
303c   temps de depart et de fin:
304c   --------------------------
305
306      itau = 0
307      iday = dayref
308      time = time_0
309         IF(time.GT.1.) THEN
310          time = time-1.
311          iday = iday+1
312         ENDIF
313      print*,'Pas de temps physique ',dtphys
314
315c=======================================================================
316c   Initialisation des fichiers de sortie
317c=======================================================================
318
319      print*,'Initialisation de wrgras'
320      file='sol'
321      call inigrads(1,iip1
322     s  ,rlonv,180./pi,-190.,190.,jjp1,rlatu,-90.,90.,180./pi
323     s  ,llm,presnivs,1.
324     s  ,dtphys,file,'gcmq2 ')
325      file='atm'
326      print*,'Pas de temps physique ',dtphys
327      call inigrads(2,iip1
328     s  ,rlonv,180./pi,-190.,190.,jjp1,rlatu,-90.,90.,180./pi
329     s  ,llm,presnivs,1.
330     s  ,dtphys,file,'gcmq2 ')
331
332      print*,'dtphys ap inigrads ',dtphys
333
334      print*,'Ecriture du fichier de conditions aux limites'
335      open (98,file='limit',form='unformatted')
336      write(98) float(im),float(jm),float(nq),float(nday)
337      write(98) (rlonv(i)*180./pi,i=1,iip1)
338      write(98) (rlatu(j)*180./pi,j=1,jjp1)
339      write(98) ((phis(i,j)/g,i=1,iip1),j=1,jjp1)
340      close(98)
341
342      CALL dynredem0("restart.nc",day_end,anne_ini,phis,nq)
343
344c-----------------------------------------------------------------------
345c   Debut de la boucle en temps:
346c   ----------------------------
347
348      do itau=1,nday*iday_step
349
350          print*,'ITAU ITAU ITAU =',itau
351         injecte='.'
352         stoke='.'
353
354c Gestion du temps
355         if (avant) then
356            dtdyn=dtvr/float(nsplit*nsplit_dyn)
357            irec=1+(jour0-1)*iday_step+itau
358         else
359            dtdyn=-dtvr/float(nsplit*nsplit_dyn)
360            irec=1+jour0*iday_step-itau+1
361         endif
362
363         debut=itau.eq.1
364
365c-----------------------------------------------------------------------
366c   Lecture des fichiers fluxmass et physique:
367c   -----------------------------------------------------
368
369         print*,'CCCCC Appel a redecoupe au pas itau=',itau
370      if(redecoupage) then
371         call redecoupe(irec,masse,pbaru,pbarv,w,teta,phi,
372     s     nrec,avant,airefi,
373     s     zmfu, zmfd, zen_u, zde_u,zen_d, zde_d, coefkz,
374     s     yu1,yv1,ftsol,pctsrf,
375     s     frac_impa,frac_nucl)
376      else
377         call lectflux(irec,masse,pbaru,pbarv,w,teta,phi,
378     s     nrec,avant,airefi,
379     s     zmfu, zmfd, zen_u, zde_u,zen_d, zde_d, coefkz,
380     s     yu1,yv1,ftsol,pctsrf,
381     s     frac_impa,frac_nucl,phis)
382      endif
383
384c    ...  ouverture du fichier de stockage netcdf ...
385C
386      IF (debut) then
387         dynhistave_file = 'dyn_hist_ave.nc'
388         day_ini=0
389         anne_ini=0
390         t_ops = periodav * daysec
391         t_wrt = periodav * daysec
392         dtav=dtvr/float(nsplit)
393         mode=1
394         CALL initdynav(dynhistave_file,day_ini,anne_ini,dtav,
395     .              t_ops, t_wrt, nq, histaveid)
396         pi=2.*asin(1.)
397         do iq=iqmin,iqmax
398            do l=1,llm
399               do j=1,jjp1
400                  do i=1,iip1
401                     q(i,j,l,iq)=0.
402                  enddo
403               enddo
404            enddo
405            zlon=lon_stat(iq)*pi/180.
406            zlat=lat_stat(iq)*pi/180.
407            CALL coordij(zlon,zlat,i_stat(iq),j_stat(iq))
408         enddo
409      ENDIF
410
411C calcul de la pression de surface.
412      do j=1,jjp1
413         do i=1,iip1
414            smass(i,j)=0.
415         enddo
416      enddo
417      do l=1,llm
418        do j=1,jjp1
419          do i=1,iip1
420             smass(i,j)=smass(i,j)+masse(i,j,l)
421          enddo
422       enddo
423      enddo
424      do j=1,jjp1
425         do i=1,iip1
426            ps (i,j)=smass(i,j)/aire(i,j)*g
427         end do
428      end do
429C
430C calcul de la pression et de pk (fonction d'Exner)
431C
432      CALL pression ( ip1jmp1, ap, bp, ps, p       )
433      CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
434c
435C
436      do jj=1,jjp1
437         do ii=1,iip1
438           phis(ii,jj)= phi(ii,jj,1)-teta(ii,jj,1)*
439     s                             (pks(ii,jj)-pk(ii,jj,1))
440         end do
441      end do
442
443c=======================================================================
444c   TERMES SOURCES OU PUITS
445c=======================================================================
446c   terme source reparti sur npasinject pas
447c   On injecte tous les traceurs de la meme facon
448      if (itau.le.npasinject) then
449        injecte='X'
450         do iq=iqmin,iqmax
451            print*,'INJECTION AU POINT ',lon_stat(iq),lat_stat(iq)
452     s      ,'     I,J=',i_stat(iq),j_stat(iq),jjp1-j_stat(iq)+1
453            q(i_stat(iq),j_stat(iq),1,iq)=
454     s      q(i_stat(iq),j_stat(iq),1,iq)
455     s      +intensite(iq)/masse(i_stat(iq),j_stat(iq),1)/npasinject
456         enddo
457         print*,'fin des injections'
458      endif
459c=======================================================================
460c   FIN DE LA PARTIE TERMES SOURCES
461c=======================================================================
462c-----------------------------------------------------------------------
463c   Advection
464c   ----------
465
466c=======================================================================
467         do isplit=1,nsplit
468c=======================================================================
469CC MAF dans le cas du offline on ne fait pas la difference masse et massem
470            do iii=1,nsplit_dyn
471               do iq=1,nq
472                  CALL vlsplt(q(:,:,:,iq),2.,masse,w,pbaru,pbarv,dtdyn)
473               enddo
474
475c   mise à jour du champd de masse tenant compte de l'advection sur
476c   le pas de temps considere.
477
478            enddo
479
480            do iq=1,nq
481               write(str2,'(i2.2)') iq
482               if(mod(itau,1).eq.0)
483     s         call diagadv(q(:,:,:,iq),masse,'Traceur '//str2)
484            enddo
485
486
487c-----------------------------------------------------------------------
488c   Physique  (lessivage dans phytrac maintenant) MAF
489c   ---------------------------------------------------------
490
491        if(physic) then
492         if(debut) then
493         latfi(1)=rlatu(1)*180./pi
494         lonfi(1)=0.
495         DO j=2,jjm
496            DO i=1,iim
497               latfi((j-2)*iim+1+i)= rlatu(j)*180./pi
498               lonfi((j-2)*iim+1+i)= rlonv(i)*180./pi
499            ENDDO
500         ENDDO
501         latfi(ngridmx)= rlatu(jjp1)*180./pi
502         lonfi(ngridmx)= 0.
503         endif
504
505c
506C Calcul de paprs et pplay et de temp
507c   -----------------------------------------------------------------
508c     .... paprs  definis aux (llm +1) interfaces des couches  ....
509c     .... pplay  definis aux (  llm )    milieux des couches  ....
510c   -----------------------------------------------------------------
511C
512      DO l = 1, llmp1
513        paprs( 1,l ) = p(1,1,l)
514        ig0 = 2
515          DO j = 2, jjm
516             DO i =1, iim
517              paprs( ig0,l ) = p(i,j,l)
518              ig0 = ig0 +1
519             ENDDO
520          ENDDO
521        paprs( ngridmx,l ) = p(1,jjp1,l)
522      ENDDO
523c
524      unskap   = 1./ kappa
525      DO l=1,llm
526
527         pksurcp     =  pk(1,1,l) / cpp
528         pplay(1,l)  =  preff * pksurcp ** unskap
529         tempfi(1,l)   =  teta(1,1,l) *  pksurcp
530         ig0         = 2
531
532         DO j = 2, jjm
533            DO i = 1, iim
534              pksurcp        = pk(i,j,l) / cpp
535              pplay(ig0,l)   = preff * pksurcp ** unskap
536              tempfi(ig0,l)    = teta(i,j,l)  * pksurcp
537              ig0            = ig0 + 1
538            ENDDO
539         ENDDO
540
541         pksurcp       = pk(1,jjp1,l) / cpp
542         pplay(ig0,l)  = preff * pksurcp ** unskap
543         tempfi (ig0,l)  = teta(1,jjp1,l)  * pksurcp
544      ENDDO
545c On passe le traceur et le geopot de surf  sur la grille physique
546             call gr_dyn_fi(llm*nq,iip1,jjp1,ngridmx,q,qfi)
547             call gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,pphis)
548
549             do iii=1,nsplit_phy
550C
551      print*,'dtphys avant phytrac ',dtphys
552             call phytrac(rnpb,
553     I                   ecritphy, debutphy,
554     I                   nq,
555     I                   ngridmx,llm,dtphys,
556     I                   tempfi,paprs,pplay,
557     I                   zmfu, zmfd, zen_u, zde_u, zen_d, zde_d,
558     I                   coefkz,yu1,yv1,ftsol,pctsrf,latfi,
559     I                   frac_impa, frac_nucl,
560     I                   lonfi,presnivs,airefi,pphis,
561     O                   qfi)
562C
563             debutphy= .FALSE.
564C
565             enddo ! nsplit_phy
566C
567c   On passe le traceur sur la grille dynamique.
568             call gr_fi_dyn(llm*nq,ngridmx,iip1,jjp1,qfi,q)
569
570             endif
571
572             itauav=(itau-1)*nsplit+isplit
573             CALL writedynav(histaveid, nq, itauav,vcov ,
574     ,                   ucov,teta,pk,phi,q,masse,ps,phis)
575
576         enddo ! isplit=1,nsplit
577
578         call histsync
579
580c-----------------------------------------------------------------------
581c    Calcul de ucov et vcov pour les sorties
582c   -------------------------------------
583
584         call massbar(masse, massebx, masseby )
585         do l=1,llm
586            do j=1,jjp1
587               do i=1,iip1
588               ucov(i,j,l)=pbaru(i,j,l)/massebx(i,j,l)*cu(i,j)*cu(i,j)
589     s         /istdyn
590               enddo
591            enddo
592            do j=1,jjm
593               do i=1,iip1
594               vcov(i,j,l)=pbarv(i,j,l)/masseby(i,j,l)*cv(i,j)*cv(i,j)
595     s         /istdyn
596               enddo
597            enddo
598         enddo
599
600         iday= dayref+itau/iday_step
601         time=
602     s   float(itau-(iday-dayref)*iday_step)/iday_step+time_0
603         IF(time.GT.1.) THEN
604            time = time-1.
605            iday = iday+1
606         ENDIF
607
608
609
610C
611c=======================================================================
612         if(irec.ge.recstkmin.and.irec.le.recstkmax) then
613            if (mod(itau,4).eq.0) then
614            stoke='X'
615               do iq=1,nq
616                  write(str2,'(i2.2)') iq
617                  nom='q'//str2
618                  call wrgrads(1,1,q(:,:,:,iq),nom,nom)
619                  call wrgrads(2,1,q(:,:,ldec,iq),nom,nom)
620               enddo
621            endif
622         endif
623
624         print*,'Record=',irec,'   stoke=',stoke,
625     s      '      injecte=',injecte
626
627      enddo ! itau
628
629      CALL dynredem1("restart.nc",0.0,
630     .          vcov,ucov,teta,q,nq,masse,ps)
631
632c=======================================================================
633C fermeture des fichiers netcdf
634          CALL histclo
635C
636 300  FORMAT('1',/,15x,'run du pas',i7,2x,'au pas',i7,2x,
637     s 'c"est a dire dujour',i7,3x,'au jour',i7,/,/)
638
639      end
640C----------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.