Changeset 2830 for trunk/LMDZ.MARS

Timestamp:

Nov 23, 2022, 4:09:32 PM (2 years ago)

Author:

abierjon

Message:

Mars GCM:
The program util/aeroptical.F90 can now compute column-integrated optical depths
of the aerosols, in addition to the opacity profiles. The user has to specify it
('yes' or 'no') in aeroptical.def
Detailed changes :

• update of the aeroptical.def file to ask the user if the column optical depth should be computed (non-retrocompatible change)
• add in the init2 subroutine a computation of the layers' height delta_z, with adaptable method depending on the availability of some variables in the input file and on wether or not the input file has been zrecasted before. Preliminary validation shows that these different methods yield a +/-3% error on the output column optical depth tau_[aer]. The delta_z variable is also written in the output file.
• add a log file for warnings in the interpolation subroutine in aeropt_mod.F90 + add some comments in the code
• add the ouput "zzlev" (= interlayer altitudes) in libf/phymars/physiq_mod.F so that it can be directly used by aeroptical.F90 to compute delta_z

AB

Location:

trunk/LMDZ.MARS

Files:

• README (modified) (1 diff)
• libf/phymars/physiq_mod.F (modified) (1 diff)
• util/README (modified) (1 diff)
• util/aeropt_mod.F90 (modified) (6 diffs)
• util/aeroptical.F90 (modified) (19 diffs)
• util/aeroptical.def (modified) (2 diffs)

trunk/LMDZ.MARS/README

@@ -3809 +3809 @@
 == 23/11/2022 == RV
 Cleaning of phyredem0. Argument passed to the subroutines were unused.
+
+== 23/11/2022 == AB
+The program util/aeroptical.F90 can now compute column-integrated optical depths
+of the aerosols, in addition to the opacity profiles. The user has to specify it
+('yes' or 'no') in aeroptical.def
+Detailed changes :
+ - update of the aeroptical.def file to ask the user if the column optical depth 
+   should be computed (non-retrocompatible change)
+ - add in the init2 subroutine a computation of the layers' height delta_z, with 
+   adaptable method depending on the availability of some variables in the input file
+   and on wether or not the input file has been zrecasted before. Preliminary validation
+   shows that these different methods yield a +/-3% error on the output column optical
+   depth tau_[aer]. The delta_z variable is also written in the output file.
+ - add a log file for warnings in the interpolation subroutine in aeropt_mod.F90
+   + add some comments in the code
+ - add the ouput "zzlev" (= interlayer altitudes) in libf/phymars/physiq_mod.F
+   so that it can be directly used by aeroptical.F90 to compute delta_z

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

@@ -2985 +2985 @@
 c        call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",3,zplay)
 c        call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",3,zplev)
+         call WRITEDIAGFI(ngrid,"zzlev","Interlayer altitude",
+     &                    "m",3,zzlev(:,1:nlayer))
          call writediagfi(ngrid,"pphi","Geopotential","m2s-2",3,
      &                    pphi)

trunk/LMDZ.MARS/util/README

@@ -192 +192 @@
 --------------------------------------------------------------------
 
-Program to compute opacities [1/km] of aerosols at a wavelength given by the user.
+Program to compute opacity profiles opa_[aer] (1/km) of aerosols
+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.
+The user can also ask to compute the column-integrated optical depth tau_[aer].
 
 NB : this program requires to compile the module aeropt_mod.F90

trunk/LMDZ.MARS/util/aeropt_mod.F90

@@ -12 +12 @@
  real,dimension(:,:),save,allocatable :: Qext    ! Extinction efficiency coefficient [/]
  real,dimension(:,:),save,allocatable :: omeg    ! Single scattering albedo [/]
+ 
+ integer,save :: aeroptlogfileID = 100           ! Log file's unit ID
 
  CONTAINS 
@@ -46 +48 @@
   integer :: read_ok                         ! to know if the line is well read
   integer :: iwvl,isize                      ! Wavelength and Particle size loop indices
+  
+  !==============================================================================
+  ! aeropt_mod module variables used here:
+  !==============================================================================
+  ! nwvl,nsize,wvl,radiusdyn,Qext,omeg
 
   !==========================================================================
@@ -228 +235 @@
   real,dimension(2) :: reffQext         ! Qext after reff interpolation
   real,dimension(2) :: reffomeg         ! omeg after reff interpolation
+  
+  !==============================================================================
+  ! aeropt_mod module variables used here:
+  !==============================================================================
+  ! nwvl,nsize,wvl,radiusdyn,Qext,omeg,aeroptlogfileID
 
   !==========================================================================
@@ -262 +274 @@
   ENDDO
 
