source: trunk/LMDZ.MARS/util/hrecast.F90 @ 1009

Last change on this file since 1009 was 860, checked in by tnavarro, 12 years ago

Modification to create large files, that is to say with at least 2 variables > 2GiB. In that case, it is necessary to have an unlimited dimension.
File size: 48.0 KB
Line 
1
2
3program hrecast
4
5! This program reads fields from GCM output files
6! (diagfi.nc, stats.nc, concat.nc, ...)
7! and recasts it on a new horizontal grid specified by the user.
8! The output file name is automatically generated from input file name:
9! if input file name is 'input.nc', then output file will be 'input_h.nc'
10!
11! NB: It is OK if the output grid is not like typical GCM outputs grids,
12!     but longitudes range must be something like from -180 to 180 (although
13!     you may omit either of these endpoints), and latitudes must range
14!     from 90 to -90 (again, you may omit the endpoints).
15! EM 09/2009
16! TN 01/2013 : Adapted for large output files with at least 2 variables > 2 GiB
17
18implicit none
19
20include "netcdf.inc" ! NetCDF definitions
21
22character(len=128) :: infile ! input file name (diagfi.nc or stats.nc format)
23character(len=128) :: outfile ! output file name
24character(len=64) :: text ! to store some text
25character(len=64) :: tmpvarname ! temporarily store a variable name
26integer :: infid ! NetCDF input file ID
27integer :: outfid ! NetCDF output file ID
28integer :: tmpvarid ! temporarily store a variable ID
29integer :: tmpdimid ! temporarily store a dimension ID
30integer :: tmpndims ! temporarily store # of dimensions of a variable
31integer :: nbvarinfile ! # of variables in input file
32integer :: nbattr ! # of attributes of a given variable in input file
33integer :: nbvar3dinfile ! # of 3D (lon,lat,time) variables in input file
34integer :: nbvar4dinfile ! # of 4D (lon,lat,alt,time) variables in input file
35
36integer :: i,j
37integer :: lon_dimid,lat_dimid,alt_dimid,time_dimid ! NetCDF dimension IDs
38integer :: lon_varid,lat_varid,alt_varid,time_varid ! NetCDF variable IDs
39!integer gcm_layers_dimid ! NetCDF dimension ID for # of layers in GCM
40integer :: sigma_varid,aps_varid,bps_varid
41integer :: phisinit_varid
42integer,dimension(4) :: datashape ! shape of 4D datasets
43integer :: ierr ! NetCDF routines return code
44character (len=64), dimension(:), allocatable :: var
45! var(): names of variables that will be processed
46integer nbvar ! # of variables to process
47integer,dimension(:),allocatable :: var_id ! IDs of variables var() (in outfile)
48real :: miss_val=-9.99e+33 ! special "missing value" to specify missing data
49real,parameter :: miss_val_def=-9.99e+33 ! default value for "missing value"
50real,dimension(:),allocatable :: inlat ! input latitude
51integer :: inlatlength  ! # of elements in input latitude
52real,dimension(:),allocatable :: inlon ! input longitude
53integer :: inlonlength ! # of elements in input longitude
54real,dimension(:),allocatable :: alt ! altitude
55integer :: altlength ! # of elements in altitude
56real,dimension(:),allocatable :: time ! input time
57integer :: timelength ! # of elements in time(:)
58real,dimension(:),allocatable :: aps,bps ! hybrid vertical coordinates
59real,dimension(:),allocatable :: sigma ! sigma levels
60real,dimension(:,:),allocatable :: inphisinit ! input ground geopotential
61real,dimension(:),allocatable :: lon ! output longitude
62integer :: lonlength ! # of elements in lon()
63real,dimension(:),allocatable :: wklon ! work longitude (with modulo element)
64integer :: wklonlength ! # of elements in wklon()
65real,dimension(:),allocatable :: lat ! output latitude
66integer :: latlength ! # of elements in lat()
67real,dimension(:),allocatable :: wklat ! work latitude (includes poles)
68integer :: wklatlength ! # of elements in wklat()
69!real,dimension(:),allocatable :: lon_bound ! output longitude boundaries
70!real,dimension(:),allocatable :: lat_bound ! output latitude boundaries
71real,dimension(:),allocatable :: in_lon_bound ! input longitude boundaries
72real,dimension(:),allocatable :: in_lat_bound ! input latitude boundaries
73real,dimension(:),allocatable :: wk_lon_bound ! work longitude boundaries
74real,dimension(:),allocatable :: wk_lat_bound ! work latitude boundaries
75real,dimension(:,:),allocatable :: in_2d_data ! input 2D (lon-lat) dataset
76real,dimension(:,:),allocatable :: wk_2d_data ! work 2D dataset
77real,dimension(:,:),allocatable :: out_2d_data ! output 2D dataset
78real,dimension(:,:,:),allocatable :: in_3d_data ! intput 3D dataset
79real,dimension(:,:,:),allocatable :: wk_3d_data ! work 3D dataset
80real,dimension(:,:,:),allocatable :: out_3d_data ! output 3D dataset
81real,dimension(:,:,:,:),allocatable :: in_4d_data ! intput 4D dataset
82real,dimension(:,:,:,:),allocatable :: wk_4D_data ! work 4D dataset
83real,dimension(:,:,:,:),allocatable :: out_4d_data ! output 4D dataset
84
85real :: pi ! =3.14...
86logical :: have_sigma ! Flag: true if sigma levels are known (false if hybrid
87                      !       coordinates are used)
88logical :: have_geopot ! Flag: true if input file contains ground geopotential
89                       ! phisinit()
90logical :: out_mod_lon ! Flag: true if output grid has modulo longitude (ie:
91                       ! first and last point are in fact at same longitude)
92logical :: out_has_poles ! Flag: true if output grid includes North and South
93                         ! poles
94                         
95integer, dimension(4) :: edges,corner ! needed to write variables for big files
96
97!===============================================================================
98! 1. Input parameters
99!===============================================================================
100pi=2.*asin(1.)
101
102!===============================================================================
103! 1.1 Input file
104!===============================================================================
105
106write(*,*) ""
107write(*,*) " Program valid for diagfi.nc, concatnc.nc and stats.nc files"
108write(*,*) "Enter input file name:"
109
110read(*,'(a128)') infile
111write(*,*) ""
112
113! open input file
114
115ierr = NF_OPEN(infile,NF_NOWRITE,infid)
116if (ierr.ne.NF_NOERR) then
117   write(*,*) 'ERROR: Pb opening file ',trim(infile)
118   stop ""
119endif
120
121!===============================================================================
122! 1.2 Get # and names of variables in input file
123!===============================================================================
124
125ierr=NF_INQ_NVARS(infid,nbvarinfile)
126if (ierr.ne.NF_NOERR) then
127  write(*,*) 'ERROR: Failed geting number of variables from file'
128  stop
129endif
130
131write(*,*)" The following variables have been found:"
132nbvar3dinfile=0
133nbvar4dinfile=0
134do i=1,nbvarinfile
135  ! get name of variable # i
136  ierr=NF_INQ_VARNAME(infid,i,tmpvarname)
137  ! check if it is a 3D variable
138  ierr=NF_INQ_VARNDIMS(infid,i,tmpndims)
139  if (tmpndims.eq.3) then
140    nbvar3dinfile=nbvar3dinfile+1
141    write(*,*) trim(tmpvarname)
142  endif
143  ! check if it is a 4D variable
144  ierr=NF_INQ_VARNDIMS(infid,i,tmpndims)
145  if (tmpndims.eq.4) then
146    nbvar4dinfile=nbvar4dinfile+1
147    write(*,*) trim(tmpvarname)
148  endif
149enddo
150
151allocate(var(nbvar3dinfile+nbvar4dinfile))
152
153write(*,*) ""
154write(*,*) "Which variable do you want to keep?"
155write(*,*) "all or list of <variables> (separated by <Return>s)"
156write(*,*) "(an empty line , i.e: just <Return>, implies end of list)"
157nbvar=0
158read(*,'(a64)') tmpvarname
159do while ((tmpvarname.ne.' ').and.(trim(tmpvarname).ne.'all'))
160  ! check if tmpvarname is valid
161  ierr=NF_INQ_VARID(infid,tmpvarname,tmpvarid)
162  if (ierr.eq.NF_NOERR) then ! valid name
163    nbvar=nbvar+1
164    var(nbvar)=tmpvarname
165  else ! invalid name
166    write(*,*) 'Error: ',trim(tmpvarname),' is not a valid name'
167    write(*,*) '       (we''ll skip that one)'
168  endif
169  read(*,'(a64)') tmpvarname
170enddo
171
172! handle "all" case
173if (tmpvarname.eq.'all') then
174  nbvar=0
175  do i=1,nbvarinfile
176    ! look for 4D variables
177    ierr=NF_INQ_VARNDIMS(infid,i,tmpndims)
178    if (tmpndims.eq.4) then
179      nbvar=nbvar+1
180      ! get the corresponding name
181      ierr=NF_INQ_VARNAME(infid,i,tmpvarname)
182      var(nbvar)=tmpvarname
183    endif
184  enddo
185endif
186
187! Check that there is at least 1 variable to process
188if (nbvar.eq.0) then
189  write(*,*) 'No variables to process !?'
190  write(*,*) 'Might as well stop here'
191  stop ""
192else
193  write(*,*) ""
194  write(*,*) 'OK, the following variables will be processed:'
195  do i=1,nbvar
196    write(*,*) var(i)
197  enddo
198endif
199
200!===============================================================================
201! 1.3 Get input grids in lon,lat,alt,time,
202!     as well as hybrid coordinates aps() and bps() (or sigma levels sigma())
203!     and eventually phisinit() from input file
204!===============================================================================
205
206! 1.3.1 input longitude, latitude, altitude and time
207
208! latitude
209ierr=NF_INQ_DIMID(infid,"latitude",tmpdimid)
210if (ierr.ne.NF_NOERR) then
211  stop "Error: Failed to get latitude dimension ID"
212else
213  ierr=NF_INQ_VARID(infid,"latitude",tmpvarid)
214  if (ierr.ne.NF_NOERR) then
215    stop "Error: Failed to get latitude ID"
216  else
217    ierr=NF_INQ_DIMLEN(infid,tmpdimid,inlatlength)
218    if (ierr.ne.NF_NOERR) then
219      stop "Error: Failed to get latitude length"
220    else
221      allocate(inlat(inlatlength))
222      ierr=NF_GET_VAR_REAL(infid,tmpvarid,inlat)
223      if (ierr.ne.NF_NOERR) then
224        stop "Error: Failed reading latitude"
225      endif
226    endif
227  endif
228endif
229
230! check that these latitudes are 'GCM' latitudes (i.e. poles are included)
231if ((abs(inlat(1)-90.0)).gt.0.001) then
232  write(*,*) "Error: Input latitudes should include north pole, but"
233  write(*,*) " lat(1)=",inlat(1)
234  stop
235endif
236if ((abs(inlat(inlatlength)+90.0)).gt.0.001) then
237  write(*,*) "Error: Input latitudes should include south pole, but"
238  write(*,*) " lat(inlatlength)=",inlat(inlatlength)
239  stop
240endif
241
242
243! longitude
244ierr=NF_INQ_DIMID(infid,"longitude",tmpdimid)
245if (ierr.ne.NF_NOERR) then
246  stop "Error: Failed to get longitude dimension ID"
247else
248  ierr=NF_INQ_VARID(infid,"longitude",tmpvarid)
249  if (ierr.ne.NF_NOERR) then
250    stop "Error: Failed to get longitude ID"
251  else
252    ierr=NF_INQ_DIMLEN(infid,tmpdimid,inlonlength)
253    if (ierr.ne.NF_NOERR) then
254      stop "Error: Failed to get longitude length"
255    else
256      allocate(inlon(inlonlength))
257      ierr=NF_GET_VAR_REAL(infid,tmpvarid,inlon)
258      if (ierr.ne.NF_NOERR) then
259        stop "Error: Failed reading longitude"
260      endif
261    endif
262  endif
263endif
264
265! check that these longitudes are 'GCM' longitudes (i.e range from -180 to 180)
266if (abs(180.+inlon(1)).gt.0.001) then
267  write(*,*) "Error: Input latitudes should start at -180, but"
268  write(*,*) " lon(1)=",inlon(1)
269  stop
270endif
271if (abs(inlon(inlonlength)-180).gt.0.001) then
272  write(*,*) "Error: Input latitudes should end at 180, but"
273  write(*,*) " lon(inlonlength)=",inlon(inlonlength)
274  stop
275endif
276
277
278! altitude
279ierr=NF_INQ_DIMID(infid,"altitude",tmpdimid)
280if (ierr.ne.NF_NOERR) then
281  stop "Error: Failed to get altitude dimension ID"
282else
283  ierr=NF_INQ_VARID(infid,"altitude",tmpvarid)
284  if (ierr.ne.NF_NOERR) then
285    stop "Error: Failed to get altitude ID"
286  else
287    ierr=NF_INQ_DIMLEN(infid,tmpdimid,altlength)
288    if (ierr.ne.NF_NOERR) then
289      stop "Error: Failed to get altitude length"
290    else
291      allocate(alt(altlength))
292      ierr=NF_GET_VAR_REAL(infid,tmpvarid,alt)
293      if (ierr.ne.NF_NOERR) then
294        stop "Error: Failed reading altitude"
295      endif
296    endif
297  endif
298endif
299
300! time
301ierr=NF_INQ_DIMID(infid,"Time",tmpdimid)
302if (ierr.ne.NF_NOERR) then
303  stop "Error: Failed to get Time dimension ID"
304else
305  ierr=NF_INQ_VARID(infid,"Time",tmpvarid)
306  if (ierr.ne.NF_NOERR) then
307    stop "Error: Failed to get Time ID"
308  else
309    ierr=NF_INQ_DIMLEN(infid,tmpdimid,timelength)
310    if (ierr.ne.NF_NOERR) then
311      stop "Error: Failed to get Time length"
312    else
313      allocate(time(timelength))
314      ierr=NF_GET_VAR_REAL(infid,tmpvarid,time)
315      if (ierr.ne.NF_NOERR) then
316        stop "Error: Failed reading Time"
317      endif
318    endif
319  endif
320endif
321
322! 1.3.2 Get hybrid coordinates (or sigma levels)
323
324! start by looking for sigma levels
325ierr=NF_INQ_VARID(infid,"sigma",tmpvarid)
326if (ierr.ne.NF_NOERR) then
327  have_sigma=.false.
328  write(*,*) "Could not find sigma levels... will look for hybrid coordinates"
329else
330  have_sigma=.true.
331  allocate(sigma(altlength))
332  ierr=NF_GET_VAR_REAL(infid,tmpvarid,sigma)
333  if (ierr.ne.NF_NOERR) then
334    stop "Error: Failed reading sigma"
335  endif
336endif
337
338! if no sigma levels, look for hybrid coordinates
339if (.not.have_sigma) then
340  ! hybrid coordinate aps
341  ierr=NF_INQ_VARID(infid,"aps",tmpvarid)
342  if (ierr.ne.NF_NOERR) then
343    stop "Error: Failed to get aps ID"
344  else
345    allocate(aps(altlength))
346    ierr=NF_GET_VAR_REAL(infid,tmpvarid,aps)
347    if (ierr.ne.NF_NOERR) then
348      stop "Error: Failed reading aps"
349    endif
350  endif
351
352  ! hybrid coordinate bps
353  ierr=NF_INQ_VARID(infid,"bps",tmpvarid)
354  if (ierr.ne.NF_NOERR) then
355    stop "Error: Failed to get bps ID"
356  else
357    allocate(bps(altlength))
358    ierr=NF_GET_VAR_REAL(infid,tmpvarid,bps)
359    if (ierr.ne.NF_NOERR) then
360      stop "Error: Failed reading bps"
361    endif
362  endif
363endif !of if (.not.have_sigma)
364
365! 1.3.3 Get ground geopotential phisinit, if available
366
367! look for 'phisinit' in current file
368ierr=NF_INQ_VARID(infid,"phisinit",tmpvarid)
369if (ierr.ne.NF_NOERR) then
370  write(*,*) "Warning: Failed to get phisinit ID from file ",trim(infile)
371  write(*,*) "  ...will not store geopotential in output... "
372  have_geopot=.false.
373else
374  have_geopot=.true.
375  ! Get input physinit
376  allocate(inphisinit(inlonlength,inlatlength))
377  ierr=NF_GET_VAR_REAL(infid,tmpvarid,inphisinit)
378  if (ierr.ne.NF_NOERR) then
379    stop "Error: Failed reading phisinit"
380  endif
381endif
382
383! 1.3.4 Create input longitude and latitude boundaries
384
385! build input longitude boundaries (in radians)
386allocate(in_lon_bound(inlonlength))
387do i=1,inlonlength-1
388  in_lon_bound(i)=0.5*(inlon(i+1)+inlon(i))*pi/180.0
389enddo
390! we have inlon(1)=inlon(inlonlength) modulo 360
391! and thus in_lon_bound(inlonlength)=in_lon_bound(1)+360
392in_lon_bound(inlonlength)=2.*pi+in_lon_bound(1)
393!do i=1,inlonlength
394!  write(*,*) "i=",i,"180/pi*in_lon_bound(i)=",(180./pi)*in_lon_bound(i)
395!enddo
396
397! build input latitude boundaries (in radians)
398allocate(in_lat_bound(inlatlength-1))
399do i=1,inlatlength-1
400  in_lat_bound(i)=0.5*(inlat(i)+inlat(i+1))*pi/180.0
401enddo
402!do i=1,inlatlength-1
403!  write(*,*) "i=",i,"180/pi*in_lat_bound(i)",(180./pi)*in_lat_bound(i)
404!enddo
405
406!===============================================================================
407! 1.4 Get output longitude and latitude coordinates
408!===============================================================================
409
410write(*,*) ""
411write(*,*) "Output horizontal grid:"
412write(*,*) "Number of grid points in longitude?"
413read(*,*) lonlength
414write(*,*) "Enter longitudes (degrees, in [-180:180]), "
415write(*,*) " in increasing order (one per line):"
416allocate(lon(lonlength))
417do i=1,lonlength
418  read(*,*) lon(i)
419enddo
420
421!! build 'work' longitude (which must be a modulo axis; i.e. first and
422!! last points e.g. -180 and 180 are included)
423if (abs((lon(1)+360.-lon(lonlength))).le.0.01) then
424  ! the axis already has modulo endpoints
425  out_mod_lon=.true.
426  wklonlength=lonlength
427  allocate(wklon(wklonlength))
428  wklon(1:lonlength)=lon(1:lonlength)
429else
430  ! add an extra point
431  out_mod_lon=.false.
432  wklonlength=lonlength+1
433  allocate(wklon(wklonlength))
434  wklon(1:lonlength)=lon(1:lonlength)
435  wklon(wklonlength)=wklon(1)+360.0
436endif
437
438! build work longitude boundaries (in radians)
439allocate(wk_lon_bound(wklonlength))
440do i=1,lonlength-1
441  wk_lon_bound(i)=0.5*(wklon(i+1)+wklon(i))*pi/180.0
442enddo
443if (out_mod_lon) then
444  ! we have lon(1)=lon(lonlength) modulo 360
445  ! and thus lon_bound(lonlength)=lon_bound(1)+360
446  wk_lon_bound(wklonlength)=2.*pi+wk_lon_bound(1)
447else
448  wk_lon_bound(wklonlength-1)=0.5*(wklon(wklonlength)+wklon(wklonlength-1))*pi/180.0
449  wk_lon_bound(wklonlength)=2.*pi+wk_lon_bound(1)
450endif
451
452
453write(*,*) "Number of grid points in latitude?"
454read(*,*) latlength
455write(*,*) "Enter latitudes (degrees), in decreasing order (from northernmost"
456write(*,*) " to southernmost), one per line:"
457allocate(lat(latlength))
458do i=1,latlength
459  read(*,*) lat(i)
460enddo
461
462! build 'work' latitude (which must include poles, ie lat=90 and -90)
463if (abs(lat(1)-90.0).le.0.001) then
464  out_has_poles=.true.
465  wklatlength=latlength
466  allocate(wklat(wklatlength))
467  wklat(1:latlength)=lat(1:latlength)
468else
469  out_has_poles=.false.
470  ! add poles
471  wklatlength=latlength+2
472  allocate(wklat(wklatlength))
473  wklat(1)=90
474  wklat(2:latlength+1)=lat(1:latlength)
475  wklat(wklatlength)=-90
476endif
477
478! build work latitude boundaries (in radians)
479allocate(wk_lat_bound(wklatlength-1))
480if (out_has_poles) then
481  do i=1,wklatlength-1
482    wk_lat_bound(i)=0.5*(wklat(i)+wklat(i+1))*pi/180.0
483  enddo
484else
485  ! put northermost boundary near pole
486  wk_lat_bound(1)=(90-0.01*(90.-lat(1)))*pi/180.0
487  do i=2,wklatlength-2
488    wk_lat_bound(i)=0.5*(wklat(i)+wklat(i+1))*pi/180.0
489  enddo
490  ! put southernmost boundary near pole
491  wk_lat_bound(wklatlength-1)=(-90.0-0.01*(-90.-lat(latlength)))*pi/180.0
492endif
493
494!do i=1,wklatlength-1
495!  write(*,*) "i=",i,"180/pi*wk_lat_bound(i)",(180./pi)*wk_lat_bound(i)
496!enddo
497
498!===============================================================================
499! 1.5 Output file
500!===============================================================================
501write(*,*) ""
502outfile=infile(1:len_trim(infile)-3)//"_h.nc"
503write(*,*) "Output file name is: ",trim(outfile)
504
505
506!===============================================================================
507! 2. Create output file and initialize definitions of variables and dimensions
508!===============================================================================
509
510!===============================================================================
511! 2.1. Output file
512!===============================================================================
513
514! Create output file
515ierr=NF_CREATE(outfile,IOR(NF_CLOBBER,NF_64BIT_OFFSET),outfid)
516if (ierr.ne.NF_NOERR) then
517  write(*,*)"Error: could not create file ",outfile
518  stop
519endif
520
521!===============================================================================
522! 2.2. Define dimensions
523!===============================================================================
524! longitude
525ierr=NF_DEF_DIM(outfid,"longitude",lonlength,lon_dimid)
526if (ierr.ne.NF_NOERR) then
527  stop "Error: Could not define longitude dimension"
528endif
529
530! latitude
531ierr=NF_DEF_DIM(outfid,"latitude",latlength,lat_dimid)
532if (ierr.ne.NF_NOERR) then
533  stop "Error: Could not define latitude dimension"
534endif
535
536! altitude
537ierr=NF_DEF_DIM(outfid,"altitude",altlength,alt_dimid)
538if (ierr.ne.NF_NOERR) then
539  stop "Error: Could not define altitude dimension"
540endif
541
542! time
543ierr=NF_DEF_DIM(outfid,"Time",NF_UNLIMITED,time_dimid)
544if (ierr.ne.NF_NOERR) then
545  stop "Error: Could not define latitude dimension"
546endif
547
548!===============================================================================
549! 2.3. Define variables and their attributes
550!===============================================================================
551
552! 2.3.1 Define 1D variables
553
554! longitude
555datashape(1)=lon_dimid
556ierr=NF_DEF_VAR(outfid,"longitude",NF_REAL,1,datashape(1),lon_varid)
557if (ierr.ne.NF_NOERR) then
558  stop "Error: Could not define longitude variable"
559endif
560
561! longitude attributes
562text='east longitude'
563ierr=NF_PUT_ATT_TEXT(outfid,lon_varid,'long_name',len_trim(text),text)
564if (ierr.ne.NF_NOERR) then
565  stop "Error: Problem writing long_name for longitude"
566endif
567text='degrees_east'
568ierr=NF_PUT_ATT_TEXT(outfid,lon_varid,'units',len_trim(text),text)
569if (ierr.ne.NF_NOERR) then
570  stop "Error: Problem writing units for longitude"
571endif
572
573! latitude
574datashape(2)=lat_dimid
575ierr=NF_DEF_VAR(outfid,"latitude",NF_REAL,1,datashape(2),lat_varid)
576if (ierr.ne.NF_NOERR) then
577  stop "Error: Could not define latitude variable"
578endif
579
580! latitude attributes
581text='north latitude'
582ierr=NF_PUT_ATT_TEXT(outfid,lat_varid,'long_name',len_trim(text),text)
583if (ierr.ne.NF_NOERR) then
584  stop "Error: Problem writing long_name for latitude"
585endif
586text='degrees_north'
587ierr=NF_PUT_ATT_TEXT(outfid,lat_varid,'units',len_trim(text),text)
588if (ierr.ne.NF_NOERR) then
589  stop "Error: Problem writing units for latitude"
590endif
591
592! altitude
593datashape(3)=alt_dimid
594ierr=NF_DEF_VAR(outfid,"altitude",NF_REAL,1,datashape(3),alt_varid)
595if (ierr.ne.NF_NOERR) then
596  stop "Error: Could not define altitude variable"
597endif
598
599! altitude attributes
600! preliminary stuff, get input file altitude variable ID
601ierr=NF_INQ_VARID(infid,"altitude",tmpvarid)
602
603! look for a 'long_name' attribute
604text=" "
605ierr=NF_GET_ATT_TEXT(infid,tmpvarid,"long_name",text)
606if (ierr.eq.NF_NOERR) then
607  ! found the attribute; write it to output file
608  ierr=NF_PUT_ATT_TEXT(outfid,alt_varid,'long_name',len_trim(text),text)
609endif
610
611! look for a 'unit' attribute
612text=" "
613ierr=NF_GET_ATT_TEXT(infid,tmpvarid,"units",text)
614if (ierr.eq.NF_NOERR) then
615  ! found the attribute; write it to output file
616  ierr=NF_PUT_ATT_TEXT(outfid,alt_varid,'units',len_trim(text),text)
617endif
618
619! look for a 'positive' attribute
620text=" "
621ierr=NF_GET_ATT_TEXT(infid,tmpvarid,"positive",text)
622if (ierr.eq.NF_NOERR) then
623  ! found the attribute; write it to output file
624  ierr=NF_PUT_ATT_TEXT(outfid,alt_varid,'positive',len_trim(text),text)
625endif
626
627! sigma levels or hybrid coordinates
628if (have_sigma) then
629  ierr=NF_DEF_VAR(outfid,"sigma",NF_REAL,1,alt_dimid,sigma_varid)
630  if (ierr.ne.NF_NOERR) then
631    stop "Error: Could not define sigma variable"
632  endif
633else ! hybrid coordinates
634  ierr=NF_DEF_VAR(outfid,"aps",NF_REAL,1,alt_dimid,aps_varid)
635  if (ierr.ne.NF_NOERR) then
636    stop "Error: Could not define aps variable"
637  endif
638  ierr=NF_DEF_VAR(outfid,"bps",NF_REAL,1,alt_dimid,bps_varid)
639  if (ierr.ne.NF_NOERR) then
640    stop "Error: Could not define bps variable"
641  endif
642endif
643
644! sigma levels (or hybrid coordinates) attributes
645if (have_sigma) then
646  text="sigma levels"
647  ierr=NF_PUT_ATT_TEXT(outfid,sigma_varid,'long_name',len_trim(text),text)
648  if (ierr.ne.NF_NOERR) then
649    stop "Error: Problem writing long_name for sigma"
650  endif
651else ! hybrid coordinates
652  text="hybrid pressure at midlayers"
653  ierr=NF_PUT_ATT_TEXT(outfid,aps_varid,'long_name',len_trim(text),text)
654  if (ierr.ne.NF_NOERR) then
655    stop "Error: Problem writing long_name for aps"
656  endif
657  text="hybrid sigma at midlayers"
658  ierr=NF_PUT_ATT_TEXT(outfid,bps_varid,'long_name',len_trim(text),text)
659  if (ierr.ne.NF_NOERR) then
660    stop "Error: Problem writing long_name for bps"
661  endif
662endif ! of if (have_sigma)
663
664! time
665datashape(4)=time_dimid
666ierr=NF_DEF_VAR(outfid,"Time",NF_REAL,1,datashape(4),time_varid)
667if (ierr.ne.NF_NOERR) then
668  stop "Error: Could not define Time variable"
669endif
670
671! time attributes
672! preliminary stuff, get input file time variable ID
673ierr=NF_INQ_VARID(infid,"Time",tmpvarid)
674! look for a 'unit' attribute
675text=" "
676ierr=NF_GET_ATT_TEXT(infid,tmpvarid,"units",text)
677if (ierr.eq.NF_NOERR) then
678  ! found the attribute; write it to output file
679  ierr=NF_PUT_ATT_TEXT(outfid,time_varid,'units',len_trim(text),text)
680else
681  ! write something not too stupid
682  text='days since 0000-01-1 00:00:00'
683  ierr=NF_PUT_ATT_TEXT(outfid,time_varid,'units',len_trim(text),text)
684  if (ierr.ne.NF_NOERR) then
685    stop "Error: Problem writing units for Time"
686  endif
687endif
688
689! 2.3.2 Define 2D (lon-lat, time independent) variables
690
691! ground geopotential
692if (have_geopot) then
693  if (have_geopot) then
694    ierr=NF_DEF_VAR(outfid,"phisinit",NF_REAL,2,datashape,phisinit_varid)
695    if (ierr.ne.NF_NOERR) then
696      stop "Error: Could not define phisinit variable"
697    endif
698  endif
699  ! ground geopotential attributes
700  text='Geopotential at the surface'
701  ierr=NF_PUT_ATT_TEXT(outfid,phisinit_varid,'long_name',len_trim(text),text)
702  if (ierr.ne.NF_NOERR) then
703    stop "Error: Problem writing long_name for phisinit"
704  endif
705endif ! of if (have_geopot)
706
707! 2.3.3 Define 3D&4D variables (ie: surface/volume+time variables)
708
709! variables requested by user
710allocate(var_id(nbvar))
711do i=1,nbvar
712  ! Get the (input file) ID of the variable (stored in tmpvarid)
713  ierr=NF_INQ_VARID(infid,var(i),tmpvarid)
714  if (ierr.ne.NF_NOERR) then
715    write(*,*) 'Error, failed to get ID for input variable ',trim(var(i))
716    stop ""
717  endif
718 
719  ! Get # of dimensions of the variable
720  ! and define the variable in output file
721  write(*,*) ""
722  write(*,*) "Creating ",trim(var(i))
723  ierr=NF_INQ_VARNDIMS(infid,tmpvarid,tmpndims)
724  if (tmpndims.eq.3) then
725    datashape(3)=time_dimid
726    ierr=NF_DEF_VAR(outfid,var(i),NF_REAL,3,datashape,var_id(i))
727  else
728    if (tmpndims.eq.4) then
729      datashape(3)=alt_dimid
730      datashape(4)=time_dimid
731      ierr=NF_DEF_VAR(outfid,var(i),NF_REAL,4,datashape,var_id(i))
732    else
733      write(*,*) "Error: number of dimensions of input variable ",trim(var(i))
734      write(*,*) "       is ",tmpndims
735      stop
736    endif
737  endif
738  if (ierr.ne.NF_NOERR) then
739    write(*,*) 'Error, could not define variable ',trim(var(i))
740    stop ""
741  endif
742 
743  ! Get the (input file) ID for the variable and
744  ! the # of attributes associated to that variable
745  ierr=NF_INQ_VARNATTS(infid,tmpvarid,nbattr)
746  if (ierr.ne.NF_NOERR) then
747    write(*,*) 'Error, could not get number of attributes for variable ',&
748               trim(var(i))
749    stop ""
750  endif
751  ! inititialize j == number of attributes written to output
752  j=0
753 
754  ! look for a "long_name" attribute (or eventually a 'title' attribute)
755  text='   '
756  ierr=NF_GET_ATT_TEXT(infid,tmpvarid,'long_name',text)
757  if (ierr.ne.NF_NOERR) then ! no long_name attribute
758    ! try to find an equivalent 'title' attribute
759    text='   '
760    ierr=NF_GET_ATT_TEXT(infid,tmpvarid,'title',text)
761    if (ierr.eq.NF_NOERR) then ! found 'title' attribute
762      write(*,*) "Found title ",trim(text)
763      j=j+1
764      ! write it as a 'long_name' attribute
765      ierr=NF_PUT_ATT_TEXT(outfid,var_id(i),'long_name',len_trim(text),text)
766      if (ierr.ne.NF_NOERR) then
767        write(*,*) "Error failed to copy title attribute:",trim(text)
768      stop ""
769      endif
770    endif
771  else ! found long_name; write it to outfile
772    write(*,*) "Found long_name ",trim(text)
773    j=j+1
774    ierr=NF_PUT_ATT_TEXT(outfid,var_id(i),'long_name',len_trim(text),text)
775    if (ierr.ne.NF_NOERR) then
776      write(*,*) "Error failed to copy long_name attribute:",trim(text)
777      stop""
778    endif
779  endif
780 
781  ! look for a "units" attribute
782  text='   '
783  ierr=NF_GET_ATT_TEXT(infid,tmpvarid,'units',text)
784  if (ierr.eq.NF_NOERR) then ! found 'units' attribute
785    write(*,*) "Found units ",trim(text)
786    j=j+1
787    ! write it to output
788    ierr=NF_PUT_ATT_TEXT(outfid,var_id(i),'units',len_trim(text),text)
789    if (ierr.ne.NF_NOERR) then
790      write(*,*) "Error failed to copy units attribute:",trim(text)
791      stop""
792    endif
793  endif
794 
795  ! look for a "missing_value" attribute
796  ierr=NF_GET_ATT_REAL(infid,tmpvarid,"missing_value",miss_val)
797  if (ierr.eq.NF_NOERR) then ! found 'missing_value' attribute
798    write(*,*) "Found missing_value ",miss_val
799    j=j+1
800  else ! no 'missing_value' attribute, set miss_val to default
801    miss_val=miss_val_def
802  endif
803 
804  ! write the missing_value attribute to output
805  ierr=NF_PUT_ATT_REAL(outfid,var_id(i),'missing_value',NF_REAL,1,miss_val)
806  if (ierr.ne.NF_NOERR) then
807    stop "Error, failed to write missing_value attribute"
808  endif
809 
810  ! warn if some attributes were missed
811  if (j.ne.nbattr) then
812    write(*,*)'Warning, it seems some attributes of variable ',trim(var(i))
813    write(*,*)"were not transfered to the new file"
814    write(*,*)'nbattr:',nbattr,' j:',j
815  endif
816
817enddo ! of do i=1,nbvar
818
819!===============================================================================
820! 2.4. Write dimensions (and time-independent variables)
821!===============================================================================
822! Switch out of NetCDF define mode
823ierr=NF_ENDDEF(outfid)
824if (ierr.ne.NF_NOERR) then
825  stop "Error: Could not switch out of define mode"
826endif
827
828! Write longitude
829ierr=NF_PUT_VAR_REAL(outfid,lon_varid,lon)
830if (ierr.ne.NF_NOERR) then
831  stop "Error: Could not write longitude data to output file"
832endif
833
834! Write latitude
835ierr=NF_PUT_VAR_REAL(outfid,lat_varid,lat)
836if (ierr.ne.NF_NOERR) then
837  stop "Error: Could not write latitude data to output file"
838endif
839
840! write altitude
841ierr=NF_PUT_VAR_REAL(outfid,alt_varid,alt)
842if (ierr.ne.NF_NOERR) then
843  stop "Error: Could not write altitude data to output file"
844endif
845
846! Write sigma levels (or hybrid coordinates)
847if (have_sigma) then
848  ierr=NF_PUT_VAR_REAL(outfid,sigma_varid,sigma)
849  if (ierr.ne.NF_NOERR) then
850    stop "Error: Could not write sigma data to output file"
851  endif
852else ! hybrid coordinates
853  ierr=NF_PUT_VAR_REAL(outfid,aps_varid,aps)
854  if (ierr.ne.NF_NOERR) then
855    stop "Error: Could not write aps data to output file"
856  endif
857  ierr=NF_PUT_VAR_REAL(outfid,bps_varid,bps)
858  if (ierr.ne.NF_NOERR) then
859    stop "Error: Could not write bps data to output file"
860  endif
861endif
862
863! write time
864ierr=NF_PUT_VARA_REAL(outfid,time_varid,1,timelength,time)
865if (ierr.ne.NF_NOERR) then
866  stop "Error: Could not write Time data to output file"
867endif
868
869!===============================================================================
870! 3. Interpolate and write 2D variables
871!===============================================================================
872write(*,*) "interpolate 2D variables"
873allocate(in_2d_data(inlonlength,inlatlength))
874allocate(wk_2d_data(wklonlength,wklatlength))
875allocate(out_2d_data(lonlength,latlength))
876
877! ground geopotential
878if (have_geopot) then
879  ! load input geopotential: get ID for input phisinit
880  ierr=NF_INQ_VARID(infid,"phisinit",tmpvarid)
881  if (ierr.ne.NF_NOERR) then
882    stop "Error: Failed to get phisinit ID"
883  endif
884  ! Get physinit
885  ierr=NF_GET_VAR_REAL(infid,tmpvarid,in_2d_data)
886  if (ierr.ne.NF_NOERR) then
887    stop "Error: Failed reading input phisinit"
888  endif
889 
890  ! interpolate onto new grid
891  call interp_horiz(in_2d_data,wk_2d_data,inlonlength-1,inlatlength-1,&
892                    wklonlength-1,wklatlength-1,1,&
893                    in_lon_bound,in_lat_bound,wk_lon_bound,wk_lat_bound)
894 
895  ! copy (and possibly reshape) data to output grid
896  if (out_has_poles) then
897    out_2d_data(1:lonlength,1:latlength)=wk_2d_data(1:lonlength,1:latlength)
898  else
899    out_2d_data(1:lonlength,1:latlength)=wk_2d_data(1:lonlength,2:latlength+1)
900  endif
901 
902  ! write interpolated phisinit to output
903  ierr=NF_PUT_VAR_REAL(outfid,phisinit_varid,out_2d_data)
904  if (ierr.ne.NF_NOERR) then
905    stop "Error: Could not write phisinit data to output file"
906  endif
907endif ! of if (have_geopot)
908
909!===============================================================================
910! 4. Interpolate and write 3D (surface+time) variables
911!===============================================================================
912write(*,*) "interpolate 3D variables"
913allocate(in_3d_data(inlonlength,inlatlength,timelength))
914allocate(wk_3d_data(wklonlength,wklatlength,timelength))
915allocate(out_3d_data(lonlength,latlength,timelength))
916
917do i=1,nbvar ! loop on all selected 3D&4D variables
918  ! get input file ID for this variable
919  ierr=NF_INQ_VARID(infid,var(i),tmpvarid)
920  if (ierr.ne.NF_NOERR) then
921    write(*,*) "Error: Failed to get input ID for ",trim(var(i))
922    stop
923  endif
924  ! get # of dimensions for this variable
925  ierr=NF_INQ_VARNDIMS(infid,tmpvarid,tmpndims)
926
927  if (tmpndims.eq.3) then ! if it is indeed at 3D variable
928    ! get the data
929    ierr=NF_GET_VAR_REAL(infid,tmpvarid,in_3d_data)
930    if (ierr.ne.NF_NOERR) then
931      write(*,*) "Error: Failed reading input ",trim(var(i))
932      stop
933    endif
934    ! interpolate data
935    do j=1,timelength
936      call interp_horiz(in_3d_data(1,1,j),wk_3d_data(1,1,j),&
937                        inlonlength-1,inlatlength-1,&
938                        wklonlength-1,wklatlength-1,1,&
939                        in_lon_bound,in_lat_bound,wk_lon_bound,wk_lat_bound)
940      ! copy (and possibly reshape) data to output grid
941      if (out_has_poles) then
942        out_3d_data(1:lonlength,1:latlength,j)=&
943                        wk_3d_data(1:lonlength,1:latlength,j)
944      else
945        out_3d_data(1:lonlength,1:latlength,j)=&
946                        wk_3d_data(1:lonlength,2:latlength+1,j)
947      endif
948    enddo
949    ! write interpolated data to output
950    corner(:)=1
951    edges(1)=lonlength
952    edges(2)=latlength
953    edges(3)=timelength
954    ierr=NF_PUT_VARA_REAL(outfid,var_id(i),corner(1:3),edges(1:3),out_3d_data)
955    if (ierr.ne.NF_NOERR) then
956      write(*,*) "Error: Could not write ",trim(var(i))," data to output file"
957      stop
958    else
959      write(*,*) "wrote ",trim(var(i))
960    endif
961  endif ! of if (tmpndims.eq.3)
962
963enddo ! of do i=1,nbvar
964
965!===============================================================================
966! 5. Interpolate and write 4D variables
967!===============================================================================
968allocate(in_4d_data(inlonlength,inlatlength,altlength,timelength))
969allocate(wk_4d_data(wklonlength,wklatlength,altlength,timelength))
970allocate(out_4d_data(lonlength,latlength,altlength,timelength))
971
972do i=1,nbvar ! loop on all selected 3D&4D variables
973  ! get input file ID for this variable
974  ierr=NF_INQ_VARID(infid,var(i),tmpvarid)
975  if (ierr.ne.NF_NOERR) then
976    write(*,*) "Error: Failed to get input ID for ",trim(var(i))
977    stop
978  endif
979  ! get # of dimensions for this variable
980  ierr=NF_INQ_VARNDIMS(infid,tmpvarid,tmpndims)
981
982  if (tmpndims.eq.4) then ! if it is indeed at 4D variable
983    ! get the data
984    ierr=NF_GET_VAR_REAL(infid,tmpvarid,in_4d_data)
985    if (ierr.ne.NF_NOERR) then
986      write(*,*) "Error: Failed reading input ",trim(var(i))
987      stop
988    endif
989    ! interpolate data
990    do j=1,timelength
991      call interp_horiz(in_4d_data(1,1,1,j),wk_4d_data(1,1,1,j),&
992                        inlonlength-1,inlatlength-1,&
993                        wklonlength-1,wklatlength-1,altlength,&
994                        in_lon_bound,in_lat_bound,wk_lon_bound,wk_lat_bound)
995      ! copy (and possibly reshape) data to output grid
996      if (out_has_poles) then
997        out_4d_data(1:lonlength,1:latlength,1:altlength,j)=&
998                        wk_4d_data(1:lonlength,1:latlength,1:altlength,j)
999      else
1000        out_4d_data(1:lonlength,1:latlength,1:altlength,j)=&
1001                        wk_4d_data(1:lonlength,2:latlength+1,1:altlength,j)
1002      endif
1003    enddo
1004    ! write interpolated data to output
1005    corner(:)=1
1006    edges(1)=lonlength
1007    edges(2)=latlength
1008    edges(3)=altlength
1009    edges(4)=timelength
1010    ierr=NF_PUT_VARA_REAL(outfid,var_id(i),corner,edges,out_4d_data)
1011    if (ierr.ne.NF_NOERR) then
1012      write(*,*) "Error: Could not write ",trim(var(i))," data to output file"
1013      stop
1014    else
1015      write(*,*) "wrote ",trim(var(i))
1016    endif
1017  endif ! of if (tmpndims.eq.3)
1018
1019enddo ! of do i=1,nbvar
1020
1021!===============================================================================
1022! 6. Close output file
1023!===============================================================================
1024ierr=NF_CLOSE(outfid)
1025if (ierr.ne.NF_NOERR) then
1026  write(*,*) 'Error, failed to close output file ',outfile
1027endif
1028
1029end program hrecast
1030
1031
1032
1033!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1034      subroutine interp_horiz (varo,varn,imo,jmo,imn,jmn,lm,            &
1035     &  rlonuo,rlatvo,rlonun,rlatvn) 
1036
1037!===========================================================
1038!  Interpolation Horizontales des variables d'une grille LMDZ
1039! (des points SCALAIRES au point SCALAIRES)
1040!  dans une autre grille LMDZ en conservant la quantite
1041!  totale pour les variables intensives (/m2) : ex : Pression au sol
1042!
1043! Francois Forget (01/1995)
1044!===========================================================
1045
1046      IMPLICIT NONE
1047
1048!   Declarations:
1049! ==============
1050!
1051!  ARGUMENTS
1052!  """""""""
1053       
1054       INTEGER imo, jmo ! dimensions ancienne grille (input)
1055       INTEGER imn,jmn  ! dimensions nouvelle grille (input)
1056
1057       REAL rlonuo(imo+1)     !  Latitude et
1058       REAL rlatvo(jmo)       !  longitude des
1059       REAL rlonun(imn+1)     !  bord des
1060       REAL rlatvn(jmn)     !  cases "scalaires" (input)
1061
1062       INTEGER lm ! dimension verticale (input)
1063       REAL varo (imo+1, jmo+1,lm) ! var dans l'ancienne grille (input)
1064       REAL varn (imn+1,jmn+1,lm) ! var dans la nouvelle grille (output)
1065
1066! Autres variables
1067! """"""""""""""""
1068       INTEGER imnmx2,jmnmx2
1069!       parameter (imnmx2=190,jmnmx2=100)
1070       parameter (imnmx2=360,jmnmx2=190)
1071       REAL airetest(imnmx2+1,jmnmx2+1)
1072       INTEGER ii,jj,l
1073
1074       REAL airen ((imnmx2+1)*(jmnmx2+1)) ! aire dans la nouvelle grille
1075       REAL airentotn   ! aire totale pole nord dans la nouvelle grille
1076       REAL airentots   ! aire totale pole sud dans la nouvelle grille
1077!    Info sur les ktotal intersection entre les cases new/old grille
1078
1079! kmax: le nombre  max  d'intersections entre les 2 grilles horizontales
1080! On fixe kmax a la taille de la grille des donnees martiennes (360x179)
1081! + des pouiemes (cas ou une maille est a cheval sur 2 ou 4 mailles)
1082!  Il y a un test dans iniinterp_h pour s'assurer que ktotal < kmax
1083       INTEGER kmax, k, ktotal
1084       parameter (kmax = 360*179 + 200000)
1085!      parameter (kmax = 360*179 + 40000)
1086
1087       INTEGER iik(kmax), jjk(kmax),jk(kmax),ik(kmax)
1088       REAL intersec(kmax)
1089       REAL R
1090       REAL totn, tots
1091       integer prev_sumdim
1092       save prev_sumdim
1093       data prev_sumdim /0/
1094       
1095
1096       logical firsttest, aire_ok
1097       save firsttest
1098       data firsttest /.true./
1099       data aire_ok /.true./
1100
1101       integer imoS,jmoS,imnS,jmnS
1102       save imoS,jmoS,imnS,jmnS
1103       save ktotal,iik,jjk,jk,ik,intersec,airen
1104!       REAL pi
1105
1106! Test dimensions imnmx2 jmnmx2
1107!""""""""""""""""""""""""""""""
1108! test dimensionnement tableau airetest
1109      if (imn.GT.imnmx2.OR.jmn.GT.jmnmx2) then
1110         write(*,*) 'STOP pb dimensionnement tableau airetest'
1111         write(*,*) 'il faut imn < imnmx2 et jmn < jmnmx2'
1112         write(*,*) 'imn imnmx2', imn,imnmx2
1113         write(*,*) 'jmn jmnmx2', jmn,jmnmx2
1114         call exit(1)
1115      endif
1116
1117! initialisation
1118! --------------
1119! Si c'est le premier appel,  on prepare l'interpolation
1120! en calculant pour chaque case autour d'un point scalaire de la
1121! nouvelle grille, la surface  de intersection avec chaque
1122!    case de l'ancienne grille.
1123
1124!  This must also be done if we change the dimension
1125      if (imo+jmo+imn+jmn.ne.prev_sumdim) then
1126          firsttest=.true.
1127          prev_sumdim=imo+jmo+imn+jmn
1128      end if       
1129
1130      if (firsttest) then
1131        call iniinterp_h(imo,jmo,imn,jmn ,kmax,                         &
1132     &       rlonuo,rlatvo,rlonun,rlatvn,                               &
1133     &          ktotal,iik,jjk,jk,ik,intersec,airen)
1134       imoS=imo
1135       jmoS=jmo
1136       imnS=imn
1137       jmnS=jmn
1138      else
1139       if(imo.NE.imoS.OR.jmo.NE.jmoS.OR.imn.NE.imnS.OR.jmn.NE.jmnS) then
1140        call iniinterp_h(imo,jmo,imn,jmn ,kmax,                         &
1141     &       rlonuo,rlatvo,rlonun,rlatvn,                               &
1142     &          ktotal,iik,jjk,jk,ik,intersec,airen)
1143       imoS=imo
1144       jmoS=jmo
1145       imnS=imn
1146       jmnS=jmn
1147       end if
1148      end if
1149
1150      do l=1,lm
1151       do jj =1 , jmn+1
1152        do ii=1, imn+1
1153          varn(ii,jj,l) =0.
1154        end do
1155       end do
1156      end do
1157       
1158! Interpolation horizontale
1159! -------------------------
1160! boucle sur toute les ktotal intersections entre les cases
1161! de l'ancienne et la  nouvelle grille
1162!
1163     
1164      do k=1,ktotal
1165        do l=1,lm
1166         varn(iik(k),jjk(k),l) = varn(iik(k),jjk(k),l)                  &
1167     &   + varo(ik(k), jk(k),l)*intersec(k)/airen(iik(k)                &
1168     &   +(jjk(k)-1)*(imn+1))
1169        end do
1170      end do
1171
1172! Une seule valeur au pole pour les variables ! :
1173! -----------------------------------------------
1174      DO l=1, lm
1175         totn =0.
1176         tots =0.
1177
1178
1179! moyenne du champ au poles (ponderee par les aires)
1180!"""""""""""""""""""""""""""""""
1181         airentotn=0.
1182         airentots=0.
1183
1184         do ii =1, imn+1
1185            totn = totn + varn(ii,1,l)*airen(ii)
1186            tots = tots + varn (ii,jmn+1,l)*airen(jmn*(imn+1)+ii)
1187            airentotn=airentotn + airen(ii)
1188            airentots=airentots + airen(jmn*(imn+1)+ii)
1189         end do
1190
1191         do ii =1, imn+1
1192            varn(ii,1,l) = totn/airentotn
1193            varn(ii,jmn+1,l) = tots/airentots
1194         end do
1195
1196      ENDDO
1197           
1198
1199!---------------------------------------------------------------
1200!  TEST  TEST  TEST  TEST  TEST  TEST  TEST  TEST  TEST  TEST
1201      if (firsttest) then
1202!      pi=2.*asin(1.)
1203      firsttest = .false.
1204!      write (*,*) 'INTERP. HORIZ. : TEST SUR LES AIRES:'
1205
1206      do jj =1 , jmn+1
1207        do ii=1, imn+1
1208          airetest(ii,jj) =0.
1209        end do
1210      end do
1211      do k=1,ktotal
1212         airetest(iik(k),jjk(k))= airetest(iik(k),jjk(k)) +intersec(k)
1213      end do
1214      do jj =1 , jmn+1
1215       do ii=1, imn+1
1216         r = airen(ii+(jj-1)*(imn+1))/airetest(ii,jj)
1217         if ((r.gt.1.001).or.(r.lt.0.999)) then
1218             write (*,*) '********** PROBLEME D'' AIRES !!!',           &
1219     &                   ' DANS L''INTERPOLATION HORIZONTALE'
1220             write(*,*)'ii,jj,airen,airetest',                          &
1221     &          ii,jj,airen(ii+(jj-1)*(imn+1)),airetest(ii,jj)
1222             aire_ok = .false.
1223         end if
1224       end do
1225      end do
1226!      if (aire_ok) write(*,*) 'INTERP. HORIZ. : AIRES OK'
1227      endif
1228
1229! FIN TEST  FIN TEST  FIN TEST  FIN TEST  FIN TEST  FIN TEST  FIN TEST
1230! --------------------------------------------------------------------
1231
1232
1233        return
1234        end
1235!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1236      subroutine iniinterp_h (imo,jmo,imn,jmn ,kmax,                    &
1237     &       rlonuo,rlatvo,rlonun,rlatvn,                               &
1238     &       ktotal,iik,jjk,jk,ik,intersec,airen)
1239   
1240      implicit none
1241
1242
1243
1244! ---------------------------------------------------------
1245! Prepare l' interpolation des variables d'une grille LMDZ
1246!  dans une autre grille LMDZ en conservant la quantite
1247!  totale pour les variables intensives (/m2) : ex : Pression au sol
1248!
1249!   (Pour chaque case autour d'un point scalaire de la nouvelle
1250!    grille, on calcule la surface (en m2)en intersection avec chaque
1251!    case de l'ancienne grille , pour la future interpolation)
1252!
1253! on calcule aussi l' aire dans la nouvelle grille
1254!
1255!
1256!   Auteur:  F.Forget 01/1995
1257!   -------
1258!
1259! ---------------------------------------------------------
1260!   Declarations:
1261! ==============
1262!
1263!  ARGUMENTS
1264!  """""""""
1265! INPUT
1266       integer imo, jmo ! dimensions ancienne grille
1267       integer imn,jmn  ! dimensions nouvelle grille
1268       integer kmax ! taille du tableau des intersections
1269       real rlonuo(imo+1)     !  Latitude et
1270       real rlatvo(jmo)       !  longitude des
1271       real rlonun(imn+1)     !  bord des
1272       real rlatvn(jmn)     !  cases "scalaires" (input)
1273
1274! OUTPUT
1275       integer ktotal ! nombre totale d''intersections reperees
1276       integer iik(kmax), jjk(kmax),jk(kmax),ik(kmax)
1277       real intersec(kmax)  ! surface des intersections (m2)
1278       real airen (imn+1,jmn+1) ! aire dans la nouvelle grille
1279
1280
1281       
1282 
1283! Autres variables
1284! """"""""""""""""
1285       integer i,j,ii,jj
1286       integer imomx,jmomx,imnmx1,jmnmx1
1287!       parameter (imomx=361,jmomx=180,imnmx1=190,jmnmx1=100)
1288       parameter (imomx=361,jmomx=180,imnmx1=360,jmnmx1=190)
1289       real a(imomx+1),b(imomx+1),c(jmomx+1),d(jmomx+1)
1290       real an(imnmx1+1),bn(imnmx1+1)
1291       real cn(jmnmx1+1),dn(jmnmx1+1)
1292       real aa, bb,cc,dd
1293       real pi
1294
1295       pi      = 2.*ASIN( 1. )
1296
1297! Test dimensions imnmx1 jmnmx1
1298!""""""""""""""""""""""""""""""
1299! test dimensionnement tableau airetest
1300      if (imn.GT.imnmx1.OR.jmn.GT.jmnmx1) then
1301        write(*,*) 'STOP pb dimensionnement'
1302        write(*,*) 'il faut imn < imnmx1 et jmn < jmnmx1'
1303        write(*,*) 'imn imnmx1', imn,imnmx1
1304        write(*,*) 'jmn jmnmx1', jmn,jmnmx1
1305        call exit(1)
1306      endif
1307
1308      if (imo.GT.imomx.OR.jmo.GT.jmomx) then
1309        write(*,*) 'STOP pb dimensionnement'
1310        write(*,*) 'il faut imo < imomx et jmo < jmomx'
1311        write(*,*) 'imo imomx', imo,imomx
1312        write(*,*) 'jmo jmomx', jmo,jmomx
1313        call exit(1)
1314      endif
1315
1316! On repere les frontieres des cases :
1317! ===================================
1318!
1319! Attention, on ruse avec des latitudes = 90 deg au pole.
1320
1321
1322!  Ancienne grile
1323!  """"""""""""""
1324      a(1) =   - rlonuo(imo+1)
1325      b(1) = rlonuo(1)
1326      do i=2,imo+1
1327         a(i) = rlonuo(i-1)
1328         b(i) =  rlonuo(i)
1329      end do
1330
1331      d(1) = pi/2
1332      do j=1,jmo
1333         c(j) = rlatvo(j)
1334         d(j+1) = rlatvo(j)
1335      end do
1336      c(jmo+1) = -pi/2
1337     
1338
1339!  Nouvelle grille
1340!  """""""""""""""
1341      an(1) =  - rlonun(imn+1)
1342      bn(1) = rlonun(1)
1343      do i=2,imn+1
1344         an(i) = rlonun(i-1)
1345         bn(i) =  rlonun(i)
1346      end do
1347
1348      dn(1) = pi/2
1349      do j=1,jmn
1350         cn(j) = rlatvn(j)
1351         dn(j+1) = rlatvn(j)
1352      end do
1353      cn(jmn+1) = -pi/2
1354
1355! Calcul de la surface des cases scalaires de la nouvelle grille
1356! ==============================================================
1357      do ii=1,imn + 1
1358        do jj = 1,jmn+1
1359           airen(ii,jj) = (bn(ii)-an(ii))*(sin(dn(jj))-sin(cn(jj)))
1360        end do
1361      end do
1362
1363! Calcul de la surface des intersections
1364! ======================================
1365
1366!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1367! test dimenssion kmax (calcul de ktotal)
1368!"""""""""""""""""""""""""""""""""""""""
1369! Calcul de ktotal, mais ralonge beaucoup en temps => pour debug
1370!      write(*,*)
1371!      write(*,*) 'DEBUT DU TEST KMAX'
1372      ktotal = 0
1373      do jj = 1,jmn+1
1374       do j=1, jmo+1
1375          if((cn(jj).lt.d(j)).and.(dn(jj).gt.c(j)))then
1376              do ii=1,imn + 1
1377                do i=1, imo +1
1378                    if (  ((an(ii).lt.b(i)).and.(bn(ii).gt.a(i)))       &
1379     &        .or. ((an(ii).lt.b(i)-2*pi).and.(bn(ii).gt.a(i)-2*pi)     &
1380     &             .and.(b(i)-2*pi.lt.-pi) )                            &
1381     &        .or. ((an(ii).lt.b(i)+2*pi).and.(bn(ii).gt.a(i)+2*pi)     &
1382     &             .and.(a(i)+2*pi.gt.pi) )                             &
1383     &                     )then
1384                      ktotal = ktotal +1
1385                     end if
1386                enddo
1387              enddo
1388           end if
1389        enddo
1390      enddo
1391
1392      if (kmax.LT.ktotal) then
1393         write(*,*)
1394         write(*,*) '******** ATTENTION ********'
1395         write(*,*) 'kmax =',kmax
1396         write(*,*) 'ktotal =',ktotal
1397         write(*,*) 'Changer la valeur de kmax dans interp_horiz.F '
1398         write(*,*) 'avec kmax >= ktotal'
1399         write(*,*) 'EXIT dans iniinterp_h'
1400         call exit(1)
1401      else
1402!         write(*,*) 'kmax =',kmax
1403!         write(*,*) 'ktotal =',ktotal
1404      end if
1405!      write(*,*) 'FIN DU TEST KMAX'
1406!      write(*,*)
1407!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1408
1409!     boucle sur la nouvelle grille
1410!     """"""""""""""""""""""""""""
1411      ktotal = 0
1412      do jj = 1,jmn+1
1413       do j=1, jmo+1
1414          if((cn(jj).lt.d(j)).and.(dn(jj).gt.c(j)))then
1415              do ii=1,imn + 1
1416                do i=1, imo +1
1417                    if (  ((an(ii).lt.b(i)).and.(bn(ii).gt.a(i)))       &
1418     &        .or. ((an(ii).lt.b(i)-2*pi).and.(bn(ii).gt.a(i)-2*pi)     &
1419     &             .and.(b(i)-2*pi.lt.-pi) )                            &
1420     &        .or. ((an(ii).lt.b(i)+2*pi).and.(bn(ii).gt.a(i)+2*pi)     &
1421     &             .and.(a(i)+2*pi.gt.pi) )                             &
1422     &                     )then
1423                      ktotal = ktotal +1
1424                      iik(ktotal) =ii
1425                      jjk(ktotal) =jj
1426                      ik(ktotal) =i
1427                      jk(ktotal) =j
1428
1429                      dd = min(d(j), dn(jj))
1430                      cc = cn(jj)
1431                      if (cc.lt. c(j))cc=c(j)
1432                      if((an(ii).lt.b(i)-2*pi).and.                     &
1433     &                  (bn(ii).gt.a(i)-2*pi)) then
1434                          bb = min(b(i)-2*pi,bn(ii))
1435                          aa = an(ii)
1436                          if (aa.lt.a(i)-2*pi) aa=a(i)-2*pi
1437                      else if((an(ii).lt.b(i)+2*pi).and.                &
1438     &                       (bn(ii).gt.a(i)+2*pi)) then
1439                          bb = min(b(i)+2*pi,bn(ii))
1440                          aa = an(ii)
1441                          if (aa.lt.a(i)+2*pi) aa=a(i)+2*pi
1442                      else
1443                          bb = min(b(i),bn(ii))
1444                          aa = an(ii)
1445                          if (aa.lt.a(i)) aa=a(i)
1446                      end if
1447                      intersec(ktotal)=(bb-aa)*(sin(dd)-sin(cc))
1448                     end if
1449                end do
1450               end do
1451             end if
1452         end do
1453       end do       
1454
1455
1456!     TEST  INFO
1457!     DO k=1,ktotal
1458!      ii = iik(k)
1459!      jj = jjk(k)
1460!      i = ik(k)
1461!      j = jk(k)
1462!      if ((ii.eq.10).and.(jj.eq.10).and.(i.eq.10).and.(j.eq.10))then
1463!      if (jj.eq.2.and.(ii.eq.1))then
1464!          write(*,*) '**************** jj=',jj,'ii=',ii
1465!          write(*,*) 'i,j =',i,j
1466!          write(*,*) 'an,bn,cn,dn', an(ii), bn(ii), cn(jj),dn(jj)
1467!          write(*,*) 'a,b,c,d', a(i), b(i), c(j),d(j)
1468!          write(*,*) 'intersec(k)',intersec(k)
1469!          write(*,*) 'airen(ii,jj)=',airen(ii,jj)
1470!      end if
1471!     END DO
1472
1473      return
1474      end
Note: See TracBrowser for help on using the repository browser.