Changeset 2443 for trunk

Timestamp:

Jan 15, 2021, 5:41:46 PM (4 years ago)

Author:

abierjon

Message:

Mars GCM:
Changes in aeroptical.F90 :

• add the flag nf90_64bit_offset when creating the output file. Without it, the program couldn't deal with large files (concat) with multiple aerosols
• modify the code structure to make it a bit less memory-intensive
• add topdust in computed aerosol opacities
• add the possibility to compute absorption opacities, instead of extinction : the user must specify 'ext' or 'abs' in aeroptical.def
• update the description of aeroptical in util/README

+ adapt simu_MCS to take into account aeroptical output variables when choosing

the binning reference

+ add outputs for topdust in physiq_mod

AB

Location:

trunk/LMDZ.MARS

Files:

• README (modified) (1 diff)
• libf/phymars/physiq_mod.F (modified) (2 diffs)
• util/README (modified) (2 diffs)
• util/aeroptical.F90 (modified) (38 diffs)
• util/aeroptical.def (modified) (2 diffs)
• util/simu_MCS.F90 (modified) (1 diff)

trunk/LMDZ.MARS/README

@@ -3252 +3252 @@
 with Solar longitude as the Time variable (for localtime, the check on long_name
 attribute is still compliant with files from previous versions)
+
+== 15/01/2021 == AB
+Changes in aeroptical.F90 :
+- add the flag nf90_64bit_offset when creating the output file. Without it, the 
+  program couldn't deal with large files (concat) with multiple aerosols
+- modify the code structure to make it a bit less memory-intensive
+- add topdust in computed aerosol opacities
+- add the possibility to compute absorption opacities, instead of extinction :
+  the user must specify 'ext' or 'abs' in aeroptical.def
+- update the description of aeroptical in util/README
+
++ adapt simu_MCS to take into account aeroptical output variables when choosing
+  the binning reference
++ add outputs for topdust in physiq_mod

trunk/LMDZ.MARS/libf/phymars/physiq_mod.F

@@ -3167 +3167 @@
 
            if (rdstorm) then  ! writediagfi tendencies stormdust tracers
-             call WRITEDIAGFI(ngrid,'reffdust','reffdust',
-     &                        'm',3,rdust*ref_r0)
              call WRITEDIAGFI(ngrid,'reffstormdust','reffstormdust',
      &                        'm',3,rstormdust*ref_r0)
@@ -3226 +3224 @@
 
            if (slpwind) then
+             call WRITEDIAGFI(ngrid,'refftopdust','refftopdust',
+     &                        'm',3,rtopdust*ref_r0)
+             call WRITEDIAGFI(ngrid,'topdustq','top Dust mass mr',
+     &                        'kg/kg',3,pq(:,:,igcm_topdust_mass))
+             call WRITEDIAGFI(ngrid,'topdustN','top Dust number',
+     &                        'part/kg',3,pq(:,:,igcm_topdust_number))
              select case (trim(dustiropacity))
               case ("tes")

trunk/LMDZ.MARS/util/README

@@ -191 +191 @@
 --------------------------------------------------------------------
 
-Program to compute opacities [1/km] of all aerosols from dust, water ice
-and stormdust present in the input file, at a wavelength given by the user.
+Program to compute opacities [1/km] of all aerosols from dust, water ice,
+stormdust and topdust present in the input file, at a wavelength given by the user.
 Computation is made from the aerosol mass mixing ratios and effective radius,
 associated with air density (rho) and files containing the aerosol 
 optical properties (generally present in the GCM datadir/ directory).
+The user have to precise the type of opacity : extinction or absorption.
 
 NB : water ice effective radius is known to be a bit inaccurate in the GCM,
@@ -202 +203 @@
 input : diagfi.nc  / concat.nc / stats.nc kind of files
 
-output file is : name_of_input_file_OPA.nc
+output file is :
+ name_of_input_file_OPAext.nc for extinction opacities
+ name_of_input_file_OPAabs.nc for absorption opacities
+
 aeroptical.e < aeroptical.def

trunk/LMDZ.MARS/util/aeroptical.F90

@@ -8 +8 @@
 ! This program takes as inputs a GCM output file (diagfi,stats,concat) that
 ! contains:
-!      - the Mass Mixing Ratios of dust (dustq), water ice (h2o_ice)
-!                                 & stormdust (rdsdustq)
-!      - their effective radius (reffdust, reffice(*), reffstormdust)
+!      - the Mass Mixing Ratios of dust (dustq), water ice (h2o_ice),
+!                               stormdust (rdsdustq) & topdust (topdustq)
+!      - their effective radius (reffdust, reffice(*),
+!                           reffstormdust, refftopdust)
 !      - the atmospheric density (rho)
 ! as well as a wavelength to calibrate the opacities, and the directory
