SUBROUTINE suaer IMPLICIT NONE !================================================================== ! Purpose. ! -------- ! initialize yomaer, the common that contains the ! radiative characteristics of the aerosols ! ! AUTHORS. ! -------- ! Richard Fournier (1996) Francois Forget (1996) ! Frederic Hourdin ! Jean-jacques morcrette *ECMWF* ! MODIF Francois Forget (2000) ! MODIF Franck Montmessin (add water ice) ! MODIF J.-B. Madeleine 2008W27 ! - Optical properties read in ASCII files ! - Add varying radius of the particules ! ! Summary. ! -------- ! ! Read the optical properties -> Mean -> Variable assignment ! (ASCII files) (see yomaer.h) ! wvl(nwvl) longsun ! ep(nwvl) epav QVISsQREF(nsun) ! omeg(nwvl) omegav omegavis(nsun) ! gfactor(nwvl) gav gvis(nsun) ! !================================================================== ! Includes: #include "callkeys.h" #include "datafile.h" #include "dimensions.h" #include "dimphys.h" #include "dimradmars.h" #include "yomaer.h" #include "aerkind.h" ! Optical properties (read in external ASCII files) INTEGER :: nwvl ! Number of wavelengths in ! the domain (VIS or IR) REAL, DIMENSION(:),& ALLOCATABLE, SAVE :: wvl ! Wavelength axis REAL, DIMENSION(:),& ALLOCATABLE, SAVE :: radiusdyn ! Particle size axis REAL, DIMENSION(:,:),& ALLOCATABLE, SAVE :: ep,& ! Extinction coefficient Qext omeg,& ! Single Scattering Albedo gfactor ! Assymetry Factor ! Local variables: INTEGER :: iaer ! Aerosol index INTEGER :: idomain ! Domain index (1=VIS,2=IR) INTEGER :: iir ! IR channel index ! iir=1: 15um CO2 bands ! iir=2 : CO2 band wings ! iir=3 : 9 um band ! iir=4 : Far IR INTEGER :: isun ! Solar band index INTEGER :: isize ! Particle size index INTEGER :: jfile ! ASCII file scan index INTEGER :: file_unit = 60 LOGICAL :: file_ok, endwhile CHARACTER(LEN=132) :: scanline ! ASCII file scanning line INTEGER :: read_ok ! I/O of "aerave" (subroutine averaging spectrally ! sing.scat.parameters) REAL tsun ! Sun brightness temperature (for SW) REAL tsol ! Surface reference brightness temp (LW) REAL longref ! reference wavelengths REAL longsun(nsun+1) ! solar band boundaries REAL longir(nir+1) ! IR band boundaries REAL epref ! reference extinction ep ! at wavelength "longref" REAL epav(nir) ! average ep ! (= /Qext(longref) if epref=1) REAL omegav(nir) ! Average sing.scat.albedo REAL gav(nir) ! Average assymetry parameter ! Local saved variables: CHARACTER(LEN=30), DIMENSION(naerkind,2), SAVE :: file_id !================================================================== !---- Please indicate the names of the optical property files below ! Please also choose the reference wavelengths of each aerosol DO iaer = 1, naerkind ! Loop on aerosol kind aerkind: SELECT CASE (name_iaer(iaer)) !================================================================== CASE("dust_conrath") aerkind ! Typical dust profile !================================================================== ! Visible domain: file_id(iaer,1) = 'optprop_dustvis_TM.dat' !M.Wolff TM ! file_id(iaer,1) = 'optprop_dustvis_clancy.dat' !Clancy-Lee ! file_id(iaer,1) = 'optprop_dustvis_ockert.dat' !Ockert-Bell ! Infrared domain: file_id(iaer,2) = 'optprop_dustir_TM.dat' !M.Wolff TM ! Toon-Forget + solsir=2 using Clancy-Lee ! file_id(iaer,2) = 'optprop_dustir_clancy.dat' ! Toon-Forget + solsir=2 using Ockert-Bell ! file_id(iaer,2) = 'optprop_dustir_ockert.dat' ! Reference wavelength in the visible: longrefvis(iaer)=0.67E-6 ! For dust: change readtesassim accordingly; ! Reference wavelength in the infrared: longrefir(iaer)=9.3E-6 !================================================================== CASE("dust_doubleq") aerkind! Two-moment scheme for dust !================================================================== ! Visible domain: file_id(iaer,1) = 'optprop_dustvis_TM_n50.dat' !T-Matrix ! file_id(iaer,1) = 'optprop_dustvis_n50.dat' !Mie ! Infrared domain: file_id(iaer,2) = 'optprop_dustir_n50.dat' !Mie ! Reference wavelength in the visible: longrefvis(iaer)=0.67E-6 ! If not equal to 0.67e-6 -> change readtesassim accordingly; ! Reference wavelength in the infrared: longrefir(iaer)=9.3E-6 !================================================================== CASE("dust_submicron") aerkind ! Small dust population !================================================================== ! Visible domain: file_id(iaer,1) = 'optprop_dustvis_01um_TM.dat' !M.Wolff ! Infrared domain: file_id(iaer,2) = 'optprop_dustir_01um_TM.dat' !M.Wolff ! Reference wavelength in the visible: longrefvis(iaer)=0.67E-6 ! If not equal to 0.67e-6 -> change readtesassim accordingly; ! Reference wavelength in the infrared: longrefir(iaer)=9.3E-6 !================================================================== CASE("h2o_ice") aerkind ! Water ice crystals !================================================================== ! Visible domain: file_id(iaer,1) = 'optprop_icevis_n30.dat' !Warren ! file_id(iaer,1) = 'optprop_icevis.dat' !Warren ! Infrared domain: file_id(iaer,2) = 'optprop_iceir_n30.dat' !Warren ! file_id(iaer,2) = 'optprop_iceir.dat' !Warren ! Reference wavelength in the visible: longrefvis(iaer)=0.67E-6 ! 1.5um OMEGA/MEx ! Reference wavelength in the infrared: longrefir(iaer)=12.1E-6 ! 825cm-1 TES/MGS !================================================================== END SELECT aerkind !================================================================== WRITE(*,*) "Scatterer: ",trim(name_iaer(iaer)) WRITE(*,*) " corresponding files: " WRITE(*,*) "VIS: ",trim(file_id(iaer,1)) WRITE(*,*) "IR : ",trim(file_id(iaer,2)) !================================================================== ENDDO ! iaer (loop on aerosol kind) ! Initializations: radiustab(1:naerkind,1:2,1:nsizemax)=0 gvis(1:nsun,1:naerkind,1:nsizemax)=0 omegavis(1:nsun,1:naerkind,1:nsizemax)=0 QVISsQREF(1:nsun,1:naerkind,1:nsizemax)=0 gIR(1:nir,1:naerkind,1:nsizemax)=0 omegaIR(1:nir,1:naerkind,1:nsizemax)=0 QIRsQREF(1:nir,1:naerkind,1:nsizemax)=0 QREFvis(1:naerkind,1:nsizemax)=0 QREFir(1:naerkind,1:nsizemax)=0 omegaREFvis(1:naerkind,1:nsizemax)=0 omegaREFir(1:naerkind,1:nsizemax)=0 DO iaer = 1, naerkind ! Loop on aerosol kind DO idomain = 1, 2 ! Loop on radiation domain (VIS or IR) !================================================================== ! 1. READ OPTICAL PROPERTIES !================================================================== ! 1.1 Open the ASCII file INQUIRE(FILE=datafile(1:LEN_TRIM(datafile))//& '/'//file_id(iaer,idomain)(1:LEN_TRIM(file_id(iaer,idomain))),& EXIST=file_ok) IF(.NOT.file_ok) THEN write(*,*)'Problem opening ',& file_id(iaer,idomain)(1:LEN_TRIM(file_id(iaer,idomain))) write(*,*)'It should be in: ',& datafile(1:LEN_TRIM(datafile)) write(*,*)'1) You can change this directory address in ' write(*,*)' file phymars/datafile.h' write(*,*)'2) If ',& file_id(iaer,idomain)(1:LEN_TRIM(file_id(iaer,idomain))),& ' is a LMD reference datafile, it' write(*,*)' can be obtained online on:' write(*,*)' http://www.lmd.jussieu.fr/',& '~forget/datagcm/datafile' write(*,*)'3) If the name of the file is wrong, you can' write(*,*)' change it in file phymars/suaer.F90. Just' write(*,*)' modify the variable called file_id.' CALL ABORT ENDIF OPEN(UNIT=file_unit,& FILE=datafile(1:LEN_TRIM(datafile))//& '/'//file_id(iaer,idomain)(1:LEN_TRIM(file_id(iaer,idomain))),& FORM='formatted') ! 1.2 Allocate the optical property table jfile = 1 endwhile = .false. DO WHILE (.NOT.endwhile) READ(file_unit,*,iostat=read_ok) scanline if (read_ok.ne.0) then write(*,*)' readoptprop: Error reading file',& trim(datafile(1:LEN_TRIM(datafile))//& '/'//file_id(iaer,idomain)(1:LEN_TRIM(file_id(iaer,idomain)))) call abort endif IF ((scanline(1:1) .ne. '#').and.& (scanline(1:1) .ne. ' ')) THEN BACKSPACE(file_unit) reading1_seq: SELECT CASE (jfile) ! ==================== CASE(1) reading1_seq ! nwvl ---------------------------- read(file_unit,*,iostat=read_ok) nwvl if (read_ok.ne.0) then write(*,*)' readoptprop: Error while reading line:',& trim(scanline) write(*,*)' of file',& trim(datafile(1:LEN_TRIM(datafile))//& '/'//file_id(iaer,idomain)(1:LEN_TRIM(file_id(iaer,idomain)))) call abort endif jfile = jfile+1 CASE(2) reading1_seq ! nsize --------------------------- read(file_unit,*,iostat=read_ok) nsize(iaer,idomain) if (read_ok.ne.0) then write(*,*)' readoptprop: Error while reading line:',& trim(scanline) write(*,*)' of file',& trim(datafile(1:LEN_TRIM(datafile))//& '/'//file_id(iaer,idomain)(1:LEN_TRIM(file_id(iaer,idomain)))) call abort endif endwhile = .true. CASE DEFAULT reading1_seq ! ---------------------------- WRITE(*,*) 'readoptprop: ',& 'Error while loading optical properties.' CALL ABORT END SELECT reading1_seq ! ============================== ENDIF ENDDO ALLOCATE(wvl(nwvl)) ! wvl ALLOCATE(radiusdyn(nsize(iaer,idomain))) ! radiusdyn ALLOCATE(ep(nwvl,nsize(iaer,idomain))) ! ep ALLOCATE(omeg(nwvl,nsize(iaer,idomain))) ! omeg ALLOCATE(gfactor(nwvl,nsize(iaer,idomain))) ! g ! 1.3 Read the data jfile = 1 endwhile = .false. DO WHILE (.NOT.endwhile) READ(file_unit,*) scanline IF ((scanline(1:1) .ne. '#').and.& (scanline(1:1) .ne. ' ')) THEN BACKSPACE(file_unit) reading2_seq: SELECT CASE (jfile) ! ==================== CASE(1) reading2_seq ! wvl ----------------------------- read(file_unit,*) wvl jfile = jfile+1 CASE(2) reading2_seq ! radiusdyn ----------------------- read(file_unit,*) radiusdyn jfile = jfile+1 CASE(3) reading2_seq ! ep ------------------------------ isize = 1 DO WHILE (isize .le. nsize(iaer,idomain)) READ(file_unit,*) scanline IF ((scanline(1:1) .ne. '#').and.& (scanline(1:1) .ne. ' ')) THEN BACKSPACE(file_unit) read(file_unit,*) ep(:,isize) isize = isize + 1 ENDIF ENDDO jfile = jfile+1 CASE(4) reading2_seq ! omeg ---------------------------- isize = 1 DO WHILE (isize .le. nsize(iaer,idomain)) READ(file_unit,*) scanline IF ((scanline(1:1) .ne. '#').and.& (scanline(1:1) .ne. ' ')) THEN BACKSPACE(file_unit) read(file_unit,*) omeg(:,isize) isize = isize + 1 ENDIF ENDDO jfile = jfile+1 CASE(5) reading2_seq ! gfactor ------------------------- isize = 1 DO WHILE (isize .le. nsize(iaer,idomain)) READ(file_unit,*) scanline IF ((scanline(1:1) .ne. '#').and.& (scanline(1:1) .ne. ' ')) THEN BACKSPACE(file_unit) read(file_unit,*) gfactor(:,isize) isize = isize + 1 ENDIF ENDDO endwhile = .true. CASE DEFAULT reading2_seq ! ---------------------------- WRITE(*,*) 'suaer.F90: ',& 'Error while loading optical properties.' CALL ABORT END SELECT reading2_seq ! ============================== ENDIF ENDDO ! 1.4 Close the file CLOSE(file_unit) !================================================================== ! 2. AVERAGED PROPERTIES AND VARIABLE ASSIGNMENTS !================================================================== domain: SELECT CASE (idomain) !================================================================== CASE(1) domain ! VISIBLE DOMAIN (idomain=1) !================================================================== ! 2.1 Parameters tsun=6000.E+0 longsun(1)=long1vis longsun(2)=long2vis longsun(3)=long3vis longref=longrefvis(iaer) epref=1.E+0 DO isize=1,nsize(iaer,idomain) ! test that there is enough room to store the data if (isize.gt.nsizemax) then write(*,*) "suaer: Error ! nsizemax is too small!" write(*,*) " nsizemax=",nsizemax write(*,*) " you must increase the value of nsizemax" write(*,*) " in dimradmars.h !" stop endif ! ------------------------------------------------ ! 2.2 Save the particle sizes radiustab(iaer,idomain,isize)=radiusdyn(isize) ! 2.3 Averaged optical properties (GCM channels) ! Notice: Aerave also computes the extinction coefficient and ! single scattering albedo at reference wavelength ! (called QREFvis and OMEGAREFvis, same in the IR, ! and not epref, which is a different parameter); ! Reference wavelengths are defined for each aerosol in ! dimradmars.h. CALL aerave ( nwvl,& wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),& longref,epref,tsun,& nsun,longsun, epav,omegav,gav,& QREFvis(iaer,isize),omegaREFvis(iaer,isize) ) ! 2.4 Variable assignements (declared by yomaer.h) DO isun=1,nsun QVISsQREF(isun,iaer,isize)=epav(isun) gvis(isun,iaer,isize)=gav(isun) omegavis(isun,iaer,isize)=omegav(isun) END DO ! 2.5 Output display WRITE(*,*) 'Les donnees spectrales :' WRITE(*,*) 'Solaire (SW) ---->' WRITE(*,*) 'Aerosol number: ', iaer WRITE(*,*) 'Rayon aerosol: ', radiustab(iaer,idomain,isize) WRITE(*,*) '/Qext(longrefvis) ; omega ; g' DO isun=1,nsun WRITE(*,*) QVISsQREF(isun,iaer,isize),& omegavis(isun,iaer,isize),& gvis(isun,iaer,isize) ENDDO WRITE(*,*) 'QREFvis(',iaer,isize,') = ',QREFvis(iaer,isize) WRITE(*,*) 'omegaREFvis(',iaer,isize,') = ',& omegaREFvis(iaer,isize) ! ------------------------------------------------ ENDDO !================================================================== CASE(2) domain ! INFRARED DOMAIN (idomain=2) !================================================================== DO isize=1,nsize(iaer,idomain) ! ---------------------------------- ! 2.1 solsir is not used anymore; division of Qext(IR) by solsir ! has to be done in the input ASCII files (if necessary). ! 2.2 Save the particle sizes radiustab(iaer,idomain,isize)=radiusdyn(isize) ! 2.3 Parameters tsol=215.D+0 longir(1)=long1ir longir(2)=long1co2 longir(3)=long2co2 longir(4)=long2ir longref=longrefir(iaer) epref=1.E+0 ! 2.4 Averaged optical properties (GCM channels) ! epav is /Qext(longrefir) since epref=1 ! Notice: Aerave also Computes the extinction coefficient at ! reference wavelength (called QREFvis or QREFir, ! and not epref, which is a different parameter); ! Reference wavelengths are defined for each aerosol in ! dimradmars.h. CALL aerave ( nwvl,& wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),& longref,epref,tsol,& nir-1,longir,epav,omegav,gav,& QREFir(iaer,isize),omegaREFir(iaer,isize) ) WRITE(*,*) 'QREFir(',iaer,isize,') = ',QREFir(iaer,isize) WRITE(*,*) 'omegaREFir(',iaer,isize,') = ',& omegaREFir(iaer,isize) ! 2.5 Computing /Qext(longrefvis) DO iir=1,nir-1 WRITE(*,*) 'QIRsQREFir Channel ',iir,': ',epav(iir) epav(iir)= epav(iir) * QREFir(iaer,isize) / & QREFvis(iaer,isize) ENDDO WRITE(*,*) 'Aerosol number', iaer WRITE(*,*) 'Particle size: ',radiustab(iaer,idomain,isize) WRITE(*,*) 'Rapport Solaire/IR:',& QREFvis(iaer,isize) / QREFir(iaer,isize) ! 2.6 Variable assignements ! (variables are declared by yomaer.h) ! Single scattering properties ! in each of the "nir" bands ! (cf. dimradmars.h) ! iir=1 : central 15um CO2 bands QIRsQREF(1,iaer,isize)=epav(2) omegaIR(1,iaer,isize)=omegav(2) gIR(1,iaer,isize)=gav(2) ! iir=2 : CO2 band wings ! (same properties than for central part) QIRsQREF(2,iaer,isize)=epav(2) omegaIR(2,iaer,isize)=omegav(2) gIR(2,iaer,isize)=gav(2) ! iir=3 : 9 um band [long1ir - long1co2] QIRsQREF(3,iaer,isize)=epav(1) omegaIR(3,iaer,isize)=omegav(1) gIR(3,iaer,isize)=gav(1) ! iir=4 : Far IR [long2co2 - long2ir] QIRsQREF(4,iaer,isize)=epav(3) omegaIR(4,iaer,isize)=omegav(3) gIR(4,iaer,isize)=gav(3) ! 2.7 Output display WRITE(*,*) 'AEROSOL PROPERTIES: Number ',iaer WRITE(*,*) 'Thermal IR (LW) ---->' WRITE(*,*) 'Particle size: ',radiustab(iaer,idomain,isize) WRITE(*,*) '/Qext(longrefvis) ; omega ; g' DO iir=1,nir WRITE(*,*) QIRsQREF(iir,iaer,isize),omegaIR(iir,iaer,isize),& gIR(iir,iaer,isize) ENDDO WRITE(*,*) 'CO2: /Qext(longrefvis) = ',& QIRsQREF(1,iaer,isize)*(1-omegaIR(1,iaer,isize)) WRITE(*,*) '' ENDDO ! isize (particle size) ------------------------------------- END SELECT domain !================================================================== ! 3. Deallocate temporary variables (read in ASCII files) !================================================================== DEALLOCATE(wvl) ! wvl DEALLOCATE(radiusdyn) ! radiusdyn DEALLOCATE(ep) ! ep DEALLOCATE(omeg) ! omeg DEALLOCATE(gfactor) ! g END DO ! Loop on iaer END DO ! Loop on idomain !================================================================== RETURN END