-  if ((radiusdyn(nsize).lt.reff_val).or.(radiusdyn(1).gt.reff_val)) then
-    write(*,*) "WARNING: the effective radius (",reff_val,"m) is out of the bounds"
-    write(*,*) "of the optical properties' tables"
-    write(*,*) "(reff_inf=",radiusdyn(1),"m ; reff_sup=",radiusdyn(nsize),"m)"
-    write(*,*) "A missing value will be written for opacity."
-
-    ! NB: this should especially handle cases when reff=0
+  if (reff_val.eq.missval) then
+    ! if the effective radius is a missing value,
+    ! put a missing value for the optical properties
+    Qext_val = missval
+    omeg_val = missval
+!    write(*,*) "reff missing value detected"
+    
+  else if ((radiusdyn(nsize).lt.reff_val).or.(radiusdyn(1).gt.reff_val)) then
+    write(aeroptlogfileID,*) "WARNING: the effective radius (",reff_val,"m) is out of the bounds"
+    write(aeroptlogfileID,*) "of the optical properties' tables"
+    write(aeroptlogfileID,*) "(reff_inf=",radiusdyn(1),"m ; reff_sup=",radiusdyn(nsize),"m)"
+    write(aeroptlogfileID,*) "A missing value will be written for optical properties."
+
     Qext_val = missval
     omeg_val = missval
@@ -333 +351 @@
 
   implicit none
+  
+  !==============================================================================
+  ! aeropt_mod module variables used here:
+  !==============================================================================
+  ! wvl,radiusdyn,Qext,omeg
 
   if (allocated(wvl)) deallocate(wvl)
@@ -342 +365 @@
  
 !*******************************************************************************
+ SUBROUTINE create_logfile(filename)
+  !==============================================================================
+  ! Purpose:
+  ! Deallocate module variables
+  !==============================================================================
+
+  implicit none
+  
+  !==============================================================================
+  ! Arguments:
+  !==============================================================================
+  character(len=128),intent(in) :: filename ! name of the aeropt log file
+  !==============================================================================
+  ! aeropt_mod module variables used here:
+  !==============================================================================
+  ! aeroptlogfileID
+  
+  ! Create output file
+  write(*,*) "Creating aeropt log file: ",trim(filename)
+  write(*,*) "It will overwrite any existing file with the same name."
+  OPEN(UNIT=aeroptlogfileID,FILE=trim(filename),FORM='formatted',STATUS='replace')
+  ! NB: STATUS='replace' will overwrite any existing file with the same name
+  
+ END SUBROUTINE create_logfile 
+ 
+!*******************************************************************************
 
 END MODULE aeropt_mod

trunk/LMDZ.MARS/util/aeroptical.F90

@@ -41 +41 @@
 real :: wvl_val                                                      ! reference wavelength of the output opacities [m]
 character(len=3) :: opatype                                          ! opacity type : extinction (ext) or absorption (abs)
+character(len=3) :: compute_col                                      ! does the user want to compute the column-integrated optical depth [yes/no]
 character(len=50),dimension(20) :: name_aer,mmrname_aer,reffname_aer ! aerosols' name, MMR varname, reff varname (or value)
 integer,dimension(20) :: reffval_ok                                  ! to know if the user gave a variable name or a value for reff_aer
@@ -50 +51 @@
 
 ! GCM file
-character(len=128) :: gcmfile                 ! name of the netcdf GCM input file
-integer :: gcmfid                             ! [netcdf] GCM input file ID number
-integer :: ierr                               ! [netcdf] subroutine returned error code
-character(len=256) :: error_text              ! text to display in case of error
-integer :: lonlen,latlen,altlen,timelen       ! nb of grid points along longitude,latitude,altitude,time
-integer :: GCM_layers                         ! number of GCM layers
-integer :: layerdimout,interlayerdimout       ! "GCM_layers" and "GCM_layers+1" IDs
-
-logical,dimension(:),allocatable :: aerok     ! to know if the needed fields are in the GCM file
-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,dimension(:,:,:,:),allocatable :: rho    ! atmospheric density [kg.m-3]
-integer :: iaer                               ! Loop index for aerosols
-integer :: ilon,ilat,ialt,it                  ! Loop indices for coordinates
+character(len=128) :: gcmfile                  ! name of the netcdf GCM input file
+integer :: gcmfid                              ! [netcdf] GCM input file ID number
+integer :: ierr                                ! [netcdf] subroutine returned error code
+character(len=256) :: error_text               ! text to display in case of error
+integer :: lonlen,latlen,altlen,timelen        ! nb of grid points along longitude,latitude,altitude,time
+integer :: GCM_layers                          ! number of GCM layers
+integer :: layerdimout,interlayerdimout        ! "GCM_layers" and "GCM_layers+1" IDs
+
+logical,dimension(:),allocatable :: aerok      ! to know if the needed fields are in the GCM file
+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,dimension(:,:,:,:),allocatable :: rho     ! atmospheric density [kg.m-3]
+integer :: iaer                                ! Loop index for aerosols
+integer :: ilon,ilat,ialt,it                   ! Loop indices for coordinates
 
 ! Opacity computation
+character(len=128) :: aeroptlogfile            ! name of the aeropt_mod log file
 real :: Qext_val                               ! Qext after interpolations
 real :: omeg_val                               ! omeg after interpolations
@@ -84 +86 @@
 character(len=16) :: tmpvarname                ! temporarily stores a variable name
 character(len=128) :: tmplongname              ! temporarily stores a variable long_name attribute
