source: trunk/LMDZ.MARS/util/gencol.F90

Last change on this file was 3795, checked in by jmauxion, 12 days ago

Mars PCM:
Add a utility to compute the column density from concentration, mass and volume ratios of various tracers.
JM
File size: 38.7 KB
Line 
1program genCol
2
3! ----------------------------------------------------------------------------
4! This program can re-compute the column density (kg m-2)
5! of any tracer from its mmr, vmr or concentration
6! input : diagfi.nc  / concat.nc / stats.nc kind of files
7! output : input_col.nc including lon,lat,alt,time and all
8! specified tracer_col
9! author: F. Forget, M. Alhossani, J.Mauxion
10! ----------------------------------------------------------------------------
11
12implicit none
13
14include "netcdf.inc" ! NetCDF definitions
15
16!!! I/O variables
17character (len=128) :: infile ! input file name (diagfi.nc or stats.nc format)
18character (len=128) :: infile2 ! second input file (may be needed for 'aire' and 'cv')
19integer infid ! NetCDF input file ID (of diagfi.nc or stats.nc format)
20integer infid2 ! NetCDF input file which contains 'phisini' dataset (diagfi.nc)
21character (len=100) :: filename
22integer :: ierr ! ierr: [netcdf] subroutine returned error code
23integer :: nout,latdimout,londimout,altdimout,timedimout
24! nout: [netcdf] output file ID
25! latdimout: [netcdf] output latitude (dimension) ID
26! londimout: [netcdf] output longitude (dimension) ID
27! altdimout: [netcdf] output altitude (dimension) ID
28! timedimout: [netcdf] output time (dimension) ID
29
30!!! Temporary variables
31character (len=64) :: tmpvarname ! temporary store a variable name
32character (len=64) :: tmpunit ! temporary store a unit
33integer :: tmpvarid ! varid: temporary store a variable ID #
34integer :: tmpdimid ! temporarily store a dimension ID
35integer tmpndims ! temporarily store # of dimensions of a variable
36real :: tmpmm ! Temporarily store the molar mass of current specie
37
38!!! Auto-detection and tracer selection
39integer nbvarinfile ! # of variables in input file
40integer nbtrinfile ! # of 4D mmr, vmr or concentration tracers in input file
41integer nbvar ! # of variables to process
42character (len=64), dimension(:), allocatable :: var
43! var(): name(s) of variable(s) that will be concatenated
44
45!!! Grid-related variables
46real, dimension(:), allocatable:: lat,lon,alt,time
47! lat(): array, stores latitude coordinates
48! lon(): array, stores longitude coordinates
49! alt(): array, stores altitude coordinates
50! time(): array, stores time coordinates
51integer lonlength ! # of grid points along longitude
52integer latlength ! # of grid points along latitude
53integer altlength ! # of grid point along altitude (of input datasets)
54integer timelength ! # of points along time
55real,dimension(:),allocatable :: ap,bp ! hybrid vertical coordinates of levels
56real,dimension(:),allocatable :: aps,bps ! hybrid vertical coordinates of layerS
57
58!!! Loops-related variables
59integer :: ilat,ilon,ialt,it,i
60! ilat: index for loops on latitude
61! ilon: index for loops on longitude
62! ialt: index for loops on altitude
63! it: index for loops on time
64! i: index for loops
65
66!!! Constants
67real :: g ! Martian gravitational acceleration
68real :: nav ! Avogadro number
69real :: rgp ! Perfect gas constant
70real :: Rmean ! Mean atmospheric molar mass
71
72!!! Booleans to handle various setup configuration
73logical :: isco, ismmr, isvmr ! is the tracer co, mmr, vmr ?
74logical :: ismmcst ! true if mm assumed constant, false otherwise
75logical :: foundmm ! to store whether we find molar mass in input file
76logical :: foundpress ! to store whether we find pressure in input file
77logical :: foundrho ! to store whether we find density in input file
78
79!!! Incrementing
80real :: incMass ! mass of air in each layer
81real :: tracerMass ! cumulated column density
82real :: cstmm ! Constant astmospheric molar mass
83
84!!! Variable to read from input file
85real,dimension(:,:,:),allocatable :: ps ! GCM surface pressure
86real,dimension(:,:,:,:),allocatable :: tracer ! vmr,mmr or concentration
87real,dimension(:,:,:,:),allocatable :: temp ! atmospheric temperature
88real,dimension(:,:,:,:),allocatable :: rho ! atmospheric density
89real,dimension(:,:,:,:),allocatable :: press ! atmospheric pressure
90real,dimension(:,:,:,:),allocatable :: locmm ! local atmospheric molar mass
91
92!!! New variable to write in output file
93integer :: tracvarout ! var id for column density
94real,dimension(:,:,:),allocatable :: tracerCol ! total tracer column
95
96!===============================================================================
97! 1.1 Input file
98!===============================================================================
99
100write(*,*) ""
101write(*,*) "Enter input file name:"
102
103read(*,'(a128)') infile
104write(*,*) ""
105
106! open input file
107
108ierr = NF_OPEN(infile,NF_NOWRITE,infid)
109if (ierr.ne.NF_NOERR) then
110   write(*,*) 'ERROR: Pb opening input file'
111   stop ""
112endif
113
114!===============================================================================
115! 1.2 Get # and names of tracers in input file
116!===============================================================================
117
118! Get number of var in input file
119ierr=NF_INQ_NVARS(infid,nbvarinfile)
120if (ierr.ne.NF_NOERR) then
121   write(*,*) 'ERROR: Failed getting number of variables from file'
122   write(*,*) NF_STRERROR(ierr)
123   stop ""
124endif
125
126! Find which var is a tracer
127write(*,*)" The following tracers have been found:"
128nbtrinfile=0
129do i=1,nbvarinfile
130   ! Get name of variable # i
131   ierr=NF_INQ_VARNAME(infid,i,tmpvarname)
132   ! Get dimension and unit
133   ierr=NF_INQ_VARNDIMS(infid,i,tmpndims)
134   ierr=NF_GET_ATT_TEXT(infid,i,'units',tmpunit)
135   ismmr = (tmpunit.eq.'kg/kg').or.(tmpunit.eq.'kg.kg-1')
136   isvmr = (tmpunit.eq.'mol/mol').or.(tmpunit.eq.'mol.mol-1')
137   isco = (tmpunit.eq.'cm-3').or.(tmpunit.eq.'m-3')
138   if ((tmpndims.eq.4).and.(ismmr.or.isvmr.or.isco)) then
139      nbtrinfile=nbtrinfile+1
140      write(*,*) trim(tmpvarname), " ", trim(tmpunit)
141   endif
142enddo
143
144! Allocate array for kept tracers
145allocate(var(nbtrinfile),stat=ierr)
146if (ierr.ne.0) then
147  write(*,*) "Error: failed memory allocation of var(nbtrinfile)"
148  stop ""
149endif
150
151! Ask for tracers to compute
152write(*,*) ""
153write(*,*) "Which variable do you want to keep?"
154write(*,*) "all or list of <variables> (separated by <Return>s)"
155write(*,*) "(an empty line , i.e: just <Return>, implies end of list)"
156nbvar=0
157read(*,'(a64)') tmpvarname
158do while ((tmpvarname.ne.' ').and.(trim(tmpvarname).ne.'all'))
159   ! check if tmpvarname is valid
160   ierr=NF_INQ_VARID(infid,tmpvarname,tmpvarid)
161   if (ierr.eq.NF_NOERR) then ! valid name
162      ! also check that it is indeed 4D tracer with consistent unit
163      ierr=NF_INQ_VARNDIMS(infid,tmpvarid,tmpndims)
164      ierr=NF_GET_ATT_TEXT(infid,tmpvarid,'units',tmpunit)
165      ismmr = (tmpunit.eq.'kg/kg').or.(tmpunit.eq.'kg.kg-1')
166      isvmr = (tmpunit.eq.'mol/mol').or.(tmpunit.eq.'mol.mol-1')
167      isco = (tmpunit.eq.'cm-3').or.(tmpunit.eq.'m-3')
168      if ((tmpndims.eq.4).and.(ismmr.or.isvmr.or.isco)) then
169         nbvar=nbvar+1
170         var(nbvar)=tmpvarname
171      else ! invalid dimension or unit
172         write(*,*) 'Error: ',trim(tmpvarname),' has wrong dimension and/or unit'
173         write(*,*) '       (we''ll skip that one)'
174      endif
175   else ! invalid name
176      write(*,*) 'Error: ',trim(tmpvarname),' is not a valid name'
177      write(*,*) '       (we''ll skip that one)'
178   endif
179   read(*,'(a64)') tmpvarname
180enddo
181
182! handle "all" case
183if (tmpvarname.eq.'all') then
184   nbvar=0
185   do i=1,nbvarinfile
186      ! look for 4D tracers with good unit
187      ierr=NF_INQ_VARNDIMS(infid,i,tmpndims)
188      ierr=NF_GET_ATT_TEXT(infid,i,'units',tmpunit)
189      ismmr = (tmpunit.eq.'kg/kg').or.(tmpunit.eq.'kg.kg-1')
190      isvmr = (tmpunit.eq.'mol/mol').or.(tmpunit.eq.'mol.mol-1')
191      isco = (tmpunit.eq.'cm-3').or.(tmpunit.eq.'m-3')
192      if ((tmpndims.eq.4).and.(ismmr.or.isvmr.or.isco)) then
193         nbvar=nbvar+1
194         ! get the corresponding name
195         ierr=NF_INQ_VARNAME(infid,i,tmpvarname)
196         var(nbvar)=tmpvarname
197      endif
198   enddo
199endif
200
201! Check that there is at least 1 variable to process
202if (nbvar.eq.0) then
203   write(*,*) 'No variables to process !?'
204   write(*,*) 'Might as well stop here'
205   stop ""
206else
207   write(*,*) ""
208   write(*,*) 'OK, the following variables will be processed:'
209   do i=1,nbvar
210      write(*,*) var(i)
211   enddo
212endif
213
214
215!===============================================================================
216! 1.3 Get grids in lon,lat,alt,time,
217!     as well as hybrid coordinates aps() and bps()
218!===============================================================================
219
220! 1.3.1 longitude, latitude, altitude and time
221
222! latitude
223ierr=NF_INQ_DIMID(infid,"latitude",tmpdimid)
224if (ierr.ne.NF_NOERR) then
225   stop "Error: Failed to get latitude dimension ID"
226else
227   ierr=NF_INQ_VARID(infid,"latitude",tmpvarid)
228   if (ierr.ne.NF_NOERR) then
229      stop "Error: Failed to get latitude ID"
230   else
231      ierr=NF_INQ_DIMLEN(infid,tmpdimid,latlength)
232      if (ierr.ne.NF_NOERR) then
233         stop "Error: Failed to get latitude length"
234      else
235         allocate(lat(latlength))
236         ierr=NF_GET_VAR_REAL(infid,tmpvarid,lat)
237         if (ierr.ne.NF_NOERR) then
238            write(*,*) "Error: Failed reading latitude"
239            stop ""
240         endif
241      endif
242   endif
243endif
244
245! longitude
246ierr=NF_INQ_DIMID(infid,"longitude",tmpdimid)
247if (ierr.ne.NF_NOERR) then
248   stop "Error: Failed to get longitude dimension ID"
249else
250   ierr=NF_INQ_VARID(infid,"longitude",tmpvarid)
251   if (ierr.ne.NF_NOERR) then
252      stop "Error: Failed to get longitude ID"
253   else
254      ierr=NF_INQ_DIMLEN(infid,tmpdimid,lonlength)
255      if (ierr.ne.NF_NOERR) then
256         stop "Error: Failed to get longitude length"
257      else
258         allocate(lon(lonlength))
259         ierr=NF_GET_VAR_REAL(infid,tmpvarid,lon)
260         if (ierr.ne.NF_NOERR) then
261            write(*,*) "Error: Failed reading longitude"
262            stop ""     
263         endif
264      endif
265   endif
266endif
267
268! time
269ierr=NF_INQ_DIMID(infid,"Time",tmpdimid)
270if (ierr.ne.NF_NOERR) then
271   stop "Error: Failed to get Time dimension ID"
272else
273   ierr=NF_INQ_VARID(infid,"Time",tmpvarid)
274   if (ierr.ne.NF_NOERR) then
275      stop "Error: Failed to get Time ID"
276   else
277      ierr=NF_INQ_DIMLEN(infid,tmpdimid,timelength)
278      if (ierr.ne.NF_NOERR) then
279         stop "Error: Failed to get Time length"
280      else
281         allocate(time(timelength))
282         ierr=NF_GET_VAR_REAL(infid,tmpvarid,time)
283         if (ierr.ne.NF_NOERR) then
284            write(*,*) "Error: Failed reading Time"
285            stop ""
286         endif
287      endif
288   endif
289endif
290
291! altitude
292ierr=NF_INQ_DIMID(infid,"altitude",tmpdimid)
293if (ierr.ne.NF_NOERR) then
294   stop "Error: Failed to get altitude dimension ID"
295else
296   ierr=NF_INQ_VARID(infid,"altitude",tmpvarid)
297   if (ierr.ne.NF_NOERR) then
298      stop "Error: Failed to get altitude length"
299   else
300      ierr=NF_INQ_DIMLEN(infid,tmpdimid,altlength)
301      if (ierr.ne.NF_NOERR) then
302         stop "Error: Failed to get altitude length"
303      else   
304         allocate(alt(altlength))
305         ierr=NF_GET_VAR_REAL(infid,tmpvarid,alt)
306         if (ierr.ne.NF_NOERR) then
307            write(*,*) "Error: Failed reading Altitude"
308            stop ""
309         endif
310      endif
311   endif
312endif
313
314! 1.3.2 Get hybrid coordinates
315
316! Look for hybrid coordinates ap and bp
317allocate(ap(altlength+1))
318allocate(bp(altlength+1))
319ierr=NF_INQ_VARID(infid,"ap",tmpvarid)
320if (ierr.ne.NF_NOERR) then
321   ! If ap not in input file, it is in a diagfi file
322   write(*,*) "Failed to get ap ID from input file ",trim(infile)
323   infile2="diagfi.nc"
324   write(*,*) "Trying file ",trim(infile2)
325   ierr=NF_OPEN(infile2,NF_NOWRITE,infid2)
326   if (ierr.ne.NF_NOERR) then
327      ! If diagfi not found, try diagfi1
328      infile2="diagfi1.nc"
329      write(*,*) "Trying file ",trim(infile2)
330      ierr=NF_OPEN(infile2,NF_NOWRITE,infid2)
331      if (ierr.ne.NF_NOERR) then
332         write(*,*) "Error: Could not open these files. I'll stop."
333         stop ""
334      endif
335   endif
336   ! Get ap and bp (assumed to be in the same file)
337   ! ap
338   ierr=NF_INQ_VARID(infid2,"ap",tmpvarid)
339   if (ierr.ne.NF_NOERR) then
340      write(*,*) "Error: Failed to get ap ID"
341      stop ""
342   endif
343   ierr=NF_GET_VAR_REAL(infid2,tmpvarid,ap)
344   if (ierr.ne.NF_NOERR) then
345      write(*,*) "Error: Failed reading ap"
346      stop ""
347   endif
348   ! bp
349   ierr=NF_INQ_VARID(infid2,"bp",tmpvarid)
350   if (ierr.ne.NF_NOERR) then
351      write(*,*) "Error: Failed to get bp ID"
352      stop ""
353   endif
354   ierr=NF_GET_VAR_REAL(infid2,tmpvarid,bp)
355   if (ierr.ne.NF_NOERR) then
356      write(*,*) "Error: Failed reading bp"
357      stop ""
358   endif
359   ! Close file
360   write(*,*) 'OK, got ap and bp'
361   ierr=NF_CLOSE(infid2)
362else
363   ! Get ap and bp (assumed to be in the same file)
364   ! ap (already got ID)
365   ierr=NF_GET_VAR_REAL(infid,tmpvarid,ap)
366   if (ierr.ne.NF_NOERR) then
367      write(*,*) "Error: Failed reading ap"
368      stop ""
369   endif
370   ! bp
371   ierr=NF_INQ_VARID(infid,"bp",tmpvarid)
372   if (ierr.ne.NF_NOERR) then
373      write(*,*) "Error: Failed to get bp ID"
374      stop ""
375   endif
376   ierr=NF_GET_VAR_REAL(infid,tmpvarid,bp)
377   if (ierr.ne.NF_NOERR) then
378      write(*,*) "Error: Failed reading bp"
379      stop ""
380   endif
381   write(*,*) 'OK, got ap and bp'
382endif
383
384! Look for hybrid coordinates aps and bps
385allocate(aps(altlength))
386allocate(bps(altlength))
387ierr=NF_INQ_VARID(infid,"aps",tmpvarid)
388if (ierr.ne.NF_NOERR) then
389   ! If ap not in input file, it is in a diagfi file
390   write(*,*) "Failed to get ap ID from input file ",trim(infile)
391   infile2="diagfi.nc"
392   write(*,*) "Trying file ",trim(infile2)
393   ierr=NF_OPEN(infile2,NF_NOWRITE,infid2)
394   if (ierr.ne.NF_NOERR) then
395      ! If diagfi not found, try diagfi1
396      infile2="diagfi1.nc"
397      write(*,*) "Trying file ",trim(infile2)
398      ierr=NF_OPEN(infile2,NF_NOWRITE,infid2)
399      if (ierr.ne.NF_NOERR) then
400         write(*,*) "Error: Could not open these files. I'll stop."
401         stop ""
402      endif
403   endif
404   ! Get aps and bps (assumed to be in the same file)
405   ! aps
406   ierr=NF_INQ_VARID(infid2,"aps",tmpvarid)
407   if (ierr.ne.NF_NOERR) then
408      write(*,*) "Error: Failed to get aps ID"
409      stop ""
410   endif
411   ierr=NF_GET_VAR_REAL(infid2,tmpvarid,aps)
412   if (ierr.ne.NF_NOERR) then
413      write(*,*) "Error: Failed reading aps"
414      stop ""
415   endif
416   ! bps
417   ierr=NF_INQ_VARID(infid2,"bps",tmpvarid)
418   if (ierr.ne.NF_NOERR) then
419      write(*,*) "Error: Failed to get bps ID"
420      stop ""
421   endif
422   ierr=NF_GET_VAR_REAL(infid2,tmpvarid,bps)
423   if (ierr.ne.NF_NOERR) then
424      write(*,*) "Error: Failed reading bps"
425      stop ""
426   endif
427   ! Close file
428   write(*,*) 'OK, got aps and bps'
429   ierr=NF_CLOSE(infid2)
430else
431   ! Get aps and bps (assumed to be in the same file)
432   ! aps (already got ID)
433   ierr=NF_GET_VAR_REAL(infid,tmpvarid,aps)
434   if (ierr.ne.NF_NOERR) then
435      write(*,*) "Error: Failed reading aps"
436      stop ""
437   endif
438   ! bps
439   ierr=NF_INQ_VARID(infid,"bps",tmpvarid)
440   if (ierr.ne.NF_NOERR) then
441      write(*,*) "Error: Failed to get bps ID"
442      stop ""
443   endif
444   ierr=NF_GET_VAR_REAL(infid,tmpvarid,bps)
445   if (ierr.ne.NF_NOERR) then
446      write(*,*) "Error: Failed reading bps"
447      stop ""
448   endif
449   write(*,*) 'OK, got aps and bps'
450endif
451
452!===============================================================================
453! 1.4 Read surface pressure for vertical integration
454!===============================================================================
455
456ierr=NF_INQ_VARID(infid,"ps",tmpvarid)
457if (ierr.ne.NF_NOERR) then
458   write(*,*) "Error: Failed to get ps ID"
459   stop ""
460else
461   allocate(ps(lonlength,latlength,timelength))
462   ierr=NF_GET_VAR_REAL(infid,tmpvarid,ps)
463   if (ierr.ne.NF_NOERR) then
464      write(*,*) "Error: Failed reading ps"
465      stop ""
466   endif
467endif
468
469!===============================================================================
470! 1.5 Read atmospheric density to convert concentration to mass
471!===============================================================================
472
473foundrho = .false.
474ierr=NF_INQ_VARID(infid,"rho",tmpvarid)
475if (ierr.ne.NF_NOERR) then
476   write(*,*) "Failed to get rho ID"
477else
478   allocate(rho(lonlength,latlength,altlength,timelength))
479   ierr=NF_GET_VAR_REAL(infid,tmpvarid,rho)
480   if (ierr.ne.NF_NOERR) then
481      write(*,*) "Failed to read rho"
482   else
483      foundrho = .true.
484   endif
485endif
486
487!===============================================================================
488! 1.6 Read/compute or assume MMEAN cst
489!     (for vmr and concentration conversion)
490!===============================================================================
491
492rgp = 8.3143 ! universal gas constant in J mol-1 K-1
493Rmean = 191 ! mean gas constant in J kg-1 K-1
494
495ismmcst = .false.
496foundmm = .false.
497foundpress = .true.
498
499! 1.5.1 Look for mmean
500ierr=NF_INQ_VARID(infid,"mmean",tmpvarid)
501if (ierr.ne.NF_NOERR) then ! can't find mean id
502   write(*,*) "Failed to get mmean ID."
503   write(*,*) "I'll try to compute it from rho, temp and pressure."
504   ! allocate locmm for later computation
505   allocate(locmm(lonlength,latlength,altlength,timelength))
506   ! Look for rho
507   if (foundrho) then
508      write(*,*) "Previously found rho, looking for temp and pressure"
509      ! Look for temp
510      ierr=NF_INQ_VARID(infid,"temp",tmpvarid)
511      if (ierr.ne.NF_NOERR) then
512         write(*,*) "Failed to get temp ID."
513         write(*,*) "Can't compute mmean, I'll assume it is constant."
514         cstmm = rgp/Rmean
515         ismmcst = .true.
516      else
517         allocate(temp(lonlength,latlength,altlength,timelength))
518         ierr=NF_GET_VAR_REAL(infid,tmpvarid,temp)
519         if (ierr.ne.NF_NOERR) then
520            write(*,*) "Failed to read temp."
521            write(*,*) "Can't compute mmean, I'll assume it is constant."
522            cstmm = rgp/Rmean
523            ismmcst = .true.
524         else
525            ! Look for pressure
526            ierr=NF_INQ_VARID(infid,"pressure",tmpvarid)
527            if (ierr.ne.NF_NOERR) then
528               write(*,*) "Failed to get pressure ID."
529               write(*,*) "I'll compute it from aps, bps and ps"
530               ! allocate press for later computation
531               allocate(press(lonlength,latlength,altlength,timelength))
532               foundpress=.false.
533            else
534               allocate(press(lonlength,latlength,altlength,timelength))
535               ierr=NF_GET_VAR_REAL(infid,tmpvarid,press)
536               if (ierr.ne.NF_NOERR) then
537                  write(*,*) "Failed to read pressure."
538                  write(*,*) "I'll compute it from aps, bps and ps."
539                  foundpress=.false.
540               endif
541            endif ! end look for pressure
542         endif
543      endif ! end look for temp
544   else ! not found rho
545      write(*,*) "Failed to read rho."
546      write(*,*) "Can't compute mmean, I'll assume it is constant."
547      cstmm = rgp/Rmean
548      ismmcst = .true.
549   endif ! end if found rho
550else ! if found mmean id, try to read it
551   allocate(locmm(lonlength,latlength,altlength,timelength))
552   ierr=NF_GET_VAR_REAL(infid,tmpvarid,locmm)
553   if (ierr.ne.NF_NOERR) then
554      write(*,*) "Failed to read mmean ID."
555      write(*,*) "I'll try to compute it from rho, temp and pressure."
556      ! Look for rho
557      if (foundrho) then
558         write(*,*) "Previously found rho, looking for temp and pressure"
559         ! Look for temp
560         ierr=NF_INQ_VARID(infid,"temp",tmpvarid)
561         if (ierr.ne.NF_NOERR) then
562            write(*,*) "Failed to get temp ID."
563            write(*,*) "Can't compute mmean, I'll assume it is constant."
564            cstmm = rgp/Rmean
565            ismmcst = .true.
566         else
567            allocate(temp(lonlength,latlength,altlength,timelength))
568            ierr=NF_GET_VAR_REAL(infid,tmpvarid,temp)
569            if (ierr.ne.NF_NOERR) then
570               write(*,*) "Failed to read temp."
571               write(*,*) "Can't compute mmean, I'll assume it is constant."
572               cstmm = rgp/Rmean
573               ismmcst = .true.
574            else
575               ! Look for pressure
576               ierr=NF_INQ_VARID(infid,"pressure",tmpvarid)
577               if (ierr.ne.NF_NOERR) then
578                  write(*,*) "Failed to get pressure ID."
579                  write(*,*) "I'll compute it from aps, bps and ps"
580                  ! allocate press for later computation
581                  allocate(press(lonlength,latlength,altlength,timelength))
582                  foundpress=.false.
583               else
584                  allocate(press(lonlength,latlength,altlength,timelength))
585                  ierr=NF_GET_VAR_REAL(infid,tmpvarid,press)
586                  if (ierr.ne.NF_NOERR) then
587                     write(*,*) "Failed to read pressure."
588                     write(*,*) "I'll compute it from aps, bps and ps."
589                     foundpress=.false.
590                  endif
591               endif ! end look for pressure
592            endif
593         endif ! end look for temp
594      else ! not found rho
595         write(*,*) "Failed to read rho."
596         write(*,*) "Can't compute mmean, I'll assume it is constant."
597         cstmm = rgp/Rmean
598         ismmcst = .true.
599      endif ! end look for rho
600   else
601      foundmm=.true.
602   endif
603endif ! end look for mmean
604
605! 1.5.2 Compute press and mm if required
606if (.not.ismmcst) then
607   if (.not.foundmm) then
608      if (.not.foundpress) then
609         ! Compute press
610         do it=1,timelength
611            do ialt=1,altlength
612               do ilat=1,latlength
613                  do ilon=1,lonlength
614                     press(ilon,ilat,ialt,it)=&
615                        aps(ialt)+bps(ialt)*ps(ilon,ialt,it)
616                  enddo
617               enddo
618            enddo
619         enddo
620         ! end compute press
621      endif ! end if not found press
622
623      ! Compute local atmospheric molecular mass
624      do it=1,timelength
625         do ialt=1,altlength
626            do ilat=1,latlength
627               do ilon=1,lonlength
628                  locmm(ilon,ilat,ialt,it)=&
629                     rgp*rho(ilon,ilat,ialt,it)*&
630                     temp(ilon,ilat,ialt,it)/&
631                     press(ilon,ilat,ialt,it)
632               enddo
633            enddo
634         enddo
635      enddo
636      ! end compute local atmospheric molecular mass
637   endif ! end if not found mm
638endif ! end if mm not cst
639
640!==============================================================================
641! 1.6. Create and initialize output file
642!==============================================================================
643
644filename=infile(1:len_trim(infile)-3)//"_col.nc"
645write(*,*)'The output file has name:'
646write(*,*) filename
647
648! Initialize output file's lat,lon,alt and time dimensions
649call initiate(filename,lat,lon,alt,time,nout,&
650     latdimout,londimout,altdimout,timedimout)
651! ! ! Initialize output file's aps,bps variables
652! call init2(infid,lonlength,latlength,altlength,&
653!      nout,londimout,latdimout,altdimout, interlayerdimout)
654
655!==============================================================================
656! 2. PROCESS DATA
657!==============================================================================
658
659!This is where you can work on data
660
661! Constants
662g = 3.72
663nav = 6.023d23
664
665!=============================================================
666allocate(tracer(lonlength,latlength,altlength,timelength))
667allocate(tracerCol(lonlength,latlength,timelength))
668
669do i=1,nbvar
670   ! Read tracer
671   ierr=NF_INQ_VARID(infid,var(i),tmpvarid)
672   if (ierr.ne.NF_NOERR) then
673      write(*,*) "Error: Failed reading ID for ", trim(var(i))
674      stop ""
675   endif
676   ierr=NF_GET_VAR_REAL(infid,tmpvarid,tracer)
677   if (ierr.ne.NF_NOERR) then
678      write(*,*) "Error: Failed reading ", trim(var(i))
679      stop ""
680   endif
681   ! check tracer unit
682   ierr=NF_GET_ATT_TEXT(infid,tmpvarid,'units',tmpunit)
683   ismmr = (tmpunit.eq.'kg/kg').or.(tmpunit.eq.'kg.kg-1')
684   isvmr = (tmpunit.eq.'mol/mol').or.(tmpunit.eq.'mol.mol-1')
685   isco = (tmpunit.eq.'cm-3').or.(tmpunit.eq.'m-3')
686   ! Compute column in mmr case
687   if (ismmr) then
688      do ilat=1,latlength
689         do ilon=1,lonlength
690            do it=1,timelength
691               tracerMass = 0.0
692               do ialt=1,altlength
693                  incMass = ((ap(ialt) - ap(ialt + 1)) + (bp(ialt) - bp(ialt + 1))*ps(ilon,ilat,it))/g
694                  tracerMass = tracerMass + incMass*tracer(ilon,ilat,ialt,it)
695               enddo
696               tracerCol(ilon,ilat,it) = tracerMass
697            enddo
698         enddo
699      enddo
700   ! Compute column in vmr case
701   else if (isvmr) then
702      write(*,*) trim(var(i)), " is a vmr tracer"
703      ! Get current specie's molecular mass
704      call getmm(var(i),tmpmm,foundmm)
705      if (.not.foundmm) then
706         write(*,*) "Cannot compute column from vmr variable ", trim(var(i))
707         write(*,*) "Check unit, name, or add it to getmm. I'll stop here"
708         stop ""
709      else
710         write(*,*) "Getting molecular mass", tmpmm, " kg/mol"
711         do ilat=1,latlength
712            do ilon=1,lonlength
713               do it=1,timelength
714                  tracerMass = 0.0
715                  do ialt=1,altlength
716                     incMass = ((ap(ialt) - ap(ialt + 1)) + (bp(ialt) - bp(ialt + 1))*ps(ilon,ilat,it))/g
717                     ! Handle cst mm case
718                     if (ismmcst) then
719                        tracerMass = tracerMass + incMass*tracer(ilon,ilat,ialt,it)*&
720                           tmpmm/cstmm
721                     ! Handle local mm case
722                     else
723                        ! mmean is stored in g/mol so we multiply by 1.e3
724                        tracerMass = tracerMass + incMass*tracer(ilon,ilat,ialt,it)*&
725                           tmpmm/locmm(ilon,ilat,ialt,it)*1.e3
726                     endif
727                  enddo
728                  tracerCol(ilon,ilat,it) = tracerMass
729               enddo
730            enddo
731         enddo
732      endif
733   ! Compute column in concentration case
734   else if(isco) then   
735      write(*,*) trim(var(i)), " is a num tracer"
736      write(*,*) "Checking rho is available..."
737      if (.not.foundrho) then
738         write(*,*) "Error: rho not available. Can't proceed"
739         stop ""
740      endif
741      ! Get current specie's molecular mass
742      call getmm(var(i),tmpmm,foundmm)
743      if (.not.foundmm) then
744         write(*,*) "Cannot compute column from concentration variable ", trim(var(i))
745         write(*,*) "Check unit, name, or add it to getmm. I'll stop here"
746         stop ""
747      else
748         write(*,*) "Getting molecular mass", tmpmm, " kg/mol"
749         do ilat=1,latlength
750            do ilon=1,lonlength
751               do it=1,timelength
752                  tracerMass = 0.0
753                  do ialt=1,altlength
754                     incMass = ((ap(ialt) - ap(ialt + 1)) + (bp(ialt) - bp(ialt + 1))*ps(ilon,ilat,it))/g
755                     tracerMass = tracerMass + incMass/rho(ilon,ilat,ialt,it)*tracer(ilon,ilat,ialt,it)*tmpmm/nav
756                  enddo
757                  ! Is it in S.I. ?
758                  if (tmpunit.eq.'m-3') then
759                     tracerCol(ilon,ilat,it) = tracerMass
760                  else
761                     tracerCol(ilon,ilat,it) = tracerMass*1.d6
762                  endif
763               enddo
764            enddo
765         enddo
766      endif
767   ! If none of previous case matches, skip
768   else
769      write(*,*) trim(var(i)), " is neither mmr, vmr or concentration"
770      write(*,*) "Will skip"
771   endif ! end of ismmr, isvmr, isco
772
773   ! Define variable to write
774   tmpvarname = trim(var(i))//"_col"
775   call def_var(nout,tmpvarname,tmpvarname,"kg/m2",3,&
776               (/londimout,latdimout,timedimout/),tracvarout,ierr)
777   if (ierr.ne.NF_NOERR) then
778      write(*,*) 'Error, could not define variable ', trim(var(i))
779      stop ""
780   endif
781   ! Write in the output file
782   ierr=NF_PUT_VAR_REAL(nout,tracvarout,tracerCol)
783   if (ierr.ne.NF_NOERR) then
784      write(*,*)'Error, Failed to write variable ', trim(var(i))
785      stop ""
786   endif
787enddo ! end of do nbvar
788
789
790!----------------------------------
791! Close input file
792ierr=NF_CLOSE(infid)
793
794! Close output file
795ierr=NF_CLOSE(nout)
796
797contains
798
799!******************************************************************************
800Subroutine initiate (filename,lat,lon,alt,time,&
801                     nout,latdimout,londimout,altdimout,timedimout)
802!==============================================================================
803! Purpose:
804! Create and initialize a data file (NetCDF format)
805! Include grid and time (lon,lat,alt,time)
806!==============================================================================
807! Remarks:
808! The NetCDF file (created in this subroutine) remains open
809!==============================================================================
810
811implicit none
812
813include "netcdf.inc" ! NetCDF definitions
814
815!==============================================================================
816! Arguments:
817!==============================================================================
818character (len=*), intent(in):: filename
819! filename(): the file's name
820real, dimension(:), intent(in):: lat
821! lat(): latitude
822real, dimension(:), intent(in):: lon
823! lon(): longitude
824real, dimension(:), intent(in):: alt
825! alt(): altitude
826real, dimension(:), intent(in):: time
827! time(): Time
828integer, intent(out):: nout
829! nout: [netcdf] file ID
830integer, intent(out):: latdimout
831! latdimout: [netcdf] lat() (i.e.: latitude)  ID
832integer, intent(out):: londimout
833! londimout: [netcdf] lon()  ID
834integer, intent(out):: altdimout
835! altdimout: [netcdf] alt()  ID
836integer, intent(out):: timedimout
837! timedimout: [netcdf] "Time"  ID
838
839!==============================================================================
840! Local variables:
841!==============================================================================
842integer :: nvarid,ierr
843! nvarid: [netcdf] ID of a variable
844! ierr: [netcdf]  return error code (from called subroutines)
845
846!==============================================================================
847! 1. Create (and open) output file
848!==============================================================================
849write(*,*) "creating "//trim(adjustl(filename))//'...'
850ierr = NF_CREATE(filename,NF_CLOBBER,nout)
851! NB: setting NF_CLOBBER mode means that it's OK to overwrite an existing file
852if (ierr.NE.NF_NOERR) then
853   WRITE(*,*)'ERROR: Impossible to create the file.'
854   stop ""
855endif
856
857!==============================================================================
858! 2. Define/write "dimensions" and get their IDs
859!==============================================================================
860
861ierr = NF_DEF_DIM(nout, "latitude", size(lat), latdimout)
862ierr = NF_DEF_DIM(nout, "longitude", size(lon), londimout)
863ierr = NF_DEF_DIM(nout, "altitude", size(alt), altdimout)
864ierr = NF_DEF_DIM(nout, "Time", size(time), timedimout)
865
866! End netcdf define mode
867ierr = NF_ENDDEF(nout)
868
869!==============================================================================
870! 3. Write "Time" (attributes)
871!==============================================================================
872
873call def_var(nout,"Time","Time","days since 0000-00-0 00:00:00",1,&
874             (/timedimout/),nvarid,ierr)
875
876#ifdef NC_DOUBLE
877ierr = NF_PUT_VAR_DOUBLE(nout,nvarid,time)
878#else
879ierr = NF_PUT_VAR_REAL(nout,nvarid,time)
880#endif
881
882!==============================================================================
883! 4. Write "latitude" (data and attributes)
884!==============================================================================
885
886call def_var(nout,"latitude","latitude","degrees_north",1,&
887             (/latdimout/),nvarid,ierr)
888
889#ifdef NC_DOUBLE
890ierr = NF_PUT_VAR_DOUBLE(nout,nvarid,lat)
891#else
892ierr = NF_PUT_VAR_REAL(nout,nvarid,lat)
893#endif
894
895!==============================================================================
896! 4. Write "longitude" (data and attributes)
897!==============================================================================
898
899call def_var(nout,"longitude","East longitude","degrees_east",1,&
900             (/londimout/),nvarid,ierr)
901
902#ifdef NC_DOUBLE
903ierr = NF_PUT_VAR_DOUBLE(nout,nvarid,lon)
904#else
905ierr = NF_PUT_VAR_REAL(nout,nvarid,lon)
906#endif
907
908!==============================================================================
909! 4. Write "altitude" (data and attributes)
910!==============================================================================
911
912call def_var(nout,"altitude","Altitude","km",1,&
913             (/altdimout/),nvarid,ierr)
914
915#ifdef NC_DOUBLE
916ierr = NF_PUT_VAR_DOUBLE(nout,nvarid,alt)
917#else
918ierr = NF_PUT_VAR_REAL(nout,nvarid,alt)
919#endif
920
921end Subroutine initiate
922
923subroutine def_var(nid,name,title,units,nbdim,dim,nvarid,ierr)
924!==============================================================================
925! Purpose: Write a variable (i.e: add a variable to a dataset)
926! called "name"; along with its attributes "title", "units"...
927! to a file (following the NetCDF format)
928!==============================================================================
929! Remarks:
930! The NetCDF file must be open
931!==============================================================================
932
933implicit none
934
935include "netcdf.inc" ! NetCDF definitions
936
937!==============================================================================
938! Arguments:
939!==============================================================================
940integer, intent(in) :: nid
941! nid: [netcdf] file ID #
942character (len=*), intent(in) :: name
943! name(): [netcdf] variable's name
944character (len=*), intent(in) :: title
945! title(): [netcdf] variable's "title" attribute
946character (len=*), intent(in) :: units
947! unit(): [netcdf] variable's "units" attribute
948integer, intent(in) :: nbdim
949! nbdim: number of dimensions of the variable
950integer, dimension(nbdim), intent(in) :: dim
951! dim(nbdim): [netcdf] dimension(s) ID(s)
952integer, intent(out) :: nvarid
953! nvarid: [netcdf] ID # of the variable
954integer, intent(out) :: ierr
955! ierr: [netcdf] subroutines returned error code
956
957! Switch to netcdf define mode
958ierr=NF_REDEF(nid)
959
960! Insert the definition of the variable
961#ifdef NC_DOUBLE
962ierr=NF_DEF_VAR(nid,adjustl(name),NF_DOUBLE,nbdim,dim,nvarid)
963#else
964ierr=NF_DEF_VAR(nid,adjustl(name),NF_FLOAT,nbdim,dim,nvarid)
965#endif
966
967! Write the attributes
968ierr=NF_PUT_ATT_TEXT(nid,nvarid,"title",len_trim(adjustl(title)),adjustl(title))
969ierr=NF_PUT_ATT_TEXT(nid,nvarid,"units",len_trim(adjustl(units)),adjustl(units))
970
971! End netcdf define mode
972ierr=NF_ENDDEF(nid)
973
974end subroutine def_var
975
976!******************************************************************************
977subroutine getmm(invarname, tmpmm, foundmm)
978!==============================================================================
979! Purpose:
980! Get the molar mass tmpmm coresponding to a given varname
981!==============================================================================
982! Remarks:
983! - varname is assumed to have the format *_name_* or *_name
984!   (based on the usual format vmr_..., num_...)
985! - Feel free to add any missing specie
986!==============================================================================
987
988implicit none
989
990!==============================================================================
991! Arguments:
992!==============================================================================
993character(len=*), intent(in) :: invarname ! name of the input variable
994real, intent(out) :: tmpmm ! molar mass to output
995logical, intent(out) :: foundmm ! did we find a correspondance ?
996
997!==============================================================================
998! Local variables
999!==============================================================================
1000integer :: first, last, pos, length, count
1001! first: beginning of specie name, last: end of specie name
1002! pos: cursor position in string varname. length: varname length
1003! count: number of encountered separator "_" (2 at most)
1004character(len=len_trim(invarname)) :: varname ! name of the specie
1005
1006!==============================================================================
1007! Find specie name
1008!==============================================================================
1009
1010length = len_trim(invarname)
1011first = 1
1012last = length
1013pos=1
1014count=0
1015! Cursor goes through varname
1016do while (pos < length)
1017   if ((invarname(pos:pos).eq.'_').and.(count.eq.0)) then
1018      ! First "_" locates the beginning of the specie's name
1019      count = count + 1
1020      first = pos + 1
1021   else if ((invarname(pos:pos).eq.'_').and.(count.eq.1)) then
1022      ! Second "_" locates the end (default is length)
1023      count = count + 1
1024      last = pos - 1
1025   endif
1026   pos = pos + 1
1027enddo
1028varname = invarname(first:last)
1029
1030!==============================================================================
1031! Check name among all possibilities
1032!==============================================================================
1033
1034foundmm=.false.
1035
1036if (varname.eq."co2") then
1037   tmpmm=44.
1038   foundmm=.true.
1039endif
1040if (varname.eq."co") then
1041   tmpmm=28.
1042   foundmm=.true.
1043endif
1044if (varname.eq."o") then
1045   tmpmm=16.
1046   foundmm=.true.
1047endif
1048if (varname.eq."o1d") then
1049   tmpmm=16.
1050   foundmm=.true.
1051endif
1052if (varname.eq."o2") then
1053   tmpmm=32.
1054   foundmm=.true.
1055endif
1056if (varname.eq."o3") then
1057   tmpmm=48.
1058   foundmm=.true.
1059endif
1060if (varname.eq."h") then
1061   tmpmm=1.
1062   foundmm=.true.
1063endif
1064if (varname.eq."h2") then
1065   tmpmm=2.
1066   foundmm=.true.
1067endif
1068if (varname.eq."oh") then
1069   tmpmm=17.
1070   foundmm=.true.
1071endif
1072if (varname.eq."ho2") then
1073   tmpmm=33.
1074   foundmm=.true.
1075endif
1076if (varname.eq."h2o2") then
1077   tmpmm=34.
1078   foundmm=.true.
1079endif
1080if (varname.eq."n2") then
1081   tmpmm=28.
1082   foundmm=.true.
1083endif
1084if (varname.eq."ch4") then
1085   tmpmm=16.
1086   foundmm=.true.
1087endif
1088if (varname.eq."ar") then
1089   tmpmm=40.
1090   foundmm=.true.
1091endif
1092if (varname.eq."n") then
1093   tmpmm=14.
1094   foundmm=.true.
1095endif
1096if (varname.eq."no") then
1097   tmpmm=30.
1098   foundmm=.true.
1099endif
1100if (varname.eq."no2") then
1101   tmpmm=46.
1102   foundmm=.true.
1103endif
1104if (varname.eq."n2d") then
1105   tmpmm=28.
1106   foundmm=.true.
1107endif
1108if (varname.eq."he") then
1109   tmpmm=4.
1110   foundmm=.true.
1111endif
1112if (varname.eq."co2plus") then
1113   tmpmm=44.
1114   foundmm=.true.
1115endif
1116if (varname.eq."oplus") then
1117   tmpmm=16.
1118   foundmm=.true.
1119endif
1120if (varname.eq."o2plus") then
1121   tmpmm=32.
1122   foundmm=.true.
1123endif
1124if (varname.eq."coplus") then
1125   tmpmm=28.
1126   foundmm=.true.
1127endif
1128if (varname.eq."cplus") then
1129   tmpmm=12.
1130   foundmm=.true.
1131endif
1132if (varname.eq."nplus") then
1133   tmpmm=14.
1134   foundmm=.true.
1135endif
1136if (varname.eq."noplus") then
1137   tmpmm=30.
1138   foundmm=.true.
1139endif
1140if (varname.eq."n2plus") then
1141   tmpmm=28.
1142   foundmm=.true.
1143endif
1144if (varname.eq."hplus") then
1145   tmpmm=1.
1146   foundmm=.true.
1147endif
1148if (varname.eq."hco2plus") then
1149   tmpmm=45.
1150   foundmm=.true.
1151endif
1152if (varname.eq."hcoplus") then
1153   tmpmm=29.
1154   foundmm=.true.
1155endif
1156if (varname.eq."h2plus") then
1157   tmpmm=18.
1158   foundmm=.true.
1159endif
1160if (varname.eq."h3oplus") then
1161   tmpmm=19.
1162   foundmm=.true.
1163endif
1164if (varname.eq."ohplus") then
1165   tmpmm=17.
1166   foundmm=.true.
1167endif
1168if (varname.eq."elec") then
1169   tmpmm=1./1822.89
1170   foundmm=.true.
1171endif
1172if (varname.eq."h2ovap") then ! tmp
1173   tmpmm=18.
1174   foundmm=.true.
1175endif
1176if (varname.eq."h2o") then ! tmp
1177   tmpmm=18.
1178   foundmm=.true.
1179endif
1180if (varname.eq."hdovap") then
1181   tmpmm=19.
1182   foundmm=.true.
1183endif
1184if (varname.eq."od") then
1185   tmpmm=18.
1186   foundmm=.true.
1187endif
1188if (varname.eq."d") then
1189   tmpmm=2.
1190   foundmm=.true.
1191endif
1192if (varname.eq."hd") then
1193   tmpmm=3.
1194   foundmm=.true.
1195endif
1196if (varname.eq."do2") then
1197   tmpmm=34.
1198   foundmm=.true.
1199endif
1200if (varname.eq."hdo2") then
1201   tmpmm=35.
1202   foundmm=.true.
1203endif
1204if (varname.eq."co2ice") then
1205   tmpmm=44.
1206   foundmm=.true.
1207endif
1208if (varname.eq."h2oice") then
1209   tmpmm=18.
1210   foundmm=.true.
1211endif
1212if (varname.eq."hdoice") then
1213   tmpmm=19.
1214   foundmm=.true.
1215endif
1216if (varname.eq."Ar_N2") then
1217   tmpmm=30.
1218   foundmm=.true.
1219endif       
1220
1221if (foundmm) then
1222   ! Back in kg
1223   tmpmm = tmpmm/1.e3
1224else
1225   write(*,*) "Found no molar mass corresponding to variable ", trim(varname)
1226endif
1227
1228end subroutine getmm
1229
1230end program
Note: See TracBrowser for help on using the repository browser.