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

Last change on this file since 1130 was 1130, checked in by emillour, 11 years ago

Mars GCM:
Series of changes to enable running in parallel (using LMDZ.COMMON dynamics);
Current LMDZ.MARS can still notheless be compiled and run in serial mode
"as previously".
Summary of main changes:

  • Main programs (newstart, start2archive, xvik) that used to be in dyn3d have been moved to phymars.
  • dyn3d/control.h is now module control_mod.F90
  • rearanged input/outputs routines everywhere to handle serial/MPI cases. physdem.F => phyredem.F90 , phyetat0.F => phyetat0.F90 ; all read/write routines for startfi files are gathered in module iostart.F90
  • added parallelism related routines init_phys_lmdz.F90, comgeomphy.F90, dimphy.F90, iniphysiq.F90, mod_grid_phy_lmdz.F90, mod_phys_lmdz_mpi_data.F90, mod_phys_lmdz_mpi_transfert.F90, mod_phys_lmdz_omp_data.F90, mod_phys_lmdz_omp_transfert.F90, mod_phys_lmdz_para.F90, mod_phys_lmdz_transfert_para.F90 in phymars and mod_const_mpi.F90 in dyn3d (for compliance with parallel case)
  • created generic routines 'planetwide_maxval' and 'planetwide_minval', in module "planetwide_mod", that enable obtaining the min and max of a field over the whole planet.

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