trunk/phylmd/test_ozonecm.f

program test_ozonecm

  ! This is a program in Fortran 2003.

  ! 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 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 dimensions, only: jjm, llm
  USE dimphy, ONLY : klon
  USE dimsoil, ONLY : nsoilmx
  use disvert_m, only: pa, disvert, ap, bp, preff, presnivs
  USE indicesol, ONLY : nbsrf
  use ozonecm_m, only: ozonecm
  use phyetat0_m, only: rlat, phyetat0
  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

  REAL pctsrf(klon, nbsrf)
  REAL ftsol(klon, nbsrf)
  REAL ftsoil(klon, nsoilmx, nbsrf)
  REAL qsurf(klon, nbsrf)
  REAL qsol(klon) ! column-density of water in soil, in kg m-2
  REAL snow(klon, nbsrf)
  REAL albe(klon, nbsrf)
  REAL rain_fall(klon)
  REAL snow_fall(klon)
  real solsw(klon)
  REAL sollw(klon)
  real fder(klon)
  REAL radsol(klon)
  REAL frugs(klon, nbsrf)
  REAL agesno(klon, nbsrf)
  REAL zmea(klon)
  REAL zstd(klon)
  REAL zsig(klon)
  REAL zgam(klon)
  REAL zthe(klon)
  REAL zpic(klon)
  REAL zval(klon)
  REAL t_ancien(klon, llm), q_ancien(klon, llm)
  LOGICAL ancien_ok
  real rnebcon(klon, llm), ratqs(klon, llm)
  REAL clwcon(klon, llm), run_off_lic_0(klon)
  real sig1(klon, llm) ! section adiabatic updraft

  real w01(klon, llm) 
  ! vertical velocity within adiabatic updraft

  integer ncid_startphy

  ! 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
  call phyetat0(pctsrf, ftsol, ftsoil, qsurf, qsol, snow, albe, rain_fall, &
       snow_fall, solsw, sollw, fder, radsol, frugs, agesno, zmea, zstd, zsig, &
       zgam, zthe, zpic, zval, t_ancien, q_ancien, ancien_ok, rnebcon, ratqs, &
       clwcon, run_off_lic_0, sig1, w01, ncid_startphy)

  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	guez	278	! This is a program in Fortran 2003.
