trunk/phylmd/test_ozonecm.f90

program test_ozonecm

  ! This is a program in Fortran 95.

  ! This program computes values of ozone abundance from Royer, on a
  ! latitude-pressure grid, and writes the values to a NetCDF file.
  ! The pressure grid is "presnivs" from "disvert".

  use dimensions, only: jjm, llm
  USE dimphy, ONLY : klon
  use disvert_m, only: pa, disvert, ap, bp, preff, presnivs
  use ozonecm_m, only: ozonecm
  use phyetat0_m, only: rlat
  use nr_util, only: arth, assert
  use netcdf95, only: nf95_create, nf95_def_dim, nf95_def_var, nf95_put_att, &
       nf95_enddef, nf95_put_var, nf95_close
  use netcdf, only: nf90_clobber, nf90_float, nf90_global
  use unit_nml_m, only: unit_nml, set_unit_nml

  implicit none

  real p(llm + 1) ! pressure at LMDZ layer interface, in Pa
  real wo(klon, llm, 360) ! column density of ozone in a cell, in kDU
  real tro3(klon, llm, 360) ! mole fraction of ozone
  integer julien, k
  real, parameter:: RG = 9.80665 ! acceleration of gravity, in m s-2
  real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2

  ! For NetCDF:
  integer ncid, dimid_time, dimid_plev, dimid_latitude
  integer varid_time, varid_plev, varid_latitude, varid_tro3

  !---------------------

  call assert(klon == jjm + 1, "test_ozonecm: iim should be 1")

  call set_unit_nml
  open(unit_nml, file="used_namelists.txt", status="replace", action="write")

  pa = 5e4
  call disvert
  p = ap + bp * preff
  rlat = arth(-90., 180. / jjm, jjm + 1)

  do julien = 1, 360
     wo(:, :, julien) = ozonecm(REAL(julien), spread(p, dim=1, ncopies=jjm+1))
  end do

  close(unit_nml)

  forall (k=1: llm)
     tro3(:, k, :) &
          = wo(:, k, :) * 1e3 * dobson_u / (p(k)-p(k+1)) * rg / 48. * 29.
  end forall

  call nf95_create("O3_Royer.nc", nf90_clobber, ncid)
  call nf95_def_dim(ncid, "time", 360, dimid_time)
  call nf95_def_dim(ncid, "plev", llm, dimid_plev)
  call nf95_def_dim(ncid, "latitude", jjm + 1, dimid_latitude)
  call nf95_def_var(ncid, "time", nf90_float, dimid_time, varid_time)
  call nf95_def_var(ncid, "plev", nf90_float, dimid_plev, varid_plev)
  call nf95_def_var(ncid, "latitude", nf90_float, dimid_latitude, &
       varid_latitude)
  call nf95_def_var(ncid, "tro3", nf90_float, &
       (/dimid_latitude, dimid_plev, dimid_time/), varid_tro3)

  call nf95_put_att(ncid, varid_time, "units", "days since 1976-1-1")
  call nf95_put_att(ncid, varid_time, "calendar", "360_day")
  call nf95_put_att(ncid, varid_time, "standard_name", "time")

  call nf95_put_att(ncid, varid_plev, "units", "mbar")
  call nf95_put_att(ncid, varid_plev, "long_name", "pressure")
  call nf95_put_att(ncid, varid_plev, "standard_name", "air_pressure")

  call nf95_put_att(ncid, varid_latitude, "units", "degrees_north")
  call nf95_put_att(ncid, varid_latitude, "standard_name", "latitude")

  call nf95_put_att(ncid, varid_tro3, "long_name", "O3 mole fraction")
  call nf95_put_att(ncid, varid_tro3, "standard_name", &
       "mole_fraction_of_ozone_in_air")

  call nf95_put_att(ncid, nf90_global, "title", &
       "ozone climatology from J.-F. Royer")

  call nf95_enddef(ncid)

  call nf95_put_var(ncid, varid_time, (/(julien - 0.5, julien = 1, 360)/))

  call nf95_put_var(ncid, varid_plev, presnivs(llm:1:-1) / 100.)
  ! sort in ascending order and convert to hPa

  call nf95_put_var(ncid, varid_latitude, rlat)

  call nf95_put_var(ncid, varid_tro3, tro3(:, llm:1:-1, :))
  ! "plev" is in ascending order whereas "p" is descending

  call nf95_close(ncid)

  print *, 'The file "O3_Royer.nc" has been created.'

