source: trunk/LMDZ.COMMON/libf/dyn3d_common/dynetat0.F @ 1436

Last change on this file since 1436 was 1422, checked in by milmd, 11 years ago

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