Index: trunk/LMDZ.GENERIC/README
===================================================================
--- trunk/LMDZ.GENERIC/README	(revision 1513)
+++ trunk/LMDZ.GENERIC/README	(revision 1514)
@@ -1116,2 +1116,3 @@
 == 10/02/2016 == MT
 - Added a generic routine (in utilities directory) to extrapolate ice fields.
+- Added a generic version of streamfunction (in utilities directory).
Index: trunk/LMDZ.GENERIC/utilities/streamfunction.F90
===================================================================
--- trunk/LMDZ.GENERIC/utilities/streamfunction.F90	(revision 1514)
+++ trunk/LMDZ.GENERIC/utilities/streamfunction.F90	(revision 1514)
@@ -0,0 +1,1152 @@
+module planet_const
+! planetary constants (set via init_planet_const routine)
+implicit none
+real :: a0 ! Mean planetary radius
+real :: g0 ! gravity
+real :: Rmean ! reduced gaz constant
+real :: omega ! rotation rate
+end module planet_const
+
+program streamfunction
+
+! ----------------------------------------------------------------------------
+! Program to calculate the stream function and the total angular momentum
+! from stats and diagfi files
+! input : diagfi.nc  / concat.nc / stats.nc kind of files
+! author: F. Gonzalez-Galindo
+!
+! Updates :
+! 
+! - 02/2016 (M. Turbet) : The code was adapted to work for any planet.
+!
+! ----------------------------------------------------------------------------
+
+use planet_const
+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 'aire' and 'cv')
+character (len=64) :: text ! to store some text
+
+character (len=50), dimension(:), allocatable :: var
+! var(): name(s) of variable(s) that will be concatenated
+character (len=100) :: filename
+! filename(): output file name
+integer :: nid,ierr,miss
+! nid: [netcdf] file ID #
+! ierr: [netcdf] subroutine returned error code
+! miss: [netcdf] subroutine returned error code
+integer :: tmpvarid
+! varid: temporary store a variable ID #
+integer tmpdimid ! temporarily store a dimension ID
+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
+real, dimension(:), allocatable:: lat,lon,alt,time
+! lat(): array, stores latitude coordinates
+! lon(): array, stores longitude coordinates
+! alt(): array, stores altitude coordinates
+! time(): array, stores time coordinates
+integer :: ilat,ilon,ialt,it,i
+! ilat: index for loops on latitude
+! ilon: index for loops on longitude 
+! ialt: index for loops on altitude
+! it: # index for loops on time
+! i: index for loops
+integer :: nout,latdimout,londimout,altdimout,timedimout,lonvarout,timevarout
+integer :: psivarout,momvarout,uvarout,vvarout,rhovarout,tempvarout
+integer :: psvarout, phisinitvarout
+! nout: [netcdf] output file ID
+! latdimout: [netcdf] output latitude (dimension) ID
+! londimout: [netcdf] output longitude (dimension) ID
+! altdimout: [netcdf] output altitude (dimension) ID
+! timedimout: [netcdf] output time (dimension) ID
+! lonvarout: [netcdf] ID of output "Time" variable
+integer lonlength ! # of grid points along longitude
+integer latlength ! # of grid points along latitude
+integer altlength ! # of grid point along altitude (of input datasets)
+integer timelength ! # of points along time
+real,dimension(:),allocatable :: aps,bps ! hybrid vertical coordinates
+real,dimension(:),allocatable :: sigma ! sigma levels
+real,dimension(:),allocatable :: lon_fake ! Fake longitude to be written
+real,dimension(:,:),allocatable :: aire ! Area of the box
+real,dimension(:,:),allocatable :: cv ! Conversion between natural and covariant v
+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 :: u ! GCM zonal wind
+real,dimension(:,:,:,:),allocatable :: v ! GCM meridional wind
+real,dimension(:,:,:,:),allocatable :: rho ! GCM atmospheric density
+real,dimension(:,:,:),allocatable :: vcont ! Covariant meridional wind
+real,dimension(:,:),allocatable :: vcontcum ! Zonal mean covariant meridional wind
+real,dimension(:,:,:),allocatable :: plev ! Pressure from hybrid coordinates
+real,dimension(:,:,:),allocatable :: mass ! Mass in a grid box
+real,dimension(:,:,:),allocatable :: dmass ! Mass variation
+real,dimension(:,:),allocatable :: psi ! Stream function
+real,dimension(:,:,:),allocatable :: mom ! Angular momentum
+real,dimension(:,:),allocatable :: momave ! Zonally averaged angular momentum
+real,dimension(:,:,:),allocatable :: ucum ! Temporally averaged zonal wind
+real,dimension(:,:,:),allocatable :: vcum ! Temporally averaged meridional wind
+real,dimension(:,:,:),allocatable :: rhocum ! Temporally averaged density
+real,dimension(:,:,:),allocatable :: tempcum ! Temporally averaged zonal wind
+real,dimension(:,:),allocatable :: pscum ! Temporally averaged zonal wind
+real,dimension(:,:),allocatable :: uzm ! Zonally averaged zonal wind
+real,dimension(:,:),allocatable :: vzm ! Zonally averaged meridional wind
+real,dimension(:,:),allocatable :: rhozm ! Zonally averaged density
+real,dimension(:,:),allocatable :: tempzm ! Zonally averaged temperature
+real,dimension(:),allocatable :: pszm ! Zonally averaged surface pressure
+real,dimension(:),allocatable :: phisinitzm ! Zonally averaged ground geopotential
+
+!===============================================================================
+! 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 input file'
+   stop ""
+endif
+
+! load planet constants (radius, gravity, ...)
+
+call init_planet_const(infid)
+
+!===============================================================================
+! 1.2 Get grids in lon,lat,alt,time,
+!     as well as hybrid coordinates aps() and bps() (or sigma levels sigma())
+!     aire() and cv() from input file
+!===============================================================================
+
+! 1.2.1 longitude, latitude, altitude and time
+
+! latitude
+ierr=NF_INQ_DIMID(infid,"latitude",tmpdimid)
+if (ierr.ne.NF_NOERR) then
+  stop "Error: Failed to get latitude dimension ID"
+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
+  stop "Error: Failed to get longitude dimension ID"
+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
+  stop "Error: Failed to get Time 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
+
+! altlength
+ierr=NF_INQ_DIMID(infid,"altitude",tmpdimid)
+if (ierr.ne.NF_NOERR) then
+  stop "Error: Failed to get altitude dimension ID"
+else
+  ierr=NF_INQ_VARID(infid,"altitude",tmpvarid)
+  if (ierr.ne.NF_NOERR) then
+     stop "Error: Failed to get altitude length"
+  else
+     ierr=NF_INQ_DIMLEN(infid,tmpdimid,altlength)
+     if (ierr.ne.NF_NOERR) then
+        stop "Error: Failed to get altitude length"
+     else   
+        allocate(alt(altlength))
+        ierr=NF_GET_VAR_REAL(infid,tmpvarid,alt)
+        if (ierr.ne.NF_NOERR) then
+           stop "Error: Failed reading Altitude"
+        endif
+     endif
+  endif
+endif
+
+! 1.2.2 Get hybrid coordinates 
+
+! look for hybrid coordinates
+
+! 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
+
+!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 ps"
+   endif
+endif
+
+
+!===============================================================================
+! 1.3 Reads 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
+
+!1.3.1 Zonal wind, meridional wind
+
+!Zonal wind
+ierr=NF_INQ_VARID(infid,"u",tmpvarid)
+if (ierr.ne.NF_NOERR) then
+  stop "Error: Failed to get u ID"
+else
+  allocate(u(lonlength,latlength,altlength,timelength))
+  ierr=NF_GET_VAR_REAL(infid,tmpvarid,u)
+  if (ierr.ne.NF_NOERR) then
+    stop "Error: Failed reading zonal wind"
+  endif
+endif
+
+!Meridional wind
+ierr=NF_INQ_VARID(infid,"v",tmpvarid)
+if (ierr.ne.NF_NOERR) then
+  stop "Error: Failed to get v ID"
+else
+  allocate(v(lonlength,latlength,altlength,timelength))
+  ierr=NF_GET_VAR_REAL(infid,tmpvarid,v)
+  if (ierr.ne.NF_NOERR) then
+    stop "Error: Failed reading zonal wind"
+  endif
+endif
+
+
+! 1.3.2 aire
+
+allocate(aire(lonlength,latlength))
+! look for 'aire' in current file
+ierr=NF_INQ_VARID(infid,"aire",tmpvarid) 
+if (ierr.ne.NF_NOERR) then
+  write(*,*) "Warning: Failed to get aire 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
+     infile2="diagfi1.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 these files"
+        stop "Might as well stop here"
+     endif
+  endif
+
+   ! Get ID for aire
+   ierr=NF_INQ_VARID(infid2,"aire",tmpvarid)
+   if (ierr.ne.NF_NOERR) then
+      stop "Error: Failed to get aire ID"
+   endif
+   ! Get aire
+   ierr=NF_GET_VAR_REAL(infid2,tmpvarid,aire)
+   if (ierr.ne.NF_NOERR) then
+      stop "Error: Failed reading aire"
+   endif
+   ! Close file
+   write(*,*) 'OK, got aire'
+   ierr=NF_CLOSE(infid2)
+else
+   ierr=NF_GET_VAR_REAL(infid,tmpvarid,aire)
+   if (ierr.ne.NF_NOERR) then
+      stop "Error: Failed reading aire"
+   endif
+endif
+
+
+! 1.3.3 phisinit
+
+allocate(phisinit(lonlength,latlength))
+! look for '' 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
+     infile2="diagfi1.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 these files"
+        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 phisinit
+   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.3.4 cv
+
+allocate(cv(lonlength,latlength))
+! look for 'cv' in current file
+ierr=NF_INQ_VARID(infid,"cv",tmpvarid) 
+if (ierr.ne.NF_NOERR) then
+  write(*,*) "Warning: Failed to get cv 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
+     infile2="diagfi1.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 these files"
+        stop "Might as well stop here"
+     endif
+  endif
+
+
+   ! Get ID for cv
+   ierr=NF_INQ_VARID(infid2,"cv",tmpvarid)
+   if (ierr.ne.NF_NOERR) then
+      stop "Error: Failed to get cv ID"
+   endif
+   ! Get cv
+   ierr=NF_GET_VAR_REAL(infid2,tmpvarid,cv)
+   if (ierr.ne.NF_NOERR) then
+      stop "Error: Failed reading cv"
+   endif
+   ! Close file
+   write(*,*) 'OK, got cv'
+   ierr=NF_CLOSE(infid2)
+else
+   ierr=NF_GET_VAR_REAL(infid,tmpvarid,cv)
+   if (ierr.ne.NF_NOERR) then
+      stop "Error: Failed reading cv"
+   endif
+endif
+
+
+!Other variables: rho, temp
+
+ierr=NF_INQ_VARID(infid,"temp",tmpvarid)
+if (ierr.ne.NF_NOERR) then
+  stop "Error: Failed to get temp ID"
+else
+  allocate(temp(lonlength,latlength,altlength,timelength))
+  ierr=NF_GET_VAR_REAL(infid,tmpvarid,temp)
+  if (ierr.ne.NF_NOERR) then
+    stop "Error: Failed reading temperature"
+  endif
+endif
+
+ierr=NF_INQ_VARID(infid,"rho",tmpvarid)
+if (ierr.ne.NF_NOERR) then
+  write(*,*) "Error: Failed to get rho ID"
+  write(*,*) "Let's compute it from temperature"
+  allocate(rho(lonlength,latlength,altlength,timelength))
+  
+  do it=1,timelength
+    do ilat=1,latlength
+      do ilon=1,lonlength
+        do ialt=1,altlength
+          rho(ilon,ilat,ialt,it)= (aps(ialt)+bps(ialt)*ps(ilon,ilat,it))/(Rmean*temp(ilon,ilat,ialt,it))
+        enddo
+      enddo
+    enddo
+  enddo
+else
+  allocate(rho(lonlength,latlength,altlength,timelength))
+  ierr=NF_GET_VAR_REAL(infid,tmpvarid,rho)
+  if (ierr.ne.NF_NOERR) then
+    stop "Error: Failed reading mass density"
+  endif
+endif
+
+!==============================================================================
+! 1.4. Output file name
+!==============================================================================
+filename=infile(1:len_trim(infile)-3)//"_stream.nc"
+write(*,*)'The output file has name:'
+write(*,*) filename
+
+
+!==============================================================================
+! 2.1. Temporal averages
+!==============================================================================
+
+allocate(ucum(lonlength,latlength,altlength)) !Time mean zonal wind
+allocate(vcum(lonlength,latlength,altlength)) !Time mean meridional wind
+allocate(tempcum(lonlength,latlength,altlength)) !Time mean temperature
+allocate(rhocum(lonlength,latlength,altlength)) !Time mean density
+allocate(pscum(lonlength,latlength)) !Time mean surface pressure
+do ilat=1,latlength
+   do ilon=1,lonlength
+      pscum(ilon,ilat)=0.
+      do ialt=1,altlength
+            ucum(ilon,ilat,ialt)=0.
+            vcum(ilon,ilat,ialt)=0.
+            tempcum(ilon,ilat,ialt)=0.
+            rhocum(ilon,ilat,ialt)=0.
+            do it=1,timelength
+               ucum(ilon,ilat,ialt)=ucum(ilon,ilat,ialt)+u(ilon,ilat,ialt,it)
+               vcum(ilon,ilat,ialt)=vcum(ilon,ilat,ialt)+v(ilon,ilat,ialt,it)
+               tempcum(ilon,ilat,ialt)=tempcum(ilon,ilat,ialt)+temp(ilon,ilat,ialt,it)
+               rhocum(ilon,ilat,ialt)=rhocum(ilon,ilat,ialt)+rho(ilon,ilat,ialt,it)
+            enddo
+            ucum(ilon,ilat,ialt)=ucum(ilon,ilat,ialt)/timelength
+            vcum(ilon,ilat,ialt)=vcum(ilon,ilat,ialt)/timelength
+            tempcum(ilon,ilat,ialt)=tempcum(ilon,ilat,ialt)/timelength
+            rhocum(ilon,ilat,ialt)=rhocum(ilon,ilat,ialt)/timelength
+      enddo
+      do it=1,timelength
+         pscum(ilon,ilat)=pscum(ilon,ilat)+ps(ilon,ilat,it)
+      enddo
+      pscum(ilon,ilat)=pscum(ilon,ilat)/timelength
+   enddo
+enddo
+
+!==============================================================================
+! 2.2. Calculation of the stream function
+!==============================================================================
+
+!Contravariant meridional wind
+allocate(vcont(lonlength,latlength,altlength))
+do ilon=1,lonlength
+   do ialt=1,altlength
+      do ilat=1,latlength-1
+         vcont(ilon,ilat,ialt)=0.5*&
+              (vcum(ilon,ilat,ialt)+vcum(ilon,ilat+1,ialt))/cv(ilon,ilat)
+      enddo
+      vcont(ilon,latlength,ialt)=0.!vcont(ilon,latlength-1,ialt)
+   enddo
+enddo
+
+!Pressure from hybrid levels
+allocate(Plev(lonlength,latlength,altlength))
+do ilon=1,lonlength
+   do ilat=1,latlength
+      do ialt=1,altlength
+         Plev(ilon,ilat,ialt)=aps(ialt)+bps(ialt)*pscum(ilon,ilat)
+      enddo
+   enddo
+enddo
+
+!Mass in the box
+allocate(mass(lonlength,latlength,altlength))
+do ilon=1,lonlength
+   do ilat=1,latlength
+      do ialt=1,altlength
+         mass(ilon,ilat,ialt)=aire(ilon,ilat)*&
+              Plev(ilon,ilat,ialt)/g0
+      enddo
+   enddo
+enddo
+
+!Mass variation in the box
+allocate(dmass(lonlength,latlength,altlength))
+do ilon=1,lonlength
+   do ilat=1,latlength
+      do ialt=1,altlength-1
+         dmass(ilon,ilat,ialt)=mass(ilon,ilat,ialt)-&
+              mass(ilon,ilat,ialt+1)
+      enddo
+      dmass(ilon,ilat,altlength)=0.!dmass(ilon,ilat,altlength-1)
+   enddo
+enddo
+
+!Stream function
+allocate(psi(latlength,altlength))
+allocate(vcontcum(latlength,altlength))
+do ilat=1,latlength-1
+   psi(ilat,altlength)=0.
+   do ialt=altlength-1,1,-1
+      psi(ilat,ialt)=psi(ilat,ialt+1)
+      vcontcum(ilat,ialt)=0.
+      do ilon=1,lonlength
+         vcontcum(ilat,ialt)=vcontcum(ilat,ialt)+vcont(ilon,ilat,ialt)/lonlength
+         psi(ilat,ialt)=psi(ilat,ialt)-(vcont(ilon,ilat,ialt)&
+              *(dmass(ilon,ilat,ialt)+dmass(ilon,ilat+1,ialt)))
+      enddo
+      psi(latlength,ialt)=0.!psi(latlength-1,ialt)
+   enddo
+enddo
+
+!==============================================================================
+! 2.3. Calculation of the total angular momentum
+!==============================================================================
+
+!Angular momentum
+allocate(mom(lonlength,latlength,altlength))
+allocate(momave(latlength,altlength))
+do ilat=1,latlength
+   do ialt=1,altlength
+      momave(ilat,ialt)=0.
+      do ilon=1,lonlength
+!            mom(ilon,ilat,ialt,it)=dmass(ilon,ilat,ialt,it)*a0*&
+         mom(ilon,ilat,ialt)=a0*cos(lat(ilat)*3.141592/180.)*&
+              (omega*a0*cos(lat(ilat)*3.141592/180.)+ucum(ilon,ilat,ialt))
+         
+         momave(ilat,ialt)=momave(ilat,ialt)+mom(ilon,ilat,ialt)
+      enddo
+      momave(ilat,ialt)=momave(ilat,ialt)/lonlength
+   enddo
+enddo
+
+!==============================================================================
+! 2.4. Zonal mean winds, temperature and density
+!==============================================================================
+
+allocate(uzm(latlength,altlength)) !Zonal mean zonal wind
+allocate(vzm(latlength,altlength)) !Zonal mean meridional wind
+allocate(tempzm(latlength,altlength)) !Zonal mean temperature
+allocate(rhozm(latlength,altlength)) !Zonal mean density
+allocate(pszm(latlength)) !Zonal mean surface pressure
+allocate(phisinitzm(latlength)) !Zonal mean ground geopotential
+do ilat=1,latlength
+   phisinitzm(ilat)=0.
+   pszm(ilat)=0.
+   do ialt=1,altlength
+      uzm(ilat,ialt)=0.
+      vzm(ilat,ialt)=0.
+      tempzm(ilat,ialt)=0.
+      rhozm(ilat,ialt)=0.
+      do ilon=1,lonlength
+         uzm(ilat,ialt)=uzm(ilat,ialt)+ucum(ilon,ilat,ialt)
+         vzm(ilat,ialt)=vzm(ilat,ialt)+vcum(ilon,ilat,ialt)
+         tempzm(ilat,ialt)=tempzm(ilat,ialt)+tempcum(ilon,ilat,ialt)
+         rhozm(ilat,ialt)=rhozm(ilat,ialt)+rhocum(ilon,ilat,ialt)
+      enddo
+      uzm(ilat,ialt)=uzm(ilat,ialt)/lonlength
+      vzm(ilat,ialt)=vzm(ilat,ialt)/lonlength
+      tempzm(ilat,ialt)=tempzm(ilat,ialt)/lonlength
+      rhozm(ilat,ialt)=rhozm(ilat,ialt)/lonlength
+   enddo
+   do ilon=1,lonlength
+      pszm(ilat)=pszm(ilat)+pscum(ilon,ilat)
+      phisinitzm(ilat)=phisinitzm(ilat)+phisinit(ilon,ilat)
+   enddo
+   pszm(ilat)=pszm(ilat)/lonlength
+   phisinitzm(ilat)=phisinitzm(ilat)/lonlength  
+enddo
+
+
+!==============================================================================
+! 2.3. Recalculation of angular momentum using zonally averaged wind
+!==============================================================================
+
+do ilat=1,latlength
+   do ialt=1,altlength
+      momave(ilat,ialt)=a0*cos(lat(ilat)*3.141592/180.)*&
+           (omega*a0*cos(lat(ilat)*3.141592/180.)+uzm(ilat,ialt))
+   enddo
+enddo
+
+
+!==============================================================================
+! 3.1 Dimensions in output file
+!==============================================================================
+
+! 3.1.1 Initialize output file's lat,lon,alt and time dimensions
+call initiate (filename,lat,alt,time,nout,&
+     latdimout,londimout,altdimout,timedimout,lonvarout,timevarout)
+! Initialize output file's aps,bps variables
+call init2(infid,lonlength,latlength,altlength,&
+     nout,londimout,latdimout,altdimout)
+
+! 3.1.2 New longitude dimension/value in output file
+allocate(lon_fake(1))
+lon_fake(1)=0.
+#ifdef NC_DOUBLE
+        ierr= NF_PUT_VARA_DOUBLE(nout,lonvarout,1,1,lon_fake)
+#else
+        ierr= NF_PUT_VARA_REAL(nout,lonvarout,1,1,lon_fake)
+#endif
+
+! 3.1.3 New time dimension/value in output file
+#ifdef NC_DOUBLE
+        ierr= NF_PUT_VARA_DOUBLE(nout,timevarout,1,1,lon_fake)
+#else
+        ierr= NF_PUT_VARA_REAL(nout,timevarout,1,1,lon_fake)
+#endif
+
+
+!==============================================================================
+! 3.2 Write variables
+!==============================================================================
+
+!Define the dimensions of the variables to be written
+
+!!3.2.1 Stream function
+call def_var(nout,"psi","Stream function","",4,&
+             (/londimout,latdimout,altdimout,timedimout/),psivarout,ierr)
+if (ierr.ne.NF_NOERR) then
+   write(*,*) 'Error, could not define variable psi'
+   stop ""
+endif
+
+!Write in the output file
+ierr=NF_PUT_VAR_REAL(nout,psivarout,psi)
+if (ierr.ne.NF_NOERR) then
+   write(*,*)'Error, Failed to write variable psi'
+   stop
+endif
+
+!3.2.2 Momentum
+call def_var(nout,"momave","Angular momentum","",4,&
+             (/londimout,latdimout,altdimout,timedimout/),momvarout,ierr)
+if (ierr.ne.NF_NOERR) then
+   write(*,*) 'Error, could not define variable momave'
+   stop ""
+endif
+
+!Write in the output file
+ierr=NF_PUT_VAR_REAL(nout,momvarout,momave)
+if (ierr.ne.NF_NOERR) then
+   write(*,*)'Error, Failed to write variable momave'
+   stop
+endif
+
+
+!3.2.3 Zonal mean zonal wind
+call def_var(nout,"u","Zonal mean zonal wind","m/s",4,&
+             (/londimout,latdimout,altdimout,timedimout/),uvarout,ierr)
+if (ierr.ne.NF_NOERR) then
+   write(*,*) 'Error, could not define variable u'
+   stop ""
+endif
+
+!Write in the output file
+ierr=NF_PUT_VAR_REAL(nout,uvarout,uzm)
+if (ierr.ne.NF_NOERR) then
+   write(*,*)'Error, Failed to write variable u'
+   stop
+endif
+
+
+!3.2.4 Zonal mean meridional wind
+call def_var(nout,"v","Zonal mean meridional wind","m/s",4,&
+             (/londimout,latdimout,altdimout,timedimout/),vvarout,ierr)
+if (ierr.ne.NF_NOERR) then
+   write(*,*) 'Error, could not define variable v'
+   stop ""
+endif
+
+!Write in the output file
+ierr=NF_PUT_VAR_REAL(nout,vvarout,vzm)
+if (ierr.ne.NF_NOERR) then
+   write(*,*)'Error, Failed to write variable v'
+   stop
+endif
+
+!3.2.5 Zonal mean density
+call def_var(nout,"rho","Zonal mean atmospheric density","",4,&
+             (/londimout,latdimout,altdimout,timedimout/),rhovarout,ierr)
+if (ierr.ne.NF_NOERR) then
+   write(*,*) 'Error, could not define variable rho'
+   stop ""
+endif
+
+!Write in the output file
+ierr=NF_PUT_VAR_REAL(nout,rhovarout,rhozm)
+if (ierr.ne.NF_NOERR) then
+   write(*,*)'Error, Failed to write variable rho'
+   stop
+endif
+
+!3.2.6 Zonal mean temperature
+call def_var(nout,"temp","Zonal mean temperature","K",4,&
+             (/londimout,latdimout,altdimout,timedimout/),tempvarout,ierr)
+if (ierr.ne.NF_NOERR) then
+   write(*,*) 'Error, could not define variable temp'
+   stop ""
+endif
+
+!Write in the output file
+ierr=NF_PUT_VAR_REAL(nout,tempvarout,tempzm)
+if (ierr.ne.NF_NOERR) then
+   write(*,*)'Error, Failed to write variable temp'
+   stop
+endif
+
+!3.2.7 Zonal mean surface pressure
+call def_var(nout,"ps","Zonal mean surface pressure","Pa",3,&
+             (/londimout,latdimout,timedimout/),psvarout,ierr)
+if (ierr.ne.NF_NOERR) then
+   write(*,*) 'Error, could not define variable ps'
+   stop ""
+endif
+
+!Write in the output file
+ierr=NF_PUT_VAR_REAL(nout,psvarout,pszm)
+if (ierr.ne.NF_NOERR) then
+   write(*,*)'Error, Failed to write variable ps'
+   stop
+endif
+
+
+!3.2.8 Zonal mean geopotential
+call def_var(nout,"phisinit","Zonal mean initial geopotential","",2,&
+             (/londimout,latdimout/),phisinitvarout,ierr)
+if (ierr.ne.NF_NOERR) then
+   write(*,*) 'Error, could not define variable phisinit'
+   stop ""
+endif
+
+!Write in the output file
+ierr=NF_PUT_VAR_REAL(nout,phisinitvarout,phisinitzm)
+if (ierr.ne.NF_NOERR) then
+   write(*,*)'Error, Failed to write variable phisinit'
+   stop
+endif
+
+! Close input file
+ierr=nf_close(nid)
+
+! Close output file
+ierr=NF_CLOSE(nout)
+
+contains
+
+!******************************************************************************
+Subroutine initiate (filename,lat,alt,time,&
+                     nout,latdimout,londimout,altdimout,timedimout,&
+                     lonvarout,timevarout)
+!==============================================================================
+! Purpose:
+! Create and initialize a data file (NetCDF format)
+!==============================================================================
+! Remarks:
+! The NetCDF file (created in this subroutine) remains open
+!==============================================================================
+
+implicit none
+
+include "netcdf.inc" ! NetCDF definitions
+
+!==============================================================================
+! Arguments:
+!==============================================================================
+character (len=*), intent(in):: filename
+! filename(): the file's name
+real, dimension(:), intent(in):: lat
+! lat(): latitude
+real, dimension(:), intent(in):: alt
+! alt(): altitude
+real, dimension(:), intent(in):: time
+! time(): Time
+integer, intent(out):: nout
+! nout: [netcdf] file ID
+integer, intent(out):: latdimout
+! latdimout: [netcdf] lat() (i.e.: latitude)  ID
+integer, intent(out):: londimout
+! londimout: [netcdf] lon()  ID
+integer, intent(out):: altdimout
+! altdimout: [netcdf] alt()  ID
+integer, intent(out):: timedimout
+! timedimout: [netcdf] "Time"  ID
+integer, intent(out):: lonvarout
+! timevarout: [netcdf] Longiture (considered as a variable) ID
+integer, intent(out):: timevarout
+! timevarout: [netcdf] Time (considered as a variable) ID
+
+!==============================================================================
+! Local variables:
+!==============================================================================
+!integer :: latdim,londim,altdim,timedim
+integer :: nvarid,ierr
+! nvarid: [netcdf] ID of a variable
+! ierr: [netcdf]  return error code (from called subroutines)
+
+!==============================================================================
+! 1. Create (and open) output file
+!==============================================================================
+write(*,*) "creating "//trim(adjustl(filename))//'...'
+ierr = NF_CREATE(filename,NF_CLOBBER,nout)
+! NB: setting NF_CLOBBER mode means that it's OK to overwrite an existing file
+if (ierr.NE.NF_NOERR) then
+   WRITE(*,*)'ERROR: Impossible to create the file.'
+   stop ""
+endif
+
+!==============================================================================
+! 2. Define/write "dimensions" and get their IDs
+!==============================================================================
+
+ierr = NF_DEF_DIM(nout, "latitude", size(lat), latdimout)
+!ierr = NF_DEF_DIM(nout, "longitude", NF_UNLIMITED, londimout)
+ierr = NF_DEF_DIM(nout, "longitude", 1, londimout)
+ierr = NF_DEF_DIM(nout, "altitude", size(alt), altdimout)
+ierr = NF_DEF_DIM(nout, "Time", 1, timedimout)
+
+! End netcdf define mode
+ierr = NF_ENDDEF(nout)
+
+!==============================================================================
+! 3. Write "Time" (attributes)
+!==============================================================================
+
+call def_var(nout,"Time","Time","years since 0000-00-0 00:00:00",1,&
+             (/timedimout/),timevarout,ierr)
+
+!==============================================================================
+! 4. Write "latitude" (data and attributes)
+!==============================================================================
+
+call def_var(nout,"latitude","latitude","degrees_north",1,&
+             (/latdimout/),nvarid,ierr)
+
+#ifdef NC_DOUBLE
+ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,lat)
+#else
+ierr = NF_PUT_VAR_REAL (nout,nvarid,lat)
+#endif
+
+!==============================================================================
+! 4. Write "longitude" (attributes)
+!==============================================================================
+
+call def_var(nout,"longitude","East longitude","degrees_east",1,&
+             (/londimout/),lonvarout,ierr)
+
+
+!==============================================================================
+! 4. Write "altitude" (data and attributes)
+!==============================================================================
+
+! Switch to netcdf define mode
+
+call def_var(nout,"altitude","Altitude","km",1,&
+             (/altdimout/),nvarid,ierr)
+
+#ifdef NC_DOUBLE
+ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,alt)
+#else
+ierr = NF_PUT_VAR_REAL (nout,nvarid,alt)
+#endif 
+
+end Subroutine initiate
+!******************************************************************************
+subroutine init2(infid,lonlen,latlen,altlen, &
+                 outfid,londimout,latdimout,altdimout)
+!==============================================================================
+! Purpose:
+! Copy aps() and bps() from input file to outpout file
+!==============================================================================
+! Remarks:
+! The NetCDF files must be open
+!==============================================================================
+
+implicit none
+
+include "netcdf.inc" ! NetCDF definitions
+
+!==============================================================================
+! Arguments:
+!==============================================================================
+integer, intent(in) :: infid  ! NetCDF output file ID
+integer, intent(in) :: lonlen ! # of grid points along longitude
+integer, intent(in) :: latlen ! # of grid points along latitude
+integer, intent(in) :: altlen ! # of grid points along latitude
+integer, intent(in) :: outfid ! NetCDF output file ID
+integer, intent(in) :: londimout ! longitude dimension ID
+integer, intent(in) :: latdimout ! latitude dimension ID
+integer, intent(in) :: altdimout ! altitude dimension ID
+!==============================================================================
+! Local variables:
+!==============================================================================
+real,dimension(:),allocatable :: aps,bps ! hybrid vertical coordinates
+integer :: apsid,bpsid
+integer :: ierr
+integer :: tmpvarid ! temporary variable ID
+logical :: aps_ok, bps_ok ! are "phisinit" "aps" "bps" available ?
+
+
+!==============================================================================
+! 1. Read data from input file
+!==============================================================================
+
+! hybrid coordinate aps
+  allocate(aps(altlen))
+ierr=NF_INQ_VARID(infid,"aps",tmpvarid)
+if (ierr.ne.NF_NOERR) then
+  write(*,*) "oops Failed to get aps ID. OK"
+  aps_ok=.false.
+else
+  ierr=NF_GET_VAR_REAL(infid,tmpvarid,aps)
+  if (ierr.ne.NF_NOERR) then
+   stop "error: Failed reading aps"
+  endif
+  aps_ok=.true.
+endif
+
+! hybrid coordinate bps
+  allocate(bps(altlen))
+ierr=NF_INQ_VARID(infid,"bps",tmpvarid)
+if (ierr.ne.NF_NOERR) then
+  write(*,*) "oops: Failed to get bps ID. OK"
+  bps_ok=.false.
+else
+  ierr=NF_GET_VAR_REAL(infid,tmpvarid,bps)
+  if (ierr.ne.NF_NOERR) then
+    stop "Error: Failed reading bps"
+  endif
+  bps_ok=.true.
+endif
+
+!==============================================================================
+! 2. Write
+!==============================================================================
+
+!==============================================================================
+! 2.2. Hybrid coordinates aps() and bps()
+!==============================================================================
+
+IF (aps_ok) then 
+   ! define aps
+   call def_var(nout,"aps","hybrid pressure at midlayers"," ",1,&
+        (/altdimout/),apsid,ierr)
+   if (ierr.ne.NF_NOERR) then
+      stop "Error: Failed to def_var aps"
+   endif
+
+   ! write aps
+#ifdef NC_DOUBLE
+   ierr=NF_PUT_VAR_DOUBLE(outfid,apsid,aps)
+#else
+   ierr=NF_PUT_VAR_REAL(outfid,apsid,aps)
+#endif
+   if (ierr.ne.NF_NOERR) then
+      stop "Error: Failed to write aps"
+   endif
+ENDIF 
+
+IF (bps_ok) then 
+   ! define bps
+   call def_var(nout,"bps","hybrid sigma at midlayers"," ",1,&
+        (/altdimout/),bpsid,ierr)
+   if (ierr.ne.NF_NOERR) then
+      stop "Error: Failed to def_var bps"
+   endif
+
+   ! write bps
+#ifdef NC_DOUBLE
+   ierr=NF_PUT_VAR_DOUBLE(outfid,bpsid,bps)
+#else
+   ierr=NF_PUT_VAR_REAL(outfid,bpsid,bps)
+#endif
+   if (ierr.ne.NF_NOERR) then
+      stop "Error: Failed to write bps"
+   endif
+ENDIF
+
+
+
+! Cleanup
+deallocate(aps)
+deallocate(bps)
+
+end subroutine init2
+!******************************************************************************
+subroutine def_var(nid,name,title,units,nbdim,dim,nvarid,ierr)
+!==============================================================================
+! Purpose: Write a variable (i.e: add a variable to a dataset)
+! called "name"; along with its attributes "title", "units"...
+! to a file (following the NetCDF format)
+!==============================================================================
+! Remarks:
+! The NetCDF file must be open
+!==============================================================================
+
+implicit none
+
+include "netcdf.inc" ! NetCDF definitions
+
+!==============================================================================
+! Arguments:
+!==============================================================================
+integer, intent(in) :: nid
+! nid: [netcdf] file ID #
+character (len=*), intent(in) :: name
+! name(): [netcdf] variable's name
+character (len=*), intent(in) :: title
+! title(): [netcdf] variable's "title" attribute
+character (len=*), intent(in) :: units
+! unit(): [netcdf] variable's "units" attribute
+integer, intent(in) :: nbdim
+! nbdim: number of dimensions of the variable
+integer, dimension(nbdim), intent(in) :: dim
+! dim(nbdim): [netcdf] dimension(s) ID(s)
+integer, intent(out) :: nvarid
+! nvarid: [netcdf] ID # of the variable
+integer, intent(out) :: ierr
+! ierr: [netcdf] subroutines returned error code
+
+! Switch to netcdf define mode
+ierr=NF_REDEF(nid)
+
+! Insert the definition of the variable
+#ifdef NC_DOUBLE
+ierr=NF_DEF_VAR(nid,adjustl(name),NF_DOUBLE,nbdim,dim,nvarid)
+#else
+ierr=NF_DEF_VAR(nid,adjustl(name),NF_FLOAT,nbdim,dim,nvarid)
+#endif
+
+! Write the attributes
+ierr=NF_PUT_ATT_TEXT(nid,nvarid,"title",len_trim(adjustl(title)),adjustl(title))
+ierr=NF_PUT_ATT_TEXT(nid,nvarid,"units",len_trim(adjustl(units)),adjustl(units))
+
+! End netcdf define mode
+ierr=NF_ENDDEF(nid)
+
+end subroutine def_var
+
+end program streamfunction
+
+subroutine init_planet_const(infid)
+! initialize planetary constants using the "controle" array in the file
+! if "controle" array not found in file, look for it in "diagfi.nc"
+use planet_const
+use netcdf
+implicit none
+
+integer,intent(in) :: infid ! input file ID
+
+! local variables
+character(len=100) :: varname
+integer :: varid ! variable ID
+integer :: status
+integer :: infid2 ! for an alternate input file
+real :: controle(100) ! to store "controle" array
+
+! look for "controle" in input file
+infid2=infid ! initialization
+varname="controle" 
+status=nf90_inq_varid(infid,trim(varname),varid)
+if (status.ne.nf90_noerr) then
+  write(*,*) "init_planet_const: Failed to find ",trim(varname)
+  write(*,*) " looking for it in file diagfi.nc"
+  status=nf90_open("diagfi.nc",NF90_NOWRITE,infid2)
+  if (status.ne.nf90_noerr) then
+    write(*,*) " no diafi.nc file ... looking for diagfi1.nc"
+    status=nf90_open("diagfi1.nc",NF90_NOWRITE,infid2)
+    if (status.ne.nf90_noerr) then
+      write(*,*) "might as well stop here..."
+      stop
+    endif
+  endif
+  status=nf90_inq_varid(infid2,trim(varname),varid)
+  if (status.ne.nf90_noerr) then
+    write(*,*) " Failed to find ",trim(varname)," in file!!"
+    stop
+  endif
+  write(*,*) "OK, found ",trim(varname)
+endif
+status=nf90_get_var(infid2,varid,controle)
+if (status.ne.nf90_noerr) then
+  write(*,*) "init_planet_const: Failed to load ",trim(varname)
+  stop
+endif
+
+! now that we have "controle"; extract relevent informations
+a0=controle(5) ! = radius of planet (m)
+g0=controle(7) ! = gravity (m.s-2) at a0
+Rmean=1000.*8.3144598/controle(8) ! controle(8)=mugaz = molar mass (g.mol-1) of atm.
+omega=controle(6) ! rotation rate (rad.s-1)
+
+end subroutine init_planet_const
