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

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

Mars GCM:

  • IMPORTANT CHANGE: Removed all reference/use of ngridmx (dimphys.h) in routines (necessary prerequisite to using parallel dynamics); in most cases this just means adding 'ngrid' as routine argument, and making local saved variables allocatable (and allocated at first call). In the process, had to convert many *.h files to equivalent modules: yomaer.h => yomaer_h.F90 , surfdat.h => surfdat_h.F90 , comsaison.h => comsaison_h.F90 , yomlw.h => yomlw_h.F90 , comdiurn.h => comdiurn_h.F90 , dimradmars.h => dimradmars_mod.F90 , comgeomfi.h => comgeomfi_h.F90, comsoil.h => comsoil_h.F90 , slope.h => slope_mod.F90
  • Also updated EOF routines, everything is now in eofdump_mod.F90
  • Removed unused routine lectfux.F (in dyn3d)

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