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

Last change on this file since 706 was 690, checked in by acolaitis, 13 years ago

Code re-organization in diverse parts of the GCM code. These are NOT cosmetic changes, but are needed for compilation of the Mesoscale model in NESTED configuration
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
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/~forget/datagcm/datafile'
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          call multscal(imd*jmdp1,zdata,1000.,zdata)
231          call multscal(imd,longitude,pi/180.,longitude)
232          call multscal(jmdp1,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          CALL dump2d(iip1,jjp1,zmea,'zmea')
238          CALL dump2d(iip1,jjp1,zstd,'zstd')
239          CALL dump2d(iip1,jjp1,zsig,'zsig')
240          CALL dump2d(iip1,jjp1,zgam,'zgam')
241          CALL dump2d(iip1,jjp1,zthe,'zthe')
242
243      endif
244
245c-----------------------------------------------------------------------
246c   Passage de zdata en grille (imdp1 jmdp1)
247c-----------------------------------------------------------------------
248      do j=1,jmdp1
249          do i=1,imd
250              zdataS(i+imdp1*(j-1)) = zdata(i+ngridmxgdat*(j-1))
251          enddo
252          zdataS(imdp1+imdp1*(j-1)) = zdata(1+ngridmxgdat*(j-1))
253      enddo
254
255c-----------------------------------------------------------------------
256c    Interpolation
257c-----------------------------------------------------------------------
258      call interp_horiz(zdataS,pfield,imd,jmd,
259     .    iim, jjm,1,rlonud,rlatvd,rlonu,rlatv)
260
261c-----------------------------------------------------------------------
262c    Periodicite   
263c-----------------------------------------------------------------------
264
265      do j=1,jjp1
266         pfield(iimp1*j) =  pfield(1+iimp1*(j-1))
267      enddo
268 
269c-----------------------------------------------------------------------
270c    Sauvegarde des champs   
271c-----------------------------------------------------------------------
272
273      if (k.eq.0) then                    ! z0
274         z0(1:iimp1*jjp1)=pfield(1:iimp1*jjp1)*.01
275         ! multiplied by 0.01 to have z0 in m
276         CALL dump2d(iimp1,jjp1,z0,'z0 in m')
277      elseif (k.eq.1) then                    ! albedo
278         do i=1,iimp1*jjp1
279              alb(i) = pfield(i)
280          enddo
281      elseif (k.eq.2) then                ! thermal
282         do i=1,iimp1*jjp1
283              ith(i) = pfield(i)
284          enddo
285      elseif (k.eq.3) then                ! relief
286        if (relief.eq.'pla') then
287              call initial0(iimp1*jjp1,phisinit)
288        else
289             do i=1,iimp1*jjp1
290                  phisinit(i) = pfield(i)
291              enddo
292        endif
293      endif
294
295      ENDDO
296
297c-----------------------------------------------------------------------
298c    Traitement Phisinit
299c-----------------------------------------------------------------------
300
301      DO i=1,iimp1*jjp1
302            phisinit(i)=1000.*phisinit(i)
303      ENDDO
304      CALL dump2d(iimp1,jjp1,phisinit,'Altitude in m')
305      CALL multscal(iimp1*jjp1,phisinit,g,phisinit)
306
307c-----------------------------------------------------------------------
308c    FIN
309c-----------------------------------------------------------------------
310
311      END
Note: See TracBrowser for help on using the repository browser.