Changeset 2432 for trunk


Ignore:
Timestamp:
Nov 14, 2020, 3:27:44 PM (4 years ago)
Author:
emillour
Message:

Mars GCM:

  • Make a single README file (get rid of REAME.exec) and update "compile" script to (hopefully) better illustrate how to compile utilities
  • Update utilities zrecast and lslin: zrecast.F90: Force the vertical interpolation of any variable starting by "rho" to be in log space (as it should for density, molecular concentration in mol.cm-3 etc..). Set missing value to 1.E20 lslin.F90: fix various problems in the output.

FF + EM

Location:
trunk/LMDZ.MARS
Files:
4 edited
1 moved

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/README

    r2429 r2432  
    31983198distance.
    31993199
     3200== 14/11/2020 == FF + EM
     3201- Make a single README file (get rid of REAME.exec) and update "compile"
     3202  script to (hopefully) better illustrate how to compile utilities
     3203- Update utilities zrecast and lslin:
     3204  zrecast.F90: Force the vertical interpolation of any variable starting
     3205  by "rho" to be in log space (as it should for density, molecular
     3206  concentration in mol.cm-3 etc..). Set missing value to 1.E20
     3207  lslin.F90: fix various problems in the output.

  • trunk/LMDZ.MARS/util/README

    r2431 r2432  
    1 This directory contains executable files (and their source) that can be
     1This directory contains source codes of utilities (programs) that can be
    22used to process the LMD Mars GCM output files (like "diagfi" or
    3 "stats"), whatever the dimension.
    4 In addition most output file from one of these program can be
    5 processes by another.
     3"stats").
     4In addition most output file from one of these program can then be
     5processed by another.
    66
    7 The executable files should run on any Linux platform.
     7See script "compile" for instructions and recommendations to compile
     8the utilities. Once adequately modified, that script can be used to
     9compile each utility; e.g. to build the zrecast.e utility program:
     10> compile zrecast
    811
    9 Inputs can be provided by (1) replying to questions on screen or
    10 (2) filling the corresponding  *.def and direct the input of these *.def
     12Note that inputs can be provided by (1) replying to questions on screen or
     13(2) filling the corresponding  *.def and redirecting the input of these *.def
    1114files instead. For example :
    1215> concatnc.e < concatnc.def
    … …  
    1417
    1518--------------------------------------------------------------------
    16 1) concatnc.e
     191) concatnc.e : concatenate successive GCM output files
    1720--------------------------------------------------------------------
    1821
    … …  
    2326software like Grads may not be able to use it. To have a linear "Ls"
    2427timescale, you can use "Ls_Linear.e" (see below).
     28You can also add a "Ls" variable corresponding solar longitude for each
     29timestep by using the "adls" option.
    2530
    2631Output file is : concat.nc
    2732
    28 MODIFICATION:
    29 07/2008 Utility concatnc.F90 (not used by the gcm): improvement in order
    30 to 1) concatenate 1D variable and 2) increase the number of input files
    31 up to 1000
    32 
    33 
    3433--------------------------------------------------------------------
    35 2) localtime.e
     342) zrecast.e : put GCM data in physical vertical coordinate
    3635--------------------------------------------------------------------
    3736
    38 Program to redistribute and interpolate the variable a the same
    39 local times everywhere (useful to mimic satellite observations, or
    40 analyse day to day variations at a given local time).
    41 input : diagfi.nc  / concat.nc / stats.nc kind of files
    42 
    43 output file is :
    44 name_of_input_file_LT.nc with pressure coordinate
    45 
    46 --------------------------------------------------------------------
    47 3) zrecast.e
    48 --------------------------------------------------------------------
     37GCM outputs are in GCM hybrid coordinate which do not correspond to any
     38physical vertical coordinate ! zrecast is NECESSARY to make any publishable
     39scientific figure.
    4940
    5041This program reads 4D (lon-lat-alt-time) fields from GCM output files
    … …  
    6556- surface pressure
    6657- atmospheric temperature
    67 - hybrid coordinates aps() and bps(), or sigma levels() (see section
    68 1.3.2)
     58- hybrid coordinates aps() and bps(), or sigma levels()
    6959- ground geopotential (in input file; if not found, it is sought
    7060  in a 'diagfi.nc' file. If not found there, it is then sought in
    … …  
    77672.2 in source file)
    7868
     69- Vertical interpolation : note that density (kg.m-3) and species density
     70(e.g. in molecules.cm-3) must be vertically interpolated in log-space. For
     71that purpose, these variables names must start with "rho" (exemple : rho for
     72density, rho_co2, rho_o2, etc....)
     73
    7974output file is :
    8075name_of_input_file_P.nc with pressure coordinate
    … …  
    8378name_of_input_file_R.nc with altitude as distance to center of planet
    8479
    85 MODIFICATION :
    86 01/2010 : correction to interpolate above surface if density is not available.
    87 03/2011 : added possibility to have output as distance to center of planet
     80--------------------------------------------------------------------
     813) localtime.e : put GCM data in local time coordinate
     82--------------------------------------------------------------------
     83
     84Program to redistribute and interpolate the variable a the same
     85local times everywhere (useful to mimic satellite observations, or
     86analyse day to day variations at a given local time).
     87input : diagfi.nc  / concat.nc / stats.nc kind of files
     88
     89output file is : name_of_input_file_LT.nc
    8890
    8991--------------------------------------------------------------------
    90 4) lslin.e
     924) lslin.e : redistribute and AVERAGE gcm output in Ls coordinate.
    9193--------------------------------------------------------------------
    9294
    9395This program has been designed to interpol data in Solar Longitude (Ls)
    9496linear time coordinate (usable with grads) from Netcdf diagfi or concatnc 
    95 files.
    96 output file is : lslin.nc
    97 lslin also create a lslin.ctl file that can be read
     97files and to AVERAGE the instantaneous data in Ls period (for comparison with
     98binned dataset, for instance).
     99
     100lslin also creates a lslin.ctl file that can be read
    98101directly by grads (>xdfopen lsllin.ctl) to plot in Ls coordinate to
    99102avoid some problem with grads when grads think that "the time interval
    100103is too small"...
    101104
    102 MODIFICATION
    103 10/2007 Utility lslin.F90 (not used by the gcm)
    104 changed evaluation of 'start_var' from hard-coded values to a computed value
    105 04/2015 Added possibility to bin data (instead of interpolating) over
    106 the time intervals
     105Output file is : name_of_input_file_LS.nc
    107106
    108107--------------------------------------------------------------------
    109 5) hrecast.e
     1085) solzenangle.e : select GCM data at a given solar zenith angle
    110109--------------------------------------------------------------------
    111110
    112 This program can interpolate GCM output on any horizontal grid (regular lat - lon) as long as it cover all the
    113 planet. The grid can be given points by points. The best way is to use the redirected input hrecast.def
     111Program to redistribute and interpolate the variable a the same
     112solar zenith angle (notably useful to mimic satellite observations and in particular
     113SOLAR OCCULTATIONS by choosing solar zenith angle = 90°)
     114The user choose between Morning and Evening side.
     115
     116input : diagfi.nc  / concat.nc / stats.nc kind of files
     117
     118On the morning side output file (1 output per sol per grid point) is :
     119 name_of_input_file_MO.nc 
     120
     121On the evening side output file is (1 output per sol per grid point):
     122 name_of_input_file_EV.nc 
     123
     124--------------------------------------------------------------------
     1256) hrecast.e : interpolate data at another horizontal grid resolution.
     126--------------------------------------------------------------------
     127
     128This program can interpolate GCM output on any horizontal grid (regular
     129lat - lon) as long as it cover all the planet.
     130The grid can be given points by points. The best way is to use the
     131redirected input hrecast.def
    114132
    115133hrecast.e < hrecast.def
    … …  
    117135
    118136--------------------------------------------------------------------
    119 6) expandstartfi.e
     1377) expandstartfi.e : to plot data in a startfi.nc file
    120138--------------------------------------------------------------------
    121139
     140In startfi.nc file, data are not plotable because the horizontal coordinate
     141in only a 1D list of all atmospheric columns on the planets.
    122142This program takes a physics start file ("startfi.nc") and recasts it
    123143on the corresponding  lonxlat grid (so it contents may easily be displayed
    124 using Grads, Ferret, etc.)
     144using Python, Grads, Ferret, etc.)
    125145
    126146Simply run expandstartfi.e as a command line with arguments:
    … …  
    132152
    133153--------------------------------------------------------------------
    134 6) extract.e
     1548) extract.e : get data at specific coordinates for comparison with observations
    135155--------------------------------------------------------------------
    136156

  • trunk/LMDZ.MARS/util/compile

    r2395 r2432  
    33# > compile concat
    44# > compile zrecast
     5## BUT first you must customize this script to your personal settings:
     6# 1) set up the correct environment; e.g. environment variable
     7#    NETCDF_HOME should point to your NetCDF distribution root directory
     8#    (and possibly you might need to "module load ..." a few things)
     9# 2) put the appropriate compiler and compiler options
     10#    in variables COMPILER and COMPILER_OPTIONS
     11# 3) Note that when you will run the executable, you might need to
     12#    also add the paths to the used libraries (e.g. $NETCDF_HOME/lib)
     13#    in environment variable LD_LIBRARY_PATH (most often the "module load ..."
     14#    command does this, so you should run it before running the executable)
    515
    6 # pgf90 -Bstatic   $1.F90 \
    7 #-I/distrib/local/netcdf/pgi_7.1-6_32/include \
    8 #-L/distrib/local/netcdf/pgi_7.1-6_32/lib -lnetcdf  -o $1.e
     16# Setup: (see at the end of this script for real world examples)
     17# possibly source some modules here and adapt variables below:
     18NETCDF_HOME="/path/to/the/NetCDF/root/directory"
     19COMPILER="gfortran"
     20COMPILER_OPTIONS="-O2"
    921
    10 #ifort $1.F90 \
    11 #pgf90 $1.F90 \
    12 ifort $1.F90 \
    13 -I$NETCDF/include \
    14 -L$NETCDF/lib -lnetcdf -o $1.e
     22# Compilation:
     23# (on some very old systems the Fortran NetCDF library is included
     24#  in the C library and "-lnetcdff" should be replaced with "-lnetcdf")
    1525
    16 # Before running that on you computer you might want to change :
    17 # 1) replace "pgf90" with the name of your favorite compiler
    18 #    (you may also add some non-agressive optimization options e.g. -O2)
    19 # 2) replace "/distrib/local/netcdf/pgi_7.1-6_32/lib" with the
    20 # address of the
    21 # directory that contains the NetCDF library (file libnetcdf.a that can
    22 # be obtained for free on
    23 # http://www.unidata.ucar.edu/packages/netcdf/index.html
    24 # (see user manual)
     26$COMPILER $COMPILER_OPTIONS $1.F90 \
     27-I$NETCDF_HOME/include \
     28-L$NETCDF_HOME/lib -lnetcdff \
     29-o $1.e
     30
    2531#
    26 # 3) Replace "/distrib/local/netcdf/pgi_7.1-6_32/lib" with the address of the
    27 # directory that contains the NetCDF  include file "netcdf.inc"
    28 # that can be obtained at the web address above.
    29 #
    30 # 4) The "-Bstatic" option is here to ensure that the executable will
    31 # work on any Linux machine (only necessary if you want to export the
    32 # executable from a machine to another).
     32# Example of a setup on a simple Linux system where the netcdf library
     33# is in a personal location /home/myacount/netcdf directory:
     34# NETCDF_HOME=/home/myaccount/netcdf
     35# COMPILER="gfortran"
     36# COMPILER_OPTIONS="-O2"
     37# And of course the LD_LIBRARY_PATH environement variable should contain
     38# path "/home/myaccount/netcdf/lib" to be able to run the executable
     39#
     40# Example of a setup on LMD CentOS7 machines using gfortran and NetCDF 4.5:
     41# module purge
     42# module load gnu/7.2.0
     43# module load netcdf4/4.5.0-gfortran72
     44# NETCDF_HOME=/opt/netcdf45/gfortran72
     45# COMPILER="gfortran"
     46# COMPILER_OPTIONS="-O2"
     47# And of course modules above need be loaded before running the executable
     48#
     49# Example of a setup on the Ciclad cluster using ifort and NetCDF 4.3
     50# module purge
     51# module load intel/15.0.6.233
     52# module load netcdf4/4.3.3.1-ifort
     53# NETCDF_HOME=/opt/netcdf43/ifort
     54# COMPILER="ifort"
     55# COMPILER_OPTIONS="-O2 -ip"
     56# And of course modules above need be loaded before running the executable
    3357
    3458
    35 
    36 

  • trunk/LMDZ.MARS/util/lslin.F90

    r1409 r2432  
    2727character (len=50) :: units_alt
    2828! var(): units of altitude which can be 'km' or 'Pa' (after zrecast)
     29character (len=50) :: altlong_name,altunits,altpositive
     30! altlong_name(): [netcdf] altitude "long_name" attribute
     31! altunits(): [netcdf] altitude "units" attribute
     32! altpositive(): [netcdf] altitude "positive" attribute
     33
    2934character (len=100) :: outfile
    3035! outfile(): output file name
    … …  
    3540character (len=1) :: average
    3641! average: the character 'y' to average within the Ls timestep
    37 integer :: nid,varid,ierr,miss
     42integer :: nid,varid,ierr,miss,validr
    3843! nid: [netcdf] ID # of input file
    3944! varid: [netcdf] ID # of a variable
    … …  
    5560! tab_cntrl(length): array, stored in the input file,
    5661!                which contains various control parameters.
    57 real, dimension(:), allocatable:: lat,lon,alt,time,lsgcm,timels
     62real, dimension(:), allocatable:: lat,lon,alt,time,lsgcm,timels,ctl
    5863! lat(): latitude coordinates (read from input file)
    5964! lon(): longitude coordinates (read from input file)
    … …  
    6267! lsgcm(): time coordinate (in unevenly spaced "Ls")
    6368! timels(): new time coordinates (evenly spaced "Ls"; written to output file)
    64 integer :: latlen,lonlen,altlen,timelen,Nls,Sls
     69! ctl(): array, stores controle array
     70integer :: latlen,lonlen,altlen,timelen,Nls,Sls,ctllen
    6571! latlen: # of elements of lat() array
    6672! lonlen: # of elements of lon() array
    … …  
    7278! nbvarfile: total # of variables (in netcdf file)
    7379! ndim: [netcdf] # of dimensions (3 or 4) of a variable
    74 integer :: latdim,londim,altdim,timedim
     80integer :: latdim,londim,altdim,timedim,ctldim
    7581! latdim: [netcdf] "latitude" dim ID
    7682! londim: [netcdf] "longitude" dim ID
    7783! altdim: [netcdf] "altdim" dim ID
    7884! timedim: [netcdf] "timedim" dim ID
    79 integer :: latvar,lonvar,altvar,timevar
     85! ctldim: [netcdf] "controle" dim ID
     86integer :: latvar,lonvar,altvar,timevar,ctlvar
    8087! latvar: [netcdf] ID of "latitude" variable
    8188! lonvar: [netcdf] ID of "longitude" variable
    … …  
    240247   else
    241248      ierr=NF_INQ_DIMLEN(nid,altdim,altlen)
    242       units_alt="                                                    "
     249      units_alt=""
    243250      ierr=NF_GET_ATT_TEXT(nid,varid,"units",units_alt)
    244251   endif
    245252   write(*,*) "altlen: ",altlen
     253! Get altitude attributes to handle files with any altitude type
     254   ierr=nf_get_att_text(nid,altvar,'long_name',altlong_name)
     255   ierr=nf_get_att_text(nid,altvar,'units',altunits)
     256   ierr=nf_get_att_text(nid,altvar,'positive',altpositive)
     257
    246258
    247259! Allocate lat(), lon() and alt()
    … …  
    251263
    252264! Read lat(),lon() and alt() from input file
    253 
    254 
    255 
    256 
    257 
    258265      ierr = NF_GET_VAR_REAL(nid,latvar,lat)
    259266      ierr = NF_GET_VAR_REAL(nid,lonvar,lon)
    … …  
    264271
    265272
     273   ierr=NF_INQ_DIMID(nid,"index",ctldim)
     274   if (ierr.NE.NF_NOERR) then
     275      write(*,*) ' Dimension <index> is missing in file '//trim(infile)
     276      ctllen=0
     277      !stop "" 
     278   endif
     279   ierr=NF_INQ_VARID(nid,"controle",ctlvar)
     280   if (ierr.NE.NF_NOERR) then
     281      write(*,*) 'Field <controle> is missing in file '//trim(infile)
     282      ctllen=0
     283      !stop ""
     284   else
     285      ierr=NF_INQ_DIMLEN(nid,ctldim,ctllen)
     286   endif
     287
     288   allocate(ctl(ctllen),stat=ierr)
     289   if (ierr.ne.0) then
     290     write(*,*) "Error: failed to allocate ctl(ctllen)"
     291     stop
     292   endif
     293
     294   if (ctllen .ne. 0) then
     295     ierr = NF_GET_VAR_REAL(nid,ctlvar,ctl)
     296     if (ierr.ne.0) then
     297       write(*,*) "Error: failed to load controle"
     298       write(*,*) NF_STRERROR(ierr)
     299       stop
     300     endif
     301   endif ! of if (ctllen .ne. 0)
     302
     303
     304
    266305!==============================================================================
    267306! 2.2. Create (and initialize latitude, longitude and altitude in) output file
    268307!==============================================================================
    269308
    270 call initiate(outfile,lat,lon,alt,nout,&
    271            latdimout,londimout,altdimout,timedimout,timevarout,units_alt)
     309  ! Initialize output file's lat,lon,alt and time dimensions
     310   call initiate(outfile,lat,lon,alt,ctl,latlen,lonlen,altlen,ctllen,&
     311         altlong_name,altunits,altpositive,&
     312         nout,latdimout,londimout,altdimout,timedimout,timevarout)
     313
     314  ! Initialize output file's aps,bps and phisinit variables
     315   call init2(nid,lonlen,latlen,altlen,altdim,&
     316              nout,londimout,latdimout,altdimout)
    272317
    273318!==============================================================================
    … …  
    328373   endif
    329374
    330    day_ini = tab_cntrl(4)                                               
     375   day_ini = nint(tab_cntrl(4))                                               
    331376   day_ini = modulo(day_ini,669)                                                   
    332377   write(*,*) 'day_ini', day_ini
    … …  
    392437
    393438! Build timels()
    394 timels(1) = 0.01*nint(100.*ls0) ! Ehouarn: what the !!!???!!!
     439timels(1) = 0.01*nint(100.*ls0)
    395440do k=2,Nls
    396441   timels(k) = timels(k-1) + deltals
    … …  
    482527! 3.3 Write this variable's definition and attributes to the output file
    483528!==============================================================================
    484    units="                                                    "
    485    title="                                                    "
    486    ierr=NF_GET_ATT_TEXT(nid,varid,"title",title)
    487    ierr=NF_GET_ATT_TEXT(nid,varid,"units",units)
    488 
    489    call def_var(nout,var(j),title,units,ndim,dim,varidout,ierr)
     529   if (ndim.ge.3) then
     530    units=""
     531    title=""
     532    ierr=NF_GET_ATT_TEXT(nid,varid,"title",title)
     533    ierr=NF_GET_ATT_TEXT(nid,varid,"units",units)
     534
     535    call def_var(nout,var(j),title,units,ndim,dim,varidout,ierr)
     536   end if
    490537
    491538!==============================================================================
    492539! 3.4 Read variable
    493 !==============================================================================
    494 
    495 
    496 
    497 
    498    ierr = NF_GET_VAR_REAL(nid,varid,var3d)
    499    miss = NF_GET_ATT_REAL(nid,varid,"missing_value",missing)
    500    miss = NF_GET_ATT_REAL(nid,varid,"valid_range",valid_range)
     540!==== ==========================================================================
     541
     542   if (ndim.ge.3) then
     543      ierr = NF_GET_VAR_REAL(nid,varid,var3d)
     544      miss=NF_GET_ATT_REAL(nid,varid,"missing_value",missing)
     545      validr=NF_GET_ATT_REAL(nid,varid,"valid_range",valid_range)
     546   end if
    501547
    502548
    … …  
    582628
    583629
    584 
    585    ierr= NF_PUT_VARA_REAL(nout,varidout,corner,edges,var3dls)
    586 
    587 
    588    if (ierr.ne.NF_NOERR) then
     630   if (ndim.ge.3) then
     631
     632     ierr= NF_PUT_VARA_REAL(nout,varidout,corner,edges,var3dls)
     633
     634
     635     if (ierr.ne.NF_NOERR) then
    589636      write(*,*) 'PUT_VAR ERROR: ',NF_STRERROR(ierr)
    590637      stop ""
    591    endif
     638     endif
    592639
    593640! In case there is a "missing value" attribute and "valid range"
    594    if (miss.eq.NF_NOERR) then
    595       call missing_value(nout,varidout,valid_range,missing) 
    596    end if
    597 
    598    deallocate(var3d)
    599    deallocate(var3dls)
    600    deallocate(var3dxy)
     641     if (miss.eq.NF_NOERR) then
     642      call missing_value(nout,varidout,missing) 
     643     end if
     644
     645     deallocate(var3d)
     646     deallocate(var3dls)
     647     deallocate(var3dxy)
     648   endif !  if (ndim.ge.3)
    601649
    602650enddo ! of do j=1,nbvar loop
    … …  
    641689contains
    642690
    643 !************************************
    644 subroutine initiate (filename,lat,lon,alt,&
    645                      nout,latdimout,londimout,altdimout,timedimout,timevarout,units_alt)
     691!******************************************************************************
     692
     693subroutine initiate(filename,lat,lon,alt,ctl,latlen,lonlen,altlen,ctllen,&
     694                     altlong_name,altunits,altpositive,&
     695                     nout,latdimout,londimout,altdimout,timedimout,timevarout)
    646696!==============================================================================
    647697! Purpose:
    … …  
    661711character (len=*), intent(in):: filename
    662712! filename(): the file's name
    663 real, dimension(:), intent(in):: lat
     713integer,intent(in) ::latlen,lonlen,altlen,ctllen
     714real, intent(in):: lat(latlen)
    664715! lat(): latitude
    665 real, dimension(:), intent(in):: lon
     716real, intent(in):: lon(lonlen)
    666717! lon(): longitude
    667 real, dimension(:), intent(in):: alt
     718real, intent(in):: alt(altlen)
    668719! alt(): altitude
     720real, intent(in):: ctl(ctllen)
     721! ctl(): controle
     722character (len=*), intent(in) :: altlong_name
     723! altlong_name(): [netcdf] altitude "long_name" attribute
     724character (len=*), intent(in) :: altunits
     725! altunits(): [netcdf] altitude "units" attribute
     726character (len=*), intent(in) :: altpositive
     727! altpositive(): [netcdf] altitude "positive" attribute
    669728integer, intent(out):: nout
    670729! nout: [netcdf] file ID
    … …  
    679738integer, intent(out):: timevarout
    680739! timevarout: [netcdf] "Time" (considered as a variable) ID
    681 character (len=50) :: units_alt
    682 ! var(): units of altitude which can be m or Pa (after zrecast)
    683 
    684740
    685741!==============================================================================
    … …  
    688744!integer :: latdim,londim,altdim,timedim
    689745integer :: nvarid,ierr
     746integer :: ctldimout
    690747! nvarid: [netcdf] ID of a variable
    691748! ierr: [netcdf]  return error code (from called subroutines)
    … …  
    698755! NB: setting NF_CLOBBER mode means that it's OK to overwrite an existing file
    699756if (ierr.NE.NF_NOERR) then
    700    WRITE(*,*)'ERROR: Impossible to create the file.'
     757   WRITE(*,*)'ERROR: Impossible to create the file ',trim(filename)
     758   write(*,*) NF_STRERROR(ierr)
    701759   stop ""
    702760endif
    … …  
    706764!==============================================================================
    707765
    708 ierr = NF_DEF_DIM(nout, "latitude", size(lat), latdimout)
    709 ierr = NF_DEF_DIM(nout, "longitude", size(lon), londimout)
    710 ierr = NF_DEF_DIM(nout, "altitude", size(alt), altdimout)
     766ierr = NF_DEF_DIM(nout, "latitude", latlen, latdimout)
     767if (ierr.NE.NF_NOERR) then
     768   WRITE(*,*)'initiate: error failed to define dimension <latitude>'
     769   write(*,*) NF_STRERROR(ierr)
     770   stop ""
     771endif
     772ierr = NF_DEF_DIM(nout, "longitude", lonlen, londimout)
     773if (ierr.NE.NF_NOERR) then
     774   WRITE(*,*)'initiate: error failed to define dimension <longitude>'
     775   write(*,*) NF_STRERROR(ierr)
     776   stop ""
     777endif
     778ierr = NF_DEF_DIM(nout, "altitude", altlen, altdimout)
     779if (ierr.NE.NF_NOERR) then
     780   WRITE(*,*)'initiate: error failed to define dimension <altitude>'
     781   write(*,*) NF_STRERROR(ierr)
     782   stop ""
     783endif
     784if (size(ctl).ne.0) then
     785  ierr = NF_DEF_DIM(nout, "index", ctllen, ctldimout)
     786  if (ierr.NE.NF_NOERR) then
     787    WRITE(*,*)'initiate: error failed to define dimension <index>'
     788    write(*,*) NF_STRERROR(ierr)
     789    stop ""
     790  endif
     791endif
    711792ierr = NF_DEF_DIM(nout, "Time", NF_UNLIMITED, timedimout)
     793if (ierr.NE.NF_NOERR) then
     794   WRITE(*,*)'initiate: error failed to define dimension <Time>'
     795   write(*,*) NF_STRERROR(ierr)
     796   stop ""
     797endif
    712798
    713799! End netcdf define mode
    714800ierr = NF_ENDDEF(nout)
     801if (ierr.NE.NF_NOERR) then
     802   WRITE(*,*)'initiate: error failed to switch out of define mode'
     803   write(*,*) NF_STRERROR(ierr)
     804   stop ""
     805endif
    715806
    716807!==============================================================================
    … …  
    718809!==============================================================================
    719810
    720 call def_var(nout,"Time","Time","years since 0000-00-0 00:00:00",1,&
     811call def_var(nout,"Time","Ls (Solar Longitude)","degrees",1,&
    721812             (/timedimout/),timevarout,ierr)
     813! Switch to netcdf define mode
     814ierr = NF_REDEF (nout)
     815ierr = NF_ENDDEF(nout)
    722816
    723817!==============================================================================
    … …  
    729823
    730824ierr = NF_PUT_VAR_REAL (nout,nvarid,lat)
     825if (ierr.NE.NF_NOERR) then
     826   WRITE(*,*)'initiate: error failed writing variable <latitude>'
     827   write(*,*) NF_STRERROR(ierr)
     828   stop ""
     829endif
    731830
    732831!==============================================================================
    … …  
    738837
    739838ierr = NF_PUT_VAR_REAL (nout,nvarid,lon)
    740 
    741 !==============================================================================
    742 ! 4. Write "altitude" (data and attributes)
     839if (ierr.NE.NF_NOERR) then
     840   WRITE(*,*)'initiate: error failed writing variable <longitude>'
     841   write(*,*) NF_STRERROR(ierr)
     842   stop ""
     843endif
     844
     845!==============================================================================
     846! 5. Write "altitude" (data and attributes)
    743847!==============================================================================
    744848
    … …  
    748852ierr = NF_DEF_VAR (nout,"altitude",NF_FLOAT,1,altdimout,nvarid)
    749853
    750 ierr = NF_PUT_ATT_TEXT (nout,nvarid,"long_name",8,"altitude")
    751 !********************* temlporaire ****************
    752 !ierr = NF_PUT_ATT_TEXT (nout,nvarid,'units',2,units_alt)
    753 !if((units_alt(1:2).eq.'Pa').or.(units_alt(1:2).eq.'pa')) then
    754 !   ierr = NF_PUT_ATT_TEXT (nout,nvarid,'positive',2,"down")
    755 !else
    756 !   ierr = NF_PUT_ATT_TEXT (nout,nvarid,'positive',2,"up")
    757 !end if
    758 !********************* temlporaire ****************
    759 ierr = NF_PUT_ATT_TEXT (nout,nvarid,'units',2,'Pa')
    760    ierr = NF_PUT_ATT_TEXT (nout,nvarid,'positive',4,'down')
    761 !********************* temlporaire ****************
     854ierr = NF_PUT_ATT_TEXT (nout,nvarid,'long_name',len_trim(adjustl(altlong_name)),adjustl(altlong_name))
     855ierr = NF_PUT_ATT_TEXT (nout,nvarid,'units',len_trim(adjustl(altunits)),adjustl(altunits))
     856ierr = NF_PUT_ATT_TEXT (nout,nvarid,'positive',len_trim(adjustl(altpositive)),adjustl(altpositive))
    762857
    763858! End netcdf define mode
    … …  
    765860
    766861ierr = NF_PUT_VAR_REAL (nout,nvarid,alt)
     862if (ierr.NE.NF_NOERR) then
     863   WRITE(*,*)'initiate: error failed writing variable <altitude>'
     864   write(*,*) NF_STRERROR(ierr)
     865   stop ""
     866endif
     867
     868!==============================================================================
     869! 6. Write "controle" (data and attributes)
     870!==============================================================================
     871
     872if (size(ctl).ne.0) then
     873   ! Switch to netcdf define mode
     874   ierr = NF_REDEF (nout)
     875
     876   ierr = NF_DEF_VAR (nout,"controle",NF_FLOAT,1,ctldimout,nvarid)
     877
     878   ierr = NF_PUT_ATT_TEXT (nout,nvarid,"title",18,"Control parameters")
     879
     880   ! End netcdf define mode
     881   ierr = NF_ENDDEF(nout)
     882
     883   ierr = NF_PUT_VAR_REAL (nout,nvarid,ctl)
     884   if (ierr.NE.NF_NOERR) then
     885      WRITE(*,*)'initiate: error failed writing variable <controle>'
     886      write(*,*) NF_STRERROR(ierr)
     887      stop ""
     888   endif
     889endif
    767890
    768891end Subroutine initiate
    769 !************************************
     892
     893!******************************************************************************
     894subroutine init2(infid,lonlen,latlen,altlen,altdimid, &
     895                 outfid,londimout,latdimout,altdimout)
     896!==============================================================================
     897! Purpose:
     898! Copy aps(), bps() and phisinit() from input file to output file
     899!==============================================================================
     900! Remarks:
     901! The NetCDF files must be open
     902!==============================================================================
     903
     904implicit none
     905
     906include "netcdf.inc" ! NetCDF definitions
     907
     908!==============================================================================
     909! Arguments:
     910!==============================================================================
     911integer, intent(in) :: infid  ! NetCDF output file ID
     912integer, intent(in) :: lonlen ! # of grid points along longitude
     913integer, intent(in) :: latlen ! # of grid points along latitude
     914integer, intent(in) :: altlen ! # of grid points along altitude
     915integer, intent(in) :: altdimid ! ID of altitude dimension
     916integer, intent(in) :: outfid ! NetCDF output file ID
     917integer, intent(in) :: londimout ! longitude dimension ID
     918integer, intent(in) :: latdimout ! latitude dimension ID
     919integer, intent(in) :: altdimout ! altitude dimension ID
     920!==============================================================================
     921! Local variables:
     922!==============================================================================
     923real,dimension(:),allocatable :: aps,bps ! hybrid vertical coordinates
     924real,dimension(:,:),allocatable :: phisinit ! Ground geopotential
     925integer :: apsid,bpsid,phisinitid
     926integer :: ierr
     927integer :: tmpdimid ! temporary dimension ID
     928integer :: tmpvarid ! temporary variable ID
     929integer :: tmplen ! temporary variable length
     930logical :: phis, aps_ok, bps_ok ! are "phisinit" "aps" "bps" available ?
     931
     932
     933!==============================================================================
     934! 1. Read data from input file
     935!==============================================================================
     936
     937! hybrid coordinate aps
     938  allocate(aps(altlen),stat=ierr)
     939  if (ierr.ne.0) then
     940    write(*,*) "init2: failed to allocate aps(altlen)"
     941    stop
     942  endif
     943
     944ierr=NF_INQ_VARID(infid,"aps",tmpvarid)
     945if (ierr.ne.NF_NOERR) then
     946  write(*,*) "oops Failed to get aps ID. OK"
     947  aps_ok=.false.
     948else
     949  ! Check that aps() is the right size (it most likely won't be if
     950  ! zrecast has be used to generate the input file)
     951  ierr=NF_INQ_VARDIMID(infid,tmpvarid,tmpdimid)
     952  ierr=NF_INQ_DIMLEN(infid,tmpdimid,tmplen)
     953  if (tmplen.ne.altlen) then
     954    write(*,*) "tmplen,altlen", tmpdimid, altdimid
     955    write(*,*) "init2: wrong dimension size for aps(), skipping it ..."
     956    aps_ok=.false.
     957  else
     958    ierr=NF_GET_VAR_REAL(infid,tmpvarid,aps)
     959    if (ierr.ne.NF_NOERR) then
     960     stop "init2 error: Failed reading aps"
     961     aps_ok=.false.
     962    endif
     963    aps_ok=.true.
     964  endif
     965endif
     966
     967! hybrid coordinate bps
     968  allocate(bps(altlen),stat=ierr)
     969  if (ierr.ne.0) then
     970    write(*,*) "init2: failed to allocate bps(altlen)"
     971    stop
     972  endif
     973
     974ierr=NF_INQ_VARID(infid,"bps",tmpvarid)
     975if (ierr.ne.NF_NOERR) then
     976  write(*,*) "oops: Failed to get bps ID. OK"
     977  bps_ok=.false.
     978else
     979  ! Check that bps() is the right size
     980  ierr=NF_INQ_VARDIMID(infid,tmpvarid,tmpdimid)
     981  ierr=NF_INQ_DIMLEN(infid,tmpdimid,tmplen)
     982  if (tmplen.ne.altlen) then
     983    write(*,*) "init2: wrong dimension size for bps(), skipping it ..."
     984    bps_ok=.false.
     985  else
     986    ierr=NF_GET_VAR_REAL(infid,tmpvarid,bps)
     987    if (ierr.ne.NF_NOERR) then
     988      stop "init2 Error: Failed reading bps"
     989      bps_ok=.false.
     990    endif
     991    bps_ok=.true.
     992  endif
     993endif
     994
     995! ground geopotential phisinit
     996allocate(phisinit(lonlen,latlen),stat=ierr)
     997if (ierr.ne.0) then
     998  write(*,*) "init2: failed to allocate phisinit(lonlen,latlen)"
     999  stop
     1000endif
     1001ierr=NF_INQ_VARID(infid,"phisinit",tmpvarid)
     1002if (ierr.ne.NF_NOERR) then
     1003  write(*,*) "Failed to get phisinit ID. OK"
     1004  phisinit = 0.
     1005  phis = .false.
     1006else
     1007  ierr=NF_GET_VAR_REAL(infid,tmpvarid,phisinit)
     1008  if (ierr.ne.NF_NOERR) then
     1009    stop "init2 Error: Failed reading phisinit"
     1010  endif
     1011  phis = .true.
     1012endif
     1013
     1014
     1015!==============================================================================
     1016! 2. Write
     1017!==============================================================================
     1018
     1019!==============================================================================
     1020! 2.2. Hybrid coordinates aps() and bps()
     1021!==============================================================================
     1022
     1023IF (aps_ok) then
     1024  ! define aps
     1025  call def_var(outfid,"aps","hybrid pressure at midlayers"," ",1,&
     1026               (/altdimout/),apsid,ierr)
     1027  if (ierr.ne.NF_NOERR) then
     1028    stop "Error: Failed to def_var aps"
     1029  endif
     1030
     1031  ! write aps
     1032  ierr=NF_PUT_VAR_REAL(outfid,apsid,aps)
     1033  if (ierr.ne.NF_NOERR) then
     1034    stop "Error: Failed to write aps"
     1035  endif
     1036ENDIF ! of IF (aps_ok)
     1037
     1038IF (bps_ok) then
     1039  ! define bps
     1040  call def_var(outfid,"bps","hybrid sigma at midlayers"," ",1,&
     1041               (/altdimout/),bpsid,ierr)
     1042  if (ierr.ne.NF_NOERR) then
     1043    stop "Error: Failed to def_var bps"
     1044  endif
     1045
     1046  ! write bps
     1047  ierr=NF_PUT_VAR_REAL(outfid,bpsid,bps)
     1048  if (ierr.ne.NF_NOERR) then
     1049    stop "Error: Failed to write bps"
     1050  endif
     1051ENDIF ! of IF (bps_ok)
     1052
     1053!==============================================================================
     1054! 2.2. phisinit()
     1055!==============================================================================
     1056
     1057IF (phis) THEN
     1058
     1059  !define phisinit
     1060  call def_var(outfid,"phisinit","Ground level geopotential"," ",2,&
     1061              (/londimout,latdimout/),phisinitid,ierr)
     1062  if (ierr.ne.NF_NOERR) then
     1063     stop "Error: Failed to def_var phisinit"
     1064  endif
     1065
     1066  ! write phisinit
     1067  ierr=NF_PUT_VAR_REAL(outfid,phisinitid,phisinit)
     1068  if (ierr.ne.NF_NOERR) then
     1069    stop "Error: Failed to write phisinit"
     1070  endif
     1071
     1072ENDIF ! of IF (phis)
     1073
     1074
     1075! Cleanup
     1076deallocate(aps)
     1077deallocate(bps)
     1078deallocate(phisinit)
     1079
     1080end subroutine init2
     1081
     1082!******************************************************************************
    7701083subroutine def_var(nid,name,title,units,nbdim,dim,nvarid,ierr)
    7711084!==============================================================================
    … …  
    8161129
    8171130end subroutine def_var
    818 !************************************
     1131
     1132!******************************************************************************
     1133subroutine  missing_value(nout,nvarid,missing)
     1134!==============================================================================
     1135! Purpose:
     1136! Write "valid_range" and "missing_value" attributes (of a given
     1137! variable) to a netcdf file
     1138!==============================================================================
     1139! Remarks:
     1140! NetCDF file must be open
     1141! Variable (nvarid) ID must be set
     1142!==============================================================================
     1143
     1144implicit none
     1145
     1146include "netcdf.inc"  ! NetCDF definitions
     1147
     1148!==============================================================================
     1149! Arguments:
     1150!==============================================================================
     1151integer, intent(in) :: nout
     1152! nout: [netcdf] file ID #
     1153integer, intent(in) :: nvarid
     1154! varid: [netcdf] variable ID #
     1155!real, dimension(2), intent(in) :: valid_range
     1156! valid_range(2): [netcdf] "valid_range" attribute (min and max)
     1157real, intent(in) :: missing
     1158! missing: [netcdf] "missing_value" attribute
     1159
     1160!==============================================================================
     1161! Local variables:
     1162!==============================================================================
     1163integer :: ierr
     1164! ierr: [netcdf] subroutine returned error code
     1165!      INTEGER nout,nvarid,ierr
     1166!      REAL missing
     1167!      REAL valid_range(2)
     1168
     1169! Switch to netcdf dataset define mode
     1170ierr = NF_REDEF (nout)
     1171if (ierr.ne.NF_NOERR) then
     1172   print*,'missing_value: '
     1173   print*, NF_STRERROR(ierr)
     1174endif
     1175
     1176
     1177!********* valid range not used in Lslin ****************
     1178! Write valid_range() attribute
     1179!ierr=NF_PUT_ATT_REAL(nout,nvarid,'valid_range',NF_FLOAT,2,valid_range)
     1180!
     1181!if (ierr.ne.NF_NOERR) then
     1182!   print*,'missing_value: valid_range attribution failed'
     1183!   print*, NF_STRERROR(ierr)
     1184!   !write(*,*) 'NF_NOERR', NF_NOERR
     1185!   !CALL abort
     1186!   stop ""
     1187!endif
     1188!*********************************************************
     1189
     1190! Write "missing_value" attribute
     1191ierr= NF_PUT_ATT_REAL(nout,nvarid,'missing_value',NF_FLOAT,1,missing)
     1192if (ierr.NE.NF_NOERR) then
     1193   print*, 'missing_value: missing value attribution failed'
     1194   print*, NF_STRERROR(ierr)
     1195!    WRITE(*,*) 'NF_NOERR', NF_NOERR
     1196!    CALL abort
     1197   stop ""
     1198endif
     1199
     1200! End netcdf dataset define mode
     1201ierr = NF_ENDDEF(nout)
     1202
     1203end subroutine  missing_value
     1204
     1205!******************************************************************************
    8191206subroutine sol2ls(sol,Ls)
    8201207!==============================================================================
    … …  
    9211308integer nd
    9221309!  internal
    923 integer i,j
     1310integer i
    9241311real y_undefined
    9251312                                                                               
    … …  
    9551342end subroutine interpolf
    9561343     
    957 end
    958 
    959 
    960 
    961 !******************************************************************************
    962 !******************************************************************************
    963 subroutine  missing_value(nout,nvarid,valid_range,missing)
    964 !==============================================================================
    965 ! Purpose:
    966 ! Write "valid_range" and "missing_value" attributes (of a given
    967 ! variable) to a netcdf file
    968 !==============================================================================
    969 ! Remarks:
    970 ! NetCDF file must be open
    971 ! Variable (nvarid) ID must be set
    972 !==============================================================================
    973 
    974 implicit none
    975 
    976 include "netcdf.inc"  ! NetCDF definitions
    977 
    978 !==============================================================================
    979 ! Arguments:
    980 !==============================================================================
    981 integer, intent(in) :: nout
    982 ! nout: [netcdf] file ID #
    983 integer, intent(in) :: nvarid
    984 ! varid: [netcdf] variable ID #
    985 real, dimension(2), intent(in) :: valid_range
    986 ! valid_range(2): [netcdf] "valid_range" attribute (min and max)
    987 real, intent(in) :: missing
    988 ! missing: [netcdf] "missing_value" attribute
    989 
    990 !==============================================================================
    991 ! Local variables:
    992 !==============================================================================
    993 integer :: ierr
    994 ! ierr: [netcdf] subroutine returned error code
    995 !      INTEGER nout,nvarid,ierr
    996 !      REAL missing
    997 !      REAL valid_range(2)
    998 
    999 ! Switch to netcdf dataset define mode
    1000 ierr = NF_REDEF (nout)
    1001 if (ierr.ne.NF_NOERR) then
    1002    print*,'missing_value: '
    1003    print*, NF_STRERROR(ierr)
    1004 endif
    1005 
    1006 ! Write valid_range() attribute
    1007 #ifdef NC_DOUBLE
    1008 ierr=NF_PUT_ATT_DOUBLE(nout,nvarid,'valid_range',NF_DOUBLE,2,valid_range)
    1009 #else
    1010 ierr=NF_PUT_ATT_REAL(nout,nvarid,'valid_range',NF_FLOAT,2,valid_range)
    1011 #endif
    1012 
    1013 if (ierr.ne.NF_NOERR) then
    1014    print*,'missing_value: valid_range attribution failed'
    1015    print*, NF_STRERROR(ierr)
    1016    !write(*,*) 'NF_NOERR', NF_NOERR
    1017    !CALL abort
    1018    stop ""
    1019 endif
    1020 
    1021 ! Write "missing_value" attribute
    1022 #ifdef NC_DOUBLE
    1023 ierr= NF_PUT_ATT_DOUBLE(nout,nvarid,'missing_value',NF_DOUBLE,1,missing)
    1024 #else
    1025 ierr= NF_PUT_ATT_REAL(nout,nvarid,'missing_value',NF_FLOAT,1,missing)
    1026 #endif
    1027 
    1028 if (ierr.NE.NF_NOERR) then
    1029    print*, 'missing_value: missing value attribution failed'
    1030    print*, NF_STRERROR(ierr)
    1031 !    WRITE(*,*) 'NF_NOERR', NF_NOERR
    1032 !    CALL abort
    1033    stop ""
    1034 endif
    1035 
    1036 ! End netcdf dataset define mode
    1037 ierr = NF_ENDDEF(nout)
    1038 
    1039 end subroutine  missing_value
    1040 !******************************************************************************
     1344end program lslin

  • trunk/LMDZ.MARS/util/zrecast.F90

    r2401 r2432  
    6161! TN 01/2013 : Adapted for large output files with at least 2 variables > 2 GiB
    6262! TN 10/2013 : Read and write controle field, if available
     63! FF 10/2020 : Force interpolation in log space for variable starting by "rho"
    6364!
    6465implicit none
    … …  
    9495integer,dimension(3) :: surfdatashape ! shape of 3D (surface+time) datasets
    9596
    96 real :: miss_val=-9.99e+33 ! special "missing value" to specify missing data
    97 real,parameter :: miss_val_def=-9.99e+33 ! default value for "missing value"
     97real :: miss_val= 1.E+20 ! special "missing value" to specify missing data
     98real,parameter :: miss_val_def= 1.E+20 ! default value for "missing value"
    9899character (len=64), dimension(:), allocatable :: var
    99100! var(): names of variables that will be processed
    … …  
    18091810  ! interpolate dataset onto new grid
    18101811  ! check if variable is "rho" (to set flag for interpolation below)
    1811   if (var(i).eq.'rho') then
     1812  if (var(i)(1:3).eq.'rho') then
    18121813    j=1
     1814    write(*,*) trim(var(i)), ' is interpolated in log scale'
    18131815  else
    18141816    j=0
    … …  
    18641866    ! interpolate dataset onto new grid
    18651867    ! check if variable is "rho" (to set flag for interpolation below)
    1866     if (var(i).eq.'rho') then
     1868    if (var(i)(1:3).eq.'rho') then
    18671869      j=1
     1870      write(*,*) trim(var(i)), ' is interpolated in log scale'
    18681871    else
    18691872      j=0
    … …  
    19191922  ! interpolate dataset onto new grid
    19201923  ! check if variable is "rho" (to set flag for interpolation below)
    1921   if (var(i).eq.'rho') then
     1924  if (var(i)(1:3).eq.'rho') then
    19221925    j=1
     1926    write(*,*) trim(var(i)), ' is interpolated in log scale'
    19231927  else
    19241928    j=0
Note: See TracChangeset for help on using the changeset viewer.