source: trunk/LMDZ.MARS/util/aeroptical.F90 @ 2651

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

Mars GCM utilities:
Minor fixes to run with picky gfortran 10.3.0 which requires one element arrays
(rather than scalars) when calling NetCDF routines, andf that stop statements
should not be followed by strings. While at it replaced tabs with spaces.
EM
File size: 59.8 KB
Line 
1program aeroptical
2! program for computation of aerosol opacities
3! author : Antoine Bierjon, April-May 2020
4!
5!===============================================================================
6!     PREFACE
7!===============================================================================
8! This program takes as inputs a GCM output file (diagfi,stats,concat) that
9! contains:
10!      - the Mass Mixing Ratios of dust (dustq), water ice (h2o_ice),
11!                               stormdust (rdsdustq) & topdust (topdustq)
12!      - their effective radius (reffdust, reffice(*),
13!                           reffstormdust, refftopdust)
14!      - the atmospheric density (rho)
15! as well as a wavelength to calibrate the opacities, and the directory
16! containing the ASCII files of the aerosols' optical properties.
17!
18! It outputs a netcdf file containing the opacities [1/km] of the dust,
19! water ice and stormdust aerosols calibrated at the prescribed wavelength.
20! The type of opacity (extinction or absorption) is given by the user.
21!
22! (*) due to a high uncertainty of reffice in the GCM, the value is asked
23! directly to the user (if the user returns 0, then the program reads the GCM
24! file to get reffice)
25!
26! NOTA BENE:
27! 1) if one wanted to add another aerosol to compute, one should look for
28! the comments + NEW AER? that are disseminated all along the code to indicate
29! the parts of the code that should be modified.
30! 2) another enhancement of this program could be the possibility to read its
31! own product files and recalibrate the opacities at another wavelength
32!===============================================================================
33
34
35use netcdf
36
37!===============================================================================
38! Variable declarations
39!===============================================================================
40
41implicit none ! for no implicitly typed variables
42
43
44! GCM file
45character(len=128) :: gcmfile                 ! name of the netcdf GCM input file
46integer :: gcmfid                             ! [netcdf] GCM input file ID number
47integer :: ierr                               ! [netcdf] subroutine returned error code
48character(len=256) :: error_text              ! text to display in case of error
49integer :: lonlen,latlen,altlen,timelen       ! nb of grid points along longitude,latitude,altitude,time
50integer :: GCM_layers                         ! number of GCM layers
51integer :: layerdimout,interlayerdimout       ! "GCM_layers" and "GCM_layers+1" IDs
52
53logical,dimension(:),allocatable :: aerok     ! to know if the needed fields are in the GCM file
54integer,dimension(:),allocatable :: mmrvarid  ! stores a MMR variable ID number
55integer,dimension(:),allocatable :: reffvarid ! stores a reff variable ID number
56integer :: tmpvarid                           ! temporarily stores a variable ID number
57real,dimension(:,:,:,:),allocatable :: mmr    ! aerosols mass mixing ratios [kg/kg]
58real,dimension(:,:,:,:),allocatable :: reff   ! aerosols effective radii [m]
59real :: reffwice_val                          ! value given by the user for reffwice (0 if read in the GCM file) [m]
60real,dimension(:,:,:,:),allocatable :: rho    ! atmospheric density [kg.m-3]
61integer :: naerkind                           ! nb of aerosols to consider
62integer :: iaer                               ! aerosol kind index (1=dust,2=water ice,3=stormdust,4=topdust) ! + NEW AER?
63integer :: ilon,ilat,ialt,it                  ! Loop indices for coordinates
64
65
66! Output file
67character(len=128) :: outfile                  ! name of the netcdf output file
68integer :: outfid                              ! [netcdf] file ID number
69integer :: latdimout,londimout,altdimout,timedimout
70                                               ! latdimout: stores latitude dimension ID number
71                                               ! londimout: stores longitude dimension ID number
72                                               ! altdimout: stores altitude dimension ID number
73                                               ! timedimout: stores time dimension ID number
74integer :: tmpvaridout                         ! temporarily stores a variable ID number
75
76real :: wvl_val                                ! reference wavelength of the output opacities (given by the user) [m]
77character(len=3) :: opatype                    ! opacity type : extinction (ext) or absorption (abs) (given by the user)
78integer :: varshape(4)                         ! stores a variable shape (of dimensions' IDs)
79character(len=16) :: tmpvarname                ! temporarily stores a variable name
80real,dimension(:,:,:,:),allocatable :: opa_aer ! Aerosols opacities [1/km]
81real,dimension(:),allocatable :: missval       ! Value to put in outfile when the reff is out of bounds
82                                               ! or when there is a mising value in input file
83
84
85! Optical properties (read in external ASCII files)
86character(len=256) :: datadir              ! directory containing the ASCII files
87character(len=128) :: optpropfile          ! name of the ASCII optical properties file
88logical :: file_ok                         ! to know if the file can be opened
89integer :: file_unit = 60                  ! ASCII file ID
90
91integer :: jfile                           ! ASCII file scan index
92logical :: endwhile                        ! Reading loop boolean
93character(len=132) :: scanline             ! ASCII file scanning line
94integer :: read_ok                         ! to know if the line is well read
95                   
96integer :: nwvl                            ! Number of wavelengths in the domain (VIS or IR)
97integer :: nsize                           ! Number of particle sizes stored in the file
98integer :: iwvl,isize                      ! Wavelength and Particle size loop indices
99real,dimension(:),allocatable :: wvl       ! Wavelength axis [m]
100real,dimension(:),allocatable :: radiusdyn ! Particle size axis [m]
101real,dimension(:,:),allocatable :: Qext    ! Extinction efficiency coefficient [/]
102real,dimension(:,:),allocatable :: omeg    ! Single scattering albedo [/]
103
104
105! Opacity computation
106integer :: iwvl1,iwvl2,isize1,isize2     ! Wavelength and Particle size indices for the interpolation
107real :: coeff                            ! Interpolation coefficient
108real,dimension(2) :: reffQext            ! Qext after reff interpolation
109real :: Qext_val                         ! Qext after both interpolations
110real,dimension(2) :: reffomeg            ! omeg after reff interpolation
111real :: omeg_val                         ! omeg after both interpolations
112real,dimension(:),allocatable :: rho_aer ! Density of the aerosols [kg.m-3]
113
114!===============================================================================
115! 0. USER INPUTS
116!===============================================================================
117write(*,*) "Welcome in the aerosol opacities' computation program !"
118write(*,*)
119
120! Ask the GCM file name
121WRITE(*,*) "Enter an input file name (GCM diagfi/stats/concat) :"
122READ(*,*) gcmfile
123
124! Ask the reffwice to the user
125write(*,*)""
126write(*,*) "The water ice effective radius computed by the GCM is known to be a bit inaccurate."
127write(*,*) "Which value do you want to use for it (in meters) ?"
128write(*,*) "(put 0 if you still want the program to read the value in "//trim(gcmfile)//")"
129READ(*,*) reffwice_val
130
131! Ask the wavelength to the user
132write(*,*)""
133WRITE(*,*) "Value of the reference wavelength for the opacities' computation (in meters) ?"
134READ(*,*) wvl_val
135
136! Ask the directory containing the optical properties files
137write(*,*)""
138WRITE(*,*) "In which directory should we look for the optical properties' files ?"
139READ(*,'(a)') datadir
140
141! Ask the type of opacity that has to be computed
142write(*,*)""
143WRITE(*,*) "Do you want to compute extinction or absorption opacities ? (ext/abs)"
144READ(*,*) opatype
145if ((trim(opatype).ne."ext").and.(trim(opatype).ne."abs")) then
146  write(*,*) "opatype = "//opatype
147  write(*,*) "Error: the opacity type must be ext or abs"
148  stop
149endif
150
151!===============================================================================
152! 1. GCM FILE & OUTPUT FILE INITIALIZATIONS
153!===============================================================================
154
155!==========================================================================
156! 1.1 GCM file opening & dimensions - Output file initializations
157!==========================================================================
158
159!==========================================================================
160! --> 1.1.1 Open the netcdf GCM file given by the user
161
162ierr=nf90_open(gcmfile,nf90_nowrite,gcmfid) ! nowrite mode=the program can only read the opened file
163error_text="Error: could not open file "//trim(gcmfile)
164call status_check(ierr,error_text)
165
166!==========================================================================
167! --> 1.1.2 Creation of the outfile
168if (opatype.eq."ext") then
169  outfile=gcmfile(1:index(gcmfile, ".nc")-1)//"_OPAext.nc"
170else ! opatype.eq."abs"
171  outfile=gcmfile(1:index(gcmfile, ".nc")-1)//"_OPAabs.nc"
172endif
173
174 
175ierr=NF90_CREATE(outfile,IOR(nf90_clobber,nf90_64bit_offset),outfid)
176! NB: clobber mode=overwrite any dataset with the same file name !
177! 64bit_offset enables the creation of large netcdf files with variables>2GB
178error_text="Error: could not create file "//trim(outfile)
179call status_check(ierr,error_text)
180write(*,*)"";WRITE(*,*)"Output file is: ",trim(outfile);write(*,*)""
181
182!==========================================================================
183! --> 1.1.3 Get the dimensions and create them in the output file
184
185! Initialize output file's lat,lon,alt,time and controle dimensions
186call inidims(gcmfile,gcmfid,outfile,outfid,&
187             lonlen,latlen,altlen,timelen,&
188             latdimout,londimout,altdimout,timedimout,&
189             GCM_layers,layerdimout,interlayerdimout)
190
191! Initialize output file's aps,bps,ap,bp and phisinit variables
192call init2(gcmfid,lonlen,latlen,altlen,GCM_layers,&
193           outfid,londimout,latdimout,altdimout,&
194           layerdimout,interlayerdimout)
195
196!==========================================================================
197! 1.2 GCM aerosols
198!==========================================================================
199
200! Number of aerosols
201naerkind = 4 ! dust, water ice, stormdust, topdust ! + NEW AER?
202
203! Variables length allocation
204allocate(mmrvarid(naerkind))
205allocate(reffvarid(naerkind))
206
207! Initialize missing value
208allocate(missval(naerkind))
209missval(:)=1.e+20
210
211! Initialize aerok to .true. for all aerosols
212allocate(aerok(naerkind))
213aerok(:)=.true.
214
215!==========================================================================
216! --> 1.2.1 DUST
217
218! DUST MASS MIXING RATIO
219! Check that the GCM file contains that variable
220ierr=nf90_inq_varid(gcmfid,"dustq",mmrvarid(1))
221error_text="Failed to find variable dustq in "//trim(gcmfile)&
222            //" We'll skip the dust aerosol."
223if (ierr.ne.nf90_noerr) then
224  write(*,*)trim(error_text)
225  aerok(1)=.false.
226endif
227
228! DUST EFFECTIVE RADIUS
229if (aerok(1)) then
230  ! Check that the GCM file contains that variable
231  ierr=nf90_inq_varid(gcmfid,"reffdust",reffvarid(1))
232  error_text="Failed to find variable reffdust in "//trim(gcmfile)&
233              //" We'll skip the dust aerosol."
234  if (ierr.ne.nf90_noerr) then
235    write(*,*)trim(error_text)
236    aerok(1)=.false.
237  endif
238endif
239
240!==========================================================================
241! --> 1.2.2 WATER ICE
242
243! WATER ICE MASS MIXING RATIO
244! Check that the GCM file contains that variable
245ierr=nf90_inq_varid(gcmfid,"h2o_ice",mmrvarid(2))
246error_text="Failed to find variable h2o_ice in "//trim(gcmfile)&
247            //" We'll skip the water ice aerosol."
248if (ierr.ne.nf90_noerr) then
249  write(*,*)trim(error_text)
250  aerok(2)=.false.
251endif
252
253! WATER ICE EFFECTIVE RADIUS
254if (aerok(2)) then 
255  IF (reffwice_val.eq.0) THEN
256    ! Check that the GCM file contains that variable
257    ierr=nf90_inq_varid(gcmfid,"reffice",reffvarid(2))
258    error_text="Failed to find variable reffice in "//trim(gcmfile)&
259                //" We'll skip the water ice aerosol."
260    if (ierr.ne.nf90_noerr) then
261      write(*,*)trim(error_text)
262      aerok(2)=.false.
263    endif
264  ENDIF
265endif
266
267!==========================================================================
268! --> 1.2.3 STORMDUST
269
270! STORMDUST MASS MIXING RATIO
271! Check that the GCM file contains that variable
272ierr=nf90_inq_varid(gcmfid,"rdsdustq",mmrvarid(3))
273error_text="Failed to find variable rdsdustq in "//trim(gcmfile)&
274            //" We'll skip the stormdust aerosol."
275if (ierr.ne.nf90_noerr) then
276  write(*,*)trim(error_text)
277  aerok(3)=.false.
278endif
279
280! STORMDUST EFFECTIVE RADIUS
281if (aerok(3)) then
282  ! Check that the GCM file contains that variable
283  ierr=nf90_inq_varid(gcmfid,"reffstormdust",reffvarid(3))
284  error_text="Failed to find variable reffstormdust in "//trim(gcmfile)&
285              //" We'll skip the stormdust aerosol."
286  if (ierr.ne.nf90_noerr) then
287    write(*,*)trim(error_text)
288    aerok(3)=.false.
289  endif
290endif
291
292!==========================================================================
293! --> 1.2.4 TOPDUST
294
295! TOPDUST MASS MIXING RATIO
296! Check that the GCM file contains that variable
297ierr=nf90_inq_varid(gcmfid,"topdustq",mmrvarid(4))
298error_text="Failed to find variable topdustq in "//trim(gcmfile)&
299            //" We'll skip the topdust aerosol."
300if (ierr.ne.nf90_noerr) then
301  write(*,*)trim(error_text)
302  aerok(4)=.false.
303endif
304
305! TOPDUST EFFECTIVE RADIUS
306if (aerok(4)) then
307  ! Check that the GCM file contains that variable
308  ierr=nf90_inq_varid(gcmfid,"refftopdust",reffvarid(4))
309  error_text="Failed to find variable refftopdust in "//trim(gcmfile)&
310              //" We'll skip the topdust aerosol."
311  if (ierr.ne.nf90_noerr) then
312    write(*,*)trim(error_text)
313    aerok(4)=.false.
314  endif
315endif
316
317!==========================================================================
318! --> 1.2.5 + NEW AER?
319
320
321! Check if there is still at least one aerosol to compute
322IF (.NOT.ANY(aerok)) THEN ! + NEW AER?
323  write(*,*) "No aerosol among [dust/water ice/stormdust/topdust] was found in the file. Better stop now..."
324  stop
325ENDIF
326
327!==========================================================================
328! 1.3 GCM ATMOSPHERIC DENSITY
329!==========================================================================
330
331! Check that the GCM file contains that variable
332ierr=nf90_inq_varid(gcmfid,"rho",tmpvarid)
333error_text="Failed to find variable rho in "//trim(gcmfile)&
334            //" We need it to compute the opacity [1/km]."
335call status_check(ierr,error_text)
336! Length allocation for each dimension of the 4D variable
337allocate(rho(lonlen,latlen,altlen,timelen))
338! Load dataset
339ierr=nf90_get_var(gcmfid,tmpvarid,rho)
340error_text="Failed to load atmospheric density"
341call status_check(ierr,error_text)
342write(*,*) "Atmospheric density rho loaded from "//trim(gcmfile)
343
344
345!==========================================================================
346! 1.4 Output variable's initializations and datasets loading
347!==========================================================================
348! Define the ID shape of the output variables
349varshape(1)=londimout
350varshape(2)=latdimout
351varshape(3)=altdimout
352varshape(4)=timedimout
353
354! Fill the aerosol density array
355allocate(rho_aer(naerkind))
356rho_aer(1)=2500. ! dust
357rho_aer(2)=920.  ! water ice
358rho_aer(3)=2500. ! stormdust
359rho_aer(4)=2500. ! topdust
360! + NEW AER?
361
362DO iaer = 1, naerkind ! Loop on aerosol kind
363  if (aerok(iaer)) then ! check if this aerosol can be computed
364    write(*,*)""
365    SELECT CASE(iaer)
366    CASE(1) ! DUST
367      write(*,*)"Computing dust opacity..."
368    CASE(2) ! WATER ICE
369      write(*,*)"Computing water ice opacity..."
370    CASE(3) ! STORMDUST
371      write(*,*)"Computing stormdust opacity..."
372    CASE(4) ! TOPDUST
373      write(*,*)"Computing topdust opacity..."
374! + NEW AER?
375    END SELECT ! iaer
376   
377    ! Length allocation of the input variables
378    ALLOCATE(mmr(lonlen,latlen,altlen,timelen))
379    ALLOCATE(reff(lonlen,latlen,altlen,timelen))
380   
381    ! Datasets loading
382     ! MMR
383    ierr=NF90_GET_VAR(gcmfid,mmrvarid(iaer),mmr(:,:,:,:))
384    error_text="Failed to load mass mixing ratio"
385    call status_check(ierr,error_text)
386    write(*,*) "Mass mixing ratio loaded from "//trim(gcmfile)
387    ! Get missing value
388    ierr=nf90_get_att(gcmfid,mmrvarid(1),"missing_value",missval(1))
389    if (ierr.ne.nf90_noerr) then
390      missval(1) = 1.e+20
391    endif
392   
393     ! REFF
394    if (iaer.ne.2) then
395      ! Load dataset
396      ierr=NF90_GET_VAR(gcmfid,reffvarid(iaer),reff(:,:,:,:))
397      error_text="Failed to load effective radius"
398      call status_check(ierr,error_text)
399      write(*,*) "Effective radius loaded from "//trim(gcmfile)
400   
401    else ! water ice special case
402      IF (reffwice_val.eq.0) THEN
403        ! Load dataset
404        ierr=NF90_GET_VAR(gcmfid,reffvarid(iaer),reff(:,:,1,:)) ! reffice is actually a GCM
405                                                              ! surface (column-averaged) variable
406        error_text="Failed to load effective radius"
407        call status_check(ierr,error_text)
408        do ialt=2,altlen
409          reff(:,:,ialt,:) = reff(:,:,1,:) ! build up along altitude axis
410        enddo
411        write(*,*) "Effective radius loaded from "//trim(gcmfile)
412      ELSE ! if reffwice_val/=0
413        reff(:,:,:,:) = reffwice_val
414        write(*,*) "Effective radius loaded from the user input"
415      ENDIF
416     
417    endif ! iaer.ne.2
418   
419 
420    ! Length allocation of the output variable
421    ALLOCATE(opa_aer(lonlen,latlen,altlen,timelen))
422   
423 
424!===============================================================================
425! 2. OPTICAL PROPERTIES
426!===============================================================================
427!==========================================================================
428! 2.1 Open the ASCII file
429!==========================================================================
430    IF (wvl_val.lt.5.e-6) THEN
431      ! "VISIBLE" DOMAIN
432     
433      SELECT CASE(iaer)
434      CASE(1) ! DUST
435        ! Dust file
436        optpropfile = "optprop_dustvis_TM_n50.dat"
437      CASE(2) ! WATER ICE
438        ! Water ice file
439        optpropfile = "optprop_icevis_n30.dat"
440      CASE(3) ! STORMDUST
441        ! Dust file
442        optpropfile = "optprop_dustvis_TM_n50.dat"
443      CASE(4) ! TOPDUST
444        ! Dust file
445        optpropfile = "optprop_dustvis_TM_n50.dat"
446! + NEW AER?
447      END SELECT ! iaer
448     
449    ELSE ! wvl_val.ge.5.e-6
450      ! "INFRARED" DOMAIN
451     
452      SELECT CASE(iaer)
453      CASE(1) ! DUST
454        ! Dust file
455        optpropfile = "optprop_dustir_n50.dat"
456      CASE(2) ! WATER ICE
457        ! Water ice file
458        optpropfile = "optprop_iceir_n30.dat"
459      CASE(3) ! STORMDUST
460        ! Dust file
461        optpropfile = "optprop_dustir_n50.dat"
462      CASE(4) ! TOPDUST
463        ! Dust file
464        optpropfile = "optprop_dustir_n50.dat"
465! + NEW AER?
466      END SELECT ! iaer
467     
468    ENDIF ! wvl_val
469
470    INQUIRE(FILE=trim(datadir)//'/'//trim(optpropfile),EXIST=file_ok)
471    if(.not.file_ok) then
472      write(*,*)'Problem opening ',trim(optpropfile)
473      write(*,*)'It should be in: ',trim(datadir)
474      stop
475    endif
476    OPEN(UNIT=file_unit,FILE=trim(datadir)//'/'//trim(optpropfile),FORM='formatted')
477
478!==========================================================================
479! 2.2 Allocate the optical property table
480!==========================================================================
481    jfile = 1
482    endwhile = .false.
483    DO WHILE (.NOT.endwhile)
484      READ(file_unit,*,iostat=read_ok) scanline
485      if (read_ok.ne.0) then
486        write(*,*)'Error reading file ',&
487                  trim(datadir)//'/'//trim(optpropfile)
488        stop
489      endif
490      IF ((scanline(1:1).ne.'#').and.(scanline(1:1).ne.' ')) THEN
491        BACKSPACE(file_unit)
492reading1_seq: SELECT CASE (jfile) ! FIRST READING SEQUENCE
493        CASE(1) reading1_seq ! nwvl ----------------------------
494            read(file_unit,*,iostat=read_ok) nwvl
495            if (read_ok.ne.0) then
496              write(*,*)'reading1_seq: Error while reading line: ',&
497                        trim(scanline)
498              write(*,*)'   of file ',&
499                        trim(datadir)//'/'//trim(optpropfile)
500              stop
501            endif
502            jfile = jfile+1
503        CASE(2) reading1_seq ! nsize ---------------------------
504            read(file_unit,*,iostat=read_ok) nsize
505            if (read_ok.ne.0) then
506              write(*,*)'reading1_seq: Error while reading line: ',&
507                        trim(scanline)
508              write(*,*)'   of file ',&
509                        trim(datadir)//'/'//trim(optpropfile)
510              stop
511            endif
512            endwhile = .true.
513        CASE DEFAULT reading1_seq ! default --------------------
514            write(*,*) 'reading1_seq: ',&
515                       'Error while loading optical properties.'
516            stop
517        END SELECT reading1_seq ! ==============================
518      ENDIF
519    ENDDO ! DO WHILE(.NOT.endwhile)
520
521    ALLOCATE(wvl(nwvl))             ! Wavelength axis
522    ALLOCATE(radiusdyn(nsize))      ! Aerosol radius axis
523    ALLOCATE(Qext(nwvl,nsize))      ! Extinction efficiency coeff
524    ALLOCATE(omeg(nwvl,nsize))      ! Single scattering albedo
525
526!==========================================================================
527! 2.3 Read the data
528!==========================================================================
529    jfile = 1
530    endwhile = .false.
531    DO WHILE (.NOT.endwhile)
532      READ(file_unit,*) scanline
533      IF ((scanline(1:1).ne.'#').and.(scanline(1:1).ne.' ')) THEN
534        BACKSPACE(file_unit)
535reading2_seq: SELECT CASE (jfile) ! SECOND READING SEQUENCE
536        CASE(1) reading2_seq ! wvl -----------------------------
537            read(file_unit,*,iostat=read_ok) wvl
538            if (read_ok.ne.0) then
539              write(*,*)'reading2_seq: Error while reading line: ',&
540                        trim(scanline)
541              write(*,*)'   of file ',&
542                        trim(datadir)//'/'//trim(optpropfile)
543              stop
544            endif
545            jfile = jfile+1
546        CASE(2) reading2_seq ! radiusdyn -----------------------
547            read(file_unit,*,iostat=read_ok) radiusdyn
548            if (read_ok.ne.0) then
549              write(*,*)'reading2_seq: Error while reading line: ',&
550                        trim(scanline)
551              write(*,*)'   of file ',&
552                        trim(datadir)//'/'//trim(optpropfile)
553              stop
554            endif
555            jfile = jfile+1
556        CASE(3) reading2_seq ! Qext ----------------------------
557            isize = 1
558            DO WHILE (isize .le. nsize)
559              READ(file_unit,*,iostat=read_ok) scanline
560              if (read_ok.ne.0) then
561                write(*,*)'reading2_seq: Error while reading line: ',&
562                          trim(scanline)
563                write(*,*)'   of file ',&
564                          trim(datadir)//'/'//trim(optpropfile)
565                stop
566              endif
567              IF ((scanline(1:1).ne.'#').and.(scanline(1:1).ne.' ')) THEN
568                BACKSPACE(file_unit)
569                read(file_unit,*) Qext(:,isize)
570                isize = isize + 1
571              ENDIF
572            ENDDO
573            jfile = jfile+1
574        CASE(4) reading2_seq ! omeg ----------------------------
575            isize = 1
576            DO WHILE (isize .le. nsize)
577              READ(file_unit,*,iostat=read_ok) scanline
578              if (read_ok.ne.0) then
579                write(*,*)'reading2_seq: Error while reading line: ',&
580                          trim(scanline)
581                write(*,*)'   of file ',&
582                          trim(datadir)//'/'//trim(optpropfile)
583                stop
584              endif
585              IF ((scanline(1:1).ne.'#').and.(scanline(1:1).ne.' ')) THEN
586                BACKSPACE(file_unit)
587                read(file_unit,*) omeg(:,isize)
588                isize = isize + 1
589              ENDIF
590            ENDDO
591            endwhile = .true.
592        CASE DEFAULT reading2_seq ! default --------------------
593            write(*,*) 'reading2_seq: ',&
594                       'Error while loading optical properties.'
595            stop
596        END SELECT reading2_seq ! ==============================
597      ENDIF
598    ENDDO
599
600    ! Close the file
601    CLOSE(file_unit)
602   
603    write(*,*)""
604    write(*,*) "Wavelength array loaded from file ",trim(datadir)//'/'//trim(optpropfile)
605    write(*,*) "ranging from ",wvl(1)," to ",wvl(nwvl)," meters"
606    write(*,*) "Effective radius array loaded from file ",trim(datadir)//'/'//trim(optpropfile)
607    write(*,*) "ranging from ",radiusdyn(1)," to ",radiusdyn(nsize)," meters"
608
609! + NEW AER? : one may adapt this part to handle the properties' file of the new aerosol
610
611!==========================================================================
612! 2.4 Get the optpropfile wavelength values encompassing the ref wavelength
613!==========================================================================
614    iwvl=1
615    DO WHILE ((iwvl.le.nwvl).and.(wvl(iwvl).lt.wvl_val))
616      iwvl=iwvl+1
617    ENDDO
618    if ((iwvl.gt.nwvl).or.(wvl(1).gt.wvl_val)) then
619      write(*,*) "ERROR: the reference wavelength is out of the bounds"
620      write(*,*) "of the file ",trim(datadir)//'/'//trim(optpropfile)
621      write(*,*) "(wvl_inf=",wvl(1),"m ; wvl_sup=",wvl(nwvl),"m)"
622      write(*,*) "Ensure you wrote it well (in meters),"
623      write(*,*) "or supply a new optical properties' file (change in aeroptical.F90 directly)"
624      stop
625    endif
626    if (wvl(iwvl).eq.wvl_val) then
627      ! if the ref wavelength is in the optpropfile
628      iwvl1 = iwvl
629      iwvl2 = iwvl
630    else ! wvl(iwvl)>wvl_val and wvl(iwvl-1)<wvl_val
631      iwvl1 = iwvl-1
632      iwvl2 = iwvl
633    endif
634
635
636!===============================================================================
637! 3. OUTPUT FILE VARIABLES
638!===============================================================================
639!==========================================================================
640! 3.1 Creation of the output individual aerosol opacities
641!==========================================================================
642    ! Switch to netcdf define mode
643    ierr=nf90_redef(outfid)
644    error_text="ERROR: Couldn't start netcdf define mode"
645    call status_check(ierr,error_text)
646    write(*,*)""
647    SELECT CASE (iaer)
648    CASE(1) ! DUST
649      ! Insert the definition of the variable
650      tmpvarname="opadust"
651      ierr=NF90_DEF_VAR(outfid,trim(tmpvarname),nf90_float,varshape,tmpvaridout)
652      error_text="ERROR: Couldn't create "//trim(tmpvarname)//" in the outfile"
653      call status_check(ierr,error_text)
654      write(*,*) trim(tmpvarname)//" has been created in the outfile"
655     
656      ! Write the attributes
657      if (opatype.eq."ext") then
658        ierr=nf90_put_att(outfid,tmpvaridout,"long_name","Dust extinction opacity at reference wavelength")
659      else ! opatype.eq."abs"
660        ierr=nf90_put_att(outfid,tmpvaridout,"long_name","Dust absorption opacity at reference wavelength")
661      endif
662     
663    CASE(2) ! WATER ICE
664      ! Insert the definition of the variable
665      tmpvarname="opawice"
666      ierr=NF90_DEF_VAR(outfid,trim(tmpvarname),nf90_float,varshape,tmpvaridout)
667      error_text="ERROR: Couldn't create "//trim(tmpvarname)//" in the outfile"
668      call status_check(ierr,error_text)
669      write(*,*) trim(tmpvarname)//" has been created in the outfile"
670     
671      ! Write the attributes
672      if (opatype.eq."ext") then
673        ierr=nf90_put_att(outfid,tmpvaridout,"long_name","Water ice extinction opacity at reference wavelength")
674      else ! opatype.eq."abs"
675        ierr=nf90_put_att(outfid,tmpvaridout,"long_name","Water ice absorption opacity at reference wavelength")
676      endif
677     
678    CASE(3) ! STORMDUST
679      ! Insert the definition of the variable
680      tmpvarname="opastormdust"
681      ierr=NF90_DEF_VAR(outfid,trim(tmpvarname),nf90_float,varshape,tmpvaridout)
682      error_text="ERROR: Couldn't create "//trim(tmpvarname)//" in the outfile"
683      call status_check(ierr,error_text)
684      write(*,*) trim(tmpvarname)//" has been created in the outfile"
685     
686      ! Write the attributes
687      if (opatype.eq."ext") then
688        ierr=nf90_put_att(outfid,tmpvaridout,"long_name","Stormdust extinction opacity at reference wavelength")
689      else ! opatype.eq."abs"
690        ierr=nf90_put_att(outfid,tmpvaridout,"long_name","Stormdust absorption opacity at reference wavelength")
691      endif
692
693    CASE(4) ! TOPDUST
694      ! Insert the definition of the variable
695      tmpvarname="opatopdust"
696      ierr=NF90_DEF_VAR(outfid,trim(tmpvarname),nf90_float,varshape,tmpvaridout)
697      error_text="ERROR: Couldn't create "//trim(tmpvarname)//" in the outfile"
698      call status_check(ierr,error_text)
699      write(*,*) trim(tmpvarname)//" has been created in the outfile"
700     
701      ! Write the attributes
702      if (opatype.eq."ext") then
703        ierr=nf90_put_att(outfid,tmpvaridout,"long_name","Topdust extinction opacity at reference wavelength")
704      else ! opatype.eq."abs"
705        ierr=nf90_put_att(outfid,tmpvaridout,"long_name","Topdust absorption opacity at reference wavelength")
706      endif
707
708! + NEW AER?   
709    END SELECT ! iaer
710   
711    ierr=nf90_put_att(outfid,tmpvaridout,"units","opacity/km")
712    ierr=nf90_put_att(outfid,tmpvaridout,"refwavelength",wvl_val)
713    ierr=nf90_put_att(outfid,tmpvaridout,"missing_value",missval(iaer))
714    write(*,*)"with missing value = ",missval(iaer)
715   
716    ! End netcdf define mode
717    ierr=nf90_enddef(outfid)
718    error_text="ERROR: Couldn't end netcdf define mode"
719    call status_check(ierr,error_text)
720     
721!==========================================================================
722! 3.2 Computation of the opacities
723!==========================================================================
724    do ilon=1,lonlen
725     do ilat=1,latlen
726      do ialt=1,altlen
727       do it=1,timelen
728         ! Get the optpropfile reff values encompassing the GCM reff
729         isize=1
730         do while((isize.le.nsize).and.(radiusdyn(isize).lt.reff(ilon,ilat,ialt,it)))
731           isize=isize+1
732         enddo
733         if ((isize.gt.nsize).or.(radiusdyn(1).gt.reff(ilon,ilat,ialt,it))) then
734!           write(*,*) "WARNING: the GCM reff (",reff(ilon,ilat,ialt,it),"m) is out of the bounds"
735!           write(*,*) "of the file ",trim(datadir)//'/'//trim(optpropfile)
736!           write(*,*) "(reff_inf=",radiusdyn(1),"m ; reff_sup=",radiusdyn(nsize),"m)"
737!           write(*,*) "A missing value will be written for opacity"
738
739           ! NB: this should especially handle cases when reff=0
740           opa_aer(ilon,ilat,ialt,it)=missval(iaer)
741           
742         else if (mmr(ilon,ilat,ialt,it).eq.missval(iaer)) then
743           ! if there is a missing value in input file
744           opa_aer(ilon,ilat,ialt,it)=missval(iaer)
745           
746         else
747           if (radiusdyn(isize).eq.reff(ilon,ilat,ialt,it)) then
748             ! if the GCM reff is exactly in the optpropfile
749             isize1 = isize
750             isize2 = isize
751           else ! radius(isize)>reff and radiusdyn(isize-1)<reff
752             isize1 = isize-1
753             isize2 = isize
754           endif
755         
756           ! Qext COMPUTATION
757           ! Linear interpolation in effective radius
758           if (isize2.ne.isize1) then
759             coeff = (reff(ilon,ilat,ialt,it)-radiusdyn(isize1))/(radiusdyn(isize2)-radiusdyn(isize1))
760           else
761             coeff = 0.
762           endif
763           reffQext(1) = Qext(iwvl1,isize1)+coeff*(Qext(iwvl1,isize2)-Qext(iwvl1,isize1))
764           reffQext(2) = Qext(iwvl2,isize1)+coeff*(Qext(iwvl2,isize2)-Qext(iwvl2,isize1))
765           ! Linear interpolation in wavelength
766           if (iwvl2.ne.iwvl1) then
767             coeff = (wvl_val-wvl(iwvl1))/(wvl(iwvl2)-wvl(iwvl1))
768           else
769             coeff = 0.
770           endif
771           Qext_val = reffQext(1)+coeff*(reffQext(2)-reffQext(1))
772           
773           ! COMPUTATION OF THE EXTINCTION OPACITY [1/km]
774           opa_aer(ilon,ilat,ialt,it)= 750.*Qext_val*mmr(ilon,ilat,ialt,it)*rho(ilon,ilat,ialt,it)&
775                                               / ( rho_aer(iaer) * reff(ilon,ilat,ialt,it) )
776             
777             
778           if (opatype.eq."abs") then
779             ! omeg COMPUTATION
780             ! Linear interpolation in effective radius
781             if (isize2.ne.isize1) then
782               coeff = (reff(ilon,ilat,ialt,it)-radiusdyn(isize1))/(radiusdyn(isize2)-radiusdyn(isize1))
783             else
784               coeff = 0.
785             endif
786             reffomeg(1) = omeg(iwvl1,isize1)+coeff*(omeg(iwvl1,isize2)-omeg(iwvl1,isize1))
787             reffomeg(2) = omeg(iwvl2,isize1)+coeff*(omeg(iwvl2,isize2)-omeg(iwvl2,isize1))
788             ! Linear interpolation in wavelength
789             if (iwvl2.ne.iwvl1) then
790               coeff = (wvl_val-wvl(iwvl1))/(wvl(iwvl2)-wvl(iwvl1))
791             else
792               coeff = 0.
793             endif
794             omeg_val = reffomeg(1)+coeff*(reffomeg(2)-reffomeg(1))
795
796             ! COMPUTATION OF THE ABSORPTION OPACITY [1/km]
797             opa_aer(ilon,ilat,ialt,it)= (1 - omeg_val) * opa_aer(ilon,ilat,ialt,it)
798           endif ! opatype.eq."abs"
799                 
800         endif
801       enddo ! it
802      enddo ! ialt
803     enddo ! ilat
804    enddo ! ilon
805
806!==========================================================================
807! 3.3 Writing in the output file
808!==========================================================================
809    ierr = NF90_PUT_VAR(outfid, tmpvaridout, opa_aer(:,:,:,:))
810    error_text="Error: could not write "//trim(tmpvarname)//" data in the outfile"
811    call status_check(ierr,error_text)
812
813!==========================================================================
814! 3.4 Prepare next loop
815!==========================================================================
816    DEALLOCATE(mmr)
817    DEALLOCATE(reff)
818    DEALLOCATE(opa_aer)
819    DEALLOCATE(wvl)
820    DEALLOCATE(radiusdyn)
821    DEALLOCATE(Qext)
822    DEALLOCATE(omeg)
823   
824  endif ! if aerok(iaer)
825
826ENDDO ! iaer
827
828!===============================================================================
829! 4. Close the files and end the program
830!===============================================================================
831ierr = nf90_close(gcmfid)
832error_text="Error: could not close file "//trim(gcmfile)
833call status_check(ierr,error_text)
834
835ierr = nf90_close(outfid)
836error_text="Error: could not close file "//trim(outfile)
837call status_check(ierr,error_text)
838
839write(*,*)"";write(*,*)"Program aeroptical completed!"
840
841end program aeroptical
842
843!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
844!! SUBROUTINES
845!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
846
847subroutine status_check(ierr,error_text)
848
849use netcdf
850implicit none
851!================================================================
852!  Arguments:
853!================================================================
854integer,intent(in) :: ierr ! NetCDF status number
855character(len=256),intent(in) :: error_text
856
857if (ierr.ne.nf90_noerr) then
858  write(*,*)trim(error_text)
859  write(*,*)trim(nf90_strerror(ierr))
860  stop
861endif
862
863end subroutine status_check
864
865
866!*******************************************************************************
867
868subroutine inidims(gcmfile,gcmfid,outfile,outfid,lonlen,latlen,altlen,timelen,&
869                   latdimout,londimout,altdimout,timedimout,&
870                   GCM_layers,layerdimout,interlayerdimout)
871
872!==============================================================================
873! Purpose:
874! Read the dimensions of the input file
875! and write them in the output file
876!==============================================================================
877! Remarks:
878! The NetCDF files must be open
879!==============================================================================
880use netcdf
881
882implicit none
883
884!==============================================================================
885! Arguments:
886!==============================================================================
887character(len=128),intent(in) :: gcmfile ! name of the netcdf GCM input file
888integer,intent(in) :: gcmfid ! [netcdf] GCM input file ID number
889character(len=128),intent(in) :: outfile ! name of the netcdf output file
890integer,intent(in) :: outfid ! [netcdf] file ID number
891
892integer,intent(out) ::  lonlen,latlen,altlen,timelen
893! nb of grid points along longitude,latitude,altitude,time
894integer,intent(out) :: latdimout,londimout,altdimout,timedimout
895! latdimout: stores latitude dimension ID number
896! londimout: stores longitude dimension ID number
897! altdimout: stores altitude dimension ID number
898! timedimout: stores time dimension ID number
899integer,intent(out) :: GCM_layers ! number of GCM layers
900integer,intent(out) :: layerdimout ! "GCM_layers" ID
901integer,intent(out) :: interlayerdimout ! "GCM_layers+1" ID
902
903!==============================================================================
904! Local variables:
905!==============================================================================
906integer :: ierr ! [netcdf] subroutine returned error code
907character(len=256) :: error_text ! text to display in case of error
908
909integer :: tmpdimid,tmpvarid ! temporary store a dimension/variable ID number
910character (len=64) :: long_name,units,positive
911! long_name(): [netcdf] long_name attribute
912! units(): [netcdf] units attribute
913! positive(): [netcdf] positive direction attribute (for altitude)
914real, dimension(:), allocatable:: lat,lon,alt,time,ctl
915! lat(): array, stores latitude coordinates
916! lon(): array, stores longitude coordinates
917! alt(): array, stores altitude coordinates
918! time(): array, stores time coordinates
919! ctl(): array, stores controle variable
920integer :: ctllen ! nb of elements in the controle array
921integer :: tmpdimidout,tmpvaridout ! temporary stores a dimension/variable ID number
922
923!==============================================================================
924! LONGITUDE
925!==============================================================================
926! Get the dimension in GCM file
927ierr=nf90_inq_dimid(gcmfid,"longitude",tmpdimid)
928error_text="Error: Dimension <longitude> is missing in file "//trim(gcmfile)
929call status_check(ierr,error_text)
930
931ierr=nf90_inquire_dimension(gcmfid,tmpdimid,len=lonlen)
932allocate(lon(lonlen))
933
934! Create the dimension in output file
935ierr=NF90_DEF_DIM(outfid,"longitude",lonlen,londimout)
936error_text="Error: could not define the longitude dimension in the outfile"
937call status_check(ierr,error_text)
938
939! Get the field in GCM file
940ierr=nf90_inq_varid(gcmfid,"longitude",tmpvarid)
941error_text="Error: Field <longitude> is missing in file "//trim(gcmfile)
942call status_check(ierr,error_text)
943
944ierr=NF90_GET_VAR(gcmfid,tmpvarid,lon)
945error_text="Failed to load longitude"
946call status_check(ierr,error_text)
947
948! Create the field in the output file
949ierr=NF90_DEF_VAR(outfid,"longitude",nf90_float,(/londimout/),tmpvaridout)
950error_text="Error: could not define the longitude variable in the outfile"
951call status_check(ierr,error_text)
952
953! Get the field attributes in the GCM file
954ierr=nf90_get_att(gcmfid,tmpvarid,"long_name",long_name)
955if (ierr.ne.nf90_noerr) then
956! if no attribute "long_name", try "title"
957  ierr=nf90_get_att(gcmfid,tmpvarid,"title",long_name)
958endif
959ierr=nf90_get_att(gcmfid,tmpvarid,"units",units)
960
961! Put the field attributes in the output file
962ierr=nf90_put_att(outfid,tmpvaridout,"long_name",long_name)
963ierr=nf90_put_att(outfid,tmpvaridout,"units",units)
964
965! Write the field values in the output file
966ierr=nf90_enddef(outfid) ! end netcdf define mode
967error_text="Error: could not end the define mode of the outfile"
968call status_check(ierr,error_text)
969
970ierr=NF90_PUT_VAR(outfid,tmpvaridout,lon)
971error_text="Error: could not write the longitude field in the outfile"
972call status_check(ierr,error_text)
973
974!==============================================================================
975! LATITUDE
976!==============================================================================
977! Switch to netcdf define mode
978ierr=nf90_redef(outfid)
979error_text="Error: could not switch to define mode in the outfile"
980call status_check(ierr,error_text)
981
982! Get the dimension in GCM file
983ierr=nf90_inq_dimid(gcmfid,"latitude",tmpdimid)
984error_text="Error: Dimension <latitude> is missing in file "//trim(gcmfile)
985call status_check(ierr,error_text)
986
987ierr=nf90_inquire_dimension(gcmfid,tmpdimid,len=latlen)
988allocate(lat(latlen))
989
990! Create the dimension in output file
991ierr=NF90_DEF_DIM(outfid,"latitude",latlen,latdimout)
992error_text="Error: could not define the latitude dimension in the outfile"
993call status_check(ierr,error_text)
994
995! Get the field in GCM file
996ierr=nf90_inq_varid(gcmfid,"latitude",tmpvarid)
997error_text="Error: Field <latitude> is missing in file "//trim(gcmfile)
998call status_check(ierr,error_text)
999
1000ierr=NF90_GET_VAR(gcmfid,tmpvarid,lat)
1001error_text="Failed to load latitude"
1002call status_check(ierr,error_text)
1003
1004! Create the field in the output file
1005ierr=NF90_DEF_VAR(outfid,"latitude",nf90_float,(/latdimout/),tmpvaridout)
1006error_text="Error: could not define the latitude variable in the outfile"
1007call status_check(ierr,error_text)
1008
1009! Get the field attributes in the GCM file
1010ierr=nf90_get_att(gcmfid,tmpvarid,"long_name",long_name)
1011if (ierr.ne.nf90_noerr) then
1012! if no attribute "long_name", try "title"
1013  ierr=nf90_get_att(gcmfid,tmpvarid,"title",long_name)
1014endif
1015ierr=nf90_get_att(gcmfid,tmpvarid,"units",units)
1016
1017! Put the field attributes in the output file
1018ierr=nf90_put_att(outfid,tmpvaridout,"long_name",long_name)
1019ierr=nf90_put_att(outfid,tmpvaridout,"units",units)
1020
1021! Write the field values in the output file
1022ierr=nf90_enddef(outfid) ! end netcdf define mode
1023error_text="Error: could not end the define mode of the outfile"
1024call status_check(ierr,error_text)
1025
1026ierr=NF90_PUT_VAR(outfid,tmpvaridout,lat)
1027error_text="Error: could not write the latitude field in the outfile"
1028call status_check(ierr,error_text)
1029
1030!==============================================================================
1031! ALTITUDE
1032!==============================================================================
1033! Switch to netcdf define mode
1034ierr=nf90_redef(outfid)
1035error_text="Error: could not switch to define mode in the outfile"
1036call status_check(ierr,error_text)
1037
1038! Get the dimension in GCM file
1039ierr=nf90_inq_dimid(gcmfid,"altitude",tmpdimid)
1040error_text="Error: Dimension <altitude> is missing in file "//trim(gcmfile)
1041call status_check(ierr,error_text)
1042
1043ierr=nf90_inquire_dimension(gcmfid,tmpdimid,len=altlen)
1044allocate(alt(altlen))
1045
1046! Create the dimension in output file
1047ierr=NF90_DEF_DIM(outfid,"altitude",altlen,altdimout)
1048error_text="Error: could not define the altitude dimension in the outfile"
1049call status_check(ierr,error_text)
1050
1051! Get the field in GCM file
1052ierr=nf90_inq_varid(gcmfid,"altitude",tmpvarid)
1053error_text="Error: Field <altitude> is missing in file "//trim(gcmfile)
1054call status_check(ierr,error_text)
1055
1056ierr=NF90_GET_VAR(gcmfid,tmpvarid,alt)
1057error_text="Failed to load altitude"
1058call status_check(ierr,error_text)
1059
1060! Create the field in the output file
1061ierr=NF90_DEF_VAR(outfid,"altitude",nf90_float,(/altdimout/),tmpvaridout)
1062error_text="Error: could not define the altitude variable in the outfile"
1063call status_check(ierr,error_text)
1064
1065! Get the field attributes in the GCM file
1066ierr=nf90_get_att(gcmfid,tmpvarid,"long_name",long_name)
1067if (ierr.ne.nf90_noerr) then
1068! if no attribute "long_name", try "title"
1069  ierr=nf90_get_att(gcmfid,tmpvarid,"title",long_name)
1070endif
1071ierr=nf90_get_att(gcmfid,tmpvarid,"units",units)
1072ierr=nf90_get_att(gcmfid,tmpvarid,"positive",positive)
1073
1074! Put the field attributes in the output file
1075ierr=nf90_put_att(outfid,tmpvaridout,"long_name",long_name)
1076ierr=nf90_put_att(outfid,tmpvaridout,"units",units)
1077ierr=nf90_put_att(outfid,tmpvaridout,"positive",positive)
1078
1079! Write the field values in the output file
1080ierr=nf90_enddef(outfid) ! end netcdf define mode
1081error_text="Error: could not end the define mode of the outfile"
1082call status_check(ierr,error_text)
1083
1084ierr=NF90_PUT_VAR(outfid,tmpvaridout,alt)
1085error_text="Error: could not write the altitude field in the outfile"
1086call status_check(ierr,error_text)
1087
1088!==============================================================================
1089! TIME
1090!==============================================================================
1091! Switch to netcdf define mode
1092ierr=nf90_redef(outfid)
1093error_text="Error: could not switch to define mode in the outfile"
1094call status_check(ierr,error_text)
1095
1096! Get the dimension in GCM file
1097ierr=nf90_inq_dimid(gcmfid,"Time",tmpdimid)
1098error_text="Error: Dimension <Time> is missing in file "//trim(gcmfile)
1099call status_check(ierr,error_text)
1100
1101ierr=nf90_inquire_dimension(gcmfid,tmpdimid,len=timelen)
1102allocate(time(timelen))
1103
1104! Create the dimension in output file
1105ierr=NF90_DEF_DIM(outfid,"Time",timelen,timedimout)
1106error_text="Error: could not define the time dimension in the outfile"
1107call status_check(ierr,error_text)
1108
1109! Get the field in GCM file
1110ierr=nf90_inq_varid(gcmfid,"Time",tmpvarid)
1111error_text="Error: Field <Time> is missing in file "//trim(gcmfile)
1112call status_check(ierr,error_text)
1113
1114ierr=NF90_GET_VAR(gcmfid,tmpvarid,time)
1115error_text="Failed to load Time"
1116call status_check(ierr,error_text)
1117
1118! Create the field in the output file
1119ierr=NF90_DEF_VAR(outfid,"Time",nf90_float,(/timedimout/),tmpvaridout)
1120error_text="Error: could not define the Time variable in the outfile"
1121call status_check(ierr,error_text)
1122
1123! Get the field attributes in the GCM file
1124ierr=nf90_get_att(gcmfid,tmpvarid,"long_name",long_name)
1125if (ierr.ne.nf90_noerr) then
1126! if no attribute "long_name", try "title"
1127  ierr=nf90_get_att(gcmfid,tmpvarid,"title",long_name)
1128endif
1129ierr=nf90_get_att(gcmfid,tmpvarid,"units",units)
1130
1131! Put the field attributes in the output file
1132ierr=nf90_put_att(outfid,tmpvaridout,"long_name",long_name)
1133ierr=nf90_put_att(outfid,tmpvaridout,"units",units)
1134
1135! Write the field values in the output file
1136ierr=nf90_enddef(outfid) ! end netcdf define mode
1137error_text="Error: could not end the define mode of the outfile"
1138call status_check(ierr,error_text)
1139
1140ierr=NF90_PUT_VAR(outfid,tmpvaridout,time)
1141error_text="Error: could not write the Time field in the outfile"
1142call status_check(ierr,error_text)
1143
1144!==============================================================================
1145! CONTROLE
1146!==============================================================================
1147! Switch to netcdf define mode
1148ierr=nf90_redef(outfid)
1149error_text="Error: could not switch to define mode in the outfile"
1150call status_check(ierr,error_text)
1151
1152! Get the dimension in GCM file
1153ierr=nf90_inq_dimid(gcmfid,"index",tmpdimid)
1154error_text="Dimension <index> is missing in file "//trim(gcmfile)&
1155           //". We'll skip that one."
1156if (ierr.ne.nf90_noerr) then
1157  write(*,*)trim(error_text)
1158  ierr=nf90_enddef(outfid) ! end netcdf define mode
1159  error_text="Error: could not end the define mode of the outfile"
1160  call status_check(ierr,error_text)
1161else
1162  ierr=nf90_inquire_dimension(gcmfid,tmpdimid,len=ctllen)
1163  allocate(ctl(ctllen))
1164
1165  ! Create the dimension in output file
1166  ierr=NF90_DEF_DIM(outfid,"index",ctllen,tmpdimidout)
1167  error_text="Error: could not define the index dimension in the outfile"
1168  call status_check(ierr,error_text)
1169
1170  ! Get the field in GCM file
1171  ierr=nf90_inq_varid(gcmfid,"controle",tmpvarid)
1172  error_text="Error: Field <controle> is missing in file "//trim(gcmfile)
1173  call status_check(ierr,error_text)
1174
1175  ierr=NF90_GET_VAR(gcmfid,tmpvarid,ctl)
1176  error_text="Failed to load ctl"
1177  call status_check(ierr,error_text)
1178
1179  ! Create the field in the output file
1180  ierr=NF90_DEF_VAR(outfid,"controle",nf90_float,(/tmpdimidout/),tmpvaridout)
1181  error_text="Error: could not define the controle variable in the outfile"
1182  call status_check(ierr,error_text)
1183
1184  ! Get the field attributes in the GCM file
1185  ierr=nf90_get_att(gcmfid,tmpvarid,"long_name",long_name)
1186  if (ierr.ne.nf90_noerr) then
1187  ! if no attribute "long_name", try "title"
1188    ierr=nf90_get_att(gcmfid,tmpvarid,"title",long_name)
1189  endif
1190
1191  ! Put the field attributes in the output file
1192  ierr=nf90_put_att(outfid,tmpvaridout,"long_name",long_name)
1193
1194  ! Write the field values in the output file
1195  ierr=nf90_enddef(outfid) ! end netcdf define mode
1196  error_text="Error: could not end the define mode of the outfile"
1197  call status_check(ierr,error_text)
1198
1199  ierr=NF90_PUT_VAR(outfid,tmpvaridout,ctl)
1200  error_text="Error: could not write the controle field in the outfile"
1201  call status_check(ierr,error_text)
1202endif
1203
1204
1205!==============================================================================
1206! Load size of aps() or sigma() (in case it is not altlen)
1207!==============================================================================
1208! Switch to netcdf define mode
1209ierr=nf90_redef(outfid)
1210error_text="Error: could not switch to define mode in the outfile"
1211call status_check(ierr,error_text)
1212
1213! Default is that GCM_layers=altlen
1214! but for outputs of zrecast, it may be a different value
1215ierr=nf90_inq_dimid(gcmfid,"GCM_layers",tmpdimid)
1216if (ierr.ne.nf90_noerr) then
1217  ! didn't find a GCM_layers dimension; therefore we have:
1218  GCM_layers=altlen 
1219else
1220  ! load value of GCM_layers
1221  ierr=nf90_inquire_dimension(gcmfid,tmpdimid,len=GCM_layers)
1222endif
1223
1224! Create the dimensions in output file
1225ierr = NF90_DEF_DIM(outfid,"GCM_layers",GCM_layers,layerdimout)
1226error_text="Error: could not define the GCM_layers dimension in the outfile"
1227call status_check(ierr,error_text)
1228ierr = NF90_DEF_DIM(outfid,"GCM_interlayers",GCM_layers+1,interlayerdimout)
1229error_text="Error: could not define the GCM_interlayers dimension in the outfile"
1230call status_check(ierr,error_text)
1231
1232! End netcdf define mode
1233ierr=nf90_enddef(outfid)
1234error_text="Error: could not end the define mode of the outfile"
1235call status_check(ierr,error_text)
1236
1237end subroutine inidims
1238
1239
1240!*******************************************************************************
1241
1242subroutine init2(gcmfid,lonlen,latlen,altlen,GCM_layers, &
1243                 outfid,londimout,latdimout,altdimout, &
1244                 layerdimout,interlayerdimout)
1245!==============================================================================
1246! Purpose:
1247! Copy ap() , bp(), aps(), bps(), aire() and phisinit()
1248! from input file to outpout file
1249!==============================================================================
1250! Remarks:
1251! The NetCDF files must be open
1252!==============================================================================
1253
1254use netcdf
1255
1256implicit none
1257
1258!==============================================================================
1259! Arguments:
1260!==============================================================================
1261integer, intent(in) :: gcmfid  ! NetCDF output file ID
1262integer, intent(in) :: lonlen ! # of grid points along longitude
1263integer, intent(in) :: latlen ! # of grid points along latitude
1264integer, intent(in) :: altlen ! # of grid points along altitude
1265integer, intent(in) :: GCM_layers ! # of GCM atmospheric layers
1266integer, intent(in) :: outfid ! NetCDF output file ID
1267integer, intent(in) :: londimout ! longitude dimension ID
1268integer, intent(in) :: latdimout ! latitude dimension ID
1269integer, intent(in) :: altdimout ! altitude dimension ID
1270integer, intent(in) :: layerdimout ! GCM_layers dimension ID
1271integer, intent(in) :: interlayerdimout ! GCM_layers+1 dimension ID
1272!==============================================================================
1273! Local variables:
1274!==============================================================================
1275real,dimension(:),allocatable :: aps,bps ! hybrid vertical coordinates
1276real,dimension(:),allocatable :: ap,bp ! hybrid vertical coordinates
1277real,dimension(:),allocatable :: sigma ! sigma levels
1278real,dimension(:,:),allocatable :: aire ! mesh areas
1279real,dimension(:,:),allocatable :: phisinit ! Ground geopotential
1280integer :: ierr
1281integer :: tmpvarid ! temporary variable ID
1282logical :: area ! is "aire" available ?
1283logical :: phis ! is "phisinit" available ?
1284logical :: hybrid ! are "aps" and "bps" available ?
1285logical :: apbp ! are "ap" and "bp" available ?
1286
1287!==============================================================================
1288! 1. Read data from input file
1289!==============================================================================
1290
1291! hybrid coordinate aps
1292!write(*,*) "aps: altlen=",altlen," GCM_layers=",GCM_layers
1293allocate(aps(GCM_layers),stat=ierr)
1294if (ierr.ne.0) then
1295  write(*,*) "init2: failed to allocate aps!"
1296  stop
1297endif
1298ierr=nf90_inq_varid(gcmfid,"aps",tmpvarid)
1299if (ierr.ne.nf90_noerr) then
1300  write(*,*) "Ooops. Failed to get aps ID. OK, will look for sigma coord."
1301  hybrid=.false.
1302else
1303  ierr=NF90_GET_VAR(gcmfid,tmpvarid,aps)
1304  hybrid=.true.
1305  if (ierr.ne.nf90_noerr) then
1306    write(*,*) "init2 Error: Failed reading aps"
1307    stop
1308  endif
1309
1310  ! hybrid coordinate bps
1311!  write(*,*) "bps: altlen=",altlen," GCM_layers=",GCM_layers
1312  allocate(bps(GCM_layers),stat=ierr)
1313  if (ierr.ne.0) then
1314    write(*,*) "init2: failed to allocate bps!"
1315    stop
1316  endif
1317  ierr=nf90_inq_varid(gcmfid,"bps",tmpvarid)
1318  if (ierr.ne.nf90_noerr) then
1319    write(*,*) "init2 Error: Failed to get bps ID."
1320    stop
1321  endif
1322  ierr=NF90_GET_VAR(gcmfid,tmpvarid,bps)
1323  if (ierr.ne.nf90_noerr) then
1324    write(*,*) "init2 Error: Failed reading bps"
1325    stop
1326  endif
1327endif
1328
1329! hybrid coordinate ap
1330allocate(ap(GCM_layers+1),stat=ierr)
1331if (ierr.ne.0) then
1332  write(*,*) "init2: failed to allocate ap!"
1333  stop
1334else
1335  ierr=nf90_inq_varid(gcmfid,"ap",tmpvarid)
1336  if (ierr.ne.nf90_noerr) then
1337    write(*,*) "Ooops. Failed to get ap ID. OK."
1338    apbp=.false.
1339  else
1340    ierr=NF90_GET_VAR(gcmfid,tmpvarid,ap)
1341    apbp=.true.
1342    if (ierr.ne.nf90_noerr) then
1343      write(*,*) "Error: Failed reading ap"
1344      stop
1345    endif
1346  endif
1347endif
1348
1349! hybrid coordinate bp
1350allocate(bp(GCM_layers+1),stat=ierr)
1351if (ierr.ne.0) then
1352  write(*,*) "init2: failed to allocate bp!"
1353  stop
1354else
1355  ierr=nf90_inq_varid(gcmfid,"bp",tmpvarid)
1356  if (ierr.ne.nf90_noerr) then
1357    write(*,*) "Ooops. Failed to get bp ID. OK."
1358    apbp=.false.
1359  else
1360    ierr=NF90_GET_VAR(gcmfid,tmpvarid,bp)
1361    apbp=.true.
1362    if (ierr.ne.nf90_noerr) then
1363      write(*,*) "Error: Failed reading bp"
1364      stop
1365    endif
1366  endif
1367endif
1368
1369! sigma levels (if any)
1370if (.not.hybrid) then
1371  allocate(sigma(GCM_layers),stat=ierr)
1372  if (ierr.ne.0) then
1373    write(*,*) "init2: failed to allocate sigma"
1374    stop
1375  endif
1376  ierr=nf90_inq_varid(gcmfid,"sigma",tmpvarid)
1377  ierr=NF90_GET_VAR(gcmfid,tmpvarid,sigma)
1378  if (ierr.ne.nf90_noerr) then
1379    write(*,*) "init2 Error: Failed reading sigma"
1380    stop
1381  endif
1382endif ! of if (.not.hybrid)
1383
1384! mesh area
1385allocate(aire(lonlen,latlen),stat=ierr)
1386if (ierr.ne.0) then
1387  write(*,*) "init2: failed to allocate aire!"
1388  stop
1389endif
1390ierr=nf90_inq_varid(gcmfid,"aire",tmpvarid)
1391if (ierr.ne.nf90_noerr) then
1392  write(*,*)"init2 warning: Failed to get aire ID."
1393  area = .false.
1394else
1395  ierr=NF90_GET_VAR(gcmfid,tmpvarid,aire)
1396  if (ierr.ne.nf90_noerr) then
1397    write(*,*) "init2 Error: Failed reading aire"
1398    stop
1399  endif
1400  area = .true.
1401endif
1402
1403! ground geopotential phisinit
1404allocate(phisinit(lonlen,latlen),stat=ierr)
1405if (ierr.ne.0) then
1406  write(*,*) "init2: failed to allocate phisinit!"
1407  stop
1408endif
1409ierr=nf90_inq_varid(gcmfid,"phisinit",tmpvarid)
1410if (ierr.ne.nf90_noerr) then
1411  write(*,*)"init2 warning: Failed to get phisinit ID."
1412  phis = .false.
1413else
1414  ierr=NF90_GET_VAR(gcmfid,tmpvarid,phisinit)
1415  if (ierr.ne.nf90_noerr) then
1416    write(*,*) "init2 Error: Failed reading phisinit"
1417    stop
1418  endif
1419  phis = .true.
1420endif
1421
1422!==============================================================================
1423! 2. Write
1424!==============================================================================
1425
1426!==============================================================================
1427! 2.1 Hybrid coordinates ap() , bp(), aps() and bps()
1428!==============================================================================
1429if(hybrid) then
1430! define aps
1431  ! Switch to netcdf define mode
1432  ierr=nf90_redef(outfid)
1433  ! Insert the definition of the variable
1434  ierr=NF90_DEF_VAR(outfid,"aps",nf90_float,(/layerdimout/),tmpvarid)
1435  if (ierr.ne.nf90_noerr) then
1436     write(*,*) "init2 Error: Failed to define the variable aps"
1437     stop
1438  endif
1439  ! Write the attributes
1440  ierr=nf90_put_att(outfid,tmpvarid,"long_name","hybrid pressure at midlayers")
1441  ierr=nf90_put_att(outfid,tmpvarid,"units"," ")
1442  ! End netcdf define mode
1443  ierr=nf90_enddef(outfid)
1444
1445! write aps
1446  ierr=NF90_PUT_VAR(outfid,tmpvarid,aps)
1447  if (ierr.ne.nf90_noerr) then
1448    write(*,*) "init2 Error: Failed to write aps"
1449    stop
1450  endif
1451
1452! define bps
1453  ! Switch to netcdf define mode
1454  ierr=nf90_redef(outfid)
1455  ! Insert the definition of the variable
1456  ierr=NF90_DEF_VAR(outfid,"bps",nf90_float,(/layerdimout/),tmpvarid)
1457  if (ierr.ne.nf90_noerr) then
1458     write(*,*) "init2 Error: Failed to define the variable bps"
1459     stop
1460  endif
1461  ! Write the attributes
1462  ierr=nf90_put_att(outfid,tmpvarid,"long_name","hybrid sigma at midlayers")
1463  ierr=nf90_put_att(outfid,tmpvarid,"units"," ")
1464  ! End netcdf define mode
1465  ierr=nf90_enddef(outfid)
1466   
1467! write bps
1468  ierr=NF90_PUT_VAR(outfid,tmpvarid,bps)
1469  if (ierr.ne.nf90_noerr) then
1470    write(*,*) "init2 Error: Failed to write bps"
1471    stop
1472  endif
1473
1474  if (apbp) then
1475! define ap
1476
1477    ! Switch to netcdf define mode
1478    ierr=nf90_redef(outfid)
1479    ! Insert the definition of the variable
1480    ierr=NF90_DEF_VAR(outfid,"ap",nf90_float,(/interlayerdimout/),tmpvarid)
1481    if (ierr.ne.nf90_noerr) then
1482       write(*,*) "init2 Error: Failed to define the variable ap"
1483       stop
1484    endif
1485    ! Write the attributes
1486    ierr=nf90_put_att(outfid,tmpvarid,"long_name","hybrid sigma at interlayers")
1487    ierr=nf90_put_att(outfid,tmpvarid,"units"," ")
1488    ! End netcdf define mode
1489    ierr=nf90_enddef(outfid)
1490
1491! write ap
1492    ierr=NF90_PUT_VAR(outfid,tmpvarid,ap)
1493    if (ierr.ne.nf90_noerr) then
1494      write(*,*) "Error: Failed to write ap"
1495      stop
1496    endif
1497
1498! define bp
1499
1500    ! Switch to netcdf define mode
1501    ierr=nf90_redef(outfid)
1502    ! Insert the definition of the variable
1503    ierr=NF90_DEF_VAR(outfid,"bp",nf90_float,(/interlayerdimout/),tmpvarid)
1504    if (ierr.ne.nf90_noerr) then
1505       write(*,*) "init2 Error: Failed to define the variable bp"
1506       stop
1507    endif
1508    ! Write the attributes
1509    ierr=nf90_put_att(outfid,tmpvarid,"long_name","hybrid sigma at interlayers")
1510    ierr=nf90_put_att(outfid,tmpvarid,"units"," ")
1511    ! End netcdf define mode
1512    ierr=nf90_enddef(outfid)
1513 
1514! write bp
1515    ierr=NF90_PUT_VAR(outfid,tmpvarid,bp)
1516    if (ierr.ne.nf90_noerr) then
1517      write(*,*) "Error: Failed to write bp"
1518      stop
1519    endif
1520  endif ! of if (apbp)
1521
1522else
1523
1524  ! Switch to netcdf define mode
1525  ierr=nf90_redef(outfid)
1526  ! Insert the definition of the variable
1527  ierr=NF90_DEF_VAR(outfid,"sigma",nf90_float,(/layerdimout/),tmpvarid)
1528  if (ierr.ne.nf90_noerr) then
1529     write(*,*) "init2 Error: Failed to define the variable sigma"
1530     stop
1531  endif
1532  ! Write the attributes
1533  ierr=nf90_put_att(outfid,tmpvarid,"long_name","sigma at midlayers")
1534  ierr=nf90_put_att(outfid,tmpvarid,"units"," ")
1535  ! End netcdf define mode
1536  ierr=nf90_enddef(outfid)
1537 
1538! write sigma
1539  ierr=NF90_PUT_VAR(outfid,tmpvarid,sigma)
1540  if (ierr.ne.nf90_noerr) then
1541    write(*,*) "init2 Error: Failed to write sigma"
1542    stop
1543  endif
1544endif ! of if (hybrid)
1545
1546!==============================================================================
1547! 2.2 aire() and phisinit()
1548!==============================================================================
1549
1550if (area) then
1551
1552  ! Switch to netcdf define mode
1553  ierr=nf90_redef(outfid)
1554  ! Insert the definition of the variable
1555  ierr=NF90_DEF_VAR(outfid,"aire",nf90_float,(/londimout,latdimout/),tmpvarid)
1556  if (ierr.ne.nf90_noerr) then
1557     write(*,*) "init2 Error: Failed to define the variable aire"
1558     stop
1559  endif
1560  ! Write the attributes
1561  ierr=nf90_put_att(outfid,tmpvarid,"long_name","Mesh area")
1562  ierr=nf90_put_att(outfid,tmpvarid,"units","m2")
1563  ! End netcdf define mode
1564  ierr=nf90_enddef(outfid)
1565 
1566  ! write aire
1567  ierr=NF90_PUT_VAR(outfid,tmpvarid,aire)
1568  if (ierr.ne.nf90_noerr) then
1569    write(*,*) "init2 Error: Failed to write aire"
1570    stop
1571  endif
1572endif ! of if (area)
1573
1574if (phis) then
1575
1576  ! Switch to netcdf define mode
1577  ierr=nf90_redef(outfid)
1578  ! Insert the definition of the variable
1579  ierr=NF90_DEF_VAR(outfid,"phisinit",nf90_float,(/londimout,latdimout/),tmpvarid)
1580  if (ierr.ne.nf90_noerr) then
1581    write(*,*) "init2 Error: Failed to define the variable phisinit"
1582    stop
1583  endif
1584  ! Write the attributes
1585  ierr=nf90_put_att(outfid,tmpvarid,"long_name","Ground level geopotential")
1586  ierr=nf90_put_att(outfid,tmpvarid,"units"," ")
1587  ! End netcdf define mode
1588  ierr=nf90_enddef(outfid)
1589
1590  ! write phisinit
1591  ierr=NF90_PUT_VAR(outfid,tmpvarid,phisinit)
1592  if (ierr.ne.nf90_noerr) then
1593    write(*,*) "init2 Error: Failed to write phisinit"
1594    stop
1595  endif
1596
1597endif ! of if (phis)
1598
1599
1600! Cleanup
1601if (allocated(aps)) deallocate(aps)
1602if (allocated(bps)) deallocate(bps)
1603if (allocated(ap)) deallocate(ap)
1604if (allocated(bp)) deallocate(bp)
1605if (allocated(sigma)) deallocate(sigma)
1606if (allocated(phisinit)) deallocate(phisinit)
1607if (allocated(aire)) deallocate(aire)
1608
1609end subroutine init2
Note: See TracBrowser for help on using the repository browser.