@@ -17 +18 @@
 ! It outputs a netcdf file containing the opacities [1/km] of the dust,
 ! water ice and stormdust aerosols calibrated at the prescribed wavelength.
-!
+! The type of opacity (extinction or absorption) is given by the user.
 !
 ! (*) due to a high uncertainty of reffice in the GCM, the value is asked
@@ -46 +47 @@
 integer :: ierr                               ! [netcdf] subroutine returned error code
 character(len=256) :: error_text              ! text to display in case of error
-integer :: tmpvarid                           ! temporarily stores a variable ID number
 integer :: lonlen,latlen,altlen,timelen       ! nb of grid points along longitude,latitude,altitude,time
 integer :: GCM_layers                         ! number of GCM layers
@@ -52 +52 @@
 
 logical,dimension(:),allocatable :: aerok     ! to know if the needed fields are in the GCM file
-real,dimension(:,:,:,:,:),allocatable :: mmr  ! aerosols mass mixing ratios [kg/kg]
-real,dimension(:,:,:,:,:),allocatable :: reff ! aerosols effective radii [m]
+integer,dimension(:),allocatable :: mmrvarid  ! stores a MMR variable ID number
+integer,dimension(:),allocatable :: reffvarid ! stores a reff variable ID number
+integer :: tmpvarid                           ! temporarily stores a variable ID number
+real,dimension(:,:,:,:),allocatable :: mmr    ! aerosols mass mixing ratios [kg/kg]
+real,dimension(:,:,:,:),allocatable :: reff   ! aerosols effective radii [m]
 real :: reffwice_val                          ! value given by the user for reffwice (0 if read in the GCM file) [m]
 real,dimension(:,:,:,:),allocatable :: rho    ! atmospheric density [kg.m-3]
 integer :: naerkind                           ! nb of aerosols to consider
-integer :: iaer                               ! aerosol kind index (1=dust, 2=water ice, 3=stormdust) ! + NEW AER?
+integer :: iaer                               ! aerosol kind index (1=dust,2=water ice,3=stormdust,4=topdust) ! + NEW AER?
 integer :: ilon,ilat,ialt,it                  ! Loop indices for coordinates
 
 
 ! Output file
-character(len=128) :: outfile                    ! name of the netcdf output file
-integer :: outfid                                ! [netcdf] file ID number
+character(len=128) :: outfile                  ! name of the netcdf output file
+integer :: outfid                              ! [netcdf] file ID number
 integer :: latdimout,londimout,altdimout,timedimout
