3.3. A Fortran 90 example

This is an example of how to use this the Fortran 90 component of Spanlib. In this example, an analysis of Sea surface Temperature anomalies using data stored in the netcdf format. The dominant oscillatory mode of the El Nino variability is then extracted and stored in a netcdf file.

Example 1. F90 example

You can run this example typing (from package directory):

<user> cd example && make

This command will first try to download the dataset, then compile the following f90 program, run it, and finally try to visualise output data. See Section 2.1, “Requirements” if it fails.

The sources of this example is highlighted below, with comments inside.

! File: example.f90
! This file is part of the SpanLib library.
! Copyright (C) 2006  Stephane Raynaud
! Contact: stephane dot raynaud at gmail dot com
! This library is free software; you can redistribute it and/or
! modify it under the terms of the GNU Lesser General Public
! License as published by the Free Software Foundation; either
! version 2.1 of the License, or (at your option) any later version.
! This library is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! Lesser General Public License for more details.
! You should have received a copy of the GNU Lesser General Public
! License along with this library; if not, write to the Free Software
! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA

program example

	! This simple example shows how to use all subroutines from this package.
	! Warning: it requires netcdf for in/outputs.
	! We start from longitude/latitude/time value of Pacific Sea Surface Temperature
	! that include the El Nino Southern Oscillation signal.
	! Input is the netcdf file data.cdf.
	! We remove land points from the initial array according
	! to the netcdf missing_value attribute of the analysed variable (data are "packed").
	! A PCA is used to reduce the degrees of freedoom before MSSA analysis.
	! Weights for PCA are computed as a fonction of latitude.
	! Then, we assume that we have already identified an oscillation (after tests).
	! This oscillation, given by a pair of MSSA modes, is then
	! reconstructed from the MSSA and PCA spaces.
	! Finally, phase composites are computed from this reconstructed oscillation.
	! The oscillation is outputed in a netcdf file (pair_1.cdf).
	! The initial data set (data.cdf):
	! - origin: updated Reynolds and Smith (1996) SST (netcdf file)
	! - origin url: data selector from http://iridl.ldeo.columbia.edu
	! - how to get it [10Mb]: http://stefdeperou.free.fr/pub/data.cdf
	! - area of study: tropical pacific [130.5E:75.5W, 29.5S:29.5N] (155x60 grid points)
	! - period of study: Jan1982:Dec2005 (288 time steps)
	! Parameters:
	! - Only the first 20 PCs are retainedand given to the MSSA
	! - A window of 7 years (84 months) is chosen for the MSSA
	! - Phase composites use 8 phases
	! - An offset of 0.4 and is used for composites
	! - The first phase of composites is set at 180 degrees (minimal value)

	use spanlib
	use netcdf

	implicit none

	! Parameters
	! ----------
	integer,parameter :: nkeep_pca=5, nwindow=84, first_mode=1, &
		& nphases=8
	real, parameter :: offset=0., first_phase=180., &
		& new_missing_value=-999.
	character(len=20), parameter :: input_nc_file="data2.cdf", &
		& output_nc_file="output.nc", var_name='ssta'

	! Other declarations
	! ------------------
	real, allocatable :: field(:,:,:), weights(:,:), &
		& lat(:), lon(:), time(:)
	real, allocatable :: reco(:,:,:), phasecomps(:,:,:)
	logical, allocatable :: mask(:,:)
	real, allocatable :: packed_field(:,:), packed_weights(:), &
		& packed_phasecomps(:,:), stphasecomps(:,:)
	real, allocatable :: eof(:,:), ev(:), pc(:,:), &
		& stpair(:,:), pair(:,:)
	real, allocatable :: steof(:,:),stpc(:,:),stev(:)
	character(len=20) :: dim_names(3), dim_name, &
		& lon_units, lat_units, var_units, &
		&	lon_name, lat_name, time_name, time_units
	integer :: ncid, dimid, dimids(4), varid, dims(3), thisdim, &
		& lonid, latid, phaseid, timeid, phcoid, recoid, origid
	integer(kind=4) :: i, nspace, nlon, nlat, ntime
	real :: pi, missing_value

	! Get the initial sst field from the netcdf file
	! ----------------------------------------------
	print*,'Reading inputs...'
	call err(nf90_open(input_nc_file, nf90_nowrite, ncid))
	call err(nf90_inq_varid(ncid, var_name, varid))
	call err(nf90_inquire_variable(ncid, varid, dimids=dimids(1:3)))
	call err(nf90_inquire_dimension(ncid, dimids(1), &
		&	name=lon_name, len=nlon))
	call err(nf90_inquire_dimension(ncid, dimids(2), &
		&	name=lat_name, len=nlat))
	call err(nf90_inquire_dimension(ncid, dimids(3), &
		&	name=time_name, len=ntime))
	call err(nf90_get_var(ncid, varid, field))
	call err(nf90_get_att(ncid, varid, 'missing_value', &
		& missing_value))
	call err(nf90_get_att(ncid, varid, 'units', var_units))
	call err(nf90_inq_varid(ncid, lon_name, varid))
	call err(nf90_get_var(ncid, varid, lon))
	call err(nf90_get_att(ncid, varid, 'units', lon_units))
	call err(nf90_inq_varid(ncid, lat_name, varid))
	call err(nf90_get_var(ncid, varid, lat))
	call err(nf90_get_att(ncid, varid, 'units', lat_units))
	call err(nf90_inq_varid(ncid, time_name, varid))
	call err(nf90_get_var(ncid, varid, time))
	call err(nf90_get_att(ncid, varid, 'units', time_units))
	call err(nf90_close(ncid))

	! Format (pack) data to have only one space dimension
	! ---------------------------------------------------

	! Compute weights proportional to grid point area
	pi = cos(-1.)
	do i=1,nlat
		weights(:,i) = cos(lat(i)*pi/180.)
	end do

	! Now pack
	mask = (field(:,:,1) /= missing_value)
	allocate(packed_field(count(mask), ntime))
	do i=1, ntime
		packed_field(:,i) = pack(field(:,:,i), mask)
	end do
	nspace = count(mask)
	packed_weights = pack(weights, mask)

	! Perform a PCA to reduce the d.o.f
	! ---------------------------------
	allocate(eof(nspace, nkeep_pca))
	allocate(pc(ntime,   nkeep_pca))
	call sl_pca(packed_field, nkeep_pca, xeof=eof, &
		&	pc=pc, weights=packed_weights)

	! We send results from PCA to MSSA
	! --------------------------------
	allocate(steof(nkeep_pca*nwindow, first_mode+1))
	allocate(stpc(ntime-nwindow+1,    first_mode+1))
	allocate(stev(                    first_mode+1))
	call sl_mssa(transpose(pc), nwindow, first_mode+1, &
		&	steof=steof, stpc=stpc, ev=stev)

	! We reconstruct modes [first_mode + first_mode+1] of MSSA
	! --------------------------------------------------------

	allocate(stpair(nkeep_pca, ntime))
	call sl_mssarec(steof, stpc, nwindow, stpair, &
		&	istart=first_mode, iend=first_mode+1)
	deallocate(steof, stpc)

	! We compute phases composites for the reconstructed oscillation
	! ---------------------------------------------------------------
	allocate(stphasecomps(nkeep_pca, nphases))
	call sl_phasecomp(stpair, nphases, stphasecomps, &
		&	weights=packed_weights, &
		&	offset=offset, firstphase=first_phase)

	! We go back to the physical space for
	! the full oscillation AND its composites
	! ---------------------------------------
	allocate(pair(nspace, ntime))
	call sl_pcarec(eof, transpose(stpair), pair)
	allocate(packed_phasecomps(nspace, nphases))
	call sl_pcarec(eof, transpose(stphasecomps), packed_phasecomps)
	deallocate(stpair, eof, stphasecomps)

	! Unpacking
	! ---------
	do i=1, ntime
		reco(:,:,i) = unpack(pair(:,i), mask, new_missing_value)
		where(mask .eq. .false.)
			field(:,:,i) = new_missing_value
		end where
	end do
	do i=1, nphases
		phasecomps(:,:,i) = unpack(packed_phasecomps(:,i), mask, &
		 & new_missing_value)
	end do

	! Write out the phase composites of the first oscillation
	! -------------------------------------------------------
	print*,'Writing out...'
	! File
	call err(nf90_create(output_nc_file, nf90_write, ncid))
	! Dimensions
	call err(nf90_def_dim(ncid, 'lon', nlon, dimids(1)))
	call err(nf90_def_dim(ncid, 'lat', nlat, dimids(2)))
	call err(nf90_def_dim(ncid, 'time', ntime, dimids(3)))
	call err(nf90_def_dim(ncid, 'phase', nphases, dimids(4)))
	! Variables
	call err(nf90_def_var(ncid, 'lon', nf90_float, dimids(1), lonid))
	call err(nf90_put_att(ncid, lonid, 'long_name', 'Longitude'))
	call err(nf90_put_att(ncid, lonid, 'units', lon_units))
	call err(nf90_def_var(ncid, 'lat', nf90_float, dimids(2), latid))
	call err(nf90_put_att(ncid, latid, 'long_name', 'Latitude'))
	call err(nf90_put_att(ncid, latid, 'units', lat_units))
	call err(nf90_def_var(ncid, 'time', nf90_float, dimids(3), &
		& timeid))
	call err(nf90_put_att(ncid, timeid, 'long_name', 'Time'))
	call err(nf90_put_att(ncid, timeid, 'units', time_units))
	call err(nf90_def_var(ncid, 'phase', nf90_float, dimids(4), &
		& phaseid))
	call err(nf90_put_att(ncid, phaseid, 'long_name', 'Phase'))
	call err(nf90_put_att(ncid, phaseid, 'units', 'level'))
	call err(nf90_def_var(ncid, 'orig', nf90_float, dimids(1:3), &
	 & origid))
	call err(nf90_put_att(ncid, origid, 'long_name', &
		& 'SST anomaly / original field'))
	call err(nf90_put_att(ncid, origid, 'units', var_units))
	call err(nf90_put_att(ncid, origid, 'missing_value', &
		& new_missing_value))
	call err(nf90_def_var(ncid, 'reco1', nf90_float, dimids(1:3), &
		& recoid))
	call err(nf90_put_att(ncid, recoid, 'long_name', &
		& 'SST anomaly / reconstruction of first pair'))
	call err(nf90_put_att(ncid, recoid, 'units', var_units))
	call err(nf90_put_att(ncid, recoid, 'missing_value', &
		& new_missing_value))
	call err(nf90_def_var(ncid, 'pair1', nf90_float, &
		& (/dimids(1),dimids(2),dimids(4)/), phcoid))
	call err(nf90_put_att(ncid, phcoid, 'long_name', &
		&'SST anomaly / phase composite of first pair'))
	call err(nf90_put_att(ncid, phcoid, 'units', var_units))
	call err(nf90_put_att(ncid, phcoid, 'missing_value', &
		& new_missing_value))
	! Values
	call err(nf90_enddef(ncid))
	call err(nf90_put_var(ncid, lonid, lon))
	call err(nf90_put_var(ncid, latid, lat))
	call err(nf90_put_var(ncid, timeid, time))
	call err(nf90_put_var(ncid, phaseid, float((/(i,i=1,nphases)/))))
	call err(nf90_put_var(ncid, origid, field))
	call err(nf90_put_var(ncid, recoid, reco))
	call err(nf90_put_var(ncid, phcoid, phasecomps))
	call err(nf90_close(ncid))

end program example

subroutine err(jstatus)

	use netcdf

	integer :: jstatus

	if (jstatus .ne. nf90_noerr) then
		print *, trim(nf90_strerror(jstatus))
	end if

end subroutine err

SourceForge.net Logo
This document was generated using xsltproc and perl.