source: trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/obsolete/zrecast.F90 @ 598

Last change on this file since 598 was 11, checked in by aslmd, 14 years ago

spiga@svn-planeto:ajoute le modele meso-echelle martien
File size: 49.1 KB
Line 
1program zrecast
2
3! This program reads 4D (lon-lat-alt-time) fields from GCM output files
4! (ie: diagfi.nc time series or concat.nc or stats.nc files) and, by
5! integrating the hydrostatic equation, recasts data along the vertical
6! direction.
7! The vertical coordinate can be either "above areoid altitudes" or
8! "pressure". Some interpolation along the vertical direction is also
9! done, following instructions given by user.
10! For "above areoid altitudes" output, Atmospheric pressure is added to
11! output dataset; for "pressure coordinate" outputs, the above areoid
12! altitude of pressure is added to output dataset.
13!
14! Minimal requirements and dependencies:
15! The dataset must include the following data:
16! - surface pressure
17! - atmospheric temperature
18! - hybrid coordinates aps() and bps(), or sigma levels() (see section 1.3.2)
19! - ground geopotential (in input file; if not found, it is sought
20!   in a 'diagfi.nc' file. If not found there, it is then sought in
21!   a 'phisinit.nc' file  (see section 1.3.3 of program)
22!
23! - When integration the hydrostatic equation, we assume that R, the molecular
24!   Gas Constant, may not be constant, so it is computed as
25!      R=P/(rho*T) (P=Pressure, rho=density, T=temperature)
26!   If 'rho' is not available, then we use a constant R (see section 2.2)
27!
28! WARNING: Asking for many points along the vertical direction quickly
29!          leads to HUGE output files.
30!
31! EM 09/2006
32! EM 10/2006 : Modified program so that it can now process 'stats.nc'
33!              files obtained from British GCM (ie: vertical coordinate
34!              given as sigma levels and geopotential read from file
35!              'phisinit.nc')
36! EM 01/2006 : Corrected a bug in vertical (log) interpolation for pressure
37!              and density
38
39implicit none
40
41include "netcdf.inc" ! NetCDF definitions
42
43character (len=128) :: infile ! input file name (diagfi.nc or stats.nc format)
44character (len=128) :: infile2 ! second input file (may be needed for 'phisini')
45character (len=128) :: outfile ! output file name
46
47character (len=64) :: text ! to store some text
48character (len=64) :: tmpvarname ! temporarily store a variable name
49integer tmpvarid ! temporarily store a variable ID
50integer tmpdimid ! temporarily store a dimension ID
51integer tmpndims ! temporarily store # of dimensions of a variable
52integer infid ! NetCDF input file ID (of diagfi.nc or stats.nc format)
53integer infid2 ! NetCDF input file which contains 'phisini' dataset (diagfi.nc)
54integer nbvarinfile ! # of variables in input file
55integer nbattr ! # of attributes of a given variable in input file
56integer nbvar4dinfile ! # of 4D (lon,lat,alt,time) variables in input file
57integer outfid ! NetCDF output file ID
58integer lon_dimid,lat_dimid,alt_dimid,time_dimid ! NetCDF dimension IDs
59integer lon_varid,lat_varid,alt_varid,time_varid
60integer za_varid,p_varid ! above areoid and pressure data IDs
61integer,dimension(4) :: datashape
62
63real :: miss_val=-9.99e+33 ! special "missing value" to specify missing data
64real,parameter :: miss_val_def=-9.99e+33 ! default value for "missing value"
65character (len=64), dimension(:), allocatable :: var
66! var(): names of variables that will be processed
67integer nbvar ! # of variables to process
68integer,dimension(:),allocatable :: var_id ! IDs of variables var() (in outfile)
69real,dimension(:),allocatable :: lon ! longitude
70integer lonlength ! # of grid points along longitude
71real,dimension(:),allocatable :: lat ! latitude
72integer latlength ! # of grid points along latitude
73integer altlength ! # of grid point along altitude (of input datasets)
74real,dimension(:),allocatable :: time ! time
75integer timelength ! # of points along time
76real,dimension(:),allocatable :: aps,bps ! hybrid vertical coordinates
77real,dimension(:),allocatable :: sigma ! sigma levels
78real,dimension(:,:),allocatable :: phisinit ! Ground geopotential
79real,dimension(:,:,:),allocatable :: ps ! GCM surface pressure
80real,dimension(:,:,:,:),allocatable :: press ! GCM atmospheric pressure
81real,dimension(:,:,:,:),allocatable :: temp ! GCM atmospheric temperature
82real,dimension(:,:,:,:),allocatable :: rho ! GCM atmospheric density
83real,dimension(:,:,:,:),allocatable :: za_gcm ! GCM above areoid levels
84real,dimension(:,:,:,:),allocatable :: indata ! to store a GCM 4D dataset
85real,dimension(:,:,:,:),allocatable :: outdata ! to store a 4D dataset
86integer ierr ! NetCDF routines return code
87integer i,j,ilon,ilat,ilev,itim ! for loops
88
89integer ztype ! Flag for vertical coordinate of output
90! ztype=1: pressure  ztype=2: above areoid
91integer nblev ! # of levels (along vertical coordinate) for output data
92real pmin,pmax ! min and max values for output pressure coordinate
93real,dimension(:),allocatable :: plevel ! Pressure levels for output data
94real zamin,zamax ! min and max values for output above areoid coordinate
95real,dimension(:),allocatable :: zareoid ! Above areoid heights for output data
96logical :: have_rho ! Flag: true if density 'rho' is available
97logical :: have_sigma ! Flag: true if sigma levels are known (false if hybrid
98                      !       coordinates are used)
99!===============================================================================
100! 1. Input parameters
101!===============================================================================
102
103!===============================================================================
104! 1.1 Input file
105!===============================================================================
106
107write(*,*) ""
108write(*,*) " Program valid for diagfi.nc, concatnc.nc and stats.nc files"
109write(*,*) "Enter input file name:"
110
111read(*,'(a128)') infile
112write(*,*) ""
113
114! open input file
115
116ierr = NF_OPEN(infile,NF_NOWRITE,infid)
117if (ierr.ne.NF_NOERR) then
118   write(*,*) 'ERROR: Pb opening file ',trim(infile)
119   stop ""
120endif
121
122!===============================================================================
123! 1.2 Get # and names of variables in input file
124!===============================================================================
125
126ierr=NF_INQ_NVARS(infid,nbvarinfile)
127if (ierr.ne.NF_NOERR) then
128  write(*,*) 'ERROR: Failed geting number of variables from file'
129  stop
130endif
131
132write(*,*)" The following variables have been found:"
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 4D variable
138  ierr=NF_INQ_VARNDIMS(infid,i,tmpndims)
139  if (tmpndims.eq.4) then
140    nbvar4dinfile=nbvar4dinfile+1
141    write(*,*) trim(tmpvarname)
142  endif
143enddo
144
145allocate(var(nbvar4dinfile))
146
147write(*,*) ""
148write(*,*) "Which variable do you want to keep?"
149write(*,*) "all or list of <variables> (separated by <Return>s)"
150write(*,*) "(an empty line , i.e: just <Return>, implies end of list)"
151nbvar=0
152read(*,'(a64)') tmpvarname
153do while ((tmpvarname.ne.' ').and.(trim(tmpvarname).ne.'all'))
154  ! check if tmpvarname is valid
155  ierr=NF_INQ_VARID(infid,tmpvarname,tmpvarid)
156  if (ierr.eq.NF_NOERR) then ! valid name
157    nbvar=nbvar+1
158    var(nbvar)=tmpvarname
159  else ! invalid name
160    write(*,*) 'Error: ',trim(tmpvarname),' is not a valid name'
161    write(*,*) '       (we''ll skip that one)'
162  endif
163  read(*,'(a64)') tmpvarname
164enddo
165
166! handle "all" case
167if (tmpvarname.eq.'all') then
168  nbvar=0
169  do i=1,nbvarinfile
170    ! look for 4D variables
171    ierr=NF_INQ_VARNDIMS(infid,i,tmpndims)
172    if (tmpndims.eq.4) then
173      nbvar=nbvar+1
174      ! get the corresponding name
175      ierr=NF_INQ_VARNAME(infid,i,tmpvarname)
176      var(nbvar)=tmpvarname
177    endif
178  enddo
179endif
180
181! Check that there is at least 1 variable to process
182if (nbvar.eq.0) then
183  write(*,*) 'No variables to process !?'
184  write(*,*) 'Might as well stop here'
185  stop ""
186else
187  write(*,*) ""
188  write(*,*) 'OK, the following variables will be processed:'
189  do i=1,nbvar
190    write(*,*) var(i)
191  enddo
192endif
193
194!===============================================================================
195! 1.3 Get grids in lon,lat,alt,time,
196!     as well as hybrid coordinates aps() and bps() (or sigma levels sigma())
197!     and phisinit() from input file
198!===============================================================================
199
200! 1.3.1 longitude, latitude, altitude and time
201
202! latitude
203ierr=NF_INQ_DIMID(infid,"latitude",tmpdimid)
204if (ierr.ne.NF_NOERR) then
205  stop "Error: Failed to get latitude dimension ID"
206else
207  ierr=NF_INQ_VARID(infid,"latitude",tmpvarid)
208  if (ierr.ne.NF_NOERR) then
209    stop "Error: Failed to get latitude ID"
210  else
211    ierr=NF_INQ_DIMLEN(infid,tmpdimid,latlength)
212    if (ierr.ne.NF_NOERR) then
213      stop "Error: Failed to get latitude length"
214    else
215      allocate(lat(latlength))
216      ierr=NF_GET_VAR_REAL(infid,tmpvarid,lat)
217      if (ierr.ne.NF_NOERR) then
218        stop "Error: Failed reading latitude"
219      endif
220    endif
221  endif
222endif
223
224! longitude
225ierr=NF_INQ_DIMID(infid,"longitude",tmpdimid)
226if (ierr.ne.NF_NOERR) then
227  stop "Error: Failed to get longitude dimension ID"
228else
229  ierr=NF_INQ_VARID(infid,"longitude",tmpvarid)
230  if (ierr.ne.NF_NOERR) then
231    stop "Error: Failed to get longitude ID"
232  else
233    ierr=NF_INQ_DIMLEN(infid,tmpdimid,lonlength)
234    if (ierr.ne.NF_NOERR) then
235      stop "Error: Failed to get longitude length"
236    else
237      allocate(lon(lonlength))
238      ierr=NF_GET_VAR_REAL(infid,tmpvarid,lon)
239      if (ierr.ne.NF_NOERR) then
240        stop "Error: Failed reading longitude"
241      endif
242    endif
243  endif
244endif
245
246! time
247ierr=NF_INQ_DIMID(infid,"Time",tmpdimid)
248if (ierr.ne.NF_NOERR) then
249  stop "Error: Failed to get Time dimension ID"
250else
251  ierr=NF_INQ_VARID(infid,"Time",tmpvarid)
252  if (ierr.ne.NF_NOERR) then
253    stop "Error: Failed to get Time ID"
254  else
255    ierr=NF_INQ_DIMLEN(infid,tmpdimid,timelength)
256    if (ierr.ne.NF_NOERR) then
257      stop "Error: Failed to get Time length"
258    else
259      allocate(time(timelength))
260      ierr=NF_GET_VAR_REAL(infid,tmpvarid,time)
261      if (ierr.ne.NF_NOERR) then
262        stop "Error: Failed reading Time"
263      endif
264    endif
265  endif
266endif
267
268! altlength
269ierr=NF_INQ_DIMID(infid,"altitude",tmpdimid)
270if (ierr.ne.NF_NOERR) then
271  stop "Error: Failed to get altitude dimension ID"
272else
273  ierr=NF_INQ_DIMLEN(infid,tmpdimid,altlength)
274  if (ierr.ne.NF_NOERR) then
275      stop "Error: Failed to get altitude length"
276  endif
277endif
278
279! 1.3.2 Get hybrid coordinates (or sigma levels)
280
281! start by looking for sigma levels
282ierr=NF_INQ_VARID(infid,"sigma",tmpvarid)
283if (ierr.ne.NF_NOERR) then
284  have_sigma=.false.
285  write(*,*) "Could not find sigma levels... will look for hybrid coordinates"
286else
287  have_sigma=.true.
288  allocate(sigma(altlength))
289  ierr=NF_GET_VAR_REAL(infid,tmpvarid,sigma)
290  if (ierr.ne.NF_NOERR) then
291    stop "Error: Failed reading sigma"
292  endif
293endif
294
295! if no sigma levels, look for hybrid coordinates
296if (.not.have_sigma) then
297  ! hybrid coordinate aps
298  ierr=NF_INQ_VARID(infid,"aps",tmpvarid)
299  if (ierr.ne.NF_NOERR) then
300    stop "Error: Failed to get aps ID"
301  else
302    allocate(aps(altlength))
303    ierr=NF_GET_VAR_REAL(infid,tmpvarid,aps)
304    if (ierr.ne.NF_NOERR) then
305      stop "Error: Failed reading aps"
306    endif
307  endif
308
309  ! hybrid coordinate bps
310  ierr=NF_INQ_VARID(infid,"bps",tmpvarid)
311  if (ierr.ne.NF_NOERR) then
312    stop "Error: Failed to get bps ID"
313  else
314    allocate(bps(altlength))
315    ierr=NF_GET_VAR_REAL(infid,tmpvarid,bps)
316    if (ierr.ne.NF_NOERR) then
317      stop "Error: Failed reading bps"
318    endif
319  endif
320endif !of if (.not.have_sigma)
321
322! 1.3.3 ground geopotential phisinit
323
324allocate(phisinit(lonlength,latlength))
325! look for 'phisinit' in current file
326ierr=NF_INQ_VARID(infid,"phisinit",tmpvarid)
327if (ierr.ne.NF_NOERR) then
328  write(*,*) "Warning: Failed to get phisinit ID from file ",trim(infile)
329  infile2="diagfi.nc"
330  write(*,*) "         Trying file ",trim(infile2)
331  ierr=NF_OPEN(infile2,NF_NOWRITE,infid2)
332  if (ierr.ne.NF_NOERR) then
333    write(*,*) "Problem: Could not find/open that file"
334    infile2="phisinit.nc"
335    write(*,*) "         Trying file ",trim(infile2)
336    ierr=NF_OPEN(infile2,NF_NOWRITE,infid2)
337    if (ierr.ne.NF_NOERR) then
338      write(*,*) "Error: Could not open that file either"
339      stop "Might as well stop here"
340    endif
341  endif
342
343  ! Get ID for phisinit
344  ierr=NF_INQ_VARID(infid2,"phisinit",tmpvarid)
345  if (ierr.ne.NF_NOERR) then
346    stop "Error: Failed to get phisinit ID"
347  endif
348  ! Get physinit
349  ierr=NF_GET_VAR_REAL(infid2,tmpvarid,phisinit)
350  if (ierr.ne.NF_NOERR) then
351    stop "Error: Failed reading phisinit"
352  endif
353  ! Close file
354  write(*,*) 'OK, got phisinit'
355  ierr=NF_CLOSE(infid2)
356else
357  ierr=NF_GET_VAR_REAL(infid,tmpvarid,phisinit)
358  if (ierr.ne.NF_NOERR) then
359    stop "Error: Failed reading phisinit"
360  endif
361endif
362
363!===============================================================================
364! 1.4 Choose and build the new vertical coordinate
365!===============================================================================
366
367write(*,*) ""
368write(*,*) "Which vertical coordinate should the output be in?"
369ztype=0
370do while ((ztype.ne.1).and.(ztype.ne.2))
371  write(*,*) "(1: pressure, 2: above areoid altitude)"
372  read(*,*)ztype
373enddo
374
375! ask for # of points and end values
376write(*,*) ""
377write(*,*) "Enter min and max of vertical coordinate (Pa or m)"
378if (ztype.eq.1) then ! pressure coordinate
379  read(*,*) pmin,pmax
380else  ! above areoid coordinate
381  read(*,*) zamin,zamax
382endif
383write(*,*) "Number of levels along vertical coordinate?"
384read(*,*) nblev
385
386! Build corresponding vertical coordinates
387if (ztype.eq.1) then ! pressure coordinate
388  allocate(plevel(nblev))
389  if (nblev.eq.1) then ! in case only one level is asked for
390    plevel(nblev)=pmin
391  else
392    do i=1,nblev
393      !    build exponentially spread layers
394      plevel(i)=exp(log(pmax)+(log(pmin)-log(pmax))* &
395                    ((real(i)-1.0)/(real(nblev)-1.0)))
396    enddo
397  endif
398else ! above areoid heights
399  allocate(zareoid(nblev))
400  if (nblev.eq.1) then ! in case only one level is asked for
401    zareoid(nblev)=zamin
402  else
403    do i=1,nblev
404      zareoid(i)=zamin+(real(i)-1.0)*((zamax-zamin)/(real(nblev)-1.0))
405    enddo
406  endif
407endif
408
409!===============================================================================
410! 1.5 Get output file name
411!===============================================================================
412write(*,*) ""
413write(*,*) "Enter output file name"
414read(*,*) outfile
415! Francois:
416!if (ztype.eq.1) then ! pressure coordinate
417!  outfile=infile(1:len_trim(infile)-3)//"_P.nc"
418!else  ! above areoid coordinate
419!  outfile=infile(1:len_trim(infile)-3)//"_A.nc"
420!endif
421
422
423!===============================================================================
424! 2.1 Build/store GCM fields which will be used later
425!===============================================================================
426
427!===============================================================================
428! 2.1.1 Surface pressure
429!===============================================================================
430ierr=NF_INQ_VARID(infid,"ps",tmpvarid)
431if (ierr.ne.NF_NOERR) then
432  stop "Error: Failed to get ps ID"
433else
434  allocate(ps(lonlength,latlength,timelength))
435  ierr=NF_GET_VAR_REAL(infid,tmpvarid,ps)
436  if (ierr.ne.NF_NOERR) then
437    stop "Error: Failed reading surface pressure"
438  endif
439endif
440
441!===============================================================================
442! 2.1.2 Atmospheric pressure
443!===============================================================================
444allocate(press(lonlength,latlength,altlength,timelength))
445
446if (have_sigma) then ! sigma coordinate
447  do itim=1,timelength
448    do ilev=1,altlength
449      do ilat=1,latlength
450        do ilon=1,lonlength
451          press(ilon,ilat,ilev,itim)=sigma(ilev)*ps(ilon,ilat,itim)
452        enddo
453      enddo
454    enddo
455  enddo
456else ! hybrid coordinates
457  do itim=1,timelength
458    do ilev=1,altlength
459      do ilat=1,latlength
460        do ilon=1,lonlength
461          press(ilon,ilat,ilev,itim)=aps(ilev)+bps(ilev)*ps(ilon,ilat,itim)
462        enddo
463      enddo
464    enddo
465  enddo
466endif
467
468!===============================================================================
469! 2.1.3 Atmospheric temperature
470!===============================================================================
471allocate(temp(lonlength,latlength,altlength,timelength))
472
473ierr=NF_INQ_VARID(infid,"temp",tmpvarid)
474if (ierr.ne.NF_NOERR) then
475  ! stop "Error: Failed to get temp ID"
476  ! try "t" for temperature
477  ierr=NF_INQ_VARID(infid,"t",tmpvarid)
478  if (ierr.ne.NF_NOERR) then
479    stop "Error: Failed to get t ID"
480  else
481    ierr=NF_GET_VAR_REAL(infid,tmpvarid,temp)
482    if (ierr.ne.NF_NOERR) then
483      stop "Error: Failed reading atmospheric temperature"
484    endif
485  endif
486else
487  ierr=NF_GET_VAR_REAL(infid,tmpvarid,temp)
488  if (ierr.ne.NF_NOERR) then
489    stop "Error: Failed reading atmospheric temperature"
490  endif
491endif
492
493!===============================================================================
494! 2.1.4 Atmospheric density
495!===============================================================================
496
497ierr=NF_INQ_VARID(infid,"rho",tmpvarid)
498if (ierr.ne.NF_NOERR) then
499  write(*,*) "Warning: Failed to get rho ID"
500  have_rho=.false.
501else
502  have_rho=.true.
503  allocate(rho(lonlength,latlength,altlength,timelength))
504  ierr=NF_GET_VAR_REAL(infid,tmpvarid,rho)
505  if (ierr.ne.NF_NOERR) then
506    stop "Error: Failed reading atmospheric density"
507  endif
508endif
509
510!===============================================================================
511! 2.2 Build GCM Above areoid altitudes of GCM nodes
512!===============================================================================
513allocate(za_gcm(lonlength,latlength,altlength,timelength))
514
515if (have_rho) then
516  call build_gcm_alt(lonlength,latlength,altlength,timelength, &
517                     phisinit,ps,press,temp,rho,za_gcm)
518else
519  write(*,*)"Warning: Using constant R to integrate hydrostatic equation"
520  call crude_gcm_alt(lonlength,latlength,altlength,timelength, &
521                     phisinit,ps,press,temp,za_gcm)
522endif
523
524!===============================================================================
525! 3. Create output file and initialize definitions of variables and dimensions
526!===============================================================================
527
528!===============================================================================
529! 3.1. Output file
530!===============================================================================
531
532! Create output file
533ierr=NF_CREATE(outfile,NF_CLOBBER,outfid)
534if (ierr.ne.NF_NOERR) then
535  write(*,*)"Error: could not create file ",outfile
536  stop
537endif
538
539!===============================================================================
540! 3.2. Define dimensions
541!===============================================================================
542! longitude
543ierr=NF_DEF_DIM(outfid,"longitude",lonlength,lon_dimid)
544if (ierr.ne.NF_NOERR) then
545  stop "Error: Could not define longitude dimension"
546endif
547
548! latitude
549ierr=NF_DEF_DIM(outfid,"latitude",latlength,lat_dimid)
550if (ierr.ne.NF_NOERR) then
551  stop "Error: Could not define latitude dimension"
552endif
553
554! altitude
555ierr=NF_DEF_DIM(outfid,"altitude",nblev,alt_dimid)
556if (ierr.ne.NF_NOERR) then
557  stop "Error: Could not define altitude dimension"
558endif
559
560! time
561ierr=NF_DEF_DIM(outfid,"Time",timelength,time_dimid)
562if (ierr.ne.NF_NOERR) then
563  stop "Error: Could not define latitude dimension"
564endif
565
566!===============================================================================
567! 3.3. Define variables and their attributes
568!===============================================================================
569! longitude
570datashape(1)=lon_dimid
571ierr=NF_DEF_VAR(outfid,"longitude",NF_REAL,1,datashape(1),lon_varid)
572if (ierr.ne.NF_NOERR) then
573  stop "Error: Could not define longitude variable"
574endif
575
576! longitude attributes
577text='east longitude'
578ierr=NF_PUT_ATT_TEXT(outfid,lon_varid,'long_name',len_trim(text),text)
579if (ierr.ne.NF_NOERR) then
580  stop "Error: Problem writing long_name for longitude"
581endif
582text='degrees_east'
583ierr=NF_PUT_ATT_TEXT(outfid,lon_varid,'units',len_trim(text),text)
584if (ierr.ne.NF_NOERR) then
585  stop "Error: Problem writing units for longitude"
586endif
587
588! latitude
589datashape(2)=lat_dimid
590ierr=NF_DEF_VAR(outfid,"latitude",NF_REAL,1,datashape(2),lat_varid)
591if (ierr.ne.NF_NOERR) then
592  stop "Error: Could not define latitude variable"
593endif
594
595! latitude attributes
596text='north latitude'
597ierr=NF_PUT_ATT_TEXT(outfid,lat_varid,'long_name',len_trim(text),text)
598if (ierr.ne.NF_NOERR) then
599  stop "Error: Problem writing long_name for latitude"
600endif
601text='degrees_north'
602ierr=NF_PUT_ATT_TEXT(outfid,lat_varid,'units',len_trim(text),text)
603if (ierr.ne.NF_NOERR) then
604  stop "Error: Problem writing units for latitude"
605endif
606
607! altitude
608datashape(3)=alt_dimid
609ierr=NF_DEF_VAR(outfid,"altitude",NF_REAL,1,datashape(3),alt_varid)
610if (ierr.ne.NF_NOERR) then
611  stop "Error: Could not define altitude variable"
612endif
613
614!altitude attributes
615if (ztype.eq.1) then ! pressure vertical coordinate
616  text='Pressure levels'
617  ierr=NF_PUT_ATT_TEXT(outfid,alt_varid,'long_name',len_trim(text),text)
618  if (ierr.ne.NF_NOERR) then
619    stop "Error: Problem writing long_name for altitude"
620  endif
621  text='Pa'
622  ierr=NF_PUT_ATT_TEXT(outfid,alt_varid,'units',len_trim(text),text)
623  if (ierr.ne.NF_NOERR) then
624    stop "Error: Problem writing units for altitude"
625  endif
626  text='down'
627  ierr=NF_PUT_ATT_TEXT(outfid,alt_varid,'positive',len_trim(text),text)
628  if (ierr.ne.NF_NOERR) then
629    stop "Error: Problem writing positive for altitude"
630  endif
631else ! above areoid vertical coordinate
632  text='Altitude above areoid'
633  ierr=NF_PUT_ATT_TEXT(outfid,alt_varid,'long_name',len_trim(text),text)
634  if (ierr.ne.NF_NOERR) then
635    stop "Error: Problem writing long_name for altitude"
636  endif
637  text='m'
638  ierr=NF_PUT_ATT_TEXT(outfid,alt_varid,'units',len_trim(text),text)
639  if (ierr.ne.NF_NOERR) then
640    stop "Error: Problem writing units for altitude"
641  endif
642  text='up'
643  ierr=NF_PUT_ATT_TEXT(outfid,alt_varid,'positive',len_trim(text),text)
644  if (ierr.ne.NF_NOERR) then
645    stop "Error: Problem writing positive for altitude"
646  endif
647endif
648
649! time
650datashape(4)=time_dimid
651ierr=NF_DEF_VAR(outfid,"Time",NF_REAL,1,datashape(4),time_varid)
652if (ierr.ne.NF_NOERR) then
653  stop "Error: Could not define Time variable"
654endif
655
656! time attributes
657text='Time'
658ierr=NF_PUT_ATT_TEXT(outfid,time_varid,'long_name',len_trim(text),text)
659if (ierr.ne.NF_NOERR) then
660  stop "Error: Problem writing long_name for Time"
661endif
662text='days since 0000-01-1 00:00:00'
663ierr=NF_PUT_ATT_TEXT(outfid,time_varid,'units',len_trim(text),text)
664if (ierr.ne.NF_NOERR) then
665  stop "Error: Problem writing units for Time"
666endif
667
668! add pressure or zareoid
669if (ztype.eq.1) then ! pressure vertical coordinate
670  ! zareoid dataset
671  ierr=NF_DEF_VAR(outfid,"zareoid",NF_REAL,4,datashape,za_varid)
672  if (ierr.ne.NF_NOERR) then
673    stop "Error: Could not define zareoid variable"
674  endif
675  ! zareoid attributes
676  text='altitude above areoid'
677  ierr=NF_PUT_ATT_TEXT(outfid,za_varid,'long_name',len_trim(text),text)
678  if (ierr.ne.NF_NOERR) then
679    stop "Error: Problem writing long_name for zareoid"
680  endif
681  text='m'
682  ierr=NF_PUT_ATT_TEXT(outfid,za_varid,'units',len_trim(text),text)
683  if (ierr.ne.NF_NOERR) then
684    stop "Error: Problem writing units for zareoid"
685  endif
686  ! zareoid missing value
687  ierr=NF_PUT_ATT_REAL(outfid,za_varid,'missing_value',NF_REAL,1,miss_val)
688  if (ierr.ne.NF_NOERR) then
689    stop "Error: Problem writing missing_value for zareoid"
690  endif
691else ! above areoid vertical coordinate
692  ! pressure dataset
693  ierr=NF_DEF_VAR(outfid,"pressure",NF_REAL,4,datashape,p_varid)
694  if (ierr.ne.NF_NOERR) then
695    stop "Error: Could not define pressure variable"
696  endif
697  ! pressure attributes
698  text='Atmospheric pressure'
699  ierr=NF_PUT_ATT_TEXT(outfid,p_varid,'long_name',len_trim(text),text)
700  if (ierr.ne.NF_NOERR) then
701    stop "Error: Problem writing long_name for pressure"
702  endif
703  text='Pa'
704  ierr=NF_PUT_ATT_TEXT(outfid,p_varid,'units',len_trim(text),text)
705  if (ierr.ne.NF_NOERR) then
706    stop "Error: Problem writing units for pressure"
707  endif
708  ! pressure missing value
709  ierr=NF_PUT_ATT_REAL(outfid,p_varid,'missing_value',NF_REAL,1,miss_val)
710  if (ierr.ne.NF_NOERR) then
711    stop "Error: Problem writing missing_value for pressure"
712  endif
713endif
714
715! Handle 4D variables
716allocate(var_id(nbvar))
717do i=1,nbvar
718  write(*,*) ""
719  write(*,*) "Processing ",trim(var(i))
720  ! define the variable
721  ierr=NF_DEF_VAR(outfid,var(i),NF_REAL,4,datashape,var_id(i))
722  if (ierr.ne.NF_NOERR) then
723    write(*,*) 'Error, could not define variable ',trim(var(i))
724    stop ""
725  endif
726
727  ! Get the (input file) ID for the variable and
728  ! the # of attributes associated to that variable
729  ierr=NF_INQ_VARID(infid,var(i),tmpvarid)
730  if (ierr.ne.NF_NOERR) then
731    write(*,*) 'Error, failed to get ID for variable ',trim(var(i))
732    stop ""
733  endif
734  ierr=NF_INQ_VARNATTS(infid,tmpvarid,nbattr)
735  if (ierr.ne.NF_NOERR) then
736    write(*,*) 'Error, could not get number of attributes for variable ',&
737               trim(var(i))
738    stop ""
739  endif
740  ! inititialize j == number of attributes written to output
741  j=0
742
743  ! look for a "long_name" attribute
744  text='   '
745  ierr=NF_GET_ATT_TEXT(infid,tmpvarid,'long_name',text)
746  if (ierr.ne.NF_NOERR) then ! no long_name attribute
747    ! try to find an equivalent 'title' attribute
748    text='   '
749    ierr=NF_GET_ATT_TEXT(infid,tmpvarid,'title',text)
750    if (ierr.eq.NF_NOERR) then ! found 'title' attribute
751      write(*,*) "Found title ",trim(text)
752      j=j+1
753      ! write it as a 'long_name' attribute
754      ierr=NF_PUT_ATT_TEXT(outfid,var_id(i),'long_name',len_trim(text),text)
755      if (ierr.ne.NF_NOERR) then
756        write(*,*) "Error failed to copy title attribute:",trim(text)
757      stop ""
758      endif
759    endif
760  else ! found long_name; write it to outfile
761    write(*,*) "Found long_name ",trim(text)
762    j=j+1
763    ierr=NF_PUT_ATT_TEXT(outfid,var_id(i),'long_name',len_trim(text),text)
764    if (ierr.ne.NF_NOERR) then
765      write(*,*) "Error failed to copy long_name attribute:",trim(text)
766      stop""
767    endif
768  endif
769 
770  ! look for a "units" attribute
771  text='   '
772  ierr=NF_GET_ATT_TEXT(infid,tmpvarid,'units',text)
773  if (ierr.eq.NF_NOERR) then ! found 'units' attribute
774    write(*,*) "Found units ",trim(text)
775    j=j+1
776    ! write it to output
777    ierr=NF_PUT_ATT_TEXT(outfid,var_id(i),'units',len_trim(text),text)
778    if (ierr.ne.NF_NOERR) then
779      write(*,*) "Error failed to copy units attribute:",trim(text)
780      stop""
781    endif
782  endif
783 
784  ! look for a "missing_value" attribute
785  ierr=NF_GET_ATT_REAL(infid,tmpvarid,"missing_value",miss_val)
786  if (ierr.eq.NF_NOERR) then ! found 'missing_value' attribute
787    write(*,*) "Found missing_value ",miss_val
788    j=j+1
789  else ! no 'missing_value' attribute, set miss_val to default
790    miss_val=miss_val_def
791  endif
792 
793  ! write the missing_value attribute to output
794  ierr=NF_PUT_ATT_REAL(outfid,var_id(i),'missing_value',NF_REAL,1,miss_val)
795  if (ierr.ne.NF_NOERR) then
796    stop "Error, failed to write missing_value attribute"
797  endif
798
799  ! warn if some attributes were missed
800  if (j.ne.nbattr) then
801    write(*,*)'Warning, it seems some attributes of variable ',trim(var(i))
802    write(*,*)"were not transfered to the new file"
803    write(*,*)'nbattr:',nbattr,' j:',j
804  endif
805   
806enddo ! of do=1,nbvar
807
808
809!===============================================================================
810! 3.4. Write dimensions
811!===============================================================================
812! Switch out of NetCDF define mode
813ierr=NF_ENDDEF(outfid)
814if (ierr.ne.NF_NOERR) then
815  stop "Error: Could not switch out of define mode"
816endif
817
818! Write longitude
819ierr=NF_PUT_VAR_REAL(outfid,lon_varid,lon)
820if (ierr.ne.NF_NOERR) then
821  stop "Error: Could not write longitude data to output file"
822endif
823
824! Write latitude
825ierr=NF_PUT_VAR_REAL(outfid,lat_varid,lat)
826if (ierr.ne.NF_NOERR) then
827  stop "Error: Could not write latitude data to output file"
828endif
829
830! Write altitude
831if (ztype.eq.1) then ! pressure vertical coordinate
832  ierr=NF_PUT_VAR_REAL(outfid,alt_varid,plevel)
833  if (ierr.ne.NF_NOERR) then
834    stop "Error: Could not write altitude data to output file"
835  endif
836else ! above areoid altitude
837  ierr=NF_PUT_VAR_REAL(outfid,alt_varid,zareoid)
838  if (ierr.ne.NF_NOERR) then
839    stop "Error: Could not write altitude data to output file"
840  endif
841endif
842
843! Write time
844ierr=NF_PUT_VAR_REAL(outfid,time_varid,time)
845if (ierr.ne.NF_NOERR) then
846  stop "Error: Could not write Time data to output file"
847endif
848
849!===============================================================================
850! 4. Interpolate and write variables
851!===============================================================================
852
853! 4.0 Allocations
854!indata() to store input (on GCM grid) data
855allocate(indata(lonlength,latlength,altlength,timelength))
856! outdata() to store output (on new vertical grid) data
857allocate(outdata(lonlength,latlength,nblev,timelength))
858
859! 4.1 If output is in pressure coordinate
860if (ztype.eq.1) then
861  do i=1,nbvar ! loop on 4D variable to process
862  ! identify and read a dataset
863  ierr=NF_INQ_VARID(infid,var(i),tmpvarid)
864  if (ierr.ne.NF_NOERR) then
865    write(*,*) 'Error, failed to get ID for variable ',var(i)
866    stop
867  endif
868  ierr=NF_GET_VAR_REAL(infid,tmpvarid,indata)
869  if (ierr.ne.NF_NOERR) then
870    write(*,*) 'Error, failed to load variable ',var(i)
871    stop
872  endif
873 
874  ! interpolate dataset onto new grid
875  call p_coord_interp(lonlength,latlength,altlength,timelength,nblev, &
876                      miss_val,ps,press,indata,plevel,outdata)
877 
878  ! write new dataset to output file
879  ierr=NF_PUT_VAR_REAL(outfid,var_id(i),outdata)
880  if (ierr.ne.NF_NOERR) then
881    write(*,*)'Error, Failed to write variable ',var(i)
882    stop
883  endif
884 
885  enddo ! of do i=1,nbvar
886 
887  ! interpolate zareoid onto new grid
888  call p_coord_interp(lonlength,latlength,altlength,timelength,nblev, &
889                      miss_val,ps,press,za_gcm,plevel,outdata)
890  ! write result to output file
891  ierr=NF_PUT_VAR_REAL(outfid,za_varid,outdata)
892  if (ierr.ne.NF_NOERR) then
893    stop "Error, Failed to write zareoid to output file"
894  endif
895endif ! of if (ztype.eq.1)
896
897! 4.2 If output is in above areoid altitude
898if (ztype.eq.2) then
899  do i=1,nbvar ! loop on 4D variable to process
900  ! identify and read a dataset
901  ierr=NF_INQ_VARID(infid,var(i),tmpvarid)
902  if (ierr.ne.NF_NOERR) then
903    write(*,*) 'Error, failed to get ID for variable ',var(i)
904    stop
905  endif
906  ierr=NF_GET_VAR_REAL(infid,tmpvarid,indata)
907  if (ierr.ne.NF_NOERR) then
908    write(*,*) 'Error, failed to load variable ',var(i)
909    stop
910  endif
911 
912  ! interpolate dataset onto new grid
913  ! check if variable is "rho" (to set flag for interpolation below)
914  if (var(i).eq.'rho') then
915    j=1
916  else
917    j=0
918  endif
919 
920  call za_coord_interp(lonlength,latlength,altlength,timelength,nblev, &
921                      miss_val,za_gcm,indata,j,zareoid,outdata)
922
923  ! write new dataset to output file
924  ierr=NF_PUT_VAR_REAL(outfid,var_id(i),outdata)
925  if (ierr.ne.NF_NOERR) then
926    write(*,*)'Error, Failed to write variable ',var(i)
927    stop
928  endif
929
930  enddo ! of do i=1,nbvar
931 
932  ! interpolate pressure onto new grid
933  call za_coord_interp(lonlength,latlength,altlength,timelength,nblev, &
934                      miss_val,za_gcm,press,1,zareoid,outdata)
935
936  ! write result to output file
937  ierr=NF_PUT_VAR_REAL(outfid,p_varid,outdata)
938  if (ierr.ne.NF_NOERR) then
939    stop "Error, Failed to write pressure to output file"
940  endif
941endif ! of if (ztype.eq.2)
942
943! 4.3 Close output file
944ierr=NF_CLOSE(outfid)
945if (ierr.ne.NF_NOERR) then
946  write(*,*) 'Error, failed to close output file ',outfile
947endif
948
949end program
950
951!===============================================================================
952
953subroutine build_gcm_alt(lonlength,latlength,altlength,timelength, &
954                         phis,ps,press,temp,rho,za_gcm)
955!==============================================================================
956! Purpose: Integrate hydrostatic equation in order to build the "above areoid
957!          altitudes" corresponding to GCM atmospheric levels
958!==============================================================================
959implicit none
960!==============================================================================
961! Arguments:
962!==============================================================================
963integer,intent(in) :: lonlength ! # of points along longitude
964integer,intent(in) :: latlength ! # of points along latitude
965integer,intent(in) :: altlength ! # of points along altitude
966integer,intent(in) :: timelength ! # of stored time steps
967real,dimension(lonlength,latlength),intent(in) :: phis
968! phis(:,:) is the ground geopotential
969real,dimension(lonlength,latlength,timelength),intent(in) :: ps
970! ps(:,:) is the surface pressure
971real,dimension(lonlength,latlength,altlength,timelength),intent(in) :: press
972! press(:,:,:,:) is atmospheric pressure
973real,dimension(lonlength,latlength,altlength,timelength),intent(in) :: temp
974! temp(:,:,:,:) is atmospheric temperature
975real,dimension(lonlength,latlength,altlength,timelength),intent(in) :: rho
976! rho(:,:,:,:) is atmospheric density
977
978real,dimension(lonlength,latlength,altlength,timelength),intent(out) :: za_gcm
979! za_gcm(:,:,:,:) are the above aroid heights of GCM levels
980
981!===============================================================================
982! Local variables:
983!===============================================================================
984real,dimension(:),allocatable :: sigma ! GCM sigma levels
985real,dimension(:,:,:,:),allocatable :: R ! molecular gas constant
986real,dimension(:,:,:,:),allocatable :: zlocal ! local above surface altitude
987real :: Rmean ! mean value of R for a given level
988real :: Tmean ! "mean" value of temperature for a given level
989integer iloop,jloop,kloop,tloop
990
991! Parameters needed to integrate hydrostatic equation:
992real :: g0 ! exact mean gravity at radius=3396.km (Lemoine et al. 2001)
993data    g0 /3.7257964/
994real :: a0 
995data    a0 /3396.E3/  ! radius=3396.km
996real :: gz
997! gz: gravitational acceleration at a given (zareoid) altitude
998
999!===============================================================================
1000! 1. Various initialisations
1001!===============================================================================
1002allocate(sigma(altlength))
1003allocate(R(lonlength,latlength,altlength,timelength))
1004allocate(zlocal(lonlength,latlength,altlength,timelength))
1005
1006!==============================================================================
1007! 2. Compute Molecular Gas Constant (J kg-1 K-1) at every grid point
1008!   --later used to integrate the hydrostatic equation--
1009!==============================================================================
1010do tloop=1,timelength
1011  do kloop=1,altlength
1012    do jloop=1,latlength
1013      do iloop=1,lonlength
1014        R(iloop,jloop,kloop,tloop)=press(iloop,jloop,kloop,tloop)/ &
1015                                            (rho(iloop,jloop,kloop,tloop)* &
1016                                             temp(iloop,jloop,kloop,tloop))
1017      enddo
1018    enddo
1019  enddo
1020enddo
1021
1022!===============================================================================
1023! 3. Integrate hydrostatic equation to compute zlocal and za_gcm
1024!===============================================================================
1025do tloop=1,timelength
1026  do jloop=1,latlength
1027    do iloop=1,lonlength
1028      ! handle case of first altitude level
1029      sigma(1)=press(iloop,jloop,1,tloop)/ps(iloop,jloop,tloop)
1030      zlocal(iloop,jloop,1,tloop)=-log(sigma(1))*R(iloop,jloop,1,tloop)* &
1031                                                 temp(iloop,jloop,1,tloop)/g0
1032      za_gcm(iloop,jloop,1,tloop)=zlocal(iloop,jloop,1,tloop)+ &
1033                                  phis(iloop,jloop)/g0
1034      do kloop=2,altlength
1035        ! compute sigma level of layer
1036        sigma(kloop)=press(iloop,jloop,kloop,tloop)/ps(iloop,jloop,tloop)
1037       
1038        ! compute "mean" temperature of the layer
1039        if(temp(iloop,jloop,kloop,tloop).eq.  &
1040           temp(iloop,jloop,kloop-1,tloop)) then
1041          Tmean=temp(iloop,jloop,kloop,tloop)
1042        else
1043          Tmean=(temp(iloop,jloop,kloop,tloop)-      &
1044                 temp(iloop,jloop,kloop-1,tloop))/   &
1045                 log(temp(iloop,jloop,kloop,tloop)/  &
1046                     temp(iloop,jloop,kloop-1,tloop))
1047        endif
1048       
1049        ! compute mean value of R of the layer
1050        Rmean=0.5*(R(iloop,jloop,kloop,tloop)+R(iloop,jloop,kloop-1,tloop))
1051       
1052        ! compute gravitational acceleration (at altitude zaeroid(kloop-1))
1053        gz=g0*(a0**2)/(a0+za_gcm(iloop,jloop,kloop-1,tloop))**2
1054       
1055        ! compute zlocal(kloop)
1056        zlocal(iloop,jloop,kloop,tloop)=zlocal(iloop,jloop,kloop-1,tloop)- &
1057                log(sigma(kloop)/sigma(kloop-1))*Rmean*Tmean/gz
1058       
1059        ! compute za_gcm(kloop)
1060        za_gcm(iloop,jloop,kloop,tloop)=zlocal(iloop,jloop,kloop,tloop)+ &
1061                                        phis(iloop,jloop)/g0
1062      enddo ! kloop
1063    enddo ! iloop
1064  enddo ! jloop
1065enddo ! tloop
1066
1067! Cleanup
1068deallocate(sigma)
1069deallocate(R)
1070deallocate(zlocal)
1071
1072end subroutine build_gcm_alt
1073
1074!===============================================================================
1075
1076subroutine crude_gcm_alt(lonlength,latlength,altlength,timelength, &
1077                         phis,ps,press,temp,za_gcm)
1078!==============================================================================
1079! Purpose: Integrate hydrostatic equation in order to build the "above areoid
1080!          altitudes" corresponding to GCM atmospheric levels
1081!==============================================================================
1082implicit none
1083!==============================================================================
1084! Arguments:
1085!==============================================================================
1086integer,intent(in) :: lonlength ! # of points along longitude
1087integer,intent(in) :: latlength ! # of points along latitude
1088integer,intent(in) :: altlength ! # of points along altitude
1089integer,intent(in) :: timelength ! # of stored time steps
1090real,dimension(lonlength,latlength),intent(in) :: phis
1091! phis(:,:) is the ground geopotential
1092real,dimension(lonlength,latlength,timelength),intent(in) :: ps
1093! ps(:,:) is the surface pressure
1094real,dimension(lonlength,latlength,altlength,timelength),intent(in) :: press
1095! press(:,:,:,:) is atmospheric pressure
1096real,dimension(lonlength,latlength,altlength,timelength),intent(in) :: temp
1097! temp(:,:,:,:) is atmospheric temperature
1098
1099real,dimension(lonlength,latlength,altlength,timelength),intent(out) :: za_gcm
1100! za_gcm(:,:,:,:) are the above aroid heights of GCM levels
1101
1102!===============================================================================
1103! Local variables:
1104!===============================================================================
1105real,dimension(:),allocatable :: sigma ! GCM sigma levels
1106real,parameter :: R=190 ! molecular gas constant
1107real,dimension(:,:,:,:),allocatable :: zlocal ! local above surface altitude
1108real :: Tmean ! "mean" value of temperature for a given level
1109integer iloop,jloop,kloop,tloop
1110
1111! Parameters needed to integrate hydrostatic equation:
1112real :: g0 ! exact mean gravity at radius=3396.km (Lemoine et al. 2001)
1113data    g0 /3.7257964/
1114real :: a0 
1115data    a0 /3396.E3/  ! radius=3396.km
1116real :: gz
1117! gz: gravitational acceleration at a given (zareoid) altitude
1118
1119!===============================================================================
1120! 1. Various initialisations
1121!===============================================================================
1122allocate(sigma(altlength))
1123allocate(zlocal(lonlength,latlength,altlength,timelength))
1124
1125!===============================================================================
1126! 2. Integrate hydrostatic equation to compute zlocal and za_gcm
1127!===============================================================================
1128do tloop=1,timelength
1129  do jloop=1,latlength
1130    do iloop=1,lonlength
1131      ! handle case of first altitude level
1132      sigma(1)=press(iloop,jloop,1,tloop)/ps(iloop,jloop,tloop)
1133      zlocal(iloop,jloop,1,tloop)=-log(sigma(1))* &
1134                                               R*temp(iloop,jloop,1,tloop)/g0
1135      za_gcm(iloop,jloop,1,tloop)=zlocal(iloop,jloop,1,tloop)+ &
1136                                  phis(iloop,jloop)/g0
1137      do kloop=2,altlength
1138        ! compute sigma level of layer
1139        sigma(kloop)=press(iloop,jloop,kloop,tloop)/ps(iloop,jloop,tloop)
1140       
1141        ! compute "mean" temperature of the layer
1142        if(temp(iloop,jloop,kloop,tloop).eq.  &
1143           temp(iloop,jloop,kloop-1,tloop)) then
1144          Tmean=temp(iloop,jloop,kloop,tloop)
1145        else
1146          Tmean=(temp(iloop,jloop,kloop,tloop)-      &
1147                 temp(iloop,jloop,kloop-1,tloop))/   &
1148                 log(temp(iloop,jloop,kloop,tloop)/  &
1149                     temp(iloop,jloop,kloop-1,tloop))
1150        endif
1151               
1152        ! compute gravitational acceleration (at altitude zaeroid(kloop-1))
1153        gz=g0*(a0**2)/(a0+za_gcm(iloop,jloop,kloop-1,tloop))**2
1154       
1155        ! compute zlocal(kloop)
1156        zlocal(iloop,jloop,kloop,tloop)=zlocal(iloop,jloop,kloop-1,tloop)- &
1157                log(sigma(kloop)/sigma(kloop-1))*R*Tmean/gz
1158       
1159        ! compute za_gcm(kloop)
1160        za_gcm(iloop,jloop,kloop,tloop)=zlocal(iloop,jloop,kloop,tloop)+ &
1161                                        phis(iloop,jloop)/g0
1162      enddo ! kloop
1163    enddo ! iloop
1164  enddo ! jloop
1165enddo ! tloop
1166
1167! Cleanup
1168deallocate(sigma)
1169deallocate(zlocal)
1170
1171end subroutine crude_gcm_alt
1172
1173!===============================================================================
1174
1175subroutine p_coord_interp(lonlen,latlen,altlen,tlen,newaltlen, &
1176                          missing,ps,press,gcmdata,plevels,newdata)
1177!==============================================================================
1178! Purpose:
1179! Recast a 4D (spatio-temporal) GCM dataset, which has a vertical coordinate
1180! in pseudo-altitude, into a dataset which has a vertical coordinate at given
1181! pressure levels
1182!==============================================================================
1183implicit none
1184!==============================================================================
1185! Arguments:
1186!==============================================================================
1187integer,intent(in) :: lonlen ! # of points along longitude (GCM dataset)
1188integer,intent(in) :: latlen ! # of points along latitude (GCM dataset)
1189integer,intent(in) :: altlen ! # of points along altitude (GCM dataset)
1190integer,intent(in) :: tlen   ! # of stored time steps (GCM dataset)
1191integer,intent(in) :: newaltlen ! # of points along altitude
1192real,intent(in) :: missing ! missing value
1193real,dimension(lonlen,latlen,tlen),intent(in) :: ps ! GCM surface pressure
1194real,dimension(lonlen,latlen,altlen,tlen),intent(in) :: press ! GCM pressure
1195real,dimension(lonlen,latlen,altlen,tlen),intent(in) :: gcmdata ! GCM dataset
1196real,dimension(newaltlen),intent(in) :: plevels
1197! plevels(:) pressure levels at which gcmdata is to be interpolated
1198
1199real,dimension(lonlen,latlen,newaltlen,tlen),intent(out) :: newdata
1200! newdata(:,:,:,:) gcmdata recasted along vertical coordinate
1201!==============================================================================
1202! Local variables:
1203!==============================================================================
1204real,dimension(:),allocatable :: lnp
1205! lnp(:): used to store log(pressure) values
1206real,dimension(:),allocatable :: q
1207! q(:): used to store values along vertical direction
1208real :: x,y
1209! x,y: temporary variables
1210integer :: iloop,jloop,kloop,tloop
1211
1212allocate(lnp(altlen))
1213allocate(q(altlen))
1214
1215do tloop=1,tlen
1216  do jloop=1,latlen
1217    do iloop=1,lonlen
1218      ! build lnp() and q() along vertical coordinate
1219      do kloop=1,altlen
1220        lnp(kloop)=-log(press(iloop,jloop,kloop,tloop))
1221        q(kloop)=gcmdata(iloop,jloop,kloop,tloop)
1222      enddo
1223     
1224      ! Interpolation 
1225      do kloop=1,newaltlen
1226        ! compute, by interpolation, value at pressure level plevels(kloop)
1227        if ((plevels(kloop).le.ps(iloop,jloop,tloop)).and. &
1228            (plevels(kloop).ge.press(iloop,jloop,altlen,tloop))) then
1229          x=-log(plevels(kloop))
1230          call interpolf(x,y,missing,lnp,q,altlen)
1231          newdata(iloop,jloop,kloop,tloop) = y
1232        else ! if plevels(kloop) is out of range,
1233          ! assign a "missing_value" at this node
1234          newdata(iloop,jloop,kloop,tloop) = missing
1235        endif
1236      enddo
1237     
1238    enddo !iloop
1239  enddo !jloop
1240enddo !tloop
1241
1242! Cleanup
1243deallocate(lnp)
1244deallocate(q)
1245
1246end subroutine p_coord_interp
1247
1248!===============================================================================
1249
1250subroutine za_coord_interp(lonlen,latlen,altlen,tlen,newaltlen, &
1251                          missing,za_gcm,gcmdata,flag,zareoid,newdata)
1252!==============================================================================
1253! Purpose:
1254! Recast a 4D (spatio-temporal) GCM dataset, for which corresponding grid
1255! values of above areoid altitude are known, along chosen above areoid
1256! altitude
1257!==============================================================================
1258implicit none
1259!==============================================================================
1260! Arguments:
1261!==============================================================================
1262integer,intent(in) :: lonlen ! # of points along longitude (GCM dataset)
1263integer,intent(in) :: latlen ! # of points along latitude (GCM dataset)
1264integer,intent(in) :: altlen ! # of points along altitude (GCM dataset)
1265integer,intent(in) :: tlen   ! # of stored time steps (GCM dataset)
1266integer,intent(in) :: newaltlen ! # of points along altitude
1267real,intent(in) :: missing ! missing value
1268real,dimension(lonlen,latlen,altlen,tlen),intent(in) :: za_gcm
1269!za_gcm(:,:,:,) above areoid altitude of GCM grid points
1270real,dimension(lonlen,latlen,altlen,tlen),intent(in) :: gcmdata ! GCM dataset
1271integer,intent(in) :: flag ! flag (==1 for 'log' interpolation)
1272! flag==0 (standard linear interpolation)
1273real,dimension(newaltlen),intent(in) :: zareoid
1274! zareoid(:) above areoid altitudes at which data must be recast
1275
1276real,dimension(lonlen,latlen,newaltlen,tlen),intent(out) :: newdata
1277! newdata(:,:,:,:) gcmdata recasted along vertical coordinate
1278
1279!==============================================================================
1280! Local variables:
1281!==============================================================================
1282real,dimension(:),allocatable :: za,q
1283real,dimension(:),allocatable :: logq
1284! za(:): used to store za_gcm at a given location
1285! q(:): used to store field values along the vertical direction
1286real :: x,y ! temporary variables
1287integer :: iloop,jloop,kloop,tloop
1288
1289allocate(za(altlen))
1290allocate(q(altlen))
1291allocate(logq(altlen))
1292
1293if (flag.eq.0) then
1294 do tloop=1,tlen
1295  do jloop=1,latlen
1296    do iloop=1,lonlen
1297      do kloop=1,altlen
1298        ! extract the vertical coordinates zaij(), or build p_za()
1299        za(kloop)=za_gcm(iloop,jloop,kloop,tloop)
1300        ! store values along altitude
1301        q(kloop)=gcmdata(iloop,jloop,kloop,tloop)
1302      enddo !kloop
1303      ! Interpolation
1304      do kloop=1,newaltlen
1305      ! Check if zareoid(kloop) is within range
1306      if ((zareoid(kloop).ge.za_gcm(iloop,jloop,1,tloop)).and. &
1307          (zareoid(kloop).le.za_gcm(iloop,jloop,altlen,tloop))) then
1308        ! zareoid(kloop) is within range
1309        x=zareoid(kloop)
1310        call interpolf(x,y,missing,za,q,altlen)
1311        newdata(iloop,jloop,kloop,tloop)=y
1312      else ! zareoid(kloop) is out of range
1313        newdata(iloop,jloop,kloop,tloop)=missing
1314      endif
1315      enddo !kloop
1316    enddo !iloop
1317  enddo !jloop
1318 enddo !tloop
1319else ! case when flag=1 (i.e.: rho)
1320 do tloop=1,tlen
1321  do jloop=1,latlen
1322    do iloop=1,lonlen
1323      do kloop=1,altlen
1324        ! extract the vertical coordinates zaij(), or build p_za()
1325        za(kloop)=za_gcm(iloop,jloop,kloop,tloop)
1326        ! extract the vertical coordinates zaij(), or build p_za()
1327        za(kloop)=za_gcm(iloop,jloop,kloop,tloop)
1328        ! store log values along altitude
1329        logq(kloop)=log(gcmdata(iloop,jloop,kloop,tloop))
1330      enddo !kloop
1331
1332      ! Interpolation
1333      do kloop=1,newaltlen
1334      ! Check if zareoid(kloop) is within range
1335      if ((zareoid(kloop).ge.za_gcm(iloop,jloop,1,tloop)).and. &
1336          (zareoid(kloop).le.za_gcm(iloop,jloop,altlen,tloop))) then
1337        ! zareoid(kloop) is within range
1338        x=zareoid(kloop)
1339        call interpolf(x,y,missing,za,logq,altlen)
1340        newdata(iloop,jloop,kloop,tloop)=exp(y)
1341      else  ! zareoid(kloop) is out of range
1342        newdata(iloop,jloop,kloop,tloop)=missing
1343      endif
1344      enddo !kloop
1345    enddo !iloop
1346  enddo !jloop
1347 enddo !tloop
1348endif
1349
1350! Cleanup
1351deallocate(za)
1352deallocate(q)
1353deallocate(logq)
1354
1355end subroutine za_coord_interp
1356
1357!===============================================================================
1358
1359subroutine interpolf(x,y,missing,xd,yd,nd)
1360!==============================================================================
1361! Purpose:
1362! Yield y=f(x), where xd() end yd() are arrays of known values,
1363! using linear interpolation
1364! If x is not included in the interval spaned by xd(), then y is set
1365! to a default value 'missing'
1366! Note:
1367! Array xd() should contain ordered (either increasing or decreasing) abscissas
1368!==============================================================================
1369implicit none
1370!==============================================================================
1371! Arguments:
1372!==============================================================================
1373real,intent(in) :: x ! abscissa at which interpolation is to be done
1374real,intent(in) :: missing ! missing value (if no interpolation is performed)
1375integer :: nd ! size of arrays
1376real,dimension(nd),intent(in) :: xd ! array of known absissas
1377real,dimension(nd),intent(in) :: yd ! array of correponding values
1378
1379real,intent(out) :: y ! interpolated value
1380!==============================================================================
1381! Local variables:
1382!==============================================================================
1383integer :: i
1384
1385! default: set y to 'missing'
1386y=missing
1387
1388   do i=1,nd-1
1389     if (((x.ge.xd(i)).and.(x.le.xd(i+1))).or.&
1390          ((x.le.xd(i)).and.(x.ge.xd(i+1)))) then
1391        y=yd(i)+(x-xd(i))*(yd(i+1)-yd(i))/(xd(i+1)-xd(i))
1392        exit
1393     endif
1394   enddo
1395
1396
1397end subroutine interpolf
Note: See TracBrowser for help on using the repository browser.