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

Last change on this file since 253 was 137, checked in by aslmd, 14 years ago

ATTENTION COMMIT MAJEUR NON-CONSERVATIF

  • CHANGEMENT ARBORESCENCE POUR SEPARATION CLAIRE DES COMPOSANTES et MODELES
  • UTILISATION DE LIENS SYMBOLIQUES POUR GARDER UNE BASE COMMUNE LMDZ.COMMON
  • VOIR LE FICHIER DOC/000-MODELS POUR EN SAVOIR PLUS !

EN CAS DE PROBLEME AVEC svn update IL EST CONSEILLE DE REVENIR A UN svn co

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