subroutine suaer_corrk ! inputs use radinc_h, only: L_NSPECTV,nsizemax,naerkind use radcommon_h, only: blamv,lamrefvis use datafile_mod, only: datadir, aerdir ! outputs use radcommon_h, only: QVISsQREF,omegavis,gvis use radcommon_h, only: radiustab,nsize,tstellar use radcommon_h, only: qrefvis use aerosol_mod 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 ! !================================================================== ! Optical properties (read in external ASCII files) INTEGER,SAVE :: nwvl ! Number of wavelengths in ! the domain (VIS or IR), read by master ! REAL :: solsir ! visible to infrared ratio ! ! (tau_0.67um/tau_9um) REAL, DIMENSION(:),& ALLOCATABLE, SAVE :: wvl ! Wavelength axis, read by master REAL, DIMENSION(:),& ALLOCATABLE, SAVE :: radiusdyn ! Particle size axis, read by master REAL, DIMENSION(:,:),& ALLOCATABLE, SAVE :: ep,& ! Extinction coefficient Qext, read by master omeg,& ! Single Scattering Albedo, read by master gfactor ! Assymetry Factor, read by master ! 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(lamrefvis) if epref=1) REAL omegavVI(L_NSPECTV) ! Average single scattering albedo REAL gavVI(L_NSPECTV) ! Average assymetry parameter logical forgetit ! use francois' data? integer iwvl, ia ! Local saved variables: character(len=100) :: message character(len=10),parameter :: subname="suaercorrk" CHARACTER(LEN=50), DIMENSION(naerkind,2), SAVE :: file_id !$OMP THREADPRIVATE(file_id) !---- Please indicate the names of the optical property files below ! Please also choose the reference wavelengths of each aerosol !--------- README TO UNDERSTAND WHAT FOLLOWS (JVO 20) ------- ! The lamref variables comes from the Martian model ! where the visible one is the one used for computing ! and the IR one is used to output scaled opacity to ! match instrumental data ... This is done (at least ! for now) in the generic, so lamrefir is dummy*! ! The important one is the VISIBLE one as it will be used ! to rescale the values in callcork.F90 assuming 'aerosol' is ! at this visible reference wavelenght. ! *Actually if you change lamrefir here there is a ! slight sensitvity in the outputs because of some ! machine precision differences that amplifys and lead ! up to 10-6 differences in the radiative balance... ! This could be good to clean the code but would require ! a lot of modifs and to take care ! !-------------------------------------------------------------- ! VENUS, ONLY VISIBLE do iaer = 1, naerkind ! Loop on aerosol kind if (iaer.eq.iaero_venus1) then print*, 'naerkind= venus1', iaer ! visible file_id(iaer,1) = 'optprop_h2so4vis_n50.dat' lamrefvis(iaer)=1.5E-6 ! no idea, must find endif if (iaer.eq.iaero_venus2) then print*, 'naerkind= venus2', iaer ! visible file_id(iaer,1) = 'optprop_h2so4vis_n50.dat' lamrefvis(iaer)=1.5E-6 ! no idea, must find endif if (iaer.eq.iaero_venus2p) then print*, 'naerkind= venus2p', iaer ! visible file_id(iaer,1) = 'optprop_h2so4vis_n50.dat' lamrefvis(iaer)=1.5E-6 ! no idea, must find endif if (iaer.eq.iaero_venus3) then print*, 'naerkind= venus3', iaer ! visible file_id(iaer,1) = 'optprop_h2so4vis_n50.dat' lamrefvis(iaer)=1.5E-6 ! no idea, must find endif if (iaer.eq.iaero_venusUV) then print*, 'naerkind= venusUV', iaer ! visible file_id(iaer,1) = 'optprop_venusUVvis.dat' lamrefvis(iaer)=3.5E-7 ! Haus et al. 2015 endif enddo !------------------------------------------------------------------ ! Initializations: radiustab(:,:,:) = 0.0 gvis(:,:,:) = 0.0 omegavis(:,:,:) = 0.0 QVISsQREF(:,:,:) = 0.0 DO iaer = 1, naerkind ! Loop on aerosol kind ! DO idomain = 1, 2 ! Loop on radiation domain (VIS or IR) idomain=1 !================================================================== ! 1. READ OPTICAL PROPERTIES !================================================================== ! 1.1 Open the ASCII file !$OMP MASTER INQUIRE(FILE=TRIM(datadir)//'/'//TRIM(aerdir)//& '/'//TRIM(file_id(iaer,idomain)),& EXIST=file_ok) IF (file_ok) THEN OPEN(UNIT=file_unit,& FILE=TRIM(datadir)//'/'//TRIM(aerdir)//& '/'//TRIM(file_id(iaer,idomain)),& FORM='formatted') ELSE ! In ye old days these files were stored in datadir; ! let's be retro-compatible INQUIRE(FILE=TRIM(datadir)//& '/'//TRIM(file_id(iaer,idomain)),& EXIST=file_ok) IF (file_ok) THEN OPEN(UNIT=file_unit,& FILE=TRIM(datadir)//& '/'//TRIM(file_id(iaer,idomain)),& FORM='formatted') ENDIF ENDIF IF(.NOT.file_ok) THEN write(*,*)'suaer_corrk: Problem 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' write(*,*)'2) If ',& TRIM(file_id(iaer,idomain)),& ' is a LMD reference datafile, it' write(*,*)' can be obtained online at:' message=' http://www.lmd.jussieu.fr/~lmdz/planets/LMDZ.GENERIC/datagcm/' call abort_physic(subname,message,1) 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 ! ---------------------------- message='readoptprop: Error while loading optical properties.' call abort_physic(subname,message,1) 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.iaero_co2).OR.(iaer.NE.iaero_co2)) THEN endwhile = .true. ENDIF CASE(6) reading2_seq endwhile = .true. CASE DEFAULT reading2_seq ! ---------------------------- message='readoptprop: Error while loading optical properties.' call abort_physic(subname,message,1) END SELECT reading2_seq ! ============================== ENDIF ENDDO ! 1.4 Close the file CLOSE(file_unit) !$OMP END MASTER !$OMP BARRIER !================================================================== ! 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)) ! 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 ! Only visible for Venus ! END SELECT domain !======================================================================== ! 3. Deallocate temporary variables that were read in the ASCII files !======================================================================== !$OMP BARRIER !$OMP MASTER IF (ALLOCATED(wvl)) DEALLOCATE(wvl) ! wvl IF (ALLOCATED(radiusdyn)) DEALLOCATE(radiusdyn) ! radiusdyn IF (ALLOCATED(ep)) DEALLOCATE(ep) ! ep IF (ALLOCATED(omeg)) DEALLOCATE(omeg) ! omeg IF (ALLOCATED(gfactor)) DEALLOCATE(gfactor) ! g !$OMP END MASTER !$OMP BARRIER END DO ! Loop on iaer ! END DO ! Loop on idomain !======================================================================== RETURN END subroutine suaer_corrk