trunk/phylmd/test_ozonecm.f90

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: disvert, ap, bp, preff, presnivs
  use dynetat0_chosen_m, only: read_serre
  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")

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