From 8624a75e20c0cc6dff42b4eb97182a138f9aa1e7 Mon Sep 17 00:00:00 2001
From: Don Morton <Don.Morton@borealscicomp.com>
Date: Thu, 24 Aug 2017 19:07:08 +0000
Subject: [PATCH] Enhancements to FPv9.3.2
Documented in Ticket #182 (as well as CTBTO ticket 357)
---
flexpart_code/fpmetbinary_mod.F90 | 581 ++++++++++++++-
flexpart_code/grib2nc4/Makefile | 41 +-
flexpart_code/grib2nc4/README | 44 +-
flexpart_code/grib2nc4/fp2nc4io_mod.F90 | 43 ++
flexpart_code/grib2nc4/grib2nc4.F90 | 74 +-
flexpart_code/grib2nc4/nc4cmp.F90 | 244 ++++++
flexpart_code/grib2nc4/processmetfields.F90 | 17 +-
.../grib2nc4/verttransform_grib2nc4_ecmwf.f90 | 701 ++++++++++++++++++
.../grib2nc4/verttransform_grib2nc4_gfs.f90 | 680 +++++++++++++++++
9 files changed, 2378 insertions(+), 47 deletions(-)
create mode 100644 flexpart_code/grib2nc4/nc4cmp.F90
create mode 100644 flexpart_code/grib2nc4/verttransform_grib2nc4_ecmwf.f90
create mode 100644 flexpart_code/grib2nc4/verttransform_grib2nc4_gfs.f90
diff --git a/flexpart_code/fpmetbinary_mod.F90 b/flexpart_code/fpmetbinary_mod.F90
index a5b6eec1..e0a65048 100644
--- a/flexpart_code/fpmetbinary_mod.F90
+++ b/flexpart_code/fpmetbinary_mod.F90
@@ -237,38 +237,82 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_def_dim(ncid, 'nxmax', nxmax, nxmax_dimid)
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","maximum dimension of wind fields in x")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
ncret = nf90_def_dim(ncid, 'nymax', nymax, nymax_dimid)
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","maximum dimension of wind fields in y")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
ncret = nf90_def_dim(ncid, 'nzmax', nzmax, nzmax_dimid)
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","maximum dimension of wind fields in z")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
ncret = nf90_def_dim(ncid, 'nuvzmax', nuvzmax, nuvzmax_dimid)
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description"," maximum dimension of (u,v) wind field in z")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
ncret = nf90_def_dim(ncid, 'nwzmax', nwzmax, nwzmax_dimid)
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description"," maximum dimension of w component of wind field in z")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
ncret = nf90_def_dim(ncid, 'maxspec', maxspec, maxspec_dimid)
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description"," maximum number of species")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
ncret = nf90_def_dim(ncid, 'numclass', numclass, numclass_dimid)
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description"," maximum number of ageclasses")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
! There are a handful of variables indexed from 0 to n, rather than 0 to n-1,
! so these dimensions handle that. What a pain.
ncret = nf90_def_dim(ncid, 'zero_to_nzmax', nzmax+1, zero_to_nzmax_dimid)
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","variable to change indexing of variables with nzmax dimension")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
ncret = nf90_def_dim(ncid, 'zero_to_maxnests', maxnests+1, zero_to_maxnests_dimid)
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","variable to change indexing of variables with maxnests dimension")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
! This is for a couple of small arrays that store polar stereographic stuff
ncret = nf90_def_dim(ncid, 'polemap_dim', 9, polemap_dimid)
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","variable to store polar stereographic indexing")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
! These two values come from conv_mod
ncret = nf90_def_dim(ncid, 'nconvlevmax_dim', nconvlevmax, nconvlevmax_dimid)
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","maximum number of levels for convection")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
ncret = nf90_def_dim(ncid, 'na_dim', na, na_dimid)
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","maximum number of levels for convection +1")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
- ! Scalar values
+ ! Scalar values
dim1dids = (/preproc_fmt_str_dimid/)
ncret = nf90_def_var(ncid, 'preproc_fmt_str', NF90_CHAR, &
@@ -287,83 +331,130 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
call handle_nf90_err(ncret)
ncret = nf90_put_var(ncid, ncvarid, nx)
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","actual dimensions of wind fields in x")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
ncret = nf90_def_var(ncid, 'ny', NF90_INT, ncvarid)
call handle_nf90_err(ncret)
ncret = nf90_put_var(ncid, ncvarid, ny)
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","actual dimensions of wind fields in y")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
ncret = nf90_def_var(ncid, 'nxmin1', NF90_INT, ncvarid)
call handle_nf90_err(ncret)
ncret = nf90_put_var(ncid, ncvarid, nxmin1)
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","actual dimensions of wind fields in x -1 (nx-1)")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
ncret = nf90_def_var(ncid, 'nymin1', NF90_INT, ncvarid)
call handle_nf90_err(ncret)
ncret = nf90_put_var(ncid, ncvarid, nymin1)
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","actual dimensions of wind fields in y -1 (ny-1)")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
ncret = nf90_def_var(ncid, 'nxfield', NF90_INT, ncvarid)
call handle_nf90_err(ncret)
ncret = nf90_put_var(ncid, ncvarid, nxfield)
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","actual dimensions of wind fields in x for limited area fields")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
ncret = nf90_def_var(ncid, 'nuvz', NF90_INT, ncvarid)
call handle_nf90_err(ncret)
ncret = nf90_put_var(ncid, ncvarid, nuvz)
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","actual dimensions of (u,v) wind fields in z direction")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
ncret = nf90_def_var(ncid, 'nwz', NF90_INT, ncvarid)
call handle_nf90_err(ncret)
ncret = nf90_put_var(ncid, ncvarid, nwz)
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","actual dimensions of w wind fields in z direction")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
ncret = nf90_def_var(ncid, 'nz', NF90_INT, ncvarid)
call handle_nf90_err(ncret)
ncret = nf90_put_var(ncid, ncvarid, nz)
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","number of vertical levels in the transformed coordinates")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
ncret = nf90_def_var(ncid, 'nmixz', NF90_INT, ncvarid)
call handle_nf90_err(ncret)
ncret = nf90_put_var(ncid, ncvarid, nmixz)
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","number of levels up to maximum PBL height (3500 m)")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
ncret = nf90_def_var(ncid, 'nlev_ec', NF90_INT, ncvarid)
call handle_nf90_err(ncret)
ncret = nf90_put_var(ncid, ncvarid, nlev_ec)
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","number of vertical levels ecmwf model")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
ncret = nf90_def_var(ncid, 'dx', NF90_FLOAT, ncvarid)
call handle_nf90_err(ncret)
ncret = nf90_put_var(ncid, ncvarid, dx)
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","grid distance in x direction")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
ncret = nf90_def_var(ncid, 'dy', NF90_FLOAT, ncvarid)
call handle_nf90_err(ncret)
ncret = nf90_put_var(ncid, ncvarid, dy)
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","grid distance in y direction")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
ncret = nf90_def_var(ncid, 'xlon0', NF90_FLOAT, ncvarid)
call handle_nf90_err(ncret)
ncret = nf90_put_var(ncid, ncvarid, xlon0)
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","geographical longitude of the lower left corner of the wind fields")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
ncret = nf90_def_var(ncid, 'ylat0', NF90_FLOAT, ncvarid)
call handle_nf90_err(ncret)
ncret = nf90_put_var(ncid, ncvarid, ylat0)
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","geographical latitude of the lower left corner of the wind fields")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
ncret = nf90_def_var(ncid, 'dxconst', NF90_FLOAT, ncvarid)
call handle_nf90_err(ncret)
ncret = nf90_put_var(ncid, ncvarid, dxconst)
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","auxiliary variable needed for utransform")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
ncret = nf90_def_var(ncid, 'dyconst', NF90_FLOAT, ncvarid)
call handle_nf90_err(ncret)
ncret = nf90_put_var(ncid, ncvarid, dyconst)
call handle_nf90_err(ncret)
-
-
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","auxiliary variable needed for vtransform")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
! Fixed fields, static in time
dim2dids = (/nxmax_dimid, nymax_dimid/)
@@ -379,6 +470,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& oro(0:nxmax-1, 0:nymax-1))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","orography of the ECMWF model")
+ ncret = nf90_put_att(ncid, ncvarid, "units","m")
+
ncret = nf90_def_var(ncid, 'excessoro', NF90_FLOAT, &
& dim2dids, ncvarid)
@@ -391,6 +486,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& excessoro(0:nxmax-1, 0:nymax-1))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","excess orography mother domain")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
ncret = nf90_def_var(ncid, 'lsm', NF90_FLOAT, &
& dim2dids, ncvarid)
@@ -403,6 +502,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& lsm(0:nxmax-1, 0:nymax-1))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","land sea mask")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
dim3dids = (/nxmax_dimid, nymax_dimid, numclass_dimid/)
! numclass comes from par_mod - number of land use classes
@@ -417,6 +520,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& xlanduse(0:nxmax-1, 0:nymax-1, 1:numclass))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","area fraction in percent")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
dim1dids = (/nzmax_dimid/)
ncret = nf90_def_var(ncid, 'height', NF90_FLOAT, &
@@ -430,7 +537,9 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& height(1:nzmax))
call handle_nf90_err(ncret)
-
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","heights of all model levels")
+ ncret = nf90_put_att(ncid, ncvarid, "units","m")
@@ -448,6 +557,9 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& uu(0:nxmax-1, 0:nymax-1, 1:nzmax, cm_index))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","U component of wind in the [horizontal] direction")
+ ncret = nf90_put_att(ncid, ncvarid, "units","m s**-1")
ncret = nf90_def_var(ncid, 'vv', NF90_FLOAT, &
& dim3dids, ncvarid)
@@ -460,6 +572,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& vv(0:nxmax-1, 0:nymax-1, 1:nzmax, cm_index))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","V component of wind in the Y[horizontal] direction")
+ ncret = nf90_put_att(ncid, ncvarid, "units","m s**-1")
+
ncret = nf90_def_var(ncid, 'uupol', NF90_FLOAT, &
& dim3dids, ncvarid)
@@ -472,6 +588,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& uupol(0:nxmax-1, 0:nymax-1, 1:nzmax, cm_index))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","horizontal component in polar stereographic projection of wind")
+ ncret = nf90_put_att(ncid, ncvarid, "units","m s**-1")
+
ncret = nf90_def_var(ncid, 'vvpol', NF90_FLOAT, &
& dim3dids, ncvarid)
@@ -484,6 +604,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& vvpol(0:nxmax-1, 0:nymax-1, 1:nzmax, cm_index))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","horizontal component in polar stereographic projection of wind")
+ ncret = nf90_put_att(ncid, ncvarid, "units","m s**-1")
+
ncret = nf90_def_var(ncid, 'ww', NF90_FLOAT, &
& dim3dids, ncvarid)
@@ -496,6 +620,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& ww(0:nxmax-1, 0:nymax-1, 1:nzmax, cm_index))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","wind component in the Z[vertical] direction")
+ ncret = nf90_put_att(ncid, ncvarid, "units","m s**-1")
+
ncret = nf90_def_var(ncid, 'tt', NF90_FLOAT, &
& dim3dids, ncvarid)
@@ -508,6 +636,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& tt(0:nxmax-1, 0:nymax-1, 1:nzmax, cm_index))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","temperature")
+ ncret = nf90_put_att(ncid, ncvarid, "units","K")
+
ncret = nf90_def_var(ncid, 'qv', NF90_FLOAT, &
& dim3dids, ncvarid)
@@ -520,6 +652,9 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& qv(0:nxmax-1, 0:nymax-1, 1:nzmax, cm_index))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","specific humidity data")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
ncret = nf90_def_var(ncid, 'pv', NF90_FLOAT, &
& dim3dids, ncvarid)
@@ -532,6 +667,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& pv(0:nxmax-1, 0:nymax-1, 1:nzmax, cm_index))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","potential vorticity")
+ ncret = nf90_put_att(ncid, ncvarid, "units","K*m**2 kg**-1 s**-1")
+
ncret = nf90_def_var(ncid, 'rho', NF90_FLOAT, &
& dim3dids, ncvarid)
@@ -544,6 +683,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& rho(0:nxmax-1, 0:nymax-1, 1:nzmax, cm_index))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","air density")
+ ncret = nf90_put_att(ncid, ncvarid, "units","kg m**-3")
+
ncret = nf90_def_var(ncid, 'drhodz', NF90_FLOAT, &
& dim3dids, ncvarid)
@@ -556,6 +699,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& drhodz(0:nxmax-1, 0:nymax-1, 1:nzmax, cm_index))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","vertical air density gradient")
+ ncret = nf90_put_att(ncid, ncvarid, "units","kg m**-2")
+
ncret = nf90_def_var(ncid, 'clouds', NF90_BYTE, &
& dim3dids, ncvarid)
@@ -568,7 +715,12 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& clouds(0:nxmax-1, 0:nymax-1, 1:nzmax, cm_index))
call handle_nf90_err(ncret)
-
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","vcloud mask, &
+ no cloud no precip = 0, cloud no precip = 1, &
+ rainout conv/lsp dominated 2/3, &
+ washout conv/lsp dominated = 4/5 ")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
! Note the change in z dimension for the following
@@ -585,6 +737,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& tth(0:nxmax-1, 0:nymax-1, 1:nuvzmax, cm_index))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","temperature in the original input model level")
+ ncret = nf90_put_att(ncid, ncvarid, "units","K")
+
ncret = nf90_def_var(ncid, 'qvh', NF90_FLOAT, &
& dim3dids, ncvarid)
@@ -597,6 +753,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& qvh(0:nxmax-1, 0:nymax-1, 1:nuvzmax, cm_index))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","specific humidity in the original input model level")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
ncret = nf90_def_var(ncid, 'pplev', NF90_FLOAT, &
& dim3dids, ncvarid)
@@ -609,6 +769,9 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& pplev(0:nxmax-1, 0:nymax-1, 1:nuvzmax, cm_index))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","pressure in pressure levels for GFS version")
+ ncret = nf90_put_att(ncid, ncvarid, "units","Pa")
dim2dids = (/nxmax_dimid, nymax_dimid/)
@@ -623,6 +786,9 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& cloudsh(0:nxmax-1, 0:nymax-1, cm_index))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","height of the cloud layer")
+ ncret = nf90_put_att(ncid, ncvarid, "units","m")
@@ -648,6 +814,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& ps(0:nxmax-1, 0:nymax-1, 1, cm_index))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","surface pressure")
+ ncret = nf90_put_att(ncid, ncvarid, "units","Pa")
+
ncret = nf90_def_var(ncid, 'sd', NF90_FLOAT, &
& dim2dids, ncvarid)
@@ -660,6 +830,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& sd(0:nxmax-1, 0:nymax-1, 1, cm_index))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","snow depth")
+ ncret = nf90_put_att(ncid, ncvarid, "units","m")
+
ncret = nf90_def_var(ncid, 'msl', NF90_FLOAT, &
& dim2dids, ncvarid)
@@ -672,6 +846,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& msl(0:nxmax-1, 0:nymax-1, 1, cm_index))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","mean sea level pressure")
+ ncret = nf90_put_att(ncid, ncvarid, "units","Pa")
+
ncret = nf90_def_var(ncid, 'tcc', NF90_FLOAT, &
& dim2dids, ncvarid)
@@ -684,6 +862,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& tcc(0:nxmax-1, 0:nymax-1, 1, cm_index))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","total cloud cover")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
ncret = nf90_def_var(ncid, 'u10', NF90_FLOAT, &
& dim2dids, ncvarid)
@@ -696,6 +878,9 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& u10(0:nxmax-1, 0:nymax-1, 1, cm_index))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","10 m u component of wind velocity")
+ ncret = nf90_put_att(ncid, ncvarid, "units","m s**-1")
ncret = nf90_def_var(ncid, 'v10', NF90_FLOAT, &
& dim2dids, ncvarid)
@@ -708,6 +893,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& v10(0:nxmax-1, 0:nymax-1, 1, cm_index))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","10 m v component of wind velocity")
+ ncret = nf90_put_att(ncid, ncvarid, "units","m s**-1")
+
ncret = nf90_def_var(ncid, 'tt2', NF90_FLOAT, &
& dim2dids, ncvarid)
@@ -720,6 +909,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& tt2(0:nxmax-1, 0:nymax-1, 1, cm_index))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","2 m temperature")
+ ncret = nf90_put_att(ncid, ncvarid, "units","K")
+
ncret = nf90_def_var(ncid, 'td2', NF90_FLOAT, &
& dim2dids, ncvarid)
@@ -732,6 +925,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& td2(0:nxmax-1, 0:nymax-1, 1, cm_index))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","2 m dew point temperature")
+ ncret = nf90_put_att(ncid, ncvarid, "units","K")
+
ncret = nf90_def_var(ncid, 'lsprec', NF90_FLOAT, &
& dim2dids, ncvarid)
@@ -744,6 +941,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& lsprec(0:nxmax-1, 0:nymax-1, 1, cm_index))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","large scale total precipitation")
+ ncret = nf90_put_att(ncid, ncvarid, "units","mm h**-1")
+
ncret = nf90_def_var(ncid, 'convprec', NF90_FLOAT, &
& dim2dids, ncvarid)
@@ -756,6 +957,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& convprec(0:nxmax-1, 0:nymax-1, 1, cm_index))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","convective precipitation")
+ ncret = nf90_put_att(ncid, ncvarid, "units","mm h**-1")
+
ncret = nf90_def_var(ncid, 'sshf', NF90_FLOAT, &
& dim2dids, ncvarid)
@@ -768,6 +973,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& sshf(0:nxmax-1, 0:nymax-1, 1, cm_index))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","surface sensible heat flux")
+ ncret = nf90_put_att(ncid, ncvarid, "units","J m**-2")
+
ncret = nf90_def_var(ncid, 'ssr', NF90_FLOAT, &
& dim2dids, ncvarid)
@@ -780,6 +989,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& ssr(0:nxmax-1, 0:nymax-1, 1, cm_index))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","surface solar radiation")
+ ncret = nf90_put_att(ncid, ncvarid, "units","J m**-2")
+
ncret = nf90_def_var(ncid, 'surfstr', NF90_FLOAT, &
& dim2dids, ncvarid)
@@ -792,6 +1005,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& surfstr(0:nxmax-1, 0:nymax-1, 1, cm_index))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","surface stress")
+ ncret = nf90_put_att(ncid, ncvarid, "units","N m**-2 s")
+
ncret = nf90_def_var(ncid, 'ustar', NF90_FLOAT, &
& dim2dids, ncvarid)
@@ -804,6 +1021,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& ustar(0:nxmax-1, 0:nymax-1, 1, cm_index))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","friction velocity")
+ ncret = nf90_put_att(ncid, ncvarid, "units","m s**-1")
+
ncret = nf90_def_var(ncid, 'wstar', NF90_FLOAT, &
& dim2dids, ncvarid)
@@ -816,6 +1037,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& wstar(0:nxmax-1, 0:nymax-1, 1, cm_index))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","convective velocity scale")
+ ncret = nf90_put_att(ncid, ncvarid, "units","m s**-1")
+
ncret = nf90_def_var(ncid, 'hmix', NF90_FLOAT, &
& dim2dids, ncvarid)
@@ -828,6 +1053,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& hmix(0:nxmax-1, 0:nymax-1, 1, cm_index))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","mixing height")
+ ncret = nf90_put_att(ncid, ncvarid, "units","m")
+
ncret = nf90_def_var(ncid, 'tropopause', NF90_FLOAT, &
& dim2dids, ncvarid)
@@ -840,6 +1069,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& tropopause(0:nxmax-1, 0:nymax-1, 1, cm_index))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","altitude of thermal tropopause")
+ ncret = nf90_put_att(ncid, ncvarid, "units","m")
+
ncret = nf90_def_var(ncid, 'oli', NF90_FLOAT, &
& dim2dids, ncvarid)
@@ -852,6 +1085,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& oli(0:nxmax-1, 0:nymax-1, 1, cm_index))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","inverse Obukhov length (1/L)")
+ ncret = nf90_put_att(ncid, ncvarid, "units","m")
+
ncret = nf90_def_var(ncid, 'diffk', NF90_FLOAT, &
& dim2dids, ncvarid)
@@ -864,7 +1101,9 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& diffk(0:nxmax-1, 0:nymax-1, 1, cm_index))
call handle_nf90_err(ncret)
-
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","diffusion coefficient at reference height")
+ ncret = nf90_put_att(ncid, ncvarid, "units","m**2 s**-1")
! !PRINT *, 'SUM(ps(0:nxmax-1, 0:nymax-1, 1, cm_index)): ', &
@@ -887,7 +1126,9 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& vdep(0:nxmax-1, 0:nymax-1, 1:maxspec, cm_index))
call handle_nf90_err(ncret)
-
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","deposition velocity")
+ ncret = nf90_put_att(ncid, ncvarid, "units","m s**-2")
! 1d fields
@@ -901,6 +1142,9 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
& deflate_level=DEF_LEVEL)
ncret = nf90_put_var(ncid, ncvarid, &
& z0(1:numclass))
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","surface roughness length")
+ ncret = nf90_put_att(ncid, ncvarid, "units","m")
dim1dids = (/nwzmax_dimid/)
@@ -916,6 +1160,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& akm(1:nwzmax))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","ECMWF vertical discretization parameter")
+ ncret = nf90_put_att(ncid, ncvarid, "units","Pa")
+
ncret = nf90_def_var(ncid, 'bkm', NF90_FLOAT, &
& dim1dids, ncvarid)
@@ -928,6 +1176,9 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& bkm(1:nwzmax))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","ECMWF vertical discretization parameter")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
dim1dids = (/nuvzmax_dimid/)
@@ -943,6 +1194,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& akz(1:nuvzmax))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","ECMWF vertical discretization parameter at center of layer")
+ ncret = nf90_put_att(ncid, ncvarid, "units","Pa")
+
ncret = nf90_def_var(ncid, 'bkz', NF90_FLOAT, &
& dim1dids, ncvarid)
@@ -955,6 +1210,9 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& bkz(1:nuvzmax))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","ECMWF vertical discretization parameter at center of layer")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
dim1dids = (/nzmax_dimid/)
@@ -970,6 +1228,9 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& aknew(1:nzmax))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","model discretization coefficient at the interpolated levels")
+ ncret = nf90_put_att(ncid, ncvarid, "units","Pa")
ncret = nf90_def_var(ncid, 'bknew', NF90_FLOAT, &
& dim1dids, ncvarid)
@@ -982,6 +1243,9 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& bknew(1:nzmax))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","model discretization coefficient at the interpolated levels")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
! PRINT *, 'SUM(bknew(1:nzmax)): ', &
@@ -1016,6 +1280,9 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& nxn(1:maxnests))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","actual dimensions of wind fields in x for the nexted domain")
+
ncret = nf90_def_var(ncid, 'nyn', NF90_INT, &
& dim1dids, ncvarid)
@@ -1028,6 +1295,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& nyn(1:maxnests))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","actual dimensions of wind fields in y for the nested domain")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
ncret = nf90_def_var(ncid, 'dxn', NF90_FLOAT, &
& dim1dids, ncvarid)
@@ -1040,6 +1311,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& dxn(1:maxnests))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","grid distance in x direction for the nested domain")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
ncret = nf90_def_var(ncid, 'dyn', NF90_FLOAT, &
& dim1dids, ncvarid)
@@ -1052,6 +1327,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& dyn(1:maxnests))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","grid distance in y direction for the nested domain")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
ncret = nf90_def_var(ncid, 'xlon0n', NF90_FLOAT, &
& dim1dids, ncvarid)
@@ -1064,6 +1343,11 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& xlon0n(1:maxnests))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","geographical longitude of the &
+ lower left corner of the nested wind fields")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
ncret = nf90_def_var(ncid, 'ylat0n', NF90_FLOAT, &
& dim1dids, ncvarid)
@@ -1076,8 +1360,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& ylat0n(1:maxnests))
call handle_nf90_err(ncret)
-
-
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","geographical latitude of &
+ the lower left corner of the nested wind fields")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
! Nested fields, static over time
@@ -1094,6 +1380,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& oron(0:nxmaxn-1, 0:nymaxn-1, 1:maxnests))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","orography of the ECMWF model for the nested domain")
+ ncret = nf90_put_att(ncid, ncvarid, "units","m")
+
ncret = nf90_def_var(ncid, 'excessoron', NF90_FLOAT, &
& dim3dids, ncvarid)
@@ -1106,6 +1396,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& excessoron(0:nxmaxn-1, 0:nymaxn-1, 1:maxnests))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","excess orography nested domain")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
ncret = nf90_def_var(ncid, 'lsmn', NF90_FLOAT, &
& dim3dids, ncvarid)
@@ -1118,6 +1412,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& lsmn(0:nxmaxn-1, 0:nymaxn-1, 1:maxnests))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","land sea mask for the nested domian")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
dim4dids = (/nxmaxn_dimid, nymaxn_dimid, numclass_dimid, maxnests_dimid/)
@@ -1132,6 +1430,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& xlandusen(0:nxmaxn-1, 0:nymaxn-1, 1:numclass, 1:maxnests))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","area fraction in percent for the nested domain")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
! PRINT *, 'SUM(oron): ', SUM(oron)
@@ -1151,6 +1453,11 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& uun(0:nxmaxn-1, 0:nymaxn-1, 1:nzmax, cm_index, 1:maxnests))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","U component of wind in the [horizontal] &
+ direction for the nested domain")
+ ncret = nf90_put_att(ncid, ncvarid, "units","m s**-1")
+
ncret = nf90_def_var(ncid, 'vvn', NF90_FLOAT, &
& dim4dids, ncvarid)
@@ -1163,6 +1470,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& vvn(0:nxmaxn-1, 0:nymaxn-1, 1:nzmax, cm_index, 1:maxnests))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","V component of wind in the Y[horizontal] &
+ direction for the nested domain")
+ ncret = nf90_put_att(ncid, ncvarid, "units","m s**-1")
ncret = nf90_def_var(ncid, 'wwn', NF90_FLOAT, &
& dim4dids, ncvarid)
@@ -1175,6 +1486,11 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& wwn(0:nxmaxn-1, 0:nymaxn-1, 1:nzmax, cm_index, 1:maxnests))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","wind componend in the Z[vertical] &
+ direction for the nested domain")
+ ncret = nf90_put_att(ncid, ncvarid, "units","m s**-1")
+
ncret = nf90_def_var(ncid, 'ttn', NF90_FLOAT, &
& dim4dids, ncvarid)
@@ -1187,6 +1503,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& ttn(0:nxmaxn-1, 0:nymaxn-1, 1:nzmax, cm_index, 1:maxnests))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","temperature for the nested domain")
+ ncret = nf90_put_att(ncid, ncvarid, "units","K")
+
ncret = nf90_def_var(ncid, 'qvn', NF90_FLOAT, &
& dim4dids, ncvarid)
@@ -1199,6 +1519,9 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& qvn(0:nxmaxn-1, 0:nymaxn-1, 1:nzmax, cm_index, 1:maxnests))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","specific humidity data for the nested domain")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
ncret = nf90_def_var(ncid, 'pvn', NF90_FLOAT, &
& dim4dids, ncvarid)
@@ -1211,6 +1534,9 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& pvn(0:nxmaxn-1, 0:nymaxn-1, 1:nzmax, cm_index, 1:maxnests))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","potential vorticity for the nested domain")
+ ncret = nf90_put_att(ncid, ncvarid, "units","K*m**2 kg**-1 s**-1")
ncret = nf90_def_var(ncid, 'rhon', NF90_FLOAT, &
& dim4dids, ncvarid)
@@ -1223,6 +1549,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& rhon(0:nxmaxn-1, 0:nymaxn-1, 1:nzmax, cm_index, 1:maxnests))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","air density for the nested domain")
+ ncret = nf90_put_att(ncid, ncvarid, "units","kg m**-3")
+
ncret = nf90_def_var(ncid, 'drhodzn', NF90_FLOAT, &
& dim4dids, ncvarid)
@@ -1235,6 +1565,9 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& drhodzn(0:nxmaxn-1, 0:nymaxn-1, 1:nzmax, cm_index, 1:maxnests))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","vertical air density gradient for the nested domain")
+ ncret = nf90_put_att(ncid, ncvarid, "units","kg m**-2")
! Note the new dimensions
@@ -1251,6 +1584,11 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& tthn(0:nxmaxn-1, 0:nymaxn-1, 1:nuvzmax, cm_index, 1:maxnests))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","temperature in the original &
+ input model level for the nested domain")
+ ncret = nf90_put_att(ncid, ncvarid, "units","K")
+
ncret = nf90_def_var(ncid, 'qvhn', NF90_FLOAT, &
& dim4dids, ncvarid)
@@ -1263,6 +1601,11 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& qvhn(0:nxmaxn-1, 0:nymaxn-1, 1:nuvzmax, cm_index, 1:maxnests))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","specific humidity in the original &
+ input model level for the nested domain")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
! Note the new dimensions
dim4dids = (/nxmaxn_dimid, nymaxn_dimid, zero_to_nzmax_dimid, maxnests_dimid/)
@@ -1278,6 +1621,13 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& cloudsn(0:nxmaxn-1, 0:nymaxn-1, 0:nzmax, cm_index, 1:maxnests))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","vcloud mask, &
+ no cloud no precip = 0, cloud no precip = 1, &
+ rainout conv/lsp dominated 2/3, &
+ washout conv/lsp dominated = 4/5 for the nested domain")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
! Note the new dimensions
dim3dids = (/nxmaxn_dimid, nymaxn_dimid, maxnests_dimid/)
@@ -1293,6 +1643,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& cloudsnh(0:nxmaxn-1, 0:nymaxn-1, cm_index, 1:maxnests))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","cloud mask at the horizontal levels")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
! PRINT *, 'SUM(uun): ', SUM(uun(:,:,:,cm_index,:))
@@ -1315,6 +1669,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& psn(0:nxmaxn-1, 0:nymaxn-1, 1, cm_index, 1:maxnests))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","surface pressure of the nested domian")
+ ncret = nf90_put_att(ncid, ncvarid, "units","Pa")
+
ncret = nf90_def_var(ncid, 'sdn', NF90_FLOAT, &
& dim3dids, ncvarid)
@@ -1327,6 +1685,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& sdn(0:nxmaxn-1, 0:nymaxn-1, 1, cm_index, 1:maxnests))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","snow depth for the nested domain")
+ ncret = nf90_put_att(ncid, ncvarid, "units","m")
+
ncret = nf90_def_var(ncid, 'msln', NF90_FLOAT, &
& dim3dids, ncvarid)
@@ -1339,6 +1701,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& msln(0:nxmaxn-1, 0:nymaxn-1, 1, cm_index, 1:maxnests))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","mean sea level pressure for the nested domain")
+ ncret = nf90_put_att(ncid, ncvarid, "units","Pa")
+
ncret = nf90_def_var(ncid, 'tccn', NF90_FLOAT, &
& dim3dids, ncvarid)
@@ -1351,6 +1717,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& tccn(0:nxmaxn-1, 0:nymaxn-1, 1, cm_index, 1:maxnests))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","total cloud cover for the nested domain")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
ncret = nf90_def_var(ncid, 'u10n', NF90_FLOAT, &
& dim3dids, ncvarid)
@@ -1363,6 +1733,9 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& u10n(0:nxmaxn-1, 0:nymaxn-1, 1, cm_index, 1:maxnests))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","10 m u component of wind velocity for the nested domain")
+ ncret = nf90_put_att(ncid, ncvarid, "units","m s**-1")
ncret = nf90_def_var(ncid, 'v10n', NF90_FLOAT, &
& dim3dids, ncvarid)
@@ -1375,6 +1748,9 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& v10n(0:nxmaxn-1, 0:nymaxn-1, 1, cm_index, 1:maxnests))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","10 m v component of wind velocity for the nested domain")
+ ncret = nf90_put_att(ncid, ncvarid, "units","m s**-1")
ncret = nf90_def_var(ncid, 'tt2n', NF90_FLOAT, &
& dim3dids, ncvarid)
@@ -1387,6 +1763,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& tt2n(0:nxmaxn-1, 0:nymaxn-1, 1, cm_index, 1:maxnests))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","2 m temperature for the nested domain")
+ ncret = nf90_put_att(ncid, ncvarid, "units","K")
+
ncret = nf90_def_var(ncid, 'td2n', NF90_FLOAT, &
& dim3dids, ncvarid)
@@ -1399,6 +1779,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& td2n(0:nxmaxn-1, 0:nymaxn-1, 1, cm_index, 1:maxnests))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","2 m dew point temperature for the nested domain")
+ ncret = nf90_put_att(ncid, ncvarid, "units","K")
+
ncret = nf90_def_var(ncid, 'lsprecn', NF90_FLOAT, &
& dim3dids, ncvarid)
@@ -1411,6 +1795,9 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& lsprecn(0:nxmaxn-1, 0:nymaxn-1, 1, cm_index, 1:maxnests))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","large scale total precipitation for the nested domain")
+ ncret = nf90_put_att(ncid, ncvarid, "units","mm h**-1")
ncret = nf90_def_var(ncid, 'convprecn', NF90_FLOAT, &
& dim3dids, ncvarid)
@@ -1423,6 +1810,9 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& convprecn(0:nxmaxn-1, 0:nymaxn-1, 1, cm_index, 1:maxnests))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","convective precipitation for the nested domain")
+ ncret = nf90_put_att(ncid, ncvarid, "units","mm h**-1")
ncret = nf90_def_var(ncid, 'sshfn', NF90_FLOAT, &
& dim3dids, ncvarid)
@@ -1435,6 +1825,9 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& sshfn(0:nxmaxn-1, 0:nymaxn-1, 1, cm_index, 1:maxnests))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","surface sensible heat flux for the nested domain")
+ ncret = nf90_put_att(ncid, ncvarid, "units","J m**-2")
ncret = nf90_def_var(ncid, 'ssrn', NF90_FLOAT, &
& dim3dids, ncvarid)
@@ -1447,6 +1840,9 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& ssrn(0:nxmaxn-1, 0:nymaxn-1, 1, cm_index, 1:maxnests))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","surface solar radiation for the nested domain")
+ ncret = nf90_put_att(ncid, ncvarid, "units","J m**-2")
ncret = nf90_def_var(ncid, 'surfstrn', NF90_FLOAT, &
& dim3dids, ncvarid)
@@ -1459,6 +1855,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& surfstrn(0:nxmaxn-1, 0:nymaxn-1, 1, cm_index, 1:maxnests))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","surface stress for the nested domain")
+ ncret = nf90_put_att(ncid, ncvarid, "units","N m**-2 s")
+
ncret = nf90_def_var(ncid, 'ustarn', NF90_FLOAT, &
& dim3dids, ncvarid)
@@ -1471,6 +1871,9 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& ustarn(0:nxmaxn-1, 0:nymaxn-1, 1, cm_index, 1:maxnests))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","friction velocity for the nested domain")
+ ncret = nf90_put_att(ncid, ncvarid, "units","m s**-1")
ncret = nf90_def_var(ncid, 'wstarn', NF90_FLOAT, &
& dim3dids, ncvarid)
@@ -1483,6 +1886,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& wstarn(0:nxmaxn-1, 0:nymaxn-1, 1, cm_index, 1:maxnests))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","convective velocity scale for the nested domain")
+ ncret = nf90_put_att(ncid, ncvarid, "units","m s**-1")
+
ncret = nf90_def_var(ncid, 'hmixn', NF90_FLOAT, &
& dim3dids, ncvarid)
@@ -1495,6 +1902,9 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& hmixn(0:nxmaxn-1, 0:nymaxn-1, 1, cm_index, 1:maxnests))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","mixing height for the nested domain")
+ ncret = nf90_put_att(ncid, ncvarid, "units","m")
ncret = nf90_def_var(ncid, 'tropopausen', NF90_FLOAT, &
& dim3dids, ncvarid)
@@ -1507,6 +1917,9 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& tropopausen(0:nxmaxn-1, 0:nymaxn-1, 1, cm_index, 1:maxnests))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","altitude of thermal tropopause for the nested domain")
+ ncret = nf90_put_att(ncid, ncvarid, "units","m")
ncret = nf90_def_var(ncid, 'olin', NF90_FLOAT, &
& dim3dids, ncvarid)
@@ -1519,6 +1932,9 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& olin(0:nxmaxn-1, 0:nymaxn-1, 1, cm_index, 1:maxnests))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","inverse Obukhov length (1/L) for the nested domain")
+ ncret = nf90_put_att(ncid, ncvarid, "units","m")
ncret = nf90_def_var(ncid, 'diffkn', NF90_FLOAT, &
& dim3dids, ncvarid)
@@ -1531,6 +1947,12 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& diffkn(0:nxmaxn-1, 0:nymaxn-1, 1, cm_index, 1:maxnests))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","diffusion coefficient at reference height &
+ for the nested domain")
+ ncret = nf90_put_att(ncid, ncvarid, "units","m**2 s**-1")
+
+
dim4dids = (/nxmaxn_dimid, nymaxn_dimid, maxspec_dimid, maxnests_dimid/)
@@ -1547,7 +1969,9 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& vdepn(0:nxmaxn-1, 0:nymaxn-1, 1:maxspec, cm_index, 1:maxnests))
call handle_nf90_err(ncret)
-
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","deposition velocity for the nested domain")
+ ncret = nf90_put_att(ncid, ncvarid, "units","m s**-2")
! PRINT *, 'SUM(psn): ', SUM(psn(:,:,:,cm_index,:))
! PRINT *, 'SUM(surfstrn): ', SUM(surfstrn(:,:,:,cm_index,:))
@@ -1569,6 +1993,11 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& xresoln(0:maxnests))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","factor by which the resolutions in the nests &
+ is enhanced compared to mother grid")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
ncret = nf90_def_var(ncid, 'yresoln', NF90_FLOAT, &
& dim1dids, ncvarid)
@@ -1581,6 +2010,11 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& yresoln(0:maxnests))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","factor by which the resolutions in the nests &
+ is enhanced compared to mother grid")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
dim1dids = (/maxnests_dimid/)
@@ -1595,6 +2029,11 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& xln(1:maxnests))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","lower left corner x point of nested grids &
+ in grid coordinates of mother grid")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
ncret = nf90_def_var(ncid, 'yln', NF90_FLOAT, &
& dim1dids, ncvarid)
@@ -1607,6 +2046,11 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& yln(1:maxnests))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","lower left corner y point of nested grids &
+ in grid coordinates of mother grid")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
ncret = nf90_def_var(ncid, 'xrn', NF90_FLOAT, &
& dim1dids, ncvarid)
@@ -1619,6 +2063,11 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& xrn(1:maxnests))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","upper right corner x point of nested grids &
+ in grid coordinates of mother grid")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
ncret = nf90_def_var(ncid, 'yrn', NF90_FLOAT, &
& dim1dids, ncvarid)
@@ -1631,6 +2080,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& yrn(1:maxnests))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","upper right corner y point of nested grids &
+ in grid coordinates of mother grid")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
! PRINT *, 'SUM(yresoln): ', SUM(yresoln)
! PRINT *, 'SUM(xrn): ', SUM(xrn)
@@ -1648,6 +2101,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
& deflate_level=DEF_LEVEL)
ncret = nf90_put_var(ncid, ncvarid, &
& southpolemap(:))
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","auxiliary variable to define stereographic projections")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
ncret = nf90_def_var(ncid, 'northpolemap', NF90_FLOAT, &
& dim1dids, ncvarid)
@@ -1657,7 +2114,9 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
& deflate_level=DEF_LEVEL)
ncret = nf90_put_var(ncid, ncvarid, &
& northpolemap(:))
-
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","auxiliary variable to define stereographic projections")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
! xglobal, sglobal, nglobal are LOGICAL vars, and need to be converted
@@ -1666,26 +2125,45 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
call handle_nf90_err(ncret)
ncret = nf90_put_var(ncid, ncvarid, logical2integer(xglobal))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","logical variable T for global fields, F for limited area fields")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
ncret = nf90_def_var(ncid, 'sglobal', NF90_INT, ncvarid)
call handle_nf90_err(ncret)
ncret = nf90_put_var(ncid, ncvarid, logical2integer(sglobal))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","logical variable T if domain extends towards south pole")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
ncret = nf90_def_var(ncid, 'nglobal', NF90_INT, ncvarid)
call handle_nf90_err(ncret)
ncret = nf90_put_var(ncid, ncvarid, logical2integer(nglobal))
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","logical variable T if domain extends towards north pole")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
ncret = nf90_def_var(ncid, 'switchnorthg', NF90_FLOAT, ncvarid)
call handle_nf90_err(ncret)
ncret = nf90_put_var(ncid, ncvarid, switchnorthg)
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","degrees for use polar stereographic grid north of switchnorth")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
ncret = nf90_def_var(ncid, 'switchsouthg', NF90_FLOAT, ncvarid)
call handle_nf90_err(ncret)
ncret = nf90_put_var(ncid, ncvarid, switchsouthg)
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","degrees for use polar stereographic grid south of switchshouth")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
@@ -1706,6 +2184,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
& deflate_level=DEF_LEVEL)
ncret = nf90_put_var(ncid, ncvarid, &
& pconv(:))
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","pressure at the input model levels for the convection calculation")
+ ncret = nf90_put_att(ncid, ncvarid, "units","Pa")
+
ncret = nf90_def_var(ncid, 'dpr', NF90_FLOAT, &
& dim1dids, ncvarid)
@@ -1715,6 +2197,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
& deflate_level=DEF_LEVEL)
ncret = nf90_put_var(ncid, ncvarid, &
& dpr(:))
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","pressure difference around tconv(k)")
+ ncret = nf90_put_att(ncid, ncvarid, "units","Pa")
+
ncret = nf90_def_var(ncid, 'pconv_hpa', NF90_FLOAT, &
& dim1dids, ncvarid)
@@ -1724,6 +2210,11 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
& deflate_level=DEF_LEVEL)
ncret = nf90_put_var(ncid, ncvarid, &
& pconv_hpa(:))
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","pressure at the input model levels &
+ in hPa for the convection calculation")
+ ncret = nf90_put_att(ncid, ncvarid, "units","hPa")
+
ncret = nf90_def_var(ncid, 'ft', NF90_FLOAT, &
& dim1dids, ncvarid)
@@ -1733,6 +2224,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
& deflate_level=DEF_LEVEL)
ncret = nf90_put_var(ncid, ncvarid, &
& ft(:))
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","auxiliary variable for the convection scheme")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
ncret = nf90_def_var(ncid, 'fq', NF90_FLOAT, &
& dim1dids, ncvarid)
@@ -1742,6 +2237,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
& deflate_level=DEF_LEVEL)
ncret = nf90_put_var(ncid, ncvarid, &
& fq(:))
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","auxiliary variable for the convection scheme")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
ncret = nf90_def_var(ncid, 'sub', NF90_FLOAT, &
& dim1dids, ncvarid)
@@ -1752,6 +2251,11 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
ncret = nf90_put_var(ncid, ncvarid, &
& sub(:))
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","auxiliary variable for the convection scheme")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
+
dim1dids = (/na_dimid/)
ncret = nf90_def_var(ncid, 'phconv', NF90_FLOAT, &
@@ -1762,6 +2266,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
& deflate_level=DEF_LEVEL)
ncret = nf90_put_var(ncid, ncvarid, &
& phconv(:))
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","pressure between levels k-1 and k")
+ ncret = nf90_put_att(ncid, ncvarid, "units","Pa")
+
ncret = nf90_def_var(ncid, 'phconv_hpa', NF90_FLOAT, &
& dim1dids, ncvarid)
@@ -1771,6 +2279,9 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
& deflate_level=DEF_LEVEL)
ncret = nf90_put_var(ncid, ncvarid, &
& phconv_hpa(:))
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","pressure between levels k-1 and k in hPa")
+ ncret = nf90_put_att(ncid, ncvarid, "units","hPa")
ncret = nf90_def_var(ncid, 'tconv', NF90_FLOAT, &
& dim1dids, ncvarid)
@@ -1780,6 +2291,11 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
& deflate_level=DEF_LEVEL)
ncret = nf90_put_var(ncid, ncvarid, &
& tconv(:))
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","temperature at the input model levels &
+ for the convection calculation")
+ ncret = nf90_put_att(ncid, ncvarid, "units","K")
+
ncret = nf90_def_var(ncid, 'qconv', NF90_FLOAT, &
& dim1dids, ncvarid)
@@ -1789,6 +2305,11 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
& deflate_level=DEF_LEVEL)
ncret = nf90_put_var(ncid, ncvarid, &
& qconv(:))
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","specific humidity at the input model &
+ levels for the convection calculation")
+ ncret = nf90_put_att(ncid, ncvarid, "units","none")
+
ncret = nf90_def_var(ncid, 'qsconv', NF90_FLOAT, &
& dim1dids, ncvarid)
@@ -1798,6 +2319,11 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
& deflate_level=DEF_LEVEL)
ncret = nf90_put_var(ncid, ncvarid, &
& qsconv(:))
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","Saturation water vapor specific humidity &
+ for the convection calculation")
+ ncret = nf90_put_att(ncid, ncvarid, "units","kg/kg")
+
! New dimensions
dim2dids = (/nconvlevmax_dimid, nconvlevmax_dimid/)
@@ -1810,6 +2336,10 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
& deflate_level=DEF_LEVEL)
ncret = nf90_put_var(ncid, ncvarid, &
& fmass(:,:))
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","mass in the grid column for the convection calculation")
+ ncret = nf90_put_att(ncid, ncvarid, "units","kg")
+
ncret = nf90_def_var(ncid, 'fmassfrac', NF90_FLOAT, &
& dim2dids, ncvarid)
@@ -1819,6 +2349,11 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
& deflate_level=DEF_LEVEL)
ncret = nf90_put_var(ncid, ncvarid, &
& fmassfrac(:,:))
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","mass fraction due to the convection &
+ displacement from level k to level k")
+ ncret = nf90_put_att(ncid, ncvarid, "units","kg")
+
! New dimensions
@@ -1832,6 +2367,9 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
& deflate_level=DEF_LEVEL)
ncret = nf90_put_var(ncid, ncvarid, &
& cbaseflux(0:nxmax-1,0:nymax-1))
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","cloudbase massflux due to the convection")
+ ncret = nf90_put_att(ncid, ncvarid, "units","(kg/m**2)/s")
! New dimensions
dim3dids = (/nxmaxn_dimid, nymaxn_dimid, maxnests_dimid/)
@@ -1844,6 +2382,9 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
& deflate_level=DEF_LEVEL)
ncret = nf90_put_var(ncid, ncvarid, &
& cbasefluxn(0:nxmaxn-1,0:nymaxn-1,1:maxnests))
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","cloudbase massflux due to the convection for the nested domain")
+ ncret = nf90_put_att(ncid, ncvarid, "units","(kg/m**2)/s")
! Scalars
@@ -1851,26 +2392,44 @@ PRINT *, 'OPENED NC4 FILE FOR READING...'
call handle_nf90_err(ncret)
ncret = nf90_put_var(ncid, ncvarid, psconv)
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","surface pressure for the convection calculation")
+ ncret = nf90_put_att(ncid, ncvarid, "units","Pa")
+
ncret = nf90_def_var(ncid, 'tt2conv', NF90_FLOAT, ncvarid)
call handle_nf90_err(ncret)
ncret = nf90_put_var(ncid, ncvarid, tt2conv)
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","2 m temperature for the convection calculation")
+ ncret = nf90_put_att(ncid, ncvarid, "units","K")
+
ncret = nf90_def_var(ncid, 'td2conv', NF90_FLOAT, ncvarid)
call handle_nf90_err(ncret)
ncret = nf90_put_var(ncid, ncvarid, td2conv)
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","2 m dew point temperature for the convection calculation")
+ ncret = nf90_put_att(ncid, ncvarid, "units","K")
+
ncret = nf90_def_var(ncid, 'nconvlev', NF90_INT, ncvarid)
call handle_nf90_err(ncret)
ncret = nf90_put_var(ncid, ncvarid, nconvlev)
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","number of levels for convection")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
ncret = nf90_def_var(ncid, 'nconvtop', NF90_INT, ncvarid)
call handle_nf90_err(ncret)
ncret = nf90_put_var(ncid, ncvarid, nconvtop)
call handle_nf90_err(ncret)
+! attributes
+ ncret = nf90_put_att(ncid, ncvarid, "description","upper levels for convection")
+ ncret = nf90_put_att(ncid, ncvarid, "units"," ")
! PRINT *, 'SUM(pconv): ', SUM(pconv)
! PRINT *, 'SUM(qconv): ', SUM(qconv)
diff --git a/flexpart_code/grib2nc4/Makefile b/flexpart_code/grib2nc4/Makefile
index e938350e..b4e672eb 100644
--- a/flexpart_code/grib2nc4/Makefile
+++ b/flexpart_code/grib2nc4/Makefile
@@ -10,37 +10,41 @@ GRIBAPI = /opt/grib-api
HDF5 = /opt/hdf5-1.8.16
NETCDFF = /opt/netcdf-fortran-4.4.3
NETCDF = /opt/netcdf-c-4.4.0
-
+#GRIBAPI = /usr/local/grib-api
+#HDF5 = /usr/local/hdf5-1.8.16
+#NETCDFF = /usr/local/netcdf-fortran-4.4.3
+#NETCDF = /usr/local/netcdf-c-4.4.0
BINARY = grib2nc4
-OBJS = processmetfields.o
-FPMODOBJS = par_mod.o com_mod.o class_vtable_mod.o cmapf_mod.o conv_mod.o
-FLXPRTOBJS = detectformat.o grib2check.o shift_field_0.o gridcheck.o \
- readwind.o readwind_nests.o calcpar.o calcpar_nests.o \
- shift_field.o pbl_profile.o scalev.o obukhov.o \
- richardson.o ew.o getvdep.o calcpv.o obukhov_gfs.o \
- richardson_gfs.o getvdep_nests.o calcpv_nests.o psim.o psih.o \
- qvsat.o caldate.o getrb.o raerod.o getrc.o partdep.o \
- verttransform.o verttransform_nests.o readwind_gfs.o \
- calcpar_gfs.o verttransform_gfs.o gridcheck_gfs.o
+CMP_BINARY = nc4cmp
+OBJS = processmetfields.o verttransform_grib2nc4_ecmwf.o verttransform_grib2nc4_gfs.o
+FPMODOBJS = ${FLEXPART_SRC}/par_mod.o ${FLEXPART_SRC}/com_mod.o ${FLEXPART_SRC}/class_vtable_mod.o ${FLEXPART_SRC}/cmapf_mod.o ${FLEXPART_SRC}/conv_mod.o
+FLXPRTOBJS = ${FLEXPART_SRC}/detectformat.o ${FLEXPART_SRC}/grib2check.o ${FLEXPART_SRC}/shift_field_0.o ${FLEXPART_SRC}/gridcheck.o \
+ ${FLEXPART_SRC}/readwind.o ${FLEXPART_SRC}/readwind_nests.o ${FLEXPART_SRC}/calcpar.o ${FLEXPART_SRC}/calcpar_nests.o \
+ ${FLEXPART_SRC}/shift_field.o ${FLEXPART_SRC}/pbl_profile.o ${FLEXPART_SRC}/scalev.o ${FLEXPART_SRC}/obukhov.o \
+ ${FLEXPART_SRC}/richardson.o ${FLEXPART_SRC}/ew.o ${FLEXPART_SRC}/getvdep.o ${FLEXPART_SRC}/calcpv.o ${FLEXPART_SRC}/obukhov_gfs.o \
+ ${FLEXPART_SRC}/richardson_gfs.o ${FLEXPART_SRC}/getvdep_nests.o ${FLEXPART_SRC}/calcpv_nests.o ${FLEXPART_SRC}/psim.o ${FLEXPART_SRC}/psih.o \
+ ${FLEXPART_SRC}/qvsat.o ${FLEXPART_SRC}/caldate.o ${FLEXPART_SRC}/getrb.o ${FLEXPART_SRC}/raerod.o ${FLEXPART_SRC}/getrc.o ${FLEXPART_SRC}/partdep.o \
+ ${FLEXPART_SRC}/verttransform.o ${FLEXPART_SRC}/verttransform_nests.o ${FLEXPART_SRC}/readwind_gfs.o \
+ ${FLEXPART_SRC}/calcpar_gfs.o ${FLEXPART_SRC}/verttransform_gfs.o ${FLEXPART_SRC}/gridcheck_gfs.o
VPATH = ${FLEXPART_SRC}
FFLAGS = -mcmodel=medium
INCLUDES_NETCDF = -I${NETCDFF}/include
-INCLUDES = -I${GRIBAPI}/include ${INCLUDES_NETCDF}
+INCLUDES = -I${GRIBAPI}/include ${INCLUDES_NETCDF} -I${FLEXPART_SRC}
### NetCDF link flags - use the first one for dynamic libs, the second
### one for static libs
-#LDFLAGS_NETCDF = -L${NETCDFF}/lib -lnetcdff -L${NETCDF}/lib -lnetcdf
-LDFLAGS_NETCDF=-static -L${NETCDFF}/lib -lnetcdff -L${NETCDF}/lib -lnetcdf -lnetcdf -L${HDF5}/lib -lhdf5_fortran -lhdf5_hl -lhdf5hl_fortran -lhdf5 -ldl -lz
+LDFLAGS_NETCDF = -L${NETCDFF}/lib -lnetcdff -L${NETCDF}/lib -lnetcdf
+#LDFLAGS_NETCDF=-static -L${NETCDFF}/lib -lnetcdff -L${NETCDF}/lib -lnetcdf -lnetcdf -L${HDF5}/lib -lhdf5_fortran -lhdf5_hl -lhdf5hl_fortran -lhdf5 -ldl -lz
-LDFLAGS = -L${GRIBAPI}/lib -lgrib_api_f90 -lgrib_api ${LDFLAGS_NETCDF} -ljasper
+LDFLAGS = -L${GRIBAPI}/lib -lgrib_api_f90 -lgrib_api ${LDFLAGS_NETCDF} -ljasper -L${FLEXPART_SRC}
#------------ Creating the binary ------------------
@@ -51,6 +55,13 @@ ${BINARY} : ${BINARY}.o fp2nc4io_mod.o ${FLXPRTOBJS} ${FPMODOBJS} ${OBJS}
${BINARY}.o : ${BINARY}.F90 fp2nc4io_mod.mod ${FPMODOBJS} Makefile
${FC} -c ${BINARY}.F90 ${FFLAGS} ${INCLUDES}
+#----------- NC4 compare ------------------------
+${CMP_BINARY} : ${CMP_BINARY}.o
+ ${FC} -o ${CMP_BINARY} ${CMP_BINARY}.o ${LDFLAGS}
+
+${CMP_BINARY}.o : ${CMP_BINARY}.F90 Makefile
+ ${FC} -c ${CMP_BINARY}.F90 ${FFLAGS} ${INCLUDES}
+
fp2nc4io_mod.mod : ${FPMODOBJS}
${FC} -c fp2nc4io_mod.F90 ${FFLAGS} ${INCLUDES}
diff --git a/flexpart_code/grib2nc4/README b/flexpart_code/grib2nc4/README
index 231f49cc..aa8b1ff3 100644
--- a/flexpart_code/grib2nc4/README
+++ b/flexpart_code/grib2nc4/README
@@ -6,7 +6,7 @@ Contact: Don Morton, don.morton@borealscicomp.com
Delia Arnold, deliona.arnold@gmail.com
M. Harustak
-Last Update: 29 May 2016
+Last Update: 05 July 2017
======================================================================
1. Introduction
@@ -19,6 +19,10 @@ all the operations that FLEXPART normally does in order to prepare the met
data for computations. It then takes the data which has been stored in
FLEXPART global arrays and selectively writes it to a NetCDF4 file.
+Optional parameter allows user to request arrays calculated for specified
+location instead of default algorithm selecting lowest left corner of meteorological
+domain with a surface pressure above 100000Pa.
+
This distribution is integrated into the FLEXPART source code, in a grib2nc4
subdirectory, but it doesn't necessarily have to reside there. The
FLEXPART_SRC variable in the Makefile allows for an alternate FP source code
@@ -41,6 +45,14 @@ write it to NC4 file, then read it back and compare with the preprocessed
data still in the FLEXPART arrays. It does not test whether preprocessing
is good, but only whether data was stored correctly in the NC4 files.
+To be able to quickly compare variables in NetCDF files, simple utility
+nc4cmp was introduced. It compares specified variable in two NetCDF files
+(with dimension up to 3) and it prints average and maximal difference between
+the values. It allows for specification of optional tolerance (in %).
+Difference is calculated as abs((value1-value2)/value1)*100.
+It is calculated for all values of variable multidimensional array.
+
+
======================================================================
2. Quick-Start
======================================================================
@@ -64,6 +76,8 @@ distribution, then there should be no need to change FLEXPART_SRC.
Simply type "make" (noting that you may see some warnings related to
variables for nests) to produce grib2nc4.
+NetCDF compare utility nc4cmp can be built using command "make nc4cmp".
+
Before trying to run, you might want to try "make test" to make sure that
the simple test works as expected.
@@ -80,6 +94,11 @@ WARNING - there is currently nothing to protect you if you get these out
of order. It is possible for you to overwrite your met file if you use it
as the second argument.
+Usage:
+
+ ./grib2nc4 <inpath> <outpath> [lat=value lon=value] [optional varnames]
+
+
Sample usage:
./grib2nc4 GD15051200 def9_nc1p0.nc4
@@ -92,6 +111,29 @@ them to the command. For example:
./grib2nc4 GD15051200 def9_nc1p0.nc4 q w
+Calculation at user specified location can be requested by adding coord parameter
+to the command line:
+
+ ./grib2nc4 GD15051200 def9_nc1p0.nc4 lon=20 lat=45
+
+To compare two NetCDF files, nc4cmp can be used.
+
+Usage:
+
+ ./nc4cmp <file1> <file2> <variable> [optional tolerance in %]
+
+Sample usage:
+
+ ./nc4cmp def9_nc1p0.nc4 def9_nc1p0_reference.nc4 U
+
+Sample usage with tolerance 1%:
+
+ ./nc4cmp def9_nc1p0.nc4 def9_nc1p0_reference.nc4 U 1.0
+
+This will extract variable U from both NetCDF files, print and compare dimensions,
+compare values and print average and maximal difference.
+
+
======================================================================
3. Installation
======================================================================
diff --git a/flexpart_code/grib2nc4/fp2nc4io_mod.F90 b/flexpart_code/grib2nc4/fp2nc4io_mod.F90
index ebda676c..6db0ce94 100644
--- a/flexpart_code/grib2nc4/fp2nc4io_mod.F90
+++ b/flexpart_code/grib2nc4/fp2nc4io_mod.F90
@@ -156,6 +156,9 @@ CONTAINS
! Write the height field - variable 'height' is defined in com_mod
ncfunc_retval = nf90_def_var(ncid, 'height', NF90_DOUBLE, &
& z_dimid, varid)
+! attributes
+ ncfunc_retval = nf90_put_att(ncid, varid, "description","height of the FLEXPART model levels")
+ ncfunc_retval = nf90_put_att(ncid, varid, "units","m a.g.l")
ncfunc_retval = nf90_def_var_deflate(ncid, varid, &
& shuffle=0, &
@@ -168,15 +171,31 @@ CONTAINS
! dx, dy, xlon0, xlat0 are all defined in com_mod
ncfunc_retval = nf90_def_var(ncid, 'dx', NF90_DOUBLE, varid)
ncfunc_retval = nf90_put_var(ncid, varid, dx)
+! attributes
+ ncfunc_retval = nf90_put_att(ncid, varid, "description","grid distance in x direction")
+ ncfunc_retval = nf90_put_att(ncid, varid, "units","degrees")
+
ncfunc_retval = nf90_def_var(ncid, 'dy', NF90_DOUBLE, varid)
ncfunc_retval = nf90_put_var(ncid, varid, dy)
+! attributes
+ ncfunc_retval = nf90_put_att(ncid, varid, "description","grid distance in y direction")
+ ncfunc_retval = nf90_put_att(ncid, varid, "units","degrees")
+
ncfunc_retval = nf90_def_var(ncid, 'xlon0', NF90_DOUBLE, varid)
ncfunc_retval = nf90_put_var(ncid, varid, xlon0)
+! attributes
+ ncfunc_retval = nf90_put_att(ncid, varid, "description","longitude of the lowest left corner")
+ ncfunc_retval = nf90_put_att(ncid, varid, "units","degrees")
+
ncfunc_retval = nf90_def_var(ncid, 'ylat0', NF90_DOUBLE, varid)
ncfunc_retval = nf90_put_var(ncid, varid, ylat0)
+! attributes
+ ncfunc_retval = nf90_put_att(ncid, varid, "description","latitude of the lowest left corner")
+ ncfunc_retval = nf90_put_att(ncid, varid, "units","degrees")
+
! All done, close the NetCDF file
ncfunc_retval = nf90_close(ncid)
@@ -295,18 +314,42 @@ CONTAINS
IF (nc_varname == 'U') THEN
ncfunc_retval = nf90_put_var(ncid, varid, &
& uu(0:nx-1, 0:ny-1, 1:nz, 1))
+! attributes
+ ncfunc_retval = nf90_put_att(ncid, varid, "description","U component of wind in the X[horizontal] direction")
+ ncfunc_retval = nf90_put_att(ncid, varid, "units","m s**-1")
+
+
ELSEIF (nc_varname == 'V') THEN
ncfunc_retval = nf90_put_var(ncid, varid, &
& vv(0:nx-1, 0:ny-1, 1:nz, 1))
+! attributes
+ ncfunc_retval = nf90_put_att(ncid, varid, "description","V component of wind in the Y[horizontal] direction")
+ ncfunc_retval = nf90_put_att(ncid, varid, "units","m s**-1")
+
+
ELSEIF (nc_varname == 'T') THEN
ncfunc_retval = nf90_put_var(ncid, varid, &
& tt(0:nx-1, 0:ny-1, 1:nz, 1))
+! attributes
+ ncfunc_retval = nf90_put_att(ncid, varid, "description","temperature")
+ ncfunc_retval = nf90_put_att(ncid, varid, "units","k")
+
ELSEIF (nc_varname == 'W') THEN
ncfunc_retval = nf90_put_var(ncid, varid, &
& ww(0:nx-1, 0:ny-1, 1:nz, 1))
+! attributes
+ ncfunc_retval = nf90_put_att(ncid, varid, "description","wind component in the Z[vertical] direction")
+ ncfunc_retval = nf90_put_att(ncid, varid, "units","m s**-1")
+
+
ELSEIF (nc_varname == 'Q') THEN
ncfunc_retval = nf90_put_var(ncid, varid, &
& qv(0:nx-1, 0:ny-1, 1:nz, 1))
+! attributes
+ ncfunc_retval = nf90_put_att(ncid, varid, "description","specific humidity")
+ ncfunc_retval = nf90_put_att(ncid, varid, "units"," ")
+
+
ELSE
PRINT *,
PRINT *, 'fp2nc4io:private_dump_3d_field() bad var: ', nc_varname
diff --git a/flexpart_code/grib2nc4/grib2nc4.F90 b/flexpart_code/grib2nc4/grib2nc4.F90
index 710d6dad..fe106d64 100644
--- a/flexpart_code/grib2nc4/grib2nc4.F90
+++ b/flexpart_code/grib2nc4/grib2nc4.F90
@@ -10,6 +10,11 @@ PROGRAM grib2nc4
! *
! May 2016 *
!*************************************************************************
+ ! M. Harustak *
+ ! -) modification to generate the output in single precission *
+ ! -) possibility to add a lat lon selection to obtain the met variables*
+ ! in the vertical levels defined in that location *
+ !*************************************************************************
USE par_mod
USE com_mod
@@ -20,20 +25,21 @@ PROGRAM grib2nc4
IMPLICIT NONE
- LOGICAL :: metfile_exists
+ LOGICAL :: metfile_exists, coordinates_provided, lat_provided, lon_provided
INTEGER :: i, j, k
INTEGER :: num_optional_vars, num_vars
INTEGER, PARAMETER :: DEFLATE_LEVEL = 2 ! Compression level (0-9)
- CHARACTER(LEN=512) :: met_filepath, netcdf4_filepath
+ CHARACTER(LEN=512) :: met_filepath, netcdf4_filepath, param_str, coord_name_str, coord_val_str
CHARACTER, DIMENSION(:), ALLOCATABLE :: vars_list ! list of variables
-
+ INTEGER :: coordX, coordY, stat
+ REAL :: coord_lat, coord_lon
INTEGER :: metdata_format = UNKNOWN_METDATA ! From FP par_mod
!--------------------------------------------------------
! Read in mandatory arguments
IF (IARGC() < 2) THEN
- PRINT *, 'Usage: grib2netcdf4 <inpath> <outpath> [optional varnames]'
+ PRINT *, 'Usage: grib2netcdf4 <inpath> <outpath> [lon=X lat=Y] [optional varnames]'
STOP
ELSE
CALL GETARG(1, met_filepath)
@@ -47,23 +53,47 @@ PROGRAM grib2nc4
! First, get the number of optional args and allocate vars_list,
! and fill the first three elements
- IF (IARGC() > 2) THEN
- num_optional_vars = IARGC() - 2
- ELSE
- num_optional_vars = 0
- ENDIF
+ coordinates_provided = .FALSE.
+ lat_provided = .FALSE.
+ lon_provided = .FALSE.
+ ALLOCATE( vars_list(IARGC()+3),stat=stat )
- num_vars = num_optional_vars + 3
- ALLOCATE( vars_list(num_vars) )
vars_list(1) = 'u'
vars_list(2) = 'v'
vars_list(3) = 't'
- ! Read in optional variable arguments, starting at element 4 of vars_list
- IF (IARGC() > 2) THEN
- DO i=1,num_optional_vars
- CALL GETARG( i+2, vars_list(i+3) )
- ENDDO
+ num_vars = 3
+ DO i=3,IARGC()
+ CALL GETARG(i,param_str)
+ param_str = TRIM(param_str)
+ j = SCAN(param_str,"=")
+ if (j>1) then
+ coord_name_str=param_str(1:j-1)
+ coord_val_str=param_str(j+1:)
+ IF ( coord_name_str == "lat" .or. coord_name_str == "LAT" ) THEN
+ read(coord_val_str,*,iostat=stat) coord_lat
+ if ( stat == 0 ) then
+ lat_provided = .TRUE.
+ else
+ print *, "Incorrect coordinates: ", coord_val_str
+ stop
+ endif
+ ELSE IF ( coord_name_str == "lon" .or. coord_name_str == "LON" ) THEN
+ read(coord_val_str,*,iostat=stat) coord_lon
+ if ( stat == 0 ) then
+ lon_provided = .TRUE.
+ else
+ print *, "Incorrect coordinates: ", coord_val_str
+ stop
+ endif
+ ENDIF
+ else
+ num_vars = num_vars + 1
+ vars_list(num_vars) = param_str
+ endif
+ ENDDO
+ IF (lat_provided .AND. lon_provided) THEN
+ coordinates_provided = .TRUE.
ENDIF
! Before proceeding, let's make sure the vars_list is good - otherwise,
@@ -125,8 +155,18 @@ PROGRAM grib2nc4
! to do so
!CALL gridcheck_nests
+ ! If the coordinates are provided, then we need to obtain
+ ! the corresponding grid indexes since this is what verttransform needs
+ if ( coordinates_provided ) then
+ coordX = (coord_lon-xlon0)/dx
+ coordY = (coord_lat-ylat0)/dy
+ print *, "Coordinates: "
+ print *, "lon: ", coord_lon, ", lat: ",coord_lat
+ print *, "x: ", coordX, ", y: ", coordY
+ endif
+
PRINT *, 'Calling processmetfields()...'
- call processmetfields( 1, metdata_format)
+ call processmetfields( 1, metdata_format, coordinates_provided, coordX, coordY)
PRINT *, 'Calling fp2nc4io_dump()...'
call fp2nc4io_dump( netcdf4_filepath, num_vars, vars_list, DEFLATE_LEVEL)
diff --git a/flexpart_code/grib2nc4/nc4cmp.F90 b/flexpart_code/grib2nc4/nc4cmp.F90
new file mode 100644
index 00000000..061ef516
--- /dev/null
+++ b/flexpart_code/grib2nc4/nc4cmp.F90
@@ -0,0 +1,244 @@
+PROGRAM nc4cmp
+
+ !*************************************************************************
+ ! This program compares variable in two NetCDF files *
+ ! (with dimension up to 3) and it prints average and maximal difference *
+ ! between the values. It allows for specification of optional tolerance *
+ ! *
+ ! Usage: *
+ ! *
+ ! ./nc4cmp <file1> <file2> <variable> [optional tolerance in %] *
+ ! *
+ ! M. Harustak *
+ ! *
+ ! March 2017 *
+ !*************************************************************************
+
+
+ USE netcdf
+
+ IMPLICIT NONE
+
+ LOGICAL :: file1_exists, file2_exists
+ CHARACTER(LEN=512) :: file1, file2, precision_str, variable
+ CHARACTER(LEN=NF90_MAX_NAME) :: returned_name
+ REAL :: tolerance
+ INTEGER :: ncid, varid, retval, xtype, numDims1, numDims2, i, j, k, stat
+ INTEGER, dimension(NF90_MAX_VAR_DIMS) :: varDimIDs1, varDimIDs2, dims1, dims2
+ double precision :: variable_array_0d_1, variable_array_0d_2
+ double precision, allocatable, dimension(:) :: variable_array_1d_1, variable_array_1d_2
+ double precision, allocatable, dimension(:,:) :: variable_array_2d_1, variable_array_2d_2
+ double precision, allocatable, dimension(:,:,:) :: variable_array_3d_1, variable_array_3d_2
+ double precision :: v1, v2
+ double precision :: sum_diff, max_diff
+ integer :: count_diff
+ logical :: different = .FALSE.
+
+ ! Read in mandatory arguments
+ IF (IARGC() < 3) THEN
+ PRINT *, 'Usage: nc4cmp <file1> <file2> <variable> [optional tolerance in %]'
+ STOP
+ ELSE
+ CALL GETARG(1, file1)
+ CALL GETARG(2, file2)
+ CALL GETARG(3, variable)
+ ENDIF
+
+ IF (IARGC() > 3) THEN
+ CALL GETARG( 4, precision_str)
+ READ(precision_str, '(f10.5)') tolerance
+ ELSE
+ tolerance = 0
+ ENDIF
+ print *,"Tolerance ",tolerance,"%"
+
+ INQUIRE( FILE=file1, EXIST=file1_exists )
+ IF ( file1_exists ) THEN
+ PRINT *, 'File 1: ', TRIM(file1)
+ ELSE
+ PRINT *, 'Unable to find file: ', TRIM(file1)
+ STOP
+ ENDIF
+ INQUIRE( FILE=file2, EXIST=file2_exists )
+ IF ( file2_exists ) THEN
+ PRINT *, 'File 2: ', TRIM(file2)
+ ELSE
+ PRINT *, 'Unable to find file: ', TRIM(file2)
+ STOP
+ ENDIF
+
+ print *, "Opening ",TRIM(file1)," for reading"
+ retval = nf90_open(file1, NF90_NOWRITE, ncid)
+ retval = nf90_inq_varid(ncid, TRIM(variable) , varid)
+ retval = nf90_inquire_variable(ncid, varid, ndims = numDims1, dimids = varDimIDs1)
+ print *, "Reading variable ",TRIM(variable), ", no of dimensions: ",numDims1
+ do i=1,numDims1
+ retval = nf90_inquire_dimension(ncid, varDimIDs1(i),len=j)
+ dims1(i) = j
+ print *," - dimension ",i,": ", j
+ enddo
+ if ( numDims1 == 0) then
+ retval = nf90_get_var(ncid, varid, variable_array_0d_1)
+ else if ( numDims1 == 1) then
+ allocate(variable_array_1d_1(dims1(1)),stat=stat)
+ retval = nf90_get_var(ncid, varid, variable_array_1d_1)
+ else if ( numDims1 == 2) then
+ allocate(variable_array_2d_1(dims1(1),dims1(2)),stat=stat)
+ retval = nf90_get_var(ncid, varid, variable_array_2d_1)
+ else if ( numDims1 == 3) then
+ allocate(variable_array_3d_1(dims1(1),dims1(2),dims1(3)),stat=stat)
+ retval = nf90_get_var(ncid, varid, variable_array_3d_1)
+ endif
+ retval = nf90_close(ncid)
+ print *, "...done"
+
+ print *, "Opening ",TRIM(file2)," for reading"
+ retval = nf90_open(file2, NF90_NOWRITE, ncid)
+ retval = nf90_inq_varid(ncid, TRIM(variable) , varid)
+ retval = nf90_inquire_variable(ncid, varid, ndims = numDims2, dimids = varDimIDs2)
+ print *, "Reading variable ",TRIM(variable), ", no of dimensions: ",numDims2
+ do i=1,numDims2
+ retval = nf90_inquire_dimension(ncid, varDimIDs2(i),len=j)
+ dims2(i) = j
+ print *," - dimension ",i,": ", j
+ enddo
+ if ( numDims2 == 0) then
+ retval = nf90_get_var(ncid, varid, variable_array_0d_2)
+ else if ( numDims2 == 1) then
+ allocate(variable_array_1d_2(dims2(1)),stat=stat)
+ retval = nf90_get_var(ncid, varid, variable_array_1d_2)
+ else if ( numDims2 == 2) then
+ allocate(variable_array_2d_2(dims2(1),dims2(2)),stat=stat)
+ retval = nf90_get_var(ncid, varid, variable_array_2d_2)
+ else if ( numDims2 == 3) then
+ allocate(variable_array_3d_2(dims2(1),dims2(2),dims2(3)),stat=stat)
+ retval = nf90_get_var(ncid, varid, variable_array_3d_2)
+ endif
+ retval = nf90_close(ncid)
+ print *, "...done"
+
+ if (numDims1 /= numDims2) THEN
+ print *, "Number of dimensions differs..."
+ STOP
+ endif
+ do i=1,numDims1
+ if (dims1(i) /= dims2(i)) then
+ print *, "Dimension ",i," differs..."
+ STOP
+ endif
+ enddo
+
+ if (numDims1>3) then
+ print *,"Only up to 3 dimensions supported..."
+ STOP
+ endif
+
+ sum_diff = 0.0
+ count_diff = 0
+
+ if (numDims1==0) then
+ v1 = variable_array_0d_1
+ v2 = variable_array_0d_2
+ print *, "Value 1: ",v1
+ print *, "Value 2: ",v2
+ if ( ( v1 == 0 .AND. v2 /= 0 ) .OR. ( v1 /= 0 .AND. v2 == 0 ) ) then
+ print *," === variables differ ==="
+ STOP
+ else if ( v1 == 0 .AND. v2 == 0) then
+ print *," === variables are within tolerance ==="
+ STOP
+ else if ( abs((v1-v2)/v1*100) > tolerance ) then
+ print *," === variables differ ==="
+ STOP
+ else
+ print *," === variables are within tolerance ==="
+ STOP
+ endif
+ endif
+ if (numDims1==1) then
+ do i=1,dims1(1)
+ v1 = variable_array_1d_1(i)
+ v2 = variable_array_1d_2(i)
+ if ( ( v1 == 0 .AND. v2 /= 0 ) .OR. ( v1 /= 0 .AND. v2 == 0 ) ) then
+ different = .TRUE.
+ else if (v1/=0 .AND. v2/=0 .AND. abs((v1-v2)/v1*100) > tolerance ) then
+ different = .TRUE.
+ endif
+ if ( v1 /= 0) then
+ sum_diff = sum_diff+abs((v1-v2)/v1*100)
+ count_diff = count_diff+1
+ if (abs((v1-v2)/v1*100)>max_diff) then
+ max_diff=abs((v1-v2)/v1*100)
+ endif
+ endif
+ enddo
+ if ( different ) then
+ print *," === variables differ ==="
+ else
+ print *," === variables are within tolerance ==="
+ endif
+ print *," === avg(diff%): ", sum_diff/count_diff
+ print *," === max(diff%): ",max_diff
+ STOP
+ endif
+ if (numDims1==2) then
+ do i=1,dims1(1)
+ do j=1,dims1(2)
+ v1 = variable_array_2d_1(i,j)
+ v2 = variable_array_2d_2(i,j)
+ if ( ( v1 == 0 .AND. v2 /= 0 ) .OR. ( v1 /= 0 .AND. v2 == 0 ) ) then
+ different = .TRUE.
+ else if (v1/=0 .AND. v2/=0 .AND. abs((v1-v2)/v1*100) > tolerance ) then
+ different = .TRUE.
+ endif
+ if (v1/=0) then
+ sum_diff = sum_diff+abs((v1-v2)/v1*100)
+ count_diff = count_diff+1
+ if (abs((v1-v2)/v1*100)>max_diff) then
+ max_diff=abs((v1-v2)/v1*100)
+ endif
+ endif
+ enddo
+ enddo
+ if ( different ) then
+ print *," === variables differ ==="
+ else
+ print *," === variables are within tolerance ==="
+ endif
+ print *," === avg(diff%): ", sum_diff/count_diff
+ print *," === max(diff%): ",max_diff
+ STOP
+ endif
+ if (numDims1==3) then
+ do i=1,dims1(1)
+ do j=1,dims1(2)
+ do k=1,dims1(3)
+ v1 = variable_array_3d_1(i,j,k)
+ v2 = variable_array_3d_2(i,j,k)
+ if ( ( v1 == 0 .AND. v2 /= 0 ) .OR. ( v1 /= 0 .AND. v2 == 0 ) ) then
+ different = .TRUE.
+ else if (v1/=0 .AND. v2/=0 .AND. abs((v1-v2)/v1*100) > tolerance ) then
+ different = .TRUE.
+ endif
+ if (v1/=0) then
+ sum_diff = sum_diff+abs((v1-v2)/v1*100)
+ count_diff = count_diff+1
+ if (abs((v1-v2)/v1*100)>max_diff) then
+ max_diff=abs((v1-v2)/v1*100)
+ endif
+ endif
+ enddo
+ enddo
+ enddo
+ if ( different ) then
+ print *," === variables differ ==="
+ else
+ print *," === variables are within tolerance ==="
+ endif
+ print *," === avg(diff%): ", sum_diff/count_diff
+ print *," === max(diff%): ",max_diff
+ STOP
+ endif
+
+
+END PROGRAM nc4cmp
diff --git a/flexpart_code/grib2nc4/processmetfields.F90 b/flexpart_code/grib2nc4/processmetfields.F90
index c9e2e2fc..ccdab7e9 100644
--- a/flexpart_code/grib2nc4/processmetfields.F90
+++ b/flexpart_code/grib2nc4/processmetfields.F90
@@ -1,4 +1,5 @@
-subroutine processmetfields(ind,metdata_format)
+subroutine processmetfields(ind,metdata_format,coordinates_provided, coordX, coordY)
+
! i o
!*****************************************************************************
! *
@@ -36,6 +37,8 @@ subroutine processmetfields(ind,metdata_format)
integer :: ind, metdata_format, i, lastSlash
integer :: dumpData
character(len=512):: fpfname, dumpPath, filename
+ integer :: coordX, coordY
+ logical :: coordinates_provided
!****************************************************************************
@@ -106,7 +109,11 @@ subroutine processmetfields(ind,metdata_format)
call readwind_nests(ind,memind(1),uuhn,vvhn,wwhn)
call calcpar_ecmwf(memind(1),uuh,vvh,pvh)
call calcpar_nests(memind(1),uuhn,vvhn,pvhn,metdata_format)
- call verttransform_ecmwf(memind(1),uuh,vvh,wwh,pvh)
+ if ( coordinates_provided ) then
+ call verttransform_grib2nc4_ecmwf(memind(1),uuh,vvh,wwh,pvh,coordX,coordY)
+ else
+ call verttransform_grib2nc4_ecmwf(memind(1),uuh,vvh,wwh,pvh,-1,-1)
+ endif
call verttransform_nests(memind(1),uuhn,vvhn,wwhn,pvhn)
memtime(1)=wftime(ind)
endif
@@ -115,7 +122,11 @@ subroutine processmetfields(ind,metdata_format)
call readwind_nests(ind,memind(1),uuhn,vvhn,wwhn)
call calcpar_gfs(memind(1),uuh,vvh,pvh)
call calcpar_nests(memind(1),uuhn,vvhn,pvhn,metdata_format)
- call verttransform_gfs(memind(1),uuh,vvh,wwh,pvh)
+ if ( coordinates_provided ) then
+ call verttransform_grib2nc4_gfs(memind(1),uuh,vvh,wwh,pvh,coordX,coordY)
+ else
+ call verttransform_grib2nc4_gfs(memind(1),uuh,vvh,wwh,pvh,-1,-1)
+ endif
call verttransform_nests(memind(1),uuhn,vvhn,wwhn,pvhn)
memtime(1)=wftime(ind)
endif
diff --git a/flexpart_code/grib2nc4/verttransform_grib2nc4_ecmwf.f90 b/flexpart_code/grib2nc4/verttransform_grib2nc4_ecmwf.f90
new file mode 100644
index 00000000..9228a91d
--- /dev/null
+++ b/flexpart_code/grib2nc4/verttransform_grib2nc4_ecmwf.f90
@@ -0,0 +1,701 @@
+!**********************************************************************
+! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 *
+! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, *
+! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann *
+! *
+! This file is part of FLEXPART. *
+! *
+! FLEXPART is free software: you can redistribute it and/or modify *
+! it under the terms of the GNU General Public License as published by*
+! the Free Software Foundation, either version 3 of the License, or *
+! (at your option) any later version. *
+! *
+! FLEXPART is distributed in the hope that it will be useful, *
+! but WITHOUT ANY WARRANTY; without even the implied warranty of *
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+! GNU General Public License for more details. *
+! *
+! You should have received a copy of the GNU General Public License *
+! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
+!**********************************************************************
+
+subroutine verttransform_grib2nc4_ecmwf(n,uuh,vvh,wwh,pvh,coordx_opt,coordy_opt)
+ ! i i i i i
+ !*****************************************************************************
+ ! *
+ ! This subroutine transforms temperature, dew point temperature and *
+ ! wind components from eta to meter coordinates. *
+ ! The vertical wind component is transformed from Pa/s to m/s using *
+ ! the conversion factor pinmconv. *
+ ! In addition, this routine calculates vertical density gradients *
+ ! needed for the parameterization of the turbulent velocities. *
+ ! *
+ ! Author: A. Stohl, G. Wotawa *
+ ! *
+ ! 12 August 1996 *
+ ! Update: 16 January 1998 *
+ ! *
+ ! Major update: 17 February 1999 *
+ ! by G. Wotawa *
+ ! *
+ ! - Vertical levels for u, v and w are put together *
+ ! - Slope correction for vertical velocity: Modification of calculation *
+ ! procedure *
+ ! *
+ !*****************************************************************************
+ ! Changes, Bernd C. Krueger, Feb. 2001:
+ ! Variables tth and qvh (on eta coordinates) from common block
+ !*****************************************************************************
+ ! Sabine Eckhardt, March 2007
+ ! added the variable cloud for use with scavenging - descr. in com_mod
+ !*****************************************************************************
+ ! Changes Arnold, D. and Morton, D. (2015): *
+ ! -- description of local and common variables *
+ !*****************************************************************************
+ ! July 2017 - M. Harustak *
+ ! In order to be used by GRIB2NC4, this routine has been modified so that *
+ ! it has the possibility to get two additional arguments, coordx_opt and *
+ ! coordy_opt. These two are used to define the location where the FLEXPART *
+ ! model levels are calculated and used in GRIB2NC4 to obtaine the met data *
+ ! *
+ ! The routine works normally if -1 (no value) is passed *
+ !*****************************************************************************
+
+
+ use par_mod
+ use com_mod
+ use cmapf_mod, only: cc2gll
+
+ implicit none
+
+ !***********************************************************************
+ ! Subroutine Parameters: *
+ ! input *
+ ! n temporal index for meteorological fields (1 to 3)*
+ ! uuh,vvh, wwh wind components in ECMWF model levels *
+ ! pvh potential vorticity in ECMWF model levels *
+ integer :: n
+ real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
+ real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
+ real :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)
+ real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
+ integer :: coordx_opt, coordy_opt
+ !***********************************************************************
+
+ !***********************************************************************
+ ! Local variables *
+ ! *
+ ! ew subroutine/function to calculate saturation *
+ ! water vaport for a given air temperature *
+ ! clouds(0:nxmax,0:nymax,0:nzmax,2) cloud field for wet deposition *
+ ! rain_cloud_above [0,1] whether there is a raining cloud or not *
+ ! kz_inv inverted indez for kz *
+ ! f_qvsat Saturation water vapor specific humidity (kg/kg) *
+ ! pressure *
+ ! rh relative humidity *
+ ! lsp large scale precip at one grid cell *
+ ! convp convective precip at one grid cell *
+ ! uvzlev height of the eta half-levels *
+ ! rhoh density in model levels *
+ ! pinmconv conversion factor *
+ ! pint pressure on model levels (using akz, bkz, ps) *
+ ! tv virtual temperature *
+ ! tvold,pold dummy variables to keep previous value *
+ ! dz1, dz2, dz differences heights model levels, pressure levels *
+ ! ui, vi *
+ ! xlon,xlat *
+ ! xlonr xlon*pi/180. *
+ ! dzdx,dzdy slope corrections *
+ ! dzdx1,dzdx2,dzdy1 slope corrections in each direction *
+ ! uuaux,vvaux,uupolaux,vvpolaux auxiliary variables for polar sterero.*
+ ! ddpol,ffpol for special treatment of poles *
+ ! wdummy *
+ ! wzlev *
+ ! uvwzlev *
+
+ integer :: rain_cloud_above,kz_inv
+ real :: f_qvsat,pressure
+ real :: rh,lsp,convp
+ real :: uvzlev(nuvzmax),rhoh(nuvzmax),pinmconv(nzmax)
+ real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
+ real :: xlon,ylat,xlonr,dzdx,dzdy
+ real :: dzdx1,dzdx2,dzdy1,dzdy2
+ real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy
+ real :: wzlev(nwzmax),uvwzlev(0:nxmax-1,0:nymax-1,nzmax)
+
+ ! Other variables:
+ ! ix,jy,kz,kzmin loop control indices in each direction *
+ ! kl,klp
+ ! ix1,jy1,ixp,jyp,ixm,jym
+ integer :: ix,jy,kz,iz,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym
+ !***********************************************************************
+
+ !***********************************************************************
+ ! Local constants *
+ real,parameter :: const=r_air/ga
+ !***********************************************************************
+
+ !***********************************************************************
+ ! Global variables *
+ ! from com_mod *
+ ! nxmin1,nymin1 nx-1, ny-1, respectively *
+ ! tt2 2 meter temperature *
+ ! ps surface pressure *
+ ! td2 2 meter dew point *
+ ! height heights of all levels *
+ ! akm, bkm coeffs. which regulate vertical discretization of ecmwf *
+ ! akz, bkz model discretization coeffizients at the centre of layers *
+ ! nuvz,nwz vertical dimension of original ECMWF data *
+ ! nx,ny,nz actual dimensions of wind fields in x,y and z*
+ ! aknew,bknew model discretization coeffs. at the interpolated levels *
+ ! uu,vv,ww [m/2] wind components in x,y and z direction *
+ ! qv specific humidity data *
+ ! pv (pvu) potential vorticity *
+ ! rho [kg/m3] air density *
+ ! drhodz [kg/m2] vertical air density gradient *
+ ! dxconst,dyconst auxiliary variables for utransform,vtransform*
+ ! ylat0 geographical latitude of lower left grid point *
+ ! xlon0 geographical longitude of lower left grid point*
+ ! uupol,vvpol [m/s] wind components in polar stereographic projection*
+ ! clouds(0:nxmax,0:nymax,0:nzmax,2) cloud field for wet deposition *
+ ! lsprec [mm/h] large scale total precipitation *
+ ! convprec [mm/h] convective precipitation *
+ ! cloudsh *
+ ! *
+ !***********************************************************************
+
+!-----------------------------------------------------------------------------
+
+
+
+ logical :: init = .true.
+
+
+ !*************************************************************************
+ ! If verttransform is called the first time, initialize heights of the *
+ ! z levels in meter. The heights are the heights of model levels, where *
+ ! u,v,T and qv are given, and of the interfaces, where w is given. So, *
+ ! the vertical resolution in the z system is doubled. As reference point,*
+ ! the lower left corner of the grid is used. *
+ ! Unlike in the eta system, no difference between heights for u,v and *
+ ! heights for w exists. *
+ !*************************************************************************
+
+
+ ! do 897 kz=1,nuvz
+ ! write (*,*) 'akz: ',akz(kz),'bkz',bkz(kz)
+ !897 continue
+
+ if (init) then
+
+ ! Search for a point with high surface pressure (i.e. not above significant topography)
+ ! Then, use this point to construct a reference z profile, to be used at all times
+ !*****************************************************************************
+
+ ! Alternatively, and to be used in GRIB2NC4, the reference point can be
+ ! provided by the user
+
+ if ( coordx_opt >= 0 .and. coordy_opt >= 0 ) then
+ print *, "Using user provided coordinates: ", coordx_opt,", ",coordy_opt
+ ixm = coordx_opt
+ jym = coordy_opt
+ else
+ print *, "Searching for coordinates"
+ do jy=0,nymin1
+ do ix=0,nxmin1
+ if (ps(ix,jy,1,n).gt.100000.) then
+ ixm=ix
+ jym=jy
+ goto 3
+ endif
+ end do
+ end do
+3 continue
+ endif
+
+
+ tvold=tt2(ixm,jym,1,n)*(1.+0.378*ew(td2(ixm,jym,1,n))/ &
+ ps(ixm,jym,1,n))
+ pold=ps(ixm,jym,1,n)
+ height(1)=0.
+
+ do kz=2,nuvz
+ pint=akz(kz)+bkz(kz)*ps(ixm,jym,1,n)
+ tv=tth(ixm,jym,kz,n)*(1.+0.608*qvh(ixm,jym,kz,n))
+
+
+ ! NOTE: In FLEXPART versions up to 4.0, the number of model levels was doubled
+ ! upon the transformation to z levels. In order to save computer memory, this is
+ ! not done anymore in the standard version. However, this option can still be
+ ! switched on by replacing the following lines with those below, that are
+ ! currently commented out.
+ ! Note that two more changes are necessary in this subroutine below.
+ ! One change is also necessary in gridcheck.f, and another one in verttransform_nests.
+ !*****************************************************************************
+
+ if (abs(tv-tvold).gt.0.2) then
+ height(kz)= &
+ height(kz-1)+const*log(pold/pint)* &
+ (tv-tvold)/log(tv/tvold)
+ else
+ height(kz)=height(kz-1)+ &
+ const*log(pold/pint)*tv
+ endif
+
+ ! Switch on following lines to use doubled vertical resolution
+ !*************************************************************
+ ! if (abs(tv-tvold).gt.0.2) then
+ ! height((kz-1)*2)=
+ ! + height(max((kz-2)*2,1))+const*log(pold/pint)*
+ ! + (tv-tvold)/log(tv/tvold)
+ ! else
+ ! height((kz-1)*2)=height(max((kz-2)*2,1))+
+ ! + const*log(pold/pint)*tv
+ ! endif
+ ! End doubled vertical resolution
+
+ tvold=tv
+ pold=pint
+ end do
+
+
+ ! Switch on following lines to use doubled vertical resolution
+ !*************************************************************
+ ! do 7 kz=3,nz-1,2
+ ! height(kz)=0.5*(height(kz-1)+height(kz+1))
+ ! height(nz)=height(nz-1)+height(nz-1)-height(nz-2)
+ ! End doubled vertical resolution
+
+
+ ! Determine highest levels that can be within PBL
+ !************************************************
+
+ do kz=1,nz
+ if (height(kz).gt.hmixmax) then
+ nmixz=kz
+ goto 9
+ endif
+ end do
+9 continue
+
+ ! Do not repeat initialization of the Cartesian z grid
+ !*****************************************************
+
+ init=.false.
+
+ endif
+
+
+ ! Loop over the whole grid
+ !*************************
+
+ do jy=0,nymin1
+ do ix=0,nxmin1
+ tvold=tt2(ix,jy,1,n)*(1.+0.378*ew(td2(ix,jy,1,n))/ &
+ ps(ix,jy,1,n))
+ pold=ps(ix,jy,1,n)
+ uvzlev(1)=0.
+ wzlev(1)=0.
+ rhoh(1)=pold/(r_air*tvold)
+
+
+ ! Compute heights of eta levels
+ !******************************
+
+ do kz=2,nuvz
+ pint=akz(kz)+bkz(kz)*ps(ix,jy,1,n)
+ tv=tth(ix,jy,kz,n)*(1.+0.608*qvh(ix,jy,kz,n))
+ rhoh(kz)=pint/(r_air*tv)
+
+ if (abs(tv-tvold).gt.0.2) then
+ uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)* &
+ (tv-tvold)/log(tv/tvold)
+ else
+ uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)*tv
+ endif
+
+ tvold=tv
+ pold=pint
+ end do
+
+
+ do kz=2,nwz-1
+ wzlev(kz)=(uvzlev(kz+1)+uvzlev(kz))/2.
+ end do
+ wzlev(nwz)=wzlev(nwz-1)+ &
+ uvzlev(nuvz)-uvzlev(nuvz-1)
+
+ uvwzlev(ix,jy,1)=0.0
+ do kz=2,nuvz
+ uvwzlev(ix,jy,kz)=uvzlev(kz)
+ end do
+
+ ! Switch on following lines to use doubled vertical resolution
+ ! Switch off the three lines above.
+ !*************************************************************
+ !22 uvwzlev(ix,jy,(kz-1)*2)=uvzlev(kz)
+ ! do 23 kz=2,nwz
+ !23 uvwzlev(ix,jy,(kz-1)*2+1)=wzlev(kz)
+ ! End doubled vertical resolution
+
+ ! pinmconv=(h2-h1)/(p2-p1)
+
+ pinmconv(1)=(uvwzlev(ix,jy,2)-uvwzlev(ix,jy,1))/ &
+ ((aknew(2)+bknew(2)*ps(ix,jy,1,n))- &
+ (aknew(1)+bknew(1)*ps(ix,jy,1,n)))
+ do kz=2,nz-1
+ pinmconv(kz)=(uvwzlev(ix,jy,kz+1)-uvwzlev(ix,jy,kz-1))/ &
+ ((aknew(kz+1)+bknew(kz+1)*ps(ix,jy,1,n))- &
+ (aknew(kz-1)+bknew(kz-1)*ps(ix,jy,1,n)))
+ end do
+ pinmconv(nz)=(uvwzlev(ix,jy,nz)-uvwzlev(ix,jy,nz-1))/ &
+ ((aknew(nz)+bknew(nz)*ps(ix,jy,1,n))- &
+ (aknew(nz-1)+bknew(nz-1)*ps(ix,jy,1,n)))
+
+ ! Levels, where u,v,t and q are given
+ !************************************
+
+ uu(ix,jy,1,n)=uuh(ix,jy,1)
+ vv(ix,jy,1,n)=vvh(ix,jy,1)
+ tt(ix,jy,1,n)=tth(ix,jy,1,n)
+ qv(ix,jy,1,n)=qvh(ix,jy,1,n)
+ pv(ix,jy,1,n)=pvh(ix,jy,1)
+ rho(ix,jy,1,n)=rhoh(1)
+ uu(ix,jy,nz,n)=uuh(ix,jy,nuvz)
+ vv(ix,jy,nz,n)=vvh(ix,jy,nuvz)
+ tt(ix,jy,nz,n)=tth(ix,jy,nuvz,n)
+ qv(ix,jy,nz,n)=qvh(ix,jy,nuvz,n)
+ pv(ix,jy,nz,n)=pvh(ix,jy,nuvz)
+ rho(ix,jy,nz,n)=rhoh(nuvz)
+ kmin=2
+ do iz=2,nz-1
+ do kz=kmin,nuvz
+ if(height(iz).gt.uvzlev(nuvz)) then
+ uu(ix,jy,iz,n)=uu(ix,jy,nz,n)
+ vv(ix,jy,iz,n)=vv(ix,jy,nz,n)
+ tt(ix,jy,iz,n)=tt(ix,jy,nz,n)
+ qv(ix,jy,iz,n)=qv(ix,jy,nz,n)
+ pv(ix,jy,iz,n)=pv(ix,jy,nz,n)
+ rho(ix,jy,iz,n)=rho(ix,jy,nz,n)
+ goto 30
+ endif
+ if ((height(iz).gt.uvzlev(kz-1)).and. &
+ (height(iz).le.uvzlev(kz))) then
+ dz1=height(iz)-uvzlev(kz-1)
+ dz2=uvzlev(kz)-height(iz)
+ dz=dz1+dz2
+ uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)*dz2+uuh(ix,jy,kz)*dz1)/dz
+ vv(ix,jy,iz,n)=(vvh(ix,jy,kz-1)*dz2+vvh(ix,jy,kz)*dz1)/dz
+ tt(ix,jy,iz,n)=(tth(ix,jy,kz-1,n)*dz2 &
+ +tth(ix,jy,kz,n)*dz1)/dz
+ qv(ix,jy,iz,n)=(qvh(ix,jy,kz-1,n)*dz2 &
+ +qvh(ix,jy,kz,n)*dz1)/dz
+ pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)*dz2+pvh(ix,jy,kz)*dz1)/dz
+ rho(ix,jy,iz,n)=(rhoh(kz-1)*dz2+rhoh(kz)*dz1)/dz
+ kmin=kz
+ goto 30
+ endif
+ end do
+30 continue
+ end do
+
+
+ ! Levels, where w is given
+ !*************************
+
+ ww(ix,jy,1,n)=wwh(ix,jy,1)*pinmconv(1)
+ ww(ix,jy,nz,n)=wwh(ix,jy,nwz)*pinmconv(nz)
+ kmin=2
+ do iz=2,nz
+ do kz=kmin,nwz
+ if ((height(iz).gt.wzlev(kz-1)).and. &
+ (height(iz).le.wzlev(kz))) then
+ dz1=height(iz)-wzlev(kz-1)
+ dz2=wzlev(kz)-height(iz)
+ dz=dz1+dz2
+ ww(ix,jy,iz,n)=(wwh(ix,jy,kz-1)*pinmconv(kz-1)*dz2 &
+ +wwh(ix,jy,kz)*pinmconv(kz)*dz1)/dz
+ kmin=kz
+ goto 40
+ endif
+ end do
+40 continue
+ end do
+
+ ! Compute density gradients at intermediate levels
+ !*************************************************
+
+ drhodz(ix,jy,1,n)=(rho(ix,jy,2,n)-rho(ix,jy,1,n))/ &
+ (height(2)-height(1))
+ do kz=2,nz-1
+ drhodz(ix,jy,kz,n)=(rho(ix,jy,kz+1,n)-rho(ix,jy,kz-1,n))/ &
+ (height(kz+1)-height(kz-1))
+ end do
+ drhodz(ix,jy,nz,n)=drhodz(ix,jy,nz-1,n)
+
+ end do
+ end do
+
+
+ !****************************************************************
+ ! Compute slope of eta levels in windward direction and resulting
+ ! vertical wind correction
+ !****************************************************************
+
+ do jy=1,ny-2
+ do ix=1,nx-2
+
+ kmin=2
+ do iz=2,nz-1
+
+ ui=uu(ix,jy,iz,n)*dxconst/cos((real(jy)*dy+ylat0)*pi180)
+ vi=vv(ix,jy,iz,n)*dyconst
+
+ do kz=kmin,nz
+ if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
+ (height(iz).le.uvwzlev(ix,jy,kz))) then
+ dz1=height(iz)-uvwzlev(ix,jy,kz-1)
+ dz2=uvwzlev(ix,jy,kz)-height(iz)
+ dz=dz1+dz2
+ kl=kz-1
+ klp=kz
+ kmin=kz
+ goto 47
+ endif
+ end do
+
+47 ix1=ix-1
+ jy1=jy-1
+ ixp=ix+1
+ jyp=jy+1
+
+ dzdx1=(uvwzlev(ixp,jy,kl)-uvwzlev(ix1,jy,kl))/2.
+ dzdx2=(uvwzlev(ixp,jy,klp)-uvwzlev(ix1,jy,klp))/2.
+ dzdx=(dzdx1*dz2+dzdx2*dz1)/dz
+
+ dzdy1=(uvwzlev(ix,jyp,kl)-uvwzlev(ix,jy1,kl))/2.
+ dzdy2=(uvwzlev(ix,jyp,klp)-uvwzlev(ix,jy1,klp))/2.
+ dzdy=(dzdy1*dz2+dzdy2*dz1)/dz
+
+ ww(ix,jy,iz,n)=ww(ix,jy,iz,n)+(dzdx*ui+dzdy*vi)
+
+ end do
+
+ end do
+ end do
+
+
+ ! If north pole is in the domain, calculate wind velocities in polar
+ ! stereographic coordinates
+ !*******************************************************************
+
+ if (nglobal) then
+ do jy=int(switchnorthg)-2,nymin1
+ ylat=ylat0+real(jy)*dy
+ do ix=0,nxmin1
+ xlon=xlon0+real(ix)*dx
+ do iz=1,nz
+ call cc2gll(northpolemap,ylat,xlon,uu(ix,jy,iz,n), &
+ vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
+ vvpol(ix,jy,iz,n))
+ end do
+ end do
+ end do
+
+
+ do iz=1,nz
+
+ ! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
+ !
+ ! AMSnauffer Nov 18 2004 Added check for case vv=0
+ !
+ xlon=xlon0+real(nx/2-1)*dx
+ xlonr=xlon*pi/180.
+ ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+ &
+ vv(nx/2-1,nymin1,iz,n)**2)
+ if (vv(nx/2-1,nymin1,iz,n).lt.0.) then
+ ddpol=atan(uu(nx/2-1,nymin1,iz,n)/ &
+ vv(nx/2-1,nymin1,iz,n))-xlonr
+ else if (vv(nx/2-1,nymin1,iz,n).gt.0.) then
+ ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/ &
+ vv(nx/2-1,nymin1,iz,n))-xlonr
+ else
+ ddpol=pi/2-xlonr
+ endif
+ if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
+ if(ddpol.gt.2.0*pi) ddpol=ddpol-2.0*pi
+
+ ! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
+ xlon=180.0
+ xlonr=xlon*pi/180.
+ ylat=90.0
+ uuaux=-ffpol*sin(xlonr+ddpol)
+ vvaux=-ffpol*cos(xlonr+ddpol)
+ call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
+ vvpolaux)
+
+ jy=nymin1
+ do ix=0,nxmin1
+ uupol(ix,jy,iz,n)=uupolaux
+ vvpol(ix,jy,iz,n)=vvpolaux
+ end do
+ end do
+
+
+ ! Fix: Set W at pole to the zonally averaged W of the next equator-
+ ! ward parallel of latitude
+
+ do iz=1,nz
+ wdummy=0.
+ jy=ny-2
+ do ix=0,nxmin1
+ wdummy=wdummy+ww(ix,jy,iz,n)
+ end do
+ wdummy=wdummy/real(nx)
+ jy=nymin1
+ do ix=0,nxmin1
+ ww(ix,jy,iz,n)=wdummy
+ end do
+ end do
+
+ endif
+
+
+ ! If south pole is in the domain, calculate wind velocities in polar
+ ! stereographic coordinates
+ !*******************************************************************
+
+ if (sglobal) then
+ do jy=0,int(switchsouthg)+3
+ ylat=ylat0+real(jy)*dy
+ do ix=0,nxmin1
+ xlon=xlon0+real(ix)*dx
+ do iz=1,nz
+ call cc2gll(southpolemap,ylat,xlon,uu(ix,jy,iz,n), &
+ vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
+ vvpol(ix,jy,iz,n))
+ end do
+ end do
+ end do
+
+ do iz=1,nz
+
+ ! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
+ !
+ ! AMSnauffer Nov 18 2004 Added check for case vv=0
+ !
+ xlon=xlon0+real(nx/2-1)*dx
+ xlonr=xlon*pi/180.
+ ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+ &
+ vv(nx/2-1,0,iz,n)**2)
+ if (vv(nx/2-1,0,iz,n).lt.0.) then
+ ddpol=atan(uu(nx/2-1,0,iz,n)/ &
+ vv(nx/2-1,0,iz,n))+xlonr
+ else if (vv(nx/2-1,0,iz,n).gt.0.) then
+ ddpol=pi+atan(uu(nx/2-1,0,iz,n)/ &
+ vv(nx/2-1,0,iz,n))+xlonr
+ else
+ ddpol=pi/2-xlonr
+ endif
+ if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
+ if(ddpol.gt.2.0*pi) ddpol=ddpol-2.0*pi
+
+ ! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
+ xlon=180.0
+ xlonr=xlon*pi/180.
+ ylat=-90.0
+ uuaux=+ffpol*sin(xlonr-ddpol)
+ vvaux=-ffpol*cos(xlonr-ddpol)
+ call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
+ vvpolaux)
+
+ jy=0
+ do ix=0,nxmin1
+ uupol(ix,jy,iz,n)=uupolaux
+ vvpol(ix,jy,iz,n)=vvpolaux
+ end do
+ end do
+
+
+ ! Fix: Set W at pole to the zonally averaged W of the next equator-
+ ! ward parallel of latitude
+
+ do iz=1,nz
+ wdummy=0.
+ jy=1
+ do ix=0,nxmin1
+ wdummy=wdummy+ww(ix,jy,iz,n)
+ end do
+ wdummy=wdummy/real(nx)
+ jy=0
+ do ix=0,nxmin1
+ ww(ix,jy,iz,n)=wdummy
+ end do
+ end do
+ endif
+
+
+ !write (*,*) 'initializing clouds, n:',n,nymin1,nxmin1,nz
+ ! create a cloud and rainout/washout field, clouds occur where rh>80%
+ ! total cloudheight is stored at level 0
+ do jy=0,nymin1
+ do ix=0,nxmin1
+ rain_cloud_above=0
+ lsp=lsprec(ix,jy,1,n)
+ convp=convprec(ix,jy,1,n)
+ cloudsh(ix,jy,n)=0
+ do kz_inv=1,nz-1
+ kz=nz-kz_inv+1
+ pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
+ rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
+ clouds(ix,jy,kz,n)=0
+ if (rh.gt.0.8) then ! in cloud
+ if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
+ rain_cloud_above=1
+ cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
+ height(kz)-height(kz-1)
+ if (lsp.ge.convp) then
+ clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
+ else
+ clouds(ix,jy,kz,n)=2 ! convp dominated rainout
+ endif
+ else ! no precipitation
+ clouds(ix,jy,kz,n)=1 ! cloud
+ endif
+ else ! no cloud
+ if (rain_cloud_above.eq.1) then ! scavenging
+ if (lsp.ge.convp) then
+ clouds(ix,jy,kz,n)=5 ! lsp dominated washout
+ else
+ clouds(ix,jy,kz,n)=4 ! convp dominated washout
+ endif
+ endif
+ endif
+ end do
+ end do
+ end do
+
+ !do 102 kz=1,nuvz
+ !write(an,'(i02)') kz+10
+ !write(*,*) nuvz,nymin1,nxmin1,'--',an,'--'
+ !open(4,file='/nilu_wrk2/sec/cloudtest/cloud'//an,form='formatted')
+ !do 101 jy=0,nymin1
+ ! write(4,*) (clouds(ix,jy,kz,n),ix=1,nxmin1)
+ !101 continue
+ ! close(4)
+ !102 continue
+
+ ! open(4,file='/nilu_wrk2/sec/cloudtest/height',form='formatted')
+ ! do 103 jy=0,nymin1
+ ! write (4,*)
+ !+ (height(kz),kz=1,nuvz)
+ !103 continue
+ ! close(4)
+
+ !open(4,file='/nilu_wrk2/sec/cloudtest/p',form='formatted')
+ ! do 104 jy=0,nymin1
+ ! write (4,*)
+ !+ (r_air*tt(ix,jy,1,n)*rho(ix,jy,1,n),ix=1,nxmin1)
+ !104 continue
+ ! close(4)
+end subroutine verttransform_grib2nc4_ecmwf
diff --git a/flexpart_code/grib2nc4/verttransform_grib2nc4_gfs.f90 b/flexpart_code/grib2nc4/verttransform_grib2nc4_gfs.f90
new file mode 100644
index 00000000..ccd94f79
--- /dev/null
+++ b/flexpart_code/grib2nc4/verttransform_grib2nc4_gfs.f90
@@ -0,0 +1,680 @@
+!**********************************************************************
+! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 *
+! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, *
+! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann *
+! *
+! This file is part of FLEXPART. *
+! *
+! FLEXPART is free software: you can redistribute it and/or modify *
+! it under the terms of the GNU General Public License as published by*
+! the Free Software Foundation, either version 3 of the License, or *
+! (at your option) any later version. *
+! *
+! FLEXPART is distributed in the hope that it will be useful, *
+! but WITHOUT ANY WARRANTY; without even the implied warranty of *
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+! GNU General Public License for more details. *
+! *
+! You should have received a copy of the GNU General Public License *
+! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
+!**********************************************************************
+
+subroutine verttransform_grib2nc4_gfs(n,uuh,vvh,wwh,pvh,coordx_opt,coordy_opt)
+ ! i i i i i
+ !*****************************************************************************
+ ! *
+ ! This subroutine transforms temperature, dew point temperature and *
+ ! wind components from eta to meter coordinates. *
+ ! The vertical wind component is transformed from Pa/s to m/s using *
+ ! the conversion factor pinmconv. *
+ ! In addition, this routine calculates vertical density gradients *
+ ! needed for the parameterization of the turbulent velocities. *
+ ! *
+ ! Author: A. Stohl, G. Wotawa *
+ ! *
+ ! 12 August 1996 *
+ ! Update: 16 January 1998 *
+ ! *
+ ! Major update: 17 February 1999 *
+ ! by G. Wotawa *
+ ! CHANGE 17/11/2005 Caroline Forster, NCEP GFS version *
+ ! *
+ ! - Vertical levels for u, v and w are put together *
+ ! - Slope correction for vertical velocity: Modification of calculation *
+ ! procedure *
+ ! *
+ !*****************************************************************************
+ ! Changes, Bernd C. Krueger, Feb. 2001:
+ ! Variables tth and qvh (on eta coordinates) from common block
+ !*****************************************************************************
+ ! Changes Arnold, D. and Morton, D. (2015): *
+ ! -- description of local and common variables *
+ !*****************************************************************************
+ ! July 2017 - M. Harustak *
+ ! In order to be used by GRIB2NC4, this routine has been modified so that *
+ ! it has the possibility to get two additional arguments, coordx_opt and *
+ ! coordy_opt. These two are used to define the location where the FLEXPART *
+ ! model levels are calculated and used in GRIB2NC4 to obtaine the met data *
+ ! *
+ ! The routine works normally if -1 (no value) is passed *
+ !*****************************************************************************
+
+
+ use par_mod
+ use com_mod
+ use cmapf_mod
+
+ implicit none
+ !***********************************************************************
+ ! Subroutine Parameters: *
+ ! input *
+ ! n temporal index for meteorological fields (1 to 3)*
+ ! uuh,vvh, wwh wind components in ECMWF model levels *
+ ! pvh potential vorticity in ECMWF model levels *
+ integer :: n
+ real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
+ real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
+ real :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)
+ real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
+ integer :: coordx_opt, coordy_opt
+ !***********************************************************************
+
+ !***********************************************************************
+ ! Local variables *
+ ! *
+ ! ew subroutine/function to calculate saturation *
+ ! water vaport for a given air temperature *
+ ! clouds(0:nxmax,0:nymax,0:nzmax,2) cloud field for wet deposition *
+ ! rain_cloud_above [0,1] whether there is a raining cloud or not *
+ ! kz_inv inverted indez for kz *
+ ! f_qvsat Saturation water vapor specific humidity (kg/kg) *
+ ! pressure *
+ ! rh relative humidity *
+ ! lsp large scale precip at one grid cell *
+ ! convp convective precip at one grid cell *
+ ! uvzlev height of the eta half-levels *
+ ! rhoh density in model levels *
+ ! pinmconv conversion factor *
+ ! pint pressure on model levels (using akz, bkz, ps) *
+ ! tv virtual temperature *
+ ! tvold,pold dummy variables to keep previous value *
+ ! dz1, dz2, dz differences heights model levels, pressure levels *
+ ! ui, vi *
+ ! xlon,xlat *
+ ! xlonr xlon*pi/180. *
+ ! dzdx,dzdy slope corrections *
+ ! dzdx1,dzdx2,dzdy1 slope corrections in each direction *
+ ! uuaux,vvaux,uupolaux,vvpolaux auxiliary variables for polar sterero.*
+ ! ddpol,ffpol for special treatment of poles *
+ ! wdummy *
+ ! wzlev *
+ ! uvwzlev *
+
+ integer :: rain_cloud_above,kz_inv
+ real :: f_qvsat,pressure
+ real :: rh,lsp,convp
+ real :: uvzlev(nuvzmax),rhoh(nuvzmax),pinmconv(nzmax)
+ real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
+ real :: xlon,ylat,xlonr,dzdx,dzdy
+ real :: dzdx1,dzdx2,dzdy1,dzdy2
+ real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy
+ real :: wzlev(nwzmax),uvwzlev(0:nxmax-1,0:nymax-1,nzmax)
+
+ ! NCEP version
+ integer :: llev, i
+
+ ! Other variables:
+ ! ix,jy,kz,kzmin loop control indices in each direction *
+ ! kl,klp
+ ! ix1,jy1,ixp,jyp,ixm,jym
+ integer :: ix,jy,kz,iz,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym
+ !***********************************************************************
+
+ !***********************************************************************
+ ! Local constants *
+ real,parameter :: const=r_air/ga
+ !***********************************************************************
+
+ !***********************************************************************
+ ! Global variables *
+ ! from par_mod and com_mod *
+ ! nxmin1,nymin1 nx-1, ny-1, respectively *
+ ! tt2 2 meter temperature *
+ ! ps surface pressure *
+ ! td2 2 meter dew point *
+ ! height heights of all levels *
+ ! akm, bkm coeffs. which regulate vertical discretization of ecmwf *
+ ! akz, bkz model discretization coeffizients at the centre of layers *
+ ! tv virtual temperature *
+ ! nuvz,nwz vertical dimension of original ECMWF data *
+ ! nx,ny,nz actual dimensions of wind fields in x,y and z*
+ ! aknew,bknew model discretization coeffs. at the interpolated levels *
+ ! uu,vv,ww [m/2] wind components in x,y and z direction *
+ ! qv specific humidity data *
+ ! pv (pvu) potential vorticity *
+ ! rho [kg/m3] air density *
+ ! drhodz [kg/m2] vertical air density gradient *
+ ! dxconst,dyconst auxiliary variables for utransform,vtransform*
+ ! ylat0 geographical latitude of lower left grid point *
+ ! xlon0 geographical longitude of lower left grid point*
+ ! uupol,vvpol [m/s] wind components in polar stereographic projection*
+ ! clouds(0:nxmax,0:nymax,0:nzmax,2) cloud field for wet deposition *
+ ! lsprec [mm/h] large scale total precipitation *
+ ! convprec [mm/h] convective precipitation *
+ ! cloudsh *
+ ! *
+ !***********************************************************************
+
+!-----------------------------------------------------------------------------
+
+
+ logical :: init = .true.
+
+
+ !*************************************************************************
+ ! If verttransform is called the first time, initialize heights of the *
+ ! z levels in meter. The heights are the heights of model levels, where *
+ ! u,v,T and qv are given, and of the interfaces, where w is given. So, *
+ ! the vertical resolution in the z system is doubled. As reference point,*
+ ! the lower left corner of the grid is used. *
+ ! Unlike in the eta system, no difference between heights for u,v and *
+ ! heights for w exists. *
+ !*************************************************************************
+
+ if (init) then
+
+ ! Search for a point with high surface pressure (i.e. not above significant topography)
+ ! Then, use this point to construct a reference z profile, to be used at all times
+ !*****************************************************************************
+
+
+ if ( coordx_opt >= 0 .and. coordy_opt>=0 ) then
+ print *, "Using user provided coordinates: ", coordx_opt,", ",coordy_opt
+ ixm = coordx_opt
+ jym = coordy_opt
+ else
+ do jy=0,nymin1
+ do ix=0,nxmin1
+ if (ps(ix,jy,1,n).gt.100000.) then
+ ixm=ix
+ jym=jy
+ goto 3
+ endif
+ end do
+ end do
+3 continue
+ endif
+
+
+ tvold=tt2(ixm,jym,1,n)*(1.+0.378*ew(td2(ixm,jym,1,n))/ &
+ ps(ixm,jym,1,n))
+ pold=ps(ixm,jym,1,n)
+ height(1)=0.
+
+ do kz=2,nuvz
+ pint=akz(kz)+bkz(kz)*ps(ixm,jym,1,n)
+ tv=tth(ixm,jym,kz,n)*(1.+0.608*qvh(ixm,jym,kz,n))
+
+
+ ! NOTE: In FLEXPART versions up to 4.0, the number of model levels was doubled
+ ! upon the transformation to z levels. In order to save computer memory, this is
+ ! not done anymore in the standard version. However, this option can still be
+ ! switched on by replacing the following lines with those below, that are
+ ! currently commented out.
+ ! Note that two more changes are necessary in this subroutine below.
+ ! One change is also necessary in gridcheck.f, and another one in verttransform_nests.
+ !*****************************************************************************
+
+ if (abs(tv-tvold).gt.0.2) then
+ height(kz)= &
+ height(kz-1)+const*log(pold/pint)* &
+ (tv-tvold)/log(tv/tvold)
+ else
+ height(kz)=height(kz-1)+ &
+ const*log(pold/pint)*tv
+ endif
+
+ ! Switch on following lines to use doubled vertical resolution
+ !*************************************************************
+ ! if (abs(tv-tvold).gt.0.2) then
+ ! height((kz-1)*2)=
+ ! + height(max((kz-2)*2,1))+const*log(pold/pint)*
+ ! + (tv-tvold)/log(tv/tvold)
+ ! else
+ ! height((kz-1)*2)=height(max((kz-2)*2,1))+
+ ! + const*log(pold/pint)*tv
+ ! endif
+ ! End doubled vertical resolution
+
+ tvold=tv
+ pold=pint
+ end do
+
+
+ ! Switch on following lines to use doubled vertical resolution
+ !*************************************************************
+ ! do 7 kz=3,nz-1,2
+ ! height(kz)=0.5*(height(kz-1)+height(kz+1))
+ ! height(nz)=height(nz-1)+height(nz-1)-height(nz-2)
+ ! End doubled vertical resolution
+
+
+ ! Determine highest levels that can be within PBL
+ !************************************************
+
+ do kz=1,nz
+ if (height(kz).gt.hmixmax) then
+ nmixz=kz
+ goto 9
+ endif
+ end do
+9 continue
+
+ ! Do not repeat initialization of the Cartesian z grid
+ !*****************************************************
+
+ init=.false.
+
+ endif
+
+
+ ! Loop over the whole grid
+ !*************************
+
+ do jy=0,nymin1
+ do ix=0,nxmin1
+
+ ! NCEP version: find first level above ground
+ llev = 0
+ do i=1,nuvz
+ if (ps(ix,jy,1,n).lt.akz(i)) llev=i
+ end do
+ llev = llev+1
+ if (llev.gt.nuvz-2) llev = nuvz-2
+ ! if (llev.eq.nuvz-2) write(*,*) 'verttransform
+ ! +WARNING: LLEV eq NUZV-2'
+ ! NCEP version
+
+
+ ! compute height of pressure levels above ground
+ !***********************************************
+
+ tvold=tth(ix,jy,llev,n)*(1.+0.608*qvh(ix,jy,llev,n))
+ pold=akz(llev)
+ uvzlev(llev)=0.
+ wzlev(llev)=0.
+ uvwzlev(ix,jy,llev)=0.
+ rhoh(llev)=pold/(r_air*tvold)
+
+ do kz=llev+1,nuvz
+ pint=akz(kz)+bkz(kz)*ps(ix,jy,1,n)
+ tv=tth(ix,jy,kz,n)*(1.+0.608*qvh(ix,jy,kz,n))
+ rhoh(kz)=pint/(r_air*tv)
+
+ if (abs(tv-tvold).gt.0.2) then
+ uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)* &
+ (tv-tvold)/log(tv/tvold)
+ else
+ uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)*tv
+ endif
+ wzlev(kz)=uvzlev(kz)
+ uvwzlev(ix,jy,kz)=uvzlev(kz)
+
+ tvold=tv
+ pold=pint
+ end do
+
+
+ ! Switch on following lines to use doubled vertical resolution
+ ! Switch off the three lines above.
+ !*************************************************************
+ !22 uvwzlev(ix,jy,(kz-1)*2)=uvzlev(kz)
+ ! do 23 kz=2,nwz
+ !23 uvwzlev(ix,jy,(kz-1)*2+1)=wzlev(kz)
+ ! End doubled vertical resolution
+
+ ! pinmconv=(h2-h1)/(p2-p1)
+
+ pinmconv(llev)=(uvwzlev(ix,jy,llev+1)-uvwzlev(ix,jy,llev))/ &
+ ((aknew(llev+1)+bknew(llev+1)*ps(ix,jy,1,n))- &
+ (aknew(llev)+bknew(llev)*ps(ix,jy,1,n)))
+ do kz=llev+1,nz-1
+ pinmconv(kz)=(uvwzlev(ix,jy,kz+1)-uvwzlev(ix,jy,kz-1))/ &
+ ((aknew(kz+1)+bknew(kz+1)*ps(ix,jy,1,n))- &
+ (aknew(kz-1)+bknew(kz-1)*ps(ix,jy,1,n)))
+ end do
+ pinmconv(nz)=(uvwzlev(ix,jy,nz)-uvwzlev(ix,jy,nz-1))/ &
+ ((aknew(nz)+bknew(nz)*ps(ix,jy,1,n))- &
+ (aknew(nz-1)+bknew(nz-1)*ps(ix,jy,1,n)))
+
+
+ ! Levels, where u,v,t and q are given
+ !************************************
+
+ uu(ix,jy,1,n)=uuh(ix,jy,llev)
+ vv(ix,jy,1,n)=vvh(ix,jy,llev)
+ tt(ix,jy,1,n)=tth(ix,jy,llev,n)
+ qv(ix,jy,1,n)=qvh(ix,jy,llev,n)
+ pv(ix,jy,1,n)=pvh(ix,jy,llev)
+ rho(ix,jy,1,n)=rhoh(llev)
+ pplev(ix,jy,1,n)=akz(llev)
+ uu(ix,jy,nz,n)=uuh(ix,jy,nuvz)
+ vv(ix,jy,nz,n)=vvh(ix,jy,nuvz)
+ tt(ix,jy,nz,n)=tth(ix,jy,nuvz,n)
+ qv(ix,jy,nz,n)=qvh(ix,jy,nuvz,n)
+ pv(ix,jy,nz,n)=pvh(ix,jy,nuvz)
+ rho(ix,jy,nz,n)=rhoh(nuvz)
+ pplev(ix,jy,nz,n)=akz(nuvz)
+ kmin=llev+1
+ do iz=2,nz-1
+ do kz=kmin,nuvz
+ if(height(iz).gt.uvzlev(nuvz)) then
+ uu(ix,jy,iz,n)=uu(ix,jy,nz,n)
+ vv(ix,jy,iz,n)=vv(ix,jy,nz,n)
+ tt(ix,jy,iz,n)=tt(ix,jy,nz,n)
+ qv(ix,jy,iz,n)=qv(ix,jy,nz,n)
+ pv(ix,jy,iz,n)=pv(ix,jy,nz,n)
+ rho(ix,jy,iz,n)=rho(ix,jy,nz,n)
+ pplev(ix,jy,iz,n)=pplev(ix,jy,nz,n)
+ goto 30
+ endif
+ if ((height(iz).gt.uvzlev(kz-1)).and. &
+ (height(iz).le.uvzlev(kz))) then
+ dz1=height(iz)-uvzlev(kz-1)
+ dz2=uvzlev(kz)-height(iz)
+ dz=dz1+dz2
+ uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)*dz2+uuh(ix,jy,kz)*dz1)/dz
+ vv(ix,jy,iz,n)=(vvh(ix,jy,kz-1)*dz2+vvh(ix,jy,kz)*dz1)/dz
+ tt(ix,jy,iz,n)=(tth(ix,jy,kz-1,n)*dz2 &
+ +tth(ix,jy,kz,n)*dz1)/dz
+ qv(ix,jy,iz,n)=(qvh(ix,jy,kz-1,n)*dz2 &
+ +qvh(ix,jy,kz,n)*dz1)/dz
+ pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)*dz2+pvh(ix,jy,kz)*dz1)/dz
+ rho(ix,jy,iz,n)=(rhoh(kz-1)*dz2+rhoh(kz)*dz1)/dz
+ pplev(ix,jy,iz,n)=(akz(kz-1)*dz2+akz(kz)*dz1)/dz
+ endif
+ end do
+30 continue
+ end do
+
+
+ ! Levels, where w is given
+ !*************************
+
+ ww(ix,jy,1,n)=wwh(ix,jy,llev)*pinmconv(llev)
+ ww(ix,jy,nz,n)=wwh(ix,jy,nwz)*pinmconv(nz)
+ kmin=llev+1
+ do iz=2,nz
+ do kz=kmin,nwz
+ if ((height(iz).gt.wzlev(kz-1)).and. &
+ (height(iz).le.wzlev(kz))) then
+ dz1=height(iz)-wzlev(kz-1)
+ dz2=wzlev(kz)-height(iz)
+ dz=dz1+dz2
+ ww(ix,jy,iz,n)=(wwh(ix,jy,kz-1)*pinmconv(kz-1)*dz2 &
+ +wwh(ix,jy,kz)*pinmconv(kz)*dz1)/dz
+
+ endif
+ end do
+ end do
+
+
+ ! Compute density gradients at intermediate levels
+ !*************************************************
+
+ drhodz(ix,jy,1,n)=(rho(ix,jy,2,n)-rho(ix,jy,1,n))/ &
+ (height(2)-height(1))
+ do kz=2,nz-1
+ drhodz(ix,jy,kz,n)=(rho(ix,jy,kz+1,n)-rho(ix,jy,kz-1,n))/ &
+ (height(kz+1)-height(kz-1))
+ end do
+ drhodz(ix,jy,nz,n)=drhodz(ix,jy,nz-1,n)
+
+ end do
+ end do
+
+
+ !****************************************************************
+ ! Compute slope of eta levels in windward direction and resulting
+ ! vertical wind correction
+ !****************************************************************
+
+ do jy=1,ny-2
+ do ix=1,nx-2
+
+ ! NCEP version: find first level above ground
+ llev = 0
+ do i=1,nuvz
+ if (ps(ix,jy,1,n).lt.akz(i)) llev=i
+ end do
+ llev = llev+1
+ if (llev.gt.nuvz-2) llev = nuvz-2
+ ! if (llev.eq.nuvz-2) write(*,*) 'verttransform
+ ! +WARNING: LLEV eq NUZV-2'
+ ! NCEP version
+
+ kmin=llev+1
+ do iz=2,nz-1
+
+ ui=uu(ix,jy,iz,n)*dxconst/cos((real(jy)*dy+ylat0)*pi180)
+ vi=vv(ix,jy,iz,n)*dyconst
+
+ do kz=kmin,nz
+ if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
+ (height(iz).le.uvwzlev(ix,jy,kz))) then
+ dz1=height(iz)-uvwzlev(ix,jy,kz-1)
+ dz2=uvwzlev(ix,jy,kz)-height(iz)
+ dz=dz1+dz2
+ kl=kz-1
+ klp=kz
+ goto 47
+ endif
+ end do
+
+47 ix1=ix-1
+ jy1=jy-1
+ ixp=ix+1
+ jyp=jy+1
+
+ dzdx1=(uvwzlev(ixp,jy,kl)-uvwzlev(ix1,jy,kl))/2.
+ dzdx2=(uvwzlev(ixp,jy,klp)-uvwzlev(ix1,jy,klp))/2.
+ dzdx=(dzdx1*dz2+dzdx2*dz1)/dz
+
+ dzdy1=(uvwzlev(ix,jyp,kl)-uvwzlev(ix,jy1,kl))/2.
+ dzdy2=(uvwzlev(ix,jyp,klp)-uvwzlev(ix,jy1,klp))/2.
+ dzdy=(dzdy1*dz2+dzdy2*dz1)/dz
+
+ ww(ix,jy,iz,n)=ww(ix,jy,iz,n)+(dzdx*ui+dzdy*vi)
+
+ end do
+
+ end do
+ end do
+
+
+ ! If north pole is in the domain, calculate wind velocities in polar
+ ! stereographic coordinates
+ !*******************************************************************
+
+ if (nglobal) then
+ do jy=int(switchnorthg)-2,nymin1
+ ylat=ylat0+real(jy)*dy
+ do ix=0,nxmin1
+ xlon=xlon0+real(ix)*dx
+ do iz=1,nz
+ call cc2gll(northpolemap,ylat,xlon,uu(ix,jy,iz,n), &
+ vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
+ vvpol(ix,jy,iz,n))
+ end do
+ end do
+ end do
+
+
+ do iz=1,nz
+
+ ! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
+ xlon=xlon0+real(nx/2-1)*dx
+ xlonr=xlon*pi/180.
+ ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+ &
+ vv(nx/2-1,nymin1,iz,n)**2)
+ if(vv(nx/2-1,nymin1,iz,n).lt.0.) then
+ ddpol=atan(uu(nx/2-1,nymin1,iz,n)/ &
+ vv(nx/2-1,nymin1,iz,n))-xlonr
+ elseif (vv(nx/2-1,nymin1,iz,n).gt.0.) then
+ ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/ &
+ vv(nx/2-1,nymin1,iz,n))-xlonr
+ else
+ ddpol=pi/2-xlonr
+ endif
+ if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
+ if(ddpol.gt.2.0*pi) ddpol=ddpol-2.0*pi
+
+ ! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
+ xlon=180.0
+ xlonr=xlon*pi/180.
+ ylat=90.0
+ uuaux=-ffpol*sin(xlonr+ddpol)
+ vvaux=-ffpol*cos(xlonr+ddpol)
+ call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
+ vvpolaux)
+
+ jy=nymin1
+ do ix=0,nxmin1
+ uupol(ix,jy,iz,n)=uupolaux
+ vvpol(ix,jy,iz,n)=vvpolaux
+ end do
+ end do
+
+
+ ! Fix: Set W at pole to the zonally averaged W of the next equator-
+ ! ward parallel of latitude
+
+ do iz=1,nz
+ wdummy=0.
+ jy=ny-2
+ do ix=0,nxmin1
+ wdummy=wdummy+ww(ix,jy,iz,n)
+ end do
+ wdummy=wdummy/real(nx)
+ jy=nymin1
+ do ix=0,nxmin1
+ ww(ix,jy,iz,n)=wdummy
+ end do
+ end do
+
+ endif
+
+
+ ! If south pole is in the domain, calculate wind velocities in polar
+ ! stereographic coordinates
+ !*******************************************************************
+
+ if (sglobal) then
+ do jy=0,int(switchsouthg)+3
+ ylat=ylat0+real(jy)*dy
+ do ix=0,nxmin1
+ xlon=xlon0+real(ix)*dx
+ do iz=1,nz
+ call cc2gll(southpolemap,ylat,xlon,uu(ix,jy,iz,n), &
+ vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
+ vvpol(ix,jy,iz,n))
+ end do
+ end do
+ end do
+
+ do iz=1,nz
+
+ ! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
+ xlon=xlon0+real(nx/2-1)*dx
+ xlonr=xlon*pi/180.
+ ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+ &
+ vv(nx/2-1,0,iz,n)**2)
+ if(vv(nx/2-1,0,iz,n).lt.0.) then
+ ddpol=atan(uu(nx/2-1,0,iz,n)/ &
+ vv(nx/2-1,0,iz,n))+xlonr
+ elseif (vv(nx/2-1,0,iz,n).gt.0.) then
+ ddpol=pi+atan(uu(nx/2-1,0,iz,n)/ &
+ vv(nx/2-1,0,iz,n))-xlonr
+ else
+ ddpol=pi/2-xlonr
+ endif
+ if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
+ if(ddpol.gt.2.0*pi) ddpol=ddpol-2.0*pi
+
+ ! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
+ xlon=180.0
+ xlonr=xlon*pi/180.
+ ylat=-90.0
+ uuaux=+ffpol*sin(xlonr-ddpol)
+ vvaux=-ffpol*cos(xlonr-ddpol)
+ call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
+ vvpolaux)
+
+ jy=0
+ do ix=0,nxmin1
+ uupol(ix,jy,iz,n)=uupolaux
+ vvpol(ix,jy,iz,n)=vvpolaux
+ end do
+ end do
+
+
+ ! Fix: Set W at pole to the zonally averaged W of the next equator-
+ ! ward parallel of latitude
+
+ do iz=1,nz
+ wdummy=0.
+ jy=1
+ do ix=0,nxmin1
+ wdummy=wdummy+ww(ix,jy,iz,n)
+ end do
+ wdummy=wdummy/real(nx)
+ jy=0
+ do ix=0,nxmin1
+ ww(ix,jy,iz,n)=wdummy
+ end do
+ end do
+ endif
+
+
+ ! write (*,*) 'initializing clouds, n:',n,nymin1,nxmin1,nz
+ ! create a cloud and rainout/washout field, clouds occur where rh>80%
+ ! total cloudheight is stored at level 0
+ do jy=0,nymin1
+ do ix=0,nxmin1
+ rain_cloud_above=0
+ lsp=lsprec(ix,jy,1,n)
+ convp=convprec(ix,jy,1,n)
+ cloudsh(ix,jy,n)=0
+ do kz_inv=1,nz-1
+ kz=nz-kz_inv+1
+ pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
+ rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
+ clouds(ix,jy,kz,n)=0
+ if (rh.gt.0.8) then ! in cloud
+ if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
+ rain_cloud_above=1
+ cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
+ height(kz)-height(kz-1)
+ if (lsp.ge.convp) then
+ clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
+ else
+ clouds(ix,jy,kz,n)=2 ! convp dominated rainout
+ endif
+ else ! no precipitation
+ clouds(ix,jy,kz,n)=1 ! cloud
+ endif
+ else ! no cloud
+ if (rain_cloud_above.eq.1) then ! scavenging
+ if (lsp.ge.convp) then
+ clouds(ix,jy,kz,n)=5 ! lsp dominated washout
+ else
+ clouds(ix,jy,kz,n)=4 ! convp dominated washout
+ endif
+ endif
+ endif
+ end do
+ end do
+ end do
+
+
+end subroutine verttransform_grib2nc4_gfs
--
GitLab