source: trunk/LMDZ.MARS/libf/phymars/datareadnc.F @ 224

Last change on this file since 224 was 224, checked in by emillour, 14 years ago

Mars GCM:

Implemented using 'z0' roughness length map (important: 'z0' reference

field is in datafile surface.nc, which has also been updated).

  • made z0 a z0(ngridmx) array and moved 'z0' from 'planete.h' to 'surfdat.h'; added a 'z0_default' (common in surfdat.h) corresponding to the 'control' array value (contole(19) in startfi.nc).
  • adapted 'tabfi.F' to use 'z0_default'.
  • adapted 'phyetat0.F' to look for a 'z0' field in startfi.nc. If 'z0' is not found in the startfi.nc file, then the uniform default value (z0_default) is used.
  • modified 'physdem1.F' to write 'z0' field to restart.nc
  • adapted use of z0() in 'physiq.F' (diagnostic computation of surface stress), 'vdifc.F' and 'vdif_cd.F'.
  • adapted 'dustdevil.F' to use 'z0_default'.
  • 'testphys1d.F' now uses 'z0_default', and the value to use can be set in run.def (with "z0=TheValueYouWant?").
  • modified 'datareadnc.F' to load reference map of 'z0' from surface.nc, and added a 'z0' option in 'newstart.F' to force a uniform value of z0. Note that the use of the z0 map is automatic when using newstart, but only when it loads a start_archive.nc file.
