Changeset 2830 for trunk/LMDZ.MARS/util


Ignore:
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/util
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • trunk/LMDZ.MARS/util/README

    r2817 r2830  
    192192--------------------------------------------------------------------
    193193
    194 Program to compute opacities [1/km] of aerosols at a wavelength given by the user.
     194Program to compute opacity profiles opa_[aer] (1/km) of aerosols
     195at a wavelength given by the user.
    195196Computation is made from the aerosol mass mixing ratios and effective radius,
    196197associated with air density (rho) and files containing the aerosol
    197198optical properties (generally present in the GCM datadir/ directory).
    198199The user have to precise the type of opacity : extinction or absorption.
     200The user can also ask to compute the column-integrated optical depth tau_[aer].
    199201
    200202NB : this program requires to compile the module aeropt_mod.F90
  • trunk/LMDZ.MARS/util/aeropt_mod.F90

    r2817 r2830  
    1212 real,dimension(:,:),save,allocatable :: Qext    ! Extinction efficiency coefficient [/]
    1313 real,dimension(:,:),save,allocatable :: omeg    ! Single scattering albedo [/]
     14 
     15 integer,save :: aeroptlogfileID = 100           ! Log file's unit ID
    1416
    1517 CONTAINS
     
    4648  integer :: read_ok                         ! to know if the line is well read
    4749  integer :: iwvl,isize                      ! Wavelength and Particle size loop indices
     50 
     51  !==============================================================================
     52  ! aeropt_mod module variables used here:
     53  !==============================================================================
     54  ! nwvl,nsize,wvl,radiusdyn,Qext,omeg
    4855
    4956  !==========================================================================
     
    228235  real,dimension(2) :: reffQext         ! Qext after reff interpolation
    229236  real,dimension(2) :: reffomeg         ! omeg after reff interpolation
     237 
     238  !==============================================================================
     239  ! aeropt_mod module variables used here:
     240  !==============================================================================
     241  ! nwvl,nsize,wvl,radiusdyn,Qext,omeg,aeroptlogfileID
    230242
    231243  !==========================================================================
     
    262274  ENDDO
    263275
    264   if ((radiusdyn(nsize).lt.reff_val).or.(radiusdyn(1).gt.reff_val)) then
    265     write(*,*) "WARNING: the effective radius (",reff_val,"m) is out of the bounds"
    266     write(*,*) "of the optical properties' tables"
    267     write(*,*) "(reff_inf=",radiusdyn(1),"m ; reff_sup=",radiusdyn(nsize),"m)"
    268     write(*,*) "A missing value will be written for opacity."
    269 
    270     ! NB: this should especially handle cases when reff=0
     276  if (reff_val.eq.missval) then
     277    ! if the effective radius is a missing value,
     278    ! put a missing value for the optical properties
     279    Qext_val = missval
     280    omeg_val = missval
     281!    write(*,*) "reff missing value detected"
     282   
     283  else if ((radiusdyn(nsize).lt.reff_val).or.(radiusdyn(1).gt.reff_val)) then
     284    write(aeroptlogfileID,*) "WARNING: the effective radius (",reff_val,"m) is out of the bounds"
     285    write(aeroptlogfileID,*) "of the optical properties' tables"
     286    write(aeroptlogfileID,*) "(reff_inf=",radiusdyn(1),"m ; reff_sup=",radiusdyn(nsize),"m)"
     287    write(aeroptlogfileID,*) "A missing value will be written for optical properties."
     288
    271289    Qext_val = missval
    272290    omeg_val = missval
     
    333351
    334352  implicit none
     353 
     354  !==============================================================================
     355  ! aeropt_mod module variables used here:
     356  !==============================================================================
     357  ! wvl,radiusdyn,Qext,omeg
    335358
    336359  if (allocated(wvl)) deallocate(wvl)
     
    342365 
    343366!*******************************************************************************
     367 SUBROUTINE create_logfile(filename)
     368  !==============================================================================
     369  ! Purpose:
     370  ! Deallocate module variables
     371  !==============================================================================
     372
     373  implicit none
     374 
     375  !==============================================================================
     376  ! Arguments:
     377  !==============================================================================
     378  character(len=128),intent(in) :: filename ! name of the aeropt log file
     379  !==============================================================================
     380  ! aeropt_mod module variables used here:
     381  !==============================================================================
     382  ! aeroptlogfileID
     383 
     384  ! Create output file
     385  write(*,*) "Creating aeropt log file: ",trim(filename)
     386  write(*,*) "It will overwrite any existing file with the same name."
     387  OPEN(UNIT=aeroptlogfileID,FILE=trim(filename),FORM='formatted',STATUS='replace')
     388  ! NB: STATUS='replace' will overwrite any existing file with the same name
     389 
     390 END SUBROUTINE create_logfile
     391 
     392!*******************************************************************************
    344393
    345394END MODULE aeropt_mod
  • trunk/LMDZ.MARS/util/aeroptical.F90

    r2817 r2830  
    4141real :: wvl_val                                                      ! reference wavelength of the output opacities [m]
    4242character(len=3) :: opatype                                          ! opacity type : extinction (ext) or absorption (abs)
     43character(len=3) :: compute_col                                      ! does the user want to compute the column-integrated optical depth [yes/no]
    4344character(len=50),dimension(20) :: name_aer,mmrname_aer,reffname_aer ! aerosols' name, MMR varname, reff varname (or value)
    4445integer,dimension(20) :: reffval_ok                                  ! to know if the user gave a variable name or a value for reff_aer
     
    5051
    5152! GCM file
    52 character(len=128) :: gcmfile                 ! name of the netcdf GCM input file
    53 integer :: gcmfid                             ! [netcdf] GCM input file ID number
    54 integer :: ierr                               ! [netcdf] subroutine returned error code
    55 character(len=256) :: error_text              ! text to display in case of error
    56 integer :: lonlen,latlen,altlen,timelen       ! nb of grid points along longitude,latitude,altitude,time
    57 integer :: GCM_layers                         ! number of GCM layers
    58 integer :: layerdimout,interlayerdimout       ! "GCM_layers" and "GCM_layers+1" IDs
    59 
    60 logical,dimension(:),allocatable :: aerok     ! to know if the needed fields are in the GCM file
    61 integer,dimension(:),allocatable :: mmrvarid  ! stores a MMR variable ID number
    62 integer,dimension(:),allocatable :: reffvarid ! stores a reff variable ID number
    63 integer :: tmpvarid                           ! temporarily stores a variable ID number
    64 real,dimension(:,:,:,:),allocatable :: mmr    ! aerosols mass mixing ratios [kg/kg]
    65 real,dimension(:,:,:,:),allocatable :: reff   ! aerosols effective radii [m]
    66 real,dimension(:,:,:,:),allocatable :: rho    ! atmospheric density [kg.m-3]
    67 integer :: iaer                               ! Loop index for aerosols
    68 integer :: ilon,ilat,ialt,it                  ! Loop indices for coordinates
     53character(len=128) :: gcmfile                  ! name of the netcdf GCM input file
     54integer :: gcmfid                              ! [netcdf] GCM input file ID number
     55integer :: ierr                                ! [netcdf] subroutine returned error code
     56character(len=256) :: error_text               ! text to display in case of error
     57integer :: lonlen,latlen,altlen,timelen        ! nb of grid points along longitude,latitude,altitude,time
     58integer :: GCM_layers                          ! number of GCM layers
     59integer :: layerdimout,interlayerdimout        ! "GCM_layers" and "GCM_layers+1" IDs
     60
     61logical,dimension(:),allocatable :: aerok      ! to know if the needed fields are in the GCM file
     62integer,dimension(:),allocatable :: mmrvarid   ! stores a MMR variable ID number
     63integer,dimension(:),allocatable :: reffvarid  ! stores a reff variable ID number
     64integer :: tmpvarid                            ! temporarily stores a variable ID number
     65real,dimension(:,:,:,:),allocatable :: mmr     ! aerosols mass mixing ratios [kg/kg]
     66real,dimension(:,:,:,:),allocatable :: reff    ! aerosols effective radii [m]
     67real,dimension(:,:,:,:),allocatable :: rho     ! atmospheric density [kg.m-3]
     68integer :: iaer                                ! Loop index for aerosols
     69integer :: ilon,ilat,ialt,it                   ! Loop indices for coordinates
    6970
    7071! Opacity computation
     72character(len=128) :: aeroptlogfile            ! name of the aeropt_mod log file
    7173real :: Qext_val                               ! Qext after interpolations
    7274real :: omeg_val                               ! omeg after interpolations
     
    8486character(len=16) :: tmpvarname                ! temporarily stores a variable name
    8587character(len=128) :: tmplongname              ! temporarily stores a variable long_name attribute
    86 real,dimension(:,:,:,:),allocatable :: opa_aer ! Aerosols opacities [1/km]
     88real,dimension(:,:,:,:),allocatable :: opa_aer ! Aerosols opacities dtau/dz [1/km]
    8789real,dimension(:),allocatable :: missval       ! Value to put in outfile when the reff is out of bounds
    8890                                               ! or when there is a mising value in input file
     91integer :: tauvarshape(3)                      ! stores the column variable shape (of dimensions' IDs)
     92real,dimension(:,:,:,:),allocatable :: delta_z ! Layers' height for column computation [km]
     93real,dimension(:,:,:),allocatable :: tau_aer   ! Aerosols column-integrated optical depth [/]
    8994
    9095
     
    96101!3) The wavelength at which the output opacities are calibrated (in meters)
    97102!4) Opacity type : extinction (ext) or absorption (abs)
    98 !5) aerosol_name,aerosol_mmr(variable name),aerosol_reff(variable name or single value in m),aerosol_rho(value in kg/m3),optpropfile_name
     103!5) Does the user want to compute the column-integrated optical depth? [yes/no]
    99104!6) aerosol_name,aerosol_mmr(variable name),aerosol_reff(variable name or single value in m),aerosol_rho(value in kg/m3),optpropfile_name
     105!7) aerosol_name,aerosol_mmr(variable name),aerosol_reff(variable name or single value in m),aerosol_rho(value in kg/m3),optpropfile_name
    100106!...
    101107!N) aerosol_name,aerosol_mmr(variable name),aerosol_reff(variable name or single value in m),aerosol_rho(value in kg/m3),optpropfile_name
     
    129135endif
    130136
    131 
    132 
    133 ! 5-N) Ask for aerosols to compute to the user
     137! 5) Computation of the column-integrated optical depth?
     138write(*,*)""
     139WRITE(*,*) "Do you want to compute the column-integrated optical depth ? (yes/no)"
     140READ(*,*) compute_col
     141if ((trim(compute_col).ne."yes").and.(trim(compute_col).ne."no")) then
     142  write(*,*) "compute_col = "//compute_col
     143  write(*,*) "Error: please answer yes or no"
     144  stop
     145endif
     146
     147
     148! 6-N) Ask for aerosols to compute to the user
    134149write(*,*)""
    135150write(*,*) "For which aerosols do you want to compute the opacity?"
     
    211226             GCM_layers,layerdimout,interlayerdimout)
    212227
     228! Length allocation for each dimension of the 4D variable delta_z
     229allocate(delta_z(lonlen,latlen,altlen,timelen))
     230
    213231! Initialize output file's aps,bps,ap,bp and phisinit variables
    214 call init2(gcmfid,lonlen,latlen,altlen,GCM_layers,&
    215            outfid,londimout,latdimout,altdimout,&
    216            layerdimout,interlayerdimout)
     232! + compute delta_z for the column integration
     233call init2(gcmfid,lonlen,latlen,altlen,timelen,GCM_layers,&
     234           outfid,londimout,latdimout,altdimout,timedimout,&
     235           layerdimout,interlayerdimout,&
     236           compute_col,delta_z)
    217237
    218238!==========================================================================
     
    297317varshape(4)=timedimout
    298318
     319tauvarshape(1)=londimout
     320tauvarshape(2)=latdimout
     321tauvarshape(3)=timedimout
     322
    299323
    300324! Loop on all aerosols
     
    333357   
    334358 
    335     ! Length allocation of the output variable
     359    ! Length allocation of the output variables
    336360    ALLOCATE(opa_aer(lonlen,latlen,altlen,timelen))
    337361   
     362    if (trim(compute_col).eq."yes") then
     363      ALLOCATE(tau_aer(lonlen,latlen,timelen))
     364      tau_aer(:,:,:) = 0.
     365    endif ! trim(compute_col).eq."yes"
     366   
     367   
    338368!==========================================================================
    339369! 2.2 Reading optical properties
    340370!==========================================================================
    341371
     372    write(*,*)""
     373    ! create the aeropt_mod log file for this aerosol
     374    write(aeroptlogfile,*) "aeropt_log_",trim(name_aer(iaer)),".txt"
     375    aeroptlogfileID = 100+iaer ! update the log file unit
     376    CALL create_logfile(aeroptlogfile)
     377
     378    write(*,*)"Reading the optical properties..."
     379    ! Read the properties tables
    342380    CALL read_optpropfile(datadir,optpropfile_aer(iaer))
    343381
     
    401439           ! out of the optpropfile boundaries
    402440           opa_aer(ilon,ilat,ialt,it)=missval(iaer)
    403            write(*,*) "(mmr = ",mmr(ilon,ilat,ialt,it),"kg/kg)"
    404            write(*,*)""
     441           write(aeroptlogfileID,*) "(mmr = ",mmr(ilon,ilat,ialt,it),"kg/kg)"
     442           write(aeroptlogfileID,*)""
    405443           
    406444         else if ((mmr(ilon,ilat,ialt,it).eq.missval(iaer))) then
     
    419457
    420458         endif
     459         
     460         if (trim(compute_col).eq."yes") then
     461           ! COMPUTATION OF THE COLUMN-INTEGRATED OPTICAL DEPTH [/]
     462           if (opa_aer(ilon,ilat,ialt,it).ne.missval(iaer)) then
     463             tau_aer(ilon,ilat,it)= tau_aer(ilon,ilat,it)+opa_aer(ilon,ilat,ialt,it)*delta_z(ilon,ilat,ialt,it)
     464           endif
     465         endif ! trim(compute_col).eq."yes"
    421466             
    422467       enddo ! ialt
     468       
     469       if (trim(compute_col).eq."yes") then
     470         ! COMPUTATION OF THE COLUMN-INTEGRATED OPTICAL DEPTH [/]
     471         if (tau_aer(ilon,ilat,it).eq.0) then
     472           tau_aer(ilon,ilat,it)= missval(iaer)
     473         endif
     474       endif ! trim(compute_col).eq."yes"
     475       
    423476      enddo ! it
    424477     enddo ! ilat
     
    426479
    427480!==========================================================================
    428 ! 3.3 Writing in the output file
     481! 3.3 Writing opacity in the output file
    429482!==========================================================================
    430483    ierr = NF90_PUT_VAR(outfid, tmpvaridout, opa_aer(:,:,:,:))
     
    433486   
    434487!==========================================================================
    435 ! 3.4 Prepare next loop
     488! 3.4 Writing column-integrated optical depth in the output file
     489!==========================================================================
     490    if (trim(compute_col).eq."yes") then
     491      ! Switch to netcdf define mode
     492      ierr=nf90_redef(outfid)
     493      error_text="ERROR: Couldn't start netcdf define mode"
     494      call status_check(ierr,error_text)
     495      write(*,*)""
     496
     497      ! Insert the definition of the variable
     498      tmpvarname="tau_"//trim(name_aer(iaer))
     499      ierr=NF90_DEF_VAR(outfid,trim(tmpvarname),nf90_float,tauvarshape,tmpvaridout)
     500      error_text="ERROR: Couldn't create "//trim(tmpvarname)//" in the outfile"
     501      call status_check(ierr,error_text)
     502      write(*,*) trim(tmpvarname)//" has been created in the outfile"
     503
     504      ! Write the attributes
     505      if (opatype.eq."ext") then
     506        tmplongname=trim(name_aer(iaer))//" extinction column-integrated optical depth at ref. wavelength"
     507        ierr=nf90_put_att(outfid,tmpvaridout,"long_name",tmplongname)
     508      else ! opatype.eq."abs"
     509        tmplongname=trim(name_aer(iaer))//" absorption column-integrated optical depth at ref. wavelength"
     510        ierr=nf90_put_att(outfid,tmpvaridout,"long_name",tmplongname)
     511      endif
     512
     513      ierr=nf90_put_att(outfid,tmpvaridout,"units","/")
     514      ierr=nf90_put_att(outfid,tmpvaridout,"refwavelength",wvl_val)
     515      ierr=nf90_put_att(outfid,tmpvaridout,"missing_value",missval(iaer))
     516      write(*,*)"with missing value = ",missval(iaer)
     517
     518      ! End netcdf define mode
     519      ierr=nf90_enddef(outfid)
     520      error_text="ERROR: Couldn't end netcdf define mode"
     521      call status_check(ierr,error_text)
     522
     523
     524      ! Writing variable in output file
     525      ierr = NF90_PUT_VAR(outfid, tmpvaridout, tau_aer(:,:,:))
     526      error_text="Error: could not write "//trim(tmpvarname)//" data in the outfile"
     527      call status_check(ierr,error_text)
     528    endif ! trim(compute_col).eq."yes"
     529   
     530   
     531!==========================================================================
     532! 3.5 Prepare next loop
    436533!==========================================================================
    437534    DEALLOCATE(mmr)
    438535    DEALLOCATE(reff)
    439536    DEALLOCATE(opa_aer)
     537    IF (allocated(tau_aer)) DEALLOCATE(tau_aer)
     538   
    440539    CALL end_aeropt_mod
    441540   
     
    837936ierr=nf90_inq_dimid(gcmfid,"GCM_layers",tmpdimid)
    838937if (ierr.ne.nf90_noerr) then
    839   ! didn't find a GCM_layers dimension; therefore we have:
    840   GCM_layers=altlen 
     938  ! didn't find a GCM_layers dimension;
     939  ! we look for interlayer dimension
     940  ierr=nf90_inq_dimid(gcmfid,"interlayer",tmpdimid)
     941  if (ierr.eq.nf90_noerr) then
     942    ! load value of GCM_layers
     943    ierr=nf90_inquire_dimension(gcmfid,tmpdimid,len=GCM_layers)
     944    GCM_layers=GCM_layers-1
     945  else
     946    ! didn't find a GCM_layers or interlayer dimension; therefore we have:
     947    GCM_layers=altlen
     948  endif
    841949else
    842950  ! load value of GCM_layers
     
    862970!*******************************************************************************
    863971
    864 subroutine init2(gcmfid,lonlen,latlen,altlen,GCM_layers, &
    865                  outfid,londimout,latdimout,altdimout, &
    866                  layerdimout,interlayerdimout)
     972subroutine init2(gcmfid,lonlen,latlen,altlen,timelen,GCM_layers, &
     973                 outfid,londimout,latdimout,altdimout,timedimout, &
     974                 layerdimout,interlayerdimout,&
     975                 compute_col,delta_z)
    867976!==============================================================================
    868977! Purpose:
    869978! Copy ap() , bp(), aps(), bps(), aire() and phisinit()
    870 ! from input file to outpout file
     979! from input file to output file
     980! + compute layers' height if needed
    871981!==============================================================================
    872982! Remarks:
     
    881991! Arguments:
    882992!==============================================================================
    883 integer, intent(in) :: gcmfid  ! NetCDF output file ID
    884 integer, intent(in) :: lonlen ! # of grid points along longitude
    885 integer, intent(in) :: latlen ! # of grid points along latitude
    886 integer, intent(in) :: altlen ! # of grid points along altitude
    887 integer, intent(in) :: GCM_layers ! # of GCM atmospheric layers
    888 integer, intent(in) :: outfid ! NetCDF output file ID
    889 integer, intent(in) :: londimout ! longitude dimension ID
    890 integer, intent(in) :: latdimout ! latitude dimension ID
    891 integer, intent(in) :: altdimout ! altitude dimension ID
    892 integer, intent(in) :: layerdimout ! GCM_layers dimension ID
    893 integer, intent(in) :: interlayerdimout ! GCM_layers+1 dimension ID
     993integer, intent(in) :: gcmfid               ! NetCDF output file ID
     994integer, intent(in) :: lonlen               ! # of grid points along longitude
     995integer, intent(in) :: latlen               ! # of grid points along latitude
     996integer, intent(in) :: altlen               ! # of grid points along altitude
     997integer, intent(in) :: timelen              ! # of grid points along time
     998integer, intent(in) :: GCM_layers           ! # of GCM atmospheric layers
     999integer, intent(in) :: outfid               ! NetCDF output file ID
     1000integer, intent(in) :: londimout            ! longitude dimension ID
     1001integer, intent(in) :: latdimout            ! latitude dimension ID
     1002integer, intent(in) :: altdimout            ! altitude dimension ID
     1003integer, intent(in) :: timedimout           ! time dimension ID
     1004integer, intent(in) :: layerdimout          ! GCM_layers dimension ID
     1005integer, intent(in) :: interlayerdimout     ! GCM_layers+1 dimension ID
     1006character(len=3), intent(in) :: compute_col ! does the user want to compute the column-integrated optical depth [yes/no]
     1007
     1008real,dimension(lonlen,latlen,altlen,timelen), intent(out) :: delta_z ! altitude difference between interlayers [km]
     1009
    8941010!==============================================================================
    8951011! Local variables:
    8961012!==============================================================================
    897 real,dimension(:),allocatable :: aps,bps ! hybrid vertical coordinates
    898 real,dimension(:),allocatable :: ap,bp ! hybrid vertical coordinates
    899 real,dimension(:),allocatable :: sigma ! sigma levels
    900 real,dimension(:,:),allocatable :: aire ! mesh areas
     1013real,dimension(:),allocatable :: aps,bps    ! hybrid vertical coordinates
     1014real,dimension(:),allocatable :: ap,bp      ! hybrid vertical coordinates
     1015real,dimension(:),allocatable :: sigma      ! sigma levels
     1016real,dimension(:,:),allocatable :: aire     ! mesh areas
    9011017real,dimension(:,:),allocatable :: phisinit ! Ground geopotential
    902 integer :: ierr
    903 integer :: tmpvarid ! temporary variable ID
    904 logical :: area ! is "aire" available ?
    905 logical :: phis ! is "phisinit" available ?
     1018
     1019integer :: ierr                  ! [netcdf] subroutine returned error code
     1020character(len=256) :: error_text ! text to display in case of error
     1021integer :: tmpvarid              ! temporary variable ID
     1022
     1023logical :: area   ! is "aire" available ?
     1024logical :: phis   ! is "phisinit" available ?
    9061025logical :: hybrid ! are "aps" and "bps" available ?
    907 logical :: apbp ! are "ap" and "bp" available ?
     1026logical :: apbp   ! are "ap" and "bp" available ? + are they usable for delta_z computation ?
     1027
     1028
     1029! for delta_z computation
     1030real,dimension(:,:,:,:),allocatable :: zlev ! altitude of the interlayers [km]
     1031logical :: is_zlev                          ! is zzlev available ?
     1032integer,dimension(4) :: zlevdimids          ! dimensions of zzlev in GCM file
     1033integer :: zlevvertlen                      ! length of the vertical dimension of zlev
     1034
     1035real,dimension(altlen) :: alt               ! altitude coordinates
     1036character (len=64) :: altlong_name,altunits ! altitude long_name and units attributes
     1037
     1038real,dimension(:,:,:,:),allocatable :: plev ! pressure of the interlayers [km]
     1039real,dimension(:,:,:),allocatable :: ps     ! surface pressure [Pa]
     1040real,dimension(:,:,:,:),allocatable :: rho  ! atmospheric density [kg/m3]
     1041real,parameter :: g0 =3.7257964             ! Mars gravity [m.s-2]
     1042
     1043real,dimension(:,:,:,:),allocatable :: zlay ! altitude of the layers [km]
     1044
     1045integer :: ilon,ilat,it,ialt                ! loop indices
     1046real :: missval                             ! GCM variables missing value attribute
    9081047
    9091048!==============================================================================
     
    10461185!==============================================================================
    10471186
    1048 !==============================================================================
    1049 ! 2.1 Hybrid coordinates ap() , bp(), aps() and bps()
    1050 !==============================================================================
     1187!=========================================================================
     1188! 2.1 Hybrid coordinates ap(), bp(), aps() and bps()
     1189!=========================================================================
    10511190if(hybrid) then
    10521191! define aps
     
    11661305endif ! of if (hybrid)
    11671306
    1168 !==============================================================================
     1307!=========================================================================
    11691308! 2.2 aire() and phisinit()
    1170 !==============================================================================
     1309!=========================================================================
    11711310
    11721311if (area) then
     
    12201359
    12211360
    1222 ! Cleanup
     1361!==============================================================================
     1362! 3. Compute delta_z for column integration
     1363!==============================================================================
     1364if (trim(compute_col).eq."yes") then
     1365  write(*,*) "Computing the layers' height delta_z..."
     1366 
     1367!=========================================================================
     1368! 3.1 CASE 1: the GCM file contains zzlev variable
     1369!=========================================================================
     1370  ! Does the field exist in GCM file?
     1371  is_zlev=.false.
     1372  ierr=nf90_inq_varid(gcmfid,"zzlev",tmpvarid)
     1373 
     1374  IF (ierr.eq.nf90_noerr) THEN
     1375    ! zzlev is present in the GCM file
     1376    is_zlev=.true.
     1377    write(*,*) "zzlev was found in GCM file"
     1378   
     1379    ! Get the field's dimensions
     1380    ierr=nf90_inquire_variable(gcmfid,tmpvarid,dimids=zlevdimids)
     1381   
     1382    ! Get the length of the vertical dimension (assuming vertical dim is #3 like other variables)
     1383    ierr=nf90_inquire_dimension(gcmfid,zlevdimids(3),len=zlevvertlen)
     1384   
     1385    ! Get missing value
     1386    ierr=nf90_get_att(gcmfid,tmpvarid,"missing_value",missval)
     1387    if (ierr.ne.nf90_noerr) then
     1388      missval = 1.e+20
     1389    endif
     1390   
     1391    if (zlevvertlen.eq.altlen+1) then
     1392      ! we have every interlayer altitudes, even the top one
     1393      allocate(zlev(lonlen,latlen,altlen+1,timelen))
     1394
     1395      ! Get zlev variable
     1396      ierr=NF90_GET_VAR(gcmfid,tmpvarid,zlev)
     1397      error_text="Failed to load zzlev"
     1398      call status_check(ierr,error_text)
     1399     
     1400      ! Compute delta_z [km]
     1401      ! NB : the unit for zzlev variable is "m"
     1402      do ilon=1,lonlen
     1403       do ilat=1,latlen
     1404        do it=1,timelen
     1405         do ialt=1,altlen
     1406         
     1407          if ((zlev(ilon,ilat,ialt+1,it).eq.missval).or.(zlev(ilon,ilat,ialt,it).eq.missval)) then
     1408            delta_z(ilon,ilat,ialt,it) = 0
     1409          else
     1410            delta_z(ilon,ilat,ialt,it) = (zlev(ilon,ilat,ialt+1,it)-zlev(ilon,ilat,ialt,it)) /1000
     1411          endif ! missval
     1412           
     1413         enddo ! ialt=1,altlen
     1414        enddo ! it=1,timelen
     1415       enddo ! ilat=1,latlen
     1416      enddo ! ilon=1,lonlen
     1417     
     1418    else if (zlevvertlen.eq.altlen) then
     1419      ! we have every interlayer altitudes, except the top one
     1420      allocate(zlev(lonlen,latlen,altlen,timelen))
     1421
     1422      ! Get zlev variable
     1423      ierr=NF90_GET_VAR(gcmfid,tmpvarid,zlev)
     1424      error_text="Failed to load zzlev"
     1425      call status_check(ierr,error_text)
     1426     
     1427      ! Compute delta_z [km]
     1428      ! NB : the unit for zrecasted altitude is "m"
     1429      do ilon=1,lonlen
     1430       do ilat=1,latlen
     1431        do it=1,timelen
     1432         do ialt=1,altlen-1
     1433         
     1434          if ((zlev(ilon,ilat,ialt+1,it).eq.missval).or.(zlev(ilon,ilat,ialt,it).eq.missval)) then
     1435            delta_z(ilon,ilat,ialt,it) = 0
     1436          else
     1437            delta_z(ilon,ilat,ialt,it) = (zlev(ilon,ilat,ialt+1,it)-zlev(ilon,ilat,ialt,it)) /1000
     1438          endif ! missval
     1439           
     1440         enddo ! ialt=1,altlen-1
     1441         
     1442         ! top of last layer : we assume the same height than the layer below
     1443         delta_z(ilon,ilat,altlen,it) = delta_z(ilon,ilat,altlen-1,it)
     1444         
     1445        enddo ! it=1,timelen
     1446       enddo ! ilat=1,latlen
     1447      enddo ! ilon=1,lonlen
     1448           
     1449     
     1450      write(*,*) "delta_z has been well computed from variable zzlev"
     1451      write(*,*)""
     1452   
     1453    else
     1454      ! weird case where the GCM file would have been zrecasted but not the zzlev variable?
     1455      ! anyway, we can't use zzlev then
     1456      write(*,*) "zzlev and altitude of GCM file don't have the same dimension."
     1457      write(*,*) "We will try to use the altitude dimension to compute the layers' height."
     1458      write(*,*)""
     1459      is_zlev=.false.
     1460
     1461    endif ! (zlevvertlen.eq.altlen+1)
     1462     
     1463  ENDIF ! (ierr.eq.nf90_noerr) = zzlev present in file
     1464
     1465!=========================================================================
     1466! 3.2 CASE 2: try to use the GCM altitude dimension
     1467!=========================================================================
     1468  IF (is_zlev.eq..false.) THEN
     1469    ! Get the field in GCM file
     1470    ierr=nf90_inq_varid(gcmfid,"altitude",tmpvarid)
     1471
     1472    ierr=NF90_GET_VAR(gcmfid,tmpvarid,alt)
     1473    error_text="Failed to load altitude"
     1474    call status_check(ierr,error_text)
     1475
     1476    ! Get the field attributes in the GCM file
     1477    ierr=nf90_get_att(gcmfid,tmpvarid,"long_name",altlong_name)
     1478    if (ierr.ne.nf90_noerr) then
     1479    ! if no attribute "long_name", try "title"
     1480      ierr=nf90_get_att(gcmfid,tmpvarid,"title",altlong_name)
     1481    endif
     1482    ierr=nf90_get_att(gcmfid,tmpvarid,"units",altunits)
     1483   
     1484    IF (trim(altlong_name).eq."pseudo-alt") THEN
     1485      ! GCM file has not been zrecasted, so we can use ap,bp
     1486      ! and compute delta_z from the layer mass
     1487      ! we need ps and rho to use this method
     1488     
     1489      ! Check that ap,bp have the same dimension length than altlen+1
     1490      if (GCM_layers+1.ne.altlen+1) then
     1491        apbp=.false.
     1492      endif 
     1493     
     1494      ! Load Psurf to reconstruct pressure
     1495      ! Does the field exist in GCM file?
     1496      ierr=nf90_inq_varid(gcmfid,"ps",tmpvarid)
     1497      if (ierr.eq.nf90_noerr) then
     1498        allocate(ps(lonlen,latlen,timelen))
     1499        ! Get ps variable
     1500        ierr=NF90_GET_VAR(gcmfid,tmpvarid,ps)
     1501        error_text="Failed to load surface pressure ps"
     1502        call status_check(ierr,error_text)
     1503      else
     1504        apbp=.false.
     1505      endif
     1506
     1507      ! Load rho to compute delta_z from the layer mass
     1508      ! Does the field exist in GCM file?
     1509      ierr=nf90_inq_varid(gcmfid,"rho",tmpvarid)
     1510      if (ierr.eq.nf90_noerr) then
     1511        allocate(rho(lonlen,latlen,altlen,timelen))
     1512        ! Get temp variable
     1513        ierr=NF90_GET_VAR(gcmfid,tmpvarid,rho)
     1514        error_text="Failed to load atmospheric rho"
     1515        call status_check(ierr,error_text)
     1516      else
     1517        apbp=.false.
     1518      endif
     1519       
     1520      if (apbp) then
     1521        allocate(plev(lonlen,latlen,altlen+1,timelen))
     1522        ! Compute delta_z [km]
     1523        do ilon=1,lonlen
     1524         do ilat=1,latlen
     1525          do it=1,timelen
     1526           ! plev at first level
     1527           plev(ilon,ilat,1,it) = ap(1)+bp(1)*ps(ilon,ilat,it)
     1528           
     1529           do ialt=1,altlen
     1530              ! plev just above
     1531              plev(ilon,ilat,ialt+1,it) = ap(ialt+1)+bp(ialt+1)*ps(ilon,ilat,it)
     1532              ! delta_z [km] from layer mass
     1533              delta_z(ilon,ilat,ialt,it) = (plev(ilon,ilat,ialt,it)-plev(ilon,ilat,ialt+1,it)) &
     1534                                           / (g0*rho(ilon,ilat,ialt,it)) /1000
     1535
     1536           enddo ! ialt=1,altlen
     1537           
     1538          enddo ! it=1,timelen
     1539         enddo ! ilat=1,latlen
     1540        enddo ! ilon=1,lonlen
     1541
     1542        write(*,*) "delta_z has been well computed from variables ap,bp,ps,rho"
     1543        write(*,*)""       
     1544       
     1545      else
     1546        ! no ap,bp in GCM file
     1547        ! => take the middle of the layer altitudes as the interlayer altitude
     1548        ! NB : the unit for "pseudo-alt" altitude (not zrecasted) is "km"
     1549        do ialt=2,altlen-1
     1550          delta_z(:,:,ialt,:) = (alt(ialt+1)-alt(ialt-1))/2
     1551        enddo
     1552        delta_z(:,:,1,:) = alt(2)/2 ! bottom at 0km (ok since pseudo-alt can't be negative)
     1553        delta_z(:,:,altlen,:) = delta_z(:,:,altlen-1,:) ! top of last layer : we assume the same height than the layer below
     1554
     1555        write(*,*) "delta_z has been well computed from dimension pseudo-altitude"
     1556        write(*,*)""
     1557       
     1558      endif ! apbp
     1559     
     1560    ELSE IF (trim(altunits).eq."m") THEN
     1561      ! GCM file has been zrecasted in altitude above surface/areoid/center of planet
     1562      ! NB : the unit for zrecasted altitude is "m"
     1563     
     1564      ! Take the arithmetic mean of the layer altitudes as the interlayer altitude
     1565      do ialt=2,altlen-1
     1566        delta_z(:,:,ialt,:) = (alt(ialt+1)-alt(ialt-1))/2 /1000
     1567      enddo
     1568     
     1569      ! Bottom layer : we can't assume a bottom of the first layer at 0km
     1570      ! since "altitude above areoid" can be negative
     1571      ! We thus assume zlev(2)-zlev(1) = zlay(2)-zlay(1)
     1572      delta_z(:,:,1,:) = (alt(2)-alt(1)) /1000
     1573     
     1574      ! Top of last layer : we assume the same height than the layer below
     1575      delta_z(:,:,altlen,:) = delta_z(:,:,altlen-1,:)
     1576
     1577      write(*,*) "delta_z has been well computed from dimension altitude"
     1578      write(*,*)""
     1579   
     1580    ELSE IF (trim(altunits).eq."Pa") THEN
     1581      ! GCM file has been zrecasted in pressure
     1582     
     1583      ! Check if "zareoid" exist in GCM file
     1584      ierr=nf90_inq_varid(gcmfid,"zareoid",tmpvarid)
     1585      if (ierr.eq.nf90_noerr) then
     1586        allocate(zlay(lonlen,latlen,altlen,timelen))
     1587        ! Get zareoid variable
     1588        ierr=NF90_GET_VAR(gcmfid,tmpvarid,zlay)
     1589        error_text="Failed to load zareoid"
     1590        call status_check(ierr,error_text)
     1591       
     1592        ! Get missing value
     1593        ierr=nf90_get_att(gcmfid,tmpvarid,"missing_value",missval)
     1594        if (ierr.ne.nf90_noerr) then
     1595          missval = 1.e+20
     1596        endif
     1597       
     1598        ! Take the arithmetic mean of the layer altitudes as the interlayer altitude
     1599        ! NB : the unit for zareoid is "m"
     1600        ! Compute delta_z [km]
     1601        do ilon=1,lonlen
     1602         do ilat=1,latlen
     1603          do it=1,timelen
     1604         
     1605           ! Bottom layer
     1606           if ((zlay(ilon,ilat,2,it).eq.missval).or.(zlay(ilon,ilat,1,it).eq.missval)) then
     1607             delta_z(ilon,ilat,1,it) = 0
     1608           else
     1609             ! We can't assume a bottom of the first layer at 0km
     1610             ! since "altitude above areoid" can be negative
     1611             ! We thus assume zlev(2)-zlev(1) = zlay(2)-zlay(1)
     1612             delta_z(:,:,1,:) = (zlay(ilon,ilat,2,it)-zlay(ilon,ilat,1,it)) /1000
     1613           endif ! missval
     1614           
     1615           ! rest of the column
     1616           do ialt=2,altlen-1
     1617            if ((zlay(ilon,ilat,ialt+1,it).eq.missval).or.(zlay(ilon,ilat,ialt-1,it).eq.missval)) then
     1618              delta_z(ilon,ilat,ialt,it) = 0
     1619            else
     1620              delta_z(ilon,ilat,ialt,it) = (zlay(ilon,ilat,ialt+1,it)-zlay(ilon,ilat,ialt-1,it))/2 /1000
     1621            endif ! missval
     1622           enddo ! ialt=2,altlen-1
     1623           
     1624           ! Top of last layer : we assume the same height than the layer below
     1625           delta_z(ilon,ilat,altlen,it) = delta_z(ilon,ilat,altlen-1,it)
     1626           
     1627          enddo ! it=1,timelen
     1628         enddo ! ilat=1,latlen
     1629        enddo ! ilon=1,lonlen
     1630
     1631        write(*,*) "delta_z has been well computed from variable zareoid"
     1632        write(*,*)""
     1633       
     1634      else
     1635     
     1636        write(*,*) "No variable zareoid present in GCM file"
     1637        write(*,*) "Computing the layers' height with the layers' pressure coordinate"
     1638        write(*,*) "is not (yet?) handled by the program."
     1639        write(*,*) "The column-integrated optical depths can thus not be computed."
     1640        write(*,*) "Better stop now..."
     1641        write(*,*) "You may retry with putting 'no' for the column computation."
     1642        stop
     1643       
     1644      endif ! (ierr.eq.nf90_noerr) = zareoid present in file
     1645     
     1646
     1647    ELSE
     1648      ! abort the program
     1649      write(*,*) "delta_z can not be computed from the GCM file's variables."
     1650      write(*,*) "The column-integrated optical depths can thus not be computed."
     1651      write(*,*) "Better stop now..."
     1652      write(*,*) "You may retry with putting 'no' for the column computation."
     1653      stop
     1654
     1655    ENDIF ! trim(altlong_name).eq."pseudo-alt"
     1656   
     1657  ENDIF ! is_zlev.eq.false
     1658 
     1659!=========================================================================
     1660! 3.3 Write delta_z in output file
     1661!=========================================================================
     1662 
     1663  ! Switch to netcdf define mode
     1664  ierr=nf90_redef(outfid)
     1665  ! Insert the definition of the variable
     1666  ierr=NF90_DEF_VAR(outfid,"delta_z",nf90_float,(/londimout,latdimout,altdimout,timedimout/),tmpvarid)
     1667  if (ierr.ne.nf90_noerr) then
     1668    write(*,*) "init2 Error: Failed to define the variable delta_z"
     1669    stop
     1670  endif
     1671  ! Write the attributes
     1672  ierr=nf90_put_att(outfid,tmpvarid,"long_name","Layer height used for column optical depth computation")
     1673  ierr=nf90_put_att(outfid,tmpvarid,"units","km")
     1674  ! End netcdf define mode
     1675  ierr=nf90_enddef(outfid)
     1676
     1677  ! write delta_z
     1678  ierr=NF90_PUT_VAR(outfid,tmpvarid,delta_z)
     1679  if (ierr.ne.nf90_noerr) then
     1680    write(*,*) "init2 Error: Failed to write delta_z"
     1681    stop
     1682  endif
     1683 
     1684endif ! trim(compute_col).eq."yes"
     1685
     1686!==============================================================================
     1687! 4. Cleanup
     1688!==============================================================================
     1689
    12231690if (allocated(aps)) deallocate(aps)
    12241691if (allocated(bps)) deallocate(bps)
     
    12301697
    12311698end subroutine init2
     1699
     1700
     1701!*******************************************************************************
  • trunk/LMDZ.MARS/util/aeroptical.def

    r2817 r2830  
    335.e-7
    44ext
     5yes
    56dust,dustq,reffdust,2500,optprop_dustvis_TM_n50.dat
    67h2o_ice,h2o_ice,5.e-6,920,optprop_icevis_n30.dat
     
    14153) The wavelength at which the output opacities are calibrated (in meters)
    15164) Opacity type : extinction (ext) or absorption (abs)
    16 5) aerosol_name,aerosol_mmr(variable name),aerosol_reff(variable name or single value in m),aerosol_rho(value in kg/m3),optpropfile_name
     175) Computation of the column-integrated optical depth? (yes/no)
     18
    17196) aerosol_name,aerosol_mmr(variable name),aerosol_reff(variable name or single value in m),aerosol_rho(value in kg/m3),optpropfile_name
     207) aerosol_name,aerosol_mmr(variable name),aerosol_reff(variable name or single value in m),aerosol_rho(value in kg/m3),optpropfile_name
    1821...
    1922N) aerosol_name,aerosol_mmr(variable name),aerosol_reff(variable name or single value in m),aerosol_rho(value in kg/m3),optpropfile_name
Note: See TracChangeset for help on using the changeset viewer.