program zrecast ! This program reads 4D (lon-lat-alt-time) fields from GCM output files ! (ie: diagfi.nc time series or concat.nc or stats.nc files) and, by ! integrating the hydrostatic equation, recasts data along the vertical ! direction. ! The vertical coordinate can be either "above areoid altitudes" or ! "pressure". Some interpolation along the vertical direction is also ! done, following instructions given by user. ! For "above areoid altitudes" output, Atmospheric pressure is added to ! output dataset; for "pressure coordinate" outputs, the above areoid ! altitude of pressure is added to output dataset. ! ! Minimal requirements and dependencies: ! The dataset must include the following data: ! - surface pressure ! - atmospheric temperature ! - hybrid coordinates aps() and bps(), or sigma levels() (see section 1.3.2) ! - ground geopotential (in input file; if not found, it is sought ! in a 'diagfi.nc' file. If not found there, it is then sought in ! a 'phisinit.nc' file (see section 1.3.3 of program) ! ! - When integration the hydrostatic equation, we assume that R, the molecular ! Gas Constant, may not be constant, so it is computed as ! R=P/(rho*T) (P=Pressure, rho=density, T=temperature) ! If 'rho' is not available, then we use a constant R (see section 2.2) ! ! WARNING: Asking for many points along the vertical direction quickly ! leads to HUGE output files. ! ! EM 01/2006 : Corrected a bug in vertical (log) interpolation for pressure ! and density ! EM 10/2006 : Modified program so that it can now process 'stats.nc' ! files obtained from British GCM (ie: vertical coordinate ! given as sigma levels and geopotential read from file ! 'phisinit.nc') ! EM 02/2007 : Changed behavior for "altitude above surface" case ! (for MCD RMS computation). Number of levels is then set as ! number of levels in initial file, ! and the new set of above surface levels follow a more elaborate ! distribution (see build_zs.F90 routine). ! EM 08/2009 : User may now specify values of each vertical level, ! or just min,max and number of levels (as before) ! EM 01/2010 : Corrected bug in 'zs_coord_interp' to correctly handle the case ! when interpolating log-wise and the density field is not ! available. implicit none include "netcdf.inc" ! NetCDF definitions character (len=128) :: infile ! input file name (diagfi.nc or stats.nc format) character (len=128) :: infile2 ! second input file (may be needed for 'phisini') character (len=128) :: outfile ! output file name character (len=64) :: text ! to store some text character (len=64) :: tmpvarname ! temporarily store a variable name integer tmpvarid ! temporarily store a variable ID integer tmpdimid ! temporarily store a dimension ID integer tmpndims ! temporarily store # of dimensions of a variable integer infid ! NetCDF input file ID (of diagfi.nc or stats.nc format) integer infid2 ! NetCDF input file which contains 'phisini' dataset (diagfi.nc) integer nbvarinfile ! # of variables in input file integer nbattr ! # of attributes of a given variable in input file integer nbvar4dinfile ! # of 4D (lon,lat,alt,time) variables in input file integer outfid ! NetCDF output file ID integer lon_dimid,lat_dimid,alt_dimid,time_dimid ! NetCDF dimension IDs integer lon_varid,lat_varid,alt_varid,time_varid integer gcm_layers_dimid ! NetCDF dimension ID for # of layers in GCM integer sigma_varid,aps_varid,bps_varid integer za_varid,p_varid ! above areoid and pressure data IDs integer ps_varid ! surface pressure data ID integer,dimension(4) :: datashape ! shape of 4D datasets integer,dimension(3) :: surfdatashape ! shape of 3D (surface+time) datasets real :: miss_val=-9.99e+33 ! special "missing value" to specify missing data real,parameter :: miss_val_def=-9.99e+33 ! default value for "missing value" character (len=64), dimension(:), allocatable :: var ! var(): names of variables that will be processed integer nbvar ! # of variables to process integer,dimension(:),allocatable :: var_id ! IDs of variables var() (in outfile) real,dimension(:),allocatable :: lon ! longitude integer lonlength ! # of grid points along longitude real,dimension(:),allocatable :: lat ! latitude integer latlength ! # of grid points along latitude integer altlength ! # of grid point along altitude (of input datasets) real,dimension(:),allocatable :: time ! time integer timelength ! # of points along time real,dimension(:),allocatable :: aps,bps ! hybrid vertical coordinates real,dimension(:),allocatable :: sigma ! sigma levels real,dimension(:,:),allocatable :: phisinit ! Ground geopotential real,dimension(:,:,:),allocatable :: ps ! GCM surface pressure real,dimension(:,:,:,:),allocatable :: press ! GCM atmospheric pressure real,dimension(:,:,:,:),allocatable :: temp ! GCM atmospheric temperature real,dimension(:,:,:,:),allocatable :: rho ! GCM atmospheric density real,dimension(:,:,:,:),allocatable :: za_gcm ! GCM above areoid levels (m) real,dimension(:,:,:,:),allocatable :: zs_gcm ! GCM above surface heights (m) real,dimension(:,:,:,:),allocatable :: indata ! to store a GCM 4D dataset real,dimension(:,:,:,:),allocatable :: outdata ! to store a 4D dataset integer ierr ! NetCDF routines return code integer i,j,ilon,ilat,ilev,itim ! for loops integer ztype ! Flag for vertical coordinate of output ! ztype=1: pressure ztype=2: above areoid ztype=3: above local surface integer nblev ! # of levels (along vertical coordinate) for output data real pmin,pmax ! min and max values for output pressure coordinate real,dimension(:),allocatable :: plevel ! Pressure levels for output data real zamin,zamax ! min and max values for output above areoid coordinate real,dimension(:),allocatable :: zareoid ! Above areoid heights for output data real,dimension(:),allocatable :: zsurface ! Above surface heights for output logical :: have_rho ! Flag: true if density 'rho' is available logical :: have_sigma ! Flag: true if sigma levels are known (false if hybrid ! coordinates are used) logical :: auto_vert_levels ! Flag: true if the positions of vertical levels ! has to be computed; false if these are given ! by the user !=============================================================================== ! 1. Input parameters !=============================================================================== !=============================================================================== ! 1.1 Input file !=============================================================================== write(*,*) "" write(*,*) " Program valid for diagfi.nc, concatnc.nc and stats.nc files" write(*,*) "Enter input file name:" read(*,'(a128)') infile write(*,*) "" ! open input file ierr = NF_OPEN(infile,NF_NOWRITE,infid) if (ierr.ne.NF_NOERR) then write(*,*) 'ERROR: Pb opening file ',trim(infile) stop "" endif !=============================================================================== ! 1.2 Get # and names of variables in input file !=============================================================================== ierr=NF_INQ_NVARS(infid,nbvarinfile) if (ierr.ne.NF_NOERR) then write(*,*) 'ERROR: Failed geting number of variables from file' stop endif write(*,*)" The following variables have been found:" nbvar4dinfile=0 do i=1,nbvarinfile ! get name of variable # i ierr=NF_INQ_VARNAME(infid,i,tmpvarname) ! check if it is a 4D variable ierr=NF_INQ_VARNDIMS(infid,i,tmpndims) if (tmpndims.eq.4) then nbvar4dinfile=nbvar4dinfile+1 write(*,*) trim(tmpvarname) endif enddo allocate(var(nbvar4dinfile)) write(*,*) "" write(*,*) "Which variable do you want to keep?" write(*,*) "all or list of (separated by s)" write(*,*) "(an empty line , i.e: just , implies end of list)" nbvar=0 read(*,'(a64)') tmpvarname do while ((tmpvarname.ne.' ').and.(trim(tmpvarname).ne.'all')) ! check if tmpvarname is valid ierr=NF_INQ_VARID(infid,tmpvarname,tmpvarid) if (ierr.eq.NF_NOERR) then ! valid name ! also check that it is indeed a 4D variable ierr=NF_INQ_VARNDIMS(infid,tmpvarid,tmpndims) if (tmpndims.eq.4) then nbvar=nbvar+1 var(nbvar)=tmpvarname else write(*,*) 'Error: ',trim(tmpvarname),' is not a 4D variable' write(*,*) ' (we''ll skip that one)' endif else ! invalid name write(*,*) 'Error: ',trim(tmpvarname),' is not a valid name' write(*,*) ' (we''ll skip that one)' endif read(*,'(a64)') tmpvarname enddo ! handle "all" case if (tmpvarname.eq.'all') then nbvar=0 do i=1,nbvarinfile ! look for 4D variables ierr=NF_INQ_VARNDIMS(infid,i,tmpndims) if (tmpndims.eq.4) then nbvar=nbvar+1 ! get the corresponding name ierr=NF_INQ_VARNAME(infid,i,tmpvarname) var(nbvar)=tmpvarname endif enddo endif ! Check that there is at least 1 variable to process if (nbvar.eq.0) then write(*,*) 'No variables to process !?' write(*,*) 'Might as well stop here' stop "" else write(*,*) "" write(*,*) 'OK, the following variables will be processed:' do i=1,nbvar write(*,*) var(i) enddo endif !=============================================================================== ! 1.3 Get grids in lon,lat,alt,time, ! as well as hybrid coordinates aps() and bps() (or sigma levels sigma()) ! and phisinit() from input file !=============================================================================== ! 1.3.1 longitude, latitude, altitude and time ! latitude ierr=NF_INQ_DIMID(infid,"latitude",tmpdimid) if (ierr.ne.NF_NOERR) then write(*,*) "Could not get latitude dimension ID" write(*,*) " looking for lat dimension instead... " ierr=NF_INQ_DIMID(infid,"lat",tmpdimid) if (ierr.ne.NF_NOERR) then stop "Error: Failed to get lat dimension ID" else ierr=NF_INQ_VARID(infid,"lat",tmpvarid) if (ierr.ne.NF_NOERR) then stop "Error: Failed to get lat ID" else ierr=NF_INQ_DIMLEN(infid,tmpdimid,latlength) if (ierr.ne.NF_NOERR) then stop "Error: Failed to get lat length" else allocate(lat(latlength)) ierr=NF_GET_VAR_REAL(infid,tmpvarid,lat) if (ierr.ne.NF_NOERR) then stop "Error: Failed reading lat" endif endif endif endif else ierr=NF_INQ_VARID(infid,"latitude",tmpvarid) if (ierr.ne.NF_NOERR) then stop "Error: Failed to get latitude ID" else ierr=NF_INQ_DIMLEN(infid,tmpdimid,latlength) if (ierr.ne.NF_NOERR) then stop "Error: Failed to get latitude length" else allocate(lat(latlength)) ierr=NF_GET_VAR_REAL(infid,tmpvarid,lat) if (ierr.ne.NF_NOERR) then stop "Error: Failed reading latitude" endif endif endif endif ! longitude ierr=NF_INQ_DIMID(infid,"longitude",tmpdimid) if (ierr.ne.NF_NOERR) then write(*,*) "Could not get longitude dimension ID" write(*,*) " looking for lon dimension instead... " ierr=NF_INQ_DIMID(infid,"lon",tmpdimid) if (ierr.ne.NF_NOERR) then stop "Error: Failed to get lon dimension ID" else ierr=NF_INQ_VARID(infid,"lon",tmpvarid) if (ierr.ne.NF_NOERR) then stop "Error: Failed to get lon ID" else ierr=NF_INQ_DIMLEN(infid,tmpdimid,lonlength) if (ierr.ne.NF_NOERR) then stop "Error: Failed to get lon length" else allocate(lon(lonlength)) ierr=NF_GET_VAR_REAL(infid,tmpvarid,lon) if (ierr.ne.NF_NOERR) then stop "Error: Failed reading longitude" endif endif endif endif else ierr=NF_INQ_VARID(infid,"longitude",tmpvarid) if (ierr.ne.NF_NOERR) then stop "Error: Failed to get longitude ID" else ierr=NF_INQ_DIMLEN(infid,tmpdimid,lonlength) if (ierr.ne.NF_NOERR) then stop "Error: Failed to get longitude length" else allocate(lon(lonlength)) ierr=NF_GET_VAR_REAL(infid,tmpvarid,lon) if (ierr.ne.NF_NOERR) then stop "Error: Failed reading longitude" endif endif endif endif ! time ierr=NF_INQ_DIMID(infid,"Time",tmpdimid) if (ierr.ne.NF_NOERR) then write(*,*) "Could not get Time dimension ID" write(*,*) " looking for time dimension instead... " ierr=NF_INQ_DIMID(infid,"time",tmpdimid) if (ierr.ne.NF_NOERR) then stop "Error: Failed to get lon dimension ID" else ierr=NF_INQ_VARID(infid,"time",tmpvarid) if (ierr.ne.NF_NOERR) then stop "Error: Failed to get time ID" else ierr=NF_INQ_DIMLEN(infid,tmpdimid,timelength) if (ierr.ne.NF_NOERR) then stop "Error: Failed to get time length" else allocate(time(timelength)) ierr=NF_GET_VAR_REAL(infid,tmpvarid,time) if (ierr.ne.NF_NOERR) then stop "Error: Failed reading time" endif endif endif endif else ierr=NF_INQ_VARID(infid,"Time",tmpvarid) if (ierr.ne.NF_NOERR) then stop "Error: Failed to get Time ID" else ierr=NF_INQ_DIMLEN(infid,tmpdimid,timelength) if (ierr.ne.NF_NOERR) then stop "Error: Failed to get Time length" else allocate(time(timelength)) ierr=NF_GET_VAR_REAL(infid,tmpvarid,time) if (ierr.ne.NF_NOERR) then stop "Error: Failed reading Time" endif endif endif endif ! altlength ierr=NF_INQ_DIMID(infid,"altitude",tmpdimid) if (ierr.ne.NF_NOERR) then write(*,*) "Could not get altitude dimension ID" write(*,*) " looking for sigma dimension instead... " ierr=NF_INQ_DIMID(infid,"sigma",tmpdimid) if (ierr.ne.NF_NOERR) then stop "Error: Failed to get sigma dimension ID" else ierr=NF_INQ_DIMLEN(infid,tmpdimid,altlength) if (ierr.ne.NF_NOERR) then stop "Error: Failed to get altitude length" endif endif else ierr=NF_INQ_DIMLEN(infid,tmpdimid,altlength) if (ierr.ne.NF_NOERR) then stop "Error: Failed to get altitude length" endif endif ! 1.3.2 Get hybrid coordinates (or sigma levels) ! start by looking for sigma levels ierr=NF_INQ_VARID(infid,"sigma",tmpvarid) if (ierr.ne.NF_NOERR) then have_sigma=.false. write(*,*) "Could not find sigma levels... will look for hybrid coordinates" else have_sigma=.true. allocate(sigma(altlength)) ierr=NF_GET_VAR_REAL(infid,tmpvarid,sigma) if (ierr.ne.NF_NOERR) then stop "Error: Failed reading sigma" endif endif ! if no sigma levels, look for hybrid coordinates if (.not.have_sigma) then ! hybrid coordinate aps ierr=NF_INQ_VARID(infid,"aps",tmpvarid) if (ierr.ne.NF_NOERR) then stop "Error: Failed to get aps ID" else allocate(aps(altlength)) ierr=NF_GET_VAR_REAL(infid,tmpvarid,aps) if (ierr.ne.NF_NOERR) then stop "Error: Failed reading aps" endif endif ! hybrid coordinate bps ierr=NF_INQ_VARID(infid,"bps",tmpvarid) if (ierr.ne.NF_NOERR) then stop "Error: Failed to get bps ID" else allocate(bps(altlength)) ierr=NF_GET_VAR_REAL(infid,tmpvarid,bps) if (ierr.ne.NF_NOERR) then stop "Error: Failed reading bps" endif endif endif !of if (.not.have_sigma) ! 1.3.3 ground geopotential phisinit allocate(phisinit(lonlength,latlength),stat=ierr) if (ierr.ne.0) then write(*,*) "Failed allocation of phisinit(lonlength,latlength) !!!" write(*,*) "lonlength=",lonlength," latlength=",latlength endif ! look for 'phisinit' in current file ierr=NF_INQ_VARID(infid,"phisinit",tmpvarid) if (ierr.ne.NF_NOERR) then write(*,*) "Warning: Failed to get phisinit ID from file ",trim(infile) infile2="diagfi.nc" write(*,*) " Trying file ",trim(infile2) ierr=NF_OPEN(infile2,NF_NOWRITE,infid2) if (ierr.ne.NF_NOERR) then write(*,*) "Problem: Could not find/open that file" infile2="phisinit.nc" write(*,*) " Trying file ",trim(infile2) ierr=NF_OPEN(infile2,NF_NOWRITE,infid2) if (ierr.ne.NF_NOERR) then write(*,*) "Error: Could not open that file either" stop "Might as well stop here" endif endif ! Get ID for phisinit ierr=NF_INQ_VARID(infid2,"phisinit",tmpvarid) if (ierr.ne.NF_NOERR) then stop "Error: Failed to get phisinit ID" endif ! Get physinit ierr=NF_GET_VAR_REAL(infid2,tmpvarid,phisinit) if (ierr.ne.NF_NOERR) then stop "Error: Failed reading phisinit" endif ! Close file write(*,*) 'OK, got phisinit' ierr=NF_CLOSE(infid2) else ierr=NF_GET_VAR_REAL(infid,tmpvarid,phisinit) if (ierr.ne.NF_NOERR) then stop "Error: Failed reading phisinit" endif endif !=============================================================================== ! 1.4 Choose and build the new vertical coordinate !=============================================================================== write(*,*) "" write(*,*) "Which vertical coordinate should the output be in?" ztype=0 do while ((ztype.lt.1).or.(ztype.gt.3)) write(*,*) "(1: pressure, 2: above areoid altitude 3: above local surface)" read(*,*)ztype enddo text="dummy" ! dummy initialization auto_vert_levels=.true. ! dummy initialization to get rid of compiler warning do while ((trim(text).ne."yes").and.(trim(text).ne."no")) write(*,*) "" write(*,*) "Automatic generation of vertical levels distribution? (yes/no)" write(*,*) "(yes: you only provide min, max and number of levels)" write(*,*) "(no: you provide values for each level)" read(*,'(a64)') text if (trim(text).eq."yes") then auto_vert_levels=.true. else auto_vert_levels=.false. endif enddo if (auto_vert_levels) then ! ask for # of points and end values for pressure or above areoid cases write(*,*) "" if (ztype.le.2) then write(*,*) "Enter min and max of vertical coordinate (Pa or m)" write(*,*) " (in that order and on the same line)" if (ztype.eq.1) then ! pressure coordinate read(*,*) pmin,pmax else ! above areoid coordinate read(*,*) zamin,zamax endif endif ! Build corresponding vertical coordinates if (ztype.eq.1) then ! pressure coordinate write(*,*) "Number of levels along vertical coordinate?" read(*,*) nblev allocate(plevel(nblev)) if (nblev.eq.1) then ! in case only one level is asked for plevel(nblev)=pmin else do i=1,nblev ! build exponentially spread layers plevel(i)=exp(log(pmax)+(log(pmin)-log(pmax))* & ((real(i)-1.0)/(real(nblev)-1.0))) enddo endif else if (ztype.eq.2) then ! above areoid heights write(*,*) "Number of levels along vertical coordinate?" read(*,*) nblev allocate(zareoid(nblev)) if (nblev.eq.1) then ! in case only one level is asked for zareoid(nblev)=zamin else do i=1,nblev zareoid(i)=zamin+(real(i)-1.0)*((zamax-zamin)/(real(nblev)-1.0)) enddo endif else ! above local surface ! set nblev to # of levels from input data files nblev=altlength allocate(zsurface(nblev)) ! build specific above local surface altitudes call build_zs(nblev,have_sigma,sigma,aps,bps,zsurface) endif else ! auto_vert_levels=.false. ; user provides values ! ask for # of points along the vertical write(*,*) "" write(*,*) "Number of levels along vertical coordinate?" read(*,*) nblev if (ztype.eq.1) then ! pressure coordinate allocate(plevel(nblev)) write(*,*) "Enter Pressure (Pa) of levels, ordered" write(*,*) " from max (near-surface) to min (top of atmosphere)," write(*,*) " (one value per line)" do i=1,nblev read(*,*) plevel(i) enddo else if (ztype.eq.2) then ! above areoid heights allocate(zareoid(nblev)) write(*,*) "Enter altitude (m) above areoid of levels, ordered" write(*,*) " from min to max (one value per line)" do i=1,nblev read(*,*) zareoid(i) enddo else ! above local surface allocate(zsurface(nblev)) write(*,*) "Enter altitude (m) above surface of levels, ordered" write(*,*) " from min to max (one value per line)" do i=1,nblev read(*,*) zsurface(i) enddo endif endif ! of if (auto_vert_levels) !=============================================================================== ! 1.5 Get output file name !=============================================================================== write(*,*) "" !write(*,*) "Enter output file name" !read(*,*) outfile if (ztype.eq.1) then ! pressure coordinate outfile=infile(1:len_trim(infile)-3)//"_P.nc" else if (ztype.eq.2) then ! above areoid coordinate outfile=infile(1:len_trim(infile)-3)//"_A.nc" else ! above local surface outfile=infile(1:len_trim(infile)-3)//"_S.nc" endif write(*,*) "Output file name is: "//trim(outfile) !=============================================================================== ! 2.1 Build/store GCM fields which will be used later !=============================================================================== !=============================================================================== ! 2.1.1 Surface pressure !=============================================================================== ierr=NF_INQ_VARID(infid,"ps",tmpvarid) if (ierr.ne.NF_NOERR) then stop "Error: Failed to get ps ID" else allocate(ps(lonlength,latlength,timelength)) ierr=NF_GET_VAR_REAL(infid,tmpvarid,ps) if (ierr.ne.NF_NOERR) then stop "Error: Failed reading surface pressure" endif endif !=============================================================================== ! 2.1.2 Atmospheric pressure !=============================================================================== allocate(press(lonlength,latlength,altlength,timelength)) if (have_sigma) then ! sigma coordinate do itim=1,timelength do ilev=1,altlength do ilat=1,latlength do ilon=1,lonlength press(ilon,ilat,ilev,itim)=sigma(ilev)*ps(ilon,ilat,itim) enddo enddo enddo enddo else ! hybrid coordinates do itim=1,timelength do ilev=1,altlength do ilat=1,latlength do ilon=1,lonlength press(ilon,ilat,ilev,itim)=aps(ilev)+bps(ilev)*ps(ilon,ilat,itim) enddo enddo enddo enddo endif !=============================================================================== ! 2.1.3 Atmospheric temperature !=============================================================================== allocate(temp(lonlength,latlength,altlength,timelength)) ierr=NF_INQ_VARID(infid,"temp",tmpvarid) if (ierr.ne.NF_NOERR) then ! stop "Error: Failed to get temp ID" ! try "t" for temperature ierr=NF_INQ_VARID(infid,"t",tmpvarid) if (ierr.ne.NF_NOERR) then stop "Error: Failed to get t ID" else ierr=NF_GET_VAR_REAL(infid,tmpvarid,temp) if (ierr.ne.NF_NOERR) then stop "Error: Failed reading atmospheric temperature" endif endif else ierr=NF_GET_VAR_REAL(infid,tmpvarid,temp) if (ierr.ne.NF_NOERR) then stop "Error: Failed reading atmospheric temperature" endif endif !=============================================================================== ! 2.1.4 Atmospheric density !=============================================================================== ierr=NF_INQ_VARID(infid,"rho",tmpvarid) if (ierr.ne.NF_NOERR) then write(*,*) "Warning: Failed to get rho ID" have_rho=.false. else have_rho=.true. allocate(rho(lonlength,latlength,altlength,timelength)) ierr=NF_GET_VAR_REAL(infid,tmpvarid,rho) if (ierr.ne.NF_NOERR) then stop "Error: Failed reading atmospheric density" endif endif !=============================================================================== ! 2.2 Build GCM Above areoid (or above surface) altitudes of GCM nodes !=============================================================================== if (have_rho) then if (ztype.le.2) then ! above areoid altitudes (also needed for ztype=1) allocate(za_gcm(lonlength,latlength,altlength,timelength)) call build_gcm_za(lonlength,latlength,altlength,timelength, & phisinit,ps,press,temp,rho,za_gcm) else ! above local surface altitudes allocate(zs_gcm(lonlength,latlength,altlength,timelength)) call build_gcm_zs(lonlength,latlength,altlength,timelength, & phisinit,ps,press,temp,rho,zs_gcm) endif else write(*,*)"Warning: Using constant R to integrate hydrostatic equation" if (ztype.le.2) then ! above areoid altitudes (also needed for ztype=1) allocate(za_gcm(lonlength,latlength,altlength,timelength)) call crude_gcm_za(lonlength,latlength,altlength,timelength, & phisinit,ps,press,temp,za_gcm) else ! above local surface altitudes allocate(zs_gcm(lonlength,latlength,altlength,timelength)) call crude_gcm_zs(lonlength,latlength,altlength,timelength, & phisinit,ps,press,temp,zs_gcm) endif endif !=============================================================================== ! 3. Create output file and initialize definitions of variables and dimensions !=============================================================================== !=============================================================================== ! 3.1. Output file !=============================================================================== ! Create output file ierr=NF_CREATE(outfile,NF_CLOBBER,outfid) if (ierr.ne.NF_NOERR) then write(*,*)"Error: could not create file ",outfile stop endif !=============================================================================== ! 3.2. Define dimensions !=============================================================================== ! longitude ierr=NF_DEF_DIM(outfid,"longitude",lonlength,lon_dimid) if (ierr.ne.NF_NOERR) then stop "Error: Could not define longitude dimension" endif ! latitude ierr=NF_DEF_DIM(outfid,"latitude",latlength,lat_dimid) if (ierr.ne.NF_NOERR) then stop "Error: Could not define latitude dimension" endif ! altitude ierr=NF_DEF_DIM(outfid,"altitude",nblev,alt_dimid) if (ierr.ne.NF_NOERR) then stop "Error: Could not define altitude dimension" endif ! time ierr=NF_DEF_DIM(outfid,"Time",timelength,time_dimid) if (ierr.ne.NF_NOERR) then stop "Error: Could not define latitude dimension" endif ! GCM layers (for sigma or aps and bps) ierr=NF_DEF_DIM(outfid,"GCM_layers",altlength,gcm_layers_dimid) if (ierr.ne.NF_NOERR) then stop "Error: Could not define GCM_layers dimension" endif !=============================================================================== ! 3.3. Define variables and their attributes !=============================================================================== ! 3.3.1 Define 1D variables ! longitude datashape(1)=lon_dimid ierr=NF_DEF_VAR(outfid,"longitude",NF_REAL,1,datashape(1),lon_varid) if (ierr.ne.NF_NOERR) then stop "Error: Could not define longitude variable" endif ! longitude attributes text='east longitude' ierr=NF_PUT_ATT_TEXT(outfid,lon_varid,'long_name',len_trim(text),text) if (ierr.ne.NF_NOERR) then stop "Error: Problem writing long_name for longitude" endif text='degrees_east' ierr=NF_PUT_ATT_TEXT(outfid,lon_varid,'units',len_trim(text),text) if (ierr.ne.NF_NOERR) then stop "Error: Problem writing units for longitude" endif ! latitude datashape(2)=lat_dimid ierr=NF_DEF_VAR(outfid,"latitude",NF_REAL,1,datashape(2),lat_varid) if (ierr.ne.NF_NOERR) then stop "Error: Could not define latitude variable" endif ! latitude attributes text='north latitude' ierr=NF_PUT_ATT_TEXT(outfid,lat_varid,'long_name',len_trim(text),text) if (ierr.ne.NF_NOERR) then stop "Error: Problem writing long_name for latitude" endif text='degrees_north' ierr=NF_PUT_ATT_TEXT(outfid,lat_varid,'units',len_trim(text),text) if (ierr.ne.NF_NOERR) then stop "Error: Problem writing units for latitude" endif ! altitude datashape(3)=alt_dimid ierr=NF_DEF_VAR(outfid,"altitude",NF_REAL,1,datashape(3),alt_varid) if (ierr.ne.NF_NOERR) then stop "Error: Could not define altitude variable" endif !altitude attributes if (ztype.eq.1) then ! pressure vertical coordinate text='Pressure levels' ierr=NF_PUT_ATT_TEXT(outfid,alt_varid,'long_name',len_trim(text),text) if (ierr.ne.NF_NOERR) then stop "Error: Problem writing long_name for altitude" endif text='Pa' ierr=NF_PUT_ATT_TEXT(outfid,alt_varid,'units',len_trim(text),text) if (ierr.ne.NF_NOERR) then stop "Error: Problem writing units for altitude" endif text='down' ierr=NF_PUT_ATT_TEXT(outfid,alt_varid,'positive',len_trim(text),text) if (ierr.ne.NF_NOERR) then stop "Error: Problem writing positive for altitude" endif else if (ztype.eq.2) then ! above areoid vertical coordinate text='Altitude above areoid' ierr=NF_PUT_ATT_TEXT(outfid,alt_varid,'long_name',len_trim(text),text) if (ierr.ne.NF_NOERR) then stop "Error: Problem writing long_name for altitude" endif text='m' ierr=NF_PUT_ATT_TEXT(outfid,alt_varid,'units',len_trim(text),text) if (ierr.ne.NF_NOERR) then stop "Error: Problem writing units for altitude" endif text='up' ierr=NF_PUT_ATT_TEXT(outfid,alt_varid,'positive',len_trim(text),text) if (ierr.ne.NF_NOERR) then stop "Error: Problem writing positive for altitude" endif else ! above surface vertical coordinate text='Altitude above local surface' ierr=NF_PUT_ATT_TEXT(outfid,alt_varid,'long_name',len_trim(text),text) if (ierr.ne.NF_NOERR) then stop "Error: Problem writing long_name for altitude" endif text='m' ierr=NF_PUT_ATT_TEXT(outfid,alt_varid,'units',len_trim(text),text) if (ierr.ne.NF_NOERR) then stop "Error: Problem writing units for altitude" endif text='up' ierr=NF_PUT_ATT_TEXT(outfid,alt_varid,'positive',len_trim(text),text) if (ierr.ne.NF_NOERR) then stop "Error: Problem writing positive for altitude" endif endif ! of if (have_sigma) ! GCM_layers !ierr=NF_DEF_VAR(outfid,"gcm_layers",NF_REAL,1,,sigma_varid) ! sigma levels or hybrid coordinates if (have_sigma) then ierr=NF_DEF_VAR(outfid,"sigma",NF_REAL,1,gcm_layers_dimid,sigma_varid) if (ierr.ne.NF_NOERR) then stop "Error: Could not define sigma variable" endif else ! hybrid coordinates ierr=NF_DEF_VAR(outfid,"aps",NF_REAL,1,gcm_layers_dimid,aps_varid) if (ierr.ne.NF_NOERR) then stop "Error: Could not define aps variable" endif ierr=NF_DEF_VAR(outfid,"bps",NF_REAL,1,gcm_layers_dimid,bps_varid) if (ierr.ne.NF_NOERR) then stop "Error: Could not define bps variable" endif endif ! sigma levels (or hybrid coordinates) attributes if (have_sigma) then text="sigma levels" ierr=NF_PUT_ATT_TEXT(outfid,sigma_varid,'long_name',len_trim(text),text) if (ierr.ne.NF_NOERR) then stop "Error: Problem writing long_name for sigma" endif else ! hybrid coordinates text="hybrid pressure at midlayers" ierr=NF_PUT_ATT_TEXT(outfid,aps_varid,'long_name',len_trim(text),text) if (ierr.ne.NF_NOERR) then stop "Error: Problem writing long_name for aps" endif text="hybrid sigma at midlayers" ierr=NF_PUT_ATT_TEXT(outfid,bps_varid,'long_name',len_trim(text),text) if (ierr.ne.NF_NOERR) then stop "Error: Problem writing long_name for bps" endif endif ! of if (have_sigma) ! time datashape(4)=time_dimid ierr=NF_DEF_VAR(outfid,"Time",NF_REAL,1,datashape(4),time_varid) if (ierr.ne.NF_NOERR) then stop "Error: Could not define Time variable" endif ! time attributes text='Time' ierr=NF_PUT_ATT_TEXT(outfid,time_varid,'long_name',len_trim(text),text) if (ierr.ne.NF_NOERR) then stop "Error: Problem writing long_name for Time" endif text='days since 0000-01-1 00:00:00' ierr=NF_PUT_ATT_TEXT(outfid,time_varid,'units',len_trim(text),text) if (ierr.ne.NF_NOERR) then stop "Error: Problem writing units for Time" endif ! 3.3.2 Define 3D variables (ie: surface+time variables) ! Surface pressure surfdatashape(1)=lon_dimid surfdatashape(2)=lat_dimid surfdatashape(3)=time_dimid ierr=NF_DEF_VAR(outfid,"ps",NF_REAL,3,surfdatashape,ps_varid) if (ierr.ne.NF_NOERR) then stop "Error: Could not define ps variable" endif ! Surface pressure attributes text='Surface pressure' ierr=NF_PUT_ATT_TEXT(outfid,ps_varid,'long_name',len_trim(text),text) if (ierr.ne.NF_NOERR) then stop "Error: Problem writing long_name for surface pressure" endif text='Pa' ierr=NF_PUT_ATT_TEXT(outfid,ps_varid,'units',len_trim(text),text) if (ierr.ne.NF_NOERR) then stop "Error: Problem writing units for surface pressure" endif ! 3.3.3 Define 4D variables ! add pressure or zareoid if (ztype.eq.1) then ! pressure vertical coordinate ! zareoid dataset ierr=NF_DEF_VAR(outfid,"zareoid",NF_REAL,4,datashape,za_varid) if (ierr.ne.NF_NOERR) then stop "Error: Could not define zareoid variable" endif ! zareoid attributes text='altitude above areoid' ierr=NF_PUT_ATT_TEXT(outfid,za_varid,'long_name',len_trim(text),text) if (ierr.ne.NF_NOERR) then stop "Error: Problem writing long_name for zareoid" endif text='m' ierr=NF_PUT_ATT_TEXT(outfid,za_varid,'units',len_trim(text),text) if (ierr.ne.NF_NOERR) then stop "Error: Problem writing units for zareoid" endif ! zareoid missing value ierr=NF_PUT_ATT_REAL(outfid,za_varid,'missing_value',NF_REAL,1,miss_val) if (ierr.ne.NF_NOERR) then stop "Error: Problem writing missing_value for zareoid" endif else ! above areoid or above local surface vertical coordinate ! pressure dataset ierr=NF_DEF_VAR(outfid,"pressure",NF_REAL,4,datashape,p_varid) if (ierr.ne.NF_NOERR) then stop "Error: Could not define pressure variable" endif ! pressure attributes text='Atmospheric pressure' ierr=NF_PUT_ATT_TEXT(outfid,p_varid,'long_name',len_trim(text),text) if (ierr.ne.NF_NOERR) then stop "Error: Problem writing long_name for pressure" endif text='Pa' ierr=NF_PUT_ATT_TEXT(outfid,p_varid,'units',len_trim(text),text) if (ierr.ne.NF_NOERR) then stop "Error: Problem writing units for pressure" endif ! pressure missing value ierr=NF_PUT_ATT_REAL(outfid,p_varid,'missing_value',NF_REAL,1,miss_val) if (ierr.ne.NF_NOERR) then stop "Error: Problem writing missing_value for pressure" endif endif ! add zs_gcm if (ztype.eq.3) then endif ! variables requested by user allocate(var_id(nbvar)) do i=1,nbvar write(*,*) "" write(*,*) "Creating ",trim(var(i)) ! define the variable ierr=NF_DEF_VAR(outfid,var(i),NF_REAL,4,datashape,var_id(i)) if (ierr.ne.NF_NOERR) then write(*,*) 'Error, could not define variable ',trim(var(i)) stop "" endif ! Get the (input file) ID for the variable and ! the # of attributes associated to that variable ierr=NF_INQ_VARID(infid,var(i),tmpvarid) if (ierr.ne.NF_NOERR) then write(*,*) 'Error, failed to get ID for variable ',trim(var(i)) stop "" endif ierr=NF_INQ_VARNATTS(infid,tmpvarid,nbattr) if (ierr.ne.NF_NOERR) then write(*,*) 'Error, could not get number of attributes for variable ',& trim(var(i)) stop "" endif ! inititialize j == number of attributes written to output j=0 ! look for a "long_name" attribute text=' ' ierr=NF_GET_ATT_TEXT(infid,tmpvarid,'long_name',text) if (ierr.ne.NF_NOERR) then ! no long_name attribute ! try to find an equivalent 'title' attribute text=' ' ierr=NF_GET_ATT_TEXT(infid,tmpvarid,'title',text) if (ierr.eq.NF_NOERR) then ! found 'title' attribute write(*,*) "Found title ",trim(text) j=j+1 ! write it as a 'long_name' attribute ierr=NF_PUT_ATT_TEXT(outfid,var_id(i),'long_name',len_trim(text),text) if (ierr.ne.NF_NOERR) then write(*,*) "Error failed to copy title attribute:",trim(text) stop "" endif else ! no 'title' attribute ! try to find a "Physics_diagnostic" attribute (UK GCM outputs) text=' ' ierr=NF_GET_ATT_TEXT(infid,tmpvarid,'Physics_diagnostic',text) if (ierr.eq.NF_NOERR) then ! found 'Physics_diagnostic' attribute write(*,*) "Found Physics_diagnostic ",trim(text) j=j+1 ! write it as a 'long_name' attribute ierr=NF_PUT_ATT_TEXT(outfid,var_id(i),'long_name',len_trim(text),text) if (ierr.ne.NF_NOERR) then write(*,*) "Error failed to copy Physics_diagnostic attribute:",trim(text) stop endif endif endif else ! found long_name; write it to outfile write(*,*) "Found long_name ",trim(text) j=j+1 ierr=NF_PUT_ATT_TEXT(outfid,var_id(i),'long_name',len_trim(text),text) if (ierr.ne.NF_NOERR) then write(*,*) "Error failed to copy long_name attribute:",trim(text) stop"" endif endif ! look for a "units" attribute text=' ' ierr=NF_GET_ATT_TEXT(infid,tmpvarid,'units',text) if (ierr.eq.NF_NOERR) then ! found 'units' attribute write(*,*) "Found units ",trim(text) j=j+1 ! write it to output ierr=NF_PUT_ATT_TEXT(outfid,var_id(i),'units',len_trim(text),text) if (ierr.ne.NF_NOERR) then write(*,*) "Error failed to copy units attribute:",trim(text) stop"" endif endif ! look for a "missing_value" attribute ierr=NF_GET_ATT_REAL(infid,tmpvarid,"missing_value",miss_val) if (ierr.eq.NF_NOERR) then ! found 'missing_value' attribute write(*,*) "Found missing_value ",miss_val j=j+1 else ! no 'missing_value' attribute, set miss_val to default miss_val=miss_val_def endif ! write the missing_value attribute to output ierr=NF_PUT_ATT_REAL(outfid,var_id(i),'missing_value',NF_REAL,1,miss_val) if (ierr.ne.NF_NOERR) then stop "Error, failed to write missing_value attribute" endif ! warn if some attributes were missed if (j.ne.nbattr) then write(*,*)'Warning, it seems some attributes of variable ',trim(var(i)) write(*,*)"were not transfered to the new file" write(*,*)'nbattr:',nbattr,' j:',j endif enddo ! of do=1,nbvar !=============================================================================== ! 3.4. Write dimensions (and 1D varaiables) !=============================================================================== ! Switch out of NetCDF define mode ierr=NF_ENDDEF(outfid) if (ierr.ne.NF_NOERR) then stop "Error: Could not switch out of define mode" endif ! Write longitude ierr=NF_PUT_VAR_REAL(outfid,lon_varid,lon) if (ierr.ne.NF_NOERR) then stop "Error: Could not write longitude data to output file" endif ! Write latitude ierr=NF_PUT_VAR_REAL(outfid,lat_varid,lat) if (ierr.ne.NF_NOERR) then stop "Error: Could not write latitude data to output file" endif ! Write altitude if (ztype.eq.1) then ! pressure vertical coordinate ierr=NF_PUT_VAR_REAL(outfid,alt_varid,plevel) if (ierr.ne.NF_NOERR) then stop "Error: Could not write altitude data to output file" endif else if (ztype.eq.2) then ! above areoid altitude ierr=NF_PUT_VAR_REAL(outfid,alt_varid,zareoid) if (ierr.ne.NF_NOERR) then stop "Error: Could not write altitude data to output file" endif else ! above local surface ierr=NF_PUT_VAR_REAL(outfid,alt_varid,zsurface) if (ierr.ne.NF_NOERR) then stop "Error: Could not write altitude data to output file" endif endif ! Write sigma levels (or hybrid coordinates) if (have_sigma) then ierr=NF_PUT_VAR_REAL(outfid,sigma_varid,sigma) if (ierr.ne.NF_NOERR) then stop "Error: Could not write sigma data to output file" endif else ! hybrid coordinates ierr=NF_PUT_VAR_REAL(outfid,aps_varid,aps) if (ierr.ne.NF_NOERR) then stop "Error: Could not write aps data to output file" endif ierr=NF_PUT_VAR_REAL(outfid,bps_varid,bps) if (ierr.ne.NF_NOERR) then stop "Error: Could not write bps data to output file" endif endif ! Write time ierr=NF_PUT_VAR_REAL(outfid,time_varid,time) if (ierr.ne.NF_NOERR) then stop "Error: Could not write Time data to output file" endif !=============================================================================== ! 3.5 Write 3D variables !=============================================================================== ! Write surface pressure ierr=NF_PUT_VAR_REAL(outfid,ps_varid,ps) if (ierr.ne.NF_NOERR) then stop "Error: Could not write ps data to output file" endif !=============================================================================== ! 4. Interpolate and write 4D variables !=============================================================================== ! 4.0 Allocations !indata() to store input (on GCM grid) data allocate(indata(lonlength,latlength,altlength,timelength)) ! outdata() to store output (on new vertical grid) data allocate(outdata(lonlength,latlength,nblev,timelength)) ! 4.1 If output is in pressure coordinate if (ztype.eq.1) then do i=1,nbvar ! loop on 4D variable to process ! identify and read a dataset ierr=NF_INQ_VARID(infid,var(i),tmpvarid) if (ierr.ne.NF_NOERR) then write(*,*) 'Error, failed to get ID for variable ',var(i) stop endif ierr=NF_GET_VAR_REAL(infid,tmpvarid,indata) if (ierr.ne.NF_NOERR) then write(*,*) 'Error, failed to load variable ',var(i) stop endif ! interpolate dataset onto new grid call p_coord_interp(lonlength,latlength,altlength,timelength,nblev, & miss_val,ps,press,indata,plevel,outdata) ! write new dataset to output file ierr=NF_PUT_VAR_REAL(outfid,var_id(i),outdata) if (ierr.ne.NF_NOERR) then write(*,*)'Error, Failed to write variable ',var(i) stop endif enddo ! of do i=1,nbvar ! interpolate zareoid onto new grid call p_coord_interp(lonlength,latlength,altlength,timelength,nblev, & miss_val,ps,press,za_gcm,plevel,outdata) ! write result to output file ierr=NF_PUT_VAR_REAL(outfid,za_varid,outdata) if (ierr.ne.NF_NOERR) then stop "Error, Failed to write zareoid to output file" endif endif ! of if (ztype.eq.1) ! 4.2 If output is in above areoid altitude if (ztype.eq.2) then do i=1,nbvar ! loop on 4D variable to process ! identify and read a dataset ierr=NF_INQ_VARID(infid,var(i),tmpvarid) if (ierr.ne.NF_NOERR) then write(*,*) 'Error, failed to get ID for variable ',var(i) stop endif ierr=NF_GET_VAR_REAL(infid,tmpvarid,indata) if (ierr.ne.NF_NOERR) then write(*,*) 'Error, failed to load variable ',var(i) stop endif ! interpolate dataset onto new grid ! check if variable is "rho" (to set flag for interpolation below) if (var(i).eq.'rho') then j=1 else j=0 endif call z_coord_interp(lonlength,latlength,altlength,timelength,nblev, & miss_val,za_gcm,indata,j,zareoid,outdata) ! write new dataset to output file ierr=NF_PUT_VAR_REAL(outfid,var_id(i),outdata) if (ierr.ne.NF_NOERR) then write(*,*)'Error, Failed to write variable ',var(i) stop endif enddo ! of do i=1,nbvar ! interpolate pressure onto new grid write(*,*) "" write(*,*) "Processing pressure" call z_coord_interp(lonlength,latlength,altlength,timelength,nblev, & miss_val,za_gcm,press,1,zareoid,outdata) ! write result to output file ierr=NF_PUT_VAR_REAL(outfid,p_varid,outdata) if (ierr.ne.NF_NOERR) then stop "Error, Failed to write pressure to output file" endif endif ! of if (ztype.eq.2) ! 4.3 If output is in above local surface altitude if (ztype.eq.3) then do i=1,nbvar ! loop on 4D variable to process write(*,*) " " write(*,*) "Processing "//trim(var(i)) ! identify and read a dataset ierr=NF_INQ_VARID(infid,var(i),tmpvarid) if (ierr.ne.NF_NOERR) then write(*,*) 'Error, failed to get ID for variable ',var(i) stop endif ierr=NF_GET_VAR_REAL(infid,tmpvarid,indata) if (ierr.ne.NF_NOERR) then write(*,*) 'Error, failed to load variable ',var(i) stop endif ! interpolate dataset onto new grid ! check if variable is "rho" (to set flag for interpolation below) if (var(i).eq.'rho') then j=1 else j=0 endif call zs_coord_interp(lonlength,latlength,altlength,timelength,nblev, & miss_val,zs_gcm,indata,phisinit,press,temp,rho,& have_rho,j,zsurface,outdata) ! write new dataset to output file ierr=NF_PUT_VAR_REAL(outfid,var_id(i),outdata) if (ierr.ne.NF_NOERR) then write(*,*)'Error, Failed to write variable ',var(i) stop endif enddo ! of do i=1,nbvar ! interpolate pressure onto new grid write(*,*) "" write(*,*) "Processing pressure" call zs_coord_interp(lonlength,latlength,altlength,timelength,nblev, & miss_val,zs_gcm,press,phisinit,press,temp,rho,& have_rho,1,zsurface,outdata) ! write result to output file ierr=NF_PUT_VAR_REAL(outfid,p_varid,outdata) if (ierr.ne.NF_NOERR) then stop "Error, Failed to write pressure to output file" endif endif ! of if (ztype.eq.3) ! 4.4 Close output file ierr=NF_CLOSE(outfid) if (ierr.ne.NF_NOERR) then write(*,*) 'Error, failed to close output file ',outfile endif end program !=============================================================================== subroutine build_gcm_zs(lonlength,latlength,altlength,timelength, & phis,ps,press,temp,rho,zs_gcm) !============================================================================== ! Purpose: Integrate hydrostatic equation in order to build the "above local ! surface altitudes" corresponding to GCM atmospheric levels !============================================================================== implicit none !============================================================================== ! Arguments: !============================================================================== integer,intent(in) :: lonlength ! # of points along longitude integer,intent(in) :: latlength ! # of points along latitude integer,intent(in) :: altlength ! # of points along altitude integer,intent(in) :: timelength ! # of stored time steps real,dimension(lonlength,latlength),intent(in) :: phis ! phis(:,:) is the ground geopotential real,dimension(lonlength,latlength,timelength),intent(in) :: ps ! ps(:,:) is the surface pressure real,dimension(lonlength,latlength,altlength,timelength),intent(in) :: press ! press(:,:,:,:) is atmospheric pressure real,dimension(lonlength,latlength,altlength,timelength),intent(in) :: temp ! temp(:,:,:,:) is atmospheric temperature real,dimension(lonlength,latlength,altlength,timelength),intent(in) :: rho ! rho(:,:,:,:) is atmospheric density real,dimension(lonlength,latlength,altlength,timelength),intent(out) :: zs_gcm ! zs_gcm(:,:,:,:) are the above local surface altitudes of GCM levels !=============================================================================== ! Local variables: !=============================================================================== real,dimension(:),allocatable :: sigma ! GCM sigma levels real,dimension(:,:,:,:),allocatable :: R ! molecular gas constant !real,dimension(:,:,:,:),allocatable :: zlocal ! local above surface altitude real :: Rmean ! mean value of R for a given level real :: Tmean ! "mean" value of temperature for a given level integer iloop,jloop,kloop,tloop ! Parameters needed to integrate hydrostatic equation: real,parameter :: g0=3.7257964 !g0: exact mean gravity at radius=3396.km (Lemoine et al. 2001) real,parameter :: a0=3396.E3 !a0: 'mean' Mars radius=3396.km real :: gz ! gz: gravitational acceleration at a given (zareoid) altitude !=============================================================================== ! 1. Various initialisations !=============================================================================== allocate(sigma(altlength)) allocate(R(lonlength,latlength,altlength,timelength)) !allocate(zlocal(lonlength,latlength,altlength,timelength)) !============================================================================== ! 2. Compute Molecular Gas Constant (J kg-1 K-1) at every grid point ! --later used to integrate the hydrostatic equation-- !============================================================================== do tloop=1,timelength do kloop=1,altlength do jloop=1,latlength do iloop=1,lonlength R(iloop,jloop,kloop,tloop)=press(iloop,jloop,kloop,tloop)/ & (rho(iloop,jloop,kloop,tloop)* & temp(iloop,jloop,kloop,tloop)) enddo enddo enddo enddo !=============================================================================== ! 3. Integrate hydrostatic equation to compute zlocal and za_gcm !=============================================================================== do tloop=1,timelength do jloop=1,latlength do iloop=1,lonlength ! handle case of first altitude level sigma(1)=press(iloop,jloop,1,tloop)/ps(iloop,jloop,tloop) zs_gcm(iloop,jloop,1,tloop)=-log(sigma(1))*R(iloop,jloop,1,tloop)* & temp(iloop,jloop,1,tloop)/g0 ! za_gcm(iloop,jloop,1,tloop)=zlocal(iloop,jloop,1,tloop)+ & ! phis(iloop,jloop)/g0 do kloop=2,altlength ! compute sigma level of layer sigma(kloop)=press(iloop,jloop,kloop,tloop)/ps(iloop,jloop,tloop) ! compute "mean" temperature of the layer if(temp(iloop,jloop,kloop,tloop).eq. & temp(iloop,jloop,kloop-1,tloop)) then Tmean=temp(iloop,jloop,kloop,tloop) else Tmean=(temp(iloop,jloop,kloop,tloop)- & temp(iloop,jloop,kloop-1,tloop))/ & log(temp(iloop,jloop,kloop,tloop)/ & temp(iloop,jloop,kloop-1,tloop)) endif ! compute mean value of R of the layer Rmean=0.5*(R(iloop,jloop,kloop,tloop)+R(iloop,jloop,kloop-1,tloop)) ! compute gravitational acceleration (at altitude zaeroid(kloop-1)) ! NB: zareoid=zsurface+phis/g0 gz=g0*(a0**2)/ & (a0+zs_gcm(iloop,jloop,kloop-1,tloop)+phis(iloop,jloop)/g0)**2 ! compute zs_gcm(iloop,jloop,kloop,tloop) zs_gcm(iloop,jloop,kloop,tloop)=zs_gcm(iloop,jloop,kloop-1,tloop)- & log(sigma(kloop)/sigma(kloop-1))*Rmean*Tmean/gz ! compute za_gcm(kloop) ! za_gcm(iloop,jloop,kloop,tloop)=zlocal(iloop,jloop,kloop,tloop)+ & ! phis(iloop,jloop)/g0 enddo ! kloop enddo ! iloop enddo ! jloop enddo ! tloop ! Cleanup deallocate(sigma) deallocate(R) !deallocate(zlocal) end subroutine build_gcm_zs !=============================================================================== !#include"build_gcm_alt.F90" subroutine build_gcm_za(lonlength,latlength,altlength,timelength, & phis,ps,press,temp,rho,za_gcm) !============================================================================== ! Purpose: Integrate hydrostatic equation in order to build the "above areoid ! altitudes" corresponding to GCM atmospheric levels !============================================================================== implicit none !============================================================================== ! Arguments: !============================================================================== integer,intent(in) :: lonlength ! # of points along longitude integer,intent(in) :: latlength ! # of points along latitude integer,intent(in) :: altlength ! # of points along altitude integer,intent(in) :: timelength ! # of stored time steps real,dimension(lonlength,latlength),intent(in) :: phis ! phis(:,:) is the ground geopotential real,dimension(lonlength,latlength,timelength),intent(in) :: ps ! ps(:,:) is the surface pressure real,dimension(lonlength,latlength,altlength,timelength),intent(in) :: press ! press(:,:,:,:) is atmospheric pressure real,dimension(lonlength,latlength,altlength,timelength),intent(in) :: temp ! temp(:,:,:,:) is atmospheric temperature real,dimension(lonlength,latlength,altlength,timelength),intent(in) :: rho ! rho(:,:,:,:) is atmospheric density real,dimension(lonlength,latlength,altlength,timelength),intent(out) :: za_gcm ! za_gcm(:,:,:,:) are the above aroid heights of GCM levels !=============================================================================== ! Local variables: !=============================================================================== real,dimension(:),allocatable :: sigma ! GCM sigma levels real,dimension(:,:,:,:),allocatable :: R ! molecular gas constant real,dimension(:,:,:,:),allocatable :: zlocal ! local above surface altitude real :: Rmean ! mean value of R for a given level real :: Tmean ! "mean" value of temperature for a given level integer iloop,jloop,kloop,tloop ! Parameters needed to integrate hydrostatic equation: real,parameter :: g0=3.7257964 !g0: exact mean gravity at radius=3396.km (Lemoine et al. 2001) real,parameter :: a0=3396.E3 !a0: 'mean' Mars radius=3396.km real :: gz ! gz: gravitational acceleration at a given (zareoid) altitude !=============================================================================== ! 1. Various initialisations !=============================================================================== allocate(sigma(altlength)) allocate(R(lonlength,latlength,altlength,timelength)) allocate(zlocal(lonlength,latlength,altlength,timelength)) !============================================================================== ! 2. Compute Molecular Gas Constant (J kg-1 K-1) at every grid point ! --later used to integrate the hydrostatic equation-- !============================================================================== do tloop=1,timelength do kloop=1,altlength do jloop=1,latlength do iloop=1,lonlength R(iloop,jloop,kloop,tloop)=press(iloop,jloop,kloop,tloop)/ & (rho(iloop,jloop,kloop,tloop)* & temp(iloop,jloop,kloop,tloop)) enddo enddo enddo enddo !=============================================================================== ! 3. Integrate hydrostatic equation to compute zlocal and za_gcm !=============================================================================== do tloop=1,timelength do jloop=1,latlength do iloop=1,lonlength ! handle case of first altitude level sigma(1)=press(iloop,jloop,1,tloop)/ps(iloop,jloop,tloop) zlocal(iloop,jloop,1,tloop)=-log(sigma(1))*R(iloop,jloop,1,tloop)* & temp(iloop,jloop,1,tloop)/g0 za_gcm(iloop,jloop,1,tloop)=zlocal(iloop,jloop,1,tloop)+ & phis(iloop,jloop)/g0 do kloop=2,altlength ! compute sigma level of layer sigma(kloop)=press(iloop,jloop,kloop,tloop)/ps(iloop,jloop,tloop) ! compute "mean" temperature of the layer if(temp(iloop,jloop,kloop,tloop).eq. & temp(iloop,jloop,kloop-1,tloop)) then Tmean=temp(iloop,jloop,kloop,tloop) else Tmean=(temp(iloop,jloop,kloop,tloop)- & temp(iloop,jloop,kloop-1,tloop))/ & log(temp(iloop,jloop,kloop,tloop)/ & temp(iloop,jloop,kloop-1,tloop)) endif ! compute mean value of R of the layer Rmean=0.5*(R(iloop,jloop,kloop,tloop)+R(iloop,jloop,kloop-1,tloop)) ! compute gravitational acceleration (at altitude zaeroid(kloop-1)) gz=g0*(a0**2)/(a0+za_gcm(iloop,jloop,kloop-1,tloop))**2 ! compute zlocal(kloop) zlocal(iloop,jloop,kloop,tloop)=zlocal(iloop,jloop,kloop-1,tloop)- & log(sigma(kloop)/sigma(kloop-1))*Rmean*Tmean/gz ! compute za_gcm(kloop) za_gcm(iloop,jloop,kloop,tloop)=zlocal(iloop,jloop,kloop,tloop)+ & phis(iloop,jloop)/g0 enddo ! kloop enddo ! iloop enddo ! jloop enddo ! tloop ! Cleanup deallocate(sigma) deallocate(R) deallocate(zlocal) end subroutine build_gcm_za !=============================================================================== subroutine build_zs(altlength,have_sigma,sigma,aps,bps,zsurface) !============================================================================== ! Build generic above surface altitudes, using either sigma levels ! or hybrid vertical coordinates. ! In order to do so, we need to use scale heights that vary with altitude. ! The scale heights distribution is then set to vary from a min. to a max. ! following a tanh() law over an imposed transition range (see below). ! Similarily, a reference mean surface pressure (610 Pa) is used to ! compute generic above surface altitudes. !============================================================================== implicit none !============================================================================== ! Arguments: !============================================================================== integer,intent(in) :: altlength ! # of points along altitude logical,intent(in) :: have_sigma ! true if sigma(:) are known ! have_sigma is false if aps() and bps(:) are given instead real,dimension(altlength),intent(in) :: sigma ! sigma levels real,dimension(altlength),intent(in) :: aps ! hybrid pressure levels real,dimension(altlength),intent(in) :: bps ! hybrid sigma levels real,dimension(altlength),intent(out) :: zsurface ! altitudes (m) !=============================================================================== ! Local variables: !=============================================================================== real,dimension(:),allocatable :: H ! scale heights real,parameter :: H_low=9650 ! scale height at low altitudes real,parameter :: H_high=15000 ! scale height at high altitudes real,parameter :: trans_window=10 ! # of levels over which H(:) goes ! from H_low to H_high real,parameter :: lev_trans=32+trans_window/2 ! lev_trans: level at which H(lev_trans)=(H_low+H_high)/2 ! N.B.: keep lev_trans and lev_trans as reals to avoid truncation issues real,parameter :: P_ref=610 ! reference surface pressure used to build zsurface integer i ! 1. Build scale heights allocate(H(altlength)) do i=1,altlength H(i)=H_low+(H_high-H_low)*0.5*(1.0+tanh(6.*(i-lev_trans)/trans_window)) enddo ! 2. Compute zsurface(:) if (have_sigma) then ! use sigma levels do i=1,altlength zsurface(i)=-H(i)*log(sigma(i)*P_ref) enddo else ! use hybrid coordinates do i=1,altlength zsurface(i)=-H(i)*log((aps(i)/P_ref)+bps(i)) enddo endif ! Cleanup deallocate(H) end subroutine build_zs !=============================================================================== subroutine crude_gcm_zs(lonlength,latlength,altlength,timelength, & phis,ps,press,temp,zs_gcm) !============================================================================== ! Purpose: Integrate hydrostatic equation in order to build the "above local ! surface altitudes" corresponding to GCM atmospheric levels !============================================================================== implicit none !============================================================================== ! Arguments: !============================================================================== integer,intent(in) :: lonlength ! # of points along longitude integer,intent(in) :: latlength ! # of points along latitude integer,intent(in) :: altlength ! # of points along altitude integer,intent(in) :: timelength ! # of stored time steps real,dimension(lonlength,latlength),intent(in) :: phis ! phis(:,:) is the ground geopotential real,dimension(lonlength,latlength,timelength),intent(in) :: ps ! ps(:,:) is the surface pressure real,dimension(lonlength,latlength,altlength,timelength),intent(in) :: press ! press(:,:,:,:) is atmospheric pressure real,dimension(lonlength,latlength,altlength,timelength),intent(in) :: temp ! temp(:,:,:,:) is atmospheric temperature real,dimension(lonlength,latlength,altlength,timelength),intent(out) :: zs_gcm ! zs_gcm(:,:,:,:) are the above local surface altitudes of GCM levels !=============================================================================== ! Local variables: !=============================================================================== real,dimension(:),allocatable :: sigma ! GCM sigma levels real,parameter :: R=191 ! molecular gas constant !real,dimension(:,:,:,:),allocatable :: zlocal ! local above surface altitude real :: Tmean ! "mean" value of temperature for a given level integer iloop,jloop,kloop,tloop ! Parameters needed to integrate hydrostatic equation: real,parameter :: g0=3.7257964 !g0: exact mean gravity at radius=3396.km (Lemoine et al. 2001) real,parameter :: a0=3396.E3 !a0: 'mean' Mars radius=3396.km real :: gz ! gz: gravitational acceleration at a given (zareoid) altitude !=============================================================================== ! 1. Various initialisations !=============================================================================== allocate(sigma(altlength)) !allocate(zlocal(lonlength,latlength,altlength,timelength)) !=============================================================================== ! 2. Integrate hydrostatic equation to compute zlocal and za_gcm !=============================================================================== do tloop=1,timelength do jloop=1,latlength do iloop=1,lonlength ! handle case of first altitude level sigma(1)=press(iloop,jloop,1,tloop)/ps(iloop,jloop,tloop) zs_gcm(iloop,jloop,1,tloop)=-log(sigma(1))* & R*temp(iloop,jloop,1,tloop)/g0 ! za_gcm(iloop,jloop,1,tloop)=zlocal(iloop,jloop,1,tloop)+ & ! phis(iloop,jloop)/g0 do kloop=2,altlength ! compute sigma level of layer sigma(kloop)=press(iloop,jloop,kloop,tloop)/ps(iloop,jloop,tloop) ! compute "mean" temperature of the layer if(temp(iloop,jloop,kloop,tloop).eq. & temp(iloop,jloop,kloop-1,tloop)) then Tmean=temp(iloop,jloop,kloop,tloop) else Tmean=(temp(iloop,jloop,kloop,tloop)- & temp(iloop,jloop,kloop-1,tloop))/ & log(temp(iloop,jloop,kloop,tloop)/ & temp(iloop,jloop,kloop-1,tloop)) endif ! compute gravitational acceleration (at altitude zaeroid(kloop-1)) ! NB: zareoid=zsurface+phis/g0 gz=g0*(a0**2)/ & (a0+zs_gcm(iloop,jloop,kloop-1,tloop)+phis(iloop,jloop)/g0)**2 ! compute zs_gcm(iloop,jloop,kloop,tloop) zs_gcm(iloop,jloop,kloop,tloop)=zs_gcm(iloop,jloop,kloop-1,tloop)- & log(sigma(kloop)/sigma(kloop-1))*R*Tmean/gz ! compute za_gcm(kloop) ! za_gcm(iloop,jloop,kloop,tloop)=zlocal(iloop,jloop,kloop,tloop)+ & ! phis(iloop,jloop)/g0 enddo ! kloop enddo ! iloop enddo ! jloop enddo ! tloop ! Cleanup deallocate(sigma) !deallocate(zlocal) end subroutine crude_gcm_zs !=============================================================================== !#include"crude_gcm_alt.F90" subroutine crude_gcm_za(lonlength,latlength,altlength,timelength, & phis,ps,press,temp,za_gcm) !============================================================================== ! Purpose: Integrate hydrostatic equation in order to build the "above areoid ! altitudes" corresponding to GCM atmospheric levels !============================================================================== implicit none !============================================================================== ! Arguments: !============================================================================== integer,intent(in) :: lonlength ! # of points along longitude integer,intent(in) :: latlength ! # of points along latitude integer,intent(in) :: altlength ! # of points along altitude integer,intent(in) :: timelength ! # of stored time steps real,dimension(lonlength,latlength),intent(in) :: phis ! phis(:,:) is the ground geopotential real,dimension(lonlength,latlength,timelength),intent(in) :: ps ! ps(:,:) is the surface pressure real,dimension(lonlength,latlength,altlength,timelength),intent(in) :: press ! press(:,:,:,:) is atmospheric pressure real,dimension(lonlength,latlength,altlength,timelength),intent(in) :: temp ! temp(:,:,:,:) is atmospheric temperature real,dimension(lonlength,latlength,altlength,timelength),intent(out) :: za_gcm ! za_gcm(:,:,:,:) are the above aroid heights of GCM levels !=============================================================================== ! Local variables: !=============================================================================== real,dimension(:),allocatable :: sigma ! GCM sigma levels real,parameter :: R=191 ! molecular gas constant real,dimension(:,:,:,:),allocatable :: zlocal ! local above surface altitude real :: Tmean ! "mean" value of temperature for a given level integer iloop,jloop,kloop,tloop ! Parameters needed to integrate hydrostatic equation: real,parameter :: g0=3.7257964 !g0: exact mean gravity at radius=3396.km (Lemoine et al. 2001) real,parameter :: a0=3396.E3 !a0: 'mean' Mars radius=3396.km real :: gz ! gz: gravitational acceleration at a given (zareoid) altitude !=============================================================================== ! 1. Various initialisations !=============================================================================== allocate(sigma(altlength)) allocate(zlocal(lonlength,latlength,altlength,timelength)) !=============================================================================== ! 2. Integrate hydrostatic equation to compute zlocal and za_gcm !=============================================================================== do tloop=1,timelength do jloop=1,latlength do iloop=1,lonlength ! handle case of first altitude level sigma(1)=press(iloop,jloop,1,tloop)/ps(iloop,jloop,tloop) zlocal(iloop,jloop,1,tloop)=-log(sigma(1))* & R*temp(iloop,jloop,1,tloop)/g0 za_gcm(iloop,jloop,1,tloop)=zlocal(iloop,jloop,1,tloop)+ & phis(iloop,jloop)/g0 do kloop=2,altlength ! compute sigma level of layer sigma(kloop)=press(iloop,jloop,kloop,tloop)/ps(iloop,jloop,tloop) ! compute "mean" temperature of the layer if(temp(iloop,jloop,kloop,tloop).eq. & temp(iloop,jloop,kloop-1,tloop)) then Tmean=temp(iloop,jloop,kloop,tloop) else Tmean=(temp(iloop,jloop,kloop,tloop)- & temp(iloop,jloop,kloop-1,tloop))/ & log(temp(iloop,jloop,kloop,tloop)/ & temp(iloop,jloop,kloop-1,tloop)) endif ! compute gravitational acceleration (at altitude zaeroid(kloop-1)) gz=g0*(a0**2)/(a0+za_gcm(iloop,jloop,kloop-1,tloop))**2 ! compute zlocal(kloop) zlocal(iloop,jloop,kloop,tloop)=zlocal(iloop,jloop,kloop-1,tloop)- & log(sigma(kloop)/sigma(kloop-1))*R*Tmean/gz ! compute za_gcm(kloop) za_gcm(iloop,jloop,kloop,tloop)=zlocal(iloop,jloop,kloop,tloop)+ & phis(iloop,jloop)/g0 enddo ! kloop enddo ! iloop enddo ! jloop enddo ! tloop ! Cleanup deallocate(sigma) deallocate(zlocal) end subroutine crude_gcm_za !=============================================================================== subroutine p_coord_interp(lonlen,latlen,altlen,tlen,newaltlen, & missing,ps,press,gcmdata,plevels,newdata) !============================================================================== ! Purpose: ! Recast a 4D (spatio-temporal) GCM dataset, which has a vertical coordinate ! in pseudo-altitude, into a dataset which has a vertical coordinate at given ! pressure levels !============================================================================== implicit none !============================================================================== ! Arguments: !============================================================================== integer,intent(in) :: lonlen ! # of points along longitude (GCM dataset) integer,intent(in) :: latlen ! # of points along latitude (GCM dataset) integer,intent(in) :: altlen ! # of points along altitude (GCM dataset) integer,intent(in) :: tlen ! # of stored time steps (GCM dataset) integer,intent(in) :: newaltlen ! # of points along altitude real,intent(in) :: missing ! missing value real,dimension(lonlen,latlen,tlen),intent(in) :: ps ! GCM surface pressure real,dimension(lonlen,latlen,altlen,tlen),intent(in) :: press ! GCM pressure real,dimension(lonlen,latlen,altlen,tlen),intent(in) :: gcmdata ! GCM dataset real,dimension(newaltlen),intent(in) :: plevels ! plevels(:) pressure levels at which gcmdata is to be interpolated real,dimension(lonlen,latlen,newaltlen,tlen),intent(out) :: newdata ! newdata(:,:,:,:) gcmdata recasted along vertical coordinate !============================================================================== ! Local variables: !============================================================================== real,dimension(:),allocatable :: lnp ! lnp(:): used to store log(pressure) values real,dimension(:),allocatable :: q ! q(:): used to store values along vertical direction real :: x,y ! x,y: temporary variables integer :: iloop,jloop,kloop,tloop allocate(lnp(altlen)) allocate(q(altlen)) do tloop=1,tlen do jloop=1,latlen do iloop=1,lonlen ! build lnp() and q() along vertical coordinate do kloop=1,altlen lnp(kloop)=-log(press(iloop,jloop,kloop,tloop)) q(kloop)=gcmdata(iloop,jloop,kloop,tloop) enddo ! Interpolation do kloop=1,newaltlen ! compute, by interpolation, value at pressure level plevels(kloop) if ((plevels(kloop).le.ps(iloop,jloop,tloop)).and. & (plevels(kloop).ge.press(iloop,jloop,altlen,tloop))) then x=-log(plevels(kloop)) call interpolf(x,y,missing,lnp,q,altlen) newdata(iloop,jloop,kloop,tloop) = y else ! if plevels(kloop) is out of range, ! assign a "missing_value" at this node newdata(iloop,jloop,kloop,tloop) = missing endif enddo enddo !iloop enddo !jloop enddo !tloop ! Cleanup deallocate(lnp) deallocate(q) end subroutine p_coord_interp !=============================================================================== !#include"za_coord_interp.F90" subroutine z_coord_interp(lonlen,latlen,altlen,tlen,newaltlen, & missing,z_gcm,gcmdata,flag,z_new,newdata) !============================================================================== ! Purpose: ! Recast a 4D (spatio-temporal) GCM dataset 'gcmdata', for which corresponding ! grid values of altitude are known 'z_gcm', onto a new altitude grid 'z_new'. ! "Altitudes" can be above areoid or above local surface altitudes, as long as ! both 'z_gcm' and 'znew' refer to altitudes of a same type. ! Note: If altitudes in 'znew' fall outside of the range of altitudes ! in 'z_gcm' then corresponding 'newdata' value is set to 'missing'. ! 'z_gcm' and 'znew' altitudes must be monotone increasing sequences. !============================================================================== implicit none !============================================================================== ! Arguments: !============================================================================== integer,intent(in) :: lonlen ! # of points along longitude (GCM dataset) integer,intent(in) :: latlen ! # of points along latitude (GCM dataset) integer,intent(in) :: altlen ! # of points along altitude (GCM dataset) integer,intent(in) :: tlen ! # of stored time steps (GCM dataset) integer,intent(in) :: newaltlen ! # of points along altitude real,intent(in) :: missing ! missing value real,dimension(lonlen,latlen,altlen,tlen),intent(in) :: z_gcm !z_gcm(:,:,:,) GCM grid points altitude (above areoid or above surface) real,dimension(lonlen,latlen,altlen,tlen),intent(in) :: gcmdata ! GCM dataset integer,intent(in) :: flag ! flag (==1 for 'log' interpolation) ! flag==0 (standard linear interpolation) real,dimension(newaltlen),intent(in) :: z_new ! z_new(:) altitudes (above areoid or surface) at which data must be recast real,dimension(lonlen,latlen,newaltlen,tlen),intent(out) :: newdata ! newdata(:,:,:,:) gcmdata recasted along vertical coordinate !============================================================================== ! Local variables: !============================================================================== real,dimension(:),allocatable :: za,q real,dimension(:),allocatable :: logq ! za(:): used to store z_gcm at a given location ! q(:): used to store field values along the vertical direction real :: x,y ! temporary variables integer :: iloop,jloop,kloop,tloop allocate(za(altlen)) allocate(q(altlen)) allocate(logq(altlen)) if (flag.eq.0) then do tloop=1,tlen do jloop=1,latlen do iloop=1,lonlen do kloop=1,altlen ! extract the vertical coordinates za(kloop)=z_gcm(iloop,jloop,kloop,tloop) ! store values along altitude q(kloop)=gcmdata(iloop,jloop,kloop,tloop) enddo !kloop ! Interpolation do kloop=1,newaltlen ! Check if z_new(kloop) is within range if ((z_new(kloop).ge.z_gcm(iloop,jloop,1,tloop)).and. & (z_new(kloop).le.z_gcm(iloop,jloop,altlen,tloop))) then ! z_new(kloop) is within range x=z_new(kloop) call interpolf(x,y,missing,za,q,altlen) newdata(iloop,jloop,kloop,tloop)=y else ! z_new(kloop) is out of range newdata(iloop,jloop,kloop,tloop)=missing endif enddo !kloop enddo !iloop enddo !jloop enddo !tloop else ! case when flag=1 (i.e.: rho) do tloop=1,tlen do jloop=1,latlen do iloop=1,lonlen do kloop=1,altlen ! extract the vertical coordinates za(kloop)=z_gcm(iloop,jloop,kloop,tloop) ! store log values along altitude logq(kloop)=log(gcmdata(iloop,jloop,kloop,tloop)) enddo !kloop ! Interpolation do kloop=1,newaltlen ! Check if z_new(kloop) is within range if ((z_new(kloop).ge.z_gcm(iloop,jloop,1,tloop)).and. & (z_new(kloop).le.z_gcm(iloop,jloop,altlen,tloop))) then ! z_new(kloop) is within range x=z_new(kloop) call interpolf(x,y,missing,za,logq,altlen) newdata(iloop,jloop,kloop,tloop)=exp(y) else ! z_new(kloop) is out of range newdata(iloop,jloop,kloop,tloop)=missing endif enddo !kloop enddo !iloop enddo !jloop enddo !tloop endif ! Cleanup deallocate(za) deallocate(q) deallocate(logq) end subroutine z_coord_interp !=============================================================================== subroutine zs_coord_interp(lonlen,latlen,altlen,tlen,newaltlen, & missing,z_gcm,gcmdata,phis,press,temp,rho, & have_rho,flag,z_new,newdata) !============================================================================== ! Purpose: ! Recast a 4D (spatio-temporal) GCM dataset 'gcmdata', for which corresponding ! grid values of altitude are known 'z_gcm', onto a new altitude grid 'z_new'. ! "Altitudes" must be above local surface altitudes. ! Notes: ! 'z_gcm' and 'znew' altitudes must be monotone increasing sequences. ! If altitudes in 'znew(i)' fall below those in 'z_gcm(:,:,1,:)', then ! 'newdata(:,:,i,:)' is set to constant value 'gcmdata(:,:,1,:)' if ! flag=0, and extrapolated (exponentially) if flag=1. ! If altitudes in 'znew(i)' are above those in 'z_gcm(:,:,altlen,:)', then ! 'newdata(:,:,i,:)' is set to constant value 'gcmdata(:,:,altlen,:)' ! if flag=0, and extrapolated (exponentially) if flag=1. !============================================================================== implicit none !============================================================================== ! Arguments: !============================================================================== integer,intent(in) :: lonlen ! # of points along longitude (GCM dataset) integer,intent(in) :: latlen ! # of points along latitude (GCM dataset) integer,intent(in) :: altlen ! # of points along altitude (GCM dataset) integer,intent(in) :: tlen ! # of stored time steps (GCM dataset) integer,intent(in) :: newaltlen ! # of points along altitude real,intent(in) :: missing ! missing value real,dimension(lonlen,latlen,altlen,tlen),intent(in) :: z_gcm !z_gcm(:,:,:,) GCM grid points altitude (above areoid or above surface) real,dimension(lonlen,latlen,altlen,tlen),intent(in) :: gcmdata ! GCM dataset real,dimension(lonlen,latlen),intent(in) :: phis ! phis(:,:) is the ground geopotential real,dimension(lonlen,latlen,altlen,tlen),intent(in) :: press ! press(:,:,:,:) is atmospheric pressure on GCM levels real,dimension(lonlen,latlen,altlen,tlen),intent(in) :: temp ! temp(:,:,:,:) is atmospheric temperature on GCM levels real,dimension(lonlen,latlen,altlen,tlen),intent(in) :: rho ! rho(:,:,:,:) is atmospheric density on GCM levels logical,intent(in) :: have_rho ! trueif we have density at hand integer,intent(in) :: flag ! flag (==1 for 'log' interpolation) ! flag==0 (standard linear interpolation) real,dimension(newaltlen),intent(in) :: z_new ! z_new(:) altitudes (above areoid or surface) at which data must be recast real,dimension(lonlen,latlen,newaltlen,tlen),intent(out) :: newdata ! newdata(:,:,:,:) gcmdata recasted along vertical coordinate !============================================================================== ! Local variables: !============================================================================== real,dimension(:),allocatable :: z,q real,dimension(:),allocatable :: logq ! z(:): used to store z_gcm at a given location ! q(:): used to store field values along the vertical direction real,dimension(:,:,:),allocatable :: Rbottom,Rtop ! values of R (gas constant) below and above GCM layers real :: x,y ! temporary variables integer :: iloop,jloop,kloop,tloop real,parameter :: g0=3.7257964 !g0: exact mean gravity at radius=3396.km (Lemoine et al. 2001) real,parameter :: a0=3396.E3 !a0: 'mean' Mars radius=3396.km real,parameter :: Rmean=191 ! molecular gas constant real :: gz ! gz: gravitational acceleration at a given (above areoid) altitude allocate(z(altlen)) allocate(q(altlen)) allocate(logq(altlen)) ! 1. Build Rbottom and Rtop (only necessary if flag=1) if (flag.eq.1) then allocate(Rbottom(lonlen,latlen,tlen)) allocate(Rtop(lonlen,latlen,tlen)) if (have_rho) then do iloop=1,lonlen do jloop=1,latlen do tloop=1,tlen Rbottom(iloop,jloop,tloop)=press(iloop,jloop,1,tloop)/ & (rho(iloop,jloop,1,tloop)* & temp(iloop,jloop,1,tloop)) Rtop(iloop,jloop,tloop)=press(iloop,jloop,altlen,tloop)/ & (rho(iloop,jloop,altlen,tloop)* & temp(iloop,jloop,altlen,tloop)) enddo enddo enddo else ! we don't have density at hand; use mean molecular gas constant value Rbottom(:,:,:)=Rmean Rtop(:,:,:)=Rmean endif endif ! if (flag.eq.1) ! 2. Interpolation if (flag.eq.0) then do tloop=1,tlen do jloop=1,latlen do iloop=1,lonlen ! preliminary stuff do kloop=1,altlen ! extract the vertical coordinates z(kloop)=z_gcm(iloop,jloop,kloop,tloop) ! store values along altitude q(kloop)=gcmdata(iloop,jloop,kloop,tloop) enddo !kloop ! Interpolation do kloop=1,newaltlen if (z_new(kloop).lt.z_gcm(iloop,jloop,1,tloop)) then ! z_new(kloop) is below lowest GCM level newdata(iloop,jloop,kloop,tloop)=gcmdata(iloop,jloop,1,tloop) else if (z_new(kloop).gt.z_gcm(iloop,jloop,altlen,tloop)) then ! z_new(kloop) is above highest GCM level newdata(iloop,jloop,kloop,tloop)=gcmdata(iloop,jloop,altlen,tloop) else ! z_new(kloop) is within range x=z_new(kloop) call interpolf(x,y,missing,z,q,altlen) newdata(iloop,jloop,kloop,tloop)=y endif enddo !kloop enddo !iloop enddo !jloop enddo !tloop else ! case when flag=1 (i.e.: rho or pressure) do tloop=1,tlen do jloop=1,latlen do iloop=1,lonlen ! preliminary stuff do kloop=1,altlen ! extract the vertical coordinates z(kloop)=z_gcm(iloop,jloop,kloop,tloop) ! store log values along altitude logq(kloop)=log(gcmdata(iloop,jloop,kloop,tloop)) enddo !kloop ! Interpolation do kloop=1,newaltlen if (z_new(kloop).lt.z_gcm(iloop,jloop,1,tloop)) then ! z_new(kloop) is below lowest GCM level newdata(iloop,jloop,kloop,tloop)=gcmdata(iloop,jloop,1,tloop)* & exp(((z_gcm(iloop,jloop,1,tloop)-z_new(kloop))*g0)/ & (Rbottom(iloop,jloop,tloop)* & temp(iloop,jloop,1,tloop))) else if (z_new(kloop).gt.z_gcm(iloop,jloop,altlen,tloop)) then ! z_new(kloop) is above highest GCM level ! NB: zareoid=zsurface+phis/g0 gz=g0*a0*a0/(a0+z_new(kloop)+phis(iloop,jloop)/g0)**2 newdata(iloop,jloop,kloop,tloop)=gcmdata(iloop,jloop,altlen,tloop)* & exp(((z_gcm(iloop,jloop,altlen,tloop)-z_new(kloop))*gz)/ & (Rtop(iloop,jloop,tloop)* & temp(iloop,jloop,altlen,tloop))) else ! z_new(kloop) is within range x=z_new(kloop) call interpolf(x,y,missing,z,logq,altlen) newdata(iloop,jloop,kloop,tloop)=exp(y) endif enddo !kloop enddo !iloop enddo !jloop enddo !tloop endif ! Cleanup deallocate(z) deallocate(q) deallocate(logq) if (flag.eq.1) then deallocate(Rbottom) deallocate(Rtop) endif end subroutine zs_coord_interp !=============================================================================== subroutine interpolf(x,y,missing,xd,yd,nd) !============================================================================== ! Purpose: ! Yield y=f(x), where xd() end yd() are arrays of known values, ! using linear interpolation ! If x is not included in the interval spaned by xd(), then y is set ! to a default value 'missing' ! Note: ! Array xd() should contain ordered (either increasing or decreasing) abscissas !============================================================================== implicit none !============================================================================== ! Arguments: !============================================================================== real,intent(in) :: x ! abscissa at which interpolation is to be done real,intent(in) :: missing ! missing value (if no interpolation is performed) integer :: nd ! size of arrays real,dimension(nd),intent(in) :: xd ! array of known absissas real,dimension(nd),intent(in) :: yd ! array of correponding values real,intent(out) :: y ! interpolated value !============================================================================== ! Local variables: !============================================================================== integer :: i ! default: set y to 'missing' y=missing do i=1,nd-1 if (((x.ge.xd(i)).and.(x.le.xd(i+1))).or.& ((x.le.xd(i)).and.(x.ge.xd(i+1)))) then y=yd(i)+(x-xd(i))*(yd(i+1)-yd(i))/(xd(i+1)-xd(i)) exit endif enddo end subroutine interpolf