File size: 10.0 KB
Line 
1c=======================================================================
2      SUBROUTINE datareadnc(relief,phisinit,alb,ith,z0,
3     &                    zmea,zstd,zsig,zgam,zthe)
4c=======================================================================
5c
6c
7c   Author: F. Hourdin      01/1997
8c   -------
9c
10c   Object: To read data from Martian surface to use in a GCM
11c   ------                from NetCDF file "surface.nc"
12c
13c
14c   Arguments:
15c   ----------
16c
17c     Inputs:
18c     ------
19c
20c     Outputs:
21c     --------
22c
23c=======================================================================
24c   donnees ALBEDO, INERTIE THERMIQUE, RELIEF:
25c
26c       Ces donnees sont au format NetCDF dans le fichier "surface.nc"
27c
28c   360 valeurs en longitude (de -179.5 a 179.5)
29c   180 valeurs en latitudes (de 89.5 a -89.5)
30c
31c   Pour les passer au format de la grille, on utilise "interp_horiz.F"
32c
33c   Il faut donc que ces donnees soient au format grille scalaire
34c               (imold+1 jmold+1)
35c       avec passage des coordonnees de la "boite" (rlonu, rlatv)
36c
37c   On prend imd (d pour donnees!)
38c           imd = 360 avec copie de la 1ere valeur sur la imd+1
39c                   (rlonud de -179 a -181)
40c           jmd = 179
41c                   (rlatvd de 89 a -89)
42c=======================================================================
43
44! to use  'getin'
45       use ioipsl_getincom
46
47      implicit none
48
49#include "dimensions.h"
50#include "paramet.h"
51#include "comgeom.h"
52#include "comconst.h"
53#include "netcdf.inc"
54#include "datafile.h"
55
56c=======================================================================
57c   Declarations:
58C=======================================================================
59
60      INTEGER    imd,jmd,imdp1,jmdp1
61      parameter    (imd=360,jmd=179,imdp1=361,jmdp1=180)
62
63      INTEGER    iimp1
64      parameter    (iimp1=iim+1-1/iim)
65
66! Arguments:
67      CHARACTER(len=3),intent(inout) :: relief
68      REAL,intent(out) :: phisinit(iimp1*jjp1)
69      REAL,intent(out) :: alb(iimp1*jjp1)
70      REAL,intent(out) :: ith(iimp1*jjp1)
71      REAL,intent(out) :: z0(iimp1*jjp1)
72      REAL,intent(out) :: zmea(imdp1*jmdp1)
73      REAL,intent(out) :: zstd(imdp1*jmdp1)
74      REAL,intent(out) :: zsig(imdp1*jmdp1)
75      REAL,intent(out) :: zgam(imdp1*jmdp1)
76      REAL,intent(out) :: zthe(imdp1*jmdp1)
77
78! Local variables:
79      REAL        zdata(imd*jmdp1)
80      REAL        zdataS(imdp1*jmdp1)
81      REAL        pfield(iimp1*jjp1)
82
83      INTEGER     ierr
84
85      INTEGER   unit,nvarid
86
87      INTEGER    i,j,k
88
89      INTEGER klatdat,ngridmxgdat
90      PARAMETER (klatdat=180,ngridmxgdat=360)
91
92c    on passe une grille en rlonu rlatv et im+1 jm a interp_horiz)
93
94      REAL longitude(imd),latitude(jmdp1) ! Pour lecture des donnees
95      REAL rlonud(imdp1),rlatvd(jmd)
96
97      CHARACTER*20 string
98      DIMENSION string(0:4)
99
100
101#include "lmdstd.h"
102#include "fxyprim.h"
103
104      pi=2.*ASIN(1.)
105
106c=======================================================================
107c    rlonud, rlatvd
108c=======================================================================
109
110c-----------------------------------------------------------------------
111c    Lecture NetCDF des donnees latitude et longitude
112c-----------------------------------------------------------------------
113      write(*,*) 'datareadnc: opening file surface.nc'
114
115      datafile="/u/forget/WWW/datagcm/datafile" ! default path to surface.nc
116      call getin("datadir",datafile) ! but users may specify another path
117     
118      ierr = NF_OPEN (trim(datafile)//'/surface.nc',
119     &  NF_NOWRITE,unit)
120      IF (ierr.NE.NF_NOERR) THEN
121        write(*,*)'Error : cannot open file surface.nc '
122        write(*,*)'(in phymars/datareadnc.F)'
123        write(*,*)'It should be in :',trim(datafile),'/'
124        write(*,*)'1) You can set this path in the callphys.def file:'
125        write(*,*)'   datadir=/path/to/the/datafiles'
126        write(*,*)'2) If necessary, surface.nc (and other datafiles)'
127        write(*,*)'   can be obtained online on:'
128        write(*,*)' http://www.lmd.jussieu.fr/~forget/datagcm/datafile'
129        CALL ABORT
130      ENDIF
131
132c
133c Lecture des latitudes (coordonnees):
134c
135      ierr = NF_INQ_VARID (unit, "latitude", nvarid)
136#ifdef NC_DOUBLE
137      ierr = NF_GET_VAR_DOUBLE(unit, nvarid, latitude)
138#else
139      ierr = NF_GET_VAR_REAL(unit, nvarid, latitude)
140#endif
141c
142c Lecture des longitudes (coordonnees):
143c
144      ierr = NF_INQ_VARID (unit, "longitude", nvarid)
145#ifdef NC_DOUBLE
146      ierr = NF_GET_VAR_DOUBLE(unit, nvarid, longitude)
147#else
148      ierr = NF_GET_VAR_REAL(unit, nvarid, longitude)
149#endif
150
151c-----------------------------------------------------------------------
152c    Passage au format boites scalaires
153c-----------------------------------------------------------------------
154
155c-----------------------------------------------------------------------
156c       longitude(imd)        -->      rlonud(imdp1)
157c-----------------------------------------------------------------------
158
159c Passage en coordonnees boites scalaires et en radian
160      do i=1,imd
161          rlonud(i)=(longitude(i)+.5)*pi/180.
162      enddo
163
164c Repetition de la valeur im+1
165      rlonud(imdp1)=rlonud(1) + 2*pi
166
167c-----------------------------------------------------------------------
168c        latitude(jmdp1)         -->        rlonvd(jmd)
169c-----------------------------------------------------------------------
170
171c Passage en coordonnees boites scalaires et en radian
172      do j=1,jmd
173          rlatvd(j)=(latitude(j)-.5)*pi/180.
174      enddo
175
176c=======================================================================
177c   lecture NetCDF de albedo, thermal, relief, zdtm (pour francois Lott)
178c=======================================================================
179
180      string(0) = 'z0'
181      string(1) = 'albedo'
182      string(2) = 'thermal'
183      if (relief.ne.'pla') then
184        write(*,*) ' MOLA topography'
185        relief = 'MOL'
186          string(3) = 'z'//relief
187      else
188          string(3) = 'zMOL'  ! pour qu''il lise qqchose sur le fichier
189                            ! remise a 0 derriere
190      endif
191      string(4) = 'zMOL'    ! lecture pour calcul topog. sous-maille
192 
193
194      DO k=0,4
195          write(*,*) 'string',k,string(k)
196         
197c-----------------------------------------------------------------------
198c    initialisation
199c-----------------------------------------------------------------------
200      call initial0(iimp1*jjp1,pfield)
201      call initial0(imd*jmdp1,zdata)
202      call initial0(imdp1*jmdp1,zdataS)
203
204c-----------------------------------------------------------------------
205c    Lecture NetCDF 
206c-----------------------------------------------------------------------
207
208      ierr = NF_INQ_VARID (unit, string(k), nvarid)
209      if (ierr.ne.nf_noerr) then
210        write(*,*) 'datareadnc error, cannot find ',trim(string(k))
211        write(*,*) ' in file ',trim(datafile),'/surface.nc'
212        stop
213      endif
214#ifdef NC_DOUBLE
215      ierr = NF_GET_VAR_DOUBLE(unit, nvarid, zdata)
216#else
217      ierr = NF_GET_VAR_REAL(unit, nvarid, zdata)
218#endif
219      if (ierr.ne.nf_noerr) then
220        write(*,*) 'datareadnc: error failed loading ',trim(string(k))
221        stop
222      endif
223
224c-----------------------------------------------------------------------
225c        Cas particulier "Francois Lott" ( k=4 ) (relief sous-maille)
226c-----------------------------------------------------------------------
227      if (k.eq.4) then
228
229          call multscal(imd*jmdp1,zdata,1000.,zdata)
230          call multscal(imd,longitude,pi/180.,longitude)
231          call multscal(jmdp1,latitude,pi/180.,latitude)
232
233          call grid_noro1(360, 180, longitude, latitude, zdata,
234     .         iim, jjp1, rlonv, rlatu, zmea,zstd,zsig,zgam,zthe)
235
236          CALL dump2d(iip1,jjp1,zmea,'zmea')
237          CALL dump2d(iip1,jjp1,zstd,'zstd')
238          CALL dump2d(iip1,jjp1,zsig,'zsig')
239          CALL dump2d(iip1,jjp1,zgam,'zgam')
240          CALL dump2d(iip1,jjp1,zthe,'zthe')
241
242      endif
243
244c-----------------------------------------------------------------------
245c   Passage de zdata en grille (imdp1 jmdp1)
246c-----------------------------------------------------------------------
247      do j=1,jmdp1
248          do i=1,imd
249              zdataS(i+imdp1*(j-1)) = zdata(i+ngridmxgdat*(j-1))
250          enddo
251          zdataS(imdp1+imdp1*(j-1)) = zdata(1+ngridmxgdat*(j-1))
252      enddo
253
254c-----------------------------------------------------------------------
255c    Interpolation
256c-----------------------------------------------------------------------
257      call interp_horiz(zdataS,pfield,imd,jmd,
258     .    iim, jjm,1,rlonud,rlatvd,rlonu,rlatv)
259
260c-----------------------------------------------------------------------
261c    Periodicite   
262c-----------------------------------------------------------------------
263
264      do j=1,jjp1
265         pfield(iimp1*j) =  pfield(1+iimp1*(j-1))
266      enddo
267 
268c-----------------------------------------------------------------------
269c    Sauvegarde des champs   
270c-----------------------------------------------------------------------
271
272      if (k.eq.0) then                    ! z0
273         z0(1:iimp1*jjp1)=pfield(1:iimp1*jjp1)*.01
274         ! multiplied by 0.01 to have z0 in m
275         CALL dump2d(iimp1,jjp1,z0,'z0 in m')
276      elseif (k.eq.1) then                    ! albedo
277         do i=1,iimp1*jjp1
278              alb(i) = pfield(i)
279          enddo
280      elseif (k.eq.2) then                ! thermal
281         do i=1,iimp1*jjp1
282              ith(i) = pfield(i)
283          enddo
284      elseif (k.eq.3) then                ! relief
285        if (relief.eq.'pla') then
286              call initial0(iimp1*jjp1,phisinit)
287        else
288             do i=1,iimp1*jjp1
289                  phisinit(i) = pfield(i)
290              enddo
291        endif
292      endif
293
294      ENDDO
295
296c-----------------------------------------------------------------------
297c    Traitement Phisinit
298c-----------------------------------------------------------------------
299
300      DO i=1,iimp1*jjp1
301            phisinit(i)=1000.*phisinit(i)
302      ENDDO
303      CALL dump2d(iimp1,jjp1,phisinit,'Altitude in m')
304      CALL multscal(iimp1*jjp1,phisinit,g,phisinit)
305
306c-----------------------------------------------------------------------
307c    FIN
308c-----------------------------------------------------------------------
309
310      END
Note: See TracBrowser for help on using the repository browser.