! DART software - Copyright 2004 - 2013 UCAR. This open source software is ! provided by UCAR, "as is", without charge, subject to all terms of use at ! http://www.image.ucar.edu/DAReS/DART/DART_download ! ! $Id: model_mod.f90 6256 2013-06-12 16:19:10Z thoar $ module model_mod !----------------------------------------------------------------------- ! ! Interface for ROSE model ! !----------------------------------------------------------------------- ! DART Modules use types_mod, only : r8, digits12, pi use time_manager_mod, only : time_type, set_time, print_time, get_time, & operator(<), operator(>), operator(+), & operator(-), operator(/), operator(*), & operator(==), operator(/=), set_time_missing 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, 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, 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 !----------------------------------------------------------------------- implicit none private public :: get_model_size, & adv_1step, & get_state_meta_data, & model_interpolate, & get_model_time_step, & end_model, & static_init_model, & init_time, & init_conditions, & nc_write_model_atts, & nc_write_model_vars, & pert_model_state, & get_close_maxdist_init, & get_close_obs_init, & get_close_obs, & ens_mean_for_model ! ROSE-specific routines public :: model_type, & prog_var_to_vector, & vector_to_prog_var, & read_ROSE_restart, & update_ROSE_restart, & update_ROSE_namelist, & read_ROSE_tref,& init_model_instance, & end_model_instance !----------------------------------------------------------------------- type model_type real(r8), pointer :: vars_3d(:,:,:,:) type(time_type) :: valid_time 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 ! Arrays to store lats, lons real(r8) :: lons(nx), lats(ny), levs(nz) ! Random sequence and init for pert_model_state logical :: first_pert_call = .true. type(random_seq_type) :: random_seq !---------------------------------------------------------------------- ! Nameslist variables with default values follow ! Namelist variables for defining state vector, and default values integer :: state_num_3d = 6 ! # of 3d fields to read from file namelist /model_nml/ state_num_3d logical :: output_prog_diag = .false. character(len=128) :: input_dir = '../input_current/' character(len=50) :: out_dir = '/ptmp/tmatsuo/rose/' character(len=30) :: restart_file = 'dyn_restart_001-1999.nc' real(kind=r8) :: amp_tune = 1. real(kind=r8) :: pha_tune = 0. real(kind=digits12) :: target_time = 0.125 ! [hr] integer :: ens_element = 1 namelist /ROSE_NML/ target_time, & input_dir, out_dir, restart_file,& output_prog_diag, & amp_tune, pha_tune, & ntime, ens_element integer, parameter :: TYPE_U0 = 0, TYPE_V0 = 1, TYPE_T0 = 2, & TYPE_U = 3, TYPE_V = 4, TYPE_T = 5, & TYPE_Q_H = 6, TYPE_Q_OH = 7, TYPE_Q_O = 8 !---------------------------------------------------------------------- ! version controlled file description for error handling, do not edit character(len=256), parameter :: source = & "$URL: https://svn-dares-dart.cgd.ucar.edu/DART/releases/classic/models/rose/model_mod.f90 $" character(len=32 ), parameter :: revision = "$Revision: 6256 $" character(len=128), parameter :: revdate = "$Date: 2013-06-12 10:19:10 -0600 (Wed, 12 Jun 2013) $" character(len = 129) :: msgstring logical, save :: module_initialized = .false. contains !================================================================== subroutine static_init_model() !------------------------------------------------------------------ ! ! Called to do one time initialization of the model. As examples, ! might define information about the model size or model timestep. ! 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, 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 ! Read the namelist ROSE_NML from the file rose.nml call find_namelist_in_file("rose.nml", "ROSE_NML", iunit) read(iunit, nml = ROSE_NML, iostat = io) call check_namelist_read(iunit, io, "ROSE_NML") 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) call check_namelist_read(iunit, io, "model_nml") if (do_nml_file()) write(nmlfileunit, nml=model_nml) if (do_nml_term()) write( * , nml=model_nml) ! Compute overall model size and put in global storage model_size = nx * ny * nz * state_num_3d ! Set the model minimum time step from the namelist seconds and days input ! ROSE control variable "ntime" specify the number of steps per hour Time_step_seconds = nint(3600.0_r8 / real(ntime)) Time_step_days = 0 write(*, *) 'time step secs days ' , Time_step_seconds, Time_step_days if (Time_step_seconds > seconds_of_day) then Time_step_days = floor(real(Time_step_seconds) / real(seconds_of_day)) Time_step_seconds = mod(Time_step_seconds, Time_step_days * seconds_of_day) endif Time_step_ROSE = set_time(Time_step_seconds, Time_step_days) if (do_output()) call print_time(Time_step_ROSE,'ROSE time step') ! lon: long_name = "geographic longitude", units = "degrees" ; d_lon = 360.0_r8/real(nx) do i = 1, nx lons(i) = (i-1)*d_lon enddo ! lat: long_name = "geographic latitude", units = "degrees" ; d_lat = 180.0_r8/real(ny) do j = 1, ny lats(j) = -90 + d_lat/2. + (j-1)*d_lat enddo ! Pressure levels ! from msetfix.f ! 1013 [mb] exp(-z/7km) with a fixed scale height (7km) ! (2.1 km increment from 16.8 km) do k = 1, nz z_m = (k-1)*dz + zbot 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 subroutine init_conditions(x) !------------------------------------------------------------------ ! subroutine init_conditions(x) ! ! Returns a model state vector, x, that is some sort of appropriate ! initial condition for starting up a long integration of the model. ! At present, this is only used if the namelist parameter ! start_from_restart is set to .false. in the program perfect_model_obs. ! If this option is not to be used in perfect_model_obs, or if no ! synthetic data experiments using perfect_model_obs are planned, ! this can be a NULL INTERFACE. real(r8), intent(out) :: x(:) if ( .not. module_initialized ) call static_init_model end subroutine init_conditions subroutine adv_1step(x, time) !------------------------------------------------------------------ ! subroutine adv_1step(x, time) ! ! Does a single timestep advance of the model. The input value of ! the vector x is the starting condition and x is updated to reflect ! the changed state after a timestep. The time argument is intent ! in and is used for models that need to know the date/time to ! compute a timestep, for instance for radiation compuations. ! This interface is only called if the namelist parameter ! async is set to 0 in perfect_model_obs of filter or if the ! program integrate_model is to be used to advance the model ! state as a separate executable. If one of these options ! is not going to be used (the model will only be advanced as ! a separate model-specific executable), this can be a ! NULL INTERFACE. real(r8), intent(inout) :: x(:) type(time_type), intent(in) :: time if ( .not. module_initialized ) call static_init_model end subroutine adv_1step function get_model_size() !------------------------------------------------------------------ ! ! Returns the size of the model as an integer. Required for all ! applications. integer :: get_model_size if ( .not. module_initialized ) call static_init_model get_model_size = model_size end function get_model_size subroutine init_time(time) !------------------------------------------------------------------ ! ! Companion interface to init_conditions. Returns a time that is somehow ! appropriate for starting up a long integration of the model. ! At present, this is only used if the namelist parameter ! start_from_restart is set to .false. in the program perfect_model_obs. ! If this option is not to be used in perfect_model_obs, or if no ! synthetic data experiments using perfect_model_obs are planned, ! this can be a NULL INTERFACE. 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 subroutine model_interpolate(x, location, itype, obs_val, istatus) !------------------------------------------------------------------ ! ! Given a state vector, a location, and a model state variable type, ! interpolates the state variable field`to that location and returns ! the value in obs_val. The istatus variable should be returned as ! 0 unless there is some problem in computing the interpolation in ! which case an alternate value should be returned. The itype variable ! is a model specific integer that specifies the type of field (for ! instance temperature, zonal wind component, etc.). In low order ! models that have no notion of types of variables, this argument can ! be ingored. For applications in which only perfect model experiments ! with identity observations (i.e. only the value of a particular ! state variable is observerd), this can be a NULL INTERFACE. real(r8), intent(in) :: x(:) type(location_type), intent(in) :: location integer, intent(in) :: itype 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() !------------------------------------------------------------------ ! ! Returns the the time step of the model; the smallest increment ! in time that the model is capable of advancing the state in a given ! implementation. This interface is required for all applications. ! ! Limited by ROSE's fixed time step (also by 2-step restart file) 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 subroutine get_state_meta_data(index_in, location, var_type) !------------------------------------------------------------------ ! ! Given an integer index into the state vector structure, returns the ! associated location. A second intent(out) optional argument kind ! can be returned if the model has more than one type of field (for ! instance temperature and zonal wind component). This interface is ! required for all filter applications as it is required for computing ! the distance between observations and state variables. integer, intent(in) :: index_in type(location_type), intent(out) :: location integer, intent(out), optional :: var_type integer :: indx, num_per_col, col_num, col_elem, lon_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 ! Compute number of items per column num_per_col = nz * state_num_3d ! What column is this index in col_num = indx / num_per_col col_elem = indx - col_num * num_per_col ! What lon and lat index for this column lon_index = col_num / ny ! ny is number of ROSE latitude grid lat_index = col_num - lon_index * ny ! Get actual lon lat values from static_init arrays ??? lon = lons(lon_index + 1) lat = lats(lat_index + 1) ! Now figure out which beast in column this is lev_index = col_elem / state_num_3d lev = levs(lev_index + 1) ! 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 = KIND_U_WIND_COMPONENT else if(var_type_temp == 1) then local_var_type = KIND_V_WIND_COMPONENT else if(var_type_temp == 2) then local_var_type = KIND_TEMPERATURE else if(var_type_temp == 3) then local_var_type = KIND_U_WIND_COMPONENT else if(var_type_temp == 4) then local_var_type = KIND_V_WIND_COMPONENT else if(var_type_temp == 5) then 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 location = set_location(lon, lat, lev, 2) ! 2 == pressure (hPa) ! If the type is wanted, return it if(present(var_type)) var_type = local_var_type end subroutine get_state_meta_data subroutine end_model() !------------------------------------------------------------------ ! ! 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 function nc_write_model_atts( ncFileID ) result (ierr) !------------------------------------------------------------------ ! Writes the model-specific attributes to a netCDF file ! TJH Jan 24 2003 ! ! TJH 29 July 2003 -- for the moment, all errors are fatal, so the ! return code is always '0 == normal', since the fatal errors stop execution. ! ! For the lorenz_96 model, each state variable is at a separate location. ! that's all the model-specific attributes I can think of ... ! ! assim_model_mod:init_diag_output uses information from the location_mod ! to define the location dimension and variable ID. All we need to do ! is query, verify, and fill ... ! ! Typical sequence for adding new dimensions,variables,attributes: ! NF90_OPEN ! open existing netCDF dataset ! NF90_redef ! put into define mode ! NF90_def_dim ! define additional dimensions (if any) ! NF90_def_var ! define variables: from name, type, and dims ! NF90_put_att ! assign attribute values ! NF90_ENDDEF ! end definitions: leave define mode ! NF90_put_var ! provide values for variable ! NF90_CLOSE ! close: save updated netCDF dataset integer, intent(in) :: ncFileID ! netCDF file identifier integer :: ierr ! return value of function integer :: nDimensions, nVariables, nAttributes, unlimitedDimID integer :: MemberDimID, TimeDimID integer :: lonDimID, latDimID, levDimID integer :: lonVarID, latVarID, levVarID integer :: u1VarID, v1VarID, t1VarID, uVarID, vVarID, tVarID ! 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 character(len=5) :: crzone ! needed by F90 DATE_AND_TIME intrinsic integer, dimension(8) :: values ! needed by F90 DATE_AND_TIME intrinsic character(len=NF90_MAX_NAME) :: str1 !------------------------------------------------------------------ if ( .not. module_initialized ) call static_init_model ierr = 0 ! assume normal termination !-------------------------------------------------------------------- ! make sure ncFileID refers to an open netCDF file !-------------------------------------------------------------------- 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 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(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 !------------------------------------------------------------------------------- ! Write Global Attributes !------------------------------------------------------------------------------- call DATE_AND_TIME(crdate,crtime,crzone,values) 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 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 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 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 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 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 !---------------------------------------------------------------------------- ! TJH NOTE: the 'other' 3D models have netcdf variables that ! are allocated ! dimids = (/ lonDimID, latDimID, levDimID, MemberDimID, unlimitedDimID /), & ! which is consistent with the Fortran convention of having the ! fastest-moving index on the left. call nc_check(nf90_def_var(ncid=ncFileID, name="u1", xtype=nf90_real, & dimids = (/ levDimID, lonDimID, latDimID, MemberDimID, unlimitedDimID /), & 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 nc_check(nf90_def_var(ncid=ncFileID, name="v1", xtype=nf90_real, & dimids = (/ levDimID, lonDimID, latDimID, MemberDimID, unlimitedDimID /), & 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 nc_check(nf90_def_var(ncid=ncFileID, name="t1", xtype=nf90_real, & dimids = (/ levDimID, lonDimID, latDimID, MemberDimID, unlimitedDimID /), & 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 nc_check(nf90_def_var(ncid=ncFileID, name="u", xtype=nf90_real, & dimids = (/ levDimID, lonDimID, latDimID, MemberDimID, unlimitedDimID /), & 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 nc_check(nf90_def_var(ncid=ncFileID, name="v", xtype=nf90_real, & dimids = (/ levDimID, lonDimID, latDimID, MemberDimID, unlimitedDimID /), & 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 nc_check(nf90_def_var(ncid=ncFileID, name="t", xtype=nf90_real, & dimids = (/ levDimID, lonDimID, latDimID, MemberDimID, unlimitedDimID /), & 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 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 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 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 nc_check(nf90_enddef(ncfileID),'nc_write_model_atts') !------------------------------------------------------------------------------- ! Fill the variables !------------------------------------------------------------------------------- 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 nc_check(nf90_sync(ncFileID),'nc_write_model_atts') write (*,*)'nc_write_model_atts: netCDF file ',ncFileID,' is synched ...' end function nc_write_model_atts function nc_write_model_vars( ncFileID, statevec, copyindex, timeindex ) result (ierr) !------------------------------------------------------------------ ! Writes the model-specific attributes to a netCDF file ! TJH 23 May 2003 ! ! TJH 29 July 2003 -- for the moment, all errors are fatal, so the ! return code is always '0 == normal', since the fatal errors stop execution. ! ! For the lorenz_96 model, each state variable is at a separate location. ! that's all the model-specific attributes I can think of ... ! ! assim_model_mod:init_diag_output uses information from the location_mod ! to define the location dimension and variable ID. All we need to do ! is query, verify, and fill ... ! ! Typical sequence for adding new dimensions,variables,attributes: ! NF90_OPEN ! open existing netCDF dataset ! NF90_redef ! put into define mode ! NF90_def_dim ! define additional dimensions (if any) ! NF90_def_var ! define variables: from name, type, and dims ! NF90_put_att ! assign attribute values ! NF90_ENDDEF ! end definitions: leave define mode ! NF90_put_var ! provide values for variable ! NF90_CLOSE ! close: save updated netCDF dataset integer, intent(in) :: ncFileID ! netCDF file identifier real(r8), dimension(:), intent(in) :: statevec integer, intent(in) :: copyindex integer, intent(in) :: timeindex integer :: ierr ! return value of function integer :: nDimensions, nVariables, nAttributes, unlimitedDimID integer :: u1VarID, v1VarID, t1VarID, uVarID, vVarID, tVarID ! integer :: qnHVarID, qnOHVarID, qnOVarID type(model_type) :: Var if ( .not. module_initialized ) call static_init_model ierr = 0 ! assume normal termination 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 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 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 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 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 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 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 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 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 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') if (do_output()) write (*,*)'Finished filling variables ...' call nc_check(nf90_sync(ncFileID),'nc_write_model_vars','sync') if (do_output()) write (*,*)'netCDF file is synched ...' call end_model_instance(Var) ! should avoid any memory leaking end function nc_write_model_vars subroutine pert_model_state(state, pert_state, interf_provided) !------------------------------------------------------------------ ! ! Perturbs a model state for generating initial ensembles. ! The perturbed state is returned in pert_state. ! A model may choose to provide a NULL INTERFACE by returning ! .false. for the inter_provided argument. This indicates to ! the filter that if it needs to generate perturbed states, it ! may do so by adding an O(0.1) magnitude perturbation to each ! model state variable independently. The interf_provided argument ! should be returned as .true. if the model want to do its own ! perturbing of states. real(r8), intent(in) :: state(:) real(r8), intent(out) :: pert_state(:) logical, intent(out) :: interf_provided integer :: i, variable_type type(location_type) :: temp_loc if ( .not. module_initialized ) call static_init_model ! An interface is provided interf_provided = .true. ! If first call initialize random sequence if(first_pert_call) then call init_random_seq(random_seq) first_pert_call = .false. endif do i = 1, get_model_size() call get_state_meta_data(i, temp_loc, variable_type) 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 end do end subroutine pert_model_state subroutine prog_var_to_vector(var, x) !======================================================================= ! based on prog_var_to_vector in model_mod for CAM. ! type(model_type), intent(in) :: var real(r8), intent(out) :: x(:) 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 ! TJH NOTE: the 'other' 3D models have variables that ! are allocated [lon,lat,lev ...] which is consistent with the Fortran ! convention of having the fastest-moving index on the left. indx = 0 do i = 1, nx !longitude do j = 1, ny !latitude ! u,v,t,q, and tracers at successively lower levels do k = 1, nz ! height do nf= 1, state_num_3d indx = indx + 1 x(indx) = var%vars_3d(k, i, j, nf) end do end do end do end do ! Temporary check if(indx /= model_size) then 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 subroutine vector_to_prog_var(x, var) !======================================================================= ! subroutine vector_to_prog_var(x, var) ! real(r8), intent(in) :: x(:) type(model_type), intent(out) :: var integer :: i, j, k, nf, indx 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 ! convention of having the fastest-moving index on the left. indx = 0 do i = 1, nx ! longitude do j = 1, ny ! latitude ! u,v,t,q and tracers at successive levels do k = 1, nz ! height do nf = 1, state_num_3d indx = indx + 1 var%vars_3d(k, i, j, nf) = x(indx) end do end do end do end do ! Temporary check if(indx /= model_size) then 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, valid_time) !======================================================================= ! subroutine init_model_instance(var) ! ! Initializes an instance of a ROSE model state variable type(model_type), intent(out) :: var type(time_type), optional, intent( in) :: valid_time 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 ! convention of having the fastest-moving index on the left. allocate(var%vars_3d(nz, nx, ny, state_num_3d)) if (present(valid_time)) then var%valid_time = valid_time else var%valid_time = set_time_missing() endif end subroutine init_model_instance subroutine end_model_instance(var) !======================================================================= ! subroutine end_model_instance(var) ! ! Ends an instance of a rose model state variable 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 character (len = *), intent(in) :: file_name type(model_type), intent(in) :: var ! ROSE restart file prognostic variables real, dimension (nz,nx,ny) :: un1, vn1, tn1, un0, vn0, tn0 ! real, dimension (nz,nx,ny,nbcon) :: qn1 ! real, dimension (nz,ny) :: q_o2, q_n2 integer :: restart_id integer :: idd integer :: ncerr, nlons, nlats, nlevs integer :: var_id integer, dimension(3) :: mtime integer :: seconds, days !==================================================================== if ( .not. module_initialized ) call static_init_model if( .not. file_exist(file_name)) then write(msgstring,*) trim(adjustl(file_name)),' not available.' call error_handler(E_ERR,'update_ROSE_restart',msgstring,source,revision,revdate) endif if (do_output()) print *, 'update_ROSE_restart: reading restart' ncerr = nf90_open( file_name, NF90_WRITE, restart_id ) call nc_check(ncerr, 'update_ROSE_restart','open') ! will die if error if (do_output()) print *, 'update_ROSE_restart: opened with '//trim(nf90_strerror(ncerr)) !... check for matching dimensions call nc_check( nf90_inq_dimid(restart_id, 'lon', idd), & 'update_rose_restart', 'inq_dimid lon') call nc_check( nf90_inquire_dimension(restart_id, idd, len=nlons), & 'update_rose_restart', 'inquire_dimension lon') call nc_check( nf90_inq_dimid(restart_id, 'lat', idd), & 'update_rose_restart', 'inq_dimid lat') call nc_check( nf90_inquire_dimension(restart_id, idd, len=nlats), & 'update_rose_restart', 'inquire_dimension lat') call nc_check( nf90_inq_dimid(restart_id, 'lev', idd), & 'update_rose_restart', 'inq_dimid lev') call nc_check( nf90_inquire_dimension(restart_id, idd, len=nlevs), & 'update_rose_restart', 'inquire_dimension lev') if ( nlons /= nx ) then write(msgstring,*) trim(file_name),' nlons = ',nlons,' DART expects ',nx call error_handler(E_ERR,'update_rose_restart',msgstring,source,revision,revdate) endif if ( nlats /= ny ) then write(msgstring,*) trim(file_name),' nlats = ',nlats,' DART expects ',ny call error_handler(E_ERR,'update_rose_restart',msgstring,source,revision,revdate) endif if ( nlevs /= nz ) then write(msgstring,*) trim(file_name),' nlevs = ',nlevs,' DART expects ',nz call error_handler(E_ERR,'update_rose_restart',msgstring,source,revision,revdate) endif ! unpack the variables into something familiar and then stuff them ! into the existing netCDF variables. call get_time(var%valid_time, seconds, days) un1 = var%vars_3d(:,:,:, 1) vn1 = var%vars_3d(:,:,:, 2) tn1 = var%vars_3d(:,:,:, 3) un0 = var%vars_3d(:,:,:, 4) vn0 = var%vars_3d(:,:,:, 5) tn0 = var%vars_3d(:,:,:, 6) !qn1(:,:,:,7) = var%vars_3d(:,:,:,7) ! H !qn1(:,:,:,8) = var%vars_3d(:,:,:,8) ! OH !qn1(:,:,:,18) = var%vars_3d(:,:,:,9) ! O call nc_check( nf90_inq_varid(restart_id, 'un1', var_id), & 'update_rose_restart', 'inq_varid un1') call nc_check( nf90_put_var(restart_id, var_id, values=un1), & 'update_rose_restart', 'put_var un1') call nc_check( nf90_inq_varid(restart_id, 'vn1', var_id), & 'update_rose_restart', 'inq_varid vn1') call nc_check( nf90_put_var(restart_id, var_id, values=vn1), & 'update_rose_restart', 'put_var vn1') call nc_check( nf90_inq_varid(restart_id, 'tn1', var_id), & 'update_rose_restart', 'inq_varid tn1') call nc_check( nf90_put_var(restart_id, var_id, values=tn1), & 'update_rose_restart', 'put_var tn1') call nc_check( nf90_inq_varid(restart_id, 'un0', var_id), & 'update_rose_restart', 'inq_varid un0') call nc_check( nf90_put_var(restart_id, var_id, values=un0), & 'update_rose_restart', 'put_var un0') call nc_check( nf90_inq_varid(restart_id, 'vn0', var_id), & 'update_rose_restart', 'inq_varid vn0') call nc_check( nf90_put_var(restart_id, var_id, values=vn0), & 'update_rose_restart', 'put_var vn0') call nc_check( nf90_inq_varid(restart_id, 'tn0', var_id), & 'update_rose_restart', 'inq_varid tn0') call nc_check( nf90_put_var(restart_id, var_id, values=tn0), & 'update_rose_restart', 'put_var tn0') ! FIXME - since we're not handling the year correctly in rose, ! we need to capture the existing one and reuse it. call nc_check( nf90_inq_varid(restart_id, 'mtime', var_id), & 'update_rose_restart','inq_varid mtime') call nc_check(nf90_get_var( restart_id, var_id, mtime) , & 'update_rose_restart','get_var mtime') ! mtime(1) = cal_year FIXME ... this is that I'm talking about mtime(2) = days mtime(3) = seconds call nc_check(nf90_put_var( restart_id, var_id, mtime) , & 'update_rose_restart','get_var mtime') if (do_output()) print *, 'update_ROSE_restart: mtime (year/doy/seconds):', mtime call nc_check( nf90_sync( restart_id), 'update_rose_restart', 'sync') call nc_check( nf90_close(restart_id), 'update_rose_restart', 'close') end subroutine update_ROSE_restart subroutine read_ROSE_restart(file_name, var, model_time) !======================================================================= ! read ROSE restart file fields that have been updated ! character (len = *), intent(in) :: file_name type(model_type), intent(out) :: var type(time_type), intent(out) :: model_time integer :: ncerr, idd integer :: var_id, restart_id, nlons, nlats, nlevs integer, dimension(3) :: mtime ! ROSE restart file prognostic variables integer :: doy ! day of year integer :: utsec integer :: cal_year ! calendar year real, dimension (nz,nx,ny) :: un1, vn1, tn1, un0, vn0, tn0 ! real, dimension (nz) :: tref ! from chem.mod.f ... constituent mixing ratios ! real, dimension (nz,nx,ny,nbcon) :: qn1 ! real, dimension (nz,ny) :: q_o2, q_n2 ! !==================================================================== if ( .not. module_initialized ) call static_init_model if( .not. file_exist(file_name)) then write(msgstring,*)trim(file_name)//' does not exist.' call error_handler(E_ERR,'read_ROSE_restart',msgstring,source,revision,revdate) endif if (do_output()) print *, 'read_ROSE_restart:reading restart:', file_name ncerr = nf90_open( file_name, NF90_NOWRITE, restart_id ) call nc_check(ncerr, 'read_ROSE_restart', 'open') if (do_output()) print *, 'read_ROSE_restart:opened with '//trim(nf90_strerror(ncerr)) !... check for matching dimensions ! call dim_check( 'lon', nx) ! call dim_check( 'lat', ny) ! call dim_check( 'lev', nz) !... check for matching dimensions call nc_check( nf90_inq_dimid(restart_id, 'lon', idd), & 'read_rose_restart', 'inq_dimid lon') call nc_check( nf90_inquire_dimension(restart_id, idd, len=nlons), & 'read_rose_restart', 'inquire_dimension lon') call nc_check( nf90_inq_dimid(restart_id, 'lat', idd), & 'read_rose_restart', 'inq_dimid lat') call nc_check( nf90_inquire_dimension(restart_id, idd, len=nlats), & 'read_rose_restart', 'inquire_dimension lat') call nc_check( nf90_inq_dimid(restart_id, 'lev', idd), & 'read_rose_restart', 'inq_dimid lev') call nc_check( nf90_inquire_dimension(restart_id, idd, len=nlevs), & 'read_rose_restart', 'inquire_dimension lev') if ( nlons /= nx ) then write(msgstring,*) trim(file_name),' nlons = ',nlons,' DART expects ',nx call error_handler(E_ERR,'read_rose_restart',msgstring,source,revision,revdate) endif if ( nlats /= ny ) then write(msgstring,*) trim(file_name),' nlats = ',nlats,' DART expects ',ny call error_handler(E_ERR,'read_rose_restart',msgstring,source,revision,revdate) endif if ( nlevs /= nz ) then write(msgstring,*) trim(file_name),' nlevs = ',nlevs,' DART expects ',nz call error_handler(E_ERR,'read_rose_restart',msgstring,source,revision,revdate) endif !... get dynamical variables (un1, vn1, tn1, un0, vn0, tn0) ! call var_read( 'un1', un1) ! call var_read( 'vn1', vn1) ! call var_read( 'tn1', tn1) ! call var_read( 'un0', un0) ! call var_read( 'vn0', vn0) ! call var_read( 'tn0', tn0) call nc_check( nf90_inq_varid(restart_id, 'un1', var_id), & 'read_rose_restart', 'inq_varid un1') call nc_check( nf90_get_var(restart_id, var_id, values=un1), & 'read_rose_restart', 'get_var un1') call nc_check( nf90_inq_varid(restart_id, 'vn1', var_id), & 'read_rose_restart', 'inq_varid vn1') call nc_check( nf90_get_var(restart_id, var_id, values=vn1), & 'read_rose_restart', 'get_var vn1') call nc_check( nf90_inq_varid(restart_id, 'tn1', var_id), & 'read_rose_restart', 'inq_varid tn1') call nc_check( nf90_get_var(restart_id, var_id, values=tn1), & 'read_rose_restart', 'get_var tn1') call nc_check( nf90_inq_varid(restart_id, 'un0', var_id), & 'read_rose_restart', 'inq_varid un0') call nc_check( nf90_get_var(restart_id, var_id, values=un0), & 'read_rose_restart', 'get_var un0') call nc_check( nf90_inq_varid(restart_id, 'vn0', var_id), & 'read_rose_restart', 'inq_varid vn0') call nc_check( nf90_get_var(restart_id, var_id, values=vn0), & 'read_rose_restart', 'get_var vn0') call nc_check( nf90_inq_varid(restart_id, 'tn0', var_id), & 'read_rose_restart', 'inq_varid tn0') call nc_check( nf90_get_var(restart_id, var_id, values=tn0), & 'read_rose_restart', 'get_var tn0') call nc_check(nf90_inq_varid( restart_id, 'mtime', var_id), & 'read_rose_restart','inq_varid mtime') call nc_check(nf90_get_var( restart_id, var_id, mtime) , & 'read_rose_restart','get_var mtime') call nc_check(nf90_inq_varid( restart_id, 'tref', var_id), & 'read_rose_restart','inq_varid tref') call nc_check(nf90_get_var(restart_id, var_id, tref) , & 'read_rose_restart','get_var tref') ! That's it ... call nc_check(nf90_close( restart_id),'read_rose_restart','close') if (do_output()) print *, 'read_ROSE_restart: mtime (year/doy/seconds):', mtime cal_year = mtime(1) doy = mtime(2) utsec = mtime(3) ! repack the variables into our structure var%vars_3d(:,:,:,1) = un1 var%vars_3d(:,:,:,2) = vn1 var%vars_3d(:,:,:,3) = tn1 var%vars_3d(:,:,:,4) = un0 ! at model_time var%vars_3d(:,:,:,5) = vn0 ! at model_time var%vars_3d(:,:,:,6) = tn0 ! at model_time !var%vars_3d(:,:,:,7) = qn1(:,:,:,7) ! H !var%vars_3d(:,:,:,8) = qn1(:,:,:,8) ! OH !var%vars_3d(:,:,:,9) = qn1(:,:,:,18) ! O ! FIXME ... this ignores years - no calendar model_time = set_time(utsec, doy) var%valid_time = model_time if (do_output()) call print_time(model_time, str=" read_ROSE_restart: model_time ") end subroutine read_ROSE_restart subroutine update_ROSE_namelist(file_name, time1, timeN, ens_member, & atune, ptune ) !======================================================================= ! Update the ROSE namelist - especially the new target_time. ! The target_time is actually an offset - so we need to calculate that. !======================================================================= character(len=*), intent(in) :: file_name type(time_type), intent(in) :: time1, timeN integer, intent(in) :: ens_member real(r8), optional, intent(in) :: atune, ptune type(time_type) :: forecast_length integer :: second, day integer :: iunit forecast_length = timeN - time1 call get_time(forecast_length, second, day) target_time = real( day,digits12)*24.0_digits12 + & real(second,digits12)/3600.0_digits12 if (do_output()) then PRINT*,'update_ROSE_namelist: forecast length [days seconds] ', & day, second, ' = hours ', target_time endif ! Update the other namelist parameters ens_element = ens_member if (present(atune)) amp_tune = atune if (present(ptune)) pha_tune = ptune iunit = open_file(file_name, action='write') write(iunit, nml=ROSE_NML) call close_file(iunit) end subroutine update_ROSE_namelist subroutine read_ROSE_tref(file_name) !======================================================================= ! read the reference temperature from a ROSE restart file !==================================================================== character (len = *), intent(in) :: file_name integer :: dim_id, dim_len integer :: var_id, restart_id if( .not. file_exist(file_name)) then write(msgstring,*) trim(adjustl(file_name)),' not available.' call error_handler(E_ERR,'read_ROSE_tref',msgstring,source,revision,revdate) endif call nc_check(nf90_open(file_name, NF90_NOWRITE, restart_id), 'read_ROSE_tref', 'open') call nc_check(nf90_inq_dimid(restart_id, 'lev', dim_id), 'read_ROSE_tref', 'inq_dimid lev') call nc_check(nf90_inquire_dimension(restart_id, dim_id, len=dim_len),'read_ROSE_tref','inquire_dimension lev') if (dim_len .ne. nz) then write(msgstring,'(a,2i5)') 'restart dimension mismatch', dim_len, nz call error_handler(E_ERR,'read_ROSE_tref',msgstring,source,revision,revdate) endif call nc_check(nf90_inq_varid(restart_id, 'tref', var_id),'read_ROSE_tref', 'inq_varid tref') call nc_check(nf90_get_var(restart_id, var_id, tref),'read_ROSE_tref', 'get_var tref') call nc_check(nf90_close(restart_id),'read_ROSE_tref', 'close') end subroutine read_ROSE_tref subroutine ens_mean_for_model(ens_mean) !------------------------------------------------------------------ ! Not used ???? real(r8), intent(in) :: ens_mean(:) if ( .not. module_initialized ) call static_init_model end subroutine ens_mean_for_model !=================================================================== ! End of model_mod !=================================================================== end module model_mod ! ! $URL: https://svn-dares-dart.cgd.ucar.edu/DART/releases/classic/models/rose/model_mod.f90 $ ! $Id: model_mod.f90 6256 2013-06-12 16:19:10Z thoar $ ! $Revision: 6256 $ ! $Date: 2013-06-12 10:19:10 -0600 (Wed, 12 Jun 2013) $