-real,dimension(:,:,:,:),allocatable :: opa_aer ! Aerosols opacities [1/km]
+real,dimension(:,:,:,:),allocatable :: opa_aer ! Aerosols opacities dtau/dz [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
+integer :: tauvarshape(3)                      ! stores the column variable shape (of dimensions' IDs)
+real,dimension(:,:,:,:),allocatable :: delta_z ! Layers' height for column computation [km]
+real,dimension(:,:,:),allocatable :: tau_aer   ! Aerosols column-integrated optical depth [/]
 
 
@@ -96 +101 @@
 !3) The wavelength at which the output opacities are calibrated (in meters)
 !4) Opacity type : extinction (ext) or absorption (abs)
-!5) aerosol_name,aerosol_mmr(variable name),aerosol_reff(variable name or single value in m),aerosol_rho(value in kg/m3),optpropfile_name
+!5) Does the user want to compute the column-integrated optical depth? [yes/no]
 !6) aerosol_name,aerosol_mmr(variable name),aerosol_reff(variable name or single value in m),aerosol_rho(value in kg/m3),optpropfile_name
+!7) aerosol_name,aerosol_mmr(variable name),aerosol_reff(variable name or single value in m),aerosol_rho(value in kg/m3),optpropfile_name
 !...
 !N) aerosol_name,aerosol_mmr(variable name),aerosol_reff(variable name or single value in m),aerosol_rho(value in kg/m3),optpropfile_name
@@ -129 +135 @@
 endif
 
-
-
-! 5-N) Ask for aerosols to compute to the user
+! 5) Computation of the column-integrated optical depth?
+write(*,*)""
+WRITE(*,*) "Do you want to compute the column-integrated optical depth ? (yes/no)"
+READ(*,*) compute_col
+if ((trim(compute_col).ne."yes").and.(trim(compute_col).ne."no")) then
+  write(*,*) "compute_col = "//compute_col
+  write(*,*) "Error: please answer yes or no"
+  stop
+endif
+
+
+! 6-N) Ask for aerosols to compute to the user
 write(*,*)""
 write(*,*) "For which aerosols do you want to compute the opacity?"
@@ -211 +226 @@
              GCM_layers,layerdimout,interlayerdimout)
 
+! Length allocation for each dimension of the 4D variable delta_z
+allocate(delta_z(lonlen,latlen,altlen,timelen))
+
 ! Initialize output file's aps,bps,ap,bp and phisinit variables
-call init2(gcmfid,lonlen,latlen,altlen,GCM_layers,&
-           outfid,londimout,latdimout,altdimout,&
-           layerdimout,interlayerdimout)
+! + compute delta_z for the column integration
+call init2(gcmfid,lonlen,latlen,altlen,timelen,GCM_layers,&
+           outfid,londimout,latdimout,altdimout,timedimout,&
+           layerdimout,interlayerdimout,&
+           compute_col,delta_z)
 
 !==========================================================================
@@ -297 +317 @@
 varshape(4)=timedimout
 
+tauvarshape(1)=londimout
+tauvarshape(2)=latdimout
+tauvarshape(3)=timedimout
+
 
 ! Loop on all aerosols
@@ -333 +357 @@
     
   
-    ! Length allocation of the output variable
+    ! Length allocation of the output variables
     ALLOCATE(opa_aer(lonlen,latlen,altlen,timelen))
     
+    if (trim(compute_col).eq."yes") then
+      ALLOCATE(tau_aer(lonlen,latlen,timelen))
+      tau_aer(:,:,:) = 0.
+    endif ! trim(compute_col).eq."yes"
+    
+    
 !==========================================================================
 ! 2.2 Reading optical properties
 !==========================================================================
 
+    write(*,*)""
+    ! create the aeropt_mod log file for this aerosol
+    write(aeroptlogfile,*) "aeropt_log_",trim(name_aer(iaer)),".txt"
+    aeroptlogfileID = 100+iaer ! update the log file unit
+    CALL create_logfile(aeroptlogfile)
+
+    write(*,*)"Reading the optical properties..."
+    ! Read the properties tables
     CALL read_optpropfile(datadir,optpropfile_aer(iaer))
 
@@ -401 +439 @@
            ! out of the optpropfile boundaries
            opa_aer(ilon,ilat,ialt,it)=missval(iaer)
-           write(*,*) "(mmr = ",mmr(ilon,ilat,ialt,it),"kg/kg)"
-           write(*,*)""
+           write(aeroptlogfileID,*) "(mmr = ",mmr(ilon,ilat,ialt,it),"kg/kg)"
+           write(aeroptlogfileID,*)""
            
          else if ((mmr(ilon,ilat,ialt,it).eq.missval(iaer))) then
@@ -419 +457 @@
 
          endif
+         
+         if (trim(compute_col).eq."yes") then
+           ! COMPUTATION OF THE COLUMN-INTEGRATED OPTICAL DEPTH [/]
+           if (opa_aer(ilon,ilat,ialt,it).ne.missval(iaer)) then
+             tau_aer(ilon,ilat,it)= tau_aer(ilon,ilat,it)+opa_aer(ilon,ilat,ialt,it)*delta_z(ilon,ilat,ialt,it)
+           endif
+         endif ! trim(compute_col).eq."yes"
              
        enddo ! ialt
