source: trunk/LMDZ.GENERIC/libf/phystd/datareadnc.F @ 832

Last change on this file since 832 was 787, checked in by aslmd, 12 years ago

LMDZ.GENERIC. (Sorry for long text but this is a quite major commit)

Paving the path for parallel computations. And moving towards a more flexible code.

Automatic allocation is used within all routines in phystd. No further mention to ngridmx and nqmx.

  1. ngridmx and nqmx are still used in LMDZ.GENERIC in the dyn3d part
  2. if the LMDZ4/LMDZ5 dynamical core is used, there is no more fixed dimensions ngridmx and nqmx --> a fully flexible parallel implementation is now possible (e.g. no need to recompile when changing numbers of processors)

The important stuff :

  • Compilation checked with ifort. OK with and without debug mode. No errors. Checked for: gcm, newstart, rcm1d, kcm1d
  • RUN GCM: Running an Earth test case. Comparison with previous revision --> debug mode : perfect match. bit by bit (diff command). checked with plots --> O1 mode : close match (checked with plots) --> O2 mode : sometimes up to 0.5 K departure.... BUT in this new version O2 and O1 are quite close while in previous version O1 and O2 differed by about, well, typically 0.5 K (pictures available on request)
  • RUN NEWSTART : perfect match (bit-by-bit) in either debug or normal mode.
  • RUN RCM1D : perfect match in normal mode.
  • RUN KCM1D : not tested (I don't know what is the use of kcm1d)

List of main changes :

  • Additional arguments to some subroutines (ngrid and nq)
  • F77 include strategy is obsolete and replaced by F90 module strategy In this new strategy arrays are allocatable and allocated once at first use This has to be done for all common featuring arrays defined with ngridmx or nqmx

surfdat.h >> surfdat_h.F90
tracer.h >> tracer_h.F90
comsaison.h >> comsaison_h.F90
comgeomfi.h >> comgeomfi_h.F90
comsoil.h >> comsoil_h.F90
comdiurn.h >> comdiurn_h.F90
fisice.h >> DELETED. was not used. probably a fossil.
watercap.h >> DELETED. variable put in surfdat_h.F90

  • F77 'save' strategy is obsolete and replaced by F90 'allocatable save' strategy (see previous point and e.g. new version of physiq.F90)
  • Suppressing any mention to advtrac.h which is a common in the dynamics and needs nqmx This was easily solved by adding an argument with tracer names, coming from the dynamics This is probably not a definitive solution, ... but this allows for generic physics to work easily with either LMDZ.GENERIC or LMDZ dynamical cores
  • Removing consistency tests between nq and nqmx ; and ngrid and ngridmx. No use now!
  • Adaptation of rcm1d, kcm1d, newstart given above-mentioned changes

A note on phyetat0 and soil_setting:

  • Now written so that a slice of horizontal size 'ngrid' starting at grid point 'cursor' is read in startfi.nc 'cursor' is defined in dimphys.h and initialized by inifis (or in newstart) this is useful for parallel computations. default behavior is the usual one : sequential runs, cursor is 1, size ngrid is the whole global domain

A note on an additional change :

  • nueffrad is now an argument to callcorrk as is the case for reffrad both are saved in physiq this is for consistency and lisibility (previously nueffrad was saved in callcorrk) ... but there is a call to a function which modifies nueffrad in physiq ... previously this was not modifying nueffrad (although it was quite cumbersome to detect this) ... to be conservative I kept this behaviour and highlighted it with an array nueffrad_dummy ... I added a comment because someone might want to change this
File size: 9.7 KB
Line 
1c=======================================================================
2      SUBROUTINE datareadnc(relief,phisinit,alb,ith,
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      use datafile_mod, only: datadir
45! to use  'getin'
46      USE ioipsl_getincom
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
55c=======================================================================
56c   Declarations:
57C=======================================================================
58
59      INTEGER    imd,jmd,imdp1,jmdp1
60      parameter    (imd=360,jmd=179,imdp1=361,jmdp1=180)
61
62      INTEGER    iimp1
63      parameter    (iimp1=iim+1-1/iim)
64
65      CHARACTER relief*3
66
67      REAL        zdata(imd*jmdp1)
68      REAL        zdataS(imdp1*jmdp1)
69      REAL        pfield(iimp1*jjp1)
70
71      REAL        alb(iimp1*jjp1)
72      REAL        ith(iimp1*jjp1)
73      REAL        phisinit(iimp1*jjp1)
74
75      REAL        zmea(imdp1*jmdp1)
76      REAL        zstd(imdp1*jmdp1)
77      REAL        zsig(imdp1*jmdp1)
78      REAL        zgam(imdp1*jmdp1)
79      REAL        zthe(imdp1*jmdp1)
80
81      INTEGER   ierr
82
83      INTEGER   unit,nvarid
84
85      INTEGER    i,j,k
86
87      INTEGER klatdat,ngridmixgdat
88      PARAMETER (klatdat=180,ngridmixgdat=360)
89
90c    on passe une grille en rlonu rlatv et im+1 jm a interp_horiz)
91
92      REAL longitude(imd),latitude(jmdp1) ! Pour lecture des donnees
93      REAL rlonud(imdp1),rlatvd(jmd)
94
95      CHARACTER*20 string
96      DIMENSION string(4)
97
98      CHARACTER*20 filename
99      CHARACTER*50 filestring
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!    First get the corret value of dtadir, if not already done:
114      ! default 'datadir' is set in "datadir_mod"
115      call getin("datadir",datadir)
116      write(*,*) 'Choice of surface data is:'
117      filestring='ls '//trim(datadir)//'/*.nc'
118      call system(filestring)
119
120      write(*,*) ''
121      write(*,*) 'Please choose the relevant file:'
122      write(*,*) '(e.g. type "surface_mars.nc")'
123      read(*,fmt='(a20)') filename
124
125      ierr = NF_OPEN (trim(datadir)//'/'//trim(adjustl(filename)),
126     &  NF_NOWRITE,unit)
127      IF (ierr.NE.NF_NOERR) THEN
128        write(*,*)'Error : cannot open file '//trim(filename)
129        write(*,*)'(in phystd/datareadnc.F)'
130        write(*,*)'It should be in :',trim(datadir),'/'
131        write(*,*)'Check that your path to datagcm:',trim(datadir)
132        write(*,*)' is correct. You can change it in callphys.def with:'
133        write(*,*)' datadir = /absolute/path/to/datagcm'
134        write(*,*)'If necessary surface.nc (and other datafiles)'
135        write(*,*)' can be obtained online on:'
136        write(*,*)' http://www.lmd.jussieu.fr/~forget/datagcm/datafile'
137        CALL ABORT
138      ENDIF
139
140c
141c Lecture des latitudes (coordonnees):
142c
143      ierr = NF_INQ_VARID (unit, "latitude", nvarid)
144#ifdef NC_DOUBLE
145      ierr = NF_GET_VAR_DOUBLE(unit, nvarid, latitude)
146#else
147      ierr = NF_GET_VAR_REAL(unit, nvarid, latitude)
148#endif
149c
150c Lecture des longitudes (coordonnees):
151c
152      ierr = NF_INQ_VARID (unit, "longitude", nvarid)
153#ifdef NC_DOUBLE
154      ierr = NF_GET_VAR_DOUBLE(unit, nvarid, longitude)
155#else
156      ierr = NF_GET_VAR_REAL(unit, nvarid, longitude)
157#endif
158
159c-----------------------------------------------------------------------
160c    Passage au format boites scalaires
161c-----------------------------------------------------------------------
162
163c-----------------------------------------------------------------------
164c       longitude(imd)        -->      rlonud(imdp1)
165c-----------------------------------------------------------------------
166
167c Passage en coordonnees boites scalaires et en radian
168      do i=1,imd
169          rlonud(i)=(longitude(i)+.5)*pi/180.
170      enddo
171
172c Repetition de la valeur im+1
173      rlonud(imdp1)=rlonud(1) + 2*pi
174
175c-----------------------------------------------------------------------
176c        latitude(jmdp1)         -->        rlonvd(jmd)
177c-----------------------------------------------------------------------
178
179c Passage en coordonnees boites scalaires et en radian
180      do j=1,jmd
181          rlatvd(j)=(latitude(j)-.5)*pi/180.
182      enddo
183
184c=======================================================================
185c   lecture NetCDF de albedo, thermal, relief, zdtm (pour francois Lott)
186c=======================================================================
187
188      string(1) = 'albedo'
189      string(2) = 'thermal'
190      if (relief.ne.'pla') then
191        write(*,*) ' La topographie est celle de MOLA'
192        relief = 'MOL'
193          string(3) = 'z'//relief
194      else
195          string(3) = 'zMOL'  ! pour qu''il lise qqchose sur le fichier
196                            ! remise a 0 derriere
197      endif
198      string(4) = 'zMOL'    ! lecture pour calcul topog. sous-maille
199 
200
201      DO k=1,4
202          write(*,*) 'string',k,string(k)
203         
204c-----------------------------------------------------------------------
205c    initialisation
206c-----------------------------------------------------------------------
207      call initial0(iimp1*jjp1,pfield)
208      call initial0(imd*jmdp1,zdata)
209      call initial0(imdp1*jmdp1,zdataS)
210
211c-----------------------------------------------------------------------
212c    Lecture NetCDF 
213c-----------------------------------------------------------------------
214
215      ierr = NF_INQ_VARID (unit, string(k), nvarid)
216#ifdef NC_DOUBLE
217      ierr = NF_GET_VAR_DOUBLE(unit, nvarid, zdata)
218#else
219      ierr = NF_GET_VAR_REAL(unit, nvarid, zdata)
220#endif
221
222c-----------------------------------------------------------------------
223c        Cas particulier "Francois Lott" ( k=4 ) (relief sous-maille)
224c-----------------------------------------------------------------------
225      if (k.eq.4) then
226
227          call multscal(imd*jmdp1,zdata,1000.,zdata)
228          call multscal(imd,longitude,pi/180.,longitude)
229          call multscal(jmdp1,latitude,pi/180.,latitude)
230
231          call grid_noro1(360, 180, longitude, latitude, zdata,
232     .         iim, jjp1, rlonv, rlatu, zmea,zstd,zsig,zgam,zthe)
233
234          CALL dump2d(iip1,jjp1,zmea,'zmea')
235          CALL dump2d(iip1,jjp1,zstd,'zstd')
236          CALL dump2d(iip1,jjp1,zsig,'zsig')
237          CALL dump2d(iip1,jjp1,zgam,'zgam')
238          CALL dump2d(iip1,jjp1,zthe,'zthe')
239
240      endif
241
242c-----------------------------------------------------------------------
243c   Passage de zdata en grille (imdp1 jmdp1)
244c-----------------------------------------------------------------------
245      do j=1,jmdp1
246          do i=1,imd
247              zdataS(i+imdp1*(j-1)) = zdata(i+ngridmixgdat*(j-1))
248          enddo
249          zdataS(imdp1+imdp1*(j-1)) = zdata(1+ngridmixgdat*(j-1))
250      enddo
251
252c-----------------------------------------------------------------------
253c    Interpolation
254c-----------------------------------------------------------------------
255      call interp_horiz(zdataS,pfield,imd,jmd,
256     .    iim, jjm,1,rlonud,rlatvd,rlonu,rlatv)
257
258c-----------------------------------------------------------------------
259c    Periodicite   
260c-----------------------------------------------------------------------
261
262      do j=1,jjp1
263         pfield(iimp1*j) =  pfield(1+iimp1*(j-1))
264      enddo
265 
266c-----------------------------------------------------------------------
267c    Sauvegarde des champs   
268c-----------------------------------------------------------------------
269
270      if (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      DO i=1,iimp1*jjp1
295            phisinit(i)=1000.*phisinit(i)
296      ENDDO
297      CALL dump2d(iimp1,jjp1,phisinit,'Altitude en m')
298      CALL multscal(iimp1*jjp1,phisinit,g,phisinit)
299
300c-----------------------------------------------------------------------
301c    FIN
302c-----------------------------------------------------------------------
303
304      END
Note: See TracBrowser for help on using the repository browser.