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

Last change on this file since 2137 was 1918, checked in by emillour, 8 years ago

Mars GCM:
Code cleanup:

  • remove "comorbit.h" since it is no longer used.
  • turn "datafile.h" into module datafile_mod.F90 (and rename variable "datafile" as "datadir" since it stores the path to the datafile directory).

EM
File size: 9.7 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      USE comconst_mod, ONLY: g,pi
47      use datafile_mod, only: datadir
48
49      implicit none
50
51      include "dimensions.h"
52      include "paramet.h"
53      include "comgeom.h"
54      include "netcdf.inc"
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      datadir="/u/lmdz/WWW/planets/mars/datadir" ! default path to surface.nc
116      call getin("datadir",datadir) ! but users may specify another path
117     
118      ierr = NF_OPEN (trim(datadir)//'/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(datadir),'/'
124        write(*,*)'1) You can set this path in the
125     & callphys.def file:'
126        write(*,*)'   datadir=/path/to/the/datafiles'
127        write(*,*)'2) If necessary, surface.nc (and other datafiles)'
128        write(*,*)'   can be obtained online on:'
129        write(*,*)'http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir'
130        CALL ABORT
131      ENDIF
132
133c
134c Lecture des latitudes (coordonnees):
135c
136      ierr = NF_INQ_VARID (unit, "latitude", nvarid)
137#ifdef NC_DOUBLE
138      ierr = NF_GET_VAR_DOUBLE(unit, nvarid, latitude)
139#else
140      ierr = NF_GET_VAR_REAL(unit, nvarid, latitude)
141#endif
142c
143c Lecture des longitudes (coordonnees):
144c
145      ierr = NF_INQ_VARID (unit, "longitude", nvarid)
146#ifdef NC_DOUBLE
147      ierr = NF_GET_VAR_DOUBLE(unit, nvarid, longitude)
148#else
149      ierr = NF_GET_VAR_REAL(unit, nvarid, longitude)
150#endif
151
152c-----------------------------------------------------------------------
153c    Passage au format boites scalaires
154c-----------------------------------------------------------------------
155
156c-----------------------------------------------------------------------
157c       longitude(imd)        -->      rlonud(imdp1)
158c-----------------------------------------------------------------------
159
160c Passage en coordonnees boites scalaires et en radian
161      do i=1,imd
162          rlonud(i)=(longitude(i)+.5)*pi/180.
163      enddo
164
165c Repetition de la valeur im+1
166      rlonud(imdp1)=rlonud(1) + 2*pi
167
168c-----------------------------------------------------------------------
169c        latitude(jmdp1)         -->        rlonvd(jmd)
170c-----------------------------------------------------------------------
171
172c Passage en coordonnees boites scalaires et en radian
173      do j=1,jmd
174          rlatvd(j)=(latitude(j)-.5)*pi/180.
175      enddo
176
177c=======================================================================
178c   lecture NetCDF de albedo, thermal, relief, zdtm (pour francois Lott)
179c=======================================================================
180
181      string(0) = 'z0'
182      string(1) = 'albedo'
183      string(2) = 'thermal'
184      if (relief.ne.'pla') then
185        write(*,*) ' MOLA topography'
186        relief = 'MOL'
187          string(3) = 'z'//relief
188      else
189          string(3) = 'zMOL'  ! pour qu''il lise qqchose sur le fichier
190                            ! remise a 0 derriere
191      endif
192      string(4) = 'zMOL'    ! lecture pour calcul topog. sous-maille
193 
194
195      DO k=0,4
196          write(*,*) 'string',k,string(k)
197         
198c-----------------------------------------------------------------------
199c    initialisation
200c-----------------------------------------------------------------------
201      call initial0(iimp1*jjp1,pfield)
202      call initial0(imd*jmdp1,zdata)
203      call initial0(imdp1*jmdp1,zdataS)
204
205c-----------------------------------------------------------------------
206c    Lecture NetCDF 
207c-----------------------------------------------------------------------
208
209      ierr = NF_INQ_VARID (unit, string(k), nvarid)
210      if (ierr.ne.nf_noerr) then
211        write(*,*) 'datareadnc error, cannot find ',trim(string(k))
212        write(*,*) ' in file ',trim(datadir),'/surface.nc'
213        stop
214      endif
215#ifdef NC_DOUBLE
216      ierr = NF_GET_VAR_DOUBLE(unit, nvarid, zdata)
217#else
218      ierr = NF_GET_VAR_REAL(unit, nvarid, zdata)
219#endif
220      if (ierr.ne.nf_noerr) then
221        write(*,*) 'datareadnc: error failed loading ',trim(string(k))
222        stop
223      endif
224
225c-----------------------------------------------------------------------
226c        Cas particulier "Francois Lott" ( k=4 ) (relief sous-maille)
227c-----------------------------------------------------------------------
228      if (k.eq.4) then
229
230          zdata(:)=1000.*zdata(:)
231          longitude(:)=(pi/180.)*longitude(:)
232          latitude(:)=(pi/180.)*latitude(:)
233
234          call grid_noro1(360, 180, longitude, latitude, zdata,
235     .         iim, jjp1, rlonv, rlatu, zmea,zstd,zsig,zgam,zthe)
236
237      endif
238
239c-----------------------------------------------------------------------
240c   Passage de zdata en grille (imdp1 jmdp1)
241c-----------------------------------------------------------------------
242      do j=1,jmdp1
243          do i=1,imd
244              zdataS(i+imdp1*(j-1)) = zdata(i+ngridmxgdat*(j-1))
245          enddo
246          zdataS(imdp1+imdp1*(j-1)) = zdata(1+ngridmxgdat*(j-1))
247      enddo
248
249c-----------------------------------------------------------------------
250c    Interpolation
251c-----------------------------------------------------------------------
252      call interp_horiz(zdataS,pfield,imd,jmd,
253     .    iim, jjm,1,rlonud,rlatvd,rlonu,rlatv)
254
255c-----------------------------------------------------------------------
256c    Periodicite   
257c-----------------------------------------------------------------------
258
259      do j=1,jjp1
260         pfield(iimp1*j) =  pfield(1+iimp1*(j-1))
261      enddo
262 
263c-----------------------------------------------------------------------
264c    Sauvegarde des champs   
265c-----------------------------------------------------------------------
266
267      if (k.eq.0) then                    ! z0
268         z0(1:iimp1*jjp1)=pfield(1:iimp1*jjp1)*.01
269         ! multiplied by 0.01 to have z0 in m
270      elseif (k.eq.1) then                    ! albedo
271         do i=1,iimp1*jjp1
272              alb(i) = pfield(i)
273          enddo
274      elseif (k.eq.2) then                ! thermal
275         do i=1,iimp1*jjp1
276              ith(i) = pfield(i)
277          enddo
278      elseif (k.eq.3) then                ! relief
279        if (relief.eq.'pla') then
280              call initial0(iimp1*jjp1,phisinit)
281        else
282             do i=1,iimp1*jjp1
283                  phisinit(i) = pfield(i)
284              enddo
285        endif
286      endif
287
288      ENDDO
289
290c-----------------------------------------------------------------------
291c    Traitement Phisinit
292c-----------------------------------------------------------------------
293
294      phisinit(1:iimp1*jjp1)=1000.*phisinit(1:iimp1*jjp1)
295      phisinit(:)=g*phisinit(:)
296
297c-----------------------------------------------------------------------
298c    FIN
299c-----------------------------------------------------------------------
300
301      END
Note: See TracBrowser for help on using the repository browser.