+       
+       if (trim(compute_col).eq."yes") then
+         ! COMPUTATION OF THE COLUMN-INTEGRATED OPTICAL DEPTH [/]
+         if (tau_aer(ilon,ilat,it).eq.0) then
+           tau_aer(ilon,ilat,it)= missval(iaer)
+         endif
+       endif ! trim(compute_col).eq."yes"
+       
       enddo ! it
      enddo ! ilat
@@ -426 +479 @@
 
 !==========================================================================
-! 3.3 Writing in the output file
+! 3.3 Writing opacity in the output file
 !==========================================================================
     ierr = NF90_PUT_VAR(outfid, tmpvaridout, opa_aer(:,:,:,:))
@@ -433 +486 @@
     
 !==========================================================================
-! 3.4 Prepare next loop
+! 3.4 Writing column-integrated optical depth in the output file
+!==========================================================================
+    if (trim(compute_col).eq."yes") then
+      ! 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(*,*)""
+
+      ! Insert the definition of the variable
+      tmpvarname="tau_"//trim(name_aer(iaer))
+      ierr=NF90_DEF_VAR(outfid,trim(tmpvarname),nf90_float,tauvarshape,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
+        tmplongname=trim(name_aer(iaer))//" extinction column-integrated optical depth at ref. wavelength"
+        ierr=nf90_put_att(outfid,tmpvaridout,"long_name",tmplongname)
+      else ! opatype.eq."abs"
+        tmplongname=trim(name_aer(iaer))//" absorption column-integrated optical depth at ref. wavelength"
+        ierr=nf90_put_att(outfid,tmpvaridout,"long_name",tmplongname)
+      endif
+
+      ierr=nf90_put_att(outfid,tmpvaridout,"units","/")
+      ierr=nf90_put_att(outfid,tmpvaridout,"refwavelength",wvl_val)
+      ierr=nf90_put_att(outfid,tmpvaridout,"missing_value",missval(iaer))
+      write(*,*)"with missing value = ",missval(iaer)
+
+      ! End netcdf define mode
+      ierr=nf90_enddef(outfid)
+      error_text="ERROR: Couldn't end netcdf define mode"
+      call status_check(ierr,error_text)
+
+
+      ! Writing variable in output file
+      ierr = NF90_PUT_VAR(outfid, tmpvaridout, tau_aer(:,:,:))
+      error_text="Error: could not write "//trim(tmpvarname)//" data in the outfile"
+      call status_check(ierr,error_text)
+    endif ! trim(compute_col).eq."yes"
+    
+    
+!==========================================================================
+! 3.5 Prepare next loop
 !==========================================================================
     DEALLOCATE(mmr)
     DEALLOCATE(reff)
     DEALLOCATE(opa_aer)
+    IF (allocated(tau_aer)) DEALLOCATE(tau_aer)
+    
     CALL end_aeropt_mod
     
@@ -837 +936 @@
 ierr=nf90_inq_dimid(gcmfid,"GCM_layers",tmpdimid)
 if (ierr.ne.nf90_noerr) then
-  ! didn't find a GCM_layers dimension; therefore we have:
-  GCM_layers=altlen  
+  ! didn't find a GCM_layers dimension;
+  ! we look for interlayer dimension
+  ierr=nf90_inq_dimid(gcmfid,"interlayer",tmpdimid)
+  if (ierr.eq.nf90_noerr) then
+    ! load value of GCM_layers
+    ierr=nf90_inquire_dimension(gcmfid,tmpdimid,len=GCM_layers)
+    GCM_layers=GCM_layers-1
+  else
+    ! didn't find a GCM_layers or interlayer dimension; therefore we have:
+    GCM_layers=altlen
+  endif 
 else
   ! load value of GCM_layers
@@ -862 +970 @@
 !*******************************************************************************
 
-subroutine init2(gcmfid,lonlen,latlen,altlen,GCM_layers, &
-                 outfid,londimout,latdimout,altdimout, &
-                 layerdimout,interlayerdimout)
+subroutine init2(gcmfid,lonlen,latlen,altlen,timelen,GCM_layers, &
+                 outfid,londimout,latdimout,altdimout,timedimout, &
+                 layerdimout,interlayerdimout,&
+                 compute_col,delta_z)
 !==============================================================================
 ! Purpose:
 ! Copy ap() , bp(), aps(), bps(), aire() and phisinit()
-! from input file to outpout file
+! from input file to output file
+! + compute layers' height if needed
 !==============================================================================
 ! Remarks:
@@ -881 +991 @@
 ! Arguments:
 !==============================================================================
-integer, intent(in) :: gcmfid  ! NetCDF output file ID
-integer, intent(in) :: lonlen ! # of grid points along longitude
-integer, intent(in) :: latlen ! # of grid points along latitude
-integer, intent(in) :: altlen ! # of grid points along altitude
-integer, intent(in) :: GCM_layers ! # of GCM atmospheric layers
-integer, intent(in) :: outfid ! NetCDF output file ID
-integer, intent(in) :: londimout ! longitude dimension ID
-integer, intent(in) :: latdimout ! latitude dimension ID
-integer, intent(in) :: altdimout ! altitude dimension ID
-integer, intent(in) :: layerdimout ! GCM_layers dimension ID
-integer, intent(in) :: interlayerdimout ! GCM_layers+1 dimension ID
+integer, intent(in) :: gcmfid               ! NetCDF output file ID
+integer, intent(in) :: lonlen               ! # of grid points along longitude
+integer, intent(in) :: latlen               ! # of grid points along latitude
+integer, intent(in) :: altlen               ! # of grid points along altitude
+integer, intent(in) :: timelen              ! # of grid points along time
+integer, intent(in) :: GCM_layers           ! # of GCM atmospheric layers
+integer, intent(in) :: outfid               ! NetCDF output file ID
+integer, intent(in) :: londimout            ! longitude dimension ID
+integer, intent(in) :: latdimout            ! latitude dimension ID
+integer, intent(in) :: altdimout            ! altitude dimension ID
+integer, intent(in) :: timedimout           ! time dimension ID
+integer, intent(in) :: layerdimout          ! GCM_layers dimension ID
+integer, intent(in) :: interlayerdimout     ! GCM_layers+1 dimension ID
+character(len=3), intent(in) :: compute_col ! does the user want to compute the column-integrated optical depth [yes/no]
+
+real,dimension(lonlen,latlen,altlen,timelen), intent(out) :: delta_z ! altitude difference between interlayers [km]
+
 !==============================================================================
 ! Local variables:
 !==============================================================================
-real,dimension(:),allocatable :: aps,bps ! hybrid vertical coordinates
-real,dimension(:),allocatable :: ap,bp ! hybrid vertical coordinates
-real,dimension(:),allocatable :: sigma ! sigma levels
-real,dimension(:,:),allocatable :: aire ! mesh areas
+real,dimension(:),allocatable :: aps,bps    ! hybrid vertical coordinates
+real,dimension(:),allocatable :: ap,bp      ! hybrid vertical coordinates
+real,dimension(:),allocatable :: sigma      ! sigma levels
+real,dimension(:,:),allocatable :: aire     ! mesh areas
 real,dimension(:,:),allocatable :: phisinit ! Ground geopotential
-integer :: ierr
-integer :: tmpvarid ! temporary variable ID
-logical :: area ! is "aire" available ?
-logical :: phis ! is "phisinit" available ?
+
+integer :: ierr                  ! [netcdf] subroutine returned error code
+character(len=256) :: error_text ! text to display in case of error
+integer :: tmpvarid              ! temporary variable ID
+
+logical :: area   ! is "aire" available ?
+logical :: phis   ! is "phisinit" available ?
 logical :: hybrid ! are "aps" and "bps" available ?
-logical :: apbp ! are "ap" and "bp" available ?
+logical :: apbp   ! are "ap" and "bp" available ? + are they usable for delta_z computation ?
+
+
+! for delta_z computation
+real,dimension(:,:,:,:),allocatable :: zlev ! altitude of the interlayers [km]
+logical :: is_zlev                          ! is zzlev available ?
+integer,dimension(4) :: zlevdimids          ! dimensions of zzlev in GCM file
+integer :: zlevvertlen                      ! length of the vertical dimension of zlev
+
+real,dimension(altlen) :: alt               ! altitude coordinates
+character (len=64) :: altlong_name,altunits ! altitude long_name and units attributes
+
+real,dimension(:,:,:,:),allocatable :: plev ! pressure of the interlayers [km]
+real,dimension(:,:,:),allocatable :: ps     ! surface pressure [Pa]
+real,dimension(:,:,:,:),allocatable :: rho  ! atmospheric density [kg/m3]
+real,parameter :: g0 =3.7257964             ! Mars gravity [m.s-2]
+
+real,dimension(:,:,:,:),allocatable :: zlay ! altitude of the layers [km]
+
+integer :: ilon,ilat,it,ialt                ! loop indices
+real :: missval                             ! GCM variables missing value attribute
 
 !==============================================================================
@@ -1046 +1185 @@
 !==============================================================================
 
-!==============================================================================
-! 2.1 Hybrid coordinates ap() , bp(), aps() and bps()
-!==============================================================================
+!=========================================================================
+! 2.1 Hybrid coordinates ap(), bp(), aps() and bps()
+!=========================================================================
 if(hybrid) then 
 ! define aps
@@ -1166 +1305 @@
 endif ! of if (hybrid)
 
-!==============================================================================
+!=========================================================================
 ! 2.2 aire() and phisinit()
-!==============================================================================
+!=========================================================================
 
 if (area) then
@@ -1220 +1359 @@
 
 
-! Cleanup
+!==============================================================================
+! 3. Compute delta_z for column integration
+!==============================================================================
+if (trim(compute_col).eq."yes") then
+  write(*,*) "Computing the layers' height delta_z..."
+  
+!=========================================================================
+! 3.1 CASE 1: the GCM file contains zzlev variable
+!=========================================================================
+  ! Does the field exist in GCM file?
+  is_zlev=.false.
+  ierr=nf90_inq_varid(gcmfid,"zzlev",tmpvarid)
+  
+  IF (ierr.eq.nf90_noerr) THEN
+    ! zzlev is present in the GCM file
+    is_zlev=.true.
+    write(*,*) "zzlev was found in GCM file"
+    
+    ! Get the field's dimensions
+    ierr=nf90_inquire_variable(gcmfid,tmpvarid,dimids=zlevdimids)
+    
+    ! Get the length of the vertical dimension (assuming vertical dim is #3 like other variables)
+    ierr=nf90_inquire_dimension(gcmfid,zlevdimids(3),len=zlevvertlen)
+    
+    ! Get missing value
+    ierr=nf90_get_att(gcmfid,tmpvarid,"missing_value",missval)
+    if (ierr.ne.nf90_noerr) then
+      missval = 1.e+20
+    endif
+    
+    if (zlevvertlen.eq.altlen+1) then
+      ! we have every interlayer altitudes, even the top one
+      allocate(zlev(lonlen,latlen,altlen+1,timelen))
+
+      ! Get zlev variable
+      ierr=NF90_GET_VAR(gcmfid,tmpvarid,zlev)
+      error_text="Failed to load zzlev"
+      call status_check(ierr,error_text)
+      
+      ! Compute delta_z [km]
+      ! NB : the unit for zzlev variable is "m"
+      do ilon=1,lonlen
+       do ilat=1,latlen
+        do it=1,timelen
+         do ialt=1,altlen
+         
+          if ((zlev(ilon,ilat,ialt+1,it).eq.missval).or.(zlev(ilon,ilat,ialt,it).eq.missval)) then
+            delta_z(ilon,ilat,ialt,it) = 0
+          else
+            delta_z(ilon,ilat,ialt,it) = (zlev(ilon,ilat,ialt+1,it)-zlev(ilon,ilat,ialt,it)) /1000
+          endif ! missval
+           
+         enddo ! ialt=1,altlen
+        enddo ! it=1,timelen
+       enddo ! ilat=1,latlen
+      enddo ! ilon=1,lonlen
+      
+    else if (zlevvertlen.eq.altlen) then
+      ! we have every interlayer altitudes, except the top one
+      allocate(zlev(lonlen,latlen,altlen,timelen))
+
+      ! Get zlev variable
+      ierr=NF90_GET_VAR(gcmfid,tmpvarid,zlev)
+      error_text="Failed to load zzlev"
+      call status_check(ierr,error_text)
+      
+      ! Compute delta_z [km]
+      ! NB : the unit for zrecasted altitude is "m"
+      do ilon=1,lonlen
+       do ilat=1,latlen
+        do it=1,timelen
+         do ialt=1,altlen-1
+         
+          if ((zlev(ilon,ilat,ialt+1,it).eq.missval).or.(zlev(ilon,ilat,ialt,it).eq.missval)) then
+            delta_z(ilon,ilat,ialt,it) = 0
+          else
+            delta_z(ilon,ilat,ialt,it) = (zlev(ilon,ilat,ialt+1,it)-zlev(ilon,ilat,ialt,it)) /1000
+          endif ! missval
+           
+         enddo ! ialt=1,altlen-1
+          
+         ! top of last layer : we assume the same height than the layer below
+         delta_z(ilon,ilat,altlen,it) = delta_z(ilon,ilat,altlen-1,it)
+          
+        enddo ! it=1,timelen
+       enddo ! ilat=1,latlen
+      enddo ! ilon=1,lonlen
+            
+      
+      write(*,*) "delta_z has been well computed from variable zzlev"
+      write(*,*)""
+    
+    else
+      ! weird case where the GCM file would have been zrecasted but not the zzlev variable?
+      ! anyway, we can't use zzlev then
+      write(*,*) "zzlev and altitude of GCM file don't have the same dimension."
+      write(*,*) "We will try to use the altitude dimension to compute the layers' height."
+      write(*,*)""
+      is_zlev=.false.
+
+    endif ! (zlevvertlen.eq.altlen+1)
+      
+  ENDIF ! (ierr.eq.nf90_noerr) = zzlev present in file
+
+!=========================================================================
+! 3.2 CASE 2: try to use the GCM altitude dimension
+!=========================================================================
+  IF (is_zlev.eq..false.) THEN
+    ! Get the field in GCM file
+    ierr=nf90_inq_varid(gcmfid,"altitude",tmpvarid)
+
+    ierr=NF90_GET_VAR(gcmfid,tmpvarid,alt)
+    error_text="Failed to load altitude"
+    call status_check(ierr,error_text)
+
+    ! Get the field attributes in the GCM file
+    ierr=nf90_get_att(gcmfid,tmpvarid,"long_name",altlong_name)
+    if (ierr.ne.nf90_noerr) then
+    ! if no attribute "long_name", try "title"
+      ierr=nf90_get_att(gcmfid,tmpvarid,"title",altlong_name)
+    endif
+    ierr=nf90_get_att(gcmfid,tmpvarid,"units",altunits)
+    
+    IF (trim(altlong_name).eq."pseudo-alt") THEN
+      ! GCM file has not been zrecasted, so we can use ap,bp
+      ! and compute delta_z from the layer mass
+      ! we need ps and rho to use this method
+      
+      ! Check that ap,bp have the same dimension length than altlen+1
+      if (GCM_layers+1.ne.altlen+1) then
+        apbp=.false.
+      endif  
+      
+      ! Load Psurf to reconstruct pressure
+      ! Does the field exist in GCM file?
+      ierr=nf90_inq_varid(gcmfid,"ps",tmpvarid)
+      if (ierr.eq.nf90_noerr) then
+        allocate(ps(lonlen,latlen,timelen))
+        ! Get ps variable
+        ierr=NF90_GET_VAR(gcmfid,tmpvarid,ps)
+        error_text="Failed to load surface pressure ps"
+        call status_check(ierr,error_text)
+      else
+        apbp=.false.
+      endif
+
+      ! Load rho to compute delta_z from the layer mass
+      ! Does the field exist in GCM file?
+      ierr=nf90_inq_varid(gcmfid,"rho",tmpvarid)
+      if (ierr.eq.nf90_noerr) then
+        allocate(rho(lonlen,latlen,altlen,timelen))
+        ! Get temp variable
+        ierr=NF90_GET_VAR(gcmfid,tmpvarid,rho)
+        error_text="Failed to load atmospheric rho"
+        call status_check(ierr,error_text)
+      else
+        apbp=.false.
+      endif
+        
+      if (apbp) then
+        allocate(plev(lonlen,latlen,altlen+1,timelen)) 
+        ! Compute delta_z [km]
+        do ilon=1,lonlen
+         do ilat=1,latlen
+          do it=1,timelen
+           ! plev at first level
+           plev(ilon,ilat,1,it) = ap(1)+bp(1)*ps(ilon,ilat,it)
+           
+           do ialt=1,altlen
+              ! plev just above
+              plev(ilon,ilat,ialt+1,it) = ap(ialt+1)+bp(ialt+1)*ps(ilon,ilat,it)
+              ! delta_z [km] from layer mass
+              delta_z(ilon,ilat,ialt,it) = (plev(ilon,ilat,ialt,it)-plev(ilon,ilat,ialt+1,it)) &
+                                           / (g0*rho(ilon,ilat,ialt,it)) /1000
+
+           enddo ! ialt=1,altlen
+           
+          enddo ! it=1,timelen
+         enddo ! ilat=1,latlen
+        enddo ! ilon=1,lonlen
+
+        write(*,*) "delta_z has been well computed from variables ap,bp,ps,rho"
+        write(*,*)""        
+        
+      else
+        ! no ap,bp in GCM file
+        ! => take the middle of the layer altitudes as the interlayer altitude
+        ! NB : the unit for "pseudo-alt" altitude (not zrecasted) is "km"
+        do ialt=2,altlen-1
+          delta_z(:,:,ialt,:) = (alt(ialt+1)-alt(ialt-1))/2
+        enddo
+        delta_z(:,:,1,:) = alt(2)/2 ! bottom at 0km (ok since pseudo-alt can't be negative)
+        delta_z(:,:,altlen,:) = delta_z(:,:,altlen-1,:) ! top of last layer : we assume the same height than the layer below
+
+        write(*,*) "delta_z has been well computed from dimension pseudo-altitude"
+        write(*,*)""
+        
+      endif ! apbp
+      
+    ELSE IF (trim(altunits).eq."m") THEN
+      ! GCM file has been zrecasted in altitude above surface/areoid/center of planet
+      ! NB : the unit for zrecasted altitude is "m"
+      
+      ! Take the arithmetic mean of the layer altitudes as the interlayer altitude
+      do ialt=2,altlen-1
+        delta_z(:,:,ialt,:) = (alt(ialt+1)-alt(ialt-1))/2 /1000
+      enddo
+      
+      ! Bottom layer : we can't assume a bottom of the first layer at 0km
+      ! since "altitude above areoid" can be negative
+      ! We thus assume zlev(2)-zlev(1) = zlay(2)-zlay(1)
+      delta_z(:,:,1,:) = (alt(2)-alt(1)) /1000
+      
+      ! Top of last layer : we assume the same height than the layer below
+      delta_z(:,:,altlen,:) = delta_z(:,:,altlen-1,:)
+
+      write(*,*) "delta_z has been well computed from dimension altitude"
+      write(*,*)""
+    
+    ELSE IF (trim(altunits).eq."Pa") THEN
+      ! GCM file has been zrecasted in pressure
+      
+      ! Check if "zareoid" exist in GCM file
+      ierr=nf90_inq_varid(gcmfid,"zareoid",tmpvarid)
+      if (ierr.eq.nf90_noerr) then
+        allocate(zlay(lonlen,latlen,altlen,timelen))
+        ! Get zareoid variable
+        ierr=NF90_GET_VAR(gcmfid,tmpvarid,zlay)
+        error_text="Failed to load zareoid"
+        call status_check(ierr,error_text)
+        
+        ! Get missing value
+        ierr=nf90_get_att(gcmfid,tmpvarid,"missing_value",missval)
+        if (ierr.ne.nf90_noerr) then
+          missval = 1.e+20
+        endif
+        
+        ! Take the arithmetic mean of the layer altitudes as the interlayer altitude
+        ! NB : the unit for zareoid is "m"
+        ! Compute delta_z [km]
+        do ilon=1,lonlen
+         do ilat=1,latlen
+          do it=1,timelen
+          
+           ! Bottom layer
+           if ((zlay(ilon,ilat,2,it).eq.missval).or.(zlay(ilon,ilat,1,it).eq.missval)) then
+             delta_z(ilon,ilat,1,it) = 0
+           else
+             ! We can't assume a bottom of the first layer at 0km
+             ! since "altitude above areoid" can be negative
+             ! We thus assume zlev(2)-zlev(1) = zlay(2)-zlay(1)
+             delta_z(:,:,1,:) = (zlay(ilon,ilat,2,it)-zlay(ilon,ilat,1,it)) /1000
+           endif ! missval
+           
+           ! rest of the column
+           do ialt=2,altlen-1
+            if ((zlay(ilon,ilat,ialt+1,it).eq.missval).or.(zlay(ilon,ilat,ialt-1,it).eq.missval)) then
+              delta_z(ilon,ilat,ialt,it) = 0
+            else
+              delta_z(ilon,ilat,ialt,it) = (zlay(ilon,ilat,ialt+1,it)-zlay(ilon,ilat,ialt-1,it))/2 /1000
+            endif ! missval
+           enddo ! ialt=2,altlen-1
+           
+           ! Top of last layer : we assume the same height than the layer below
+           delta_z(ilon,ilat,altlen,it) = delta_z(ilon,ilat,altlen-1,it)
+           
+          enddo ! it=1,timelen
+         enddo ! ilat=1,latlen
+        enddo ! ilon=1,lonlen
+
+        write(*,*) "delta_z has been well computed from variable zareoid"
+        write(*,*)""
+        
+      else
+      
+        write(*,*) "No variable zareoid present in GCM file"
+        write(*,*) "Computing the layers' height with the layers' pressure coordinate"
+        write(*,*) "is not (yet?) handled by the program."
+        write(*,*) "The column-integrated optical depths can thus not be computed."
+        write(*,*) "Better stop now..."
+        write(*,*) "You may retry with putting 'no' for the column computation."
+        stop
+        
+      endif ! (ierr.eq.nf90_noerr) = zareoid present in file
+      
+
+    ELSE
+      ! abort the program
+      write(*,*) "delta_z can not be computed from the GCM file's variables."
+      write(*,*) "The column-integrated optical depths can thus not be computed."
+      write(*,*) "Better stop now..."
+      write(*,*) "You may retry with putting 'no' for the column computation."
+      stop
+
+    ENDIF ! trim(altlong_name).eq."pseudo-alt"
+    
+  ENDIF ! is_zlev.eq.false
+  
+!=========================================================================
+! 3.3 Write delta_z in output file
+!=========================================================================
+  
+  ! Switch to netcdf define mode
+  ierr=nf90_redef(outfid)
+  ! Insert the definition of the variable
+  ierr=NF90_DEF_VAR(outfid,"delta_z",nf90_float,(/londimout,latdimout,altdimout,timedimout/),tmpvarid)
+  if (ierr.ne.nf90_noerr) then
+    write(*,*) "init2 Error: Failed to define the variable delta_z"
+    stop
+  endif
+  ! Write the attributes
+  ierr=nf90_put_att(outfid,tmpvarid,"long_name","Layer height used for column optical depth computation")
+  ierr=nf90_put_att(outfid,tmpvarid,"units","km")
+  ! End netcdf define mode
+  ierr=nf90_enddef(outfid)
+
+  ! write delta_z
+  ierr=NF90_PUT_VAR(outfid,tmpvarid,delta_z)
+  if (ierr.ne.nf90_noerr) then
+    write(*,*) "init2 Error: Failed to write delta_z"
+    stop
+  endif
+  
+endif ! trim(compute_col).eq."yes"
+
+!==============================================================================
+! 4. Cleanup
+!==============================================================================
+
 if (allocated(aps)) deallocate(aps)
 if (allocated(bps)) deallocate(bps)
@@ -1230 +1697 @@
 
 end subroutine init2
+
+
+!*******************************************************************************

trunk/LMDZ.MARS/util/aeroptical.def

@@ -3 +3 @@
 5.e-7
 ext
+yes
 dust,dustq,reffdust,2500,optprop_dustvis_TM_n50.dat
 h2o_ice,h2o_ice,5.e-6,920,optprop_icevis_n30.dat
@@ -14 +15 @@
 3) The wavelength at which the output opacities are calibrated (in meters)
 4) Opacity type : extinction (ext) or absorption (abs)
-5) aerosol_name,aerosol_mmr(variable name),aerosol_reff(variable name or single value in m),aerosol_rho(value in kg/m3),optpropfile_name
+5) Computation of the column-integrated optical depth? (yes/no)
+
 6) aerosol_name,aerosol_mmr(variable name),aerosol_reff(variable name or single value in m),aerosol_rho(value in kg/m3),optpropfile_name
+7) aerosol_name,aerosol_mmr(variable name),aerosol_reff(variable name or single value in m),aerosol_rho(value in kg/m3),optpropfile_name
 ...
 N) aerosol_name,aerosol_mmr(variable name),aerosol_reff(variable name or single value in m),aerosol_rho(value in kg/m3),optpropfile_name