-                                                 ! latdimout: stores latitude dimension ID number
-                                                 ! londimout: stores longitude dimension ID number
-                                                 ! altdimout: stores altitude dimension ID number
-                                                 ! timedimout: stores time dimension ID number
-integer :: tmpvaridout                           ! temporarily stores a variable ID number
-
-real :: wvl_val                                  ! reference wavelength of the output opacities (given by the user) [m]
-integer :: varshape(4)                           ! stores a variable shape (of dimensions' IDs)
-character(len=16) :: tmpvarname                  ! temporarily stores a variable name
-real,dimension(:,:,:,:,:),allocatable :: opa_aer ! Aerosols opacities [1/km]
-real,dimension(:),allocatable :: missval         ! Value to put in outfile when the reff is out of bounds
-                                                 ! or when there is a mising value in input file
+                                               ! latdimout: stores latitude dimension ID number
+                                               ! londimout: stores longitude dimension ID number
+                                               ! altdimout: stores altitude dimension ID number
+                                               ! timedimout: stores time dimension ID number
+integer :: tmpvaridout                         ! temporarily stores a variable ID number
+
+real :: wvl_val                                ! reference wavelength of the output opacities (given by the user) [m]
+character(len=3) :: opatype                    ! opacity type : extinction (ext) or absorption (abs) (given by the user)
+integer :: varshape(4)                         ! stores a variable shape (of dimensions' IDs)
+character(len=16) :: tmpvarname                ! temporarily stores a variable name
+real,dimension(:,:,:,:),allocatable :: opa_aer ! Aerosols opacities [1/km]
+real,dimension(:),allocatable :: missval       ! Value to put in outfile when the reff is out of bounds
+                                               ! or when there is a mising value in input file
 
 
@@ -83 +87 @@
 character(len=128) :: optpropfile          ! name of the ASCII optical properties file
 logical :: file_ok                         ! to know if the file can be opened
-integer :: file_unit = 60                  !
+integer :: file_unit = 60                  ! ASCII file ID
 
 integer :: jfile                           ! ASCII file scan index
@@ -96 +100 @@
 real,dimension(:),allocatable :: radiusdyn ! Particle size axis [m]
 real,dimension(:,:),allocatable :: Qext    ! Extinction efficiency coefficient [/]
+real,dimension(:,:),allocatable :: omeg    ! Single scattering albedo [/]
 
 
@@ -103 +108 @@
 real,dimension(2) :: reffQext            ! Qext after reff interpolation
 real :: Qext_val                         ! Qext after both interpolations
+real,dimension(2) :: reffomeg            ! omeg after reff interpolation
+real :: omeg_val                         ! omeg after both interpolations
 real,dimension(:),allocatable :: rho_aer ! Density of the aerosols [kg.m-3]
 
@@ -132 +139 @@
 READ(*,'(a)') datadir
 
+! Ask the type of opacity that has to be computed
+write(*,*)""
+WRITE(*,*) "Do you want to compute extinction or absorption opacities ? (ext/abs)"
+READ(*,*) opatype
+if ((trim(opatype).ne."ext").and.(trim(opatype).ne."abs")) then
+  write(*,*) "opatype = "//opatype
+  write(*,*) "Error: the opacity type must be ext or abs"
+  stop
+endif
+
 !===============================================================================
 ! 1. GCM FILE & OUTPUT FILE INITIALIZATIONS
@@ -149 +166 @@
 !==========================================================================
 ! --> 1.1.2 Creation of the outfile
-outfile=gcmfile(1:index(gcmfile, ".nc")-1)//"_OPA.nc"
-
-ierr=NF90_CREATE(outfile,nf90_clobber,outfid) ! NB: clobber mode=overwrite any dataset with the same file name !
+if (opatype.eq."ext") then
+  outfile=gcmfile(1:index(gcmfile, ".nc")-1)//"_OPAext.nc"
+else ! opatype.eq."abs"
+  outfile=gcmfile(1:index(gcmfile, ".nc")-1)//"_OPAabs.nc"
+endif
+
+ 
+ierr=NF90_CREATE(outfile,IOR(nf90_clobber,nf90_64bit_offset),outfid)
+! NB: clobber mode=overwrite any dataset with the same file name !
+! 64bit_offset enables the creation of large netcdf files with variables>2GB
 error_text="Error: could not create file "//trim(outfile)
 call status_check(ierr,error_text)
@@ -171 +195 @@
 
 !==========================================================================
-! 1.2 GCM file variables
+! 1.2 GCM aerosols
 !==========================================================================
 
 ! Number of aerosols
-naerkind = 3 ! dust, water ice, stormdust ! + NEW AER?
-
-! Length allocation of the variables
-allocate(mmr(naerkind,lonlen,latlen,altlen,timelen))
-allocate(reff(naerkind,lonlen,latlen,altlen,timelen))
+naerkind = 4 ! dust, water ice, stormdust, topdust ! + NEW AER?
+
+! Variables length allocation
+allocate(mmrvarid(naerkind))
+allocate(reffvarid(naerkind))
 
 ! Initialize missing value
@@ -194 +218 @@
 ! DUST MASS MIXING RATIO
 ! Check that the GCM file contains that variable
-ierr=nf90_inq_varid(gcmfid,"dustq",tmpvarid)
+ierr=nf90_inq_varid(gcmfid,"dustq",mmrvarid(1))
 error_text="Failed to find variable dustq in "//trim(gcmfile)&
             //" We'll skip the dust aerosol."
 if (ierr.ne.nf90_noerr) then
   write(*,*)trim(error_text)
-  aerok(1)=.false.
-else
-  ! Load datasets
-  ierr=NF90_GET_VAR(gcmfid,tmpvarid,mmr(1,:,:,:,:))
-  error_text="Failed to load dust mass mixing ratio"
-  call status_check(ierr,error_text)
-  write(*,*) "Dust mass mixing ratio loaded from "//trim(gcmfile)
-  ! Get missing value
-  ierr=nf90_get_att(gcmfid,tmpvarid,"missing_value",missval(1))
-  if (ierr.ne.nf90_noerr) then
-    missval(1) = 1.e+20
-  endif
+  aerok(1)=.false. 
 endif
 
@@ -216 +229 @@
 if (aerok(1)) then
   ! Check that the GCM file contains that variable
-  ierr=nf90_inq_varid(gcmfid,"reffdust",tmpvarid)
+  ierr=nf90_inq_varid(gcmfid,"reffdust",reffvarid(1))
   error_text="Failed to find variable reffdust in "//trim(gcmfile)&
               //" We'll skip the dust aerosol."
@@ -222 +235 @@
     write(*,*)trim(error_text)
     aerok(1)=.false.
-  else
-    ! Load datasets
-    ierr=NF90_GET_VAR(gcmfid,tmpvarid,reff(1,:,:,:,:))
-    error_text="Failed to load dust effective radius"
-    call status_check(ierr,error_text)
-    write(*,*) "Dust effective radius loaded from "//trim(gcmfile)
   endif
 endif
@@ -236 +243 @@
 ! WATER ICE MASS MIXING RATIO
 ! Check that the GCM file contains that variable
-ierr=nf90_inq_varid(gcmfid,"h2o_ice",tmpvarid)
+ierr=nf90_inq_varid(gcmfid,"h2o_ice",mmrvarid(2))
 error_text="Failed to find variable h2o_ice in "//trim(gcmfile)&
             //" We'll skip the water ice aerosol."
@@ -242 +249 @@
   write(*,*)trim(error_text)
   aerok(2)=.false.
-else
-  ! Load datasets
-  ierr=NF90_GET_VAR(gcmfid,tmpvarid,mmr(2,:,:,:,:))
-  error_text="Failed to load water ice mass mixing ratio"
-  call status_check(ierr,error_text)
-  write(*,*) "Water ice mass mixing ratio loaded from "//trim(gcmfile)
-  ! Get missing value
-  ierr=nf90_get_att(gcmfid,tmpvarid,"missing_value",missval(2))
-  if (ierr.ne.nf90_noerr) then
-    missval(2) = 1.e+20
-  endif
 endif
 
@@ -259 +255 @@
   IF (reffwice_val.eq.0) THEN
     ! Check that the GCM file contains that variable
-    ierr=nf90_inq_varid(gcmfid,"reffice",tmpvarid)
+    ierr=nf90_inq_varid(gcmfid,"reffice",reffvarid(2))
     error_text="Failed to find variable reffice in "//trim(gcmfile)&
                 //" We'll skip the water ice aerosol."
@@ -265 +261 @@
       write(*,*)trim(error_text)
       aerok(2)=.false.
-    else
-      ! Load datasets
-      ierr=NF90_GET_VAR(gcmfid,tmpvarid,reff(2,:,:,1,:)) ! reffice is actually a GCM 
-                                                         ! surface (column-averaged) variable
-      error_text="Failed to load water ice effective radius"
-      call status_check(ierr,error_text)
-      do ialt=2,altlen
-        reff(2,:,:,ialt,:) = reff(2,:,:,1,:) ! build up along altitude axis
-      enddo
-      write(*,*) "Water ice effective radius loaded from "//trim(gcmfile)
     endif
-  ELSE ! if reffwice_val/=0
-    reff(2,:,:,:,:) = reffwice_val
-    write(*,*) "Water ice effective radius loaded from the user input"
   ENDIF
 endif
-
-! Check if there is still at least one aerosol to compute
-IF (.NOT.ANY(aerok)) THEN
-  write(*,*) "Neither dust nor water ice are found in the file. Better stop now..."
-  stop
-ENDIF
 
 !==========================================================================
@@ -293 +270 @@
 ! STORMDUST MASS MIXING RATIO
 ! Check that the GCM file contains that variable
-ierr=nf90_inq_varid(gcmfid,"rdsdustq",tmpvarid)
+ierr=nf90_inq_varid(gcmfid,"rdsdustq",mmrvarid(3))
 error_text="Failed to find variable rdsdustq in "//trim(gcmfile)&
             //" We'll skip the stormdust aerosol."
@@ -299 +276 @@
   write(*,*)trim(error_text)
   aerok(3)=.false.
-else
-  ! Load datasets
-  ierr=NF90_GET_VAR(gcmfid,tmpvarid,mmr(3,:,:,:,:))
-  error_text="Failed to load stormdust mass mixing ratio"
-  call status_check(ierr,error_text)
-  write(*,*) "Stormdust mass mixing ratio loaded from "//trim(gcmfile)
-  ! Get missing value
-  ierr=nf90_get_att(gcmfid,tmpvarid,"missing_value",missval(3))
-  if (ierr.ne.nf90_noerr) then
-    missval(3) = 1.e+20
-  endif
 endif
 
@@ -315 +281 @@
 if (aerok(3)) then
   ! Check that the GCM file contains that variable
-  ierr=nf90_inq_varid(gcmfid,"reffstormdust",tmpvarid)
+  ierr=nf90_inq_varid(gcmfid,"reffstormdust",reffvarid(3))
   error_text="Failed to find variable reffstormdust in "//trim(gcmfile)&
-              //" We'll skip the dust aerosol."
+              //" We'll skip the stormdust aerosol."
   if (ierr.ne.nf90_noerr) then
     write(*,*)trim(error_text)
     aerok(3)=.false.
-  else
-    ! Load datasets
-    ierr=NF90_GET_VAR(gcmfid,tmpvarid,reff(3,:,:,:,:))
-    error_text="Failed to load stormdust effective radius"
-    call status_check(ierr,error_text)
-    write(*,*) "Stormdust effective radius loaded from "//trim(gcmfile)
   endif
 endif 
 
-
-!==========================================================================
-! --> 1.2.4 + NEW AER? 
-
-
-!==========================================================================
-! --> 1.2.5 ATMOSPHERIC DENSITY
+!==========================================================================
+! --> 1.2.4 TOPDUST
+
+! TOPDUST MASS MIXING RATIO
+! Check that the GCM file contains that variable
+ierr=nf90_inq_varid(gcmfid,"topdustq",mmrvarid(4))
+error_text="Failed to find variable topdustq in "//trim(gcmfile)&
+            //" We'll skip the topdust aerosol."
+if (ierr.ne.nf90_noerr) then
+  write(*,*)trim(error_text)
+  aerok(4)=.false.
+endif
+
+! TOPDUST EFFECTIVE RADIUS
+if (aerok(4)) then
+  ! Check that the GCM file contains that variable
+  ierr=nf90_inq_varid(gcmfid,"refftopdust",reffvarid(4))
+  error_text="Failed to find variable refftopdust in "//trim(gcmfile)&
+              //" We'll skip the topdust aerosol."
+  if (ierr.ne.nf90_noerr) then
+    write(*,*)trim(error_text)
+    aerok(4)=.false.
+  endif
+endif
+
+!==========================================================================
+! --> 1.2.5 + NEW AER? 
+
+
+! Check if there is still at least one aerosol to compute
+IF (.NOT.ANY(aerok)) THEN ! + NEW AER?
+  write(*,*) "No aerosol among [dust/water ice/stormdust/topdust] was found in the file. Better stop now..."
+  stop
+ENDIF
+
+!==========================================================================
+! 1.3 GCM ATMOSPHERIC DENSITY
+!==========================================================================
 
 ! Check that the GCM file contains that variable
@@ -345 +336 @@
 ! Length allocation for each dimension of the 4D variable
 allocate(rho(lonlen,latlen,altlen,timelen))
-! Load datasets
+! Load dataset
 ierr=nf90_get_var(gcmfid,tmpvarid,rho)
 error_text="Failed to load atmospheric density"
 call status_check(ierr,error_text)
-write(*,*) "Loaded atmospheric density rho from "//trim(gcmfile)
-
-!==========================================================================
-! 1.3 Some output variable's initializations
+write(*,*) "Atmospheric density rho loaded from "//trim(gcmfile)
+
+
+!==========================================================================
+! 1.4 Output variable's initializations and datasets loading
 !==========================================================================
 ! Define the ID shape of the output variables
@@ -365 +357 @@
 rho_aer(2)=920.  ! water ice
 rho_aer(3)=2500. ! stormdust
+rho_aer(4)=2500. ! topdust
 ! + NEW AER?
 
-allocate(opa_aer(naerkind,lonlen,latlen,altlen,timelen))
-
-
+DO iaer = 1, naerkind ! Loop on aerosol kind
+  if (aerok(iaer)) then ! check if this aerosol can be computed
+    write(*,*)""
+    SELECT CASE(iaer)
+    CASE(1) ! DUST
+      write(*,*)"Computing dust opacity..."
+    CASE(2) ! WATER ICE
+      write(*,*)"Computing water ice opacity..."
+    CASE(3) ! STORMDUST
+      write(*,*)"Computing stormdust opacity..."
+    CASE(4) ! TOPDUST
+      write(*,*)"Computing topdust opacity..."
+! + NEW AER?
+    END SELECT ! iaer
+    
+    ! Length allocation of the input variables
+    ALLOCATE(mmr(lonlen,latlen,altlen,timelen))
+    ALLOCATE(reff(lonlen,latlen,altlen,timelen))
+    
+    ! Datasets loading
+     ! MMR
+    ierr=NF90_GET_VAR(gcmfid,mmrvarid(iaer),mmr(:,:,:,:))
+    error_text="Failed to load mass mixing ratio"
+    call status_check(ierr,error_text)
+    write(*,*) "Mass mixing ratio loaded from "//trim(gcmfile)
+    ! Get missing value
+    ierr=nf90_get_att(gcmfid,mmrvarid(1),"missing_value",missval(1))
+    if (ierr.ne.nf90_noerr) then
+      missval(1) = 1.e+20
+    endif
+    
+     ! REFF
+    if (iaer.ne.2) then
+      ! Load dataset
+      ierr=NF90_GET_VAR(gcmfid,reffvarid(iaer),reff(:,:,:,:))
+      error_text="Failed to load effective radius"
+      call status_check(ierr,error_text)
+      write(*,*) "Effective radius loaded from "//trim(gcmfile)
+    
+    else ! water ice special case
+      IF (reffwice_val.eq.0) THEN
+        ! Load dataset
+        ierr=NF90_GET_VAR(gcmfid,reffvarid(iaer),reff(:,:,1,:)) ! reffice is actually a GCM 
+                                                              ! surface (column-averaged) variable
+        error_text="Failed to load effective radius"
+        call status_check(ierr,error_text)
+        do ialt=2,altlen
+          reff(:,:,ialt,:) = reff(:,:,1,:) ! build up along altitude axis
+        enddo
+        write(*,*) "Effective radius loaded from "//trim(gcmfile)
+      ELSE ! if reffwice_val/=0
+        reff(:,:,:,:) = reffwice_val
+        write(*,*) "Effective radius loaded from the user input"
+      ENDIF
+      
+    endif ! iaer.ne.2
+    
+  
+    ! Length allocation of the output variable
+    ALLOCATE(opa_aer(lonlen,latlen,altlen,timelen))
+    
+  
 !===============================================================================
 ! 2. OPTICAL PROPERTIES
 !===============================================================================
-
-DO iaer = 1, naerkind ! Loop on aerosol kind
-  if (aerok(iaer)) then ! check if this aerosol can be computed
-
 !==========================================================================
 ! 2.1 Open the ASCII file
@@ -391 +439 @@
         optpropfile = "optprop_icevis_n30.dat"
       CASE(3) ! STORMDUST
+        ! Dust file
+        optpropfile = "optprop_dustvis_TM_n50.dat"
+      CASE(4) ! TOPDUST
         ! Dust file
         optpropfile = "optprop_dustvis_TM_n50.dat"
@@ -407 +458 @@
         optpropfile = "optprop_iceir_n30.dat"
       CASE(3) ! STORMDUST
+        ! Dust file
+        optpropfile = "optprop_dustir_n50.dat"
+      CASE(4) ! TOPDUST
         ! Dust file
         optpropfile = "optprop_dustir_n50.dat"
@@ -468 +522 @@
     ALLOCATE(radiusdyn(nsize))      ! Aerosol radius axis
     ALLOCATE(Qext(nwvl,nsize))      ! Extinction efficiency coeff
+    ALLOCATE(omeg(nwvl,nsize))      ! Single scattering albedo
 
 !==========================================================================
@@ -516 +571 @@
               ENDIF
             ENDDO
+            jfile = jfile+1
+        CASE(4) reading2_seq ! omeg ----------------------------
+            isize = 1
+            DO WHILE (isize .le. nsize)
+              READ(file_unit,*,iostat=read_ok) scanline
+              if (read_ok.ne.0) then
+                write(*,*)'reading2_seq: Error while reading line: ',&
+                          trim(scanline)
+                write(*,*)'   of file ',&
+                          trim(datadir)//'/'//trim(optpropfile)
+                stop
+              endif
+              IF ((scanline(1:1).ne.'#').and.(scanline(1:1).ne.' ')) THEN
+                BACKSPACE(file_unit)
+                read(file_unit,*) omeg(:,isize)
+                isize = isize + 1
+              ENDIF
+            ENDDO
             endwhile = .true.
         CASE DEFAULT reading2_seq ! default --------------------
@@ -569 +642 @@
     ! Switch to netcdf define mode
     ierr=nf90_redef(outfid)
+    error_text="ERROR: Couldn't start netcdf define mode"
+    call status_check(ierr,error_text)
     write(*,*)""
     SELECT CASE (iaer)
@@ -580 +655 @@
       
       ! Write the attributes
-      ierr=nf90_put_att(outfid,tmpvaridout,"long_name","Dust extinction opacity at reference wavelength")
+      if (opatype.eq."ext") then
+        ierr=nf90_put_att(outfid,tmpvaridout,"long_name","Dust extinction opacity at reference wavelength")
+      else ! opatype.eq."abs"
+        ierr=nf90_put_att(outfid,tmpvaridout,"long_name","Dust absorption opacity at reference wavelength")
+      endif
       
     CASE(2) ! WATER ICE
@@ -591 +670 @@
       
       ! Write the attributes
-      ierr=nf90_put_att(outfid,tmpvaridout,"long_name","Water ice extinction opacity at reference wavelength")
+      if (opatype.eq."ext") then
+        ierr=nf90_put_att(outfid,tmpvaridout,"long_name","Water ice extinction opacity at reference wavelength")
+      else ! opatype.eq."abs"
+        ierr=nf90_put_att(outfid,tmpvaridout,"long_name","Water ice absorption opacity at reference wavelength")
+      endif
       
     CASE(3) ! STORMDUST
@@ -602 +685 @@
       
       ! Write the attributes
-      ierr=nf90_put_att(outfid,tmpvaridout,"long_name","Stormdust extinction opacity at reference wavelength")
+      if (opatype.eq."ext") then
+        ierr=nf90_put_att(outfid,tmpvaridout,"long_name","Stormdust extinction opacity at reference wavelength")
+      else ! opatype.eq."abs"
+        ierr=nf90_put_att(outfid,tmpvaridout,"long_name","Stormdust absorption opacity at reference wavelength")
+      endif
+
+    CASE(4) ! TOPDUST
+      ! Insert the definition of the variable
+      tmpvarname="opatopdust"
+      ierr=NF90_DEF_VAR(outfid,trim(tmpvarname),nf90_float,varshape,tmpvaridout)
+      error_text="ERROR: Couldn't create "//trim(tmpvarname)//" in the outfile"
+      call status_check(ierr,error_text)
+      write(*,*) trim(tmpvarname)//" has been created in the outfile"
+      
+      ! Write the attributes
+      if (opatype.eq."ext") then
+        ierr=nf90_put_att(outfid,tmpvaridout,"long_name","Topdust extinction opacity at reference wavelength")
+      else ! opatype.eq."abs"
+        ierr=nf90_put_att(outfid,tmpvaridout,"long_name","Topdust absorption opacity at reference wavelength")
+      endif
 
 ! + NEW AER?    
     END SELECT ! iaer
-
+    
     ierr=nf90_put_att(outfid,tmpvaridout,"units","opacity/km")
     ierr=nf90_put_att(outfid,tmpvaridout,"refwavelength",wvl_val)
@@ -614 +716 @@
     ! End netcdf define mode
     ierr=nf90_enddef(outfid)
-
+    error_text="ERROR: Couldn't end netcdf define mode"
+    call status_check(ierr,error_text)
+      
 !==========================================================================
 ! 3.2 Computation of the opacities
@@ -624 +728 @@
          ! Get the optpropfile reff values encompassing the GCM reff
          isize=1
-         do while((isize.le.nsize).and.(radiusdyn(isize).lt.reff(iaer,ilon,ilat,ialt,it)))
+         do while((isize.le.nsize).and.(radiusdyn(isize).lt.reff(ilon,ilat,ialt,it)))
            isize=isize+1
          enddo
-         if ((isize.gt.nsize).or.(radiusdyn(1).gt.reff(iaer,ilon,ilat,ialt,it))) then
-!           write(*,*) "WARNING: the GCM reff (",reff(iaer,ilon,ilat,ialt,it),"m) is out of the bounds"
+         if ((isize.gt.nsize).or.(radiusdyn(1).gt.reff(ilon,ilat,ialt,it))) then
+!           write(*,*) "WARNING: the GCM reff (",reff(ilon,ilat,ialt,it),"m) is out of the bounds"
 !           write(*,*) "of the file ",trim(datadir)//'/'//trim(optpropfile)
 !           write(*,*) "(reff_inf=",radiusdyn(1),"m ; reff_sup=",radiusdyn(nsize),"m)"
@@ -634 +738 @@
 
            ! NB: this should especially handle cases when reff=0
-           opa_aer(iaer,ilon,ilat,ialt,it)=missval(iaer)
+           opa_aer(ilon,ilat,ialt,it)=missval(iaer)
            
-         else if (mmr(iaer,ilon,ilat,ialt,it).eq.missval(iaer)) then
+         else if (mmr(ilon,ilat,ialt,it).eq.missval(iaer)) then
            ! if there is a missing value in input file
-           opa_aer(iaer,ilon,ilat,ialt,it)=missval(iaer)
+           opa_aer(ilon,ilat,ialt,it)=missval(iaer)
            
          else
-           if (radiusdyn(isize).eq.reff(iaer,ilon,ilat,ialt,it)) then
-             ! if the GCM reff is in the optpropfile
+           if (radiusdyn(isize).eq.reff(ilon,ilat,ialt,it)) then
+             ! if the GCM reff is exactly in the optpropfile
              isize1 = isize
              isize2 = isize
@@ -650 +754 @@
            endif
          
+           ! Qext COMPUTATION
            ! Linear interpolation in effective radius
            if (isize2.ne.isize1) then
-             coeff = (reff(iaer,ilon,ilat,ialt,it)-radiusdyn(isize1))/(radiusdyn(isize2)-radiusdyn(isize1))
+             coeff = (reff(ilon,ilat,ialt,it)-radiusdyn(isize1))/(radiusdyn(isize2)-radiusdyn(isize1))
            else
              coeff = 0.
@@ -658 +763 @@
            reffQext(1) = Qext(iwvl1,isize1)+coeff*(Qext(iwvl1,isize2)-Qext(iwvl1,isize1))
            reffQext(2) = Qext(iwvl2,isize1)+coeff*(Qext(iwvl2,isize2)-Qext(iwvl2,isize1))
-           
            ! Linear interpolation in wavelength
            if (iwvl2.ne.iwvl1) then
@@ -667 +771 @@
            Qext_val = reffQext(1)+coeff*(reffQext(2)-reffQext(1))
            
-           ! Computation of the opacity [1/km]
-           opa_aer(iaer,ilon,ilat,ialt,it)= 750.*Qext_val*mmr(iaer,ilon,ilat,ialt,it)*rho(ilon,ilat,ialt,it)&
-                                               / ( rho_aer(iaer) * reff(iaer,ilon,ilat,ialt,it) )
-                    
+           ! COMPUTATION OF THE EXTINCTION OPACITY [1/km]
+           opa_aer(ilon,ilat,ialt,it)= 750.*Qext_val*mmr(ilon,ilat,ialt,it)*rho(ilon,ilat,ialt,it)&
+                                               / ( rho_aer(iaer) * reff(ilon,ilat,ialt,it) )
+             
+             
+           if (opatype.eq."abs") then
+             ! omeg COMPUTATION
+             ! Linear interpolation in effective radius
+             if (isize2.ne.isize1) then
+               coeff = (reff(ilon,ilat,ialt,it)-radiusdyn(isize1))/(radiusdyn(isize2)-radiusdyn(isize1))
+             else
+               coeff = 0.
+             endif
+             reffomeg(1) = omeg(iwvl1,isize1)+coeff*(omeg(iwvl1,isize2)-omeg(iwvl1,isize1))
+             reffomeg(2) = omeg(iwvl2,isize1)+coeff*(omeg(iwvl2,isize2)-omeg(iwvl2,isize1))
+             ! Linear interpolation in wavelength
+             if (iwvl2.ne.iwvl1) then
+               coeff = (wvl_val-wvl(iwvl1))/(wvl(iwvl2)-wvl(iwvl1))
+             else
+               coeff = 0.
+             endif
+             omeg_val = reffomeg(1)+coeff*(reffomeg(2)-reffomeg(1))
+
+             ! COMPUTATION OF THE ABSORPTION OPACITY [1/km]
+             opa_aer(ilon,ilat,ialt,it)= (1 - omeg_val) * opa_aer(ilon,ilat,ialt,it)
+           endif ! opatype.eq."abs"
+                  
          endif
        enddo ! it
@@ -680 +807 @@
 ! 3.3 Writing in the output file
 !==========================================================================
-    ierr = NF90_PUT_VAR(outfid, tmpvaridout, opa_aer(iaer,:,:,:,:))
+    ierr = NF90_PUT_VAR(outfid, tmpvaridout, opa_aer(:,:,:,:))
     error_text="Error: could not write "//trim(tmpvarname)//" data in the outfile"
     call status_check(ierr,error_text)
@@ -687 +814 @@
 ! 3.4 Prepare next loop
 !==========================================================================
+    DEALLOCATE(mmr)
+    DEALLOCATE(reff)
+    DEALLOCATE(opa_aer)
     DEALLOCATE(wvl)
     DEALLOCATE(radiusdyn)
     DEALLOCATE(Qext)
+    DEALLOCATE(omeg)
     
   endif ! if aerok(iaer)

trunk/LMDZ.MARS/util/aeroptical.def

@@ -3 +3 @@
 5.e-7
 ~/datadir
-
+ext
 
 
@@ -14 +14 @@
 3) The wavelength at which the output opacities are calibrated (in meters)
 4) Directory where the optical properties files are
