source: trunk/LMDZ.MARS/libf/dyn3d/dynetat0.F @ 1127

Last change on this file since 1127 was 1036, checked in by emillour, 12 years ago

Mars GCM: (a first step towards using parallel dynamics)

  • IMPORTANT CHANGE: Implemented dynamic tracers. It is no longer necessary to compile the model with the '-t #' option, number of tracers is simply read from tracer.def file (as before). Adapted makegcm_* scripts (and co.) accordingly. Technical aspects of the switch to dynamic tracers are:
    • advtrac.h (in dyn3d) removed and replaced by module infotrac.F
    • tracer.h (in phymars) removed and replaced by module tracer_mod.F90 (which contains nqmx, the number of tracers, etc. and can be used anywhere in the physics).
  • Included some side cleanups: removed unused files (in dyn3d) anldoppler2.F, anl_mcdstats.F and anl_stats-diag.F, and all the unecessary dimensions.* files in grid/dimension.
  • Checked that changes are clean and that GCM yields identical results (in debug mode) to previous svn version.

EM
File size: 15.0 KB
RevLine 
[38]1      SUBROUTINE dynetat0(fichnom,nq,vcov,ucov,
[999]2     .                    teta,q,masse,ps,phis,time0)
[38]3     
4      use netcdf
[1036]5      use infotrac, only: tnom
[38]6     
7      IMPLICIT NONE
8
9c=======================================================================
10c
11c   Auteur:  P. Le Van / L.Fairhead
12c   -------
13c
14c   objet:
15c   ------
16c
17c   Lecture de l'etat initial
18c
19c   Modifs: Oct.2008 read in tracers by name. Ehouarn Millour
20c           Aug.2010 use NetCDF90 to load variables (enables using
21c                    r4 or r8 restarts independently of having compiled
22c                    the GCM in r4 or r8)
[999]23c           June 2013 TN : Possibility to read files with a time axis.
24c                          NOW defrun_new MUST BE CALLED BEFORE dynetat0.
[38]25c=======================================================================
26c-----------------------------------------------------------------------
27c   Declarations:
28c   -------------
29
30#include "dimensions.h"
31#include "paramet.h"
32#include "temps.h"
33#include "comconst.h"
34#include "comvert.h"
35#include "comgeom2.h"
36#include "ener.h"
37#include "description.h"
38#include "serre.h"
39#include "logic.h"
[1036]40!#include "advtrac.h"
[999]41#include "control.h"
[38]42
43c   Arguments:
44c   ----------
45
46      CHARACTER*(*) fichnom
47      INTEGER nq
48      REAL vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm),teta(iip1,jjp1,llm)
49      REAL q(iip1,jjp1,llm,nq)
50      REAL masse(iip1,jjp1,llm)
51      REAL ps(iip1,jjp1),phis(iip1,jjp1)
52
[999]53      REAL time0 ! fraction of day the fields correspond to
[38]54
55c   Variables
56c
57      REAL traceur(iip1,jjp1,llm) ! to temporarily store a tracer
58      INTEGER length,iq,i,j,l
59      PARAMETER (length = 100)
60      REAL tab_cntrl(length) ! tableau des parametres du run
61      INTEGER ierr, nid, nvarid, nqold
62      CHARACTER  str3*3,yes*1
[999]63     
64      REAL,ALLOCATABLE :: time(:) ! times stored in start
65      INTEGER timelen ! number of times stored in the file
66      INTEGER indextime ! index of selected time
67      !REAL  hour_ini ! fraction of day of stored date. Equivalent of day_ini, but 0=<hour_ini<1
[38]68
[999]69      INTEGER edges(4),corner(4)
70
[38]71c-----------------------------------------------------------------------
72
73c  Ouverture NetCDF du fichier etat initial
74
75      ierr=nf90_open(fichnom,NF90_NOWRITE,nid)
76      IF (ierr.NE.nf90_noerr) THEN
77        write(6,*)' Pb d''ouverture du fichier ',trim(fichnom)
78        write(*,*)trim(nf90_strerror(ierr))
79        CALL ABORT
80      ENDIF
81
82c
83      ierr=nf90_inq_varid(nid,"controle",nvarid)
84      IF (ierr .NE. nf90_noerr) THEN
85         PRINT*, "dynetat0: Le champ <controle> est absent"
86         write(*,*)trim(nf90_strerror(ierr))
87         CALL abort
88      ENDIF
89      ierr=nf90_get_var(nid,nvarid,tab_cntrl)
90      IF (ierr .NE. nf90_noerr) THEN
91         PRINT*, "dynetat0: Lecture echoue pour <controle>"
92         write(*,*)trim(nf90_strerror(ierr))
93         CALL abort
94      ENDIF
95
96      im         = tab_cntrl(1)
97      jm         = tab_cntrl(2)
98      lllm       = tab_cntrl(3)
99      day_ini    = tab_cntrl(4)
100      rad        = tab_cntrl(5)
101      omeg       = tab_cntrl(6)
102      g          = tab_cntrl(7)
103      cpp        = tab_cntrl(8)
104      kappa      = tab_cntrl(9)
105      daysec     = tab_cntrl(10)
106      dtvr       = tab_cntrl(11)
107      etot0      = tab_cntrl(12)
108      ptot0      = tab_cntrl(13)
109      ztot0      = tab_cntrl(14)
110      stot0      = tab_cntrl(15)
111      ang0       = tab_cntrl(16)
112      pa         = tab_cntrl(17)
113      preff      = tab_cntrl(18)
[999]114     
115      hour_ini   = tab_cntrl(29)
116
[38]117c
118      clon       = tab_cntrl(19)
119      clat       = tab_cntrl(20)
120      grossismx  = tab_cntrl(21)
121      grossismy  = tab_cntrl(22)
122c
123      IF ( tab_cntrl(23).EQ.1. )  THEN
124        fxyhypb  = . TRUE .
125        dzoomx   = tab_cntrl(24)
126        dzoomy   = tab_cntrl(25)
127        taux     = tab_cntrl(27)
128        tauy     = tab_cntrl(28)
129      ELSE
130        fxyhypb = . FALSE .
131        ysinus  = . FALSE .
132        IF( tab_cntrl(26).EQ.1. ) ysinus = . TRUE.
133      ENDIF
[999]134     
[38]135c   .................................................................
136c
137c
138      PRINT*,'dynetat0: rad,omeg,g,cpp,kappa',rad,omeg,g,cpp,kappa
139 
140      IF(   im.ne.iim           )  THEN
141          PRINT 1,im,iim
142          STOP
143      ELSE  IF( jm.ne.jjm       )  THEN
144          PRINT 2,jm,jjm
145          STOP
146      ELSE  IF( lllm.ne.llm     )  THEN
147          PRINT 3,lllm,llm
148          STOP
149      ENDIF
150
151      ierr=nf90_inq_varid(nid,"rlonu",nvarid)
152      IF (ierr .NE. nf90_noerr) THEN
153         PRINT*, "dynetat0: Le champ <rlonu> est absent"
154         write(*,*)trim(nf90_strerror(ierr))
155         CALL abort
156      ENDIF
157      ierr=nf90_get_var(nid,nvarid,rlonu)
158      IF (ierr .NE. nf90_noerr) THEN
159         PRINT*, "dynetat0: Lecture echouee pour <rlonu>"
160         write(*,*)trim(nf90_strerror(ierr))
161         CALL abort
162      ENDIF
163
164      ierr=nf90_inq_varid(nid,"rlatu",nvarid)
165      IF (ierr .NE. nf90_noerr) THEN
166         PRINT*, "dynetat0: Le champ <rlatu> est absent"
167         write(*,*)trim(nf90_strerror(ierr))
168         CALL abort
169      ENDIF
170      ierr=nf90_get_var(nid,nvarid,rlatu)
171      IF (ierr .NE. nf90_noerr) THEN
172         PRINT*, "dynetat0: Lecture echouee pour <rlatu>"
173         write(*,*)trim(nf90_strerror(ierr))
174         CALL abort
175      ENDIF
176
177      ierr=nf90_inq_varid(nid,"rlonv",nvarid)
178      IF (ierr .NE. nf90_noerr) THEN
179         PRINT*, "dynetat0: Le champ <rlonv> est absent"
180         write(*,*)trim(nf90_strerror(ierr))
181         CALL abort
182      ENDIF
183      ierr=nf90_get_var(nid,nvarid,rlonv)
184      IF (ierr .NE. nf90_noerr) THEN
185         PRINT*, "dynetat0: Lecture echouee pour <rlonv>"
186         write(*,*)trim(nf90_strerror(ierr))
187         CALL abort
188      ENDIF
189
190      ierr=nf90_inq_varid(nid,"rlatv",nvarid)
191      IF (ierr .NE. nf90_noerr) THEN
192         PRINT*, "dynetat0: Le champ <rlatv> est absent"
193         write(*,*)trim(nf90_strerror(ierr))
194         CALL abort
195      ENDIF
196      ierr=nf90_get_var(nid,nvarid,rlatv)
197      IF (ierr .NE. nf90_noerr) THEN
198         PRINT*, "dynetat0: Lecture echouee pour rlatv"
199         write(*,*)trim(nf90_strerror(ierr))
200         CALL abort
201      ENDIF
202
203      ierr=nf90_inq_varid(nid,"cu",nvarid)
204      IF (ierr .NE. nf90_noerr) THEN
205         PRINT*, "dynetat0: Le champ <cu> est absent"
206         write(*,*)trim(nf90_strerror(ierr))
207         CALL abort
208      ENDIF
209      ierr=nf90_get_var(nid,nvarid,cu)
210      IF (ierr .NE. nf90_noerr) THEN
211         PRINT*, "dynetat0: Lecture echouee pour <cu>"
212         write(*,*)trim(nf90_strerror(ierr))
213         CALL abort
214      ENDIF
215
216      ierr=nf90_inq_varid(nid,"cv",nvarid)
217      IF (ierr .NE. nf90_noerr) THEN
218         PRINT*, "dynetat0: Le champ <cv> est absent"
219         write(*,*)trim(nf90_strerror(ierr))
220         CALL abort
221      ENDIF
222      ierr=nf90_get_var(nid,nvarid,cv)
223      IF (ierr .NE. nf90_noerr) THEN
224         PRINT*, "dynetat0: Lecture echouee pour <cv>"
225         write(*,*)trim(nf90_strerror(ierr))
226         CALL abort
227      ENDIF
228
229      ierr=nf90_inq_varid(nid,"aire",nvarid)
230      IF (ierr .NE. nf90_noerr) THEN
231         PRINT*, "dynetat0: Le champ <aire> est absent"
232         write(*,*)trim(nf90_strerror(ierr))
233         CALL abort
234      ENDIF
235      ierr=nf90_get_var(nid,nvarid,aire)
236      IF (ierr .NE. nf90_noerr) THEN
237         PRINT*, "dynetat0: Lecture echouee pour <aire>"
238         write(*,*)trim(nf90_strerror(ierr))
239         CALL abort
240      ENDIF
241
242      ierr=nf90_inq_varid(nid,"phisinit",nvarid)
243      IF (ierr .NE. nf90_noerr) THEN
244         PRINT*, "dynetat0: Le champ <phisinit> est absent"
245         write(*,*)trim(nf90_strerror(ierr))
246         CALL abort
247      ENDIF
248      ierr=nf90_get_var(nid,nvarid,phis)
249      IF (ierr .NE. nf90_noerr) THEN
250         PRINT*, "dynetat0: Lecture echouee pour <phisinit>"
251         write(*,*)trim(nf90_strerror(ierr))
252         CALL abort
253      ENDIF
254
[999]255
256! read time axis
257      ierr = nf90_inq_dimid(nid,"Time",nvarid)
258      ierr = nf90_inquire_dimension(nid,nvarid,len=timelen)
259
260      ALLOCATE(time(timelen))
261
[38]262      ierr=nf90_inq_varid(nid,"Time",nvarid)
263      IF (ierr .NE. nf90_noerr) THEN
264             ierr=nf90_inq_varid(nid,"temps", nvarid)
265                 IF (ierr .NE. nf90_noerr) THEN
266           PRINT*, "dynetat0: <Time> or <temps> absent"
267           write(*,*)trim(nf90_strerror(ierr))
268           CALL abort
269         ENDIF
270      ENDIF
271      ierr=nf90_get_var(nid,nvarid,time)
272      IF (ierr .NE. nf90_noerr) THEN
273         PRINT*, "dynetat0: Lecture echouee <Time>/<temps>"
274         write(*,*)trim(nf90_strerror(ierr))
275         CALL abort
276      ENDIF
[999]277     
278! select desired time
279      IF (timestart .lt. 0) THEN  ! default: we use the last time value
280        indextime = timelen
281      ELSE  ! else we look for the desired value in the time axis
282       indextime = 0
283        DO i=1,timelen
284          IF (abs(time(i) - timestart) .lt. 0.01) THEN
285             indextime = i
286             EXIT
287          ENDIF
288        ENDDO
289        IF (indextime .eq. 0) THEN
290          PRINT*, "Time", timestart," is not in "//fichnom//"!!"
291          PRINT*, "Stored times are:"
292          DO i=1,timelen
293             PRINT*, time(i)
294          ENDDO
295          CALL abort
296        ENDIF
297      ENDIF
298     
299      ! In start the absolute date is day_ini + hour_ini + time
300      ! For now on, in the GCMdynamics, it is day_ini + time0
301      time0 = time(indextime) + hour_ini
302      day_ini = day_ini + INT(time0)
303      time0 = time0 - INT(time0) ! time0 devient le nouveau hour_ini
304      hour_ini = time0
305     
306      PRINT*, "dynetat0: Selected time ",time(indextime),
307     .        " at index ",indextime
308     
309      DEALLOCATE(time)
[38]310
[999]311
312! read vcov
313      corner(1)=1
314      corner(2)=1
315      corner(3)=1
316      corner(4)=indextime
317      edges(1)=iip1
318      edges(2)=jjm
319      edges(3)=llm
320      edges(4)=1
321      ierr=nf90_inq_varid(nid,"vcov",nvarid)
[38]322      IF (ierr .NE. nf90_noerr) THEN
[999]323         PRINT*, "dynetat0: Le champ <vcov> est absent"
[38]324         write(*,*)trim(nf90_strerror(ierr))
325         CALL abort
326      ENDIF
[999]327      !ierr=nf90_get_var(nid,nvarid,vcov)
328      ierr=nf90_get_var(nid,nvarid,vcov,corner,edges)
[38]329      IF (ierr .NE. nf90_noerr) THEN
[999]330         PRINT*, "dynetat0: Lecture echouee pour <vcov>"
[38]331         write(*,*)trim(nf90_strerror(ierr))
332         CALL abort
333      ENDIF
[999]334     
335! read ucov
336      corner(1)=1
337      corner(2)=1
338      corner(3)=1
339      corner(4)=indextime
340      edges(1)=iip1
341      edges(2)=jjp1
342      edges(3)=llm
343      edges(4)=1
344      ierr=nf90_inq_varid(nid,"ucov",nvarid)
[38]345      IF (ierr .NE. nf90_noerr) THEN
[999]346         PRINT*, "dynetat0: Le champ <ucov> est absent"
[38]347         write(*,*)trim(nf90_strerror(ierr))
348         CALL abort
349      ENDIF
[999]350      ierr=nf90_get_var(nid,nvarid,ucov,corner,edges)
[38]351      IF (ierr .NE. nf90_noerr) THEN
[999]352         PRINT*, "dynetat0: Lecture echouee pour <ucov>"
[38]353         write(*,*)trim(nf90_strerror(ierr))
354         CALL abort
355      ENDIF
356
[999]357! read teta
[38]358      ierr=nf90_inq_varid(nid,"teta",nvarid)
359      IF (ierr .NE. nf90_noerr) THEN
360         PRINT*, "dynetat0: Le champ <teta> est absent"
361         write(*,*)trim(nf90_strerror(ierr))
362         CALL abort
363      ENDIF
[999]364      ierr=nf90_get_var(nid,nvarid,teta,corner,edges)
[38]365      IF (ierr .NE. nf90_noerr) THEN
366         PRINT*, "dynetat0: Lecture echouee pour <teta>"
367         write(*,*)trim(nf90_strerror(ierr))
368         CALL abort
369      ENDIF
370
[999]371! read tracers
[38]372      IF(nq.GE.1) THEN
373        write(*,*) 'dynetat0: loading tracers'
374         IF(nq.GT.99) THEN
375            PRINT*, "Trop de traceurs"
376            CALL abort
377         ENDIF
378         nqold=nq
379         DO iq=1,nq
380!           str3(1:1)='q'
381!           WRITE(str3(2:3),'(i2.2)') iq
382!           ierr =  NF_INQ_VARID (nid, str3, nvarid)
[1036]383! NB: tracers are now read in using their name ('tnom' from infotrac)
[38]384!           write(*,*) "  loading tracer:",trim(tnom(iq))
385           ierr=nf90_inq_varid(nid,tnom(iq),nvarid)
386           IF (ierr .NE. nf90_noerr) THEN
387!              PRINT*, "dynetat0: Le champ <"//str3//"> est absent"
388              PRINT*, "dynetat0: Le champ <"//trim(tnom(iq))//
389     &                "> est absent"
390              PRINT*, "          Il est donc initialise a zero"
391              CALL initial0(ijp1llm,q(1,1,1,iq))
392              nqold=min(iq-1,nqold)
393           ELSE
394           !ierr=nf90_get_var(nid,nvarid,q(1,1,1,iq))
[999]395           ierr=nf90_get_var(nid,nvarid,traceur,corner,edges)
[38]396             IF (ierr .NE. nf90_noerr) THEN
397!                 PRINT*, "dynetat0: Lecture echouee pour "//str3
398               PRINT*, "dynetat0: Lecture echouee pour "//trim(tnom(iq))
399               CALL abort
400             ENDIF
401             q(:,:,:,iq)=traceur(:,:,:)
402           ENDIF
403         ENDDO
404         if ((nqold.lt.nq).and.(nqold.ge.1)) then   
405c        case when new tracer are added in addition to old ones
406             write(*,*)'tracers 1 to ', nqold,'were already present'
[1036]407             write(*,*)'tracers ', nqold+1,' to ', nq,'are new'
[38]408!             yes=' '
409!            do while ((yes.ne.'y').and.(yes.ne.'n'))
410!             write(*,*) 'Would you like to reindex tracer # 1 ->',nqold
[1036]411!             write(*,*) 'to #',nq-nqold+1,'->', nq,'   (y or n) ?'
[38]412!             read(*,fmt='(a)') yes
413!            end do
414!            if (yes.eq.'y') then
415!              write(*,*) 'OK, let s reindex the tracers'
416!              do l=1,llm
417!                do j=1,jjp1
418!                  do i=1,iip1
[1036]419!                    do iq=nq,nq-nqold+1,-1
420!                       q(i,j,l,iq)=q(i,j,l,iq-nq+nqold)   
[38]421!                    end do
[1036]422!                    do iq=nq-nqold,1,-1
[38]423!                       q(i,j,l,iq)= 0.
424!                    end do
425!                  end do
426!                end do
427!              end do
428!            end if
429         end if
430      ENDIF
431
[999]432! read masse
[38]433      ierr=nf90_inq_varid(nid,"masse",nvarid)
434      IF (ierr .NE. nf90_noerr) THEN
435         PRINT*, "dynetat0: Le champ <masse> est absent"
436         write(*,*)trim(nf90_strerror(ierr))
437         CALL abort
438      ENDIF
[999]439      ierr=nf90_get_var(nid,nvarid,masse,corner,edges)
[38]440      IF (ierr .NE. nf90_noerr) THEN
441         PRINT*, "dynetat0: Lecture echouee pour <masse>"
442         write(*,*)trim(nf90_strerror(ierr))
443         CALL abort
444      ENDIF
445
[999]446! read ps
447      corner(1)=1
448      corner(2)=1
449      corner(3)=indextime
450      edges(1)=iip1
451      edges(2)=jjp1
452      edges(3)=1
[38]453      ierr=nf90_inq_varid(nid,"ps",nvarid)
454      IF (ierr .NE. nf90_noerr) THEN
455         PRINT*, "dynetat0: Le champ <ps> est absent"
456         write(*,*)trim(nf90_strerror(ierr))
457         CALL abort
458      ENDIF
[999]459      ierr=nf90_get_var(nid,nvarid,ps,corner,edges)
[38]460      IF (ierr .NE. nf90_noerr) THEN
461         PRINT*, "dynetat0: Lecture echouee pour <ps>"
462         write(*,*)trim(nf90_strerror(ierr))
463         CALL abort
464      ENDIF
465
466      ierr=nf90_close(nid)
467
[791]468      ! day_ini=day_ini+INT(time) ! obsolete stuff ; 0<time<1 anyways
469      ! time=time-INT(time)
[38]470
471  1   FORMAT(//10x,'la valeur de im =',i4,2x,'lue sur le fichier de dem
472     *arrage est differente de la valeur parametree iim =',i4//)
473   2  FORMAT(//10x,'la valeur de jm =',i4,2x,'lue sur le fichier de dem
474     *arrage est differente de la valeur parametree jjm =',i4//)
475   3  FORMAT(//10x,'la valeur de lmax =',i4,2x,'lue sur le fichier dema
476     *rrage est differente de la valeur parametree llm =',i4//)
477   4  FORMAT(//10x,'la valeur de dtrv =',i4,2x,'lue sur le fichier dema
478     *rrage est differente de la valeur  dtinteg =',i4//)
479
480      RETURN
481      END
Note: See TracBrowser for help on using the repository browser.