[Dart-dev] [4305] DART/trunk/models/rose/model_mod.f90: Incorporated the model_interpolate routine (from lightning) and
nancy at ucar.edu
nancy at ucar.edu
Mon Mar 8 17:35:17 MST 2010
Revision: 4305
Author: thoar
Date: 2010-03-08 17:35:17 -0700 (Mon, 08 Mar 2010)
Log Message:
-----------
Incorporated the model_interpolate routine (from lightning) and
removed all the unused variables. Now uses the netCDF F90 interface
exclusively. Need to check the get_state_meta_data() routine ...
for var_type information.
Modified Paths:
--------------
DART/trunk/models/rose/model_mod.f90
-------------- next part --------------
Modified: DART/trunk/models/rose/model_mod.f90
===================================================================
--- DART/trunk/models/rose/model_mod.f90 2010-03-07 17:31:26 UTC (rev 4304)
+++ DART/trunk/models/rose/model_mod.f90 2010-03-09 00:35:17 UTC (rev 4305)
@@ -21,13 +21,17 @@
use time_manager_mod, only : time_type,set_time,print_time
use location_mod, only : location_type, get_close_maxdist_init, &
get_close_obs_init, get_close_obs, &
- set_location, get_location, &
- query_location, get_dist
+ set_location, get_location, query_location, &
+ get_dist, vert_is_height
use utilities_mod, only : file_exist, open_file, close_file, &
error_handler, E_ERR, E_MSG, E_WARN, nmlfileunit, &
do_output, find_namelist_in_file, check_namelist_read, &
- do_nml_file, do_nml_term
+ do_nml_file, do_nml_term, nc_check
use random_seq_mod, only : random_seq_type, init_random_seq, random_gaussian
+use obs_kind_mod, only : KIND_U_WIND_COMPONENT, KIND_V_WIND_COMPONENT,&
+ KIND_TEMPERATURE
+use typesizes
+use netcdf
! ROSE Modules
use params, only : nx, ny, nz, ntime
@@ -60,6 +64,7 @@
vector_to_prog_var, &
read_ROSE_restart, &
update_ROSE_restart, &
+ read_ROSE_tref,&
init_model_instance, &
end_model_instance
@@ -69,7 +74,19 @@
real(r8), pointer :: vars_3d(:,:,:,:)
end type model_type
+integer, parameter :: TYPE_local_U0 = 0, &
+ TYPE_local_V0 = 1, &
+ TYPE_local_T0 = 2, &
+ TYPE_local_U = 3, &
+ TYPE_local_V = 4, &
+ TYPE_local_T = 5
+
!----------------------------------------------------------------------
+
+! ROSE vertical grid variables
+real(r8) :: dz = 2100.0_r8, zbot = 16800.0_r8
+real, dimension (nz) :: tref !! auxiliary variable needed for H
+
! Global storage for describing ROSE model class
integer :: model_size
type(time_type) :: Time_step_ROSE
@@ -113,7 +130,8 @@
revision = "$Revision$", &
revdate = "$Date$"
-character(len = 129) :: errstring
+character(len = 129) :: msgstring
+logical, save :: module_initialized = .false.
contains
@@ -128,15 +146,13 @@
! In models that require pre-computed static data, for instance
! spherical harmonic weights, these would also be computed here.
! Can be a NULL INTERFACE for the simplest models.
-!
-integer :: iunit, ierr, io, i, j, k
+integer :: iunit, io, i, j, k
integer :: Time_step_seconds, Time_step_days
integer :: seconds_of_day = 86400
real(r8) :: d_lat, d_lon
real(r8) :: z_m
real(r8) :: dz = 2100.0_r8, zbot = 16800.0_r8
-character(len=129) :: nml_string
! Read the namelist rose_nml from the file rose.nml
call find_namelist_in_file("rose.nml", "rose_nml", iunit)
@@ -146,7 +162,6 @@
if (do_nml_file()) write(nmlfileunit, nml=rose_nml)
if (do_nml_term()) write( * , nml=rose_nml)
-
! Read the namelist entry for model_mod from file input.nml
call find_namelist_in_file("input.nml", "model_nml", iunit)
read(iunit, nml = model_nml, iostat = io)
@@ -192,6 +207,12 @@
levs(k) = 1013._r8*exp(-z_m/7.e3) ![hPa]
enddo
+! Read in rose_restart file for tref (global mean temperature)
+! for computation of total temperature needed for H
+call read_ROSE_tref(restart_file)
+
+module_initialized = .true.
+
end subroutine static_init_model
@@ -209,6 +230,7 @@
real(r8), intent(out) :: x(:)
+if ( .not. module_initialized ) call static_init_model
end subroutine init_conditions
@@ -233,6 +255,8 @@
real(r8), intent(inout) :: x(:)
type(time_type), intent(in) :: time
+if ( .not. module_initialized ) call static_init_model
+
end subroutine adv_1step
@@ -245,6 +269,8 @@
integer :: get_model_size
+if ( .not. module_initialized ) call static_init_model
+
get_model_size = model_size
end function get_model_size
@@ -264,6 +290,10 @@
type(time_type), intent(out) :: time
+! For now, just set to 0
+time = set_time(0, 0)
+if ( .not. module_initialized ) call static_init_model
+
end subroutine init_time
@@ -289,13 +319,190 @@
real(r8), intent(out) :: obs_val
integer, intent(out) :: istatus
+integer :: i, vstatus, which_vert
+integer :: lat_below, lat_above, lon_below, lon_above
+real(r8) :: lon_fract, temp_lon, lat_fract
+real(r8) :: lon, lat, height, lon_lat_lev(3)
+real(r8) :: bot_lon, top_lon, delta_lon, bot_lat, top_lat
+real(r8) :: val(2,2), a(2)
+if ( .not. module_initialized ) call static_init_model
+
! Default for successful return
istatus = 0
+vstatus = 0
+! Get the position, determine if it is model level or pressure in vertical
+lon_lat_lev = get_location(location)
+lon = lon_lat_lev(1); lat = lon_lat_lev(2);
+if(vert_is_height(location)) then
+ height = lon_lat_lev(3)
+else
+ which_vert = nint(query_location(location))
+ write(msgstring,*) 'vertical coordinate type:',which_vert,' cannot be handled'
+ call error_handler(E_ERR,'model_interpolate',msgstring,source,revision,revdate)
+endif
+
+! Get lon and lat grid specs, nx, ny, nz from ROSE param mod
+ bot_lon = lons(1)
+ top_lon = lons(nx)
+ delta_lon = lons(2) - lons(1)
+ bot_lat = lats(1)
+ top_lat = lats(ny)
+
+! Compute bracketing lon indices
+if(lon >= bot_lon .and. lon <= top_lon) then
+ lon_below = int((lon - bot_lon) / delta_lon) + 1
+ lon_above = lon_below + 1
+ lon_fract = (lon - ((lon_below - 1) * delta_lon + bot_lon)) / delta_lon
+else
+! At wraparound point
+ lon_below = nx
+ lon_above = 1
+ if(lon < bot_lon) then
+ temp_lon = lon + 360.0
+ else
+ temp_lon = lon
+ endif
+ lon_fract = (temp_lon - top_lon) / delta_lon
+endif
+
+! Next, compute neighboring lat rows
+! NEED TO BE VERY CAREFUL ABOUT POLES; WHAT'S BEING DONE MAY BE WRONG
+! Inefficient search used for latitudes in Gaussian grid. Might want to speed up.
+if(lat >= bot_lat .and. lat <= top_lat) then
+
+ do i = 2, ny
+ if(lat <= lats(i)) then
+ lat_above = i
+ lat_below = i - 1
+ lat_fract = (lat - lats(lat_below)) / (lats(lat_above) - lats(lat_below))
+ goto 20
+ end if
+ end do
+else if(lat <= bot_lat) then
+! South of bottom lat NEED TO DO BETTER: NOT REALLY REGULAR
+ lat_below = 1
+ lat_above = 1
+ lat_fract = 1.
+else
+! North of top lat NEED TO DO BETTER: NOT REALLY REGULAR
+ lat_below = ny
+ lat_above = ny
+ lat_fract = 1.
+endif
+
+20 continue
+
+! Now, need to find the values for the four corners
+ call get_val(val(1, 1), x, lon_below, lat_below, height, itype, vstatus)
+ if (vstatus /= 1) call get_val(val(1, 2), x, lon_below, lat_above, height, itype, vstatus)
+ if (vstatus /= 1) call get_val(val(2, 1), x, lon_above, lat_below, height, itype, vstatus)
+ if (vstatus /= 1) call get_val(val(2, 2), x, lon_above, lat_above, height, itype, vstatus)
+
+
+! kdr Guam;s
+! istatus meaning return expected obs? assimilate?
+! 0 obs and model are fine; yes yes
+! 1 fatal problem; no no
+! 2 exclude valid obs yes no
+!
+istatus = vstatus
+if(istatus /= 1) then
+ do i = 1, 2
+ a(i) = lon_fract * val(2, i) + (1.0 - lon_fract) * val(1, i)
+ end do
+ obs_val = lat_fract * a(2) + (1.0 - lat_fract) * a(1)
+else
+ obs_val = 0.
+endif
+
+!print*, 'model_interpolate', lon, lat, height,obs_val
+
end subroutine model_interpolate
+subroutine get_val(val, x, lon_index, lat_index, height, obs_kind, istatus)
+!------------------------------------------------------------------
+!
+real(r8), intent(out) :: val
+real(r8), intent(in) :: x(:), height
+integer, intent(in) :: lon_index, lat_index, obs_kind
+integer, intent(out) :: istatus
+
+integer :: var_type
+integer :: k, lev_top, lev_bottom
+real(r8) :: zgrid
+real(r8) :: val_top, val_bottom, frac_lev
+
+! No errors to start with
+istatus = 0
+
+! Pressure levels ! from msetfix.f
+! 1013 [mb] exp(-z/7km) with a fixed scale height (7km)
+! (2.5 km increment from 17.5 km)
+!levs(k) = 1013._r8*exp(-z_m/7.e3) ![hPa]
+
+height_loop:do k = 1, nz
+ zgrid = (k-1)*dz + zbot
+ if (height <= zgrid) then
+ lev_top = k
+ frac_lev = (zgrid - height)/dz
+ exit height_loop
+ endif
+enddo height_loop
+
+lev_bottom = lev_top -1
+
+if (obs_kind == KIND_TEMPERATURE) then
+
+ var_type = TYPE_local_T
+! val_top = x(get_index(lat_index, lon_index, lev_top, var_type))
+! val_bottom = x(get_index(lat_index, lon_index, lev_bottom, var_type))
+
+ val_top = x(get_index(lat_index,lon_index,lev_top,var_type)) &
+ + tref(lev_top)
+ val_bottom = x(get_index(lat_index,lon_index,lev_bottom,var_type)) &
+ + tref(lev_bottom)
+
+elseif (obs_kind == KIND_U_WIND_COMPONENT) then
+
+ var_type = TYPE_local_U
+ val_top = x(get_index(lat_index, lon_index, lev_top, var_type))
+ val_bottom = x(get_index(lat_index, lon_index, lev_bottom, var_type))
+
+elseif (obs_kind == KIND_V_WIND_COMPONENT) then
+
+ var_type = TYPE_local_V
+ val_top = x(get_index(lat_index, lon_index, lev_top, var_type))
+ val_bottom = x(get_index(lat_index, lon_index, lev_bottom, var_type))
+
+else
+ istatus = 1
+ val = 0.
+ return
+end if
+
+val = frac_lev * val_bottom + (1.0 - frac_lev) * val_top
+
+end subroutine get_val
+
+
+
+function get_index(lat_index, lon_index, lev_index, var_type)
+!------------------------------------------------------------------
+!
+integer, intent(in) :: lat_index, lon_index, lev_index, var_type
+integer :: get_index
+
+get_index = 1 + var_type + (lev_index - 1)*state_num_3d &
+ + (lat_index - 1)*state_num_3d*nz &
+ + (lon_index - 1)*state_num_3d*nz*ny
+
+end function get_index
+
+
+
function get_model_time_step()
!------------------------------------------------------------------
!
@@ -307,6 +514,8 @@
type(time_type) :: get_model_time_step
+if ( .not. module_initialized ) call static_init_model
+
get_model_time_step = Time_step_ROSE
end function get_model_time_step
@@ -328,10 +537,12 @@
integer, intent(out), optional :: var_type
integer :: indx, num_per_col, col_num, col_elem, lon_index,&
- & lat_index, lev_index
+ lat_index, lev_index
real(r8) :: lon, lat, lev
integer :: local_var_type, var_type_temp
+if ( .not. module_initialized ) call static_init_model
+
! Easier to compute with a 0 to size - 1 index
indx = index_in - 1
@@ -357,26 +568,26 @@
! Find which var_type this element is
var_type_temp = mod(col_elem, state_num_3d )
if(var_type_temp == 0) then
- local_var_type = TYPE_U0
+ local_var_type = KIND_U_WIND_COMPONENT
else if(var_type_temp == 1) then
- local_var_type = TYPE_V0
+ local_var_type = KIND_V_WIND_COMPONENT
else if(var_type_temp == 2) then
- local_var_type = TYPE_T0
+ local_var_type = KIND_TEMPERATURE
else if(var_type_temp == 3) then
- local_var_type = TYPE_U
+ local_var_type = KIND_U_WIND_COMPONENT
else if(var_type_temp == 4) then
- local_var_type = TYPE_V
+ local_var_type = KIND_V_WIND_COMPONENT
else if(var_type_temp == 5) then
- local_var_type = TYPE_T
-else if(var_type_temp == 6) then
- local_var_type = TYPE_Q_H
-else if(var_type_temp == 7) then
- local_var_type = TYPE_Q_OH
-else
- local_var_type = TYPE_Q_O
+ local_var_type = KIND_TEMPERATURE
+! else if(var_type_temp == 6) then
+! local_var_type = TYPE_Q_H
+!else if(var_type_temp == 7) then
+! local_var_type = TYPE_Q_OH
+!else
+! local_var_type = TYPE_Q_O
endif
-write(*, '(1x,3(f6.2,1x),i3)') lon, lat, lev, local_var_type
+!write(*, '(1x,3(f6.2,1x),i3)') lon, lat, lev, local_var_type
location = set_location(lon, lat, lev, 2) ! 2 == pressure (hPa)
@@ -393,6 +604,8 @@
! Does any shutdown and clean-up needed for model. Can be a NULL
! INTERFACE if the model has no need to clean up storgae, etc.
+if ( .not. module_initialized ) call static_init_model
+
end subroutine end_model
@@ -425,21 +638,17 @@
! NF90_put_var ! provide values for variable
! NF90_CLOSE ! close: save updated netCDF dataset
-use typeSizes
-use netcdf
-
integer, intent(in) :: ncFileID ! netCDF file identifier
integer :: ierr ! return value of function
integer :: nDimensions, nVariables, nAttributes, unlimitedDimID
-integer :: StateVarDimID
integer :: MemberDimID, TimeDimID
integer :: lonDimID, latDimID, levDimID
integer :: lonVarID, latVarID, levVarID
integer :: u1VarID, v1VarID, t1VarID, uVarID, vVarID, tVarID
-integer :: qnHVarID, qnOHVarID, qnOVarID
+! integer :: qnHVarID, qnOHVarID, qnOVarID
character(len=8) :: crdate ! needed by F90 DATE_AND_TIME intrinsic
character(len=10) :: crtime ! needed by F90 DATE_AND_TIME intrinsic
@@ -449,25 +658,27 @@
!------------------------------------------------------------------
+if ( .not. module_initialized ) call static_init_model
+
ierr = 0 ! assume normal termination
!--------------------------------------------------------------------
! make sure ncFileID refers to an open netCDF file
!--------------------------------------------------------------------
-call check(nf90_Inquire(ncFileID, nDimensions, nVariables, nAttributes, unlimitedDimID))
-call check(nf90_Redef(ncFileID))
+call nc_check(nf90_Inquire(ncFileID, nDimensions, nVariables, nAttributes, unlimitedDimID),'nc_write_model_atts','inquire')
+call nc_check(nf90_Redef(ncFileID),'nc_write_model_atts','redef')
!-------------------------------------------------------------------------------
! We need the dimension ID for the number of copies
!-------------------------------------------------------------------------------
-call check(nf90_inq_dimid(ncid=ncFileID, name="copy", dimid=MemberDimID))
-call check(nf90_inq_dimid(ncid=ncFileID, name="time", dimid= TimeDimID))
+call nc_check(nf90_inq_dimid(ncid=ncFileID, name="copy", dimid=MemberDimID),'nc_write_model_atts')
+call nc_check(nf90_inq_dimid(ncid=ncFileID, name="time", dimid= TimeDimID),'nc_write_model_atts')
if ( TimeDimID /= unlimitedDimId ) then
- write(errstring,*)'Time Dimension ID ',TimeDimID,' should equal Unlimited Dimension ID',unlimitedDimID
- call error_handler(E_ERR,'nc_write_model_atts', errstring, source, revision, revdate)
+ write(msgstring,*)'Time Dimension ID ',TimeDimID,' should equal Unlimited Dimension ID',unlimitedDimID
+ call error_handler(E_ERR,'nc_write_model_atts', msgstring, source, revision, revdate)
endif
!-------------------------------------------------------------------------------
@@ -478,44 +689,44 @@
write(str1,'(''YYYY MM DD HH MM SS = '',i4,5(1x,i2.2))') &
values(1), values(2), values(3), values(5), values(6), values(7)
-call check(nf90_put_att(ncFileID, NF90_GLOBAL, "creation_date",str1))
-call check(nf90_put_att(ncFileID, NF90_GLOBAL, "model_source",source))
-call check(nf90_put_att(ncFileID, NF90_GLOBAL, "model_revision",revision))
-call check(nf90_put_att(ncFileID, NF90_GLOBAL, "model_revdate",revdate))
-call check(nf90_put_att(ncFileID, NF90_GLOBAL, "model","ROSE"))
+call nc_check(nf90_put_att(ncFileID, NF90_GLOBAL, "creation_date",str1),'nc_write_model_atts')
+call nc_check(nf90_put_att(ncFileID, NF90_GLOBAL, "model_source",source),'nc_write_model_atts')
+call nc_check(nf90_put_att(ncFileID, NF90_GLOBAL, "model_revision",revision),'nc_write_model_atts')
+call nc_check(nf90_put_att(ncFileID, NF90_GLOBAL, "model_revdate",revdate),'nc_write_model_atts')
+call nc_check(nf90_put_att(ncFileID, NF90_GLOBAL, "model","ROSE"),'nc_write_model_atts')
!-------------------------------------------------------------------------------
! Define the dimensions IDs
!-------------------------------------------------------------------------------
-call check(nf90_def_dim(ncid=ncFileID, name="lon", len = nx, dimid = lonDimID))
-call check(nf90_def_dim(ncid=ncFileID, name="lat", len = ny, dimid = latDimID))
-call check(nf90_def_dim(ncid=ncFileID, name="lev", len = nz, dimid = levDimID))
+call nc_check(nf90_def_dim(ncid=ncFileID, name="lon", len = nx, dimid = lonDimID),'nc_write_model_atts')
+call nc_check(nf90_def_dim(ncid=ncFileID, name="lat", len = ny, dimid = latDimID),'nc_write_model_atts')
+call nc_check(nf90_def_dim(ncid=ncFileID, name="lev", len = nz, dimid = levDimID),'nc_write_model_atts')
!-------------------------------------------------------------------------------
! Create the (empty) Variables and the Attributes
!-------------------------------------------------------------------------------
-call check(nf90_def_var(ncFileID, name="lon", xtype=nf90_double, &
- dimids=lonDimID, varid=lonVarID) )
-call check(nf90_put_att(ncFileID, lonVarID, "long_name", "longitude"))
-call check(nf90_put_att(ncFileID, lonVarID, "cartesian_axis", "X"))
-call check(nf90_put_att(ncFileID, lonVarID, "units", "degrees_east")) ! check
-call check(nf90_put_att(ncFileID, lonVarID, "valid_range", (/ 0.0_r8, 360.0_r8 /)))
+call nc_check(nf90_def_var(ncFileID, name="lon", xtype=nf90_double, &
+ dimids=lonDimID, varid=lonVarID),'nc_write_model_atts' )
+call nc_check(nf90_put_att(ncFileID, lonVarID, "long_name", "longitude"),'nc_write_model_atts')
+call nc_check(nf90_put_att(ncFileID, lonVarID, "cartesian_axis", "X"),'nc_write_model_atts')
+call nc_check(nf90_put_att(ncFileID, lonVarID, "units", "degrees_east"),'nc_write_model_atts')
+call nc_check(nf90_put_att(ncFileID, lonVarID, "valid_range", (/ 0.0_r8, 360.0_r8 /)),'nc_write_model_atts')
-call check(nf90_def_var(ncFileID, name="lat", xtype=nf90_double, &
- dimids=latDimID, varid=latVarID) )
-call check(nf90_put_att(ncFileID, latVarID, "long_name", "latitude"))
-call check(nf90_put_att(ncFileID, latVarID, "cartesian_axis", "Y"))
-call check(nf90_put_att(ncFileID, latVarID, "units", "degrees_north")) !check
-call check(nf90_put_att(ncFileID, latVarID, "valid_range", (/ -90.0_r8, 90.0_r8 /)))
+call nc_check(nf90_def_var(ncFileID, name="lat", xtype=nf90_double, &
+ dimids=latDimID, varid=latVarID),'nc_write_model_atts' )
+call nc_check(nf90_put_att(ncFileID, latVarID, "long_name", "latitude"),'nc_write_model_atts')
+call nc_check(nf90_put_att(ncFileID, latVarID, "cartesian_axis", "Y"),'nc_write_model_atts')
+call nc_check(nf90_put_att(ncFileID, latVarID, "units", "degrees_north"),'nc_write_model_atts')
+call nc_check(nf90_put_att(ncFileID, latVarID, "valid_range", (/ -90.0_r8, 90.0_r8 /)),'nc_write_model_atts')
-call check(nf90_def_var(ncFileID, name="lev", xtype=nf90_double, &
- dimids=levDimID, varid=levVarID) )
-call check(nf90_put_att(ncFileID, levVarID, "long_name", "level"))
-call check(nf90_put_att(ncFileID, levVarID, "cartesian_axis", "Z"))
-call check(nf90_put_att(ncFileID, levVarID, "units", "hPa"))
-call check(nf90_put_att(ncFileID, levVarID, "positive", "down")) !check
+call nc_check(nf90_def_var(ncFileID, name="lev", xtype=nf90_double, &
+ dimids=levDimID, varid=levVarID),'nc_write_model_atts')
+call nc_check(nf90_put_att(ncFileID, levVarID, "long_name", "level"),'nc_write_model_atts')
+call nc_check(nf90_put_att(ncFileID, levVarID, "cartesian_axis", "Z"),'nc_write_model_atts')
+call nc_check(nf90_put_att(ncFileID, levVarID, "units", "hPa"),'nc_write_model_atts')
+call nc_check(nf90_put_att(ncFileID, levVarID, "positive", "down"),'nc_write_model_atts')
!----------------------------------------------------------------------------
! Create attributes for the state vector
@@ -526,92 +737,78 @@
! which is consistent with the Fortran convention of having the
! fastest-moving index on the left.
-call check(nf90_def_var(ncid=ncFileID, name="u1", xtype=nf90_real, &
+call nc_check(nf90_def_var(ncid=ncFileID, name="u1", xtype=nf90_real, &
dimids = (/ levDimID, lonDimID, latDimID, MemberDimID, unlimitedDimID /), &
- varid = u1VarID))
-call check(nf90_put_att(ncFileID, u1VarID, "long_name", "zonal wind component")) !check
-call check(nf90_put_att(ncFileID, u1VarID, "units", "m/s"))
+ varid = u1VarID),'nc_write_model_atts')
+call nc_check(nf90_put_att(ncFileID, u1VarID, "long_name", "zonal wind component"),'nc_write_model_atts')
+call nc_check(nf90_put_att(ncFileID, u1VarID, "units", "m/s"),'nc_write_model_atts')
-call check(nf90_def_var(ncid=ncFileID, name="v1", xtype=nf90_real, &
+call nc_check(nf90_def_var(ncid=ncFileID, name="v1", xtype=nf90_real, &
dimids = (/ levDimID, lonDimID, latDimID, MemberDimID, unlimitedDimID /), &
- varid = v1VarID))
-call check(nf90_put_att(ncFileID, v1VarID, "long_name", "meridional wind component")) !check
-call check(nf90_put_att(ncFileID, v1VarID, "units", "m/s"))
+ varid = v1VarID),'nc_write_model_atts')
+call nc_check(nf90_put_att(ncFileID, v1VarID, "long_name", "meridional wind component"),'nc_write_model_atts')
+call nc_check(nf90_put_att(ncFileID, v1VarID, "units", "m/s"),'nc_write_model_atts')
-call check(nf90_def_var(ncid=ncFileID, name="t1", xtype=nf90_real, &
+call nc_check(nf90_def_var(ncid=ncFileID, name="t1", xtype=nf90_real, &
dimids = (/ levDimID, lonDimID, latDimID, MemberDimID, unlimitedDimID /), &
- varid = t1VarID))
-call check(nf90_put_att(ncFileID, t1VarID, "long_name", "temperature")) !check
-call check(nf90_put_att(ncFileID, t1VarID, "units", "degrees Kelvin"))
+ varid = t1VarID),'nc_write_model_atts')
+call nc_check(nf90_put_att(ncFileID, t1VarID, "long_name", "temperature"),'nc_write_model_atts')
+call nc_check(nf90_put_att(ncFileID, t1VarID, "units", "degrees Kelvin"),'nc_write_model_atts')
-call check(nf90_def_var(ncid=ncFileID, name="u", xtype=nf90_real, &
+call nc_check(nf90_def_var(ncid=ncFileID, name="u", xtype=nf90_real, &
dimids = (/ levDimID, lonDimID, latDimID, MemberDimID, unlimitedDimID /), &
- varid = uVarID))
-call check(nf90_put_att(ncFileID, uVarID, "long_name", "zonal wind component"))
-call check(nf90_put_att(ncFileID, uVarID, "units", "m/s"))
+ varid = uVarID),'nc_write_model_atts')
+call nc_check(nf90_put_att(ncFileID, uVarID, "long_name", "zonal wind component"),'nc_write_model_atts')
+call nc_check(nf90_put_att(ncFileID, uVarID, "units", "m/s"),'nc_write_model_atts')
-call check(nf90_def_var(ncid=ncFileID, name="v", xtype=nf90_real, &
+call nc_check(nf90_def_var(ncid=ncFileID, name="v", xtype=nf90_real, &
dimids = (/ levDimID, lonDimID, latDimID, MemberDimID, unlimitedDimID /), &
- varid = vVarID))
-call check(nf90_put_att(ncFileID, vVarID, "long_name", "meridional wind component"))
-call check(nf90_put_att(ncFileID, vVarID, "units", "m/s"))
+ varid = vVarID),'nc_write_model_atts')
+call nc_check(nf90_put_att(ncFileID, vVarID, "long_name", "meridional wind component"),'nc_write_model_atts')
+call nc_check(nf90_put_att(ncFileID, vVarID, "units", "m/s"),'nc_write_model_atts')
-call check(nf90_def_var(ncid=ncFileID, name="t", xtype=nf90_real, &
+call nc_check(nf90_def_var(ncid=ncFileID, name="t", xtype=nf90_real, &
dimids = (/ levDimID, lonDimID, latDimID, MemberDimID, unlimitedDimID /), &
- varid = tVarID))
-call check(nf90_put_att(ncFileID, tVarID, "long_name", "temperature"))
-call check(nf90_put_att(ncFileID, tVarID, "units", "degrees Kelvin"))
+ varid = tVarID),'nc_write_model_atts')
+call nc_check(nf90_put_att(ncFileID, tVarID, "long_name", "temperature"),'nc_write_model_atts')
+call nc_check(nf90_put_att(ncFileID, tVarID, "units", "degrees Kelvin"),'nc_write_model_atts')
-call check(nf90_def_var(ncid=ncFileID, name="qnH", xtype=nf90_real, &
- dimids = (/ levDimID, lonDimID, latDimID, MemberDimID, unlimitedDimID /), &
- varid = qnHVarID))
-call check(nf90_put_att(ncFileID, qnHVarID, "long_name", "mixing ratio H"))
-call check(nf90_put_att(ncFileID, qnHVarID, "units", "?"))
+! call nc_check(nf90_def_var(ncid=ncFileID, name="qnH", xtype=nf90_real, &
+! dimids = (/ levDimID, lonDimID, latDimID, MemberDimID, unlimitedDimID /), &
+! varid = qnHVarID),'nc_write_model_atts')
+! call nc_check(nf90_put_att(ncFileID, qnHVarID, "long_name", "mixing ratio H"),'nc_write_model_atts')
+! call nc_check(nf90_put_att(ncFileID, qnHVarID, "units", "?"),'nc_write_model_atts')
-call check(nf90_def_var(ncid=ncFileID, name="qnOH", xtype=nf90_real, &
- dimids = (/ levDimID, lonDimID, latDimID, MemberDimID, unlimitedDimID /), &
- varid = qnOHVarID))
-call check(nf90_put_att(ncFileID, qnOHVarID, "long_name", "mixing ratio OH"))
-call check(nf90_put_att(ncFileID, qnOHVarID, "units", "?"))
+! call nc_check(nf90_def_var(ncid=ncFileID, name="qnOH", xtype=nf90_real, &
+! dimids = (/ levDimID, lonDimID, latDimID, MemberDimID, unlimitedDimID /), &
+! varid = qnOHVarID),'nc_write_model_atts')
+! call nc_check(nf90_put_att(ncFileID, qnOHVarID, "long_name", "mixing ratio OH"),'nc_write_model_atts')
+! call nc_check(nf90_put_att(ncFileID, qnOHVarID, "units", "?"),'nc_write_model_atts')
-call check(nf90_def_var(ncid=ncFileID, name="qnO", xtype=nf90_real, &
- dimids = (/ levDimID, lonDimID, latDimID, MemberDimID, unlimitedDimID /), &
- varid = qnOVarID))
-call check(nf90_put_att(ncFileID, qnOVarID, "long_name", "mixing ratio O"))
-call check(nf90_put_att(ncFileID, qnOVarID, "units", "?"))
+! call nc_check(nf90_def_var(ncid=ncFileID, name="qnO", xtype=nf90_real, &
+! dimids = (/ levDimID, lonDimID, latDimID, MemberDimID, unlimitedDimID /), &
+! varid = qnOVarID),'nc_write_model_atts')
+! call nc_check(nf90_put_att(ncFileID, qnOVarID, "long_name", "mixing ratio O"),'nc_write_model_atts')
+! call nc_check(nf90_put_att(ncFileID, qnOVarID, "units", "?"),'nc_write_model_atts')
-call check(nf90_enddef(ncfileID))
+call nc_check(nf90_enddef(ncfileID),'nc_write_model_atts')
!-------------------------------------------------------------------------------
! Fill the variables
!-------------------------------------------------------------------------------
-call check(nf90_put_var(ncFileID, lonVarID, lons))
-call check(nf90_put_var(ncFileID, latVarID, lats))
-call check(nf90_put_var(ncFileID, levVarID, levs))
+call nc_check(nf90_put_var(ncFileID, lonVarID, lons),'nc_write_model_atts')
+call nc_check(nf90_put_var(ncFileID, latVarID, lats),'nc_write_model_atts')
+call nc_check(nf90_put_var(ncFileID, levVarID, levs),'nc_write_model_atts')
!-------------------------------------------------------------------------------
! Flush the buffer and leave netCDF file open
!-------------------------------------------------------------------------------
-call check(nf90_sync(ncFileID))
+call nc_check(nf90_sync(ncFileID),'nc_write_model_atts')
write (*,*)'nc_write_model_atts: netCDF file ',ncFileID,' is synched ...'
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-contains
-
-! Internal subroutine - checks error status after each netcdf, prints
-! text message each time an error code is returned.
-
-subroutine check(istatus)
- integer, intent (in) :: istatus
- if(istatus /= nf90_noerr) call error_handler(E_ERR, 'nc_write_model_atts', &
- trim(nf90_strerror(istatus)), source, revision, revdate)
-end subroutine check
-
-
end function nc_write_model_atts
@@ -640,9 +837,6 @@
! NF90_put_var ! provide values for variable
! NF90_CLOSE ! close: save updated netCDF dataset
-use typeSizes
-use netcdf
-
integer, intent(in) :: ncFileID ! netCDF file identifier
real(r8), dimension(:), intent(in) :: statevec
integer, intent(in) :: copyindex
@@ -651,70 +845,62 @@
integer :: nDimensions, nVariables, nAttributes, unlimitedDimID
integer :: u1VarID, v1VarID, t1VarID, uVarID, vVarID, tVarID
-integer :: qnHVarID, qnOHVarID, qnOVarID
+! integer :: qnHVarID, qnOHVarID, qnOVarID
+
type(model_type) :: Var
+if ( .not. module_initialized ) call static_init_model
+
ierr = 0 ! assume normal termination
-call check(nf90_Inquire(ncFileID, nDimensions, nVariables, nAttributes, unlimitedDimID))
+call nc_check(nf90_Inquire(ncFileID, nDimensions, nVariables, nAttributes, unlimitedDimID),'nc_write_model_vars','inquire')
call init_model_instance(Var)
call vector_to_prog_var(statevec, Var)
-call check(NF90_inq_varid(ncFileID, "u1", u1VarID))
-call check(nf90_put_var( ncFileID, u1VarId, var%vars_3d(:,:,:, 1), &
- start=(/ 1, 1, 1, copyindex, timeindex /) ))
+call nc_check(NF90_inq_varid(ncFileID, "u1", u1VarID),'nc_write_model_vars','inq_varid u1')
+call nc_check(nf90_put_var( ncFileID, u1VarId, var%vars_3d(:,:,:, 1), &
+ start=(/ 1, 1, 1, copyindex, timeindex /) ),'nc_write_model_vars','put_var u1')
-call check(NF90_inq_varid(ncFileID, "v1", v1VarID))
-call check(nf90_put_var( ncFileID, v1VarId, var%vars_3d(:,:,:, 2), &
- start=(/ 1, 1, 1, copyindex, timeindex /) ))
+call nc_check(NF90_inq_varid(ncFileID, "v1", v1VarID),'nc_write_model_vars','inq_varid v1')
+call nc_check(nf90_put_var( ncFileID, v1VarId, var%vars_3d(:,:,:, 2), &
+ start=(/ 1, 1, 1, copyindex, timeindex /) ),'nc_write_model_vars','put_var v1')
-call check(NF90_inq_varid(ncFileID, "t1", t1VarID))
-call check(nf90_put_var( ncFileID, t1VarId, var%vars_3d(:,:,:, 3), &
- start=(/ 1, 1, 1, copyindex, timeindex /) ))
+call nc_check(NF90_inq_varid(ncFileID, "t1", t1VarID),'nc_write_model_vars','inq_varid t1')
+call nc_check(nf90_put_var( ncFileID, t1VarId, var%vars_3d(:,:,:, 3), &
+ start=(/ 1, 1, 1, copyindex, timeindex /) ),'nc_write_model_vars','put_var t1')
-call check(NF90_inq_varid(ncFileID, "u", uVarID))
-call check(nf90_put_var( ncFileID, uVarId, var%vars_3d(:,:,:, 4), &
- start=(/ 1, 1, 1, copyindex, timeindex /) ))
+call nc_check(NF90_inq_varid(ncFileID, "u", uVarID),'nc_write_model_vars','inq_ varid u')
+call nc_check(nf90_put_var( ncFileID, uVarId, var%vars_3d(:,:,:, 4), &
+ start=(/ 1, 1, 1, copyindex, timeindex /) ),'nc_write_model_vars','put_var u')
-call check(NF90_inq_varid(ncFileID, "v", vVarID))
-call check(nf90_put_var( ncFileID, vVarId, var%vars_3d(:,:,:, 5), &
- start=(/ 1, 1, 1, copyindex, timeindex /) ))
+call nc_check(NF90_inq_varid(ncFileID, "v", vVarID),'nc_write_model_vars','inq_varid v')
+call nc_check(nf90_put_var( ncFileID, vVarId, var%vars_3d(:,:,:, 5), &
+ start=(/ 1, 1, 1, copyindex, timeindex /) ),'nc_write_model_vars','put_var v')
-call check(NF90_inq_varid(ncFileID, "t", tVarID))
-call check(nf90_put_var( ncFileID, tVarId, var%vars_3d(:,:,:, 6), &
- start=(/ 1, 1, 1, copyindex, timeindex /) ))
+call nc_check(NF90_inq_varid(ncFileID, "t", tVarID),'nc_write_model_vars','inq_varid t')
+call nc_check(nf90_put_var( ncFileID, tVarId, var%vars_3d(:,:,:, 6), &
+ start=(/ 1, 1, 1, copyindex, timeindex /) ),'nc_write_model_vars','put_var t')
-call check(NF90_inq_varid(ncFileID, "qnH", qnHVarID))
-call check(nf90_put_var( ncFileID, qnHVarId, var%vars_3d(:,:,:, 7), &
- start=(/ 1, 1, 1, copyindex, timeindex /) ))
+! call nc_check(NF90_inq_varid(ncFileID, "qnH", qnHVarID),'nc_write_model_vars','inq_varid qnH')
+! call nc_check(nf90_put_var( ncFileID, qnHVarId, var%vars_3d(:,:,:, 7), &
+! start=(/ 1, 1, 1, copyindex, timeindex /) ),'nc_write_model_vars','put_var qnH')
-call check(NF90_inq_varid(ncFileID, "qnOH", qnOHVarID))
-call check(nf90_put_var( ncFileID, qnOHVarId, var%vars_3d(:,:,:, 8), &
- start=(/ 1, 1, 1, copyindex, timeindex /) ))
+! call nc_check(NF90_inq_varid(ncFileID, "qnOH", qnOHVarID),'nc_write_model_vars','inq_varid qnOH')
+! call nc_check(nf90_put_var( ncFileID, qnOHVarId, var%vars_3d(:,:,:, 8), &
+! start=(/ 1, 1, 1, copyindex, timeindex /) ),'nc_write_model_vars','put_var qnOH')
-call check(NF90_inq_varid(ncFileID, "qnH", qnOVarID))
-call check(nf90_put_var( ncFileID, qnOVarId, var%vars_3d(:,:,:, 9), &
- start=(/ 1, 1, 1, copyindex, timeindex /) ))
+! call nc_check(NF90_inq_varid(ncFileID, "qnH", qnOVarID),'nc_write_model_vars','inq_varid qnH')
+! call nc_check(nf90_put_var( ncFileID, qnOVarId, var%vars_3d(:,:,:, 9), &
+! start=(/ 1, 1, 1, copyindex, timeindex /) ),'nc_write_model_vars','put_var qnH')
write (*,*)'Finished filling variables ...'
-call check(nf90_sync(ncFileID))
+call nc_check(nf90_sync(ncFileID),'nc_write_model_vars','sync')
write (*,*)'netCDF file is synched ...'
call end_model_instance(Var) ! should avoid any memory leaking
-contains
-
- ! Internal subroutine - checks error status after each netcdf, prints
- ! text message each time an error code is returned.
- subroutine check(istatus)
- integer, intent ( in) :: istatus
- if(istatus /= nf90_noerr) call error_handler(E_ERR, 'nc_write_model_vars', &
- trim(nf90_strerror(istatus)), source, revision, revdate)
-
- end subroutine check
-
end function nc_write_model_vars
@@ -738,6 +924,8 @@
integer :: i, variable_type
type(location_type) :: temp_loc
+if ( .not. module_initialized ) call static_init_model
+
! An interface is provided
interf_provided = .true.
@@ -749,9 +937,12 @@
do i = 1, get_model_size()
call get_state_meta_data(i, temp_loc, variable_type)
- if(variable_type == TYPE_U .or. variable_type == TYPE_V .or. variable_type == TYPE_T .or. &
- variable_type == TYPE_U0 .or. variable_type == TYPE_V0 .or. variable_type == TYPE_T0) then
- pert_state(i) = random_gaussian(random_seq, state(i), 0.5_r8)
+ if(variable_type == KIND_U_WIND_COMPONENT .or. &
+ variable_type == KIND_V_WIND_COMPONENT .or. &
+ variable_type == KIND_TEMPERATURE ) then
+ pert_state(i) = &
+! & random_gaussian(random_seq, state(i), state(i)*0.05)
+ & state(i)
else
pert_state(i) = state(i)
endif
@@ -770,6 +961,8 @@
integer :: i, j, k, nf, indx
+if ( .not. module_initialized ) call static_init_model
+
! Do order as ps, t, u, v, q, tracers to be consistent with b-grid
! Start copying fields to straight vector
@@ -779,7 +972,7 @@
indx = 0
do i = 1, nx !longitude
do j = 1, ny !latitude
-! u,v,t,q, and tracers at successively lower levels
+ ! u,v,t,q, and tracers at successively lower levels
do k = 1, nz ! height
do nf= 1, state_num_3d
indx = indx + 1
@@ -791,8 +984,8 @@
! Temporary check
if(indx /= model_size) then
- write(errstring, *) 'indx ',indx,' model_size ',model_size,' must be equal '
- call error_handler(E_ERR, 'prog_var_to_vector', errstring, source, revision, revdate)
+ write(msgstring, *) 'indx ',indx,' model_size ',model_size,' must be equal '
+ call error_handler(E_ERR, 'prog_var_to_vector', msgstring, source, revision, revdate)
endif
end subroutine prog_var_to_vector
@@ -807,8 +1000,9 @@
type(model_type), intent(out) :: var
integer :: i, j, k, nf, indx
-character(len=129) :: errstring
+if ( .not. module_initialized ) call static_init_model
+
! Start copying fields from straight vector
! TJH NOTE: the 'other' 3D models have variables that
! are allocated [lon,lat,lev ...] which is consistent with the Fortran
@@ -816,7 +1010,7 @@
indx = 0
do i = 1, nx ! longitude
do j = 1, ny ! latitude
-! u,v,t,q and tracers at successive levels
+ ! u,v,t,q and tracers at successive levels
do k = 1, nz ! height
do nf = 1, state_num_3d
indx = indx + 1
@@ -828,12 +1022,14 @@
! Temporary check
if(indx /= model_size) then
- write(errstring, *) 'indx ',indx,' model_size ',model_size,' must be equal '
- call error_handler(E_ERR, 'vector_to_prog_var', errstring, source, revision, revdate)
+ write(msgstring, *) 'indx ',indx,' model_size ',model_size,' must be equal '
+ call error_handler(E_ERR, 'vector_to_prog_var', msgstring, source, revision, revdate)
endif
end subroutine vector_to_prog_var
+
+
subroutine init_model_instance(var)
!=======================================================================
! subroutine init_model_instance(var)
@@ -842,6 +1038,8 @@
type(model_type), intent(out) :: var
+if ( .not. module_initialized ) call static_init_model
+
! Initialize the storage space and return
! TJH NOTE: the 'other' 3D models have variables that
! are allocated [lon,lat,lev ...] which is consistent with the Fortran
@@ -849,11 +1047,10 @@
allocate(var%vars_3d(nz, nx, ny, state_num_3d))
+end subroutine init_model_instance
-end subroutine init_model_instance
-
subroutine end_model_instance(var)
!=======================================================================
! subroutine end_model_instance(var)
@@ -862,121 +1059,86 @@
type(model_type), intent(inout) :: var
+if ( .not. module_initialized ) call static_init_model
+
deallocate(var%vars_3d)
end subroutine end_model_instance
+
subroutine update_ROSE_restart(file_name, var)
!=======================================================================
! update ROSE restart file fields
!
- use restart, only : var_read, dim_check
-
- use ncdf
- include 'netcdf.inc'
-
character (len = *), intent(in) :: file_name
- type(model_type), intent(in) :: var
+ type(model_type), intent(in) :: var
! ROSE restart file prognostic variables
-! from a3d.chem.f
- integer :: doy ! day of year
- integer :: utsec ! universal time of model in seconds
- integer :: cal_year
-! from dynam.mod.f
real, dimension (nz,nx,ny) :: un1, vn1, tn1, un0, vn0, tn0
- real, dimension (nz) :: tref
- real :: treflb, trefub
-! from params.mod.f
- integer :: tot_tstep ! total timesteps
- ! since beginning of 1st segment
- integer :: year0, day0, ut0, nseg
-! from chem.mod.f ... constituent mixing ratios
+
! real, dimension (nz,nx,ny,nbcon) :: qn1
! real, dimension (nz,ny) :: q_o2, q_n2
-! from restart.mod.f
integer :: restart_id
- integer :: nstp_id
- integer :: time_id
- integer :: mtime_id
- integer, dimension(6) :: rvar_ids
- character (LEN=6), dimension(6) :: var_names = &
- & (/ 'un0 ', 'vn0 ', 'tn0 ', &
- & 'un1 ', 'vn1 ', 'tn1 ' /)
- character (LEN=6), dimension(6) :: var_units = &
- & (/ 'm/s ', 'm/s ', 'deg k ' ,&
- & 'm/s ', 'm/s ', 'deg k ' /)
- real, dimension(nz, nx, ny, 6) :: vars
- integer :: iret, idv, idd
- integer :: ncerr
+ integer :: idd
+ integer :: ncerr, nlons, nlats, nlevs
integer :: var_id
- integer, dimension(3) :: mtime
!====================================================================
-if(file_exist(file_name)) then
+if ( .not. module_initialized ) call static_init_model
+if( .not. file_exist(file_name)) then
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list