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

Last change on this file since 3493 was 2567, checked in by emillour, 3 years ago

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