+5) Opacity type : extinction (ext) or absorption (abs)
 

trunk/LMDZ.MARS/util/simu_MCS.F90

@@ -1206 +1206 @@
 !===============================================================================
     ! Observer variables to extract :
-    IF ((trim(GCMvarname).eq."dsodust").or.(trim(GCMvarname).eq."dso") &
-      .or.(trim(GCMvarname).eq."dustq").or.(trim(GCMvarname).eq."dust") &
-      .or.(trim(GCMvarname).eq."reffdust")) THEN
+    IF ((index(GCMvarname,"dust").ne.0).or.(index(GCMvarname,"dso").ne.0)) THEN
+      ! if the variable name contains "dust" or "dso". Especially for the targeted variables :
+      ! dustq,dustN,dsodust,reffdust,opadust, and their equivalents for stormdust & topdust
       OBSvarname  = "dust"
       numbin_name = "numbindust"
     ELSE IF ((trim(GCMvarname).eq."h2o_ice").or.(trim(GCMvarname).eq."wice") &
-           .or.(trim(GCMvarname).eq."reffice")) THEN
+           .or.(trim(GCMvarname).eq."reffice").or.(trim(GCMvarname).eq."opawice")) THEN
       OBSvarname  = "wice"
       numbin_name = "numbinwice"