end program test_ozonecm
1	guez	123	program test_ozonecm
2
3			! This is a program in Fortran 95.
4
5			! This program computes values of ozone abundance from Royer, on a
6			! latitude-pressure grid, and writes the values to a NetCDF file.
7			! The pressure grid is "presnivs" from "disvert".
8
9	guez	265	use dimensions, only: jjm, llm
10	guez	123	USE dimphy, ONLY : klon
11			use disvert_m, only: pa, disvert, ap, bp, preff, presnivs
12			use ozonecm_m, only: ozonecm
13			use phyetat0_m, only: rlat
14			use nr_util, only: arth, assert
15			use netcdf95, only: nf95_create, nf95_def_dim, nf95_def_var, nf95_put_att, &
16			nf95_enddef, nf95_put_var, nf95_close
17			use netcdf, only: nf90_clobber, nf90_float, nf90_global
18	guez	276	use unit_nml_m, only: unit_nml, set_unit_nml
19	guez	123
20			implicit none
21
22			real p(llm + 1) ! pressure at LMDZ layer interface, in Pa
23			real wo(klon, llm, 360) ! column density of ozone in a cell, in kDU
24			real tro3(klon, llm, 360) ! mole fraction of ozone
25			integer julien, k
26			real, parameter:: RG = 9.80665 ! acceleration of gravity, in m s-2
27			real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
28
29			! For NetCDF:
30			integer ncid, dimid_time, dimid_plev, dimid_latitude
31			integer varid_time, varid_plev, varid_latitude, varid_tro3
32
33			!---------------------
34
35			call assert(klon == jjm + 1, "test_ozonecm: iim should be 1")
36
37	guez	276	call set_unit_nml
38	guez	123	open(unit_nml, file="used_namelists.txt", status="replace", action="write")
39
40			pa = 5e4
41			call disvert
42			p = ap + bp * preff
43			rlat = arth(-90., 180. / jjm, jjm + 1)
44
45			do julien = 1, 360
46			wo(:, :, julien) = ozonecm(REAL(julien), spread(p, dim=1, ncopies=jjm+1))
47			end do
48
49			close(unit_nml)
50
51			forall (k=1: llm)
52			tro3(:, k, :) &
53			= wo(:, k, :) * 1e3 * dobson_u / (p(k)-p(k+1)) * rg / 48. * 29.
54			end forall
55
56			call nf95_create("O3_Royer.nc", nf90_clobber, ncid)
57			call nf95_def_dim(ncid, "time", 360, dimid_time)
58			call nf95_def_dim(ncid, "plev", llm, dimid_plev)
59			call nf95_def_dim(ncid, "latitude", jjm + 1, dimid_latitude)
60			call nf95_def_var(ncid, "time", nf90_float, dimid_time, varid_time)
61			call nf95_def_var(ncid, "plev", nf90_float, dimid_plev, varid_plev)
62			call nf95_def_var(ncid, "latitude", nf90_float, dimid_latitude, &
63			varid_latitude)
64			call nf95_def_var(ncid, "tro3", nf90_float, &
65			(/dimid_latitude, dimid_plev, dimid_time/), varid_tro3)
66
67			call nf95_put_att(ncid, varid_time, "units", "days since 1976-1-1")
68			call nf95_put_att(ncid, varid_time, "calendar", "360_day")
69			call nf95_put_att(ncid, varid_time, "standard_name", "time")
70
71			call nf95_put_att(ncid, varid_plev, "units", "mbar")
72			call nf95_put_att(ncid, varid_plev, "long_name", "pressure")
73			call nf95_put_att(ncid, varid_plev, "standard_name", "air_pressure")
74
75			call nf95_put_att(ncid, varid_latitude, "units", "degrees_north")
76			call nf95_put_att(ncid, varid_latitude, "standard_name", "latitude")
77
78			call nf95_put_att(ncid, varid_tro3, "long_name", "O3 mole fraction")
79			call nf95_put_att(ncid, varid_tro3, "standard_name", &
80			"mole_fraction_of_ozone_in_air")
81
82			call nf95_put_att(ncid, nf90_global, "title", &
83			"ozone climatology from J.-F. Royer")
84
85			call nf95_enddef(ncid)
86
87			call nf95_put_var(ncid, varid_time, (/(julien - 0.5, julien = 1, 360)/))
88
89			call nf95_put_var(ncid, varid_plev, presnivs(llm:1:-1) / 100.)
90			! sort in ascending order and convert to hPa
91
92			call nf95_put_var(ncid, varid_latitude, rlat)
93
94			call nf95_put_var(ncid, varid_tro3, tro3(:, llm:1:-1, :))
95			! "plev" is in ascending order whereas "p" is descending
96
97			call nf95_close(ncid)
98
99			print *, 'The file "O3_Royer.nc" has been created.'
100
101			end program test_ozonecm