source: trunk/LMDZ.MARS/libf/phymars/phyetat0.F @ 1009

Last change on this file since 1009 was 999, checked in by tnavarro, 11 years ago

Possibility to store multiple initial states in one start/startfi. This is RETROCOMPATIBLE. New option ecrithist in run.def to write data in start/startfi every ecrithist dynamical timestep. New option timestart in run.def to initialize the GCM with the time timestart stored in start
File size: 21.2 KB
RevLine 
[38]1      SUBROUTINE phyetat0 (fichnom,tab0,Lmodif,nsoil,nq,
[999]2     .           day_ini,time0,
[38]3     .           tsurf,tsoil,emis,q2,qsurf,co2ice)
4
5      use netcdf
6
7      implicit none
8c======================================================================
9c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
10c  Adaptation à Mars : Yann Wanherdrick
11c Objet: Lecture de l etat initial pour la physique
12! Modifs: Aug.2010 EM : use NetCDF90 to load variables (enables using
13!                      r4 or r8 restarts independently of having compiled
14!                      the GCM in r4 or r8)
[999]15!         June 2013 TN : Possibility to read files with a time axis
[38]16c======================================================================
17!#include "netcdf.inc"
18#include "dimensions.h"
19#include "dimphys.h"
[697]20!#include "comgeomfi.h"
[38]21#include "surfdat.h"
22#include "planete.h"
23#include "dimradmars.h"
24#include "yomaer.h"
25#include "comcstfi.h"
26!#include "tracer.h"
[999]27#include "advtrac.h"
28#include "control.h"
[38]29c======================================================================
30      INTEGER nbsrf !Mars nbsrf a 1 au lieu de 4
31      PARAMETER (nbsrf=1) ! nombre de sous-fractions pour une maille
32!======================================================================
33!  Arguments:
34!  ---------
35!  inputs:
36      character*(*) fichnom ! "startfi.nc" file
37      integer tab0
38      integer Lmodif
39      integer nsoil ! # of soil layers
40      integer nq
41      integer day_ini
[999]42      real time0
[38]43
44!  outputs:
45      real tsurf(ngridmx) ! surface temperature
46      real tsoil(ngridmx,nsoil) ! soil temperature
47      real emis(ngridmx) ! surface emissivity
48      real q2(ngridmx, llm+1) !
49      real qsurf(ngridmx,nq) ! tracers on surface
50      real co2ice(ngridmx) ! co2 ice cover
51
52!======================================================================
53!  Local variables:
54
55      real surffield(ngridmx) ! to temporarily store a surface field
56      real xmin,xmax ! to display min and max of a field
57c
58      INTEGER ig,iq,lmax
59      INTEGER nid, nvarid
60      INTEGER ierr, i, nsrf
61!      integer isoil
62!      INTEGER length
63!      PARAMETER (length=100)
64      CHARACTER*7 str7
65      CHARACTER*2 str2
66      CHARACTER*1 yes
67c
68      REAL p_rad,p_omeg,p_g,p_mugaz,p_daysec
69      INTEGER nqold
70
71! flag which identifies if 'startfi.nc' file is using old names (qsurf01,...)
72      logical :: oldtracernames=.false.
73      integer :: count
74      character(len=30) :: txt ! to store some text
75
[999]76! specific for time
77      REAL,ALLOCATABLE :: time(:) ! times stored in start
78      INTEGER timelen ! number of times stored in the file
79      INTEGER indextime ! index of selected time
80
81      INTEGER edges(3),corner(3)
82
83
[38]84c
85c Ouvrir le fichier contenant l etat initial:
86c
87
88      ierr=nf90_open(fichnom,NF90_NOWRITE,nid)
89      IF (ierr.NE.nf90_noerr) THEN
90        write(6,*)' Pb d''ouverture du fichier ',trim(fichnom)
91        write(*,*)trim(nf90_strerror(ierr))
92        CALL ABORT
93      ENDIF
94
95! Preliminary stuff: check if tracers follow old naming convention (qsurf01,
96!                    qsurf02, ...)
97      count=0
98      do iq=1,nqmx
99        txt= " "
100        write(txt,'(a5,i2.2)')'qsurf',iq
101        ierr=nf90_inq_varid(nid,txt,nvarid)
102        if (ierr.ne.nf90_noerr) then
103          ! did not find old tracer name
104          exit ! might as well stop here
105        else
106          ! found old tracer name
107          count=count+1
108        endif
109      enddo
110      if (count.eq.nqmx) then
111        write(*,*) "phyetat0:tracers seem to follow old naming ",
112     &             "convention (qsurf01,qsurf02,...)"
113        write(*,*) "   => will work for now ... "
114        write(*,*) "      but you should run newstart to rename them"
115        oldtracernames=.true.
116      endif
117
118c modifications possibles des variables de tab_cntrl
119      PRINT*
120      write(*,*) 'TABFI de phyeta0',Lmodif,tab0
121      call tabfi (nid,Lmodif,tab0,day_ini,lmax,p_rad,
[999]122     .              p_omeg,p_g,p_mugaz,p_daysec,time0)
[38]123c
[697]124c Read latitudes (coordinates): No need, these are provided by the dynamics
[38]125c
[697]126!      ierr=nf90_inq_varid(nid,"latitude",nvarid)
127!      IF (ierr.NE.nf90_noerr) THEN
128!         PRINT*, 'phyetat0: Le champ <latitude> est absent'
129!         write(*,*)trim(nf90_strerror(ierr))
130!         CALL abort
131!      ENDIF
132!      ierr=nf90_get_var(nid,nvarid,lati)
133!      IF (ierr.NE.nf90_noerr) THEN
134!         PRINT*, 'phyetat0: Lecture echouee pour <latitude>'
135!         write(*,*)trim(nf90_strerror(ierr))
136!         CALL abort
137!      ENDIF
[38]138c
[697]139c read longitudes (coordinates): No need, these are provided by the dynamics
[38]140c
[697]141!      ierr=nf90_inq_varid(nid,"longitude",nvarid)
142!      IF (ierr.NE.nf90_noerr) THEN
143!         PRINT*, 'phyetat0: Le champ <longitude> est absent'
144!         write(*,*)trim(nf90_strerror(ierr))
145!         CALL abort
146!      ENDIF
147!      ierr=nf90_get_var(nid,nvarid,long)
148!      IF (ierr.NE.nf90_noerr) THEN
149!         PRINT*, 'phyetat0: Lecture echouee pour <longitude>'
150!         write(*,*)trim(nf90_strerror(ierr))
151!         CALL abort
152!      ENDIF
[38]153c
[697]154c Read areas of meshes: No need, these are provided by the dynamics
[38]155c
[697]156!      ierr=nf90_inq_varid(nid,"area",nvarid)
157!      IF (ierr.NE.nf90_noerr) THEN
158!         PRINT*, 'phyetat0: Le champ <area> est absent'
159!         write(*,*)trim(nf90_strerror(ierr))
160!         CALL abort
161!      ENDIF
162!      ierr=nf90_get_var(nid,nvarid,area)
163!      IF (ierr.NE.nf90_noerr) THEN
164!         PRINT*, 'phyetat0: Lecture echouee pour <area>'
165!         write(*,*)trim(nf90_strerror(ierr))
166!         CALL abort
167!      ENDIF
168!      xmin = 1.0E+20
169!      xmax = -1.0E+20
170!      xmin = MINVAL(area)
171!      xmax = MAXVAL(area)
172!      PRINT*,'Aires des mailles <area>:', xmin, xmax
[38]173c
174c Lecture du geopotentiel au sol:
175c
176      ierr=nf90_inq_varid(nid,"phisfi",nvarid)
177      IF (ierr.NE.nf90_noerr) THEN
178         PRINT*, 'phyetat0: Le champ <phisfi> est absent'
179         write(*,*)trim(nf90_strerror(ierr))
180         CALL abort
181      ENDIF
182      ierr=nf90_get_var(nid,nvarid,phisfi)
183      IF (ierr.NE.nf90_noerr) THEN
184         PRINT*, 'phyetat0: Lecture echouee pour <phisfi>'
185         write(*,*)trim(nf90_strerror(ierr))
186         CALL abort
187      ENDIF
188      xmin = 1.0E+20
189      xmax = -1.0E+20
190      xmin = MINVAL(phisfi)
191      xmax = MAXVAL(phisfi)
192      PRINT*,'Geopotentiel au sol <phisfi>:', xmin, xmax
193c
194c Lecture de l''albedo du sol nu:
195c
196      ierr=nf90_inq_varid(nid,"albedodat",nvarid)
197      IF (ierr.NE.nf90_noerr) THEN
198         PRINT*, 'phyetat0: Le champ <albedodat> est absent'
199         write(*,*)trim(nf90_strerror(ierr))
200         CALL abort
201      ENDIF
202      ierr=nf90_get_var(nid,nvarid,albedodat)
203      IF (ierr.NE.nf90_noerr) THEN
204         PRINT*, 'phyetat0: Lecture echouee pour <albedodat>'
205         CALL abort
206      ENDIF
207      xmin = 1.0E+20
208      xmax = -1.0E+20
209      xmin = MINVAL(albedodat)
210      xmax = MAXVAL(albedodat)
211      PRINT*,'Albedo du sol nu <albedodat>:', xmin, xmax
212c
213c Lecture de l''inertie thermique du sol:
214c
215!      ierr = NF_INQ_VARID (nid, "inertiedat", nvarid)
216!      IF (ierr.NE.nf90_noerr) THEN
217!         PRINT*, 'phyetat0: Le champ <inertiedat> est absent'
218!         CALL abort
219!      ENDIF
220!#ifdef NC_DOUBLE
221!      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, inertiedat)
222!#else
223!      ierr = NF_GET_VAR_REAL(nid, nvarid, inertiedat)
224!#endif
225!      IF (ierr.NE.nf90_noerr) THEN
226!         PRINT*, 'phyetat0: Lecture echouee pour <inertiedat>'
227!         CALL abort
228!      ENDIF
229!      xmin = 1.0E+20
230!      xmax = -1.0E+20
231!      xmin = MINVAL(inertiedat)
232!      xmax = MAXVAL(inertiedat)
233!      PRINT*,'Inertie thermique du sol <inertiedat>:', xmin, xmax
234c
235c ZMEA
236c
237      ierr=nf90_inq_varid(nid,"ZMEA",nvarid)
238      IF (ierr.NE.nf90_noerr) THEN
239         PRINT*, 'phyetat0: Le champ <ZMEA> est absent'
240         write(*,*)trim(nf90_strerror(ierr))
241         CALL abort
242      ENDIF
243      ierr=nf90_get_var(nid,nvarid,zmea)
244      IF (ierr.NE.nf90_noerr) THEN
245         PRINT*, 'phyetat0: Lecture echouee pour <ZMEA>'
246         write(*,*)trim(nf90_strerror(ierr))
247         CALL abort
248      ENDIF
249      xmin = 1.0E+20
250      xmax = -1.0E+20
251      DO i = 1, ngridmx
252         xmin = MIN(zmea(i),xmin)
253         xmax = MAX(zmea(i),xmax)
254      ENDDO
255      PRINT*,'<zmea>:', xmin, xmax
256c
257c ZSTD
258c
259      ierr=nf90_inq_varid(nid,"ZSTD",nvarid)
260      IF (ierr.NE.nf90_noerr) THEN
261         PRINT*, 'phyetat0: Le champ <ZSTD> est absent'
262         write(*,*)trim(nf90_strerror(ierr))
263         CALL abort
264      ENDIF
265      ierr=nf90_get_var(nid,nvarid,zstd)
266      IF (ierr.NE.nf90_noerr) THEN
267         PRINT*, 'phyetat0: Lecture echouee pour <ZSTD>'
268         write(*,*)trim(nf90_strerror(ierr))
269         CALL abort
270      ENDIF
271      xmin = 1.0E+20
272      xmax = -1.0E+20
273      DO i = 1, ngridmx
274         xmin = MIN(zstd(i),xmin)
275         xmax = MAX(zstd(i),xmax)
276      ENDDO
277      PRINT*,'<zstd>:', xmin, xmax
278c
279c ZSIG
280c
281      ierr=nf90_inq_varid(nid,"ZSIG",nvarid)
282      IF (ierr.NE.nf90_noerr) THEN
283         PRINT*, 'phyetat0: Le champ <ZSIG> est absent'
284         write(*,*)trim(nf90_strerror(ierr))
285         CALL abort
286      ENDIF
287      ierr=nf90_get_var(nid,nvarid,zsig)
288      IF (ierr.NE.nf90_noerr) THEN
289         PRINT*, 'phyetat0: Lecture echouee pour <ZSIG>'
290         write(*,*)trim(nf90_strerror(ierr))
291         CALL abort
292      ENDIF
293      xmin = 1.0E+20
294      xmax = -1.0E+20
295      DO i = 1, ngridmx
296         xmin = MIN(zsig(i),xmin)
297         xmax = MAX(zsig(i),xmax)
298      ENDDO
299      PRINT*,'<zsig>:', xmin, xmax
300c
301c ZGAM
302c
303      ierr=nf90_inq_varid(nid,"ZGAM",nvarid)
304      IF (ierr.NE.nf90_noerr) THEN
305         PRINT*, 'phyetat0: Le champ <ZGAM> est absent'
306         write(*,*)trim(nf90_strerror(ierr))
307         CALL abort
308      ENDIF
309      ierr=nf90_get_var(nid,nvarid,zgam)
310      IF (ierr.NE.nf90_noerr) THEN
311         PRINT*, 'phyetat0: Lecture echouee pour <ZGAM>'
312         write(*,*)trim(nf90_strerror(ierr))
313         CALL abort
314      ENDIF
315      xmin = 1.0E+20
316      xmax = -1.0E+20
317      DO i = 1, ngridmx
318         xmin = MIN(zgam(i),xmin)
319         xmax = MAX(zgam(i),xmax)
320      ENDDO
321      PRINT*,'<zgam>:', xmin, xmax
322c
323c ZTHE
324c
325      ierr=nf90_inq_varid(nid,"ZTHE",nvarid)
326      IF (ierr.NE.nf90_noerr) THEN
327         PRINT*, 'phyetat0: Le champ <ZTHE> est absent'
328         write(*,*)trim(nf90_strerror(ierr))
329         CALL abort
330      ENDIF
331      ierr=nf90_get_var(nid,nvarid,zthe)
332      IF (ierr.NE.nf90_noerr) THEN
333         PRINT*, 'phyetat0: Lecture echouee pour <ZTHE>'
334         write(*,*)trim(nf90_strerror(ierr))
335         CALL abort
336      ENDIF
337      xmin = 1.0E+20
338      xmax = -1.0E+20
339      DO i = 1, ngridmx
340         xmin = MIN(zthe(i),xmin)
341         xmax = MAX(zthe(i),xmax)
342      ENDDO
343      PRINT*,'<zthe>:', xmin, xmax
[999]344     
[38]345c
[999]346c Time axis
347c
348      ierr = nf90_inq_dimid(nid,"Time",nvarid)
349
350      IF (ierr .NE. nf90_noerr) THEN
351        indextime = 1
352        print*, "phyetat0: No time axis found in "//fichnom
353     
354      ELSE
355        print*, "phyetat0: Time axis found in "//fichnom
356       
357        ierr = nf90_inquire_dimension(nid,nvarid,len=timelen)
358       
359        ALLOCATE(time(timelen))
360 
361        ierr=nf90_inq_varid(nid,"Time",nvarid)
362        IF (ierr .NE. nf90_noerr) THEN
363             ierr=nf90_inq_varid(nid,"temps", nvarid)
364                 IF (ierr .NE. nf90_noerr) THEN
365             PRINT*, "dynetat0: <Time> or <temps> absent"
366             write(*,*)trim(nf90_strerror(ierr))
367             CALL abort
368           ENDIF
369        ENDIF
370        ierr=nf90_get_var(nid,nvarid,time)
371        IF (ierr .NE. nf90_noerr) THEN
372           PRINT*, "dynetat0: Lecture echouee <Time>/<temps>"
373           write(*,*)trim(nf90_strerror(ierr))
374           CALL abort
375        ENDIF
376c     
377c Select desired time
378c
379        IF (timestart .lt. 0) THEN  ! default: we use the last time value
380          indextime = timelen
381        ELSE  ! else we look for the desired value in the time axis
382         indextime = 0
383          DO i=1,timelen
384            IF (abs(time(i) - timestart) .lt. 0.01) THEN
385               indextime = i
386               EXIT
387            ENDIF
388          ENDDO
389          IF (indextime .eq. 0) THEN
390            PRINT*, "Time", timestart," is not in "//fichnom//"!!"
391            PRINT*, "Stored times are:"
392            DO i=1,timelen
393               PRINT*, time(i)
394            ENDDO
395            CALL abort
396          ENDIF
397        ENDIF
398       
399        ! In startfi the absolute date is day_ini + time0 + time
400        ! For now on, in the GCM physics, it is day_ini + time0
401        time0 = time(indextime) + time0
402        day_ini = day_ini + INT(time0)
403        time0 = time0 - INT(time0)
404       
405        PRINT*, "phyetat0: Selected time ",time(indextime),
406     .          " at index ",indextime
407     
408        DEALLOCATE(time)
409       
410      ENDIF ! of IF (Time not in file)
411     
412     
413c
[38]414c CO2 ice cover
415c
[999]416      corner(1)=1
417      corner(2)=indextime
418      edges(1)=ngridmx
419      edges(2)=1
[38]420      ierr=nf90_inq_varid(nid,"co2ice",nvarid)
421      IF (ierr.NE.nf90_noerr) THEN
422         PRINT*, 'phyetat0: Le champ <co2ice> est absent'
423         write(*,*)trim(nf90_strerror(ierr))
424         CALL abort
425      ENDIF
[999]426      ierr=nf90_get_var(nid,nvarid,co2ice,corner,edges)
[38]427      IF (ierr.NE.nf90_noerr) THEN
428         PRINT*, 'phyetat0: Lecture echouee pour <co2ice>'
429         write(*,*)trim(nf90_strerror(ierr))
430         CALL abort
431      ENDIF
432      xmin = 1.0E+20
433      xmax = -1.0E+20
434      xmin = MINVAL(co2ice)
435      xmax = MAXVAL(co2ice)
436      PRINT*,'CO2 ice cover <co2ice>:', xmin, xmax
437c
438c Lecture des temperatures du sol:
439c
440      ierr=nf90_inq_varid(nid,"tsurf",nvarid)
441      IF (ierr.NE.nf90_noerr) THEN
442         PRINT*, 'phyetat0: Le champ <tsurf> est absent'
443! Ehouarn: dump this part (causes problems with netcdf90 interface)
444!         PRINT*, '          Mais je vais essayer de lire TS**'
445!         IF (nbsrf.GT.99) THEN
446!            PRINT*, "Trop de sous-mailles"
447!            CALL abort
448!         ENDIF
449!         DO nsrf = 1, nbsrf
450!           WRITE(str2,'(i2.2)') nsrf
451!           ierr=nf90_inq_varid(nid,"TS"//str2,nvarid)
452!           IF (ierr.NE.nf90_noerr) THEN
453!              PRINT*, "phyetat0: Le champ <TS"//str2//"> est absent"
454!              write(*,*)trim(nf90_strerror(ierr))
455!              CALL abort
456!           ENDIF
457!           ierr=nf90_get_var(nid,nvarid,tsurf(1,nsrf))
458!           IF (ierr.NE.nf90_noerr) THEN
459!             PRINT*, "phyetat0: Lecture echouee pour <TS"//str2//">"
460!             write(*,*)trim(nf90_strerror(ierr))
461!             CALL abort
462!           ENDIF
463!           xmin = 1.0E+20
464!           xmax = -1.0E+20
465!           xmin = MINVAL(tsurf)
466!           xmax = MAXVAL(tsurf)
467!           PRINT*,'Temperature du sol TS**:', nsrf, xmin, xmax
468!         ENDDO
469      ELSE
470         PRINT*, 'phyetat0: Le champ <tsurf> est present'
471         PRINT*, '          J ignore donc les autres temperatures TS**'
472         !ierr=nf90_get_var(nid,nvarid,tsurf(1,1))
[999]473         ierr=nf90_get_var(nid,nvarid,tsurf,corner,edges)
[38]474         IF (ierr.NE.nf90_noerr) THEN
475            PRINT*, "phyetat0: Lecture echouee pour <TSURF>"
476            write(*,*)trim(nf90_strerror(ierr))
477            CALL abort
478         ENDIF
479         xmin = 1.0E+20
480         xmax = -1.0E+20
481         xmin = MINVAL(tsurf)
482         xmax = MAXVAL(tsurf)
483         PRINT*,'Temperature du sol <tsurf>', xmin, xmax
484! Ehouarn: Mars => nbsurf=1
485!         IF (nbsrf >= 2) THEN
486!            DO nsrf = 2, nbsrf
487!               DO i = 1, ngridmx
488!                  tsurf(i,nsrf) = tsurf(i,1)
489!               ENDDO
490!            ENDDO
491!         ENDIF
492      ENDIF
493c
494c Lecture des temperatures du sol profond:
495c
496!      IF (nsoil.GT.99 .OR. nbsrf.GT.99) THEN
497!         PRINT*, "Trop de couches ou sous-mailles"
498!         CALL abort
499!      ENDIF
500!      DO nsrf = 1, nbsrf
501!         DO isoil=1, nsoil
502!            WRITE(str7,'(i2.2,"srf",i2.2)') isoil, nsrf
503!            ierr = NF_INQ_VARID (nid, 'tsoil', nvarid)
504!            IF (ierr.NE.nf90_noerr) THEN
505!               PRINT*, "phyetat0: Le champ <tsoil> est absent"
506!               PRINT*, "          Il prend donc la valeur de surface"
507!               DO i=1, ngridmx
508!                  tsoil(i,isoil,nsrf)=tsurf(i,nsrf)
509!               ENDDO
510!            ELSE
511!#ifdef NC_DOUBLE
512!              ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tsoil(1,1,nsrf))
513!#else
514!              ierr = NF_GET_VAR_REAL(nid, nvarid, tsoil(1,1,nsrf))
515!#endif
516!               IF (ierr.NE.nf90_noerr) THEN
517!                  PRINT*, "Lecture echouee pour <tsoil>"
518!                  CALL abort
519!               ENDIF
520!            ENDIF
521!         ENDDO
522!      ENDDO
523!      xmin = 1.0E+20
524!      xmax = -1.0E+20
525!      xmin = MINVAL(tsoil)
526!      xmax = MAXVAL(tsoil)
527!      PRINT*,'Temperatures du sol profond <tsoil>', xmin, xmax
528c
529c Surface emissivity
530c
531      ierr=nf90_inq_varid(nid,"emis",nvarid)
532      IF (ierr.NE.nf90_noerr) THEN
533         PRINT*, 'phyetat0: Le champ <emis> est absent'
534         write(*,*)trim(nf90_strerror(ierr))
535         CALL abort
536      ENDIF
[999]537      ierr=nf90_get_var(nid,nvarid,emis,corner,edges)
[38]538      IF (ierr.NE.nf90_noerr) THEN
539         PRINT*, 'phyetat0: Lecture echouee pour <emis>'
540         write(*,*)trim(nf90_strerror(ierr))
541         CALL abort
542      ENDIF
543      xmin = 1.0E+20
544      xmax = -1.0E+20
545      xmin = MINVAL(emis)
546      xmax = MAXVAL(emis)
547      PRINT*,'Surface emissivity <emis>:', xmin, xmax
[224]548
549!
550! surface roughness length (NB: z0 is a common in surfdat.h)
551!
552      ierr=nf90_inq_varid(nid,"z0",nvarid)
553      if (ierr.ne.nf90_noerr) then
554        write(*,*) 'phyetat0: did not find z0 field!'
555        write(*,*) 'will use constant value of z0_default instead'
556        z0(:)=z0_default
557      else
558        ierr=nf90_get_var(nid,nvarid,z0)
559        if (ierr.ne.nf90_noerr) then
560          write(*,*) 'phyetat0: Failed loading <z0>'
561          write(*,*)trim(nf90_strerror(ierr))
562          call abort
563        endif
564        xmin=1.0E+20
565        xmax=-1.0E+20
566        xmin=minval(z0)
567        xmax=maxval(z0)
568        write(*,*)'Surface roughness <z0>:',xmin,xmax
569      endif
570
571
[38]572c
573c pbl wind variance
574c
[999]575      corner(1)=1
576      corner(2)=1
577      corner(3)=indextime
578      edges(1)=ngridmx
579      edges(2)=llm+1
580      edges(3)=1
[38]581      ierr=nf90_inq_varid(nid,"q2",nvarid)
582      IF (ierr.NE.nf90_noerr) THEN
583         PRINT*, 'phyetat0: Le champ <q2> est absent'
584         write(*,*)trim(nf90_strerror(ierr))
585         CALL abort
586      ENDIF
[999]587      ierr=nf90_get_var(nid,nvarid,q2,corner,edges)
[38]588      IF (ierr.NE.nf90_noerr) THEN
589         PRINT*, 'phyetat0: Lecture echouee pour <q2>'
590         write(*,*)trim(nf90_strerror(ierr))
591         CALL abort
592      ENDIF
593      xmin = 1.0E+20
594      xmax = -1.0E+20
595      xmin = MINVAL(q2)
596      xmax = MAXVAL(q2)
597      PRINT*,'pbl wind variance <q2>:', xmin, xmax
598c
599c tracer on surface
600c
[999]601      corner(1)=1
602      corner(2)=indextime
603      edges(1)=ngridmx
604      edges(2)=1
[38]605      IF(nq.GE.1) THEN
606         nqold=nq
607         DO iq=1,nq
608!            str7(1:5)='qsurf'
609!            WRITE(str7(6:7),'(i2.2)') iq
610!            ierr = NF_INQ_VARID (nid,str7,nvarid)
611           IF (oldtracernames) THEN
612             txt=" "
613             write(txt,'(a5,i2.2)')'qsurf',iq
614           ELSE
615             txt=tnom(iq)
616             if (txt.eq."h2o_vap") then
617               ! There is no surface tracer for h2o_vap;
618               ! "h2o_ice" should be loaded instead
619               txt="h2o_ice"
620               write(*,*) 'phyetat0: loading surface tracer',
621     &                     ' h2o_ice instead of h2o_vap'
622             endif
623           ENDIF ! of IF (oldtracernames) THEN
624           ierr=nf90_inq_varid(nid,txt,nvarid)
625           IF (ierr.NE.nf90_noerr) THEN
626             write(*,*) 'PHYETAT0: WARNING : surface tracer',trim(txt),
627     &                  ' not found in file'
628             write(*,*) trim(txt), ' set to 0'
629             do ig=1,ngridmx
630               qsurf(ig,iq)=0.
631             end do
632             nqold=min(iq-1,nqold)
633           ELSE
634           !ierr=nf90_get_var(nid,nvarid,qsurf(1,iq))
[999]635           ierr=nf90_get_var(nid,nvarid,surffield,corner,edges)
[38]636           qsurf(1:ngridmx,iq)=surffield(1:ngridmx)
637             IF (ierr.NE.nf90_noerr) THEN
638               PRINT*, 'phyetat0: Lecture echouee pour <',trim(txt),'>'
639               write(*,*)trim(nf90_strerror(ierr))
640               CALL abort
641             ENDIF
642           ENDIF
643           xmin = 1.0E+20
644           xmax = -1.0E+20
645           xmin = MINVAL(qsurf(1:ngridmx,iq))
646           xmax = MAXVAL(qsurf(1:ngridmx,iq))
647           PRINT*,'tracer on surface <',trim(txt),'>:',xmin,xmax
648         ENDDO
649         if ((nqold.lt.nq).and.(nqold.ge.1)) then
650c        case when new tracer are added in addition to old ones
651             write(*,*)'qsurf 1 to ', nqold,'were already present'
652             write(*,*)'qsurf ', nqold+1,' to ', nqmx,'are new'
653!            yes=' '
654!            do while ((yes.ne.'y').and.(yes.ne.'n'))
655!             write(*,*) 'Would you like to reindex qsurf # 1 ->',nqold
656!             write(*,*) 'to #',nqmx-nqold+1,'->', nqmx,'   (y or n) ?'
657!             read(*,fmt='(a)') yes
658!            end do
659!            if (yes.eq.'y') then
660!              write(*,*) 'OK, let s reindex qsurf'
661!                 do ig=1,ngridmx
662!                    do iq=nqmx,nqmx-nqold+1,-1
663!                       qsurf(ig,iq)=qsurf(ig,iq-nqmx+nqold)
664!                    end do
665!                    do iq=nqmx-nqold,1,-1
666!                       qsurf(ig,iq)= 0.
667!                    end do
668!                 end do
669!            end if
670         end if
671      ENDIF
672
673! Call to soil_settings, in order to read soil temperatures,
674! as well as thermal inertia and volumetric heat capacity
675
[999]676      call soil_settings(nid,ngridmx,nsoil,tsurf,tsoil,indextime)
[38]677c
678c Fermer le fichier:
679c
680      ierr=nf90_close(nid)
681c
682      RETURN
683      END
Note: See TracBrowser for help on using the repository browser.