program energy ! SL 12/2009: ! This program reads 4D (lon-lat-alt-time) fields directly from LMD outputs without regrid : histmth ! ! it computes: ! dmass -- 4D -- mass of each cell ! sek -- 4D -- specific kinetic energy ! ek -- 1D -- integrated kinetic energy ! sep -- 4D -- specific potential energy ! ep -- 1D -- integrated potential energy ! ! Minimal requirements and dependencies: ! The dataset must include the following data: ! - surface pressure ! - atmospheric temperature ! - zonal and meridional winds ! VERTICAL WIND SPEED IS NEGLECTED IN KINETIC ENERGY implicit none include "netcdf.inc" ! NetCDF definitions character (len=128) :: infile ! input file name (name_P.nc) character (len=128) :: outfile ! output file name character (len=64) :: text ! to store some text integer infid ! NetCDF input file ID 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 :: datashape1d ! shape of 1D datasets integer,dimension(2) :: datashape2d ! shape of 2D datasets integer,dimension(3) :: datashape3d ! shape of 3D datasets integer,dimension(4) :: datashape4d ! shape of 4D datasets real :: miss_val ! special "missing value" to specify missing data real,parameter :: miss_val_def=-9.99e+33 ! default value for "missing value" real :: pi real,dimension(:),allocatable :: lon ! longitude integer lonlength ! # of grid points along longitude real,dimension(:),allocatable :: lat ! latitude real,dimension(:),allocatable :: coslat ! cos of latitude integer latlength ! # of grid points along latitude real,dimension(:),allocatable :: plev ! Pressure levels (Pa) integer altlength ! # of grid point along presnivs (of input datasets) real,dimension(:),allocatable :: time ! time integer timelength ! # of points along time real,dimension(:,:,:),allocatable :: ps ! surface pressure real,dimension(:,:,:,:),allocatable :: temp ! atmospheric temperature real,dimension(:,:,:,:),allocatable :: vitu ! zonal wind (in m/s) real,dimension(:,:,:,:),allocatable :: vitv ! meridional wind (in m/s) real,dimension(:,:,:,:),allocatable :: rayon ! distance to center (m) real,dimension(:,:,:,:),allocatable :: grav ! gravity field (m s-2) real,dimension(:,:,:,:),allocatable :: dmass ! mass in cell (kg) real,dimension(:,:,:,:),allocatable :: sek ! specific kinetic energy real,dimension(:),allocatable :: ek ! total kinetic energy real,dimension(:,:,:,:),allocatable :: sep ! specific potential energy real,dimension(:),allocatable :: ep ! total potential energy integer ierr,ierr1,ierr2 ! NetCDF routines return codes integer i,j,ilon,ilat,ilev,itim ! for loops real :: deltalat,deltalon ! lat and lon intervals in radians real,dimension(:,:,:,:),allocatable :: deltap ! pressure thickness of each layer (Pa) real :: tmpp ! temporary pressure real :: dz ! altitude diff real :: signe ! orientation of lon axis for mountain torque computation logical :: lmdflag real :: cpdet include "planet.h" !=============================================================================== ! 1. Input parameters !=============================================================================== pi = 2.*asin(1.) miss_val = miss_val_def write(*,*) "" write(*,*) "You are working on the atmosphere of ",planet !=============================================================================== ! 1.1 Input file !=============================================================================== write(*,*) "" write(*,*) " Program valid for Venus or Titan LMD, or Venus CAM output 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 grids in lon,lat,alt(pressure),time !=============================================================================== call get_iddim(infid,lat_varid,latlength,lon_varid,lonlength,& alt_varid,altlength,time_varid,timelength,lmdflag ) allocate(lon(lonlength)) ierr=NF_GET_VAR_REAL(infid,lon_varid,lon) if (ierr.ne.NF_NOERR) stop "Error: Failed reading longitude" if(lon(1).gt.lon(2)) then signe=-1. else signe=1. endif allocate(lat(latlength)) ierr=NF_GET_VAR_REAL(infid,lat_varid,lat) if (ierr.ne.NF_NOERR) stop "Error: Failed reading lat" allocate(coslat(latlength)) ! Beware of rounding problems at poles... coslat(:) = max(0.,cos(lat(:)*pi/180.)) ! Lat, lon pressure intervals deltalat = abs(lat(2)-lat(1))*pi/180. deltalon = abs(lon(2)-lon(1))*pi/180. allocate(plev(altlength)) ierr=NF_GET_VAR_REAL(infid,alt_varid,plev) if (ierr.ne.NF_NOERR) stop "Error: Failed reading presnivs (ie pressure levels)" if(.not.lmdflag) then ! in CAM files, pressure in mbar and reversed... call reverselev(altlength,plev) plev=plev*100. ! convertion to Pa endif allocate(time(timelength)) ierr=NF_GET_VAR_REAL(infid,time_varid,time) if (ierr.ne.NF_NOERR) stop "Error: Failed reading time" ! Time axis IN PLANET DAYS if(.not.lmdflag) then ! in CAM files, time in Earth days... ! => seconds time=time*86400. endif time=time/localday !=============================================================================== ! 1.3 Get output file name !=============================================================================== write(*,*) "" !write(*,*) "Enter output file name" !read(*,*) outfile outfile=infile(1:len_trim(infile)-3)//"_NRG.nc" write(*,*) "Output file name is: "//trim(outfile) !=============================================================================== ! 2.1 Store needed fields !=============================================================================== !=============================================================================== ! 2.1.1 Surface pressure !=============================================================================== allocate(ps(lonlength,latlength,timelength)) text="ps" call get_var3d(infid,lonlength,latlength,timelength,text,ps,ierr1,ierr2) if (ierr1.ne.NF_NOERR) then write(*,*) " looking for psol instead... " text="psol" call get_var3d(infid,lonlength,latlength,timelength,text,ps,ierr1,ierr2) if (ierr1.ne.NF_NOERR) stop "Error: Failed to get psol ID" endif if (ierr2.ne.NF_NOERR) stop "Error: Failed reading surface pressure" if((.not.lmdflag).and.(planet.eq."Venus")) call reverse3d(lonlength,latlength,timelength,ps) !=============================================================================== ! 2.1.2 Atmospheric temperature !=============================================================================== allocate(temp(lonlength,latlength,altlength,timelength)) if(lmdflag) then text="temp" else text="T" endif call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,miss_val,ierr1,ierr2) if (ierr1.ne.NF_NOERR) then write(*,*) " looking for t instead... " text="t" call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,miss_val,ierr1,ierr2) if (ierr1.ne.NF_NOERR) stop "Error: Failed to get temperature ID" endif if (ierr2.ne.NF_NOERR) stop "Error: Failed reading temperature" if(.not.lmdflag) call reverse4d(lonlength,latlength,altlength,timelength,temp) !=============================================================================== ! 2.1.3 Winds !=============================================================================== allocate(vitu(lonlength,latlength,altlength,timelength)) allocate(vitv(lonlength,latlength,altlength,timelength)) ! zonal wind vitu (in m/s) if(lmdflag) then text="vitu" else text="U" endif call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,miss_val,ierr1,ierr2) if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitu ID" if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zonal wind" if(.not.lmdflag) call reverse4d(lonlength,latlength,altlength,timelength,vitu) ! meridional wind vitv (in m/s) if(lmdflag) then text="vitv" else text="V" endif call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,miss_val,ierr1,ierr2) if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitv ID" if (ierr2.ne.NF_NOERR) stop "Error: Failed reading meridional wind" if(.not.lmdflag) call reverse4d(lonlength,latlength,altlength,timelength,vitv) !=============================================================================== ! 2.1.4 Altitude above areoide !=============================================================================== ! Only needed if g(z) on Titan... !allocate(za(lonlength,latlength,altlength,timelength)) !text="zareoid" !call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,miss_val,ierr1,ierr2) !if (ierr1.ne.NF_NOERR) stop "Error: Failed to get za ID" !if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zareoid" !=============================================================================== ! 2.2 Computations !=============================================================================== !=============================================================================== ! 2.2.2 Mass in cells !=============================================================================== allocate(rayon(lonlength,latlength,altlength,timelength)) allocate(grav(lonlength,latlength,altlength,timelength)) allocate(dmass(lonlength,latlength,altlength,timelength)) do itim=1,timelength do ilon=1,lonlength do ilat=1,latlength do ilev=1,altlength ! Need to be consistent with GCM computations ! if (za(ilon,ilat,ilev,itim).ne.miss_val) then rayon(ilon,ilat,ilev,itim) = a0 ! rayon(ilon,ilat,ilev,itim) = za(ilon,ilat,ilev,itim) + a0 grav(ilon,ilat,ilev,itim) = g0*a0*a0 & /(rayon(ilon,ilat,ilev,itim)*rayon(ilon,ilat,ilev,itim)) ! else ! rayon(ilon,ilat,ilev,itim) = miss_val ! grav(ilon,ilat,ilev,itim) = miss_val ! endif enddo enddo enddo enddo ! timelength call cellmass(infid,latlength,lonlength,altlength,timelength,lmdflag, & miss_val,deltalon,deltalat,coslat,plev,ps,grav,rayon, & dmass ) !=============================================================================== ! 2.2.6 Specific energies !=============================================================================== allocate(sek(lonlength,latlength,altlength,timelength)) allocate(sep(lonlength,latlength,altlength,timelength)) do itim=1,timelength do ilon=1,lonlength do ilat=1,latlength do ilev=1,altlength if (rayon(ilon,ilat,ilev,itim).ne.miss_val) then if ((vitu(ilon,ilat,ilev,itim).lt.miss_val) & .and.(vitv(ilon,ilat,ilev,itim).lt.miss_val)) then sek(ilon,ilat,ilev,itim) = 0.5 * & ( vitu(ilon,ilat,ilev,itim)*vitu(ilon,ilat,ilev,itim) & + vitv(ilon,ilat,ilev,itim)*vitv(ilon,ilat,ilev,itim) ) else sek(ilon,ilat,ilev,itim) = miss_val endif if (temp(ilon,ilat,ilev,itim).ne.miss_val) then sep(ilon,ilat,ilev,itim) = temp(ilon,ilat,ilev,itim) & * cpdet(temp(ilon,ilat,ilev,itim)) else sep(ilon,ilat,ilev,itim) = miss_val endif else sek(ilon,ilat,ilev,itim) = miss_val sep(ilon,ilat,ilev,itim) = miss_val endif enddo enddo enddo enddo ! timelength !=============================================================================== ! 2.2.7 Total energies !=============================================================================== allocate(ek(timelength)) allocate(ep(timelength)) do itim=1,timelength ek(itim) = 0. ep(itim) = 0. do ilon=1,lonlength do ilat=1,latlength do ilev=1,altlength if (sek(ilon,ilat,ilev,itim).ne.miss_val) then ek(itim) = ek(itim) & + sek(ilon,ilat,ilev,itim) * dmass(ilon,ilat,ilev,itim) endif if (sep(ilon,ilat,ilev,itim).ne.miss_val) then ep(itim) = ep(itim) & + sep(ilon,ilat,ilev,itim) * dmass(ilon,ilat,ilev,itim) endif enddo enddo enddo if (ek(itim).eq.0.) then ek(itim) = miss_val ep(itim) = miss_val endif enddo ! timelength print*,"End of computations" !=============================================================================== ! 3. Create 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.1. Define and write dimensions !=============================================================================== call write_dim(outfid,lonlength,latlength,altlength,timelength, & lon,lat,plev,time,lon_dimid,lat_dimid,alt_dimid,time_dimid) !=============================================================================== ! 3.2. Define and write variables !=============================================================================== ! 1D Variables datashape1d=time_dimid call write_var1d(outfid,datashape1d,timelength,& "ek ", "Total kinetic energy","J ",miss_val,& ek ) call write_var1d(outfid,datashape1d,timelength,& "ep ", "Total pot energy ","J ",miss_val,& ep ) ! 3D variables datashape3d(1)=lon_dimid datashape3d(2)=lat_dimid datashape3d(3)=time_dimid call write_var3d(outfid,datashape3d,lonlength,latlength,timelength,& "ps ", "Surface pressure ","Pa ",miss_val,& ps ) ! 4D variables datashape4d(1)=lon_dimid datashape4d(2)=lat_dimid datashape4d(3)=alt_dimid datashape4d(4)=time_dimid call write_var4d(outfid,datashape4d,lonlength,latlength,altlength,timelength,& "dmass ", "Mass ","kg ",miss_val,& dmass ) call write_var4d(outfid,datashape4d,lonlength,latlength,altlength,timelength,& "sek ", "Specific kin energy ","J/kg ",miss_val,& sek ) call write_var4d(outfid,datashape4d,lonlength,latlength,altlength,timelength,& "sep ", "Specific pot energy ","J/kg ",miss_val,& sep ) !!!! Close output file ierr=NF_CLOSE(outfid) if (ierr.ne.NF_NOERR) then write(*,*) 'Error, failed to close output file ',outfile endif end program