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
RevLine 
[1]1!
2! $Id $
3!
4      SUBROUTINE dynetat0(fichnom,vcov,ucov,
[1107]5     .                    teta,q,masse,ps,phis,time0)
[1]6
[1107]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
[841]11
[1107]12      use control_mod, only : planet_type, timestart
[1422]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
[841]21
[1]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"
[776]41#include "comgeom2.h"
[1]42#include "netcdf.inc"
43#include "iniprint.h"
44
45c   Arguments:
46c   ----------
47
[1107]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
[1]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
[907]65      character(len=12) :: start_file_type="earth" ! default start file type
[841]66      INTEGER idecal
67
[1107]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
[1]77c-----------------------------------------------------------------------
78
79c  Ouverture NetCDF du fichier etat initial
80
[1107]81      ierr=nf90_open(fichnom,NF90_NOWRITE,nid)
82      IF (ierr.NE.nf90_noerr) THEN
[1]83        write(lunout,*)'dynetat0: Pb d''ouverture du fichier start.nc'
[1107]84        write(lunout,*)trim(nf90_strerror(ierr))
[1300]85        CALL ABORT_gcm("dynetat0", "", 1)
[1]86      ENDIF
87
88c
[1107]89      ierr = nf90_inq_varid (nid, "controle", nvarid)
90      IF (ierr .NE. nf90_noerr) THEN
[1]91         write(lunout,*)"dynetat0: Le champ <controle> est absent"
[1107]92         write(lunout,*)trim(nf90_strerror(ierr))
[1300]93         CALL ABORT_gcm("dynetat0", "", 1)
[1]94      ENDIF
[776]95      ierr = nf90_get_var(nid, nvarid, tab_cntrl)
[1107]96      IF (ierr .NE. nf90_noerr) THEN
[1]97         write(lunout,*)"dynetat0: Lecture echoue pour <controle>"
[1107]98         write(lunout,*)trim(nf90_strerror(ierr))
[1300]99         CALL ABORT_gcm("dynetat0", "", 1)
[1]100      ENDIF
101
[841]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'
[907]106          start_file_type="planeto"
[841]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
[1]116      im         = tab_cntrl(1)
117      jm         = tab_cntrl(2)
118      lllm       = tab_cntrl(3)
[907]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
[841]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)
[1]139c
[841]140      clon       = tab_cntrl(idecal+15)
141      clat       = tab_cntrl(idecal+16)
142      grossismx  = tab_cntrl(idecal+17)
143      grossismy  = tab_cntrl(idecal+18)
[1]144c
[841]145      IF ( tab_cntrl(idecal+19).EQ.1. )  THEN
[1]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 .
[841]154        IF( tab_cntrl(idecal+22).EQ.1. ) ysinus = . TRUE.
[1]155      ENDIF
156
[1107]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
[907]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
[1]172c   .................................................................
173c
174c
[1107]175      write(lunout,*)'dynetat0: rad,omeg,g,cpp,kappa ',
[1]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
[1107]189      ierr=nf90_inq_varid(nid, "rlonu", nvarid)
190      IF (ierr .NE. nf90_noerr) THEN
[1]191         write(lunout,*)"dynetat0: Le champ <rlonu> est absent"
[1107]192         write(lunout,*)trim(nf90_strerror(ierr))
[1300]193         CALL ABORT_gcm("dynetat0", "", 1)
[1]194      ENDIF
[776]195      ierr = nf90_get_var(nid, nvarid, rlonu)
[1107]196      IF (ierr .NE. nf90_noerr) THEN
[1]197         write(lunout,*)"dynetat0: Lecture echouee pour <rlonu>"
[1107]198         write(lunout,*)trim(nf90_strerror(ierr))
[1300]199         CALL ABORT_gcm("dynetat0", "", 1)
[1]200      ENDIF
201
[1107]202      ierr = nf90_inq_varid (nid, "rlatu", nvarid)
203      IF (ierr .NE. nf90_noerr) THEN
[1]204         write(lunout,*)"dynetat0: Le champ <rlatu> est absent"
[1107]205         write(lunout,*)trim(nf90_strerror(ierr))
[1300]206         CALL ABORT_gcm("dynetat0", "", 1)
[1]207      ENDIF
[776]208      ierr = nf90_get_var(nid, nvarid, rlatu)
[1107]209      IF (ierr .NE. nf90_noerr) THEN
[1]210         write(lunout,*)"dynetat0: Lecture echouee pour <rlatu>"
[1107]211         write(lunout,*)trim(nf90_strerror(ierr))
[1300]212         CALL ABORT_gcm("dynetat0", "", 1)
[1]213      ENDIF
214
[1107]215      ierr = nf90_inq_varid (nid, "rlonv", nvarid)
216      IF (ierr .NE. nf90_noerr) THEN
[1]217         write(lunout,*)"dynetat0: Le champ <rlonv> est absent"
[1107]218         write(lunout,*)trim(nf90_strerror(ierr))
[1300]219         CALL ABORT_gcm("dynetat0", "", 1)
[1]220      ENDIF
[776]221      ierr = nf90_get_var(nid, nvarid, rlonv)
[1107]222      IF (ierr .NE. nf90_noerr) THEN
[1]223         write(lunout,*)"dynetat0: Lecture echouee pour <rlonv>"
[1107]224         write(lunout,*)trim(nf90_strerror(ierr))
[1300]225         CALL ABORT_gcm("dynetat0", "", 1)
[1]226      ENDIF
227
[1107]228      ierr = nf90_inq_varid (nid, "rlatv", nvarid)
229      IF (ierr .NE. nf90_noerr) THEN
[1]230         write(lunout,*)"dynetat0: Le champ <rlatv> est absent"
[1107]231         write(lunout,*)trim(nf90_strerror(ierr))
[1300]232         CALL ABORT_gcm("dynetat0", "", 1)
[1]233      ENDIF
[776]234      ierr = nf90_get_var(nid, nvarid, rlatv)
[1107]235      IF (ierr .NE. nf90_noerr) THEN
[1]236         write(lunout,*)"dynetat0: Lecture echouee pour rlatv"
[1107]237         write(lunout,*)trim(nf90_strerror(ierr))
[1300]238         CALL ABORT_gcm("dynetat0", "", 1)
[1]239      ENDIF
240
[1107]241      ierr = nf90_inq_varid (nid, "cu", nvarid)
242      IF (ierr .NE. nf90_noerr) THEN
[1]243         write(lunout,*)"dynetat0: Le champ <cu> est absent"
[1107]244         write(lunout,*)trim(nf90_strerror(ierr))
[1300]245         CALL ABORT_gcm("dynetat0", "", 1)
[1]246      ENDIF
[776]247      ierr = nf90_get_var(nid, nvarid, cu)
[1107]248      IF (ierr .NE. nf90_noerr) THEN
[1]249         write(lunout,*)"dynetat0: Lecture echouee pour <cu>"
[1107]250         write(lunout,*)trim(nf90_strerror(ierr))
[1300]251         CALL ABORT_gcm("dynetat0", "", 1)
[1]252      ENDIF
253
[1107]254      ierr = nf90_inq_varid (nid, "cv", nvarid)
255      IF (ierr .NE. nf90_noerr) THEN
[1]256         write(lunout,*)"dynetat0: Le champ <cv> est absent"
[1107]257         write(lunout,*)trim(nf90_strerror(ierr))
[1300]258         CALL ABORT_gcm("dynetat0", "", 1)
[1]259      ENDIF
[776]260      ierr = nf90_get_var(nid, nvarid, cv)
[1107]261      IF (ierr .NE. nf90_noerr) THEN
[1]262         write(lunout,*)"dynetat0: Lecture echouee pour <cv>"
[1107]263         write(lunout,*)trim(nf90_strerror(ierr))
[1300]264         CALL ABORT_gcm("dynetat0", "", 1)
[1]265      ENDIF
266
[1107]267      ierr = nf90_inq_varid (nid, "aire", nvarid)
268      IF (ierr .NE. nf90_noerr) THEN
[1]269         write(lunout,*)"dynetat0: Le champ <aire> est absent"
[1107]270         write(lunout,*)trim(nf90_strerror(ierr))
[1300]271         CALL ABORT_gcm("dynetat0", "", 1)
[1]272      ENDIF
[776]273      ierr = nf90_get_var(nid, nvarid, aire)
[1107]274      IF (ierr .NE. nf90_noerr) THEN
[1]275         write(lunout,*)"dynetat0: Lecture echouee pour <aire>"
[1107]276         write(lunout,*)trim(nf90_strerror(ierr))
[1300]277         CALL ABORT_gcm("dynetat0", "", 1)
[1]278      ENDIF
279
[1107]280      ierr = nf90_inq_varid (nid, "phisinit", nvarid)
281      IF (ierr .NE. nf90_noerr) THEN
[1]282         write(lunout,*)"dynetat0: Le champ <phisinit> est absent"
[1107]283         write(lunout,*)trim(nf90_strerror(ierr))
[1300]284         CALL ABORT_gcm("dynetat0", "", 1)
[1]285      ENDIF
[776]286      ierr = nf90_get_var(nid, nvarid, phis)
[1107]287      IF (ierr .NE. nf90_noerr) THEN
[1]288         write(lunout,*)"dynetat0: Lecture echouee pour <phisinit>"
[1107]289         write(lunout,*)trim(nf90_strerror(ierr))
[1300]290         CALL ABORT_gcm("dynetat0", "", 1)
[1]291      ENDIF
292
[1107]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))
[1300]302           CALL ABORT_gcm("dynetat0", "", 1)
[1107]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))
[1300]314           CALL ABORT_gcm("dynetat0", "", 1)
[1107]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))
[1300]327           CALL ABORT_gcm("dynetat0", "", 1)
[1107]328        ENDIF
[1]329      ENDIF
[1107]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
[1300]349          CALL ABORT_gcm("dynetat0", "", 1)
[1107]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))
[1300]382         CALL ABORT_gcm("dynetat0", "", 1)
[1]383      ENDIF
[1107]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))
[1300]388         CALL ABORT_gcm("dynetat0", "", 1)
[1107]389      ENDIF
[1]390
[1107]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
[1]402         write(lunout,*)"dynetat0: Le champ <ucov> est absent"
[1107]403         write(lunout,*)trim(nf90_strerror(ierr))
[1300]404         CALL ABORT_gcm("dynetat0", "", 1)
[1]405      ENDIF
[1107]406      ierr=nf90_get_var(nid,nvarid,ucov,corner,edges)
407      IF (ierr .NE. nf90_noerr) THEN
[1]408         write(lunout,*)"dynetat0: Lecture echouee pour <ucov>"
[1107]409         write(lunout,*)trim(nf90_strerror(ierr))
[1300]410         CALL ABORT_gcm("dynetat0", "", 1)
[1]411      ENDIF
412 
[1107]413! read teta (same corner/edges as ucov)
414      ierr=nf90_inq_varid(nid,"teta",nvarid)
415      IF (ierr .NE. nf90_noerr) THEN
[1]416         write(lunout,*)"dynetat0: Le champ <teta> est absent"
[1107]417         write(lunout,*)trim(nf90_strerror(ierr))
[1300]418         CALL ABORT_gcm("dynetat0", "", 1)
[1]419      ENDIF
[1107]420      ierr=nf90_get_var(nid,nvarid,teta,corner,edges)
421      IF (ierr .NE. nf90_noerr) THEN
[1]422         write(lunout,*)"dynetat0: Lecture echouee pour <teta>"
[1107]423         write(lunout,*)trim(nf90_strerror(ierr))
[1300]424         CALL ABORT_gcm("dynetat0", "", 1)
[1]425      ENDIF
426
[1107]427! read tracers (same corner/edges as ucov)
[1]428      IF(nqtot.GE.1) THEN
429      DO iq=1,nqtot
[1107]430        ierr= nf90_inq_varid(nid,tname(iq),nvarid)
431        IF (ierr .NE. nf90_noerr) THEN
[1]432           write(lunout,*)"dynetat0: Le traceur <"//trim(tname(iq))//
433     &                    "> est absent"
434           write(lunout,*)"          Il est donc initialise a zero"
[776]435           q(:,:,:,iq)=0.
[1]436        ELSE
[1107]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))
[1300]442            CALL ABORT_gcm("dynetat0", "", 1)
[1]443          ENDIF
444        ENDIF
445      ENDDO
446      ENDIF
447
[1107]448!read masse (same corner/edges as ucov)
449      ierr=nf90_inq_varid(nid,"masse",nvarid)
450      IF (ierr .NE. nf90_noerr) THEN
[1]451         write(lunout,*)"dynetat0: Le champ <masse> est absent"
[1107]452         write(lunout,*)trim(nf90_strerror(ierr))
[1300]453         CALL ABORT_gcm("dynetat0", "", 1)
[1]454      ENDIF
[1107]455      ierr=nf90_get_var(nid,nvarid,masse,corner,edges)
456      IF (ierr .NE. nf90_noerr) THEN
[1]457         write(lunout,*)"dynetat0: Lecture echouee pour <masse>"
[1107]458         write(lunout,*)trim(nf90_strerror(ierr))
[1300]459         CALL ABORT_gcm("dynetat0", "", 1)
[1]460      ENDIF
461
[1107]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
[1]471         write(lunout,*)"dynetat0: Le champ <ps> est absent"
[1107]472         write(lunout,*)trim(nf90_strerror(ierr))
[1300]473         CALL ABORT_gcm("dynetat0", "", 1)
[1]474      ENDIF
[1107]475      ierr=nf90_get_var(nid,nvarid,ps,corner,edges)
476      IF (ierr .NE. nf90_noerr) THEN
[1]477         write(lunout,*)"dynetat0: Lecture echouee pour <ps>"
[1107]478         write(lunout,*)trim(nf90_strerror(ierr))
[1300]479         CALL ABORT_gcm("dynetat0", "", 1)
[1]480      ENDIF
481
[1107]482      ierr=nf90_close(nid)
[1]483
[1107]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
[1]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.