subroutine suaer_corrk ! inputs use radinc_h, only: L_NSPECTI,L_NSPECTV,naerkind,nsizemax use radcommon_h, only: blamv,blami,lamrefir,lamrefvis use datafile_mod, only: datadir, aerdir, hazeprop use aerosol_mod ! outputs use radcommon_h, only: QVISsQREF,omegavis,gvis,QIRsQREF,omegair,gir use radcommon_h, only: radiustab,nsize,tstellar use radcommon_h, only: qrefvis,qrefir,omegarefvis,omegarefir implicit none !================================================================== ! Purpose ! ------- ! Initialize all radiative aerosol properties ! ! Notes ! ----- ! Reads the optical properties -> Mean -> Variable assignment ! (ASCII files) (see radcommon_h.F90) ! wvl(nwvl) longsun ! ep(nwvl) epav QVISsQREF(L_NSPECTV) ! omeg(nwvl) omegav omegavis(L_NSPECTV) ! gfactor(nwvl) gav gvis(L_NSPECTV) ! ! 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 ! - Add water ice clouds ! MODIF R. Wordsworth (2009) ! - generalisation to arbitrary spectral bands ! MODIF Tanguy BERTRAND (2017) : adaptation Pluto !================================================================== #include "dimensions.h" #include "dimphys.h" #include "callkeys.h" #include "tracer.h" ! Optical properties (read in external ASCII files) INTEGER,SAVE :: 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 :: ilw ! longwave index INTEGER :: isw ! shortwave 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 ! I/O of "aerave" (subroutine that spectrally averages ! the single scattering parameters) REAL lamref ! reference wavelengths REAL epref ! reference extinction ep REAL epavVI(L_NSPECTV) ! Average ep (= /Qext(lamref) if epref=1) REAL omegavVI(L_NSPECTV) ! Average single scattering albedo REAL gavVI(L_NSPECTV) ! Average assymetry parameter REAL epavIR(L_NSPECTI) ! Average ep (= /Qext(lamref) if epref=1) REAL omegavIR(L_NSPECTI) ! Average single scattering albedo REAL gavIR(L_NSPECTI) ! Average assymetry parameter logical forgetit ! use francois' data? integer iwvl ! Local saved variables: CHARACTER(LEN=60), 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 ! NAERKIND dans radinc_h : = 1 ! Flag haze doit etre dans callcorrk... write(*,*) 'naerkind=',naerkind do iaer=1,naerkind ! naerkind=1: haze if (iaer.eq.1) then ! Only one table of optical properties : write(*,*)'Suaer haze optical properties, using: ', & TRIM(hazeprop) ! visible file_id(iaer,1)=TRIM(hazeprop) ! infrared file_id(iaer,2)=file_id(iaer,1) lamrefvis(iaer)=0.607E-6 ! reference wavelength for opacity vis pivot wavelength Cheng et al 2017 lamrefir(iaer)=2.E-6 ! reference wavelength for opacity IR (in the LEISA range) endif enddo IF (naerkind .gt. 1) THEN print*,'naerkind = ',naerkind print*,'but we only have data for 1 type, exiting.' call abort ENDIF !------------------------------------------------------------------ ! Initializations: radiustab(:,:,:) = 0.0 gvis(:,:,:) = 0.0 omegavis(:,:,:) = 0.0 QVISsQREF(:,:,:) = 0.0 gir(:,:,:) = 0.0 omegair(:,:,:) = 0.0 QIRsQREF(:,:,:) = 0.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=TRIM(datadir)//'/'//TRIM(aerdir)//& '/'//TRIM(file_id(iaer,idomain)),& EXIST=file_ok) write(*,*)' opening file=',TRIM(datadir)//'/'//TRIM(aerdir)//& '/'//TRIM(file_id(iaer,idomain)) IF (file_ok) THEN OPEN(UNIT=file_unit,& FILE=TRIM(datadir)//'/'//TRIM(aerdir)//& '/'//TRIM(file_id(iaer,idomain)),& FORM='formatted') ENDIF IF(.NOT.file_ok) THEN write(*,*)'suaer_corrk: Pb opening ',& TRIM(file_id(iaer,idomain)) write(*,*)'It should be in: ',& TRIM(datadir)//'/'//TRIM(aerdir) write(*,*)'1) You can set this directory address ',& 'in callphys.def with:' write(*,*)' datadir = /absolute/path/to/datagcm' CALL ABORT ENDIF ! 1.2 Allocate the optical property table 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) reading1_seq: SELECT CASE (jfile) ! ==================== CASE(1) reading1_seq ! nwvl ---------------------------- read(file_unit,*) nwvl jfile = jfile+1 CASE(2) reading1_seq ! nsize --------------------------- read(file_unit,*) nsize(iaer,idomain) 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 jfile = jfile+1 IF ((idomain.NE.1).OR.(iaer.NE.1).OR.(jfile.EQ.6)) THEN endwhile = .true. ENDIF CASE(6) reading2_seq endwhile = .true. CASE DEFAULT reading2_seq ! ---------------------------- WRITE(*,*) 'readoptprop: ',& 'Error while loading optical properties.' CALL ABORT END SELECT reading2_seq ! ============================== ENDIF ENDDO ! 1.4 Close the file CLOSE(file_unit) ! 1.5 Convert wvl to metres (needed in aerave_new) ! not needed if already the case ! do iwvl=1,nwvl ! wvl(iwvl)=wvl(iwvl)*1.e-6 ! enddo !================================================================== ! 2. AVERAGED PROPERTIES AND VARIABLE ASSIGNMENTS !================================================================== domain: SELECT CASE (idomain) !================================================================== CASE(1) domain ! VISIBLE DOMAIN (idomain=1) !================================================================== lamref=lamrefvis(iaer) epref=1.E+0 DO isize=1,nsize(iaer,idomain) ! Save the particle sizes radiustab(iaer,idomain,isize)=radiusdyn(isize) ! Averaged optical properties (GCM channels) CALL aerave_new ( nwvl,& wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),& lamref,epref,tstellar,& L_NSPECTV,blamv,epavVI,& omegavVI,gavVI,QREFvis(iaer,isize),omegaREFir(iaer,isize)) ! Variable assignments (declared in radcommon) DO isw=1,L_NSPECTV QVISsQREF(isw,iaer,isize)=epavVI(isw) gvis(isw,iaer,isize)=gavVI(isw) omegavis(isw,iaer,isize)=omegavVI(isw) END DO ENDDO !================================================================== CASE(2) domain ! INFRARED DOMAIN (idomain=2) !================================================================== DO isize=1,nsize(iaer,idomain) ! ---------------------------------- lamref=lamrefir(iaer) epref=1.E+0 ! Save the particle sizes radiustab(iaer,idomain,isize)=radiusdyn(isize) ! radiustab is the table of radius from the table (iaer,idomain,isize) ! Averaged optical properties (GCM channels) ! epav is /Qext(lamrefir) since epref=1 ! Note: aerave also computes the extinction coefficient at ! the reference wavelength. This is called QREFvis or QREFir ! (not epref, which is a different parameter). ! Reference wavelengths SHOULD BE defined for each aerosol in ! radcommon_h.F90 CALL aerave_new ( nwvl,& wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),& lamref,epref,tplanet,& L_NSPECTI,blami,epavIR,& omegavIR,gavIR,QREFir(iaer,isize),omegaREFir(iaer,isize)) ! Variable assignments (declared in radcommon) DO ilw=1,L_NSPECTI QIRsQREF(ilw,iaer,isize)=epavIR(ilw) gir(ilw,iaer,isize)=gavIR(ilw) omegair(ilw,iaer,isize)=omegavIR(ilw) END DO ENDDO ! isize (particle size) ------------------------------------- END SELECT domain !======================================================================== ! 3. Deallocate temporary variables that were read in the 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 subroutine suaer_corrk