source: LMDZ.3.3/branches/rel-1-0-patch/libf/dyn3d/offlinenc.F @ 1176

Last change on this file since 1176 was 253, checked in by (none), 23 years ago

This commit was manufactured by cvs2svn to create branch
'rel-1-0-patch'.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 27.5 KB
Line 
1c
2c $Header$
3c
4      PROGRAM offlinenc
5      USE ioipsl
6      IMPLICIT NONE
7
8#include "dimensions.h"
9#include "paramet.h"
10#include "comconst.h"
11#include "comdissnew.h"
12#include "comvert.h"
13#include "comgeom2.h"
14#include "logic.h"
15#include "temps.h"
16#include "control.h"
17#include "ener.h"
18#include "description.h"
19#include "serre.h"
20#include "tracstoke.h"
21      integer ngridmx,nbsrf
22      PARAMETER (ngridmx=iim*(jjm-1)+2,nbsrf=4)
23c-----------------------------------------------------------------------
24
25      INTEGER*4  iday ! jour julien
26      REAL       time ! Heure de la journee en fraction d'1 jour
27      REAL dtdyn,dtav
28      INTEGER iday_step
29
30      INTEGER         longcles
31      PARAMETER     ( longcles = 20 )
32      REAL  clesphy0( longcles )
33      logical clefinit
34
35c   variables dynamiques
36      REAL vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants
37      REAL q(iip1,jjp1,llm,nqmx)               ! champs advectes
38        real qmoy(iip1,jjp1,nqmx)
39      REAL pks(iip1,jjp1)                      ! exner (f pour filtre)
40      REAL phis(iip1,jjp1)                     ! geopotentiel au sol
41      REAL phi(iip1,jjp1,llm)                  ! geopotentiel
42      REAL w(iip1,jjp1,llm)                    ! vitesse verticale
43C
44C   variables physiques
45      real zmfd(ngridmx,llm),zde_d(ngridmx,llm),zen_d(ngridmx,llm)
46      REAL qfi(ngridmx,llm,nqmx)               ! champs advectes
47      real zmfu(ngridmx,llm),zde_u(ngridmx,llm),zen_u(ngridmx,llm)
48        real t(ngridmx,llm)
49      real coefkz(ngridmx,llm)
50      real frac_impa(ngridmx,llm),frac_nucl(ngridmx,llm)
51      real pphis (ngridmx)
52      REAL airefi(ngridmx)
53      REAL latfi(ngridmx),lonfi(ngridmx) ! latitude et long pour chaque point
54      real yu1(ngridmx), yv1(ngridmx)
55      real ftsol(ngridmx,nbsrf),pctsrf(ngridmx,nbsrf)
56C
57      character*10 file
58c variables dynamiques intermediaire pour le transport
59      REAL pbaru(iip1,jjp1,llm),pbarv(iip1,jjm,llm) !flux de masse
60      REAL masse(iip1,jjp1,llm)
61      REAL massebx(iip1,jjp1,llm),masseby(iip1,jjm,llm)
62      REAL teta(iip1,jjp1,llm)
63      INTEGER itau,irec,nrec,ldec,itauav
64      integer recmin,recmax,recstkmin,recstkmax
65      character*1 injecte,stoke
66      character*10 nom
67
68        real tinj2,tinj1
69        integer iiinj
70      integer jour0,isplit,nsplit_dyn,nsplit_phy,nsplit
71      logical debut,lectstart,rnpb,lafin
72
73      EXTERNAL inidissip,iniconst,inifilr
74      EXTERNAL inigeom
75      EXTERNAL exner_hyb,pression
76      EXTERNAL defrun_new
77
78      REAL time_0 ! temps de debut de simulation (% journee)
79
80      character*2 str2
81      INTEGER iq,j,i,l,nq,mode
82      real paprs(ngridmx,llm+1),pplay(ngridmx,llm),tempfi(ngridmx,llm)
83      real p(iip1,jjp1,llmp1),pk(iip1,jjp1,llm)
84      real lon_stat(nqmx),lat_stat(nqmx),zlon,zlat
85      integer i_stat(nqmx),j_stat(nqmx)
86      real intensite(nqmx)
87      integer iqmin,iqmax
88      integer npasinject
89      data npasinject/1/
90      integer histid, histvid, histaveid
91      real time_step, t_wrt, t_ops
92      integer itime_step
93      character*80 dynhist_file, dynhistave_file
94
95      logical avant
96      logical debutphy,redecoupage
97      DATA debutphy/.true./
98
99      real unskap,pksurcp,smass(iip1,jjp1)
100      integer ig0,ii,jj,iii,ij
101
102      REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
103      REAL ps(iip1,jjp1),pkf(ip1jmp1,llm)
104
105        integer ntau, iinjec2,iij
106        real zttt
107c-----------------------------------------------------------------------
108c   Initialisations:
109c   ----------------
110      descript = 'Run GCM LMDZ'
111        tinj1=0.
112        tinj2=0.
113      rad = 6400000
114      omeg = 7.272205e-05
115      g = 9.8
116      kappa = 0.285716
117      daysec = 86400
118      cpp = 1004.70885
119
120Cmaf a verifier....
121      preff = 101325.
122      pa= 50000.
123C
124      print*,'PARAMETRES A ENTRER DANS offline.def'
125      print*,'Nombre effectif de traceurs'
126      read(*,*) nq
127      write(*,*) nq
128      print*,'Splitting du pas de temps dynamique'
129      read(*,*) nsplit_dyn
130      write(*,*) nsplit_dyn
131      print*,'Splitting du pas de temps physique'
132      read(*,*) nsplit_phy
133      write(*,*) nsplit_phy
134      print*,'jour0, par rapport au debut du fichier fluxmass'
135      read(*,*) jour0
136      write(*,*) jour0
137      print*,'redecoupage (T) ou standard (F)'
138      read(*,'(l1)') redecoupage
139      write(*,*) redecoupage
140      print*,'Lecture fichier start (T)?'
141      read(*,'(l1)') lectstart
142      write(*,'(l1)') lectstart
143      print*,'Traceurs 1 et 2 = Rn222 et Pb'
144      read(*,'(l1)') rnpb
145      write(*,'(l1)') rnpb
146      print*,'Avant T ou arriere F?'
147      read(*,'(l1)') avant
148      write(*,'(l1)') avant
149      print*,'Index min et max des traceurs a initaliser'
150      read(*,*) iqmin,iqmax
151      write(*,*) iqmin,iqmax
152      print*,'localisation (lon, lat, intensite)'
153
154      if(iqmax.le.0) then
155         iqmin=iqmax+1
156      endif
157
158      do iq=iqmin,iqmax
159         read(*,*) lon_stat(iq),lat_stat(iq),intensite(iq)
160         write(*,*) lon_stat(iq),lat_stat(iq),intensite(iq)
161      enddo
162
163      print*,'Records min et max pour le stokage et couche '
164     s    ,'verticale (specifique reseaux), def 0 0 1'
165      read(*,*) recstkmin,recstkmax,ldec
166      write(*,*) recstkmin,recstkmax,ldec
167
168c   redefinition du splitting
169
170      if(mod(nsplit_dyn,nsplit_phy).eq.0) then
171         nsplit=nsplit_phy
172      elseif (mod(nsplit_phy,nsplit_dyn).eq.0) then
173         nsplit=nsplit_dyn
174      else
175         stop'prendre un des deux nsplit doit diviser l autre'
176      endif
177
178      nsplit_dyn=nsplit_dyn/nsplit
179      nsplit_phy=nsplit_phy/nsplit
180
181      print*,'nsplit,nsplit_dyn,nsplit_phy'
182     s       ,nsplit,nsplit_dyn,nsplit_phy
183
184c   Lecture eventuelle
185      time_0=0.
186      if(nq.gt.nqmx) stop'surdimensioner nqmx'
187      print*,'av initial0'
188      call initial0(nq*ijp1llm,q)
189      print*,'av lectstart'
190
191      if(lectstart) then
192         CALL dynetat0("start.nc",nq,vcov,ucov,
193     .              teta,q,masse,ps,phis, time_0)
194         print*,'Lecture du start'
195c  On zappe le radon et le plomb.
196         q(:,:,:,1)=q(:,:,:,3)
197         q(:,:,:,2)=q(:,:,:,4)
198c Initialisation d'un traceur ` 1 pour tester l'impact du lessivage.
199         q(:,:,:,3)=1.
200      else
201         day_ini=0
202         anne_ini=0
203      endif
204
205      day_end=day_ini+nday
206C      print*,'av defrun'
207c     open (99,file='run.def',status='old',form='formatted')
208      clefinit=.true.
209      CALL defrun_new( 99, clefinit, clesphy0 )
210c     close(99)
211C      print*,'av iniconst'
212c   lecture du jour de demarrage
213c    premiere cnitialisation, eventuellement bidon
214      cALL iniconst
215      CALL inigeom
216C
217      print*,'ENTREE DANS redecoupenc ou lectfluxnc
218     s        pour irec=0'
219
220c  Lecture des entetes des fichiers de flux
221
222      if(redecoupage) then
223         call redecoupenc(0,masse,pbaru,pbarv,w,teta,phi,
224     s     nrec,avant,airefi,pphis,
225     s     t,zmfu, zmfd, zen_u, zde_u,zen_d, zde_d, coefkz,
226     s     yu1,yv1,ftsol,pctsrf,
227     s     frac_impa,frac_nucl,phis)
228      else
229         call lectfluxnc(0,masse,pbaru,pbarv,w,teta,phi,
230     s     nrec,avant,airefi,pphis,
231     s     t,zmfu, zmfd, zen_u, zde_u,zen_d, zde_d, coefkz,
232     s     yu1,yv1,ftsol,pctsrf,
233     s     frac_impa,frac_nucl,phis)
234         call inigeom
235      endif
236
237c   controle sur les records
238
239      if(avant) then
240         recmin=2+(jour0-1)*iday_step
241         recmax=1+(jour0+nday-1)*iday_step
242      else
243         recmin=(jour0-nday)*iday_step+2
244         recmax=jour0*iday_step+1
245      endif
246
247      if(recmin.le.1.or.recmax.ge.nrec) then
248         print*,'Probleme de lecture '
249         print*,'recmin,recmax,nrec',recmin,recmax,nrec
250         stop
251      endif
252
253C on recalcule le nombre de pas de temps dans une journee
254c   istdyn est la frequence de stokage en "pas de temps dynamique" dans
255c   le modele on-line. dtvr est le pas de temps du meme modele en s.
256c iday-step est le nombre de lectures par jour (grand pas de temps).
257
258      iday_step=daysec/(istdyn*dtvr)
259C
260      CALL iniconst
261c  on calcule le pas de temps
262      dtvr    = daysec/FLOAT(iday_step)
263      dtphys=dtvr/(nsplit*nsplit_phy) ! a mettre apres iniconst!
264      print*,'Pas de temps physique ',dtphys
265C
266      call gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
267      print *,'Pas de temps dynamique', dtvr
268      print*,'Pas de temps physique ',dtphys
269
270c-----------------------------------------------------------------------
271
272C maf est ce utile ?
273      call suphec
274      CALL inifilr
275      print*,'Pas de temps physique ',dtphys
276      CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   ,
277     *                tetagdiv, tetagrot , tetatemp              )
278      print*,'Pas de temps physique ',dtphys
279
280c-----------------------------------------------------------------------
281c   temps de depart et de fin:
282c   --------------------------
283      itau = 0
284      iday = dayref
285      time = time_0
286         IF(time.GT.1.) THEN
287          time = time-1.
288          iday = iday+1
289         ENDIF
290      print*,'Pas de temps physique ',dtphys
291c=======================================================================
292c   Initialisation des fichiers de sortie
293c=======================================================================
294      print*,'Initialisation de wrgras'
295        do i=1,iip1
296        print*,'rlonu(',i,')=',rlonu(i)
297        enddo
298        do i=1,jjp1
299        print*,'rlatu(',i,')=',rlatu(i)
300        enddo
301        do i=1,iip1
302        print*,'rlonv(',i,')=',rlonv(i)
303        enddo
304        do i=1,jjm
305        print*,'rlatv(',i,')=',rlatv(i)
306        enddo
307
308      file='sol'
309      call inigrads(1,iip1
310     s  ,rlonv,180./pi,-180.,180.,jjp1,rlatu,-90.,90.,180./pi
311     s  ,llm,presnivs,1.
312     s  ,dtphys,file,'gcmq2 ')
313
314      file='atm'
315      print*,'Pas de temps physique ',dtphys
316      call inigrads(2,iip1
317     s  ,rlonv,180./pi,-180.,180.,jjp1,rlatu,-90.,90.,180./pi
318     s  ,llm,presnivs,1.
319     s  ,dtphys,file,'gcmq2 ')
320
321      file='pbu'
322      call inigrads(3,iip1
323     s  ,rlonu,180./pi,-180.,180.,jjp1,rlatu,-90.,90.,180./pi
324     s  ,llm,presnivs,1.
325     s  ,dtphys,file,'gcmq2 ')
326
327      file='pbv'
328      call inigrads(4,iip1
329     s  ,rlonv,180./pi,-180.,180.,jjm,rlatv,-90.,90.,180./pi
330     s  ,llm,presnivs,1.
331     s  ,dtphys,file,'gcmq2 ')
332
333      file='qadv'
334      call inigrads(5,iip1
335     s  ,rlonv,180./pi,-180.,180.,jjp1,rlatu,-90.,90.,180./pi
336     s  ,llm,presnivs,1.
337     s  ,dtphys,file,'gcmq2 ')
338
339      print*,'dtphys ap inigrads ',dtphys
340
341      print*,'Ecriture du fichier de conditions aux limites'
342      open (98,file='limit',form='unformatted')
343      write(98) float(im),float(jm),float(nq),float(nday)
344      write(98) (rlonv(i)*180./pi,i=1,iip1)
345      write(98) (rlatu(j)*180./pi,j=1,jjp1)
346      write(98) ((phis(i,j)/g,i=1,iip1),j=1,jjp1)
347      close(98)
348
349      CALL dynredem0("restart.nc",day_end,anne_ini,phis,nq)
350
351c-----------------------------------------------------------------------
352c   Debut de la boucle en temps:
353c   ----------------------------
354
355      do itau=1,nday*iday_step
356          print*,'ITAU ITAU ITAU =',itau
357         injecte='.'
358         stoke='.'
359c Gestion du temps
360         if (avant) then
361            dtdyn=dtvr/float(nsplit*nsplit_dyn)
362c            irec=1+(jour0-1)*iday_step+itau
363             irec=(jour0-1)*iday_step+itau
364         else
365            dtdyn=-dtvr/float(nsplit*nsplit_dyn)
366            irec=1+jour0*iday_step-itau+1
367         endif
368
369         debut=itau.eq.1
370
371c-----------------------------------------------------------------------
372c   Lecture des fichiers fluxmass et physique:
373c   -----------------------------------------------------
374         print*,'CCCCC Appel a REDECOUPENC ou LECTFLUXNC
375     s           au pas itau=',itau
376      if(redecoupage) then
377         call redecoupenc(irec,masse,pbaru,pbarv,w,teta,phi,
378     s     nrec,avant,airefi,pphis,
379     s     t,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      else
383         call lectfluxnc(irec,masse,pbaru,pbarv,w,teta,phi,
384     s     nrec,avant,airefi,pphis,
385     s     t,zmfu, zmfd, zen_u, zde_u,zen_d, zde_d, coefkz,
386     s     yu1,yv1,ftsol,pctsrf,
387     s     frac_impa,frac_nucl,phis)
388      endif
389
390
391      print*,'TESTPHYS: ON PREND LA PUIS ',1./float(nsplit*nsplit_phy)
392     s   ,'  DES FRAC A IT=',itau
393      frac_impa(:,:)=frac_impa(:,:)**(1./float(nsplit*nsplit_phy))
394      frac_nucl(:,:)=frac_nucl(:,:)**(1./float(nsplit*nsplit_phy))
395c    ...  ouverture du fichier de stockage netcdf ...
396C
397      IF (debut) then
398         dynhistave_file = 'histmoy.nc'
399         day_ini=0
400         anne_ini=0
401         t_ops =(1./48.)*daysec
402         t_wrt =(1./48.)*daysec
403         dtav=dtvr/float(nsplit)
404         mode=1
405
406c        CALL initdynav(dynhistave_file,day_ini,anne_ini,dtav,
407c    .              t_ops, t_wrt, nq,mode, histaveid)
408
409         pi=2.*asin(1.)
410         do iq=iqmin,iqmax
411            do l=1,llm
412               do j=1,jjp1
413                  do i=1,iip1
414                     q(i,j,l,iq)=0.
415                  enddo
416               enddo
417            enddo
418
419            zlon=lon_stat(iq)*pi/180.
420            zlat=lat_stat(iq)*pi/180.
421
422            CALL coordij(zlon,zlat,i_stat(iq),j_stat(iq))
423
424         enddo
425      ENDIF
426
427C calcul de la pression de surface.
428
429      do j=1,jjp1
430         do i=1,iip1
431            smass(i,j)=0.
432         enddo
433      enddo
434
435      do l=1,llm
436        do j=1,jjp1
437          do i=1,iip1
438             smass(i,j)=smass(i,j)+masse(i,j,l)
439          enddo
440       enddo
441      enddo
442
443      do j=1,jjp1
444         do i=1,iip1
445            ps (i,j)=smass(i,j)/aire(i,j)*g
446         end do
447      end do
448C
449C calcul de la pression et de pk (fonction d'Exner)
450C
451      CALL pression ( ip1jmp1, ap, bp, ps, p       )
452      CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
453C
454      do jj=1,jjp1
455         do ii=1,iip1
456           phis(ii,jj)= phi(ii,jj,1)-teta(ii,jj,1)*
457     s                             (pks(ii,jj)-pk(ii,jj,1))
458         end do
459      end do
460c=======================================================================
461c   TERMES SOURCES OU PUITS
462         do isplit=1,nsplit
463c=======================================================================
464c injection par Fred
465c       Print*,Verifier en fct de iapp_tracvl la boucle inj
466       goto 333
467c Inj pour le stokage ttes les 0.5h
468        if (itau.gt.32.and.itau.le.56)then
469        print*,'Masse de la maille inj :',
470     s          masse(i_stat(1),j_stat(1),1)
471        do iiinj=1,20
472        q(i_stat(1),j_stat(1),1,1)=
473     s      q(i_stat(1),j_stat(1),1,1)
474     s  +intensite(1)/masse(i_stat(1),j_stat(1),1)/
475     s   (float(nsplit)*24.*20.)
476
477        tinj1=tinj1+intensite(1)/(float(nsplit)*24.*20.)
478        enddo
479
480           print*,'Apres injection a itau ',itau
481
482            do iq=1,nq
483               write(str2,'(i2.2)') iq
484              call diagadv(q(:,:,:,iq),masse,'Traceur '//str2)
485            enddo
486
487        print*,'QT TOTALE INJECTEE TR1',tinj1
488        print*,'fin des injections'
489
490        endif
491
492c Inj pour le stokage ttes les 1h
493
494        if (itau.gt.16.and.itau.le.28)then
495        print*,'Masse de la maille inj :',
496     s          masse(i_stat(1),j_stat(1),1)
497        do iiinj=1,20
498        q(i_stat(1),j_stat(1),1,1)=
499     s      q(i_stat(1),j_stat(1),1,1)
500     s  +intensite(1)/masse(i_stat(1),j_stat(1),1)/
501     s   (float(nsplit)*12.*20.)
502
503        tinj1=tinj1+intensite(1)/(float(nsplit)*12.*20.)
504
505        enddo
506
507           print*,'Apres injection a itau ',itau
508
509            do iq=1,nq
510               write(str2,'(i2.2)') iq
511              call diagadv(q(:,:,:,iq),masse,'Traceur '//str2)
512            enddo
513
514        print*,'QT TOTALE INJECTEE TR1',tinj1
515        print*,'fin des injections'
516
517        endif
518
519c Inj pour le stokage ttes les 1.5h
520
521        if (itau.eq.11.and.isplit.gt.2)then
522        do iiinj=1,20
523        q(i_stat(1),j_stat(1),1,1)=
524     s      q(i_stat(1),j_stat(1),1,1)
525     s  +intensite(1)/masse(i_stat(1),j_stat(1),1)/
526     s   (24.*20.)
527        enddo
528        endif
529
530        if (itau.ge.12.and.itau.le.18)then
531        do iiinj=1,20
532        q(i_stat(1),j_stat(1),1,1)=
533     s      q(i_stat(1),j_stat(1),1,1)
534     s  +intensite(1)/masse(i_stat(1),j_stat(1),1)/
535     s   (24.*20.)
536        enddo
537        endif
538
539        if (itau.eq.19.and.isplit.le.2)then
540        do iiinj=1,20
541        q(i_stat(1),j_stat(1),1,1)=
542     s      q(i_stat(1),j_stat(1),1,1)
543     s  +intensite(1)/masse(i_stat(1),j_stat(1),1)/
544     s   (24.*20.)
545        enddo
546
547        endif
548
549c Inj pour le stokage ttes les 4h
550
551        if (itau.ge.5.and.itau.le.7)then
552        do iiinj=1,20
553        q(i_stat(1),j_stat(1),1,1)=
554     s      q(i_stat(1),j_stat(1),1,1)
555     s  +intensite(1)/masse(i_stat(1),j_stat(1),1)/
556     s   (24.*20.)
557        enddo
558        endif
559
560c Inj pour le stokage ttes les 3h
561
562        if (itau.eq.6.and.isplit.gt.2)then
563        do iiinj=1,20
564        q(i_stat(1),j_stat(1),1,1)=
565     s      q(i_stat(1),j_stat(1),1,1)
566     s  +intensite(1)/masse(i_stat(1),j_stat(1),1)/
567     s   (24.*20.)
568        enddo
569        endif
570
571        if (itau.ge.7.and.itau.le.9)then
572        do iiinj=1,20
573        q(i_stat(1),j_stat(1),1,1)=
574     s      q(i_stat(1),j_stat(1),1,1)
575     s  +intensite(1)/masse(i_stat(1),j_stat(1),1)/
576     s   (24.*20.)
577        enddo
578        endif
579
580        if (itau.eq.10.and.isplit.le.2)then
581        do iiinj=1,20
582        q(i_stat(1),j_stat(1),1,1)=
583     s      q(i_stat(1),j_stat(1),1,1)
584     s  +intensite(1)/masse(i_stat(1),j_stat(1),1)/
585     s   (24.*20.)
586        enddo
587        endif
588
589c  Inj pour le stokage ttes les 6h
590
591        if (itau.eq.3.and.isplit.gt.8)then
592        do iiinj=1,20
593        q(i_stat(1),j_stat(1),1,1)=
594     s      q(i_stat(1),j_stat(1),1,1)
595     s  +intensite(1)/masse(i_stat(1),j_stat(1),1)/
596     s   (24.*20.)
597        enddo
598        endif
599
600        if (itau.eq.4)then
601        do iiinj=1,20
602        q(i_stat(1),j_stat(1),1,1)=
603     s      q(i_stat(1),j_stat(1),1,1)
604     s  +intensite(1)/masse(i_stat(1),j_stat(1),1)/
605     s   (24.*20.)
606        enddo
607        endif
608
609        if (itau.eq.5.and.isplit.le.8)then
610        do iiinj=1,20
611        q(i_stat(1),j_stat(1),1,1)=
612     s      q(i_stat(1),j_stat(1),1,1)
613     s  +intensite(1)/masse(i_stat(1),j_stat(1),1)/
614     s   (24.*20.)
615        enddo
616        endif
617
618c Inj sur les 4 1ere couches pour le stokage ttes les 3h
619
620        if (itau.eq.6.and.isplit.gt.2)then
621        do l=1,4
622        q(i_stat(1),j_stat(1),l,1)=
623     s      q(i_stat(1),j_stat(1),l,1)
624     s  +intensite(1)/masse(i_stat(1),j_stat(1),l)/
625     s   (24.*4.)
626        enddo
627        endif
628
629        if (itau.ge.7.and.itau.le.9)then
630        do l=1,4
631        q(i_stat(1),j_stat(1),l,1)=
632     s      q(i_stat(1),j_stat(1),l,1)
633     s  +intensite(1)/masse(i_stat(1),j_stat(1),l)/
634     s   (24.*4.)
635        enddo
636        endif
637
638        if (itau.eq.10.and.isplit.le.2)then
639        do l=1,4
640        q(i_stat(1),j_stat(1),l,1)=
641     s      q(i_stat(1),j_stat(1),l,1)
642     s  +intensite(1)/masse(i_stat(1),j_stat(1),l)/
643     s   (24.*4.)
644        enddo
645        endif
646
647c Inj pour le stokage ttes les 12h
648
649        if (itau.eq.2.and.isplit.gt.8)then
650        do iiinj=1,20
651        q(i_stat(1),j_stat(1),1,1)=
652     s      q(i_stat(1),j_stat(1),1,1)
653     s  +intensite(1)/masse(i_stat(1),j_stat(1),1)/
654     s   (24.*20.)
655        enddo
656        endif
657
658        if (itau.eq.3.and.isplit.le.8)then
659        do iiinj=1,20
660        q(i_stat(1),j_stat(1),1,1)=
661     s      q(i_stat(1),j_stat(1),1,1)
662     s  +intensite(1)/masse(i_stat(1),j_stat(1),1)/
663     s   (24.*20.)
664        enddo
665        endif
666
667c Inj pour le stokage ttes les 24h
668
669        if (itau.eq.1.and.isplit.gt.32)then
670        q(i_stat(1),j_stat(1),1,1)=
671     s      q(i_stat(1),j_stat(1),1,1)
672     s  +intensite(1)/masse(i_stat(1),j_stat(1),1)/24.
673        endif
674
675        if (itau.eq.2.and.isplit.le.8)then
676        q(i_stat(1),j_stat(1),1,1)=
677     s      q(i_stat(1),j_stat(1),1,1)
678     s  +intensite(1)/masse(i_stat(1),j_stat(1),1)/24.
679        endif     
680
681c Inj pour le stokage ttes les 3h
682
683        if (itau.eq.6.and.isplit.gt.2)then
684        do iiinj=1,20
685        q(i_stat(1),j_stat(1),1,1)=
686     s      q(i_stat(1),j_stat(1),1,1)
687     s  +intensite(1)/masse(i_stat(1),j_stat(1),1)/
688     s   (24.*20.)
689        enddo
690        endif
691
692        if (itau.ge.7.and.itau.le.9)then
693        do iiinj=1,20
694        q(i_stat(1),j_stat(1),1,1)=
695     s      q(i_stat(1),j_stat(1),1,1)
696     s  +intensite(1)/masse(i_stat(1),j_stat(1),1)/
697     s   (24.*20.)
698        enddo
699        endif
700
701        if (itau.eq.10.and.isplit.le.2)then
702        do iiinj=1,20
703        q(i_stat(1),j_stat(1),1,1)=
704     s      q(i_stat(1),j_stat(1),1,1)
705     s  +intensite(1)/masse(i_stat(1),j_stat(1),1)/
706     s   (24.*20.)
707        enddo
708        endif
709
710333     continue
711
712        if (itau.eq.6.and.isplit.gt.2)then
713
714        do l=1,2
715        q(i_stat(1),j_stat(1),l,1)=
716     s      q(i_stat(1),j_stat(1),l,1)
717     s  +intensite(1)/masse(i_stat(1),j_stat(1),l)/
718     s   (24.*2.)
719        enddo
720        endif
721
722        if (itau.ge.7.and.itau.le.9)then
723        do l=1,2
724        q(i_stat(1),j_stat(1),l,1)=
725     s      q(i_stat(1),j_stat(1),l,1)
726     s  +intensite(1)/masse(i_stat(1),j_stat(1),l)/
727     s   (24.*2.)
728        enddo
729        endif
730
731        if (itau.eq.10.and.isplit.le.2)then
732        do l=1,2
733        q(i_stat(1),j_stat(1),l,1)=
734     s      q(i_stat(1),j_stat(1),l,1)
735     s  +intensite(1)/masse(i_stat(1),j_stat(1),l)/
736     s   (24.*2.)
737        enddo
738        endif
739
740c=======================================================================
741c   FIN DE LA PARTIE TERMES SOURCES
742c=======================================================================
743c-----------------------------------------------------------------------
744c   Advection
745c   ----------
746        print*,'ENTREE DANS VLSPLT'
747c=======================================================================
748c         do isplit=1,nsplit
749c=======================================================================
750CC MAF dans le cas du offline on ne fait pas la difference masse et massem
751            do iii=1,nsplit_dyn
752                  nom='pbaru'
753                  call wrgrads(3,llm,pbaru(:,:,1),nom,nom)
754                  nom='pbarv'
755                  call wrgrads(4,llm,pbarv(:,:,1),nom,nom)
756                  nom='masse'
757                  call wrgrads(3,llm,masse(:,:,1),nom,nom)
758                  nom='w'
759                  call wrgrads(3,llm,w(:,:,1),nom,nom)
760           print*,'Entree dans vlsplt a itau et au temps',itau,time
761           print*,'Pas advect :',dtdyn
762           print*,'Avant vlsplt'
763            do iq=1,nq
764               write(str2,'(i2.2)') iq
765              call diagadv(q(:,:,:,iq),masse,'Traceur '//str2)
766            enddo
767               do iq=1,nq
768ctest abder   CALL vlsplt(q(1,1,1,iq),2.,masse,w,pbaru,pbarv,dtdyn)
769              CALL vlsplt(q(:,:,:,iq),2.,masse,w,pbaru,pbarv,dtdyn)
770               enddo
771
772               do iq=iqmin,iqmax
773                  write(str2,'(i2.2)') iq
774                  nom='q'//str2
775                  call wrgrads(5,1,q(:,:,1,iq),nom,nom)
776               enddo
777            print*,'Apres vlsplt'
778            do iq=1,nq
779               write(str2,'(i2.2)') iq
780              call diagadv(q(:,:,:,iq),masse,'Traceur '//str2)
781            enddo
782
783c   mise a jour du champ de masse tenant compte de l'advection sur
784c   le pas de temps considere.
785
786            enddo
787
788c-----------------------------------------------------------------------
789c   Physique  (lessivage dans phytrac maintenant) MAF
790c   ---------------------------------------------------------
791        if(physic) then
792         if(debut) then
793         latfi(1)=rlatu(1)*180./pi
794         lonfi(1)=0.
795         DO j=2,jjm
796            DO i=1,iim
797               latfi((j-2)*iim+1+i)= rlatu(j)*180./pi
798               lonfi((j-2)*iim+1+i)= rlonv(i)*180./pi
799            ENDDO
800         ENDDO
801         latfi(ngridmx)= rlatu(jjp1)*180./pi
802         lonfi(ngridmx)= 0.
803         endif
804c
805C Calcul de paprs et pplay et de temp
806c   -----------------------------------------------------------------
807c     .... paprs  definis aux (llm +1) interfaces des couches  ....
808c     .... pplay  definis aux (  llm )    milieux des couches  ....
809c   -----------------------------------------------------------------
810C
811      DO l = 1, llmp1
812        paprs( 1,l ) = p(1,1,l)
813        ig0 = 2
814          DO j = 2, jjm
815             DO i =1, iim
816              paprs( ig0,l ) = p(i,j,l)
817              ig0 = ig0 +1
818             ENDDO
819          ENDDO
820        paprs( ngridmx,l ) = p(1,jjp1,l)
821      ENDDO
822c
823      unskap   = 1./ kappa
824      DO l=1,llm
825         pksurcp     =  pk(1,1,l) / cpp
826         pplay(1,l)  =  preff * pksurcp ** unskap
827         tempfi(1,l)   =  teta(1,1,l) *  pksurcp
828         ig0         = 2
829         DO j = 2, jjm
830            DO i = 1, iim
831              pksurcp        = pk(i,j,l) / cpp
832              pplay(ig0,l)   = preff * pksurcp ** unskap
833              tempfi(ig0,l)    = teta(i,j,l)  * pksurcp
834              ig0            = ig0 + 1
835            ENDDO
836         ENDDO
837         pksurcp       = pk(1,jjp1,l) / cpp
838         pplay(ig0,l)  = preff * pksurcp ** unskap
839         tempfi (ig0,l)  = teta(1,jjp1,l)  * pksurcp
840      ENDDO
841
842c On passe le traceur et le geopot de surf  sur la grille physique
843             call gr_dyn_fi(llm*nq,iip1,jjp1,ngridmx,q,qfi)
844c             call gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,pphis)
845             do iii=1,nsplit_phy
846C
847      lafin=.false. ! en attendant mieux.
848      print*,'dtphys avant phytrac ',dtphys
849      print*,'TESTPHYS: APPEL A PHYTRAC IT=',itau
850
851             call phytrac(rnpb,
852     I                   debutphy, lafin,
853     I                   nq,
854     I                   ngridmx,llm,dtphys,
855     I                   t,paprs,pplay,
856     I                   zmfu, zmfd, zen_u, zde_u, zen_d, zde_d,
857     I                   coefkz,yu1,yv1,ftsol,pctsrf,latfi,
858     I                   frac_impa, frac_nucl,
859     I                   lonfi,presnivs,airefi,pphis,
860     O                   qfi)
861C
862             debutphy= .FALSE.
863C
864             enddo ! nsplit_phy
865C
866c   On passe le traceur sur la grille dynamique.
867             call gr_fi_dyn(llm*nq,ngridmx,iip1,jjp1,qfi,q)
868
869             endif
870
871             itauav=(itau-1)*nsplit+isplit
872c             itauav=itau*nsplit+isplit
873c            CALL writedynav(histaveid, nq,mode, itauav,vcov ,
874c    ,                   ucov,teta,pk,phi,q,masse,ps,phis)
875c            qmoy(:,:,:)=qmoy(:,:,:)+q(:,:,1,:)
876         do iq=1,nq
877            do j=1,jjp1
878               do i=1,iip1
879                  qmoy(i,j,iq)=qmoy(i,j,iq)+q(i,j,1,iq)
880               enddo
881            enddo
882         enddo
883         enddo ! isplit=1,nsplit
884
885         call histsync
886
887c-----------------------------------------------------------------------
888c    Calcul de ucov et vcov pour les sorties
889c   -------------------------------------
890
891         call massbar(masse, massebx, masseby )
892         do l=1,llm
893            do j=1,jjp1
894               do i=1,iip1
895               ucov(i,j,l)=pbaru(i,j,l)/massebx(i,j,l)*cu(i,j)*cu(i,j)
896     s         /istdyn
897               enddo
898            enddo
899
900            do j=1,jjm
901               do i=1,iip1
902               vcov(i,j,l)=pbarv(i,j,l)/masseby(i,j,l)*cv(i,j)*cv(i,j)
903     s         /istdyn
904               enddo
905            enddo
906         enddo
907         iday= dayref+itau/iday_step
908         time=
909     s   float(itau-(iday-dayref)*iday_step)/iday_step+time_0
910         IF(time.GT.1.) THEN
911            time = time-1.
912            iday = iday+1
913         ENDIF
914C
915c=======================================================================
916               do iq=iqmin,iqmax
917                  write(str2,'(i2.2)') iq
918                  nom='q'//str2
919                  call wrgrads(1,1,qmoy(:,:,iq),nom,nom)
920                  call wrgrads(2,1,q(:,:,ldec,iq),nom,nom)
921               enddo
922
923         do iq=1,nq
924            do j=1,jjp1
925               do i=1,iip1
926                  qmoy(i,j,iq)=0.
927               enddo
928            enddo
929         enddo
930
931         print*,'Record=',irec,'   stoke=',stoke,
932     s      '      injecte=',injecte
933
934      enddo ! itau
935
936      CALL dynredem1("restart.nc",0.0,
937     .          vcov,ucov,teta,q,nq,masse,ps)
938
939c=======================================================================
940C fermeture des fichiers netcdf
941          CALL histclo
942C
943 300  FORMAT('1',/,15x,'run du pas',i7,2x,'au pas',i7,2x,
944     s 'c"est a dire du jour',i7,3x,'au jour',i7,/,/)
945C----------------------------------------------------------------------
946
947      stop
948      end
Note: See TracBrowser for help on using the repository browser.