program simu_MCS ! program for a MCS observer simulator of GCM data ! author : Antoine Bierjon, 2019-2020 ! contact : antoine.bierjon@lmd.jussieu.fr ! !=================================================================================================== ! PREFACE !=================================================================================================== ! This program loads on one hand, a GCM file, which can be : ! - a GCM diagfi.nc or concat.nc that covers the observation period, ! - a GCM stats file (stats.nc) which gets extended to cover the observation period ! ! On the other hand, it loads a MRO/MCS data file (binned by Luca Montabone). ! This MCS file serves then as a reference for the interpolation and the binning of the ! GCM variables contained in the GCM file and given as inputs in the simu_MCS.def ! ! Since the MCS data is binned by intervals of Ls, the GCM simulation's data is interpolated at ! the MCS spatial coordinates, at every GCM sol value contained in each interval of Ls, before ! being averaged to create output bins for each interval. The binning in sols also takes into ! account the variability of Local Time in the bin and tries to represent it (with a number of LT ! values in each sol bin equal to the number of values in each associated MCS bin : numbintemp/dust/wice, ! and centered around the MCS bin LT average : timeave). ! ! There are also specific GCM variables that the program looks for in order to make them ! comparable with their equivalent in MCS files. These variables are : ! GCM MCS ! # temp ---> dtemp,ntemp ! # dso(/dsodust/qdust)+rho ---> ddust,ndust (dust opacity/km) ! # h2o_ice+rho ---> dwice,nwice (water ice opacity/km) ! ! For dust especially, if a dust opacity variable has been created in the output file, the program also ! computes the ratio of the Column-integrated Dust Optical Depths from the binned GCM file and the ! MCS file, that can then be used to normalize the dust opacity profiles when doing comparisons. ! ! Eventually, the program outputs a netcdf file, filled with the GCM binned data for dayside ! and nightside and edited with the same format than the MCS file (and with the same ! dimensions' declaration order). ! ! Minimal requirements : ! - the MCS file must contain the variables : ! dtimeave,dtimemax,dtimemin, ! dtemp,dnumbintemp, ! ddust,dnumbindust, ! dwice,dnumbinwice, ! ntimeave,ntimemax,ntimemin, ! ntemp,nnumbintemp, ! ndust,nnumbindust, ! nwice,nnumbinwice ! - the GCM file must : ! # have the same altitude type as the MCS file ; ! # cover completely the observation period ! ! See also NOTA BENE in section 2.2 ! ! ! ! Algorithm : ! 0. Variable declarations ! 1. OPENING OF THE FILES AND INITIALIZATION OF THE DIMENSIONS ! 1.1 MCS data file : obsfile ! 1.1.1 Open the Observer data file ! 1.1.2 Get dimensions lon,lat,alt,time from the observation file ! 1.2. GCM simulation file : gcmfile ! 1.2.1 Open the GCM simulation file ! 1.2.2 Get dimensions lon,lat,alt,time from the GCM file ! 1.3 Create the output file & initialize the coordinates ! 2. VARIABLES MANAGEMENT ! 2.1 List of the GCM variables to be interpolated ! 2.1.1 Read the GCM variables ! 2.1.2 Handle dust and wice opacities (first set-ups) ! 2.2 Definition of LT variables from obsfile to outfile ![day/night loop begins... ! 2.2.1 Average of Local Time in the OBS bins ! 2.2.2 Maximum of Local Time in the OBS bins ! 2.2.3 Minimum of Local Time in the OBS bins ! 2.3 Definition of numbin variables from obsfile to outfile ! 2.3.1 Number of values in the OBS "temp" bins ! 2.3.2 Number of values in the OBS "dust" bins ! 2.3.3 Number of values in the OBS "wice" bins ! 2.4 Opening of the GCM variables ! [var loop begins... ! 2.4.1 Generic reading of the variable ! 2.4.2 Handle dust and wice opacities (second part) ! 2.5 Opening of the associated MCS variables ! 2.5.1 MCS reference variable (for the missing values) ! 2.5.2 Number of values in the OBS bin (for the sol binning) ! 2.6 Definition of GCM variables in outfile ! 3. EXTRACTION OF THE VARIABLE ! [coordinates loop begins... ! 3.1 Do some checks and preparations before the extraction ! 3.2 Compute GCM sol date corresponding to Observer Ls (via m_(min/max)sol) ! and LT (via OBSLT(min/max)) ! 3.3 Do the interpolation and binning for the given location ! ..coordinates loop ends] ! 3.4 Write the data in the netcdf output file ! ..var loop ends] ! 4. CDOD RATIO (tau_GCM/tau_MCS) ! 4.1 Preliminary checks ! 4.2 tau_ratio computation ! 4.3 Write tau_ratio in the netcdf output file ! 5. END OF THE DAY/NIGHT LOOP !..day/night loop ends] ! 6. CLOSE THE FILES ! ! Subroutines ! extraction ! inidim ! ls2sol ! gen_sol_list ! status_check ! LTmod !=================================================================================================== use netcdf !=================================================================================== ! 0. Variable declarations !=================================================================================== implicit none ! for no implicitly typed variables !------------------------ ! Files: character(len=256) :: gcmfile ! GCM simulation file logical :: is_stats = .false. ! to check if the GCM file is a stats.nc or a diagfi.nc file character(len=256) :: obsfile ! observation file = MCS/Observer data file character(len=256) :: outfile ! output file !------------------------ ! NetCDF stuff integer :: status ! NetCDF routines return code character (len=256) :: error_text ! to store the error text integer :: gcmfid ! NetCDF gcm file ID integer :: obsfid ! NetCDF observation file ID integer :: outfid ! NetCDF output file ID integer :: GCMvarid, OBSvarid, LT_id, numbin_id, outvarid ! to store the ID of a variable integer :: lat_dimid_obs,lon_dimid_obs,alt_dimid_obs,time_dimid_obs ! dimensions' ID in OBS and output files integer :: lat_dimid_gcm,lon_dimid_gcm,alt_dimid_gcm,time_dimid_gcm ! dimensions' ID in GCM file integer :: GCMvarshape(4), OBSvarshape(4), LTshape(3),numbinshape(4) ! to store a variable's coordinates order !------------------------ ! Dimensions real,dimension(:),allocatable :: GCMlon, OBSlon ! longitude in the GCM & the Observer files integer GCMlonlen, OBSlonlen ! # of grid points along GCMlon & OBSlon real,dimension(:),allocatable :: GCMlat, OBSlat ! latitude in the GCM & the Observer files integer GCMlatlen, OBSlatlen ! # of grid points along GCMlat & OBSlat real,dimension(:),allocatable :: GCMalt, OBSalt ! altitude/pressure in the GCM & the Observer files integer GCMaltlen, OBSaltlen ! # of grid point along GCMalt & OBSalt character(len=1) :: GCMalttype, OBSalttype ! altitude coord. type:'z' (altitude, m) 'p' (pressure, Pa) real,dimension(:),allocatable :: GCMtime, OBSLs ! time in the GCM diagfi (sols) & the Observer files (Ls) real,dimension(:),allocatable :: GCMstatstime ! time in the GCM stats file (LT at lon 0°) integer :: GCMtimelen, GCMstatstimelen, OBSLslen ! # of points along GCMtime, GCMstatstime, OBSLs real :: starttimeoffset=0. ! offset (in sols) wrt Ls=0 of sol 0 in GCM file !------------------------ ! Variables character(len=10) :: dayornight ! are we in the "dayside" or "nightside" loop? character(len=64),dimension(:),allocatable :: gcm_vars ! list of GCM variables to interpolate character(len=10),dimension(15) :: notprocessed ! names of the (15) variables that won't be processed integer :: nbvarfile,Nnotprocessed ! nbs of variables to deal with the non-processed ones integer :: nbvar ! nb of variables that will be processed logical :: var_ok ! is this variable to be processed? logical :: dustok1,dustok2,dustok3,wiceok ! is it possible to compute opacities and how? character(len=64) :: GCMvarname,OBSvarname,outvarname ! name of the variables integer :: nbdim ! nb of dimensions of a variable real,dimension(:,:,:,:),allocatable :: GCM_var,OBS_var ! the 4D variable extracted from GCM & OBS files real,dimension(:,:,:,:),allocatable :: GCM_rho ! atmospheric density for opacities' computation character(len=64) :: long_name,units,comment ! netcdf attributes of the variable character(len=64) :: OBSLTave_name,OBSLTmax_name,OBSLTmin_name ! names of the average, max and min of LT real,dimension(:,:,:),allocatable :: OBSLT,OBSLTmax,OBSLTmin ! 3D variables extracted from obsfile (ave, max and min of LT in a bin) character (len=64) :: numbin_name ! name of the nb of values in an OBS bin real,dimension(:,:,:,:),allocatable :: numbin ! nb of values in an OBS temp/dust/wice bin real :: GCMmiss_val, OBSmiss_val, LTmiss_val ! value to denote non-existant data real :: extr_value ! result of the extraction subroutine real, dimension(:,:,:,:), allocatable :: outvar ! outvar(,,,): 4D array to store the output variable's data !------------------------ ! Time binning management real :: OBSdeltaLs ! difference of Ls between each observation bin real :: sol, maxsol, minsol ! sol date corresponding to Observed Ls and LT integer :: m_maxsol, m_minsol ! indexes of the maximum and minimum GCM sol taken in a bin for interpolation external LTmod ! declaration of the function LTmod real :: LTmod ! declaration of the type of the function LTmod integer :: LTcount ! nb of LT samples on which the interpolation is performed, ! for every LT interval (= numbin[lon,lat,alt,Ls]) integer :: solcount ! number of GCM sol integer values in one Ls interval real,dimension(:),allocatable :: int_sol_list ! list of the integer values of GCM sol real,dimension(:),allocatable :: sol_list ! list of the sol values used for the interpolation integer :: solerrcount ! nb of GCM missing values during interpolation, removed for the binning integer :: errcount = 0 ! total number of GCM missing values real :: solbinned_value ! extracted value averaged on sol samples, which is finally put in the output bin !------------------------ ! Extraction & loop indices real :: lon_val, lat_val, alt_val, Ls_val, LT_val ! where and when the output file is written at integer :: i,j,k,l ! loop iteration indices integer :: m ! sol binning loops iteration index integer :: v,vnot ! variable loops indices !------------------------ ! CDOD ratio (tau_GCM/tau_MCS) logical :: tau_ratio_ok ! true if it is possible to compute the CDOD ratio character(len=64)::out_dztauname,OBS_dztauname,OBS_tempname ! input variable names in files integer :: out_dztauid,OBS_dztauid,OBS_tempid ! input variable ID in files real,dimension(:,:,:,:),allocatable :: out_dztau,OBS_dztau ! the 4D dust opacities from outfile and obsfile real,dimension(:,:,:,:),allocatable :: OBS_temp ! temperature from obsfile real,dimension(:),allocatable :: OBS_plev ! pressure levels encompassing each Observer pressure levels real :: OBS_rho ! temporary OBS atm density computed from OBS temp and pressure real,parameter :: r_atm=191. ! Mars atmosphere specific gas constant real,parameter :: g=3.72 ! Mars gravity real :: out_tau,OBS_tau ! integrated dust columns from GCM (outfile) and OBS real,dimension(:,:,:),allocatable :: tau_ratio ! ratio out_tau/OBS_tau !=================================================================================== ! 1. OPENING OF THE FILES AND INITIALIZATION OF THE DIMENSIONS !=================================================================================== write(*,*) "Welcome in the MRO/MCS Observer Simulator program !" !=============================================================================== ! 1.1 MCS data file : obsfile !=============================================================================== !================================================================ ! 1.1.1 Open the Observer data file !================================================================ ! Ask the user to give a netcdf observation file WRITE(*,*) "-> Enter observation file name :" READ(*,*) obsfile status=NF90_OPEN(obsfile,nf90_nowrite,obsfid) ! nowrite mode=the program can only read the opened file error_text="Error: could not open file "//trim(obsfile) call status_check(status,error_text) !================================================================ ! 1.1.2 Get dimensions lon,lat,alt,time from the observation file !================================================================ ! OBS Latitude !-------------- status=nf90_inq_dimid(obsfid,"latitude",lat_dimid_obs) error_text="Failed to find Observer latitude dimension" call status_check(status,error_text) status=nf90_inquire_dimension(obsfid,lat_dimid_obs,len=OBSlatlen) error_text="Failed to find Observer latitude length" call status_check(status,error_text) allocate(OBSlat(OBSlatlen)) status=nf90_inq_varid(obsfid,"latitude",OBSvarid) error_text="Failed to find Observer latitude ID" call status_check(status,error_text) ! Read OBSlat status=NF90_GET_VAR(obsfid,OBSvarid,OBSlat) error_text="Failed to load OBSlat" call status_check(status,error_text) ! OBS Longitude !-------------- status=nf90_inq_dimid(obsfid,"longitude",lon_dimid_obs) error_text="Failed to find Observer longitude dimension" call status_check(status,error_text) status=nf90_inquire_dimension(obsfid,lon_dimid_obs,len=OBSlonlen) error_text="Failed to find Observer longitude length" call status_check(status,error_text) allocate(OBSlon(OBSlonlen)) status=nf90_inq_varid(obsfid,"longitude",OBSvarid) error_text="Failed to find Observer longitude ID" call status_check(status,error_text) ! Read OBSlon status=NF90_GET_VAR(obsfid,OBSvarid,OBSlon) error_text="Failed to load OBSlon" call status_check(status,error_text) ! OBS Time (Ls) !-------------- status=nf90_inq_dimid(obsfid,"time",time_dimid_obs) error_text="Failed to find Observer time (Ls) dimension" call status_check(status,error_text) status=nf90_inquire_dimension(obsfid,time_dimid_obs,len=OBSLslen) error_text="Failed to find Observer time (Ls) length" call status_check(status,error_text) allocate(OBSLs(OBSLslen)) status=nf90_inq_varid(obsfid,"time",OBSvarid) error_text="Failed to find Observer time (Ls) ID" call status_check(status,error_text) ! Read OBSLs status=NF90_GET_VAR(obsfid,OBSvarid,OBSLs) error_text="Failed to load OBSLs" call status_check(status,error_text) ! Get the observation timestep between bins OBSdeltaLs=OBSLs(2)-OBSLs(1) ! OBS Altitude !-------------- status=nf90_inq_dimid(obsfid,"altitude",alt_dimid_obs) error_text="Failed to find Observer altitude dimension" call status_check(status,error_text) status=nf90_inquire_dimension(obsfid,alt_dimid_obs,len=OBSaltlen) error_text="Failed to find Observer altitude length" call status_check(status,error_text) allocate(OBSalt(OBSaltlen)) status=nf90_inq_varid(obsfid,"altitude",OBSvarid) error_text="Failed to find Observer altitude ID" call status_check(status,error_text) ! Read OBSalt status=NF90_GET_VAR(obsfid,OBSvarid,OBSalt) error_text="Failed to load OBSalt" call status_check(status,error_text) ! Check altitude attribute "units" to find out altitude type and compare with the GCM file status=nf90_get_att(obsfid,OBSvarid,"units",units) error_text="Failed to load Observer altitude units attribute" call status_check(status,error_text) ! an unknown and invisible character is placed just after the unit's ! characters in the Observer file so we only take the first characters ! corresponding to the sought unit if (trim(units(1:2)).eq."Pa") then units="Pa" OBSalttype='p' ! pressure coordinate else if (trim(units(1:2)).eq."m") then units="m" OBSalttype='z' ! altitude coordinate else write(*,*)" I do not understand this unit ",trim(units)," for Observer altitude!" stop endif !=============================================================================== ! 1.2. GCM simulation file : gcmfile !=============================================================================== !================================================================ ! 1.2.1 Open the GCM simulation file !================================================================ ! Ask the user to give a netcdf input file write(*,*)"";WRITE(*,*) "-> Enter input file name (GCM simulation) :" READ(*,*) gcmfile ! Open GCM file status=NF90_OPEN(gcmfile,nf90_nowrite,gcmfid) ! nowrite mode=the program can only read the opened file error_text="Failed to open datafile "//trim(gcmfile) call status_check(status,error_text) !================================================================ ! 1.2.2 Get dimensions lon,lat,alt,time from the GCM file !================================================================ ! GCM Latitude !-------------- status=nf90_inq_dimid(gcmfid,"latitude",lat_dimid_gcm) error_text="Failed to find GCM latitude dimension" call status_check(status,error_text) status=nf90_inquire_dimension(gcmfid,lat_dimid_gcm,len=GCMlatlen) error_text="Failed to find GCM latitude length" call status_check(status,error_text) allocate(GCMlat(GCMlatlen)) status=nf90_inq_varid(gcmfid,"latitude",GCMvarid) error_text="Failed to find GCM latitude ID" call status_check(status,error_text) ! Read GCMlat status=NF90_GET_VAR(gcmfid,GCMvarid,GCMlat) error_text="Failed to load GCMlat" call status_check(status,error_text) ! GCM Longitude !-------------- status=nf90_inq_dimid(gcmfid,"longitude",lon_dimid_gcm) error_text="Failed to find GCM longitude dimension" call status_check(status,error_text) status=nf90_inquire_dimension(gcmfid,lon_dimid_gcm,len=GCMlonlen) error_text="Failed to find GCM longitude length" call status_check(status,error_text) allocate(GCMlon(GCMlonlen)) status=nf90_inq_varid(gcmfid,"longitude",GCMvarid) error_text="Failed to find GCM longitude ID" call status_check(status,error_text) ! Read GCMlon status=NF90_GET_VAR(gcmfid,GCMvarid,GCMlon) error_text="Failed to load GCMlon" call status_check(status,error_text) ! GCM Time !-------------- status=nf90_inq_dimid(gcmfid,"Time",time_dimid_gcm) error_text="Failed to find GCM time dimension" call status_check(status,error_text) status=nf90_inquire_dimension(gcmfid,time_dimid_gcm,len=GCMtimelen) error_text="Failed to find GCM time length" call status_check(status,error_text) allocate(GCMtime(GCMtimelen)) status=nf90_inq_varid(gcmfid,"Time",GCMvarid) error_text="Failed to find GCM time ID" call status_check(status,error_text) status=NF90_GET_VAR(gcmfid,GCMvarid,GCMtime) error_text="Failed to load GCMtime" call status_check(status,error_text) ! is_stats ? IF ((GCMtimelen.eq.12).and.(GCMtime(1).eq.2.).and.(GCMtime(GCMtimelen).eq.24.)) then ! if GCM file is a stats, time is in LT at longitude 0° and not in sols at longitude 0° write(*,*)"The GCM file is recognized as a stats file." is_stats = .true. deallocate(GCMtime) GCMstatstimelen = GCMtimelen allocate(GCMstatstime(GCMstatstimelen)) status=NF90_GET_VAR(gcmfid,GCMvarid,GCMstatstime) error_text="Failed to load GCMstatstime (LT at lon 0)" call status_check(status,error_text) ELSE write(*,*)"The GCM file is recognized as a diagfi/concatnc file." ENDIF ! Simulation time offset management WRITE(*,*) "Beginning date of the simulation file?" WRITE(*,*) "(i.e. number of sols since Ls=0 at the Time=0.0 in the GCM file)" READ(*,*) starttimeoffset if (.not.is_stats) then ! Add the offset to GCMtime(:) if the file is not a stats file GCMtime(:)=GCMtime(:)+starttimeoffset endif ! Check of temporal coherence between gcmfile & obsfile call ls2sol(OBSLs(OBSLslen),maxsol) ! maximum date considered call ls2sol(OBSLs(1),minsol) ! minimum date considered IF (.not.is_stats) then ! if it is a diagfi, we check the time coherence between the 2 files if ((maxsol.gt.maxval(GCMtime)).or.(minsol.lt.minval(GCMtime))) then write(*,*)"Error : obsfile temporal bounds exceed the GCM simulation bounds." write(*,*)"Please use a GCM file whose time interval contains the observation period." stop else write(*,*)"Both files are temporally coherent. The program continues..." endif ELSE ! if it is a stats, we create the array GCMtime array (in sols) covering the observation period ! and filled with the mean GCM day stored in stats.nc GCMtimelen = ((ceiling(maxsol)-floor(minsol)+1)+2) ! we add 2 days in the beginning and the end ! to be sure we cover the observation period allocate(GCMtime(GCMstatstimelen * GCMtimelen)) do l=1,GCMtimelen do m=1,GCMstatstimelen GCMtime(m+(l-1)*GCMstatstimelen) = (floor(minsol)-1) + (l-1) + GCMstatstime(m)/24. enddo enddo GCMtimelen = GCMstatstimelen * GCMtimelen write(*,*)"GCMtime has been created from the stats.nc time and the observation period. The program continues..." ENDIF ! GCM Altitude !-------------- status=nf90_inq_dimid(gcmfid,"altitude",alt_dimid_gcm) error_text="Failed to find GCM altitude dimension" call status_check(status,error_text) status=nf90_inquire_dimension(gcmfid,alt_dimid_gcm,len=GCMaltlen) error_text="Failed to find GCM altitude length" call status_check(status,error_text) allocate(GCMalt(GCMaltlen)) status=nf90_inq_varid(gcmfid,"altitude",GCMvarid) error_text="Failed to find GCM altitude ID" call status_check(status,error_text) ! Read GCMalt status=NF90_GET_VAR(gcmfid,GCMvarid,GCMalt) error_text="Failed to load GCMalt" call status_check(status,error_text) ! Check altitude attribute "units" to find out altitude type status=nf90_get_att(gcmfid,GCMvarid,"units",units) error_text="Failed to load GCM altitude units attribute" call status_check(status,error_text) if (trim(units).eq."Pa") then GCMalttype='p' ! pressure coordinate else if (trim(units).eq."m") then GCMalttype='z' ! altitude coordinate else write(*,*)"I do not understand this unit ",trim(units)," for GCM altitude!" if (OBSalttype.eq.'p') then write(*,*)"Please use zrecast to put the altitude in the same type as the MCS file (pressure in Pa)" else if (OBSalttype.eq.'z') then write(*,*)"Please use zrecast to put the altitude in the same type as the MCS file (altitude in m)" endif stop endif IF(OBSalttype.ne.GCMalttype) then write(*,*)"Observer altitude type (", OBSalttype,") and ", & "GCM altitude type (",GCMalttype,") don't match!" stop ENDIF !=============================================================================== ! 1.3 Create the output file & initialize the coordinates !=============================================================================== ! Name of the outfile IF (.not.is_stats) then outfile=obsfile(1:index(obsfile, ".nc")-1)//"_GCMdiagfi.nc" ELSE outfile=obsfile(1:index(obsfile, ".nc")-1)//"_GCMstats.nc" ENDIF ! Creation of the outfile status=NF90_CREATE(outfile,nf90_clobber,outfid)!NB: clobber mode=overwrite any dataset with the same file name ! error_text="Error: could not create file "//trim(outfile) call status_check(status,error_text) write(*,*)"";WRITE(*,*)"-> Output file is: ",trim(outfile) ! Creation of the dimensions call inidim(outfid,OBSlonlen,OBSlatlen,OBSaltlen,OBSLslen,OBSlon,OBSlat,OBSalt,OBSLs,OBSalttype,& lon_dimid_obs,lat_dimid_obs,alt_dimid_obs,time_dimid_obs) !write(*,*)"Dimensions ID in the outfile are :" !write(*,*)"lon_dimid=",lon_dimid_obs !write(*,*)"lat_dimid=",lat_dimid_obs !write(*,*)"alt_dimid=",alt_dimid_obs !write(*,*)"time_dimid=",time_dimid_obs !=================================================================================== ! 2. VARIABLES MANAGEMENT !=================================================================================== !=============================================================================== ! 2.1 List of the GCM variables to be interpolated !=============================================================================== !================================================================ ! 2.1.1 Read the GCM variables !================================================================ ! Initialize logicals dustok1 = .false. dustok2 = .false. dustok3 = .false. wiceok = .false. ! Get nbvarfile (total number of variables in the GCM file) status=NF90_INQUIRE(gcmfid,nVariables=nbvarfile) error_text="Error : Pb with nf90_inquire(gcmfid,nVariables=nbvarfile)" call status_check(status,error_text) ! List of variables that should not be processed notprocessed(1)='Time' notprocessed(2)='controle' notprocessed(3)='rlonu' notprocessed(4)='latitude' notprocessed(5)='longitude' notprocessed(6)='altitude' notprocessed(7)='rlatv' notprocessed(8)='aps' notprocessed(9)='bps' notprocessed(10)='ap' notprocessed(11)='bp' notprocessed(12)='cu' notprocessed(13)='cv' notprocessed(14)='aire' notprocessed(15)='phisinit' ! List of variables in the GCM file write(*,*)"" write(*,*)"List of variables in the GCM file :" Nnotprocessed=0 do v=1,nbvarfile status=NF90_INQUIRE_VARIABLE(gcmfid,v,name=GCMvarname) ! GCMvarname now contains the "name" of variable of ID # v var_ok=.true. do vnot=1,15 if (GCMvarname.eq.notprocessed(vnot)) then var_ok=.false. Nnotprocessed=Nnotprocessed+1 endif enddo if (var_ok) write(*,*) trim(GCMvarname) ! Detect if we can compute dust and wice opacities if (trim(GCMvarname).eq."dso") then dustok1 = .true. else if (trim(GCMvarname).eq."dsodust") then dustok2 = .true. else if (trim(GCMvarname).eq."dustq") then dustok3 = .true. else if (trim(GCMvarname).eq."h2o_ice") then wiceok = .true. endif enddo ! Nnotprocessed: # of variables that won't be processed ! nbvarfile: total # of variables in file ! +2: the dust and wice opacities allocate(gcm_vars(nbvarfile-Nnotprocessed+2),stat=status) if (status.ne.0) then write(*,*) "Error: failed allocation of gcm_vars(nbvarfile-Nnotprocessed+2)" write(*,*) " nbvarfile=",nbvarfile write(*,*) " Nnotprocessed=",Nnotprocessed stop endif ! List of variables to process write(*,*) write(*,*) "Which variables do you want to redistribute ?" write(*,*) "list of (separated by s)" write(*,*) "(an empty line , i.e: just , implies end of list)" write(*,*) "NB: this program handles only 4D netcdf variables for now" nbvar=0 read(*,'(a50)') GCMvarname do while ((GCMvarname.ne.' ').AND.(trim(GCMvarname).ne."all")) nbvar=nbvar+1 gcm_vars(nbvar)=GCMvarname read(*,'(a50)') GCMvarname enddo if (GCMvarname.eq."all") then nbvar=nbvarfile-Nnotprocessed do v=Nnotprocessed+1,nbvarfile status=nf90_inquire_variable(gcmfid,v,name=gcm_vars(v-Nnotprocessed)) enddo ! Variables names from the file are stored in gcm_vars() nbvar=nbvarfile-Nnotprocessed do v=1,nbvar status=nf90_inquire_variable(gcmfid,v+Nnotprocessed,name=gcm_vars(v)) write(*,'(a9,1x,i2,1x,a1,1x,a64)') "variable ",v,":",gcm_vars(v) enddo else if(nbvar==0) then write(*,*) "No variables to process in the GCM file... program stopped" stop endif ! of if (GCMvarname.eq."all") !================================================================ ! 2.1.2 Handle dust and wice opacities (first set-ups) !================================================================ ! 2nd part is in section 2.4.2 write(*,*) ! Load atmospheric density "rho" if (dustok1.or.dustok2.or.dustok3.or.wiceok) then ! Check that the GCM file contains that variable status=nf90_inq_varid(gcmfid,"rho",GCMvarid) if (status.ne.nf90_noerr) then write(*,*) "Failed to find variable rho in "//trim(gcmfile) write(*,*) "No computation of opacities..." dustok1 =.false. dustok2 =.false. dustok3 =.false. wiceok =.false. else ! Length allocation for each dimension of the 4D variable allocate(GCM_rho(GCMlonlen,GCMlatlen,GCMaltlen,GCMtimelen)) ! Load datasets if (.not.is_stats) then status=NF90_GET_VAR(gcmfid,GCMvarid,GCM_rho) error_text="Failed to load rho" call status_check(status,error_text) else ! if it is a stats file, we load only the first sol, and then copy it to all the other sols status=NF90_GET_VAR(gcmfid,GCMvarid,GCM_rho(:,:,:,1:GCMstatstimelen)) error_text="Failed to load rho" call status_check(status,error_text) ! write(*,*)"GCMstatstimelen = ", GCMstatstimelen ! write(*,*)"GCMtimelen = ", GCMtimelen do l=(GCMstatstimelen+1),GCMtimelen if (modulo(l,GCMstatstimelen).ne.0) then GCM_rho(:,:,:,l) = GCM_rho(:,:,:,modulo(l,GCMstatstimelen)) else ! if l is a multiple of GCMstatstimelen, since the index modulo(l,GCMstatstimelen)=0 ! doesn't exist, we make a special case GCM_rho(:,:,:,l) = GCM_rho(:,:,:,GCMstatstimelen) endif enddo endif write(*,*) "Variable rho loaded from the GCM file" endif endif ! dustok1.or.dustok2.or.dustok3.or.wiceok ! Dust and wice opacity booleans if (dustok1.or.dustok2.or.dustok3) then nbvar=nbvar+1 gcm_vars(nbvar)="dust" endif if (wiceok) then nbvar=nbvar+1 gcm_vars(nbvar)="wice" endif !write(*,*) "gcm_vars retrieved : ",gcm_vars(1:nbvar) !=============================================================================== ! 2.2 Definition of LT variables from obsfile to outfile !=============================================================================== ! --> the day/night loop begins here !******************** NOTA BENE (cf sections 2.2 and 4)************************* ! We execute the program a first time with the daytime values, and then a second ! time with the nighttime values. !******************************************************************************* dayornight = "dayside" ! we begin with daytime temperature write(*,*)"" ; write(*,*) "Beginning the 1st loop, on daytime values"; write(*,*)"" DAY_OR_NIGHT: DO ! (the end of the loop is in section 4.) SELECT CASE (dayornight) CASE ("dayside") OBSLTave_name = "dtimeave" OBSLTmax_name = "dtimemax" OBSLTmin_name = "dtimemin" CASE ("nightside") OBSLTave_name = "ntimeave" OBSLTmax_name = "ntimemax" OBSLTmin_name = "ntimemin" END SELECT !================================================================ ! 2.2.1 Average of Local Time in the OBS bins !================================================================ ! Read the OBS file !------------------ status=nf90_inq_varid(obsfid,trim(OBSLTave_name),LT_id) error_text="Failed to find Observer local time ("//trim(OBSLTave_name)//") ID in "//trim(obsfile) call status_check(status,error_text) status=nf90_inquire_variable(obsfid,LT_id,dimids=LTshape) error_text="Failed to get the dim shape of variable "//trim(OBSLTave_name) call status_check(status,error_text) ! Length allocation for each dimension of the 3D variable allocate(OBSLT(OBSlonlen,OBSlatlen,OBSLslen)) ! Load datasets status=NF90_GET_VAR(obsfid,LT_id,OBSLT) error_text="Failed to load "//trim(OBSLTave_name)//" from the obsfile" call status_check(status,error_text) write(*,*) trim(OBSLTave_name)," loaded from the obsfile" ! Get LT missing_value attribute status=nf90_get_att(obsfid,LT_id,"_FillValue",LTmiss_val) error_text="Failed to load missing_value attribute" call status_check(status,error_text) ! Create the variable in the outfile !----------------------------------- ! Switch to netcdf define mode status=nf90_redef(outfid) error_text="Error: could not switch to define mode in the outfile" call status_check(status,error_text) ! Definition of the variable status=NF90_DEF_VAR(outfid,trim(OBSLTave_name),nf90_float,LTshape,LT_id) error_text="Error: could not define the variable "//trim(OBSLTave_name)//" in the outfile" call status_check(status,error_text) ! Write the attributes select case (dayornight) case ("dayside") status=nf90_put_att(outfid,LT_id,"long_name","Average local time in bin - day side [6h, 18h]") case ("nightside") status=nf90_put_att(outfid,LT_id,"long_name","Average local time in bin - night side [18h, 6h]") end select status=nf90_put_att(outfid,LT_id,"units","hours") status=nf90_put_att(outfid,LT_id,"_FillValue",LTmiss_val) ! End the netcdf define mode (and thus enter the "data writing" mode) status=nf90_enddef(outfid) error_text="Error: could not close the define mode of the outfile" call status_check(status,error_text) ! Write the data in the output file status = NF90_PUT_VAR(outfid, LT_id, OBSLT) ! write the MCS d/ntimeave as the output Local Time error_text="Error: could not write "//trim(OBSLTave_name)//" data in the outfile" call status_check(status,error_text) write(*,*)"Local Time (",trim(OBSLTave_name),") has been created in the outfile" write(*,'(" with missing_value attribute : ",1pe12.5)')LTmiss_val !================================================================ ! 2.2.2 Maximum of Local Time in the OBS bins !================================================================ ! Read the OBS file !------------------ status=nf90_inq_varid(obsfid,trim(OBSLTmax_name),LT_id) error_text="Failed to find Observer max local time ("//trim(OBSLTmax_name)//") ID in "//trim(obsfile) call status_check(status,error_text) status=nf90_inquire_variable(obsfid,LT_id,dimids=LTshape) error_text="Failed to get the dim shape of variable "//trim(OBSLTmax_name) call status_check(status,error_text) ! Length allocation for each dimension of the 3D variable allocate(OBSLTmax(OBSlonlen,OBSlatlen,OBSLslen)) ! Load datasets status=NF90_GET_VAR(obsfid,LT_id,OBSLTmax) error_text="Failed to load "//trim(OBSLTmax_name)//" from the obsfile" call status_check(status,error_text) write(*,*) trim(OBSLTmax_name)," loaded from the obsfile" !================================================================ ! 2.2.3 Minimum of Local Time in the OBS bins !================================================================ ! Read the OBS file !------------------ status=nf90_inq_varid(obsfid,trim(OBSLTmin_name),LT_id) error_text="Failed to find Observer min local time ("//trim(OBSLTmin_name)//") ID in "//trim(obsfile) call status_check(status,error_text) status=nf90_inquire_variable(obsfid,LT_id,dimids=LTshape) error_text="Failed to obtain information on variable "//trim(OBSLTmin_name) call status_check(status,error_text) ! Length allocation for each dimension of the 3D variable allocate(OBSLTmin(OBSlonlen,OBSlatlen,OBSLslen)) ! Load datasets status=NF90_GET_VAR(obsfid,LT_id,OBSLTmin) error_text="Failed to load "//trim(OBSLTmin_name)//" from the obsfile" call status_check(status,error_text) write(*,*) trim(OBSLTmin_name)," loaded from the obsfile" write(*,*)"" !=============================================================================== ! 2.3 Definition of numbin variables from obsfile to outfile !=============================================================================== !================================================================ ! 2.3.1 Number of values in the OBS "temp" bins !================================================================ SELECT CASE (dayornight) CASE ("dayside") numbin_name = "dnumbintemp" CASE ("nightside") numbin_name = "nnumbintemp" END SELECT ! Read the OBS file !------------------ status=nf90_inq_varid(obsfid,trim(numbin_name),numbin_id) error_text="Failed to find Observer nb of temp values in bin ("//trim(numbin_name)//")'s ID in "//trim(obsfile) call status_check(status,error_text) status=nf90_inquire_variable(obsfid,numbin_id,dimids=numbinshape) error_text="Failed to obtain information on variable "//trim(numbin_name) call status_check(status,error_text) ! Length allocation for each dimension of the 4D variable allocate(numbin(OBSlonlen,OBSlatlen,OBSaltlen,OBSLslen)) ! Load datasets status=NF90_GET_VAR(obsfid,numbin_id,numbin) error_text="Failed to load "//trim(numbin_name)//" from the obsfile" call status_check(status,error_text) write(*,*) trim(numbin_name)," loaded from the obsfile" ! Create the variable in the outfile !----------------------------------- ! Switch to netcdf define mode status=nf90_redef(outfid) error_text="Error: could not switch to define mode in the outfile" call status_check(status,error_text) ! Definition of the variable status=NF90_DEF_VAR(outfid,trim(numbin_name),nf90_float,numbinshape,numbin_id) error_text="Error: could not define the variable "//trim(numbin_name)//" in the outfile" call status_check(status,error_text) ! Write the attributes select case (dayornight) case ("dayside") status=nf90_put_att(outfid,numbin_id,"long_name","Number of temp values in bin - day side") case ("nightside") status=nf90_put_att(outfid,numbin_id,"long_name","Number of temp values in bin - night side") end select ! End the netcdf define mode (and thus enter the "data writing" mode) status=nf90_enddef(outfid) error_text="Error: could not close the define mode of the outfile" call status_check(status,error_text) ! Write the data in the output file status = NF90_PUT_VAR(outfid, numbin_id, numbin) error_text="Error: could not write "//trim(numbin_name)//" data in the outfile" call status_check(status,error_text) write(*,*)"Number of temp values in bin (",trim(numbin_name),") has been created in the outfile" write(*,*)"" deallocate(numbin) !================================================================ ! 2.3.2 Number of values in the OBS "dust" bins !================================================================ SELECT CASE (dayornight) CASE ("dayside") numbin_name = "dnumbindust" CASE ("nightside") numbin_name = "nnumbindust" END SELECT ! Read the OBS file !------------------ status=nf90_inq_varid(obsfid,trim(numbin_name),numbin_id) error_text="Failed to find Observer nb of dust values in bin ("//trim(numbin_name)//")'s ID in "//trim(obsfile) call status_check(status,error_text) status=nf90_inquire_variable(obsfid,numbin_id,dimids=numbinshape) error_text="Failed to obtain information on variable "//trim(numbin_name) call status_check(status,error_text) ! Length allocation for each dimension of the 4D variable allocate(numbin(OBSlonlen,OBSlatlen,OBSaltlen,OBSLslen)) ! Load datasets status=NF90_GET_VAR(obsfid,numbin_id,numbin) error_text="Failed to load "//trim(numbin_name)//" from the obsfile" call status_check(status,error_text) write(*,*) trim(numbin_name)," loaded from the obsfile" ! Create the variable in the outfile !----------------------------------- ! Switch to netcdf define mode status=nf90_redef(outfid) error_text="Error: could not switch to define mode in the outfile" call status_check(status,error_text) ! Definition of the variable status=NF90_DEF_VAR(outfid,trim(numbin_name),nf90_float,numbinshape,numbin_id) error_text="Error: could not define the variable "//trim(numbin_name)//" in the outfile" call status_check(status,error_text) ! Write the attributes select case (dayornight) case ("dayside") status=nf90_put_att(outfid,numbin_id,"long_name","Number of dust values in bin - day side") case ("nightside") status=nf90_put_att(outfid,numbin_id,"long_name","Number of dust values in bin - night side") end select ! End the netcdf define mode (and thus enter the "data writing" mode) status=nf90_enddef(outfid) error_text="Error: could not close the define mode of the outfile" call status_check(status,error_text) ! Write the data in the output file status = NF90_PUT_VAR(outfid, numbin_id, numbin) error_text="Error: could not write "//trim(numbin_name)//" data in the outfile" call status_check(status,error_text) write(*,*)"Number of dust values in bin (",trim(numbin_name),") has been created in the outfile" write(*,*)"" deallocate(numbin) !================================================================ ! 2.3.3 Number of values in the OBS "wice" bins !================================================================ SELECT CASE (dayornight) CASE ("dayside") numbin_name = "dnumbinwice" CASE ("nightside") numbin_name = "nnumbinwice" END SELECT ! Read the OBS file !------------------ status=nf90_inq_varid(obsfid,trim(numbin_name),numbin_id) error_text="Failed to find Observer nb of wice values in bin ("//trim(numbin_name)//")'s ID in "//trim(obsfile) call status_check(status,error_text) status=nf90_inquire_variable(obsfid,numbin_id,dimids=numbinshape) error_text="Failed to obtain information on variable "//trim(numbin_name) call status_check(status,error_text) ! Length allocation for each dimension of the 4D variable allocate(numbin(OBSlonlen,OBSlatlen,OBSaltlen,OBSLslen)) ! Load datasets status=NF90_GET_VAR(obsfid,numbin_id,numbin) error_text="Failed to load "//trim(numbin_name)//" from the obsfile" call status_check(status,error_text) write(*,*) trim(numbin_name)," loaded from the obsfile" ! Create the variable in the outfile !----------------------------------- ! Switch to netcdf define mode status=nf90_redef(outfid) error_text="Error: could not switch to define mode in the outfile" call status_check(status,error_text) ! Definition of the variable status=NF90_DEF_VAR(outfid,trim(numbin_name),nf90_float,numbinshape,numbin_id) error_text="Error: could not define the variable "//trim(numbin_name)//" in the outfile" call status_check(status,error_text) ! Write the attributes select case (dayornight) case ("dayside") status=nf90_put_att(outfid,numbin_id,"long_name","Number of wice values in bin - day side") case ("nightside") status=nf90_put_att(outfid,numbin_id,"long_name","Number of wice values in bin - night side") end select ! End the netcdf define mode (and thus enter the "data writing" mode) status=nf90_enddef(outfid) error_text="Error: could not close the define mode of the outfile" call status_check(status,error_text) ! Write the data in the output file status = NF90_PUT_VAR(outfid, numbin_id, numbin) error_text="Error: could not write "//trim(numbin_name)//" data in the outfile" call status_check(status,error_text) write(*,*)"Number of wice values in bin (",trim(numbin_name),") has been created in the outfile" write(*,*)"" deallocate(numbin) !=============================================================================== ! 2.4 Opening of the GCM variables !=============================================================================== ! --> the var loop begins here VAR: DO v=1,nbvar ! LOOP ON ALL THE GCM VARIABLES TO PROCESS ! (the end of the loop is in section 3.4) GCMvarname = gcm_vars(v) ! Detect the dust and wice opacities special cases if (trim(GCMvarname).eq."dust") then ! Default : dso > dsodust > dustq if (dustok1) then ! "dso" is detected in gcmfile GCMvarname="dso" dustok2=.false. dustok3=.false. else if (dustok2) then ! "dsodust" is detected in gcmfile GCMvarname="dsodust" dustok3=.false. else if (dustok3) then ! "dustq" is detected in gcmfile GCMvarname="dustq" endif write(*,*) "Computing dust opacity..." endif if (trim(GCMvarname).eq."wice") then ! "h2o_ice" detected in gcmfile GCMvarname="h2o_ice" write(*,*) "Computing water ice opacity..." endif !================================================================ ! 2.4.1 Generic reading of the variable !================================================================ ! Check that the GCM file contains that variable status=nf90_inq_varid(gcmfid,trim(GCMvarname),GCMvarid) if (status.ne.nf90_noerr) then write(*,*) "Failed to find variable "//trim(GCMvarname)//" in "//trim(gcmfile) write(*,*) "We'll skip that variable..." CYCLE VAR ! go directly to the next variable endif ! Sanity checks on the variable status=nf90_inquire_variable(gcmfid,GCMvarid,ndims=nbdim,dimids=GCMvarshape) error_text="Failed to obtain information on variable "//trim(GCMvarname) call status_check(status,error_text) ! Check that it is a 4D variable if (nbdim.ne.4) then write(*,*) "Error:",trim(GCMvarname)," is not a 4D variable" write(*,*) "We'll skip that variable...";write(*,*)"" CYCLE VAR ! go directly to the next variable endif ! Check that its dimensions are indeed lon,lat,alt,time (in the right order) if (GCMvarshape(1).ne.lon_dimid_gcm) then write(*,*) "Error, expected first dimension of ",trim(GCMvarname)," to be longitude!" write(*,*) "We'll skip that variable..." CYCLE VAR ! go directly to the next variable endif if (GCMvarshape(2).ne.lat_dimid_gcm) then write(*,*) "Error, expected second dimension of ",trim(GCMvarname)," to be latitude!" write(*,*) "We'll skip that variable..." CYCLE VAR ! go directly to the next variable endif if (GCMvarshape(3).ne.alt_dimid_gcm) then write(*,*) "Error, expected third dimension of ",trim(GCMvarname)," to be altitude!" write(*,*) "We'll skip that variable..." CYCLE VAR ! go directly to the next variable endif if (GCMvarshape(4).ne.time_dimid_gcm) then write(*,*) "Error, expected fourth dimension of ",trim(GCMvarname)," to be time!" write(*,*) "We'll skip that variable..." CYCLE VAR ! go directly to the next variable endif ! Length allocation for each dimension of the 4D variable allocate(GCM_var(GCMlonlen,GCMlatlen,GCMaltlen,GCMtimelen)) ! Load datasets if (.not.is_stats) then status=NF90_GET_VAR(gcmfid,GCMvarid,GCM_var) error_text="Failed to load "//trim(GCMvarname) call status_check(status,error_text) else ! if it is a stats file, we load only the first sol, and then copy it to all the other sols status=NF90_GET_VAR(gcmfid,GCMvarid,GCM_var(:,:,:,1:GCMstatstimelen)) error_text="Failed to load "//trim(GCMvarname) call status_check(status,error_text) ! write(*,*)"GCMstatstimelen = ", GCMstatstimelen ! write(*,*)"GCMtimelen = ", GCMtimelen do l=(GCMstatstimelen+1),GCMtimelen if (modulo(l,GCMstatstimelen).ne.0) then GCM_var(:,:,:,l) = GCM_var(:,:,:,modulo(l,GCMstatstimelen)) else ! if l is a multiple of GCMstatstimelen, since the index modulo(l,GCMstatstimelen)=0 ! doesn't exist, we make a special case GCM_var(:,:,:,l) = GCM_var(:,:,:,GCMstatstimelen) endif enddo endif write(*,*) "Variable ",trim(GCMvarname)," loaded from the GCM file" ! Get dataset's missing_value attribute status=nf90_get_att(gcmfid,GCMvarid,"missing_value",GCMmiss_val) error_text="Failed to load missing_value attribute" call status_check(status,error_text) ! Get other variable's attributes status=nf90_get_att(gcmfid,GCMvarid,"long_name",long_name) if (status.ne.nf90_noerr) then ! if no attribute "long_name", try "title" status=nf90_get_att(gcmfid,GCMvarid,"title",long_name) endif status=nf90_get_att(gcmfid,GCMvarid,"units",units) !================================================================ ! 2.4.2 Handle dust and wice opacities (second part) !================================================================ ! DUST !----- if (trim(gcm_vars(v)).eq."dust") then IF (dustok1.or.dustok2) THEN ! Dust opacity computed from its density-scaled opacity ! write(*,*)long_name(index(long_name,"(")+1:index(long_name,")")-1) do i=1,GCMlonlen do j=1,GCMlatlen do k=1,GCMaltlen do l=1,GCMtimelen if (GCM_var(i,j,k,l).ne.GCMmiss_val) then ! Multiply by rho to have opacity [1/km] GCM_var(i,j,k,l) = GCM_var(i,j,k,l) * GCM_rho(i,j,k,l) *1000. if (long_name(index(long_name,"(")+1:index(long_name,")")-1).eq."TES") then ! The density-scaled opacity was calibrated on TES wavelength (9.3um) ! so we must recalibrate it to MCS wavelength (21.6um) using recalibration ! coefficients from Montabone et al. 2015, section 2.3 GCM_var(i,j,k,l) = 1.3/2.7 * GCM_var(i,j,k,l) endif endif enddo enddo enddo enddo long_name = "IR Dust opacity (from DSO)" ELSE IF (dustok3) THEN ! Dust opacity computed from its mass mixing ratio do i=1,GCMlonlen do j=1,GCMlatlen do k=1,GCMaltlen do l=1,GCMtimelen if (GCM_var(i,j,k,l).ne.GCMmiss_val) then ! Opacity is computed from the equation of Heavens et al. 2014, section 2.3, ! assuming a rho_dust=3000 kg/m3 and an effective radius reff=1.06 microns ! + the opacity is here in 1/km and the MMR in kg/kg GCM_var(i,j,k,l) = GCM_var(i,j,k,l) * GCM_rho(i,j,k,l) / 0.012 * 1000 endif enddo enddo enddo enddo long_name = "IR Dust opacity (from MMR)" ENDIF GCMvarname = gcm_vars(v) ! reput the right name in GCMvarname units = "opacity/km" endif ! trim(gcm_vars(v)).eq."dust" ! WICE !----- if (trim(gcm_vars(v)).eq."wice") then ! Water ice opacity computed from its mass mixing ratio do i=1,GCMlonlen do j=1,GCMlatlen do k=1,GCMaltlen do l=1,GCMtimelen if (GCM_var(i,j,k,l).ne.GCMmiss_val) then ! Opacity at MCS wavelength (11.9um) is computed from an equation ! similar to the one of Heavens et al. 2014, section 2.3. ! We assume a rho_wice=920 kg/m3, an effective radius reff=3um, ! an extinction coefficient Qext(wvl,reff)=1.54471 GCM_var(i,j,k,l) = 750*1.54471* GCM_var(i,j,k,l) * GCM_rho(i,j,k,l) / (920*3e-6) endif enddo enddo enddo enddo long_name = "IR Water ice opacity (from MMR)" GCMvarname = gcm_vars(v) ! reput the right name in GCMvarname units = "opacity/km" endif ! trim(gcm_vars(v)).eq."wice" !=============================================================================== ! 2.5 Opening of the associated MCS variables !=============================================================================== ! Observer variables to extract : 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").or.(trim(GCMvarname).eq."opawice")) THEN OBSvarname = "wice" numbin_name = "numbinwice" ELSE ! default case is temp binning, since it contains the most values OBSvarname = "temp" numbin_name = "numbintemp" ENDIF SELECT CASE (dayornight) CASE ("dayside") OBSvarname = "d"//OBSvarname numbin_name = "d"//numbin_name CASE ("nightside") OBSvarname = "n"//OBSvarname numbin_name = "n"//numbin_name END SELECT !================================================================ ! 2.5.1 MCS reference variable (for the missing values) !================================================================ ! Check that the observation file contains that variable status=nf90_inq_varid(obsfid,trim(OBSvarname),OBSvarid) error_text="Failed to find variable "//trim(OBSvarname)//" in "//trim(obsfile) call status_check(status,error_text) ! Sanity checks on the variable status=nf90_inquire_variable(obsfid,OBSvarid,ndims=nbdim,dimids=OBSvarshape) error_text="Failed to obtain information on variable "//trim(OBSvarname) call status_check(status,error_text) ! Check that it is a 4D variable if (nbdim.ne.4) then write(*,*) "Error, expected a 4D (lon-lat-alt-time) variable for ",trim(OBSvarname) stop endif ! Check that its dimensions are indeed lon,lat,alt,time (in the right order) if (OBSvarshape(1).ne.lon_dimid_obs) then write(*,*) "Error, expected first dimension of ",trim(OBSvarname)," to be longitude!" stop endif if (OBSvarshape(2).ne.lat_dimid_obs) then write(*,*) "Error, expected second dimension of ",trim(OBSvarname)," to be latitude!" stop endif if (OBSvarshape(3).ne.alt_dimid_obs) then write(*,*) "Error, expected third dimension of ",trim(OBSvarname)," to be altitude!" stop endif if (OBSvarshape(4).ne.time_dimid_obs) then write(*,*) "Error, expected fourth dimension of ",trim(OBSvarname)," to be time!" stop endif ! Length allocation for each dimension of the 4D variable allocate(OBS_var(OBSlonlen,OBSlatlen,OBSaltlen,OBSLslen)) ! Load datasets status=NF90_GET_VAR(obsfid,OBSvarid,OBS_var) error_text="Failed to load "//trim(OBSvarname)//" from the obsfile" call status_check(status,error_text) write(*,*) trim(OBSvarname)," loaded from the obsfile as reference variable" ! Get OBS_var missing_value attribute status=nf90_get_att(obsfid,OBSvarid,"_FillValue",OBSmiss_val) error_text="Failed to load missing_value attribute" call status_check(status,error_text) !================================================================ ! 2.5.2 Number of values in the OBS bin (for the sol binning) !================================================================ ! Check that the observation file contains that variable status=nf90_inq_varid(obsfid,trim(numbin_name),numbin_id) ! Checks have already been done in section 2.3 ! Length allocation for each dimension of the 4D variable allocate(numbin(OBSlonlen,OBSlatlen,OBSaltlen,OBSLslen)) ! Load datasets status=NF90_GET_VAR(obsfid,numbin_id,numbin) !=============================================================================== ! 2.6 Definition of GCM variables in outfile !=============================================================================== ! Switch to netcdf define mode status=nf90_redef(outfid) error_text="Error: could not switch to define mode in the outfile" call status_check(status,error_text) ! Definition of the variable SELECT CASE (dayornight) CASE ("dayside") outvarname = "d"//GCMvarname CASE ("nightside") outvarname = "n"//GCMvarname END SELECT status=NF90_DEF_VAR(outfid,trim(outvarname),nf90_float,OBSvarshape,outvarid) error_text="Error: could not define the variable "//trim(outvarname)//" in the outfile" call status_check(status,error_text) ! Write the attributes SELECT CASE (dayornight) CASE ("dayside") long_name = trim(long_name)//" - day side" status=nf90_put_att(outfid,outvarid,"long_name",long_name) CASE ("nightside") long_name = trim(long_name)//" - night side" status=nf90_put_att(outfid,outvarid,"long_name",long_name) END SELECT status=nf90_put_att(outfid,outvarid,"units",units) status=nf90_put_att(outfid,outvarid,"_FillValue",OBSmiss_val) comment = "Reference numbin: "//trim(numbin_name) status=nf90_put_att(outfid,outvarid,"comment",comment) write(*,*)trim(outvarname)," has been created in the outfile" write(*,'(" with missing_value attribute : ",1pe12.5)')OBSmiss_val write(*,*)"" ! End the netcdf define mode (and thus enter the "data writing" mode) status=nf90_enddef(outfid) error_text="Error: could not close the define mode of the outfile" call status_check(status,error_text) allocate(outvar(OBSlonlen,OBSlatlen,OBSaltlen,OBSLslen)) !=================================================================================== ! 3. EXTRACTION OF THE VARIABLE !=================================================================================== !=============================================================================== ! 3.1 Do some checks and preparations before the extraction !=============================================================================== Ls: do l=1,OBSLslen ! loop on all observed Solar longitudes Ls_val=OBSLs(l) ! Ls_val=center of the output bin if ((Ls_val.lt.0.).or.(Ls_val.gt.360.)) then write(*,*) "Unexpected value for OBSLs: ",Ls_val stop endif ! Convert the Ls bin into a sol interval on which the binning is done : !---------------------------------------------------------------------- !-> get the index of the maximum value of GCM sol (m_maxsol) that is lower than Ls bin's superior bound (maxsol) call ls2sol(Ls_val+OBSdeltaLs/2.,maxsol) m_maxsol=1 do while (((GCMtime(m_maxsol)+0.5).lt.maxsol).AND.(m_maxsol.le.(GCMtimelen-1))) !! The +0.5 is there to take into account the whole planet (lon=[-180°;180°]) and not just the lon=0° from the GCM m_maxsol=m_maxsol+1 enddo !-> get the index of the minimum value of GCM sol (m_minsol) that is greater than Ls bin's inferior bound (minsol) call ls2sol(Ls_val-OBSdeltaLs/2.,minsol) m_minsol=1 do while (((GCMtime(m_minsol)-0.5).le.minsol).AND.(m_minsol.le.(GCMtimelen-1))) !! Same comment for the -0.5 m_minsol=m_minsol+1 enddo if (m_minsol.gt.m_maxsol) then write(*,*) "No value in gcmfile between sol=",minsol," and sol=",maxsol," (Ls=",Ls_val,"°)" ! Write a missing_value to output outvar(:,:,:,l) = OBSmiss_val CYCLE Ls ! go directly to the next Ls endif ! Get all the integer values of GCM sols that fit in this interval solcount=floor(GCMtime(m_maxsol))-ceiling(GCMtime(m_minsol))+1 ! sols that are not fully in the interval are not counted allocate(int_sol_list(solcount)) ! write(*,*) "GCMminsol=", GCMtime(m_minsol) ! write(*,*)"GCMmaxsol=", GCMtime(m_maxsol) do m=1,solcount int_sol_list(m)=ceiling(GCMtime(m_minsol)) + m-1 enddo ! write(*,*)"int_sol_list=",int_sol_list latitude: do j=1,OBSlatlen ! loop on all observed latitudes lat_val=OBSlat(j) if ((lat_val.lt.-90.).or.(lat_val.gt.90.)) then write(*,*) "Unexpected value for OBSlat: ",lat_val stop endif longitude: do i=1,OBSlonlen ! loop on all observed longitudes lon_val=OBSlon(i) if ((lon_val.lt.-360.).or.(lon_val.gt.360.)) then write(*,*) "Unexpected value for lon_val: ",lon_val stop endif ! We want lon_val in [-180:180] for the subroutine extraction if (lon_val.lt.-180.) lon_val=lon_val+360. if (lon_val.gt.180.) lon_val=lon_val-360. LT_val=OBSLT(i,j,l) ! find the Observer average LT value at bin(lon_val, lat_val, Ls_val) if ((LT_val.lt.0.).or.(LT_val.gt.24.)) then if (LT_val.eq.LTmiss_val) then ! write(*,*) "Missing value in obsfile for LT_val" ! Write a missing_value to output outvar(i,j,:,l) = OBSmiss_val CYCLE longitude ! go directly to the next longitude else write(*,*) "Unexpected value for LT_val: ",LT_val stop endif endif altitude: do k=1,OBSaltlen ! loop on all observed altitudes alt_val=OBSalt(k) if (OBS_var(i,j,k,l).eq.OBSmiss_val) then ! write(*,*) "Missing value in obsfile for ",OBSvarname ! Write a missing_value to output outvar(i,j,k,l) = OBSmiss_val CYCLE altitude ! go directly to the next altitude endif !=============================================================================== ! 3.2 Compute GCM sol date corresponding to Observer Ls (via m_(min/max)sol) ! and LT (via OBSLT(min/max)) !=============================================================================== LTcount=floor(numbin(i,j,k,l)) ! find the Observer number of temp values ! at bin(lon_val,lat_val,alt_val,Ls_val) if (LTcount.eq.0.) then ! Write a missing_value to output outvar(i,j,k,l) = OBSmiss_val CYCLE altitude ! go directly to the next altitude endif if (LTcount.lt.0.) then write(*,*) "Unexpected value for LTcount: ",LTcount stop endif ! Generate the sol list for the interpolation allocate(sol_list(solcount*LTcount)) call gen_sol_list(solcount,int_sol_list,LTcount,LT_val,OBSLTmax(i,j,l),& OBSLTmin(i,j,l),lon_val,LTmod,dayornight,& sol_list) solerrcount=0 solbinned_value=0 sol_bin: do m=1,solcount*LTcount ! loop on all GCM sols of the bin sol=sol_list(m) ! write(*,*)"sol=",sol !=============================================================================== ! 3.3 Do the interpolation and binning for the given location !=============================================================================== call extraction(lon_val,lat_val,alt_val,sol,& GCMlonlen,GCMlatlen,GCMaltlen,GCMtimelen,& GCMlon,GCMlat,GCMalt,GCMtime,& GCM_var,GCMmiss_val,GCMalttype,GCMvarname,extr_value) if (extr_value.eq.GCMmiss_val) then ! write(*,*) "Missing value in gcmfile at lon=",lon_val,"; lat=",lat_val,"; alt=",alt_val,"; sol=",sol solerrcount=solerrcount+1 CYCLE sol_bin ! go directly to the next GCM sol of the bin endif solbinned_value=solbinned_value+extr_value enddo sol_bin ! end loop on all GCM sols of the bin if ((solcount*LTcount-solerrcount).ne.0) then solbinned_value=solbinned_value/(solcount*LTcount-solerrcount) else ! write(*,*)"No GCM value in this sol bin" solbinned_value=OBSmiss_val endif ! Write value to output outvar(i,j,k,l)=solbinned_value errcount=errcount+solerrcount deallocate(sol_list) enddo altitude ! end loop on observed altitudes enddo longitude ! end loop on observed longitudes enddo latitude ! end loop on observed latitudes deallocate(int_sol_list) enddo Ls ! end loop on observed Solar longitudes ! write(*,*)"Nb of GCM missing values :",errcount !=============================================================================== ! 3.4 Write the data in the netcdf output file !=============================================================================== status = nf90_put_var(outfid, outvarid, outvar) error_text="Error: could not write "//trim(outvarname)//" data in the outfile" call status_check(status,error_text) ! Deallocations before going to the next GCM variable deallocate(GCM_var) deallocate(OBS_var) deallocate(numbin) deallocate(outvar) ENDDO VAR ! end loop on variables !=================================================================================== ! 4. CDOD RATIO (tau_GCM/tau_MCS) !=================================================================================== !=============================================================================== ! 4.1 Preliminary checks !=============================================================================== if (OBSalttype.eq.'p') then tau_ratio_ok = .true. else tau_ratio_ok = .false. endif IF (dayornight.EQ."dayside") THEN out_dztauname = "d" OBS_dztauname = "ddust" OBS_tempname = "dtemp" ELSE ! i.e. dayornight="nightside" out_dztauname = "n" OBS_dztauname = "ndust" OBS_tempname = "ntemp" ENDIF ! Check that the output file contains dust opacity ! NB: a dust opacity in outfile requires the dust ! opacity variable to be present in the obsfile, ! so we don't have to check it out if (tau_ratio_ok) then out_dztauname = out_dztauname(1:1)//"dust" status=nf90_inq_varid(outfid,trim(out_dztauname),out_dztauid) if (status.ne.nf90_noerr) then ! if no "d/ndust" in outfile, look for "d/nopadust" out_dztauname = out_dztauname(1:1)//"opadust" status=nf90_inq_varid(outfid,trim(out_dztauname),out_dztauid) if (status.ne.nf90_noerr) then tau_ratio_ok = .false. write(*,*)"No dust opacity in the outfile, we'll pass the computation of the CDOD ratio" endif endif endif ! Check that the obsfile contains temperature (for rho computation) if (tau_ratio_ok) then status=nf90_inq_varid(obsfid,trim(OBS_tempname),OBS_tempid) if (status.ne.nf90_noerr) then tau_ratio_ok = .false. write(*,*)"No temperature in the outfile, we'll pass the computation of the CDOD ratio" endif endif !=============================================================================== ! 4.2 tau_ratio computation !=============================================================================== if (tau_ratio_ok) then write(*,*)"Computing the Column Dust Optical Depths ratio (tau_ratio = tau_GCM/tau_MCS)..." write(*,*)"(it is recommended to use it to normalize the GCM dust profiles when comparing with MCS dust profiles)" ! Creation of the OBS_plev table IF (dayornight.EQ."dayside") THEN ! only do it once allocate(OBS_plev(OBSaltlen+1)) DO k=2,OBSaltlen ! mid-altitude pressure between 2 OBS pressure layers (OBSalt) ! NB: OBS_plev(k) > OBSalt(k) > OBS_plev(k+1) OBS_plev(k) = sqrt(OBSalt(k-1)*OBSalt(k)) ENDDO OBS_plev(1)=2000. OBS_plev(OBSaltlen+1)=0 write(*,*)"OBS_plev = ",OBS_plev ENDIF ! Arrays allocation allocate(out_dztau(OBSlonlen,OBSlatlen,OBSaltlen,OBSLslen)) allocate(OBS_dztau(OBSlonlen,OBSlatlen,OBSaltlen,OBSLslen)) allocate(OBS_temp(OBSlonlen,OBSlatlen,OBSaltlen,OBSLslen)) allocate(tau_ratio(OBSlonlen,OBSlatlen,OBSLslen)) tau_ratio(:,:,:) = 0 ! Load GCM opacity status=NF90_GET_VAR(outfid,out_dztauid,out_dztau) error_text="Failed to load outfile dust opacity" call status_check(status,error_text) ! Load OBS opacity status=nf90_inq_varid(obsfid,trim(OBS_dztauname),OBS_dztauid) status=NF90_GET_VAR(OBSfid,OBS_dztauid,OBS_dztau) error_text="Failed to load obsfile dust opacity" call status_check(status,error_text) ! Load OBS temperature status=NF90_GET_VAR(OBSfid,OBS_tempid,OBS_temp) error_text="Failed to load obsfile temperature" call status_check(status,error_text) do l=1,OBSLslen do j=1,OBSlatlen do i=1,OBSlonlen ! initialize both columns to 0 out_tau = 0 OBS_tau = 0 do k=1,OBSlatlen ! Compute optical depths where both GCM and OBS have non-missing values ! (especially above the highest (in altitude) of the lowest (in altitude) valid values in the GCM and Observer profile) IF ((out_dztau(i,j,k,l).ne.OBSmiss_val).and.(OBS_dztau(i,j,k,l).ne.OBSmiss_val).and.(OBS_temp(i,j,k,l).ne.OBSmiss_val)) then OBS_rho = OBSalt(k) / (r_atm*OBS_temp(i,j,k,l)) OBS_tau = OBS_tau + OBS_dztau(i,j,k,l)*1.e-3 / OBS_rho /g *(OBS_plev(k)-OBS_plev(k+1)) out_tau = out_tau + out_dztau(i,j,k,l)*1.e-3 / OBS_rho /g *(OBS_plev(k)-OBS_plev(k+1)) ENDIF enddo !alt ! Compute the Column Dust Optical Depth ratio if ((out_tau.eq.0).or.(OBS_tau.eq.0)) then tau_ratio(i,j,l) = LTmiss_val else tau_ratio(i,j,l) = out_tau/OBS_tau endif enddo !lon enddo !lat enddo !Ls !=============================================================================== ! 4.3 Write tau_ratio in the netcdf output file !=============================================================================== ! Switch to netcdf define mode status=nf90_redef(outfid) error_text="Error: could not switch to define mode in the outfile" call status_check(status,error_text) ! Definition of the variable SELECT CASE (dayornight) CASE ("dayside") outvarname = "dtau_ratio" CASE ("nightside") outvarname = "ntau_ratio" END SELECT status=NF90_DEF_VAR(outfid,trim(outvarname),nf90_float,LTshape,outvarid) error_text="Error: could not define the variable "//trim(outvarname)//" in the outfile" call status_check(status,error_text) ! Write the attributes long_name = "Integrated optical depths ratio tau_GCM/tau_MCS" SELECT CASE (dayornight) CASE ("dayside") long_name = trim(long_name)//" - day side" status=nf90_put_att(outfid,outvarid,"long_name",long_name) CASE ("nightside") long_name = trim(long_name)//" - night side" status=nf90_put_att(outfid,outvarid,"long_name",long_name) END SELECT status=nf90_put_att(outfid,outvarid,"units","/") status=nf90_put_att(outfid,outvarid,"_FillValue",LTmiss_val) comment = "Reference numbin: dust+temp" status=nf90_put_att(outfid,outvarid,"comment",comment) write(*,*)trim(outvarname)," has been created in the outfile" write(*,'(" with missing_value attribute : ",1pe12.5)')LTmiss_val write(*,*)"" ! End the netcdf define mode (and thus enter the "data writing" mode) status=nf90_enddef(outfid) error_text="Error: could not close the define mode of the outfile" call status_check(status,error_text) ! Write data status = nf90_put_var(outfid, outvarid, tau_ratio) error_text="Error: could not write "//trim(outvarname)//" data in the outfile" call status_check(status,error_text) ! Deallocations deallocate(out_dztau) deallocate(OBS_dztau) deallocate(OBS_temp) deallocate(tau_ratio) endif !tau_ratio_ok !=================================================================================== ! 5. END OF THE DAY/NIGHT LOOP !=================================================================================== IF (dayornight.EQ."dayside") THEN ! this is the end of the first loop (on daytime values) ! and we still have to loop on nighttime values dayornight="nightside" deallocate(OBSLT) deallocate(OBSLTmax) deallocate(OBSLTmin) write(*,*)"" write(*,*) "Beginning the 2nd loop, on nighttime values"; write(*,*)"" CYCLE DAY_OR_NIGHT ELSE ! i.e. dayornight="nightside" ! this is the end of the second loop (on nighttime values) ! and thus the end of the program EXIT DAY_OR_NIGHT ENDIF ENDDO DAY_OR_NIGHT ! end of the day/night loop that begins in section 1.3 !=================================================================================== ! 6. CLOSE THE FILES !=================================================================================== status = nf90_close(gcmfid) error_text="Error: could not close file "//trim(gcmfile) call status_check(status,error_text) status = nf90_close(obsfid) error_text="Error: could not close file "//trim(obsfile) call status_check(status,error_text) status = nf90_close(outfid) error_text="Error: could not close file "//trim(outfile) call status_check(status,error_text) write(*,*)"";write(*,*)"-> Program simu_MCS completed!" end program simu_MCS !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !=================================================================================================== ! Subroutines !=================================================================================================== subroutine extraction(lon,lat,alt,sol,& lonlen,latlen,altlen,timelen,& longitude,latitude,altitude,time,& field,missing_value,alttype,varname,value) implicit none !================================================================ ! Arguments: !================================================================ real,intent(in) :: lon ! sought longitude (deg, in [-180:180]) real,intent(in) :: lat ! sought latitude (deg, in [-90:90]) real,intent(in) :: alt ! sought altitude (m or Pa) real,intent(in) :: sol ! sought date (sols) integer,intent(in) :: lonlen integer,intent(in) :: latlen integer,intent(in) :: altlen integer,intent(in) :: timelen real,intent(in) :: longitude(lonlen) real,intent(in) :: latitude(latlen) real,intent(in) :: altitude(altlen) real,intent(in) :: time(timelen) real,intent(in) :: field(lonlen,latlen,altlen,timelen) real,intent(in) :: missing_value ! default value for "no data" character,intent(in) :: alttype ! altitude coord. type:'z' (altitude, m) ! 'p' (pressure, Pa) character(len=*),intent(in) :: varname ! variable name (in GCM file) real,intent(out) :: value !================================================================ ! Local variables: !================================================================ real,save :: prev_lon=-99 ! previous value of 'lon' routine was called with real,save :: prev_lat=-99 ! previous value of 'lat' routine was called with real,save :: prev_alt=-9.e20 ! ! previous value of 'alt' real,save :: prev_sol=-99 ! previous value of 'sol' routine was called with ! encompasing indexes: integer,save :: ilon_inf=-1,ilon_sup=-1 ! along longitude integer,save :: ilat_inf=-1,ilat_sup=-1 ! along latitude integer,save :: ialt_inf=-1,ialt_sup=-1 ! along altitude integer,save :: itim_inf=-1,itim_sup=-1 ! along time ! intermediate interpolated values real :: t_interp(2,2,2) ! after time interpolation real :: zt_interp(2,2) ! after altitude interpolation real :: yzt_interp(2) ! after latitude interpolation real :: coeff ! interpolation coefficient integer :: i ! By default, set value to missing_value value=missing_value !================================================================ ! 1. Find encompassing indexes !================================================================ if (lon.ne.prev_lon) then do i=1,lonlen-1 if (longitude(i).le.lon) then ilon_inf=i else exit endif enddo ilon_sup=ilon_inf+1 endif ! of if (lon.ne.prev_lon) !write(*,*) 'ilon_inf=',ilon_inf," longitude(ilon_inf)=",longitude(ilon_inf) if (lat.ne.prev_lat) then if (latitude(1).gt.latitude(2)) then ! decreasing latitudes, from 90N to (-)90S (LMDZ-like) do i=1,latlen-1 if (latitude(i).ge.lat) then ilat_inf=i else exit endif enddo else ! increasing latitudes, from (-)90S to 90N (DYNAMICO-like) do i=1,latlen-1 if (latitude(i).le.lat) then ilat_inf=i else exit endif enddo endif ilat_sup=ilat_inf+1 endif ! of if (lat.ne.prev_lat) !write(*,*) 'ilat_inf=',ilat_inf," latitude(ilat_inf)=",latitude(ilat_inf) if (alt.ne.prev_alt) then if (alttype.eq.'p') then ! pressures are ordered from max to min !handle special case for alt not in the altitude(1:altlen) interval if ((alt.gt.altitude(1)).or.(alt.lt.altitude(altlen))) then ialt_inf=-1 ialt_sup=-1 ! return to main program (with value=missing_value) ! write(*,*)"Problem in extraction : GCM alt p" return else ! general case do i=1,altlen-1 if (altitude(i).ge.alt) then ialt_inf=i else exit endif enddo ialt_sup=ialt_inf+1 endif ! of if ((alt.gt.altitude(1)).or.(alt.lt.altitude(altlen))) else ! altitudes (m) are ordered from min to max !handle special case for alt not in the altitude(1:altlen) interval if ((alt.lt.altitude(1)).or.(alt.gt.altitude(altlen))) then ialt_inf=-1 ialt_sup=-1 ! return to main program (with value=missing_value) ! write(*,*)"Problem in extraction : GCM alt z" return else ! general case do i=1,altlen-1 if (altitude(i).le.alt) then ialt_inf=i else exit endif enddo ialt_sup=ialt_inf+1 endif ! of if ((alt.lt.altitude(1)).or.(alt.gt.altitude(altlen))) endif ! of if (alttype.eq.'p') endif ! of if (alt.ne.prev_alt) !write(*,*) 'ialt_inf=',ialt_inf," altitude(ialt_inf)=",altitude(ialt_inf) if (sol.ne.prev_sol) then !handle special case for sol not in the time(1):time(timenlen) interval if ((sol.lt.time(1)).or.(sol.gt.time(timelen))) then itim_inf=-1 itim_sup=-1 ! return to main program (with value=missing_value) ! write(*,*)"Problem in extraction : GCM sol" return else ! general case do i=1,timelen-1 if (time(i).le.sol) then itim_inf=i else exit endif enddo itim_sup=itim_inf+1 endif ! of if ((sol.lt.time(1)).or.(sol.gt.time(timenlen))) endif ! of if (sol.ne.prev_sol) !write(*,*) 'itim_inf=',itim_inf," time(itim_inf)=",time(itim_inf) !write(*,*) 'itim_sup=',itim_sup," time(itim_sup)=",time(itim_sup) !================================================================ ! 2. Interpolate !================================================================ ! check that there are no "missing_value" in the field() elements we need ! otherwise return to main program (with value=missing_value) if (field(ilon_inf,ilat_inf,ialt_inf,itim_inf).eq.missing_value) then ! write(*,*)"Error 1 in extraction" return endif if (field(ilon_inf,ilat_inf,ialt_inf,itim_sup).eq.missing_value) then ! write(*,*)"Error 2 in extraction" return endif if (field(ilon_inf,ilat_inf,ialt_sup,itim_inf).eq.missing_value) then ! write(*,*)"Error 3 in extraction" return endif if (field(ilon_inf,ilat_inf,ialt_sup,itim_sup).eq.missing_value) then ! write(*,*)"Error 4 in extraction" return endif if (field(ilon_inf,ilat_sup,ialt_inf,itim_inf).eq.missing_value) then ! write(*,*)"Error 5 in extraction" return endif if (field(ilon_inf,ilat_sup,ialt_inf,itim_sup).eq.missing_value) then ! write(*,*)"Error 6 in extraction" return endif if (field(ilon_inf,ilat_sup,ialt_sup,itim_inf).eq.missing_value) then ! write(*,*)"Error 7 in extraction" return endif if (field(ilon_inf,ilat_sup,ialt_sup,itim_sup).eq.missing_value) then ! write(*,*)"Error 8 in extraction" return endif if (field(ilon_sup,ilat_inf,ialt_inf,itim_inf).eq.missing_value) then ! write(*,*)"Error 9 in extraction" return endif if (field(ilon_sup,ilat_inf,ialt_inf,itim_sup).eq.missing_value) then ! write(*,*)"Error 10 in extraction" return endif if (field(ilon_sup,ilat_inf,ialt_sup,itim_inf).eq.missing_value) then ! write(*,*)"Error 11 in extraction" return endif if (field(ilon_sup,ilat_inf,ialt_sup,itim_sup).eq.missing_value) then ! write(*,*)"Error 12 in extraction" return endif if (field(ilon_sup,ilat_sup,ialt_inf,itim_inf).eq.missing_value) then ! write(*,*)"Error 13 in extraction" return endif if (field(ilon_sup,ilat_sup,ialt_inf,itim_sup).eq.missing_value) then ! write(*,*)"Error 14 in extraction" return endif if (field(ilon_sup,ilat_sup,ialt_sup,itim_inf).eq.missing_value) then ! write(*,*)"Error 15 in extraction" return endif if (field(ilon_sup,ilat_sup,ialt_sup,itim_sup).eq.missing_value) then ! write(*,*)"Error 16 in extraction" return endif !================================================================ ! 2.1 Linear interpolation in time !================================================================ coeff=(sol-time(itim_inf))/(time(itim_sup)-time(itim_inf)) t_interp(1,1,1)=field(ilon_inf,ilat_inf,ialt_inf,itim_inf)+ & coeff*(field(ilon_inf,ilat_inf,ialt_inf,itim_sup)- & field(ilon_inf,ilat_inf,ialt_inf,itim_inf)) t_interp(1,1,2)=field(ilon_inf,ilat_inf,ialt_sup,itim_inf)+ & coeff*(field(ilon_inf,ilat_inf,ialt_sup,itim_sup)- & field(ilon_inf,ilat_inf,ialt_sup,itim_inf)) t_interp(1,2,1)=field(ilon_inf,ilat_sup,ialt_inf,itim_inf)+ & coeff*(field(ilon_inf,ilat_sup,ialt_inf,itim_sup)- & field(ilon_inf,ilat_sup,ialt_inf,itim_inf)) t_interp(1,2,2)=field(ilon_inf,ilat_sup,ialt_sup,itim_inf)+ & coeff*(field(ilon_inf,ilat_sup,ialt_sup,itim_sup)- & field(ilon_inf,ilat_sup,ialt_sup,itim_inf)) t_interp(2,1,1)=field(ilon_sup,ilat_inf,ialt_inf,itim_inf)+ & coeff*(field(ilon_sup,ilat_inf,ialt_inf,itim_sup)- & field(ilon_sup,ilat_inf,ialt_inf,itim_inf)) t_interp(2,1,2)=field(ilon_sup,ilat_inf,ialt_sup,itim_inf)+ & coeff*(field(ilon_sup,ilat_inf,ialt_sup,itim_sup)- & field(ilon_sup,ilat_inf,ialt_sup,itim_inf)) t_interp(2,2,1)=field(ilon_sup,ilat_sup,ialt_inf,itim_inf)+ & coeff*(field(ilon_sup,ilat_sup,ialt_inf,itim_sup)- & field(ilon_sup,ilat_sup,ialt_inf,itim_inf)) t_interp(2,2,2)=field(ilon_sup,ilat_sup,ialt_sup,itim_inf)+ & coeff*(field(ilon_sup,ilat_sup,ialt_sup,itim_sup)- & field(ilon_sup,ilat_sup,ialt_sup,itim_inf)) !================================================================ ! 2.2 Vertical interpolation !================================================================ if (((varname=='rho').or.(varname=='pressure')).and.(alttype=='z')) then ! do the interpolation on the log of the quantity coeff=(alt-altitude(ialt_inf))/(altitude(ialt_sup)-altitude(ialt_inf)) zt_interp(1,1)=log(t_interp(1,1,1))+coeff* & (log(t_interp(1,1,2))-log(t_interp(1,1,1))) zt_interp(1,2)=log(t_interp(1,2,1))+coeff* & (log(t_interp(1,2,2))-log(t_interp(1,2,1))) zt_interp(2,1)=log(t_interp(2,1,1))+coeff* & (log(t_interp(2,1,2))-log(t_interp(2,1,1))) zt_interp(2,2)=log(t_interp(2,2,1))+coeff* & (log(t_interp(2,2,2))-log(t_interp(2,2,1))) zt_interp(1:2,1:2)=exp(zt_interp(1:2,1:2)) else ! general case coeff=(alt-altitude(ialt_inf))/(altitude(ialt_sup)-altitude(ialt_inf)) zt_interp(1,1)=t_interp(1,1,1)+coeff*(t_interp(1,1,2)-t_interp(1,1,1)) zt_interp(1,2)=t_interp(1,2,1)+coeff*(t_interp(1,2,2)-t_interp(1,2,1)) zt_interp(2,1)=t_interp(2,1,1)+coeff*(t_interp(2,1,2)-t_interp(2,1,1)) zt_interp(2,2)=t_interp(2,2,1)+coeff*(t_interp(2,2,2)-t_interp(2,2,1)) endif !================================================================ ! 2.3 Latitudinal interpolation !================================================================ coeff=(lat-latitude(ilat_inf))/(latitude(ilat_sup)-latitude(ilat_inf)) yzt_interp(1)=zt_interp(1,1)+coeff*(zt_interp(1,2)-zt_interp(1,1)) yzt_interp(2)=zt_interp(2,1)+coeff*(zt_interp(2,2)-zt_interp(2,1)) !================================================================ ! 2.4 longitudinal interpolation !================================================================ coeff=(lon-longitude(ilon_inf))/(longitude(ilon_sup)-longitude(ilon_inf)) value=yzt_interp(1)+coeff*(yzt_interp(2)-yzt_interp(1)) end subroutine extraction !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine inidim(outfid,lonlen,latlen,altlen,timelen,lon,lat,alt,time,units_alt,& londimid,latdimid,altdimid,timedimid) !================================================================ ! Purpose: ! Initialize a data file (NetCDF format) !================================================================ ! Remarks: ! The NetCDF file remains open !================================================================ use netcdf ! NetCDF definitions implicit none !================================================================ ! Arguments: !================================================================ integer, intent(in):: outfid ! outfid: [netcdf] file ID integer, intent(in):: lonlen ! lonlen: longitude length (# of longitude values) integer, intent(in):: latlen ! latlen: latitude length (# of latitude values) integer, intent(in):: altlen ! altlen: altitude length (# of altitude values) integer, intent(in):: timelen ! timelen: time length (# of time values) real, intent(in):: lon(lonlen) ! lon(): longitude real, intent(in):: lat(latlen) ! lat(): latitude real, intent(in):: alt(altlen) ! alt(): altitude real, intent(in):: time(timelen) ! time(): time (Ls) character(len=1), intent(in) :: units_alt ! units_alt: altitude coord. type:'z' (altitude, m) 'p' (pressure, Pa) integer,intent(inout):: londimid ! londimid: [netcdf] lon() (i.e.: longitude) ID in MCS and output files (they are the same) integer,intent(inout) :: latdimid ! latdimid: [netcdf] lat() ID in MCS and output files (they are the same) integer,intent(inout):: altdimid ! altdimid: [netcdf] alt() ID in MCS and output files (they are the same) integer,intent(inout):: timedimid ! timedimid: [netcdf] time() ID in MCS and output files (they are the same) !================================================================ ! Local variables: !================================================================ integer :: lonvarid,latvarid,altvarid,timevarid,status ! *varid: [netcdf] ID of a variable ! status: [netcdf] return error code (from called subroutines) character(len=256) :: error_text integer :: d ! loop index on dimensions ID !=============================================================== ! 1. Define/write "dimensions" in the same order than MCS file ! and get their IDs !================================================================ do d=1,4 if (altdimid.eq.d) then status=nf90_def_dim(outfid, "altitude", altlen, altdimid) error_text="Error: could not define the altitude dimension in the outfile" call status_check(status,error_text) else if (timedimid.eq.d) then status=nf90_def_dim(outfid, "time", timelen, timedimid) error_text="Error: could not define the time dimension in the outfile" call status_check(status,error_text) else if (latdimid.eq.d) then status=nf90_def_dim(outfid, "latitude", latlen, latdimid) error_text="Error: could not define the latitude dimension in the outfile" call status_check(status,error_text) else if (londimid.eq.d) then status=nf90_def_dim(outfid, "longitude", lonlen, londimid) error_text="Error: could not define the longitude dimension in the outfile" call status_check(status,error_text) endif enddo !================================================================ ! 2. Define "variables" and their attributes !================================================================ !================================================================ ! 2.1 Write "longitude" (data and attributes) !================================================================ ! Insert the definition of the variable status=nf90_def_var(outfid,"longitude",nf90_float,(/londimid/),lonvarid) error_text="Error: could not define the longitude variable in the outfile" call status_check(status,error_text) ! Write the attributes status=nf90_put_att(outfid,lonvarid,"long_name","longitude") error_text="Error: could not put the long_name attribute for the longitude variable in the outfile" call status_check(status,error_text) status=nf90_put_att(outfid,lonvarid,"units","degrees_east") error_text="Error: could not put the units attribute for the longitude variable in the outfile" call status_check(status,error_text) !================================================================ ! 2.2 "latitude" !================================================================ ! Insert the definition of the variable status=nf90_def_var(outfid,"latitude",nf90_float,(/latdimid/),latvarid) error_text="Error: could not define the latitude variable in the outfile" call status_check(status,error_text) ! Write the attributes status=nf90_put_att(outfid,latvarid,"long_name","latitude") error_text="Error: could not put the long_name attribute for the latitude variable in the outfile" call status_check(status,error_text) status=nf90_put_att(outfid,latvarid,"units","degrees_north") error_text="Error: could not put the units attribute for the latitude variable in the outfile" call status_check(status,error_text) !================================================================ ! 2.3 Write "altitude" (data and attributes) !================================================================ ! Insert the definition of the variable status=nf90_def_var(outfid,"altitude",nf90_float,(/altdimid/),altvarid) error_text="Error: could not define the altitude variable in the outfile" call status_check(status,error_text) ! Write the attributes if (units_alt.eq.'p') then ! pressure coordinate status=nf90_put_att(outfid,altvarid,"long_name","pressure") error_text="Error: could not put the long_name attribute for the altitude variable in the outfile" call status_check(status,error_text) status=nf90_put_att(outfid,altvarid,'units',"Pa") error_text="Error: could not put the units attribute for the altitude variable in the outfile" call status_check(status,error_text) status=nf90_put_att(outfid,altvarid,'positive',"down") error_text="Error: could not put the positive attribute for the altitude variable in the outfile" call status_check(status,error_text) status=nf90_put_att(outfid,altvarid,'comment',& "The vertical variable is in fact pressure, not altitude. We just call it altitude for easy reading in Ferret and Grads.") error_text="Error: could not put the comment attribute for the altitude variable in the outfile" call status_check(status,error_text) else if (units_alt.eq.'z') then ! altitude coordinate status=nf90_put_att(outfid,altvarid,"long_name","altitude") error_text="Error: could not put the long_name attribute for the altitude variable in the outfile" call status_check(status,error_text) status=nf90_put_att(outfid,altvarid,'units',"m") error_text="Error: could not put the units attribute for the altitude variable in the outfile" call status_check(status,error_text) status=nf90_put_att(outfid,altvarid,'positive',"up") error_text="Error: could not put the positive attribute for the altitude variable in the outfile" call status_check(status,error_text) else write(*,*)"I do not understand this unit type ",units_alt," for outfile altitude!" stop end if !================================================================ ! 2.4 "time" !================================================================ ! Insert the definition of the variable status=nf90_def_var(outfid,"time",nf90_float,(/timedimid/),timevarid) error_text="Error: could not define the time variable in the outfile" call status_check(status,error_text) ! Write the attributes status=nf90_put_att(outfid,timevarid,"long_name","Solar longitude") error_text="Error: could not put the long_name attribute for the time variable in the outfile" call status_check(status,error_text) status=nf90_put_att(outfid,timevarid,"units","days since 0000-01-1 00:00:00") error_text="Error: could not put the units attribute for the time variable in the outfile" call status_check(status,error_text) status=nf90_put_att(outfid,timevarid,"comment",& "Units is in fact degrees, but set to a dummy value of days for compatibility with Ferret and Grads.") error_text="Error: could not put the comment attribute for the time variable in the outfile" call status_check(status,error_text) !================================================================ ! 2.5 End netcdf define mode !================================================================ status=nf90_enddef(outfid) error_text="Error: could not end the define mode of the outfile" call status_check(status,error_text) !================================================================ ! 3. Write "variables" (data of the dimension variables) !================================================================ ! "time" status=nf90_put_var(outfid,timevarid,time) error_text="Error: could not write the time variable's data in the outfile" call status_check(status,error_text) ! "latitude" status=nf90_put_var(outfid,latvarid,lat) error_text="Error: could not write the latitude variable's data in the outfile" call status_check(status,error_text) ! "longitude" status=nf90_put_var(outfid,lonvarid,lon) error_text="Error: could not write the longitude variable's data in the outfile" call status_check(status,error_text) ! "altitude" status=nf90_put_var(outfid,altvarid,alt) error_text="Error: could not write the altitude variable's data in the outfile" call status_check(status,error_text) end subroutine inidim !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine ls2sol(ls,sol) implicit none !================================================================ ! Arguments: !================================================================ real,intent(in) :: ls real,intent(out) :: sol !================================================================ ! Local: !================================================================ double precision xref,zx0,zteta,zz !xref: mean anomaly, zteta: true anomaly, zx0: eccentric anomaly double precision year_day double precision peri_day,timeperi,e_elips double precision pi,degrad parameter (year_day=668.6d0) ! number of sols in a martian year parameter (peri_day=485.35d0) ! date (in sols) of perihelion !timeperi: 2*pi*( 1 - Ls(perihelion)/ 360 ); Ls(perihelion)=250.99 parameter (timeperi=1.90258341759902d0) parameter (e_elips=0.0934d0) ! eccentricity of orbit parameter (pi=3.14159265358979d0) parameter (degrad=57.2957795130823d0) if (abs(ls).lt.1.0e-5) then if (ls.ge.0.0) then sol = 0.0 else sol = real(year_day) end if return end if zteta = ls/degrad + timeperi zx0 = 2.0*datan(dtan(0.5*zteta)/dsqrt((1.+e_elips)/(1.-e_elips))) xref = zx0-e_elips*dsin(zx0) zz = xref/(2.*pi) sol = real(zz*year_day + peri_day) if (sol.lt.0.0) sol = sol + real(year_day) if (sol.ge.year_day) sol = sol - real(year_day) end subroutine ls2sol !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine gen_sol_list(solcount,int_sol_list,LTcount,LTave,LTmax,LTmin,lon,f_LT,dayornight,& sol_list) !================================================================ ! Purpose: ! Generate a list that is the combination of two lists : ! - int_sol_list, which is a list of integer values of sol ! - LT_list, which contains a number LTcount of LT values in the ! interval [LTmin;LTmax] which are evenly distributed around ! LTave (so that their average is LTave) !================================================================ implicit none !================================================================ ! Arguments: !================================================================ integer, intent(in) :: solcount ! nb of integer values of sol real, intent(in) :: int_sol_list(solcount) ! list of the integer values of sol integer, intent(in) :: LTcount ! nb of LT per sol to be interpolated real, intent(in) :: LTave ! average of LT real, intent(in) :: LTmax, LTmin ! bounds of the LT interval for the interpolation real, intent(in) :: lon ! longitude value for the transformation into a sol value at lon=0° external f_LT ! function that changes LT interval for night LT real f_LT character (len=10), intent(in) :: dayornight ! to know if we have day or night values real, intent(out) :: sol_list(solcount*LTcount) ! all the sol values at lon=0° to interpolate !================================================================ ! Local variables: !================================================================ integer :: N integer :: a,b,c ! loop iteration indexes real :: LT_list(LTcount) N = floor(LTcount/2.) !================================================================ ! 1. Filling of LT_list !================================================================ SELECT CASE (dayornight) CASE ("dayside") if (abs(LTave-LTmax).le.abs(LTave-LTmin)) then ! LTave is closer to LTmax than to LTmin do a=1,N LT_list(a)=LTave+a*abs(LTave-LTmax)/(N+1) LT_list(N+a)=LTave-a*abs(LTave-LTmax)/(N+1) enddo else ! LTave is closer to LTmin than to LTmax do a=1,N LT_list(a)=LTave+a*abs(LTave-LTmin)/(N+1) LT_list(N+a)=LTave-a*abs(LTave-LTmin)/(N+1) enddo endif CASE ("nightside") if (abs(f_LT(LTave)-f_LT(LTmax)).le.abs(f_LT(LTave)-f_LT(LTmin))) then ! LTave is closer to LTmax than to LTmin do a=1,N LT_list(a)=LTave+a*abs(f_LT(LTave)-f_LT(LTmax))/(N+1) LT_list(a)=mod(LT_list(a),24.) ! reput the LT in a [0;24[ interval if needed LT_list(N+a)=LTave-a*abs(f_LT(LTave)-f_LT(LTmax))/(N+1) LT_list(N+a)=mod(LT_list(N+a),24.) enddo else ! LTave is closer to LTmin than to LTmax do a=1,N LT_list(a)=LTave+a*abs(f_LT(LTave)-f_LT(LTmin))/(N+1) LT_list(a)=mod(LT_list(a),24.) LT_list(N+a)=LTave-a*abs(f_LT(LTave)-f_LT(LTmin))/(N+1) LT_list(N+a)=mod(LT_list(N+a),24.) enddo endif END SELECT if (mod(LTcount,2).eq.1) then ! if LTcount is an odd number LT_list(LTcount)=LTave ! add LTave to the list endif !================================================================ ! 2. Combination of int_sol_list & LT_list into sol_list !================================================================ c=1 do a=1,solcount do b=1,LTcount sol_list(c)=int_sol_list(a)+(LT_list(b)-lon/15.)/24. c=c+1 enddo enddo end subroutine gen_sol_list !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine status_check(status,error_text) use netcdf implicit none !================================================================ ! Arguments: !================================================================ integer,intent(in) :: status character(len=256),intent(in) :: error_text if (status.ne.nf90_noerr) then write(*,*)trim(error_text) write(*,*)trim(nf90_strerror(status)) stop endif end subroutine status_check !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! real function LTmod(LT) !================================================================ ! Purpose: ! For night local times management, replace hours ! from a [0;24[ interval to a [-12;12[ interval in which ! midnight = 0 (in order to ensure continuity when comparing ! night local times) !================================================================ implicit none real,intent(in) :: LT LTmod = 2*mod(LT,12.)-LT return end function LTmod