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

Last change on this file since 1564 was 1422, checked in by milmd, 10 years ago

In GENERIC, MARS and COMMON models replace some include files by modules (usefull for decoupling physics with dynamics).
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
48      implicit none
49
50#include "dimensions.h"
51#include "paramet.h"
52#include "comgeom.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/lmdz/WWW/planets/mars/datadir" ! 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
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(datafile),'/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.