4	guez	123
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	278	use nr_util, only: arth, assert
10			use netcdf95, only: nf95_create, nf95_def_dim, nf95_def_var, nf95_put_att, &
11			nf95_enddef, nf95_put_var, nf95_close
12			use netcdf, only: nf90_clobber, nf90_float, nf90_global
13
14	guez	265	use dimensions, only: jjm, llm
15	guez	123	USE dimphy, ONLY : klon
16	guez	278	USE dimsoil, ONLY : nsoilmx
17	guez	123	use disvert_m, only: pa, disvert, ap, bp, preff, presnivs
18	guez	278	USE indicesol, ONLY : nbsrf
19	guez	123	use ozonecm_m, only: ozonecm
20	guez	278	use phyetat0_m, only: rlat, phyetat0
21	guez	276	use unit_nml_m, only: unit_nml, set_unit_nml
22	guez	123
23			implicit none
24
25			real p(llm + 1) ! pressure at LMDZ layer interface, in Pa
26			real wo(klon, llm, 360) ! column density of ozone in a cell, in kDU
27			real tro3(klon, llm, 360) ! mole fraction of ozone
28			integer julien, k
29			real, parameter:: RG = 9.80665 ! acceleration of gravity, in m s-2
30			real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
31
32	guez	278	REAL pctsrf(klon, nbsrf)
33			REAL ftsol(klon, nbsrf)
34			REAL ftsoil(klon, nsoilmx, nbsrf)
35			REAL qsurf(klon, nbsrf)
36			REAL qsol(klon) ! column-density of water in soil, in kg m-2
37			REAL snow(klon, nbsrf)
38			REAL albe(klon, nbsrf)
39			REAL rain_fall(klon)
40			REAL snow_fall(klon)
41			real solsw(klon)
42			REAL sollw(klon)
43			real fder(klon)
44			REAL radsol(klon)
45			REAL frugs(klon, nbsrf)
46			REAL agesno(klon, nbsrf)
47			REAL zmea(klon)
48			REAL zstd(klon)
49			REAL zsig(klon)
50			REAL zgam(klon)
51			REAL zthe(klon)
52			REAL zpic(klon)
53			REAL zval(klon)
54			REAL t_ancien(klon, llm), q_ancien(klon, llm)
55			LOGICAL ancien_ok
56			real rnebcon(klon, llm), ratqs(klon, llm)
57			REAL clwcon(klon, llm), run_off_lic_0(klon)
58			real sig1(klon, llm) ! section adiabatic updraft
59
60			real w01(klon, llm)
61			! vertical velocity within adiabatic updraft
62
63			integer ncid_startphy
64
65	guez	123	! For NetCDF:
66			integer ncid, dimid_time, dimid_plev, dimid_latitude
67			integer varid_time, varid_plev, varid_latitude, varid_tro3
68
69			!---------------------
70
71			call assert(klon == jjm + 1, "test_ozonecm: iim should be 1")
72
73	guez	276	call set_unit_nml
74	guez	123	open(unit_nml, file="used_namelists.txt", status="replace", action="write")
75
76			pa = 5e4
77			call disvert
78			p = ap + bp * preff
79	guez	304	call phyetat0(pctsrf, ftsol, ftsoil, qsurf, qsol, snow, albe, rain_fall, &
80			snow_fall, solsw, sollw, fder, radsol, frugs, agesno, zmea, zstd, zsig, &
81			zgam, zthe, zpic, zval, t_ancien, q_ancien, ancien_ok, rnebcon, ratqs, &
82			clwcon, run_off_lic_0, sig1, w01, ncid_startphy)
83	guez	123
84			do julien = 1, 360
85			wo(:, :, julien) = ozonecm(REAL(julien), spread(p, dim=1, ncopies=jjm+1))
86			end do
87
88			close(unit_nml)
89
90			forall (k=1: llm)
91			tro3(:, k, :) &
92			= wo(:, k, :) * 1e3 * dobson_u / (p(k)-p(k+1)) * rg / 48. * 29.
93			end forall
94
95			call nf95_create("O3_Royer.nc", nf90_clobber, ncid)
96			call nf95_def_dim(ncid, "time", 360, dimid_time)
97			call nf95_def_dim(ncid, "plev", llm, dimid_plev)
98			call nf95_def_dim(ncid, "latitude", jjm + 1, dimid_latitude)
99			call nf95_def_var(ncid, "time", nf90_float, dimid_time, varid_time)
100			call nf95_def_var(ncid, "plev", nf90_float, dimid_plev, varid_plev)
101			call nf95_def_var(ncid, "latitude", nf90_float, dimid_latitude, &
102			varid_latitude)
103			call nf95_def_var(ncid, "tro3", nf90_float, &
104			(/dimid_latitude, dimid_plev, dimid_time/), varid_tro3)
105
106			call nf95_put_att(ncid, varid_time, "units", "days since 1976-1-1")
107			call nf95_put_att(ncid, varid_time, "calendar", "360_day")
108			call nf95_put_att(ncid, varid_time, "standard_name", "time")
109
110			call nf95_put_att(ncid, varid_plev, "units", "mbar")
111			call nf95_put_att(ncid, varid_plev, "long_name", "pressure")
112			call nf95_put_att(ncid, varid_plev, "standard_name", "air_pressure")
113
114			call nf95_put_att(ncid, varid_latitude, "units", "degrees_north")
115			call nf95_put_att(ncid, varid_latitude, "standard_name", "latitude")
116
117			call nf95_put_att(ncid, varid_tro3, "long_name", "O3 mole fraction")
118			call nf95_put_att(ncid, varid_tro3, "standard_name", &
119			"mole_fraction_of_ozone_in_air")
120
121			call nf95_put_att(ncid, nf90_global, "title", &
122			"ozone climatology from J.-F. Royer")
123
124			call nf95_enddef(ncid)
125
126			call nf95_put_var(ncid, varid_time, (/(julien - 0.5, julien = 1, 360)/))
127
128			call nf95_put_var(ncid, varid_plev, presnivs(llm:1:-1) / 100.)
129			! sort in ascending order and convert to hPa
130
131			call nf95_put_var(ncid, varid_latitude, rlat)
132
133			call nf95_put_var(ncid, varid_tro3, tro3(:, llm:1:-1, :))
134			! "plev" is in ascending order whereas "p" is descending
135
136			call nf95_close(ncid)
137
138			print *, 'The file "O3_Royer.nc" has been created.'
139
140			end program test_ozonecm