diff --git a/.gitignore b/.gitignore index f92426caa9..5d1f27c778 100644 --- a/.gitignore +++ b/.gitignore @@ -90,6 +90,8 @@ gitm_to_netcdf netcdf_to_gitm_blocks streamflow_obs_diag cam_dart_obs_preprocessor +aether_to_dart +dart_to_aether # Observation converter exectutables convert_aviso diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 84bd662946..5ad6ff2d62 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -22,6 +22,18 @@ individual files. The changes are now listed with the most recent at the top. +**March 12 2024 :: MITgcm/N-BLING with Compressed Staggered Grids. Tag v11.3.0** + +- The DART-MITgcm code now supports compressed grids, especially suited for areas like + the Red Sea where land occupies more than 90% of the domain. + Built upon work *contributed by Jiachen Liu*. +- Allows writing the BGC fields into MITgcm's pickup files. +- Allows different compression for the regular and staggered grids. + +**March 12 2024 :: Aether lat-lon. Tag v11.2.0** + +- Aether lat-lon interface added to DART. + **March 11 2024 :: SEIR model for infectious diseases. Tag v11.1.0** - Added SEIR model which simulates the spread of infectious diseases, for example COVID-19. diff --git a/assimilation_code/modules/utilities/netcdf_utilities_mod.f90 b/assimilation_code/modules/utilities/netcdf_utilities_mod.f90 index 74301bd6f4..667138cfe5 100644 --- a/assimilation_code/modules/utilities/netcdf_utilities_mod.f90 +++ b/assimilation_code/modules/utilities/netcdf_utilities_mod.f90 @@ -63,7 +63,7 @@ module netcdf_utilities_mod nc_begin_define_mode, & nc_end_define_mode, & nc_synchronize_file, & - NF90_MAX_NAME, NF90_MAX_VAR_DIMS + NF90_MAX_NAME, NF90_MAX_VAR_DIMS, NF90_FILL_REAL ! note here that you only need to distinguish between diff --git a/assimilation_code/programs/model_mod_check/model_mod_check.f90 b/assimilation_code/programs/model_mod_check/model_mod_check.f90 index ee4eea51a0..f56c7617a1 100644 --- a/assimilation_code/programs/model_mod_check/model_mod_check.f90 +++ b/assimilation_code/programs/model_mod_check/model_mod_check.f90 @@ -416,7 +416,7 @@ subroutine check_meta_data( iloc ) kind_index=qty_index, & kind_string=qty_string) -write(string1,'("index ",i11," is i,j,k",3(1x,i4)," and is in domain ",i2)') & +write(string1,'("index ",i11," is i,j,k",3(1x,i10)," and is in domain ",i2)') & iloc, ix, iy, iz, dom_id write(string2,'("is quantity ", I4,", ",A)') var_type, trim(qty_string)//' at location' call write_location(0,loc,charstring=string3) @@ -556,9 +556,12 @@ subroutine check_all_meta_data() kind_string=qty_string) ! CLM has (potentially many) columns and needs i7 ish precision - write(string1,'(i11,1x,''i,j,k'',3(1x,i7),'' domain '',i2)') & +! write(string1,'(i11,1x,''i,j,k'',3(1x,i7),'' domain '',i2)') & +! iloc, ix, iy, iz, dom_id + ! EL: integer to short for the new I/O method + ! Change to long int to avoid problems + write(string1,'(i21,1x,''i,j,k'',3(1x,i21),'' domain '',i2)') & iloc, ix, iy, iz, dom_id - call get_state_meta_data(iloc, loc, var_type) metadata_qty_string = trim(get_name_for_quantity(var_type)) diff --git a/conf.py b/conf.py index f55af3ec11..0dd2606f18 100644 --- a/conf.py +++ b/conf.py @@ -21,7 +21,7 @@ author = 'Data Assimilation Research Section' # The full version, including alpha/beta/rc tags -release = '11.1.0' +release = '11.3.0' root_doc = 'index' # -- General configuration --------------------------------------------------- diff --git a/index.rst b/index.rst index f04390488f..d0f3490564 100644 --- a/index.rst +++ b/index.rst @@ -439,6 +439,7 @@ References :hidden: models/9var/readme + models/aether_lat-lon/readme models/am2/readme models/bgrid_solo/readme models/cam-fv/readme diff --git a/models/MITgcm_ocean/model_mod.f90 b/models/MITgcm_ocean/model_mod.f90 index 5ef1900bc6..85d635e520 100644 --- a/models/MITgcm_ocean/model_mod.f90 +++ b/models/MITgcm_ocean/model_mod.f90 @@ -20,7 +20,7 @@ module model_mod get_close_state, get_close_obs, set_location, & VERTISHEIGHT, get_location, is_vertical, & convert_vertical_obs, convert_vertical_state - +! EL use only nc_check was here, deleted for now for testing use utilities_mod, only : error_handler, E_ERR, E_WARN, E_MSG, & logfileunit, get_unit, do_output, to_upper, & find_namelist_in_file, check_namelist_read, & @@ -56,7 +56,10 @@ module model_mod get_index_start, get_index_end, & get_dart_vector_index, get_num_variables, & get_domain_size, & - get_io_clamping_minval + get_io_clamping_minval, get_kind_index + +use netcdf_utilities_mod, only : nc_open_file_readonly, nc_get_variable, & + nc_get_dimension_size, nc_close_file use netcdf @@ -255,9 +258,16 @@ module model_mod ! standard MITgcm namelist and filled in here. integer :: Nx=-1, Ny=-1, Nz=-1 ! grid counts for each field +integer :: comp2d = -1, comp3d=-1, comp3dU = -1, comp3dV = -1 ! size of commpressed variables ! locations of cell centers (C) and edges (G) for each axis. real(r8), allocatable :: XC(:), XG(:), YC(:), YG(:), ZC(:), ZG(:) +real(r4), allocatable :: XC_sq(:), YC_sq(:), XG_sq(:), YG_sq(:) +real(r8), allocatable :: ZC_sq(:) + +integer, allocatable :: Xc_Ti(:), Yc_Ti(:), Zc_Ti(:) +integer, allocatable :: Xc_Ui(:), Yc_Ui(:), Zc_Ui(:) +integer, allocatable :: Xc_Vi(:), Yc_Vi(:), Zc_Vi(:) real(r8) :: ocean_dynamics_timestep = 900.0_r4 integer :: timestepcount = 0 @@ -277,7 +287,6 @@ module model_mod integer, parameter :: NUM_STATE_TABLE_COLUMNS = 5 character(len=vtablenamelength) :: mitgcm_variables(NUM_STATE_TABLE_COLUMNS, MAX_STATE_VARIABLES ) = ' ' - character(len=256) :: model_shape_file = ' ' integer :: assimilation_period_days = 7 integer :: assimilation_period_seconds = 0 @@ -292,8 +301,9 @@ module model_mod logical :: go_to_dart = .false. logical :: do_bgc = .false. logical :: log_transform = .false. +logical :: compress = .false. -namelist /trans_mitdart_nml/ go_to_dart, do_bgc, log_transform +namelist /trans_mitdart_nml/ go_to_dart, do_bgc, log_transform, compress ! /pkg/mdsio/mdsio_write_meta.F writes the .meta files type MIT_meta_type @@ -327,6 +337,7 @@ subroutine static_init_model() integer :: i, iunit, io integer :: ss, dd +integer :: ncid ! for reading compressed coordinates ! The Plan: ! @@ -528,13 +539,62 @@ subroutine static_init_model() domain_id = add_domain(model_shape_file, nvars, & var_names, quantity_list, clamp_vals, update_list ) +if (compress) then ! read in compressed coordinates + + ncid = nc_open_file_readonly(model_shape_file) + comp2d = nc_get_dimension_size(ncid, 'comp2d' , 'static_init_model', model_shape_file) + comp3d = nc_get_dimension_size(ncid, 'comp3d' , 'static_init_model', model_shape_file) + comp3dU = nc_get_dimension_size(ncid, 'comp3dU', 'static_init_model', model_shape_file) + comp3dV = nc_get_dimension_size(ncid, 'comp3dV', 'static_init_model', model_shape_file) + + allocate(XC_sq(comp3d)) + allocate(YC_sq(comp3d)) + allocate(ZC_sq(comp3d)) ! ZC is r8 + + allocate(XG_sq(comp3d)) + allocate(YG_sq(comp3d)) + + allocate(Xc_Ti(comp3d)) + allocate(Yc_Ti(comp3d)) + allocate(Zc_Ti(comp3d)) + + allocate(Xc_Ui(comp3dU)) + allocate(Yc_Ui(comp3dU)) + allocate(Zc_Ui(comp3dU)) + + allocate(Xc_Vi(comp3dV)) + allocate(Yc_Vi(comp3dV)) + allocate(Zc_Vi(comp3dV)) + + call nc_get_variable(ncid, 'XCcomp', XC_sq) + call nc_get_variable(ncid, 'YCcomp', YC_sq) + call nc_get_variable(ncid, 'ZCcomp', ZC_sq) + + call nc_get_variable(ncid, 'XGcomp', XG_sq) + call nc_get_variable(ncid, 'YGcomp', YG_sq) + + call nc_get_variable(ncid, 'Xcomp_ind', Xc_Ti) + call nc_get_variable(ncid, 'Ycomp_ind', Yc_Ti) + call nc_get_variable(ncid, 'Zcomp_ind', Zc_Ti) + + call nc_get_variable(ncid, 'Xcomp_indU', Xc_Ui) + call nc_get_variable(ncid, 'Ycomp_indU', Yc_Ui) + call nc_get_variable(ncid, 'Zcomp_indU', Zc_Ui) + + call nc_get_variable(ncid, 'Xcomp_indV', Xc_Vi) + call nc_get_variable(ncid, 'Ycomp_indV', Yc_Vi) + call nc_get_variable(ncid, 'Zcomp_indV', Zc_Vi) + + call nc_close_file(ncid) + +endif + model_size = get_domain_size(domain_id) if (do_output()) write(*,*) 'model_size = ', model_size end subroutine static_init_model - function get_model_size() !------------------------------------------------------------------ ! @@ -954,6 +1014,63 @@ function lon_dist(lon1, lon2) end function lon_dist +function get_compressed_dart_vector_index(iloc, jloc, kloc, dom_id, var_id) +!======================================================================= +! + +! returns the dart vector index for the compressed state + +integer, intent(in) :: iloc, jloc, kloc +integer, intent(in) :: dom_id, var_id +integer(i8) :: get_compressed_dart_vector_index + +integer :: i ! loop counter +integer :: qty +integer(i8) :: offset + +offset = get_index_start(dom_id, var_id) + +qty = get_kind_index(dom_id, var_id) + +get_compressed_dart_vector_index = -1 + +! MEG: Using the already established compressed indices +! +! 2D compressed variables +if (qty == QTY_SEA_SURFACE_HEIGHT .or. qty == QTY_SURFACE_CHLOROPHYLL ) then + do i = 1, comp2d + if (Xc_Ti(i) == iloc .and. Yc_Ti(i) == jloc .and. Zc_Ti(i) == 1) then + get_compressed_dart_vector_index = offset + i - 1 + endif + enddo + return +endif + +! 3D compressed variables +if (qty == QTY_U_CURRENT_COMPONENT) then + do i = 1, comp3dU + if (Xc_Ui(i) == iloc .and. Yc_Ui(i) == jloc .and. Zc_Ui(i) == kloc) then + get_compressed_dart_vector_index = offset + i - 1 + endif + enddo +elseif (qty == QTY_V_CURRENT_COMPONENT) then + do i = 1, comp3dV + if (Xc_Vi(i) == iloc .and. Yc_Vi(i) == jloc .and. Zc_Vi(i) == kloc) then + get_compressed_dart_vector_index = offset + i - 1 + endif + enddo +else + do i = 1, comp3d + if (Xc_Ti(i) == iloc .and. Yc_Ti(i) == jloc .and. Zc_Ti(i) == kloc) then + get_compressed_dart_vector_index = offset + i - 1 + endif + enddo +endif + + +end function get_compressed_dart_vector_index + + function get_val(lon_index, lat_index, level, var_id, state_handle,ens_size, masked) !======================================================================= ! @@ -971,24 +1088,28 @@ function get_val(lon_index, lat_index, level, var_id, state_handle,ens_size, mas if ( .not. module_initialized ) call static_init_model -state_index = get_dart_vector_index(lon_index, lat_index, level, domain_id, var_id) -get_val = get_state(state_index,state_handle) +masked = .false. -! Masked returns false if the value is masked -! A grid variable is assumed to be masked if its value is FVAL. -! Just to maintain legacy, we also assume that A grid variable is assumed -! to be masked if its value is exactly 0. -! See discussion in lat_lon_interpolate. +if (compress) then -! MEG CAUTION: THE ABOVE STATEMENT IS INCORRECT -! trans_mitdart already looks for 0.0 and makes them FVAL -! So, in the condition below we don't need to check for zeros -! The only mask is FVAL -masked = .false. -do i=1,ens_size -! if(get_val(i) == FVAL .or. get_val(i) == 0.0_r8 ) masked = .true. - if(get_val(i) == FVAL) masked = .true. -enddo + state_index = get_compressed_dart_vector_index(lon_index, lat_index, level, domain_id, var_id) + + if (state_index .ne. -1) then + get_val = get_state(state_index,state_handle) + else + masked = .true. + endif + +else + + state_index = get_dart_vector_index(lon_index, lat_index, level, domain_id, var_id) + get_val = get_state(state_index,state_handle) + + do i=1,ens_size ! HK this is checking the whole ensemble, can you have different masks for each ensemble member? + if(get_val(i) == FVAL) masked = .true. + enddo + +endif end function get_val @@ -1079,16 +1200,28 @@ subroutine get_state_meta_data(index_in, location, qty) call get_model_variable_indices(index_in, iloc, jloc, kloc, kind_index = qty) -lon = XC(iloc) -lat = YC(jloc) -depth = ZC(kloc) +if (compress) then ! all variables ae 1D + lon = XC_sq(iloc) + lat = YC_sq(iloc) + depth = ZC_sq(iloc) + ! Acounting for variables those on staggered grids + if (qty == QTY_U_CURRENT_COMPONENT) lon = XG_sq(iloc) + if (qty == QTY_V_CURRENT_COMPONENT) lat = YG_sq(iloc) +else + + lon = XC(iloc) + lat = YC(jloc) + depth = ZC(kloc) + + ! Acounting for variables those on staggered grids + if (qty == QTY_U_CURRENT_COMPONENT) lon = XG(iloc) + if (qty == QTY_V_CURRENT_COMPONENT) lat = YG(jloc) + +endif -! Acounting for surface variables and those on staggered grids ! MEG: check chl's depth here if (qty == QTY_SEA_SURFACE_HEIGHT .or. & qty == QTY_SURFACE_CHLOROPHYLL) depth = 0.0_r8 -if (qty == QTY_U_CURRENT_COMPONENT) lon = XG(iloc) -if (qty == QTY_V_CURRENT_COMPONENT) lat = YG(jloc) location = set_location(lon, lat, depth, VERTISHEIGHT) @@ -1297,6 +1430,8 @@ end subroutine nc_write_model_atts !------------------------------------------------------------------ ! Create an ensemble of states from a single state. +! Note if you perturb a compressed state, this will not be bitwise +! with perturbing a non-compressed state. subroutine pert_model_copies(state_ens_handle, ens_size, pert_amp, interf_provided) type(ensemble_type), intent(inout) :: state_ens_handle diff --git a/models/MITgcm_ocean/model_mod.nml b/models/MITgcm_ocean/model_mod.nml index ef64b88caa..03b81505ae 100644 --- a/models/MITgcm_ocean/model_mod.nml +++ b/models/MITgcm_ocean/model_mod.nml @@ -2,5 +2,6 @@ assimilation_period_days = 7 assimilation_period_seconds = 0 model_perturbation_amplitude = 0.2 + model_shape_file = 'mem01_reduced.nc' / diff --git a/models/MITgcm_ocean/readme.rst b/models/MITgcm_ocean/readme.rst index d63e56fea0..c6b7dc3dfd 100644 --- a/models/MITgcm_ocean/readme.rst +++ b/models/MITgcm_ocean/readme.rst @@ -34,8 +34,14 @@ can be set in the ``&trans_mitdart_nml`` namelist in ``input.nml``. &trans_mitdart_nml do_bgc = .false. ! change to .true. if doing bio-geo-chemistry log_transform = .false. ! change to .true. if using log_transform + compress = .false. ! change to .true. to compress the state vector / +``compress = .true.`` can be used to generate netcdf files for use with DART which has missing values (land) removed. +For some datasets this reduces the state vector size significantly. For example, the state vector size is +reduced by approximately 90% for the Red Sea. The program ``expand_netcdf`` can be used to uncompress the netcdf +file to view the data in a convenient form. + .. Warning:: diff --git a/models/MITgcm_ocean/trans_mitdart_mod.f90 b/models/MITgcm_ocean/trans_mitdart_mod.f90 index 15861b8164..8ce45307df 100644 --- a/models/MITgcm_ocean/trans_mitdart_mod.f90 +++ b/models/MITgcm_ocean/trans_mitdart_mod.f90 @@ -9,6 +9,7 @@ module trans_mitdart_mod use utilities_mod, only: initialize_utilities, register_module, & get_unit, find_namelist_in_file, file_exist, & check_namelist_read +use netcdf_utilities_mod, only : nc_get_variable, nc_get_dimension_size use netcdf implicit none @@ -20,9 +21,16 @@ module trans_mitdart_mod integer :: io, iunit logical :: do_bgc = .false. -logical :: log_transform = .false. +logical :: log_transform = .false. +logical :: compress = .false. +! set compress = .true. remove missing values from state +logical :: output_chl_data = .false. +! CHL.data is not written to mit .data files by default -namelist /trans_mitdart_nml/ do_bgc, log_transform +namelist /trans_mitdart_nml/ do_bgc, log_transform, compress + +real(r4), parameter :: FVAL=-999.0_r4 ! may put this as a namelist option +real(r4), parameter :: binary_fill=0.0_r4 !------------------------------------------------------------------ ! @@ -43,7 +51,7 @@ module trans_mitdart_mod integer :: recl3d integer :: recl2d -!-- Gridding parameters variable declarations +!-- Gridding parameters variable declarations logical :: usingCartesianGrid, usingCylindricalGrid, & usingSphericalPolarGrid, usingCurvilinearGrid, & deepAtmosphere @@ -71,14 +79,49 @@ module trans_mitdart_mod ! standard MITgcm namelist and filled in here. integer :: Nx=-1, Ny=-1, Nz=-1 ! grid counts for each field +integer :: ncomp2 = -1 ! length of 2D compressed dim +integer :: ncomp3 = -1, ncomp3U = -1, ncomp3V = -1 ! length of 3D compressed dim + +integer, parameter :: MITgcm_3D_FIELD = 1 +integer, parameter :: MITgcm_3D_FIELD_U = 2 +integer, parameter :: MITgcm_3D_FIELD_V = 3 ! locations of cell centers (C) and edges (G) for each axis. real(r8), allocatable :: XC(:), XG(:), YC(:), YG(:), ZC(:), ZG(:) +real(r8), allocatable :: XCcomp(:), XGcomp(:), YCcomp(:), YGcomp(:), ZCcomp(:), ZGcomp(:) + +integer, allocatable :: Xcomp_ind(:), Ycomp_ind(:), Zcomp_ind(:) !HK are the staggered grids compressed the same? +!MEG: For staggered grids +integer, allocatable :: Xcomp_indU(:), Ycomp_indU(:), Zcomp_indU(:) +integer, allocatable :: Xcomp_indV(:), Ycomp_indV(:), Zcomp_indV(:) + +! 3D variables, 3 grids: +! +! XC, YC, ZC 1 PSAL, PTMP, NO3, PO4, O2, PHY, ALK, DIC, DOP, DON, FET +! XC, YC, ZG 2 UVEL +! XC, YG, ZC 3 VVEL + +! MEG: For compression, especially if we're doing Arakawa C-grid, +! we will need 3 different compressions for the above variables + +! 2D variables, 1 grid: +! +! YC, XC ETA, CHL private public :: static_init_trans, mit2dart, dart2mit +interface write_compressed + module procedure write_compressed_2d + module procedure write_compressed_3d +end interface write_compressed + +interface read_compressed + module procedure read_compressed_2d + module procedure read_compressed_3d +end interface read_compressed + contains !================================================================== @@ -100,7 +143,6 @@ subroutine static_init_trans() read(iunit, nml = trans_mitdart_nml, iostat = io) call check_namelist_read(iunit, io, 'trans_mitdart_nml') - ! Grid-related variables are in PARM04 delX(:) = 0.0_r4 delY(:) = 0.0_r4 @@ -196,13 +238,6 @@ subroutine static_init_trans() recl3d = Nx*Ny*Nz*4 recl2d = Nx*Ny*4 -! MEG Better have that as inout namelist parameter -! Are we also doing bgc on top of physics? -! If we found nitrate then the rest of the binaries (for the -! remaining 9 variables) should be also there. -! TODO may also enhance this functionality -! if (file_exist('NO3.data')) do_bgc = .true. - end subroutine static_init_trans !------------------------------------------------------------------ @@ -210,11 +245,17 @@ end subroutine static_init_trans subroutine mit2dart() -integer :: ncid, iunit +integer :: ncid ! for the dimensions and coordinate variables integer :: XGDimID, XCDimID, YGDimID, YCDimID, ZGDimID, ZCDimID integer :: XGVarID, XCVarID, YGVarID, YCVarID, ZGVarID, ZCVarID +integer :: comp2ID, comp3ID, comp3UD, comp3VD ! compressed dim +integer :: XGcompVarID, XCcompVarID, YGcompVarID, YCcompVarID, ZGcompVarID, ZCcompVarID +integer :: XindID, YindID, ZindID +integer :: XindUD, YindUD, ZindUD +integer :: XindVD, YindVD, ZindVD +integer :: all_dimids(9) ! store the 9 dimension ids that are used ! for the prognostic variables integer :: SVarID, TVarID, UVarID, VVarID, EtaVarID @@ -224,27 +265,47 @@ subroutine mit2dart() ! diagnostic variable integer :: chl_varid -real(r4), allocatable :: data_3d(:,:,:), data_2d(:,:) - -real(r4) :: FVAL - if (.not. module_initialized) call static_init_trans -FVAL=-999.0_r4 - -allocate(data_3d(Nx,Ny,Nz)) -allocate(data_2d(Nx,Ny)) - call check(nf90_create(path="OUTPUT.nc",cmode=or(nf90_clobber,nf90_64bit_offset),ncid=ncid)) ! Define the new dimensions IDs - -call check(nf90_def_dim(ncid=ncid, name="XG", len = Nx, dimid = XGDimID)) + call check(nf90_def_dim(ncid=ncid, name="XC", len = Nx, dimid = XCDimID)) -call check(nf90_def_dim(ncid=ncid, name="YG", len = Ny, dimid = YGDimID)) call check(nf90_def_dim(ncid=ncid, name="YC", len = Ny, dimid = YCDimID)) call check(nf90_def_dim(ncid=ncid, name="ZC", len = Nz, dimid = ZCDimID)) - + +call check(nf90_def_dim(ncid=ncid, name="XG", len = Nx, dimid = XGDimID)) +call check(nf90_def_dim(ncid=ncid, name="YG", len = Ny, dimid = YGDimID)) + +print *, '' + +if (compress) then + ncomp2 = get_compressed_size_2d() + + write(*, '(A, I12, A, I8)') '2D: ', Nx*Ny, ', COMP2D: ', ncomp2 + + ncomp3 = get_compressed_size_3d(MITgcm_3D_FIELD) + ncomp3U = get_compressed_size_3d(MITgcm_3D_FIELD_U) + ncomp3V = get_compressed_size_3d(MITgcm_3D_FIELD_V) + + write(*, '(A, I12, A, 3I8)') '3D: ', Nx*Ny*Nz, ', COMP3D [T-S, U, V]: ', ncomp3, ncomp3U, ncomp3V + + ! Put the compressed dimensions in the restart file + call check(nf90_def_dim(ncid=ncid, name="comp2d", len = ncomp2, dimid = comp2ID)) + call check(nf90_def_dim(ncid=ncid, name="comp3d", len = ncomp3, dimid = comp3ID)) + call check(nf90_def_dim(ncid=ncid, name="comp3dU", len = ncomp3U, dimid = comp3UD)) + call check(nf90_def_dim(ncid=ncid, name="comp3dV", len = ncomp3V, dimid = comp3VD)) +else + comp2ID = -1 + comp3ID = -1 +endif + +all_dimids = (/XCDimID, YCDimID, ZCDimID, XGDimID, YGDimID, & + comp2ID, comp3ID, comp3UD, comp3VD/) + +print *, '' + ! Create the (empty) Coordinate Variables and the Attributes ! U Grid Longitudes @@ -290,142 +351,89 @@ subroutine mit2dart() call check(nf90_put_att(ncid, ZCVarID, "axis", "Z")) call check(nf90_put_att(ncid, ZCVarID, "standard_name", "depth")) +! Compressed grid variables +if (compress) then + call check(nf90_def_var(ncid,name="XGcomp",xtype=nf90_real,dimids=comp3ID,varid=XGcompVarID)) + call check(nf90_def_var(ncid,name="XCcomp",xtype=nf90_real,dimids=comp3ID,varid=XCcompVarID)) + call check(nf90_def_var(ncid,name="YGcomp",xtype=nf90_real,dimids=comp3ID,varid=YGcompVarID)) + call check(nf90_def_var(ncid,name="YCcomp",xtype=nf90_real,dimids=comp3ID,varid=YCcompVarID)) + call check(nf90_def_var(ncid,name="ZCcomp",xtype=nf90_double,dimids=comp3ID,varid=ZCcompVarID)) + + call check(nf90_def_var(ncid,name="Xcomp_ind",xtype=nf90_int,dimids=comp3ID,varid=XindID)) + call check(nf90_def_var(ncid,name="Ycomp_ind",xtype=nf90_int,dimids=comp3ID,varid=YindID)) + call check(nf90_def_var(ncid,name="Zcomp_ind",xtype=nf90_int,dimids=comp3ID,varid=ZindID)) + + call check(nf90_def_var(ncid,name="Xcomp_indU",xtype=nf90_int,dimids=comp3UD,varid=XindUD)) + call check(nf90_def_var(ncid,name="Ycomp_indU",xtype=nf90_int,dimids=comp3UD,varid=YindUD)) + call check(nf90_def_var(ncid,name="Zcomp_indU",xtype=nf90_int,dimids=comp3UD,varid=ZindUD)) + + call check(nf90_def_var(ncid,name="Xcomp_indV",xtype=nf90_int,dimids=comp3VD,varid=XindVD)) + call check(nf90_def_var(ncid,name="Ycomp_indV",xtype=nf90_int,dimids=comp3VD,varid=YindVD)) + call check(nf90_def_var(ncid,name="Zcomp_indV",xtype=nf90_int,dimids=comp3VD,varid=ZindVD)) +endif + +! The size of these variables will depend on the compression ! Create the (empty) Prognostic Variables and the Attributes -call check(nf90_def_var(ncid=ncid, name="PSAL", xtype=nf90_real, & - dimids = (/XCDimID,YCDimID,ZCDimID/),varid=SVarID)) -call check(nf90_put_att(ncid, SVarID, "long_name", "potential salinity")) -call check(nf90_put_att(ncid, SVarID, "missing_value", FVAL)) -call check(nf90_put_att(ncid, SVarID, "_FillValue", FVAL)) -call check(nf90_put_att(ncid, SVarID, "units", "psu")) -call check(nf90_put_att(ncid, SVarID, "units_long_name", "practical salinity units")) - -call check(nf90_def_var(ncid=ncid, name="PTMP", xtype=nf90_real, & - dimids=(/XCDimID,YCDimID,ZCDimID/),varid=TVarID)) -call check(nf90_put_att(ncid, TVarID, "long_name", "Potential Temperature")) -call check(nf90_put_att(ncid, TVarID, "missing_value", FVAL)) -call check(nf90_put_att(ncid, TVarID, "_FillValue", FVAL)) -call check(nf90_put_att(ncid, TVarID, "units", "C")) -call check(nf90_put_att(ncid, TVarID, "units_long_name", "degrees celsius")) - -call check(nf90_def_var(ncid=ncid, name="UVEL", xtype=nf90_real, & - dimids=(/XGDimID,YCDimID,ZCDimID/),varid=UVarID)) -call check(nf90_put_att(ncid, UVarID, "long_name", "Zonal Velocity")) -call check(nf90_put_att(ncid, UVarID, "mssing_value", FVAL)) -call check(nf90_put_att(ncid, UVarID, "_FillValue", FVAL)) -call check(nf90_put_att(ncid, UVarID, "units", "m/s")) -call check(nf90_put_att(ncid, UVarID, "units_long_name", "meters per second")) - -call check(nf90_def_var(ncid=ncid, name="VVEL", xtype=nf90_real, & - dimids=(/XCDimID,YGDimID,ZCDimID/),varid=VVarID)) -call check(nf90_put_att(ncid, VVarID, "long_name", "Meridional Velocity")) -call check(nf90_put_att(ncid, VVarID, "missing_value", FVAL)) -call check(nf90_put_att(ncid, VVarID, "_FillValue", FVAL)) -call check(nf90_put_att(ncid, VVarID, "units", "m/s")) -call check(nf90_put_att(ncid, VVarID, "units_long_name", "meters per second")) - -call check(nf90_def_var(ncid=ncid, name="ETA", xtype=nf90_real, & - dimids=(/XCDimID,YCDimID/),varid=EtaVarID)) -call check(nf90_put_att(ncid, EtaVarID, "long_name", "sea surface height")) -call check(nf90_put_att(ncid, EtaVarID, "missing_value", FVAL)) -call check(nf90_put_att(ncid, EtaVarID, "_FillValue", FVAL)) -call check(nf90_put_att(ncid, EtaVarID, "units", "m")) -call check(nf90_put_att(ncid, EtaVarID, "units_long_name", "meters")) - -!> Add BLING data: +SVarID = define_variable(ncid,"PSAL", nf90_real, all_dimids, MITgcm_3D_FIELD) +call add_attributes_to_variable(ncid, SVarID, "potential salinity", "psu", "practical salinity units") + +TVarID = define_variable(ncid,"PTMP", nf90_real, all_dimids, MITgcm_3D_FIELD) +call add_attributes_to_variable(ncid, TVarID, "Potential Temperature", "C", "degrees celsius") + +UVarID = define_variable(ncid,"UVEL", nf90_real, all_dimids, MITgcm_3D_FIELD_U) +call add_attributes_to_variable(ncid, UVarID, "Zonal Velocity", "m/s", "meters per second") + +VVarID = define_variable(ncid,"VVEL", nf90_real, all_dimids, MITgcm_3D_FIELD_V) +call add_attributes_to_variable(ncid, VVarID, "Meridional Velocity", "m/s", "meters per second") + +EtaVarID = define_variable_2d(ncid,"ETA", nf90_real, all_dimids) +call add_attributes_to_variable(ncid, EtaVarID, "sea surface height", "m", "meters") + +! Create the BLING netcdf variables: if (do_bgc) then ! 1. BLING tracer: nitrate NO3 - call check(nf90_def_var(ncid=ncid, name="NO3", xtype=nf90_real, & - dimids=(/XCDimID,YCDimID,ZCDimID/),varid=no3_varid)) - call check(nf90_put_att(ncid, no3_varid, "long_name" , "Nitrate")) - call check(nf90_put_att(ncid, no3_varid, "missing_value" , FVAL)) - call check(nf90_put_att(ncid, no3_varid, "_FillValue" , FVAL)) - call check(nf90_put_att(ncid, no3_varid, "units" , "mol N/m3")) - call check(nf90_put_att(ncid, no3_varid, "units_long_name", "moles Nitrogen per cubic meters")) - + no3_varid = define_variable(ncid,"NO3", nf90_real, all_dimids, MITgcm_3D_FIELD) + call add_attributes_to_variable(ncid, no3_varid, "Nitrate", "mol N/m3", "moles Nitrogen per cubic meters") + ! 2. BLING tracer: phosphate PO4 - call check(nf90_def_var(ncid=ncid, name="PO4", xtype=nf90_real, & - dimids=(/XCDimID,YCDimID,ZCDimID/),varid=po4_varid)) - call check(nf90_put_att(ncid, po4_varid, "long_name" , "Phosphate")) - call check(nf90_put_att(ncid, po4_varid, "missing_value" , FVAL)) - call check(nf90_put_att(ncid, po4_varid, "_FillValue" , FVAL)) - call check(nf90_put_att(ncid, po4_varid, "units" , "mol P/m3")) - call check(nf90_put_att(ncid, po4_varid, "units_long_name", "moles Phosphorus per cubic meters")) - + po4_varid = define_variable(ncid,"PO4", nf90_real, all_dimids, MITgcm_3D_FIELD) + call add_attributes_to_variable(ncid, po4_varid, "Phosphate", "mol P/m3", "moles Phosphorus per cubic meters") + ! 3. BLING tracer: oxygen O2 - call check(nf90_def_var(ncid=ncid, name="O2", xtype=nf90_real, & - dimids=(/XCDimID,YCDimID,ZCDimID/),varid=o2_varid)) - call check(nf90_put_att(ncid, o2_varid, "long_name" , "Dissolved Oxygen")) - call check(nf90_put_att(ncid, o2_varid, "missing_value" , FVAL)) - call check(nf90_put_att(ncid, o2_varid, "_FillValue" , FVAL)) - call check(nf90_put_att(ncid, o2_varid, "units" , "mol O/m3")) - call check(nf90_put_att(ncid, o2_varid, "units_long_name", "moles Oxygen per cubic meters")) - + o2_varid = define_variable(ncid,"O2", nf90_real, all_dimids, MITgcm_3D_FIELD) + call add_attributes_to_variable(ncid, o2_varid, "Dissolved Oxygen", "mol O/m3", "moles Oxygen per cubic meters") + ! 4. BLING tracer: phytoplankton PHY - call check(nf90_def_var(ncid=ncid, name="PHY", xtype=nf90_real, & - dimids=(/XCDimID,YCDimID,ZCDimID/),varid=phy_varid)) - call check(nf90_put_att(ncid, phy_varid, "long_name" , "Phytoplankton Biomass")) - call check(nf90_put_att(ncid, phy_varid, "missing_value" , FVAL)) - call check(nf90_put_att(ncid, phy_varid, "_FillValue" , FVAL)) - call check(nf90_put_att(ncid, phy_varid, "units" , "mol C/m3")) - call check(nf90_put_att(ncid, phy_varid, "units_long_name", "moles Carbon per cubic meters")) + phy_varid = define_variable(ncid,"PHY", nf90_real, all_dimids, MITgcm_3D_FIELD) + call add_attributes_to_variable(ncid, phy_varid, "Phytoplankton Biomass", "mol C/m3", "moles Carbon per cubic meters") ! 5. BLING tracer: alkalinity ALK - call check(nf90_def_var(ncid=ncid, name="ALK", xtype=nf90_real, & - dimids=(/XCDimID,YCDimID,ZCDimID/),varid=alk_varid)) - call check(nf90_put_att(ncid, alk_varid, "long_name" , "Alkalinity")) - call check(nf90_put_att(ncid, alk_varid, "missing_value" , FVAL)) - call check(nf90_put_att(ncid, alk_varid, "_FillValue" , FVAL)) - call check(nf90_put_att(ncid, alk_varid, "units" , "mol eq/m3")) - call check(nf90_put_att(ncid, alk_varid, "units_long_name", "moles equivalent per cubic meters")) - + alk_varid = define_variable(ncid,"ALK", nf90_real, all_dimids, MITgcm_3D_FIELD) + call add_attributes_to_variable(ncid, alk_varid, "Alkalinity", "mol eq/m3", "moles equivalent per cubic meters") + ! 6. BLING tracer: dissolved inorganic carbon DIC - call check(nf90_def_var(ncid=ncid, name="DIC", xtype=nf90_real, & - dimids=(/XCDimID,YCDimID,ZCDimID/),varid=dic_varid)) - call check(nf90_put_att(ncid, dic_varid, "long_name" , "Dissolved Inorganic Carbon")) - call check(nf90_put_att(ncid, dic_varid, "missing_value" , FVAL)) - call check(nf90_put_att(ncid, dic_varid, "_FillValue" , FVAL)) - call check(nf90_put_att(ncid, dic_varid, "units" , "mol C/m3")) - call check(nf90_put_att(ncid, dic_varid, "units_long_name", "moles Carbon per cubic meters")) - - ! 7. BLING tracer: dissolved organic phosphorus DOP - call check(nf90_def_var(ncid=ncid, name="DOP", xtype=nf90_real, & - dimids=(/XCDimID,YCDimID,ZCDimID/),varid=dop_varid)) - call check(nf90_put_att(ncid, dop_varid, "long_name" , "Dissolved Organic Phosphorus")) - call check(nf90_put_att(ncid, dop_varid, "missing_value" , FVAL)) - call check(nf90_put_att(ncid, dop_varid, "_FillValue" , FVAL)) - call check(nf90_put_att(ncid, dop_varid, "units" , "mol P/m3")) - call check(nf90_put_att(ncid, dop_varid, "units_long_name", "moles Phosphorus per cubic meters")) + dic_varid = define_variable(ncid,"DIC", nf90_real, all_dimids, MITgcm_3D_FIELD) + call add_attributes_to_variable(ncid, dic_varid, "Dissolved Inorganic Carbon", "mol C/m3", "moles Carbon per cubic meters") + + ! 7. BLING tracer: dissolved organic phosphorus DOP + dop_varid = define_variable(ncid,"DOP", nf90_real, all_dimids, MITgcm_3D_FIELD) + call add_attributes_to_variable(ncid, dop_varid, "Dissolved Organic Phosphorus", "mol P/m3", "moles Phosphorus per cubic meters") ! 8. BLING tracer: dissolved organic nitrogen DON - call check(nf90_def_var(ncid=ncid, name="DON", xtype=nf90_real, & - dimids=(/XCDimID,YCDimID,ZCDimID/),varid=don_varid)) - call check(nf90_put_att(ncid, don_varid, "long_name" , "Dissolved Organic Nitrogen")) - call check(nf90_put_att(ncid, don_varid, "missing_value" , FVAL)) - call check(nf90_put_att(ncid, don_varid, "_FillValue" , FVAL)) - call check(nf90_put_att(ncid, don_varid, "units" , "mol N/m3")) - call check(nf90_put_att(ncid, don_varid, "units_long_name", "moles Nitrogen per cubic meters")) + don_varid = define_variable(ncid,"DON", nf90_real, all_dimids, MITgcm_3D_FIELD) + call add_attributes_to_variable(ncid, don_varid, "Dissolved Organic Nitrogen", "mol N/m3", "moles Nitrogen per cubic meters") ! 9. BLING tracer: dissolved inorganic iron FET - call check(nf90_def_var(ncid=ncid, name="FET", xtype=nf90_real, & - dimids=(/XCDimID,YCDimID,ZCDimID/),varid=fet_varid)) - call check(nf90_put_att(ncid, fet_varid, "long_name" , "Dissolved Inorganic Iron")) - call check(nf90_put_att(ncid, fet_varid, "missing_value" , FVAL)) - call check(nf90_put_att(ncid, fet_varid, "_FillValue" , FVAL)) - call check(nf90_put_att(ncid, fet_varid, "units" , "mol Fe/m3")) - call check(nf90_put_att(ncid, fet_varid, "units_long_name", "moles Iron per cubic meters")) - + fet_varid = define_variable(ncid,"FET", nf90_real, all_dimids, MITgcm_3D_FIELD) + call add_attributes_to_variable(ncid, fet_varid, "Dissolved Inorganic Iron", "mol Fe/m3", "moles Iron per cubic meters") + ! 10. BLING tracer: Surface Chlorophyl CHL - call check(nf90_def_var(ncid=ncid, name="CHL", xtype=nf90_real, & - dimids=(/XCDimID,YCDimID/),varid=chl_varid)) - call check(nf90_put_att(ncid, chl_varid, "long_name" , "Surface Chlorophyll")) - call check(nf90_put_att(ncid, chl_varid, "missing_value" , FVAL)) - call check(nf90_put_att(ncid, chl_varid, "_FillValue" , FVAL)) - call check(nf90_put_att(ncid, chl_varid, "units" , "mg/m3")) - call check(nf90_put_att(ncid, chl_varid, "units_long_name", "milligram per cubic meters")) -endif + chl_varid = define_variable_2d(ncid,"CHL", nf90_real, all_dimids) + call add_attributes_to_variable(ncid, chl_varid, "Surface Chlorophyll", "mg/m3", "milligram per cubic meters" ) +endif ! Finished with dimension/variable definitions, must end 'define' mode to fill. @@ -439,125 +447,72 @@ subroutine mit2dart() call check(nf90_put_var(ncid, YCVarID, YC )) call check(nf90_put_var(ncid, ZCVarID, ZC )) -! Fill the data - -iunit = get_unit() -open(iunit, file='PSAL.data', form='UNFORMATTED', status='OLD', & - access='DIRECT', recl=recl3d, convert='BIG_ENDIAN') -read(iunit,rec=1)data_3d -close(iunit) -where (data_3d == 0.0_r4) data_3d = FVAL -call check(nf90_put_var(ncid,SVarID,data_3d,start=(/1,1,1/))) - -open(iunit, file='PTMP.data', form='UNFORMATTED', status='OLD', & - access='DIRECT', recl=recl3d, convert='BIG_ENDIAN') -read(iunit,rec=1)data_3d -close(iunit) -where (data_3d == 0.0_r4) data_3d = FVAL -call check(nf90_put_var(ncid,TVarID,data_3d,start=(/1,1,1/))) - -open(iunit, file='UVEL.data', form='UNFORMATTED', status='OLD', & - access='DIRECT', recl=recl3d, convert='BIG_ENDIAN') -read(iunit,rec=1)data_3d -close(iunit) -where (data_3d == 0.0_r4) data_3d = FVAL -call check(nf90_put_var(ncid,UVarID,data_3d,start=(/1,1,1/))) - -open(iunit, file='VVEL.data', form='UNFORMATTED', status='OLD', & - access='DIRECT', recl=recl3d, convert='BIG_ENDIAN') -read(iunit,rec=1)data_3d -close(iunit) -where (data_3d == 0.0_r4) data_3d = FVAL -call check(nf90_put_var(ncid,VVarID,data_3d,start=(/1,1,1/))) - -open(iunit, file='ETA.data', form='UNFORMATTED', status='OLD', & - access='DIRECT', recl=recl2d, convert='BIG_ENDIAN') -read(iunit,rec=1)data_2d -close(iunit) -where (data_2d == 0.0_r4) data_2d = FVAL -call check(nf90_put_var(ncid,EtaVarID,data_2d,start=(/1,1/))) +if (compress) then + allocate(XCcomp(ncomp3)) + allocate(XGcomp(ncomp3)) + allocate(YCcomp(ncomp3)) + allocate(YGcomp(ncomp3)) + allocate(ZCcomp(ncomp3)) + allocate(ZGcomp(ncomp3)) + allocate(Xcomp_ind(ncomp3)) + allocate(Ycomp_ind(ncomp3)) + allocate(Zcomp_ind(ncomp3)) + + allocate(Xcomp_indU(ncomp3U)) + allocate(Ycomp_indU(ncomp3U)) + allocate(Zcomp_indU(ncomp3U)) + + allocate(Xcomp_indV(ncomp3V)) + allocate(Ycomp_indV(ncomp3V)) + allocate(Zcomp_indV(ncomp3V)) + + call fill_compressed_coords() + + call check(nf90_put_var(ncid, XGcompVarID, XGcomp )) + call check(nf90_put_var(ncid, XCcompVarID, XCcomp )) + call check(nf90_put_var(ncid, YGcompVarID, YGcomp )) + call check(nf90_put_var(ncid, YCcompVarID, YCcomp )) + call check(nf90_put_var(ncid, ZCcompVarID, ZCcomp )) + + call check(nf90_put_var(ncid, XindID, Xcomp_ind )) + call check(nf90_put_var(ncid, YindID, Ycomp_ind )) + call check(nf90_put_var(ncid, ZindID, Zcomp_ind )) + + call check(nf90_put_var(ncid, XindUD, Xcomp_indU )) + call check(nf90_put_var(ncid, YindUD, Ycomp_indU )) + call check(nf90_put_var(ncid, ZindUD, Zcomp_indU )) + + call check(nf90_put_var(ncid, XindVD, Xcomp_indV )) + call check(nf90_put_var(ncid, YindVD, Ycomp_indV )) + call check(nf90_put_var(ncid, ZindVD, Zcomp_indV )) +endif -if (do_bgc) then - open(iunit, file='NO3.data', form='UNFORMATTED', status='OLD', & - access='DIRECT', recl=recl3d, convert='BIG_ENDIAN') - read(iunit,rec=1)data_3d - close(iunit) - call fill_var_md(data_3d, FVAL) - call check(nf90_put_var(ncid,no3_varid,data_3d,start=(/1,1,1/))) - - open(iunit, file='PO4.data', form='UNFORMATTED', status='OLD', & - access='DIRECT', recl=recl3d, convert='BIG_ENDIAN') - read(iunit,rec=1)data_3d - close(iunit) - call fill_var_md(data_3d, FVAL) - call check(nf90_put_var(ncid,po4_varid,data_3d,start=(/1,1,1/))) - - open(iunit, file='O2.data', form='UNFORMATTED', status='OLD', & - access='DIRECT', recl=recl3d, convert='BIG_ENDIAN') - read(iunit,rec=1)data_3d - close(iunit) - call fill_var_md(data_3d, FVAL) - call check(nf90_put_var(ncid,o2_varid,data_3d,start=(/1,1,1/))) - - open(iunit, file='PHY.data', form='UNFORMATTED', status='OLD', & - access='DIRECT', recl=recl3d, convert='BIG_ENDIAN') - read(iunit,rec=1)data_3d - close(iunit) - call fill_var_md(data_3d, FVAL) - call check(nf90_put_var(ncid,phy_varid,data_3d,start=(/1,1,1/))) - - open(iunit, file='ALK.data', form='UNFORMATTED', status='OLD', & - access='DIRECT', recl=recl3d, convert='BIG_ENDIAN') - read(iunit,rec=1)data_3d - close(iunit) - call fill_var_md(data_3d, FVAL) - call check(nf90_put_var(ncid,alk_varid,data_3d,start=(/1,1,1/))) - - open(iunit, file='DIC.data', form='UNFORMATTED', status='OLD', & - access='DIRECT', recl=recl3d, convert='BIG_ENDIAN') - read(iunit,rec=1)data_3d - close(iunit) - call fill_var_md(data_3d, FVAL) - call check(nf90_put_var(ncid,dic_varid,data_3d,start=(/1,1,1/))) - - open(iunit, file='DOP.data', form='UNFORMATTED', status='OLD', & - access='DIRECT', recl=recl3d, convert='BIG_ENDIAN') - read(iunit,rec=1)data_3d - close(iunit) - call fill_var_md(data_3d, FVAL) - call check(nf90_put_var(ncid,dop_varid,data_3d,start=(/1,1,1/))) - - open(iunit, file='DON.data', form='UNFORMATTED', status='OLD', & - access='DIRECT', recl=recl3d, convert='BIG_ENDIAN') - read(iunit,rec=1)data_3d - close(iunit) - call fill_var_md(data_3d, FVAL) - call check(nf90_put_var(ncid,don_varid,data_3d,start=(/1,1,1/))) - - open(iunit, file='FET.data', form='UNFORMATTED', status='OLD', & - access='DIRECT', recl=recl3d, convert='BIG_ENDIAN') - read(iunit,rec=1)data_3d - close(iunit) - call fill_var_md(data_3d, FVAL) - call check(nf90_put_var(ncid,fet_varid,data_3d,start=(/1,1,1/))) - - open(iunit, file='CHL.data', form='UNFORMATTED', status='OLD', & - access='DIRECT', recl=recl2d, convert='BIG_ENDIAN') - read(iunit,rec=1)data_2d - close(iunit) - where (data_2d == 0.0_r4) - data_2d = FVAL - elsewhere - data_2d = log10(data_2d) - endwhere - call check(nf90_put_var(ncid,chl_varid,data_2d,start=(/1,1/))) +! Fill the netcdf variables +call from_mit_to_netcdf_3d('PSAL.data', ncid, SVarID, MITgcm_3D_FIELD) +call from_mit_to_netcdf_3d('PTMP.data', ncid, TVarID, MITgcm_3D_FIELD) +call from_mit_to_netcdf_3d('UVEL.data', ncid, UVarID, MITgcm_3D_FIELD_U) +call from_mit_to_netcdf_3d('VVEL.data', ncid, VVarID, MITgcm_3D_FIELD_V) +call from_mit_to_netcdf_2d('ETA.data' , ncid, EtaVarID) + +print *, 'Done writing physical variables' + +if (do_bgc) then + call from_mit_to_netcdf_tracer_3d('NO3.data', ncid, no3_varid) + call from_mit_to_netcdf_tracer_3d('PO4.data', ncid, po4_varid) + call from_mit_to_netcdf_tracer_3d('O2.data' , ncid, o2_varid) + call from_mit_to_netcdf_tracer_3d('PHY.data', ncid, phy_varid) + call from_mit_to_netcdf_tracer_3d('ALK.data', ncid, alk_varid) + call from_mit_to_netcdf_tracer_3d('DIC.data', ncid, dic_varid) + call from_mit_to_netcdf_tracer_3d('DOP.data', ncid, dop_varid) + call from_mit_to_netcdf_tracer_3d('DON.data', ncid, don_varid) + call from_mit_to_netcdf_tracer_3d('FET.data', ncid, fet_varid) + call from_mit_to_netcdf_tracer_2d('CHL.data', ncid, chl_varid) + + print *, 'Done writing biogeochemical variables' endif call check(nf90_close(ncid)) -deallocate(data_3d) -deallocate(data_2d) - end subroutine mit2dart !------------------------------------------------------------------ @@ -565,165 +520,74 @@ end subroutine mit2dart subroutine dart2mit() -integer :: ncid, varid, iunit -real(r4), allocatable :: data_3d(:,:,:),data_2d(:,:) -real(r4) :: FVAL - -allocate(data_3d(Nx,Ny,Nz)) -allocate(data_2d(Nx,Ny)) +integer :: ncid +recl2d = Nx*Ny*8 if (.not. module_initialized) call static_init_trans -iunit = get_unit() call check(nf90_open("INPUT.nc",NF90_NOWRITE,ncid)) -!Fill the data -call check( NF90_INQ_VARID(ncid,'PSAL',varid) ) -call check( NF90_GET_VAR(ncid,varid,data_3d)) -call check( nf90_get_att(ncid,varid,"_FillValue",FVAL)) -where (data_3d == FVAL) data_3d = 0.0_r4 - -open(iunit, file='PSAL.data', form="UNFORMATTED", status='UNKNOWN', & - access='DIRECT', recl=recl3d, convert='BIG_ENDIAN') -write(iunit,rec=1)data_3d -close(iunit) - -call check( NF90_INQ_VARID(ncid,'PTMP',varid) ) -call check( NF90_GET_VAR(ncid,varid,data_3d)) -call check( nf90_get_att(ncid,varid,"_FillValue",FVAL)) -where (data_3d == FVAL) data_3d = 0.0_r4 +if (compress) then + ncomp2 = nc_get_dimension_size(ncid,'comp2d') -open(iunit, file='PTMP.data', form="UNFORMATTED", status='UNKNOWN', & - access='DIRECT', recl=recl3d, convert='BIG_ENDIAN') -write(iunit,rec=1)data_3d -close(iunit) - -call check( NF90_INQ_VARID(ncid,'UVEL',varid) ) -call check( NF90_GET_VAR(ncid,varid,data_3d)) -call check( nf90_get_att(ncid,varid,"_FillValue",FVAL)) -where (data_3d == FVAL) data_3d = 0.0_r4 + ncomp3 = nc_get_dimension_size(ncid,'comp3d') + ncomp3U = nc_get_dimension_size(ncid,'comp3dU') + ncomp3V = nc_get_dimension_size(ncid,'comp3dV') -open(iunit, file='UVEL.data', form="UNFORMATTED", status='UNKNOWN', & - access='DIRECT', recl=recl3d, convert='BIG_ENDIAN') -write(iunit,rec=1)data_3d -close(iunit) + allocate(Xcomp_ind(ncomp3)) + allocate(Ycomp_ind(ncomp3)) + allocate(Zcomp_ind(ncomp3)) + + allocate(Xcomp_indU(ncomp3U)) + allocate(Ycomp_indU(ncomp3U)) + allocate(Zcomp_indU(ncomp3U)) -call check( NF90_INQ_VARID(ncid,'VVEL',varid) ) -call check( NF90_GET_VAR(ncid,varid,data_3d)) -call check( nf90_get_att(ncid,varid,"_FillValue",FVAL)) -where (data_3d == FVAL) data_3d = 0.0_r4 + allocate(Xcomp_indV(ncomp3V)) + allocate(Ycomp_indV(ncomp3V)) + allocate(Zcomp_indV(ncomp3V)) + + call nc_get_variable(ncid, 'Xcomp_ind', Xcomp_ind) + call nc_get_variable(ncid, 'Ycomp_ind', Ycomp_ind) + call nc_get_variable(ncid, 'Zcomp_ind', Zcomp_ind) -open(iunit, file='VVEL.data', form="UNFORMATTED", status='UNKNOWN', & - access='DIRECT', recl=recl3d, convert='BIG_ENDIAN') -write(iunit,rec=1)data_3d -close(iunit) + call nc_get_variable(ncid, 'Xcomp_indU', Xcomp_indU) + call nc_get_variable(ncid, 'Ycomp_indU', Ycomp_indU) + call nc_get_variable(ncid, 'Zcomp_indU', Zcomp_indU) -call check( NF90_INQ_VARID(ncid,'ETA',varid) ) -call check( NF90_GET_VAR(ncid,varid,data_2d)) -call check( nf90_get_att(ncid,varid,"_FillValue",FVAL)) -where (data_2d == FVAL) data_2d = 0.0_r4 + call nc_get_variable(ncid, 'Xcomp_indV', Xcomp_indV) + call nc_get_variable(ncid, 'Ycomp_indV', Ycomp_indV) + call nc_get_variable(ncid, 'Zcomp_indV', Zcomp_indV) +endif -open(iunit, file='ETA.data', form="UNFORMATTED", status='UNKNOWN', & +!Fill the data +iunit = get_unit() +open(iunit, file='PICKUP.OUTPUT', form="UNFORMATTED", status='UNKNOWN', & access='DIRECT', recl=recl2d, convert='BIG_ENDIAN') -write(iunit,rec=1)data_2d -close(iunit) -if (do_bgc) then - call check( NF90_INQ_VARID(ncid,'NO3',varid) ) - call check( NF90_GET_VAR(ncid,varid,data_3d)) - call check( nf90_get_att(ncid,varid,"_FillValue",FVAL)) - call fill_var_dm(data_3d, FVAL) - - open(iunit, file='NO3.data', form="UNFORMATTED", status='UNKNOWN', & - access='DIRECT', recl=recl3d, convert='BIG_ENDIAN') - write(iunit,rec=1)data_3d - close(iunit) - - call check( NF90_INQ_VARID(ncid,'PO4',varid) ) - call check( NF90_GET_VAR(ncid,varid,data_3d)) - call check( nf90_get_att(ncid,varid,"_FillValue",FVAL)) - call fill_var_dm(data_3d, FVAL) - - open(iunit, file='PO4.data', form="UNFORMATTED", status='UNKNOWN', & - access='DIRECT', recl=recl3d, convert='BIG_ENDIAN') - write(iunit,rec=1)data_3d - close(iunit) - - call check( NF90_INQ_VARID(ncid,'O2',varid) ) - call check( NF90_GET_VAR(ncid,varid,data_3d)) - call check( nf90_get_att(ncid,varid,"_FillValue",FVAL)) - call fill_var_dm(data_3d, FVAL) - - open(iunit, file='O2.data', form="UNFORMATTED", status='UNKNOWN', & - access='DIRECT', recl=recl3d, convert='BIG_ENDIAN') - write(iunit,rec=1)data_3d - close(iunit) - - call check( NF90_INQ_VARID(ncid,'PHY',varid) ) - call check( NF90_GET_VAR(ncid,varid,data_3d)) - call check( nf90_get_att(ncid,varid,"_FillValue",FVAL)) - call fill_var_dm(data_3d, FVAL) - - open(iunit, file='PHY.data', form="UNFORMATTED", status='UNKNOWN', & - access='DIRECT', recl=recl3d, convert='BIG_ENDIAN') - write(iunit,rec=1)data_3d - close(iunit) - - call check( NF90_INQ_VARID(ncid,'ALK',varid) ) - call check( NF90_GET_VAR(ncid,varid,data_3d)) - call check( nf90_get_att(ncid,varid,"_FillValue",FVAL)) - call fill_var_dm(data_3d, FVAL) - - open(iunit, file='ALK.data', form="UNFORMATTED", status='UNKNOWN', & - access='DIRECT', recl=recl3d, convert='BIG_ENDIAN') - write(iunit,rec=1)data_3d - close(iunit) - - call check( NF90_INQ_VARID(ncid,'DIC',varid) ) - call check( NF90_GET_VAR(ncid,varid,data_3d)) - call check( nf90_get_att(ncid,varid,"_FillValue",FVAL)) - call fill_var_dm(data_3d, FVAL) - - open(iunit, file='DIC.data', form="UNFORMATTED", status='UNKNOWN', & - access='DIRECT', recl=recl3d, convert='BIG_ENDIAN') - write(iunit,rec=1)data_3d - close(iunit) - - call check( NF90_INQ_VARID(ncid,'DOP',varid) ) - call check( NF90_GET_VAR(ncid,varid,data_3d)) - call check( nf90_get_att(ncid,varid,"_FillValue",FVAL)) - call fill_var_dm(data_3d, FVAL) - - open(iunit, file='DOP.data', form="UNFORMATTED", status='UNKNOWN', & - access='DIRECT', recl=recl3d, convert='BIG_ENDIAN') - write(iunit,rec=1)data_3d - close(iunit) - - call check( NF90_INQ_VARID(ncid,'DON',varid) ) - call check( NF90_GET_VAR(ncid,varid,data_3d)) - call check( nf90_get_att(ncid,varid,"_FillValue",FVAL)) - call fill_var_dm(data_3d, FVAL) - - open(iunit, file='DON.data', form="UNFORMATTED", status='UNKNOWN', & - access='DIRECT', recl=recl3d, convert='BIG_ENDIAN') - write(iunit,rec=1)data_3d - close(iunit) - - call check( NF90_INQ_VARID(ncid,'FET',varid) ) - call check( NF90_GET_VAR(ncid,varid,data_3d)) - call check( nf90_get_att(ncid,varid,"_FillValue",FVAL)) - call fill_var_dm(data_3d, FVAL) - - open(iunit, file='FET.data', form="UNFORMATTED", status='UNKNOWN', & - access='DIRECT', recl=recl3d, convert='BIG_ENDIAN') - write(iunit,rec=1)data_3d - close(iunit) +call from_netcdf_to_mit_3d_pickup(ncid, 'UVEL', 1, MITgcm_3D_FIELD_U) +call from_netcdf_to_mit_3d_pickup(ncid, 'VVEL', 2, MITgcm_3D_FIELD_V) +call from_netcdf_to_mit_3d_pickup(ncid, 'PTMP', 3, MITgcm_3D_FIELD) +call from_netcdf_to_mit_3d_pickup(ncid, 'PSAL', 4, MITgcm_3D_FIELD) +call from_netcdf_to_mit_2d_pickup(ncid, 'ETA') + +print *, 'Done writing physical variables into model binary files' + +if (do_bgc) then + call from_netcdf_to_mit_tracer_pickup(ncid, 'DIC', 1) + call from_netcdf_to_mit_tracer_pickup(ncid, 'ALK', 2) + call from_netcdf_to_mit_tracer_pickup(ncid, 'O2' , 3) + call from_netcdf_to_mit_tracer_pickup(ncid, 'NO3', 4) + call from_netcdf_to_mit_tracer_pickup(ncid, 'PO4', 5) + call from_netcdf_to_mit_tracer_pickup(ncid, 'FET', 6) + call from_netcdf_to_mit_tracer_pickup(ncid, 'DON', 7) + call from_netcdf_to_mit_tracer_pickup(ncid, 'DOP', 8) + call from_netcdf_to_mit_tracer_pickup(ncid, 'PHY', 9) + print *, 'Done writing biogeochemical variables into model binary files' endif call check( NF90_CLOSE(ncid) ) -deallocate(data_3d) -deallocate(data_2d) +if (compress) deallocate(Xcomp_ind, Ycomp_ind, Zcomp_ind) end subroutine dart2mit @@ -742,61 +606,820 @@ subroutine check(status) end subroutine check - !=============================================================================== -!> Check the tracer variables after reading from the binaries -!> Make sure they are non-negative -!> Do the transform if requested -!> md: mit2dart; dm: dart2mit +! 3D variable +function define_variable(ncid, VARname, nc_type, all_dimids, field) result(varid) + +integer, intent(in) :: ncid +character(len=*), intent(in) :: VARname ! variable name +integer, intent(in) :: nc_type +integer, intent(in) :: all_dimids(9) ! possible dimension ids +integer, intent(in) :: field +integer :: varid ! netcdf variable id + +integer :: dimids(3) + +if (compress) then + if (field == MITgcm_3D_FIELD) then + call check(nf90_def_var(ncid=ncid, name=VARname, xtype=nc_type, & + dimids=all_dimids(7),varid=varid)) + elseif (field == MITgcm_3D_FIELD_U) then + call check(nf90_def_var(ncid=ncid, name=VARname, xtype=nc_type, & + dimids=all_dimids(8),varid=varid)) + elseif (field == MITgcm_3D_FIELD_V) then + call check(nf90_def_var(ncid=ncid, name=VARname, xtype=nc_type, & + dimids=all_dimids(9),varid=varid)) + endif +else + dimids = which_dims(VARname, all_dimids) + call check(nf90_def_var(ncid=ncid, name=VARname, xtype=nc_type, & + dimids=dimids, varid=varid)) +endif + +end function define_variable + +!------------------------------------------------------------------ +! For the non-compressed variables, X,Y,Z dimesnions vary +! depending on the variable +function which_dims(VARname, all_dimids) result(dimids) + +character(len=*), intent(in) :: VARname ! variable name +integer, intent(in) :: all_dimids(9) +integer :: dimids(3) +! 3D variables, 3 grids: +! XC, YC, ZC 1 PSAL, PTMP, NO3, PO4, O2, PHY, ALK, DIC, DOP, DON, FET +! XG, YC, ZC 2 UVEL +! XC, YG, ZC 3 VVEL + +if (VARname == 'UVEL') then + dimids = (/all_dimids(4),all_dimids(2),all_dimids(3)/) + return +endif +if (VARname == 'VVEL') then + dimids = (/all_dimids(1),all_dimids(5),all_dimids(3)/) + return +endif + +dimids = (/all_dimids(1),all_dimids(2),all_dimids(3)/) + +end function + +!------------------------------------------------------------------ +! 2D variable +function define_variable_2d(ncid, name, nc_type, all_dimids) result(varid) + +integer, intent(in) :: ncid +character(len=*), intent(in) :: name ! variable name +integer, intent(in) :: nc_type +integer, intent(in) :: all_dimids(9) +integer :: varid ! netcdf variable id -subroutine fill_var_md(var, fillval) +! 2D variables, 1 grid: +! YC, XC 1 ETA, CHL + +if (compress) then + call check(nf90_def_var(ncid=ncid, name=name, xtype=nc_type, & + dimids = (/all_dimids(6)/),varid=varid)) +else + call check(nf90_def_var(ncid=ncid, name=name, xtype=nc_type, & + dimids = (/all_dimids(1),all_dimids(2)/),varid=varid)) +endif + +end function define_variable_2d + +!------------------------------------------------------------------ +subroutine add_attributes_to_variable(ncid, varid, long_name, units, units_long_name) + +integer, intent(in) :: ncid, varid ! which file, which variable +character(len=*), intent(in) :: long_name, units, units_long_name + +call check(nf90_put_att(ncid, varid, "long_name" , long_name)) +call check(nf90_put_att(ncid, varid, "missing_value" , FVAL)) +call check(nf90_put_att(ncid, varid, "_FillValue" , FVAL)) +call check(nf90_put_att(ncid, varid, "units" , units)) +call check(nf90_put_att(ncid, varid, "units_long_name", units_long_name)) + +end subroutine + +!------------------------------------------------------------------ +subroutine from_mit_to_netcdf_3d(mitfile, ncid, varid, field) + +character(len=*), intent(in) :: mitfile +integer, intent(in) :: ncid, varid, field ! which file, which variable, grid type + +integer :: iunit +real(r4) :: var_data(Nx,Ny,Nz) + +iunit = get_unit() +! HK are the mit files big endian by default? +open(iunit, file=mitfile, form='UNFORMATTED', status='OLD', & + access='DIRECT', recl=recl3d, convert='BIG_ENDIAN') +read(iunit,rec=1) var_data +close(iunit) -real(r4), intent(inout) :: var(:, :, :) -real(r4), intent(in) :: fillval +where (var_data == binary_fill) var_data = FVAL !HK do we also need a check for nans here? +if (compress) then + call write_compressed(ncid, varid, var_data, field) +else + call check(nf90_put_var(ncid,varid,var_data)) +endif + +end subroutine from_mit_to_netcdf_3d + +!------------------------------------------------------------------ +subroutine from_mit_to_netcdf_2d(mitfile, ncid, varid) + +character(len=*), intent(in) :: mitfile +integer, intent(in) :: ncid, varid ! which file, which variable + +integer :: iunit +real(r4) :: var_data(Nx,Ny), var_T_data(Nx,Ny,Nz) + +iunit = get_unit() +! HK are the mit files big endian by default? +open(iunit, file=mitfile, form='UNFORMATTED', status='OLD', & + access='DIRECT', recl=recl2d, convert='BIG_ENDIAN') +read(iunit,rec=1) var_data +close(iunit) + +! Manually get PTMP surface layer +open(iunit, file='PTMP.data', form='UNFORMATTED', status='OLD', & + access='DIRECT', recl=recl3d, convert='BIG_ENDIAN') +read(iunit,rec=1) var_T_data +close(iunit) + +where (var_T_data(:,:,1) == binary_fill) var_data = FVAL !HK do we also need a check for nans here? + +if (compress) then + call write_compressed(ncid, varid, var_data) +else + call check(nf90_put_var(ncid,varid,var_data)) +endif + +end subroutine from_mit_to_netcdf_2d + + +!------------------------------------------------------------------ +subroutine from_mit_to_netcdf_tracer_3d(mitfile, ncid, varid) + +character(len=*), intent(in) :: mitfile +integer, intent(in) :: ncid, varid ! which file, which variable + +integer :: iunit +real(r4) :: var_data(Nx,Ny,Nz) real(r4) :: low_conc -if (.not. module_initialized) call static_init_trans +low_conc = 1.0e-12 + +iunit = get_unit() +! HK are the mit files big endian by default? +open(iunit, file=mitfile, form='UNFORMATTED', status='OLD', & + access='DIRECT', recl=recl3d, convert='BIG_ENDIAN') +read(iunit,rec=1) var_data +close(iunit) + +! CHL is treated differently - HK CHL is 2d so you will not enter this +if (mitfile=='CHL.data') then + where (var_data == binary_fill) + var_data = FVAL + elsewhere + var_data = log10(var_data) + endwhere +else + ! Make sure the tracer concentration is positive + where(var_data < binary_fill) var_data = low_conc + + if (log_transform) then + where (var_data == binary_fill) + var_data = FVAL + elsewhere + var_data = log(var_data) + endwhere + else + where (var_data == binary_fill) var_data = FVAL + endif +endif + +if (compress) then + call write_compressed(ncid, varid, var_data, MITgcm_3D_FIELD) +else + call check(nf90_put_var(ncid,varid,var_data)) +endif + +end subroutine from_mit_to_netcdf_tracer_3d + +!------------------------------------------------------------------ +subroutine from_mit_to_netcdf_tracer_2d(mitfile, ncid, varid) + +character(len=*), intent(in) :: mitfile +integer, intent(in) :: ncid, varid ! which file, which variable + +integer :: iunit +real(r4) :: var_data(Nx,Ny) +real(r4) :: low_conc low_conc = 1.0e-12 -! Make sure the tracer concentration is positive -where(var < 0.0_r4) var = low_conc +iunit = get_unit() +! HK are the mit files big endian by default? +open(iunit, file=mitfile, form='UNFORMATTED', status='OLD', & + access='DIRECT', recl=recl3d, convert='BIG_ENDIAN') +read(iunit,rec=1) var_data +close(iunit) + +! CHL is treated differently +if (mitfile=='CHL.data') then + where (var_data == binary_fill) + var_data = FVAL + elsewhere + var_data = log10(var_data) + endwhere +else + ! Make sure the tracer concentration is positive + where(var_data < binary_fill) var_data = low_conc + + if (log_transform) then + where (var_data == binary_fill) + var_data = FVAL + elsewhere + var_data = log(var_data) + endwhere + else + where (var_data == binary_fill) var_data = FVAL + endif +endif + +if (compress) then + call write_compressed(ncid, varid, var_data) +else + call check(nf90_put_var(ncid,varid,var_data)) +endif + +end subroutine from_mit_to_netcdf_tracer_2d + +!------------------------------------------------------------------ +subroutine from_netcdf_to_mit_2d(ncid, name) + +integer, intent(in) :: ncid ! which file, +character(len=*), intent(in) :: name ! which variable + +integer :: iunit +real(r4) :: var(Nx,Ny) +integer :: varid +real(r4) :: local_fval + +call check( NF90_INQ_VARID(ncid,name,varid) ) +call check( nf90_get_att(ncid,varid,"_FillValue",local_fval)) +! initialize var to netcdf fill value +var(:,:) = local_fval + +if (compress) then + call read_compressed(ncid, varid, var) +else + call check(nf90_get_var(ncid,varid,var)) +endif + +where (var == local_fval) var = binary_fill + +iunit = get_unit() +open(iunit, file=trim(name)//'.data', form="UNFORMATTED", status='UNKNOWN', & + access='DIRECT', recl=recl2d, convert='BIG_ENDIAN') +write(iunit,rec=1)var +close(iunit) + +end subroutine from_netcdf_to_mit_2d + +!------------------------------------------------------------------ +subroutine from_netcdf_to_mit_3d(ncid, name, field) + +integer, intent(in) :: ncid ! which file, +character(len=*), intent(in) :: name ! which variable + +integer :: iunit, field +real(r4) :: var(Nx,Ny,Nz) +integer :: varid +real(r4) :: local_fval + +call check( NF90_INQ_VARID(ncid,name,varid) ) +call check( nf90_get_att(ncid,varid,"_FillValue",local_fval)) +! initialize var to netcdf fill value +var(:,:,:) = local_fval + +if (compress) then + call read_compressed(ncid, varid, var, field) +else + call check(nf90_get_var(ncid,varid,var)) +endif + +where (var == local_fval) var = binary_fill + +iunit = get_unit() +open(iunit, file=trim(name)//'.data', form="UNFORMATTED", status='UNKNOWN', & + access='DIRECT', recl=recl3d, convert='BIG_ENDIAN') +write(iunit,rec=1)var +close(iunit) + +end subroutine from_netcdf_to_mit_3d + +!------------------------------------------------------------------ +subroutine from_netcdf_to_mit_2d_pickup(ncid, name) + +integer, intent(in) :: ncid ! which file, +character(len=*), intent(in) :: name ! which variable + +integer :: iunit +real(r4) :: var(Nx,Ny) +real(r8) :: var8(Nx,Ny) +integer :: varid +real(r4) :: local_fval + + + +call check( NF90_INQ_VARID(ncid,name,varid) ) +call check( nf90_get_att(ncid,varid,"_FillValue",local_fval)) + +! initialize var to netcdf fill value +var(:,:) = local_fval + +if (compress) then + call read_compressed(ncid, varid, var) +else + call check(nf90_get_var(ncid,varid,var)) +endif + +where (var == local_fval) var = binary_fill +var8 = var + + +if (do_bgc) then + write(iunit,rec=401) var8 +else + write(iunit,rec=481) var8 +endif +close(iunit) + +end subroutine from_netcdf_to_mit_2d_pickup + +!------------------------------------------------------------------ +subroutine from_netcdf_to_mit_3d_pickup(ncid, name, lev, field) + +integer, intent(in) :: ncid ! which file, +character(len=*), intent(in) :: name ! which variable + +integer :: iunit, lev, field +real(r4) :: var(Nx,Ny,Nz) +real(r8) :: var8(Nx,Ny,Nz) +integer :: varid, i +real(r4) :: local_fval +integer :: LB, RB, RF + + +call check( NF90_INQ_VARID(ncid,name,varid) ) +call check( nf90_get_att(ncid,varid,"_FillValue",local_fval)) + +! initialize var to netcdf fill value +var(:,:,:) = local_fval + +if (compress) then + call read_compressed(ncid, varid, var, field) +else + call check(nf90_get_var(ncid,varid,var)) +endif + +where (var == local_fval) var = binary_fill +var8 = var + + + +LB = Nz * (lev-1) + 1 +RB = Nz * lev +RF = Nz * (lev-1) +do i = LB, RB + write(iunit,rec=i) var8(:, :, i - RF) +enddo +close(iunit) + +end subroutine from_netcdf_to_mit_3d_pickup + +!------------------------------------------------------------------ +subroutine from_netcdf_to_mit_tracer(ncid, name) + +integer, intent(in) :: ncid ! which file +character(len=*), intent(in) :: name ! which variable + +integer :: iunit +real(r4) :: var(Nx,Ny,Nz) +integer :: varid +real(r4) :: local_fval + +call check( NF90_INQ_VARID(ncid,name,varid) ) +call check( nf90_get_att(ncid,varid,"_FillValue",local_fval)) +! initialize var to netcdf fill value +var(:,:,:) = local_fval + +if (compress) then + call read_compressed(ncid, varid, var, MITgcm_3D_FIELD) +else + call check(nf90_get_var(ncid,varid,var)) +endif if (log_transform) then - where (var == 0.0_r4) - var = fillval + where (var == local_fval) + var = binary_fill elsewhere - var = log(var) + var = exp(var) endwhere else - where (var == 0.0_r4) var = fillval + where (var == local_fval) var = binary_fill endif -end subroutine +iunit = get_unit() +open(iunit, file=trim(name)//'.data', form="UNFORMATTED", status='UNKNOWN', & + access='DIRECT', recl=recl3d, convert='BIG_ENDIAN') +write(iunit,rec=1)var +close(iunit) + +end subroutine from_netcdf_to_mit_tracer !------------------------------------------------------------------ +subroutine from_netcdf_to_mit_tracer_pickup(ncid, name, lev) -subroutine fill_var_dm(var, fillval) +integer, intent(in) :: ncid ! which file +character(len=*), intent(in) :: name ! which variable -real(r4), intent(inout) :: var(:, :, :) -real(r4), intent(in) :: fillval +integer :: iunit, lev +real(r4) :: var(Nx,Ny,Nz) +real(r8) :: var8(Nx,Ny,Nz) +integer :: varid +real(r4) :: local_fval +real(r4) :: low_conc, large_conc = 5.0 ! From Siva's old code -if (.not. module_initialized) call static_init_trans +low_conc = 1.0e-12 -if (log_transform) then - where (var == fillval) - var = 0.0_r4 +call check( NF90_INQ_VARID(ncid,name,varid) ) +call check( nf90_get_att(ncid,varid,"_FillValue",local_fval)) + +! initialize var to netcdf fill value +var(:,:,:) = local_fval + +if (compress) then + call read_compressed(ncid, varid, var, MITgcm_3D_FIELD) +else + call check(nf90_get_var(ncid,varid,var)) +endif + +if (log_transform) then + where (var == local_fval) + var = binary_fill elsewhere var = exp(var) endwhere else - where (var == fillval) var = 0.0_r4 + where (var == local_fval) var = binary_fill + where (var > large_conc) var = low_conc endif -end subroutine +var8 = var + +iunit = get_unit() +open(iunit, file='PICKUP_PTRACERS.OUTPUT', form="UNFORMATTED", status='UNKNOWN', & + access='DIRECT', recl=recl3d, convert='BIG_ENDIAN') +write(iunit,rec=lev) var8 +close(iunit) + +end subroutine from_netcdf_to_mit_tracer_pickup + +!------------------------------------------------------------------ +subroutine from_netcdf_to_mit_tracer_chl(ncid, name) + +integer, intent(in) :: ncid ! which file +character(len=*), intent(in) :: name ! which variable + +integer :: iunit +real(r4) :: var(Nx,Ny) +integer :: varid +real(r4) :: local_fval + +call check( NF90_INQ_VARID(ncid,name,varid) ) +call check( nf90_get_att(ncid,varid,"_FillValue",local_fval)) +! initialize var to netcdf fill value +var(:,:) = local_fval + +if (compress) then + call read_compressed(ncid, varid, var) +else + call check(nf90_get_var(ncid,varid,var)) +endif + +where (var == local_fval) + var = binary_fill +elsewhere + var = 10**(var) +endwhere + + +iunit = get_unit() +open(iunit, file=trim(name)//'.data', form="UNFORMATTED", status='UNKNOWN', & + access='DIRECT', recl=recl2d, convert='BIG_ENDIAN') +write(iunit,rec=1)var +close(iunit) + +end subroutine from_netcdf_to_mit_tracer_chl + + +!------------------------------------------------------------------ +! Assumes all 3D variables are masked in the +! same location +function get_compressed_size_3d(field) result(n3) + +integer :: n3, field +integer :: iunit +real(r4) :: var3d(NX,NY,NZ) +integer :: i, j, k +character(len=MAX_LEN_FNAM) :: source + +if (field == MITgcm_3D_FIELD) source = 'PSAL.data' +if (field == MITgcm_3D_FIELD_U) source = 'UVEL.data' +if (field == MITgcm_3D_FIELD_V) source = 'VVEL.data' + +iunit = get_unit() +open(iunit, file=trim(source), form='UNFORMATTED', status='OLD', & + access='DIRECT', recl=recl3d, convert='BIG_ENDIAN') +read(iunit,rec=1) var3d +close(iunit) + +n3 = 0 + +! Get compressed size +do i=1,NX + do j=1,NY + do k=1,NZ + if (var3d(i,j,k) /= binary_fill) then !HK also NaN? + n3 = n3 + 1 + endif + enddo + enddo +enddo + +end function get_compressed_size_3d !------------------------------------------------------------------ +! Assumes all 2D variables are masked in the +! same location +function get_compressed_size_2d() result(n2) + +integer :: n2 +integer :: iunit +real(r4) :: var3d(NX,NY,NZ) +integer :: i,j + +iunit = get_unit() +open(iunit, file='PTMP.data', form='UNFORMATTED', status='OLD', & + access='DIRECT', recl=recl3d, convert='BIG_ENDIAN') +read(iunit,rec=1) var3d +close(iunit) + +n2 = 0 + +! Get compressed size +do i=1,NX + do j=1,NY + if (var3d(i,j,1) /= binary_fill) then !HK also NaN? + n2 = n2 + 1 + endif + enddo +enddo + +end function get_compressed_size_2d + +!------------------------------------------------------------------ +subroutine fill_compressed_coords() + +!XG,etc read from PARAM04 in static_init_trans +real(r4) :: var3d(NX,NY,NZ) +integer :: n, i, j, k + +iunit = get_unit() +open(iunit, file='PSAL.data', form='UNFORMATTED', status='OLD', & + access='DIRECT', recl=recl3d, convert='BIG_ENDIAN') +read(iunit,rec=1) var3d +close(iunit) + +n = 1 + +do k=1,NZ ! k first so 2d is first + do i=1,NX + do j=1,NY + if (var3d(i,j,k) /= binary_fill) then !HK also NaN? + XCcomp(n) = XC(i) + YCcomp(n) = YC(j) + ZCcomp(n) = ZC(k) + XGcomp(n) = XG(i) + YGcomp(n) = YG(j) + ZGcomp(n) = ZG(k) + + Xcomp_ind(n) = i ! Assuming grids are compressed the same + Ycomp_ind(n) = j + Zcomp_ind(n) = k + + n = n + 1 + endif + enddo + enddo +enddo + +! UVEL: +iunit = get_unit() +open(iunit, file='UVEL.data', form='UNFORMATTED', status='OLD', & + access='DIRECT', recl=recl3d, convert='BIG_ENDIAN') +read(iunit,rec=1) var3d +close(iunit) + +n = 1 + +do k=1,NZ ! k first so 2d is first + do i=1,NX + do j=1,NY + if (var3d(i,j,k) /= binary_fill) then !HK also NaN? + Xcomp_indU(n) = i + Ycomp_indU(n) = j + Zcomp_indU(n) = k + + n = n + 1 + endif + enddo + enddo +enddo + +! VVEL: +iunit = get_unit() +open(iunit, file='VVEL.data', form='UNFORMATTED', status='OLD', & + access='DIRECT', recl=recl3d, convert='BIG_ENDIAN') +read(iunit,rec=1) var3d +close(iunit) + +n = 1 + +do k=1,NZ ! k first so 2d is first + do i=1,NX + do j=1,NY + if (var3d(i,j,k) /= binary_fill) then !HK also NaN? + Xcomp_indV(n) = i + Ycomp_indV(n) = j + Zcomp_indV(n) = k + + n = n + 1 + endif + enddo + enddo +enddo + +end subroutine fill_compressed_coords + +!------------------------------------------------------------------ +subroutine write_compressed_2d(ncid, varid, var_data) + +integer, intent(in) :: ncid, varid +real(r4), intent(in) :: var_data(Nx,Ny) + +real(r4) :: comp_var(ncomp2) +integer :: n +integer :: i,j ! loop variables + +n = 1 +do i = 1, NX + do j = 1, NY + if (var_data(i,j) /= FVAL) then + comp_var(n) = var_data(i,j) + n = n + 1 + endif + enddo +enddo + +call check(nf90_put_var(ncid,varid,comp_var)) + +end subroutine write_compressed_2d + +!------------------------------------------------------------------ +subroutine write_compressed_3d(ncid, varid, var_data, field) + +integer, intent(in) :: ncid, varid, field +real(r4), intent(in) :: var_data(Nx,Ny,Nz) + +real(r4), allocatable :: comp_var(:) +integer :: n +integer :: i,j,k ! loop variables + +if (field == MITgcm_3D_FIELD_U) then + allocate(comp_var(ncomp3U)) + do i = 1,ncomp3U + comp_var(i) = var_data(Xcomp_indU(i), Ycomp_indU(i), Zcomp_indU(i)) + enddo + +elseif (field == MITgcm_3D_FIELD_V) then + allocate(comp_var(ncomp3V)) + do i = 1,ncomp3V + comp_var(i) = var_data(Xcomp_indV(i), Ycomp_indV(i), Zcomp_indV(i)) + enddo + +else + allocate(comp_var(ncomp3)) + do i = 1,ncomp3 + comp_var(i) = var_data(Xcomp_ind(i), Ycomp_ind(i), Zcomp_ind(i)) + enddo +endif + +!n = 1 +!do k = 1, NZ !k first so 2d is first +! do i = 1, NX +! do j = 1, NY +! if (var_data(i,j,k) /= FVAL) then +! print *, 'n: ', n, ', var_data(i,j,k): ', var_data(i,j,k) +! comp_var(n) = var_data(i,j,k) +! n = n + 1 +! endif +! enddo +! enddo +!enddo + +call check(nf90_put_var(ncid,varid,comp_var)) + +deallocate(comp_var) + +end subroutine write_compressed_3d + +!------------------------------------------------------------------ +subroutine read_compressed_2d(ncid, varid, var) + +integer, intent(in) :: ncid, varid +real(r4), intent(inout) :: var(NX,NY) + +real(r4) :: comp_var(ncomp2) +integer :: n ! loop variable +integer :: i,j,k ! x,y,z +integer :: c + +c = 1 + +call check(nf90_get_var(ncid,varid,comp_var)) + +do n = 1, ncomp3 + i = Xcomp_ind(n) + j = Ycomp_ind(n) + k = Zcomp_ind(n) + if (k == 1) then + var(i,j) = comp_var(c) + c = c + 1 + endif +enddo + +end subroutine read_compressed_2d + +!------------------------------------------------------------------ +subroutine read_compressed_3d(ncid, varid, var, field) + +integer, intent(in) :: ncid, varid, field +real(r4), intent(inout) :: var(NX,NY,NZ) + +real(r4), allocatable :: comp_var(:) +integer :: n ! loop variable +integer :: i,j,k ! x,y,k + +if (field == MITgcm_3D_FIELD_U) then + allocate(comp_var(ncomp3U)) + call check(nf90_get_var(ncid,varid,comp_var)) + do n = 1, ncomp3U + i = Xcomp_indU(n) + j = Ycomp_indU(n) + k = Zcomp_indU(n) + var(i,j,k) = comp_var(n) + enddo + +elseif (field == MITgcm_3D_FIELD_V) then + allocate(comp_var(ncomp3V)) + call check(nf90_get_var(ncid,varid,comp_var)) + do n = 1, ncomp3V + i = Xcomp_indV(n) + j = Ycomp_indV(n) + k = Zcomp_indV(n) + var(i,j,k) = comp_var(n) + enddo + +else + allocate(comp_var(ncomp3)) + call check(nf90_get_var(ncid,varid,comp_var)) + do n = 1, ncomp3 + i = Xcomp_ind(n) + j = Ycomp_ind(n) + k = Zcomp_ind(n) + var(i,j,k) = comp_var(n) + enddo +endif + +deallocate(comp_var) + +end subroutine read_compressed_3d end module trans_mitdart_mod diff --git a/models/MITgcm_ocean/work/input.nml b/models/MITgcm_ocean/work/input.nml index f92ab19c39..0c462c2266 100644 --- a/models/MITgcm_ocean/work/input.nml +++ b/models/MITgcm_ocean/work/input.nml @@ -457,7 +457,7 @@ # quantity_of_interest = 'QTY_DENSITY' &model_mod_check_nml - input_state_files = 'OUTPUT.nc' + input_state_files = 'mem01_reduced.nc' output_state_files = 'check_me' verbose = .TRUE. test1thru = 0 diff --git a/models/aether_lat-lon/aether_to_dart.f90 b/models/aether_lat-lon/aether_to_dart.f90 new file mode 100644 index 0000000000..dca1ef1c6d --- /dev/null +++ b/models/aether_lat-lon/aether_to_dart.f90 @@ -0,0 +1,472 @@ +! DART software - Copyright 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 +! + +program aether_to_dart + +!---------------------------------------------------------------------- +! purpose: Transform the Aether model restarts into a DART filter_input.nc. +! +! method: Read aether "restart" files of model state (multiple files, +! one block per aether mpi task) +! Reform fields into a DART netcdf file +! +! USAGE: The aether restart dirname and output filename are read from +! the aether_to_dart_nml namelist. +! +!---------------------------------------------------------------------- +! Converts Aether restart files to a netCDF file + +use types_mod, only : r4, MISSING_I, vtablenamelength + +use time_manager_mod, only: time_type + +use utilities_mod, only : & + finalize_utilities, error_handler, E_ERR, E_MSG, E_WARN, & + initialize_utilities, do_output + +use default_model_mod, only : write_model_time + +use transform_state_mod, only : & + static_init_blocks, aether_name_to_dart, & + nghost, open_block_file, aether_restart_dirname, & + VT_ORIGININDX, VT_VARNAMEINDX, nvar_neutral, nvar_ion, & + nx_per_block, ny_per_block, nz_per_block, & + nblocks_lon, nblocks_lat, variables, & + lats, levs, lons, debug, state_time, & + block_file_name, nlat, nlon, nlev, purge_chars + +use netcdf_utilities_mod, only : & + nc_create_file, nc_close_file, & + nc_begin_define_mode, nc_end_define_mode, & + nc_define_dimension, & + nc_add_global_attribute, nc_add_global_creation_time, & + nc_get_attribute_from_variable, nc_add_attribute_to_variable, & + nc_define_real_variable, nc_define_real_scalar, & + nc_get_variable, nc_put_variable, & + nc_synchronize_file + +implicit none + +!---------------------------------------------------------------------- +! global storage +!---------------------------------------------------------------------- + +integer :: member = MISSING_I, & + num_args, ncid +character(len=3) :: char_mem +character(len=31) :: filter_io_root = 'filter_input' +character(len=64) :: filter_io_file = '' +character(len=512) :: error_string_1, error_string_2 +character(len=31), parameter :: progname = 'aether_to_dart' +character(len=256), parameter :: source = 'aether_lat-lon/aether_to_dart.f90' + +character(len=4), parameter :: LEV_DIM_NAME = 'alt' +character(len=4), parameter :: LAT_DIM_NAME = 'lat' +character(len=4), parameter :: LON_DIM_NAME = 'lon' +character(len=4), parameter :: TIME_DIM_NAME = 'time' + +character(len=4), parameter :: LEV_VAR_NAME = 'alt' +character(len=4), parameter :: LAT_VAR_NAME = 'lat' +character(len=4), parameter :: LON_VAR_NAME = 'lon' +character(len=4), parameter :: TIME_VAR_NAME = 'time' + +!====================================================================== + +call initialize_utilities(progname) + +!---------------------------------------------------------------------- +! Get the ensemble member +!---------------------------------------------------------------------- +num_args = command_argument_count() +if (num_args == 0) then + write(error_string_1,*) 'Usage: ./aether_to_dart member_number (0-based)' + call error_handler(E_ERR, progname, error_string_1) +endif + +call get_command_argument(1,char_mem) +read(char_mem,'(I3)') member + +!---------------------------------------------------------------------- +! Convert the files +!---------------------------------------------------------------------- + +call static_init_blocks(member) + +! Must be after static_init_blocks, which provides filter_io_root from the namelist. +write(filter_io_file,'(2A, I0.4, A3)') trim(filter_io_root),'_', member + 1,'.nc' +call error_handler(E_MSG, '', '') +write(error_string_1,'(A,I3,2A)') 'Converting Aether member ',member, & + ' restart files to the NetCDF file ', trim(filter_io_file) +write(error_string_2,'(3A)') ' in directory ', trim(aether_restart_dirname) +call error_handler(E_MSG, progname, error_string_1, text2=error_string_2) +call error_handler(E_MSG, '', '') + +! nc_create_file does not leave define mode. +ncid = nc_create_file(trim(aether_restart_dirname)//'/'//trim(filter_io_file)) +! def_fill_dimvars does leave define mode. +call def_fill_dimvars(ncid) + +! Write_model_time will make a time variable, if needed, which it is not. +! state_time is read in transform_state_mod and is available by the use statement. +call write_model_time(ncid, state_time) + +! Define (non-time) variables +call restarts_to_filter(ncid, member, define=.true.) + +! Read and convert (non-time) variables +call restarts_to_filter(ncid, member, define=.false.) +! subr. called by this routine closes the file only if define = .true. +call nc_close_file(ncid) + +call error_handler(E_MSG, '', '') +write(error_string_1,'(3A)') 'Successfully converted the Aether restart files to ', & + "'"//trim(filter_io_file)//"'" +call error_handler(E_MSG, progname, error_string_1) +call error_handler(E_MSG, '', '') + +! end - close the log, etc +call finalize_utilities() + +!----------------------------------------------------------------------- +contains + +!----------------------------------------------------------------------- +! Open all restart files (blocks x {neutrals,ions}) for 1 member +! and transfer the requested variable contents to the filter input file. +! This is called with 'define' = +! .true. define variables in the file or +! .false. transfer the data from restart files to a filter_inpu.nc file. + +subroutine restarts_to_filter(ncid_output, member, define) + +integer, intent(in) :: ncid_output, member +logical, intent(in) :: define + +integer :: ib, jb, ib_loop, jb_loop + +if (define) then + ! if define, run one block. + ! the block_to_filter_io call defines the variables in the whole domain netCDF file. + ib_loop = 1 + jb_loop = 1 + call nc_begin_define_mode(ncid_output) +else + ! if not define, run all blocks. + ! the block_to_filter_io call adds the (ib,jb) block to the netCDF variables + ! in order to make a file containing the data for all the blocks. + ib_loop = nblocks_lon + jb_loop = nblocks_lat +end if + +do jb = 1, jb_loop + do ib = 1, ib_loop + call block_to_filter_io(ncid_output, ib, jb, member, define) + enddo +enddo + +if (define) then + call nc_end_define_mode(ncid_output) +endif + +end subroutine restarts_to_filter + +!----------------------------------------------------------------------- +! Transfer variable data from a block restart file to the filter_input.nc file. +! It's called with 2 modes: +! define = .true. define the NC variables in the filter_input.nc +! define = .false. write the data from a block to the NC file using write_filter_io. + +subroutine block_to_filter_io(ncid_output, ib, jb, member, define) + +integer, intent(in) :: ncid_output +integer, intent(in) :: ib, jb +integer, intent(in) :: member +logical, intent(in) :: define + +real(r4), allocatable :: temp1d(:), temp2d(:,:), temp3d(:,:,:) +! real(r4), allocatable :: alt1d(:), density_ion_e(:,:,:) +integer :: ivar, nb, ncid_input +! TEC? integer :: maxsize +! logical :: no_idensity +! real(r4) :: temp0d +character(len=32) :: att_val +character(len=128) :: file_root +character(len=256) :: filename +character(len=vtablenamelength) :: varname, dart_varname + +character(len=*), parameter :: routine = 'block_to_filter_io' + +! The block number, as counted in Aether. +! Lower left is 0, increase to the East, then 1 row farther north, West to East. +nb = (jb - 1) * nblocks_lon + ib - 1 + +! a temp array large enough to hold any of the +! Lon,Lat or Alt array from a block plus ghost cells +allocate(temp1d(1-nghost:max(nx_per_block, ny_per_block, nz_per_block) + nghost)) + +! treat alt specially since we want to derive TEC here +! TODO: See density_ion_e too. +! allocate( alt1d(1-nghost:max(nx_per_block, ny_per_block, nz_per_block) + nghost)) + +! temp array large enough to hold any 2D field +allocate(temp2d(1-nghost:ny_per_block+nghost, & + 1-nghost:nx_per_block+nghost)) + +! TODO: We need all altitudes, but there might be vertical blocks in the future. +! But there would be no vertical halos. +! Make transform_state_mod: zcount adapt to whether there are blocks. +! Temp needs to have C-ordering, which is what the restart files have. +! temp array large enough to hold 1 species, temperature, etc +allocate(temp3d(1:nz_per_block, & + 1-nghost:ny_per_block+nghost, & + 1-nghost:nx_per_block+nghost)) + +! TODO: Waiting for e- guidance from Aaron. +! save density_ion_e to compute TEC +! allocate(density_ion_e(1:nz_per_block, & +! 1-nghost:ny_per_block+nghost, & +! 1-nghost:nx_per_block+nghost)) + +! TODO: Aether gives a unique name to each (of 6) velocity components. +! Do we want to use a temp4d array to handle them? +! They are independent variables in the block files (and state). +! ! temp array large enough to hold velocity vect, etc +! maxsize = max(3, nvar_ion) +! allocate(temp4d(1-nghost:nx_per_block+nghost, & +! 1-nghost:ny_per_block+nghost, & +! 1-nghost:nz_per_block+nghost, maxsize)) + + +! TODO; Does Aether need a replacement for these Density fields? Yes. +! But they are probably read by the loops below. +! Don't need to fetch index because Aether has NetCDF restarts, +! so just loop over the field names to read. +! +! ! assume we could not find the electron density for VTEC calculations +! no_idensity = .true. +! +! if (inum > 0) then +! ! one or more items in the state vector need to replace the +! ! data in the output file. loop over the index list in order. +! j = 1 +! ! TODO: electron density is not in the restart files, but it's needed for TEC +! In Aether they will be from an ions file, but now only from an output file (2023-10-30). +! Can that be handled like the neutrals and ions files, using variables(VT_ORIGININDX,:) +! to build an output file name? Are outputs in block form? +! ! save the electron density for TEC computation +! density_ion_e(:,:,:) = temp3d(:,:,:) + +! Handle the 2 restart file types (ions and neutrals). +! Each field has a file type associated with it: variables(VT_ORIGININDX,f_index) + +file_root = variables(VT_ORIGININDX,1) +write(filename,'(A,"/",A)') trim(aether_restart_dirname), & + trim(block_file_name(trim(file_root), member, nb)) +ncid_input = open_block_file(filename, 'read') + +do ivar = 1, nvar_neutral + ! The nf90 functions cannot read the variable names with the '\'s in them. + varname = purge_chars(trim(variables(VT_VARNAMEINDX,ivar)), '\', plus_minus=.false.) + if (debug >= 100 .and. do_output()) print*, routine,'varname = ', varname + ! Translate the Aether field name into a CF-compliant DART field name. + dart_varname = aether_name_to_dart(varname) + + ! TODO: Given the subroutine name, perhaps these definition sections should be + ! one call higher up, with the same loop around it. + if (define) then + ! Define the variable in the filter_input.nc file (the output from this program). + ! The calling routine entered define mode. + + if (debug > 10 .and. do_output()) then + write(error_string_1,'(A,I0,2A)') 'Defining ivar = ', ivar,':', dart_varname + call error_handler(E_MSG, routine, error_string_1, source) + end if + + call nc_define_real_variable(ncid_output, dart_varname, & + (/ LEV_DIM_NAME, LAT_DIM_NAME, LON_DIM_NAME/) ) + call nc_get_attribute_from_variable(ncid_input, varname, 'units', att_val, routine) + call nc_add_attribute_to_variable(ncid_output, dart_varname, 'units', att_val, routine) + + else if (file_root == 'neutrals') then + ! Read 3D array and extract the non-halo data of this block. + call nc_get_variable(ncid_input, varname, temp3d, context=routine) + call write_filter_io(temp3d, dart_varname, ib, jb, ncid_output) + else + write(error_string_1,'(A,I3,A)') 'Trying to read neutrals, but variables(', & + VT_ORIGININDX,ivar , ') /= "neutrals"' + call error_handler(E_ERR, routine, error_string_1, source) + endif + +enddo +call nc_close_file(ncid_input) + +file_root = variables(VT_ORIGININDX,nvar_neutral+1) +write(filename,'(A,"/",A)') trim(aether_restart_dirname), & + trim(block_file_name(trim(file_root), member, nb)) +ncid_input = open_block_file(filename, 'read') + +do ivar = nvar_neutral +1, nvar_neutral + nvar_ion + ! Purging \ from aether name. + varname = purge_chars(trim(variables(VT_VARNAMEINDX,ivar)), '\', plus_minus=.false.) + dart_varname = aether_name_to_dart(varname) + + if (define) then + + if (debug > 10 .and. do_output()) then + write(error_string_1,'(A,I0,2A)') 'Defining ivar = ', ivar,':', dart_varname + call error_handler(E_MSG, routine, error_string_1, source) + end if + + call nc_define_real_variable(ncid_output, dart_varname, & + (/ LEV_DIM_NAME, LAT_DIM_NAME, LON_DIM_NAME/) ) + call nc_get_attribute_from_variable(ncid_input, varname, 'units', att_val, routine) + call nc_add_attribute_to_variable(ncid_output, dart_varname, 'units', att_val, routine) + print*, routine,': defined ivar, dart_varname, att = ', & + ivar, trim(dart_varname), trim(att_val) + + else if (file_root == 'ions') then + call nc_get_variable(ncid_input, varname, temp3d, context=routine) + call write_filter_io(temp3d, dart_varname, ib, jb, ncid_output) + else + write(error_string_1,'(A,I3,A)') 'Trying to read ions, but variables(', & + VT_ORIGININDX,ivar , ') /= "ions"' + call error_handler(E_ERR, routine, error_string_1, source) + endif + +enddo + +! Leave file open if fields were just added (define = .false.), +! so that time can be added. +if (define) call nc_close_file(ncid_input) + +! TODO: Does Aether need TEC to be calculated? Yes +! ! add the VTEC as an extended-state variable +! ! NOTE: This variable will *not* be written out to the Aether restart files +! +! if (no_idensity) then +! write(error_string_1,*) 'Cannot compute the VTEC without the electron density' +! call error_handler(E_ERR, routine, error_string_1, source) +! end if +! +! temp2d = 0._r8 +! ! compute the TEC integral +! do i =1,nz_per_block-1 ! approximate the integral over the altitude as a sum of trapezoids +! ! area of a trapezoid: A = (h2-h1) * (f2+f1)/2 +! temp2d(:,:) = temp2d(:,:) + ( alt1d(i+1)-alt1d(i) ) * & +! ( density_ion_e(:,:,i+1)+density_ion_e(:,:,i) ) /2.0_r8 +! end do +! ! convert temp2d to TEC units +! temp2d = temp2d/1e16_r8 +! call write_block_to_filter2d(temp2d, ivals(1), block, ncid, define) + +! TODO: Does Aether need f10_7 to be calculated or processed? Yes +! !gitm_index = get_index_start(domain_id, 'VerticalVelocity') +! call get_index_from_gitm_varname('f107', inum, ivals) +! if (inum > 0) then +! call write_block_to_filter0d(temp0d, ivals(1), ncid, define) !see comments in the body of the subroutine +! endif +! + +deallocate(temp1d, temp2d, temp3d) +! deallocate(alt1d, density_ion_e) + +end subroutine block_to_filter_io + +!----------------------------------------------------------------------- +! Open all restart files (neutrals,ions) for a block and read in the requested data items. +! The write_filter_io calls will write the data to the filter_input.nc. + +subroutine write_filter_io(data3d, varname, ib, jb, ncid) + +real(r4), intent(in) :: data3d(1:nz_per_block, & + 1-nghost:ny_per_block+nghost, & + 1-nghost:nx_per_block+nghost) + +character(len=vtablenamelength), intent(in) :: varname +integer, intent(in) :: ib, jb +integer, intent(in) :: ncid + +integer :: starts(3) + +character(len=*), parameter :: routine = 'write_filter_io' + +! write(varname,'(A)') trim(variables(VT_VARNAMEINDX,ivar)) + +! to compute the start, consider (ib-1)*nx_per_block+1 +starts(1) = 1 +starts(2) = (jb-1) * ny_per_block + 1 +starts(3) = (ib-1) * nx_per_block + 1 + +call nc_put_variable(ncid, varname, & + data3d(1:nz_per_block,1:ny_per_block,1:nx_per_block), & + context=routine, nc_start=starts, & + nc_count=(/nz_per_block,ny_per_block,nx_per_block/)) + +end subroutine write_filter_io + +!----------------------------------------------------------------------- +! Add dimension variable contents to the file. + +subroutine def_fill_dimvars(ncid) + +integer, intent(in) :: ncid + +character(len=*), parameter :: routine = 'def_fill_dimvars' + +! File is still in define mode from nc_create_file +! call nc_begin_define_mode(ncid) + +! Global atts for aether_to_dart and dart_to_aether. +call nc_add_global_creation_time(ncid, routine) +call nc_add_global_attribute(ncid, "model_source", source, routine) +call nc_add_global_attribute(ncid, "model", "aether", routine) + +! define grid dimensions +call nc_define_dimension(ncid, trim(LEV_DIM_NAME), nlev, routine) +call nc_define_dimension(ncid, trim(LAT_DIM_NAME), nlat, routine) +call nc_define_dimension(ncid, trim(LON_DIM_NAME), nlon, routine) + +! define grid variables +! z +call nc_define_real_variable( ncid, trim(LEV_VAR_NAME), (/ trim(LEV_DIM_NAME) /), routine) +call nc_add_attribute_to_variable(ncid, trim(LEV_VAR_NAME), 'units', 'm', routine) +call nc_add_attribute_to_variable & + (ncid, trim(LEV_VAR_NAME), 'long_name', 'height above mean sea level', routine) + +! latitude +call nc_define_real_variable( ncid, trim(LAT_VAR_NAME), (/ trim(LAT_DIM_NAME) /), routine) +call nc_add_attribute_to_variable(ncid, trim(LAT_VAR_NAME), 'units', 'degrees_north', routine) +call nc_add_attribute_to_variable(ncid, trim(LAT_VAR_NAME), 'long_name', 'latitude', routine) + +! longitude +call nc_define_real_variable( ncid, trim(LON_VAR_NAME), (/ trim(LON_VAR_NAME) /), routine) +call nc_add_attribute_to_variable(ncid, trim(LON_VAR_NAME), 'units', 'degrees_east', routine) +call nc_add_attribute_to_variable(ncid, trim(LON_VAR_NAME), 'long_name', 'longitude', routine) + +! Dimension 'time' will no longer be created by write_model_time, +! or by nc_define_unlimited_dimension. It will be a scalar variable. +! time +call nc_define_real_scalar( ncid, trim(TIME_VAR_NAME), routine) +call nc_add_attribute_to_variable(ncid, trim(TIME_VAR_NAME), 'calendar', 'gregorian', routine) +call nc_add_attribute_to_variable & + (ncid, trim(TIME_VAR_NAME), 'units', 'days since 1601-01-01 00:00:00', routine) +call nc_add_attribute_to_variable & + (ncid, trim(TIME_VAR_NAME), 'long_name', 'gregorian_days', routine) + +call nc_end_define_mode(ncid) + +call nc_put_variable(ncid, trim(LEV_VAR_NAME), levs, routine) +call nc_put_variable(ncid, trim(LAT_VAR_NAME), lats, routine) +call nc_put_variable(ncid, trim(LON_VAR_NAME), lons, routine) +! time will be written elsewhere. + +! Flush the buffer and leave netCDF file open +call nc_synchronize_file(ncid) + +end subroutine def_fill_dimvars + +!----------------------------------------------------------------------- +end program aether_to_dart + diff --git a/models/aether_lat-lon/dart_to_aether.f90 b/models/aether_lat-lon/dart_to_aether.f90 new file mode 100644 index 0000000000..c8aab30dcd --- /dev/null +++ b/models/aether_lat-lon/dart_to_aether.f90 @@ -0,0 +1,395 @@ +! DART software - Copyright 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 +! + +program dart_to_aether + +!---------------------------------------------------------------------- +! purpose: Transform a DART filter_output.nc into the Aether model restarts. +! +! method: Read DART state netcdf file and overwrite values in Aether restart files. +! +! this version assumes that the DART grid is global and the data needs to be +! blocked into one block per Aether mpi task. there is a different converter +! for when Aether only needs a single input/output file. +! +!---------------------------------------------------------------------- + +use types_mod, only : r4, MISSING_I, MISSING_R4, vtablenamelength + +use utilities_mod, only : & + finalize_utilities, error_handler, E_ERR, E_MSG, E_WARN, & + initialize_utilities, do_output + +use default_model_mod, only : write_model_time + +use transform_state_mod, only : & + debug, aether_restart_dirname, nblocks_lat, & + nblocks_lon, nghost, nlat, nlon, nlev, & + nx_per_block, ny_per_block, nz_per_block, & + nvar_ion, nvar_neutral, VT_ORIGININDX, VT_VARNAMEINDX, & + block_file_name, open_block_file, aether_name_to_dart, & + variables, purge_chars, static_init_blocks + +use netcdf_utilities_mod, only : & + nc_open_file_readonly, nc_close_file, & + nc_begin_define_mode, nc_end_define_mode, & + nc_define_dimension, & + nc_add_global_attribute, nc_add_global_creation_time, & + nc_get_attribute_from_variable, nc_add_attribute_to_variable, & + nc_define_real_variable, nc_define_real_scalar, & + nc_get_variable, nc_put_variable, nc_variable_exists, & + nc_synchronize_file, NF90_FILL_REAL + +implicit none + +!---------------------------------------------------------------------- +! global storage +!---------------------------------------------------------------------- + +integer :: member = MISSING_I, & + num_args, ncid +character(len=3) :: char_mem +character(len=31) :: filter_io_root = 'filter_input' +character(len=64) :: filter_io_file = '' +character(len=512) :: error_string_1, error_string_2 +character(len=31), parameter :: progname = 'dart_to_aether' +character(len=256), parameter :: source = 'aether_lat-lon/dart_to_aether.f90' + +!====================================================================== + +call initialize_utilities(progname) + +!---------------------------------------------------------------------- +! Get the ensemble member +!---------------------------------------------------------------------- +num_args = command_argument_count() +if (num_args == 0) then + write(error_string_1,*) 'Usage: ./dart_to_aether member_number (0-based)' + call error_handler(E_ERR, progname, error_string_1) +endif + +call get_command_argument(1,char_mem) +read(char_mem,'(I3)') member + +!---------------------------------------------------------------------- +! Convert the files +!---------------------------------------------------------------------- + +call static_init_blocks(member) + +write(filter_io_file,'(2A,I0.4,A3)') trim(filter_io_root),'_',member + 1,'.nc' + +call error_handler(E_MSG, source, '', '') +write(error_string_1,'(3A)') 'Extracting fields from DART file ',trim(filter_io_file) +write(error_string_2,'(A,I3,2A)') 'into Aether restart member ',member, & + ' in directory ', trim(aether_restart_dirname) +call error_handler(E_MSG, progname, error_string_1, text2=error_string_2) +call error_handler(E_MSG, '', '') + +ncid = nc_open_file_readonly(trim(aether_restart_dirname)//'/'//trim(filter_io_file), source) + +call filter_to_restarts(ncid, member) + +!---------------------------------------------------------------------- +! Log what we think we're doing, and exit. +!---------------------------------------------------------------------- +call error_handler(E_MSG, source,'','') +write(error_string_1,'(3A)') 'Successfully converted to the Aether restart files in directory' +write(error_string_2,'(3A)') "'"//trim(aether_restart_dirname)//"'" +call error_handler(E_MSG, source, error_string_1, source, text2=error_string_2) + +call nc_close_file(ncid) + +! end - close the log, etc +call finalize_utilities() + +!----------------------------------------------------------------------- +contains +!----------------------------------------------------------------------- +! Extract (updated) variables from a filter_output.nc file +! and write to existing block restart files. + +subroutine filter_to_restarts(ncid, member) + +integer, intent(in) :: member, ncid + +real(r4), allocatable :: fulldom3d(:,:,:) +character(len=64) :: file_root +integer :: ivar +character(len=vtablenamelength) :: varname, dart_varname + +character(len=*), parameter :: routine = 'filter_to_restarts' + +! Space for full domain field (read from filter_output.nc) +! and halo around the full domain +allocate(fulldom3d(1:nlev, & + 1-nghost:nlat+nghost, & + 1-nghost:nlon+nghost)) + +! get the dirname, construct the filenames inside open_block_file + +! Not all fields have halos suitable for calculating gradients. +! These do (2023-11-8): neutrals; temperature, O, O2, N2, and the horizontal winds. +! ions; none. +! The current transform_state will fill all neutral halos anyway, +! since that's simpler and won't break the model. +! TODO: add an attribute to the variables (?) to denote whether a field +! should have its halo filled? +do ivar = 1, nvar_neutral + varname = purge_chars(trim(variables(VT_VARNAMEINDX,ivar)), '\', plus_minus=.false.) + if (debug >= 0 .and. do_output()) then + write(error_string_1,'("varname = ",A)') trim(varname) + call error_handler(E_MSG, routine, error_string_1, source) + endif + dart_varname = aether_name_to_dart(varname) + + file_root = trim(variables(VT_ORIGININDX,ivar)) + if (trim(file_root) == 'neutrals') then + ! This parameter is available through the `use netcdf` command. + fulldom3d = NF90_FILL_REAL + + call nc_get_variable(ncid, dart_varname, fulldom3d(1:nlev,1:nlat,1:nlon), & + context=routine) + ! Copy updated field values to full domain halo. + ! Block domains+halos will be easily read from this. + call add_halo_fulldom3d(fulldom3d) + + call filter_io_to_blocks(fulldom3d, varname, file_root, member) + else + write(error_string_1,'(3A)') "file_root of varname = ",trim(varname), & + ' expected to be "neutrals"' + call error_handler(E_ERR, routine, error_string_1, source) + endif + +enddo + +do ivar = nvar_neutral + 1, nvar_neutral + nvar_ion + varname = purge_chars(trim(variables(VT_VARNAMEINDX,ivar)), '\', plus_minus=.false.) + dart_varname = aether_name_to_dart(varname) + + file_root = trim(variables(VT_ORIGININDX,ivar)) + if (debug >= 0 .and. do_output()) then + write(error_string_1,'("varname, dart_varname, file_root = ",3(2x,A))') & + trim(varname), trim(dart_varname), trim(file_root) + call error_handler(E_MSG, routine, error_string_1, source) + endif + + if (trim(file_root) == 'ions') then + fulldom3d = NF90_FILL_REAL + call nc_get_variable(ncid, dart_varname, fulldom3d(1:nlev,1:nlat,1:nlon), & + context=routine) + ! 2023-11: ions do not have real or used data in their halos. + ! Make this clear by leaving the halos filled with MISSING_R4 + ! TODO: Will this be translated into NetCDF missing_value? + ! call add_halo_fulldom3d(fulldom3d) + + call filter_io_to_blocks(fulldom3d, varname, file_root, member) + + else + write(error_string_1,'(3A)') "file_root of varname = ",trim(varname), & + ' expected to be "ions"' + call error_handler(E_ERR, routine, error_string_1, source) + endif +enddo + +deallocate(fulldom3d) + +end subroutine filter_to_restarts + + +!----------------------------------------------------------------------- +! Copy updated data from the full domain into the halo regions, +! in preparation for extracting haloed blocks into the block restart files. +! First, the halos past the East and West edges are taken from the wrap-around points. +! Then, the halos beyond the edge latitudes in the North and South +! are taken by reaching over the pole to a longitude that's half way around the globe. +! This is independent of the number of blocks. + +subroutine add_halo_fulldom3d(fulldom3d) + +! Space for full domain field (read from filter_output.nc) +! and halo around the full domain +real(r4), intent(inout) :: fulldom3d(1:nz_per_block, & + 1-nghost:nlat+nghost, & + 1-nghost:nlon+nghost) + +integer :: g, i, j, haflat, haflon +real(r4), allocatable :: normed(:,:) +character(len=16) :: debug_format + +character(len=*), parameter :: routine = 'add_halo_fulldom3d' + +! An array for debugging by renormalizing an altitude of fulldom3d. +allocate(normed(1-nghost:nlat+nghost, & + 1-nghost:nlon+nghost)) + +haflat = nlat / 2 +haflon = nlon / 2 + +do g = 1,nghost + ! left; reach around the date line. + ! There's no data at the ends of the halos for this copy. + fulldom3d (:,1:nlat, 1-g) & + = fulldom3d(:,1:nlat,nlon+1-g) + + ! right + fulldom3d (:,1:nlat,nlon+g) & + = fulldom3d(:,1:nlat,g) + + ! bottom; reach over the S Pole for halo values. + ! There is data at the ends of the halos for these.) + + fulldom3d (:, 1-g ,1-nghost :haflon) & + = fulldom3d(:, g ,1-nghost+haflon:nlon) + fulldom3d (:, 1-g ,haflon+1:nlon) & + = fulldom3d(:, g , 1:haflon) + ! Last 2 (halo) points on the right edge (at the bottom) + fulldom3d (:, 1-g , nlon+1: nlon+nghost) & + = fulldom3d(:, g ,haflon+1:haflon+nghost) + + ! top + fulldom3d (:, nlat +g, 1-nghost :haflon) & + = fulldom3d(:, nlat+1-g, 1-nghost+haflon:nlon) + fulldom3d (:, nlat +g, haflon+1:nlon) & + = fulldom3d(:, nlat+1-g, 1:haflon) + ! Last 2 (halo) points on the right edge (at the top) + fulldom3d (:, nlat +g, nlon+1: nlon+nghost) & + = fulldom3d(:, nlat+1-g, haflon+1:haflon+nghost) +enddo + +if (any(fulldom3d == MISSING_R4)) then + error_string_1 = 'ERROR: some fulldom3d contain MISSING_R4 after halos' + call error_handler(E_ERR, routine, error_string_1, source) +endif + +! TODO: Keep halo corners check for future use? +! Add more robust rescaling. +! Print the 4x4 arrays (corners & middle) to see whether values are copied correctly. +! Level 44 values range from 800-eps to 805. I don't want to see the 80. +! For O+ range from 0 to 7e+11, but are close to 1.1082e+10 near the corners. +! 2023-12-20; Aaron sent new files with 54 levels. +if (debug >= 100 .and. do_output()) then + if (fulldom3d(54,10,10) > 1.e+10_r4) then + normed = fulldom3d(54,:,:) - 1.1092e+10_r4 + debug_format = '(3(4E10.4,2X))' + else if (fulldom3d(54,10,10) < 1000._r4) then + normed = fulldom3d(54,:,:) - 800._r4 + debug_format = '(3(4F10.5,2X))' + endif + + ! Debug HDF5 + write(error_string_1,'("normed_field(10,nlat+1,nlon+2) = ",3(1x,i5))') normed(nlat+1,nlon+2) + call error_handler(E_MSG, routine, error_string_1, source) + + ! 17 format debug_format + print*,'top' + do j = nlat+2, nlat-1, -1 + write(*,debug_format) (normed(j,i), i= -1, 2), & + (normed(j,i), i=haflon-1,haflon+2), & + (normed(j,i), i= nlon-1, nlon+2) + enddo + print*,'middle' + do j = haflat+2, haflat-1 , -1 + write(*,debug_format) (normed(j,i), i= -1, 2), & + (normed(j,i), i=haflon-1,haflon+2), & + (normed(j,i), i= nlon-1, nlon+2) + enddo + print*,'bottom' + do j = 2,-1, -1 + write(*,debug_format) (normed(j,i), i= -1, 2), & + (normed(j,i), i=haflon-1,haflon+2), & + (normed(j,i), i= nlon-1, nlon+2) + enddo +endif + +deallocate(normed) + +end subroutine add_halo_fulldom3d + +!----------------------------------------------------------------------- +! Transfer part of the full field into a block restart file. + +subroutine filter_io_to_blocks(fulldom3d, varname, file_root, member) + +real(r4), intent(in) :: fulldom3d(1:nz_per_block, & + 1-nghost:nlat+nghost, & + 1-nghost:nlon+nghost ) +character(len=*), intent(in) :: varname +character(len=*), intent(in) :: file_root +integer, intent(in) :: member + +! Don't collect velocity components (6 of them) +! real(r4) :: temp0d +! , temp1d(:) ? +integer :: ncid_output +integer :: ib, jb, nb +integer :: starts(3), ends(3), xcount, ycount, zcount +character(len=256) :: block_file + +character(len=*), parameter :: routine = 'filter_io_to_blocks' + +! a temp array large enough to hold any of the +! Lon,Lat or Alt array from a block plus ghost cells +! allocate(temp1d(1-nghost:max(nx_per_block,ny_per_block,nz_per_block)+nghost)) + +zcount = nz_per_block +ycount = ny_per_block + (2 * nghost) +xcount = nx_per_block + (2 * nghost) + +if (debug > 0 .and. do_output()) then + write(error_string_1,'(A,I0,A,I0,A)') 'Now putting the data for ', nblocks_lon, & + ' blocks lon by ',nblocks_lat,' blocks lat' + call error_handler(E_MSG, routine, error_string_1, source) +end if + +starts(1) = 1 +ends(1) = nz_per_block + +do jb = 1, nblocks_lat + starts(2) = (jb - 1) * ny_per_block - nghost + 1 + ends(2) = jb * ny_per_block + nghost + + do ib = 1, nblocks_lon + starts(3) = (ib - 1) * nx_per_block - nghost + 1 + ends(3) = ib * nx_per_block + nghost + + nb = (jb - 1) * nblocks_lon + ib - 1 + + write(block_file,'(A,"/",A)') trim(aether_restart_dirname), & + trim(block_file_name(trim(file_root), member, nb)) + ncid_output = open_block_file(block_file, 'readwrite') + if (.not.nc_variable_exists(ncid_output,varname)) then + write(error_string_1,'(4A)') 'variable ', varname, ' does not exist in ',block_file + call error_handler(E_ERR, routine, error_string_1, source) + endif + + if ( debug > 0 .and. do_output()) then + write(error_string_1,'(A,3(2X,i5))') "block, ib, jb = ", nb, ib, jb + call error_handler(E_MSG, routine, error_string_1, source) + write(error_string_1,'(3(A,3i5))') & + 'starts = ',starts, 'ends = ',ends, '[xyz]counts = ',xcount,ycount,zcount + call error_handler(E_MSG, routine, error_string_1, source) + endif + + call nc_put_variable(ncid_output, trim(varname), & + fulldom3d(starts(1):ends(1), starts(2):ends(2), starts(3):ends(3)), & + context=routine, nc_count=(/ zcount,ycount,xcount /) ) + + call nc_close_file(ncid_output) + + enddo +enddo + +! +! TODO: ? Add f107 and Rho to the restart files +! call read_filter_io_block0d(ncid, ivals(1), data0d) +! if (data0d < 0.0_r8) data0d = 60.0_r8 !alex +! write(ounit) data0d + +end subroutine filter_io_to_blocks + +!----------------------------------------------------------------------- +end program dart_to_aether + diff --git a/models/aether_lat-lon/model_mod.f90 b/models/aether_lat-lon/model_mod.f90 new file mode 100644 index 0000000000..87ac6c0e9c --- /dev/null +++ b/models/aether_lat-lon/model_mod.f90 @@ -0,0 +1,622 @@ +! DART software - Copyright 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 +! + +module model_mod + +!----------------------------------------------------------------------- +! +! Interface for Aether +! +!----------------------------------------------------------------------- + +use types_mod, only : & + r8, i8, MISSING_R8, vtablenamelength + +use time_manager_mod, only : & + time_type, set_time, set_calendar_type + +use location_mod, only : & + location_type, get_close_type, & + get_close_obs, get_close_state, & + is_vertical, set_location, & + VERTISHEIGHT, query_location, get_location + +use utilities_mod, only : & + open_file, close_file, & + error_handler, E_ERR, E_MSG, E_WARN, & + nmlfileunit, do_nml_file, do_nml_term, & + find_namelist_in_file, check_namelist_read, to_upper, & + find_enclosing_indices + +use obs_kind_mod, only : get_index_for_quantity + +use netcdf_utilities_mod, only : & + nc_add_global_attribute, nc_synchronize_file, & + nc_add_global_creation_time, & + nc_begin_define_mode, nc_end_define_mode, & + nc_open_file_readonly, nc_get_dimension_size, nc_create_file, & + nc_get_variable + +use quad_utils_mod, only : & + quad_interp_handle, init_quad_interp, set_quad_coords, & + quad_lon_lat_locate, quad_lon_lat_evaluate, & + GRID_QUAD_FULLY_REGULAR, QUAD_LOCATED_CELL_CENTERS + +use state_structure_mod, only : & + add_domain, get_dart_vector_index, get_domain_size, & + get_model_variable_indices, get_varid_from_kind, & + state_structure_info + +use distributed_state_mod, only : get_state + +use ensemble_manager_mod, only : ensemble_type + +! These routines are passed through from default_model_mod. +! To write model specific versions of these routines +! remove the routine from this use statement and add your code to +! this the file. +use default_model_mod, only : & + pert_model_copies, read_model_time, write_model_time, & + init_time => fail_init_time, & + init_conditions => fail_init_conditions, & + convert_vertical_obs, convert_vertical_state, adv_1step + +implicit none +private + +! routines required by DART code - will be called from filter and other DART executables. +public :: get_model_size, & + get_state_meta_data, & + model_interpolate, & + end_model, & + static_init_model, & + nc_write_model_atts, & + get_close_obs, & + get_close_state, & + pert_model_copies, & + convert_vertical_obs, & + convert_vertical_state, & + read_model_time, & + adv_1step, & + init_time, & + init_conditions, & + shortest_time_between_assimilations, & + write_model_time + +character(len=256), parameter :: source = 'aether_lat-lon/model_mod.f90' + +logical :: module_initialized = .false. +integer :: dom_id ! used to access the state structure +type(time_type) :: assimilation_time_step + +!----------------------------------------------------------------------- +! Default values for namelist +character(len=256) :: template_file = 'filter_input_0001.nc' +integer :: time_step_days = 0 +integer :: time_step_seconds = 3600 + +integer, parameter :: MAX_STATE_VARIABLES = 100 +integer, parameter :: NUM_STATE_TABLE_COLUMNS = 5 +character(len=vtablenamelength) :: variables(NUM_STATE_TABLE_COLUMNS,MAX_STATE_VARIABLES) = '' + +type :: var_type + integer :: count + character(len=64), allocatable :: names(:) + integer, allocatable :: qtys(:) + real(r8), allocatable :: clamp_values(:, :) + logical, allocatable :: updates(:) +end type var_type + +namelist /model_nml/ template_file, time_step_days, time_step_seconds, variables + +!----------------------------------------------------------------------- +! Dimensions + +character(len=4), parameter :: LEV_DIM_NAME = 'alt' +character(len=4), parameter :: LAT_DIM_NAME = 'lat' +character(len=4), parameter :: LON_DIM_NAME = 'lon' +character(len=4), parameter :: TIME_DIM_NAME = 'time' + +character(len=4), parameter :: LEV_VAR_NAME = 'alt' +character(len=4), parameter :: LAT_VAR_NAME = 'lat' +character(len=4), parameter :: LON_VAR_NAME = 'lon' +character(len=4), parameter :: TIME_VAR_NAME = 'time' + +! Filter +! To be assigned in assign_dimensions (for filter) +! or get_grid_from_blocks (aether_to_dart, dart_to_aether). +real(r8), allocatable :: levs(:), lats(:), lons(:) + +integer :: nlev, nlat, nlon +real(r8) :: lon_start, lon_delta, lat_start, lat_delta, lat_end + +!----------------------------------------------------------------------- +! to be assigned in the verify_variables subroutine + +type(quad_interp_handle) :: quad_interp + +integer, parameter :: GENERAL_ERROR_CODE = 99 +integer, parameter :: INVALID_VERT_COORD_ERROR_CODE = 15 +integer, parameter :: INVALID_LATLON_VAL_ERROR_CODE = 16 +integer, parameter :: INVALID_ALTITUDE_VAL_ERROR_CODE = 17 +integer, parameter :: UNKNOWN_OBS_QTY_ERROR_CODE = 20 + +type(time_type) :: state_time ! module-storage declaration of current model time +character(len=512) :: error_string_1, error_string_2 + +contains + +!----------------------------------------------------------------------- +! 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. + +subroutine static_init_model() + +integer :: iunit, io +type(var_type) :: var + +module_initialized = .true. + +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") + +call set_calendar_type('GREGORIAN') + +! Record the namelist values used for the run +if (do_nml_file()) write(nmlfileunit, nml=model_nml) +if (do_nml_term()) write( * , nml=model_nml) + +call assign_dimensions() + +! Dimension start and deltas needed for set_quad_coords +lon_start = lons(1) +lon_delta = lons(2) - lons(1) +lat_start = lats(1) +lat_delta = lats(2) - lats(1) + +var = assign_var(variables, MAX_STATE_VARIABLES) + +! This time is both the minimum time you can ask the model to advance +! (for models that can be advanced by filter) and it sets the assimilation +! window. All observations within +/- 1/2 this interval from the current +! model time will be assimilated. If this is not settable at runtime +! feel free to hardcode it and remove from the namelist. +assimilation_time_step = set_time(time_step_seconds, time_step_days) + +! Define which variables are in the model state +! This is using add_domain_from_file (arg list matches) +dom_id = add_domain(template_file, var%count, var%names, var%qtys, & + var%clamp_values, var%updates) + +call state_structure_info(dom_id) + + +call init_quad_interp(GRID_QUAD_FULLY_REGULAR, nlon, nlat, & + QUAD_LOCATED_CELL_CENTERS, & + global=.true., spans_lon_zero=.true., pole_wrap=.true., & + interp_handle=quad_interp) + +call set_quad_coords(quad_interp, lon_start, lon_delta, lat_start, lat_delta) + +end subroutine static_init_model + +!----------------------------------------------------------------------- +! Returns the number of items in the state vector as an integer. + +function get_model_size() + +integer(i8) :: get_model_size + +if ( .not. module_initialized ) call static_init_model + +get_model_size = get_domain_size(dom_id) + +end function get_model_size + +!----------------------------------------------------------------------- +! Use quad_utils_mod to interpalate the ensemble to the ob location. + +subroutine model_interpolate(state_handle, ens_size, location, qty, expected_obs, istatus) + +type(ensemble_type), intent(in) :: state_handle +integer, intent(in) :: ens_size +type(location_type), intent(in) :: location +integer, intent(in) :: qty +real(r8), intent(out) :: expected_obs(ens_size) +integer, intent(out) :: istatus(ens_size) + +! Local storage + +character(len=*), parameter :: routine = 'model_interpolate' + +real(r8) :: loc_array(3), llon, llat, lvert, lon_fract, lat_fract +integer :: four_lons(4), four_lats(4) +integer :: status1, which_vert, varid +real(r8) :: quad_vals(4, ens_size) + +if ( .not. module_initialized ) call static_init_model + +! Assume failure. Set return val to missing, then the code can +! just set istatus to something indicating why it failed, and return. +! If the interpolation is good, expected_obs will be set to the +! good values, and the last line here sets istatus to 0. +! make any error codes set here be in the 10s + +expected_obs = MISSING_R8 ! the DART bad value flag +istatus = GENERAL_ERROR_CODE ! unknown error + +! Get the individual locations values + +loc_array = get_location(location) +llon = loc_array(1) +llat = loc_array(2) +lvert = loc_array(3) +which_vert = nint(query_location(location)) + +! Only height and level for vertical location type is supported at this point +if (.not. is_vertical(location, "HEIGHT") .and. .not. is_vertical(location, "LEVEL")) THEN + istatus = INVALID_VERT_COORD_ERROR_CODE + return +endif + +! See if the state contains the obs quantity +varid = get_varid_from_kind(dom_id, qty) + +if (varid > 0) then + istatus = 0 +else + istatus = UNKNOWN_OBS_QTY_ERROR_CODE +endif + +! get the indices for the 4 corners of the quad in the horizontal, plus +! the fraction across the quad for the obs location +call quad_lon_lat_locate(quad_interp, llon, llat, & + four_lons, four_lats, lon_fract, lat_fract, status1) +if (status1 /= 0) then + istatus(:) = INVALID_LATLON_VAL_ERROR_CODE ! cannot locate enclosing horizontal quad + return +endif + +call get_quad_vals(state_handle, ens_size, varid, four_lons, four_lats, & + loc_array, which_vert, quad_vals, istatus) +if (any(istatus /= 0)) return + +! do the horizontal interpolation for each ensemble member +call quad_lon_lat_evaluate(quad_interp, lon_fract, lat_fract, ens_size, & + quad_vals, expected_obs, istatus) + +! All good. +istatus(:) = 0 + +end subroutine model_interpolate + +!----------------------------------------------------------------------- +! Returns the smallest increment in time that the model is capable +! of advancing the state in a given implementation, or the shortest +! time you want the model to advance between assimilations. + +function shortest_time_between_assimilations() + +type(time_type) :: shortest_time_between_assimilations + +if ( .not. module_initialized ) call static_init_model + +shortest_time_between_assimilations = assimilation_time_step + +end function shortest_time_between_assimilations + + + +!----------------------------------------------------------------------- +! Given an integer index into the state vector, returns the +! associated location and optionally the physical quantity. + +subroutine get_state_meta_data(index_in, location, qty) + +integer(i8), intent(in) :: index_in +type(location_type), intent(out) :: location +integer, optional , intent(out) :: qty + +! Local variables + +integer :: lat_index, lon_index, lev_index +integer :: my_var_id, my_qty + +if ( .not. module_initialized ) call static_init_model + +! Restart data is ordered (lev,lat,lon) (translated from C to fortran). +call get_model_variable_indices(index_in, lev_index, lat_index, lon_index, & + var_id=my_var_id, kind_index=my_qty) + +! should be set to the actual location using set_location() +location = set_location(lons(lon_index), lats(lat_index), levs(lev_index), VERTISHEIGHT) + +! should be set to the physical quantity, e.g. QTY_TEMPERATURE +if (present(qty)) qty = my_qty + +end subroutine get_state_meta_data + +!----------------------------------------------------------------------- +! Does any shutdown and clean-up needed for model. Can be a NULL +! INTERFACE if the model has no need to clean up storage, etc. + +subroutine end_model() + +end subroutine end_model + +!----------------------------------------------------------------------- +! write any additional attributes to the output and diagnostic files + +subroutine nc_write_model_atts(ncid, domain_id) + +integer, intent(in) :: ncid ! netCDF file identifier +integer, intent(in) :: domain_id + +character(len=*), parameter :: routine = 'nc_write_model_atts' + +if ( .not. module_initialized ) call static_init_model + +! It is already in define mode from nc_create_file. +! OR NOT, if called by create_and_open_state_output +call nc_begin_define_mode(ncid) + +! nc_write_model_atts is called by create_and_open_state_output, +! which calls nf90_enddef before it. +call nc_add_global_creation_time(ncid, routine) + +call nc_add_global_attribute(ncid, "model_source", source, routine) +call nc_add_global_attribute(ncid, "model", "aether", routine) + +call nc_end_define_mode(ncid) + +! Flush the buffer and leave netCDF file open +call nc_synchronize_file(ncid) + +end subroutine nc_write_model_atts + +!----------------------------------------------------------------------- +! Read dimension information from the template file and use +! it to assign values to variables. + +subroutine assign_dimensions() + +integer :: ncid +character(len=24), parameter :: routine = 'assign_dimensions' + +call error_handler(E_MSG, routine, 'reading filter input ['//trim(template_file)//']') + +ncid = nc_open_file_readonly(template_file, routine) + +! levels +nlev = nc_get_dimension_size(ncid, trim(LEV_DIM_NAME), routine) +allocate(levs(nlev)) +call nc_get_variable(ncid, trim(LEV_VAR_NAME), levs, routine) + +! latitiude +nlat = nc_get_dimension_size(ncid, trim(LAT_DIM_NAME), routine) +allocate(lats(nlat)) +call nc_get_variable(ncid, trim(LAT_VAR_NAME), lats, routine) + +! longitude +nlon = nc_get_dimension_size(ncid, trim(LON_DIM_NAME), routine) +allocate(lons(nlon)) +call nc_get_variable(ncid, trim(LON_VAR_NAME), lons, routine) + +end subroutine assign_dimensions + +!----------------------------------------------------------------------- +! Parse the table of variables characteristics into arrays for easier access. + +function assign_var(variables, MAX_STATE_VARIABLES) result(var) + +character(len=vtablenamelength), intent(in) :: variables(:, :) +integer, intent(in) :: MAX_STATE_VARIABLES + +type(var_type) :: var +integer :: ivar +character(len=vtablenamelength) :: table_entry + +!----------------------------------------------------------------------- +! Codes for interpreting the NUM_STATE_TABLE_COLUMNS of the variables table +integer, parameter :: NAME_INDEX = 1 ! ... variable name +integer, parameter :: QTY_INDEX = 2 ! ... DART qty +integer, parameter :: MIN_VAL_INDEX = 3 ! ... minimum value if any +integer, parameter :: MAX_VAL_INDEX = 4 ! ... maximum value if any +integer, parameter :: UPDATE_INDEX = 5 ! ... update (state) or not + +! Loop through the variables array to get the actual count of the number of variables +do ivar = 1, MAX_STATE_VARIABLES + ! If the element is an empty string, the loop has exceeded the extent of the variables + if (variables(1, ivar) == '') then + var%count = ivar-1 + exit + endif +enddo + +! Allocate the arrays in the var derived type +allocate(var%names(var%count), var%qtys(var%count), var%clamp_values(var%count, 2), var%updates(var%count)) + +do ivar = 1, var%count + + var%names(ivar) = trim(variables(NAME_INDEX, ivar)) + + table_entry = variables(QTY_INDEX, ivar) + call to_upper(table_entry) + + var%qtys(ivar) = get_index_for_quantity(table_entry) + + if (variables(MIN_VAL_INDEX, ivar) /= 'NA') then + read(variables(MIN_VAL_INDEX, ivar), '(d16.8)') var%clamp_values(ivar,1) + else + var%clamp_values(ivar,1) = MISSING_R8 + endif + + if (variables(MAX_VAL_INDEX, ivar) /= 'NA') then + read(variables(MAX_VAL_INDEX, ivar), '(d16.8)') var%clamp_values(ivar,2) + else + var%clamp_values(ivar,2) = MISSING_R8 + endif + + table_entry = variables(UPDATE_INDEX, ivar) + call to_upper(table_entry) + + if (table_entry == 'UPDATE') then + var%updates(ivar) = .true. + else + var%updates(ivar) = .false. + endif + +enddo + +end function assign_var + +!----------------------------------------------------------------------- +! Extract state values needed by the interpolation from all ensemble members. + +subroutine get_quad_vals(state_handle, ens_size, varid, four_lons, four_lats, & + lon_lat_vert, which_vert, quad_vals, istatus) + +type(ensemble_type), intent(in) :: state_handle +integer, intent(in) :: ens_size +integer, intent(in) :: varid +integer, intent(in) :: four_lons(4), four_lats(4) +real(r8), intent(in) :: lon_lat_vert(3) +integer, intent(in) :: which_vert +real(r8), intent(out) :: quad_vals(4, ens_size) +integer, intent(out) :: istatus(ens_size) + +integer :: lev1, lev2, stat +real(r8) :: vert_val, vert_fract +character(len=512) :: error_string_1 + +character(len=*), parameter :: routine = 'get_quad_vals' + +quad_vals(:,:) = MISSING_R8 +istatus(:) = GENERAL_ERROR_CODE + +vert_val = lon_lat_vert(3) + +if ( which_vert == VERTISHEIGHT ) then + call find_enclosing_indices(nlev, levs(:), vert_val, lev1, lev2, & + vert_fract, stat, log_scale = .false.) + + if (stat /= 0) then + istatus = INVALID_ALTITUDE_VAL_ERROR_CODE + end if +else + istatus(:) = INVALID_VERT_COORD_ERROR_CODE + write(error_string_1, *) 'unsupported vertical type: ', which_vert + call error_handler(E_ERR, routine, error_string_1, source) +endif + +! we have all the indices and fractions we could ever want. +! now get the data values at the bottom levels, the top levels, +! and do vertical interpolation to get the 4 values in the columns. +! the final horizontal interpolation will happen later. + +if (varid > 0) then + + call get_four_state_values(state_handle, ens_size, four_lons, four_lats, & + lev1, lev2, vert_fract, varid, quad_vals, istatus) +else + write(error_string_1, *) 'unsupported variable: ', varid + call error_handler(E_ERR, routine, error_string_1, source) +endif + +if (any(istatus /= 0)) return + +! when you get here, istatus() was set either by passing it to a +! subroutine, or setting it explicitly here. +end subroutine get_quad_vals + +!----------------------------------------------------------------------- +! interpolate in the vertical between 2 arrays of items. + +! vert_fracts: 0 is 100% of the first level and +! 1 is 100% of the second level + +subroutine vert_interp(nitems, levs1, levs2, vert_fract, out_vals) + +integer, intent(in) :: nitems +real(r8), intent(in) :: levs1(nitems) +real(r8), intent(in) :: levs2(nitems) +real(r8), intent(in) :: vert_fract +real(r8), intent(out) :: out_vals(nitems) + +out_vals(:) = (levs1(:) * (1.0_r8 - vert_fract)) + & + (levs2(:) * vert_fract ) + +end subroutine vert_interp + +!----------------------------------------------------------------------- +! Extract the state values at the corners of the 2 quads used for interpolation. + +subroutine get_four_state_values(state_handle, ens_size, four_lons, four_lats, & + lev1, lev2, vert_fract, varid, quad_vals, istatus) + +type(ensemble_type), intent(in) :: state_handle +integer, intent(in) :: ens_size +integer, intent(in) :: four_lons(4), four_lats(4) +integer, intent(in) :: lev1, lev2 +real(r8), intent(in) :: vert_fract +integer, intent(in) :: varid +real(r8), intent(out) :: quad_vals(4, ens_size) !< array of interpolated values +integer, intent(out) :: istatus(ens_size) + +integer :: icorner +integer(i8) :: state_indx +real(r8) :: vals1(ens_size), vals2(ens_size) +real(r8) :: qvals(ens_size) + +character(len=*), parameter :: routine = 'get_four_state_values:' + +do icorner = 1, 4 + + ! Most rapidly varying dim must be first + state_indx = get_dart_vector_index(lev1, four_lats(icorner), & + four_lons(icorner), dom_id, varid) + + if (state_indx < 0) then + write(error_string_1,'(A)') 'Could not find dart state index from ' + write(error_string_2,'(A,3F15.4)') 'lon, lat, and lev1 index :', & + four_lons(icorner), four_lats(icorner), lev1 + call error_handler(E_ERR, routine, error_string_1, source, & + text2=error_string_2) + return + endif + + vals1(:) = get_state(state_indx, state_handle) ! all the ensemble members for level (i) + + state_indx = get_dart_vector_index(lev2, four_lats(icorner), & + four_lons(icorner), dom_id, varid) + + if (state_indx < 0) then + write(error_string_1,'(A)') 'Could not find dart state index from ' + write(error_string_2,'(A,3F15.4)') 'lon, lat, and lev2 index :', & + four_lons(icorner), four_lats(icorner), lev2 + call error_handler(E_ERR, routine, error_string_1, source, & + text2=error_string_2) + return + endif + + vals2(:) = get_state(state_indx, state_handle) ! all the ensemble members for level (i) + + ! if directly using quad_vals here, it would create a temporary array and give a warning + call vert_interp(ens_size, vals1, vals2, vert_fract, qvals) + quad_vals(icorner, :) = qvals +enddo + +istatus = 0 + +end subroutine get_four_state_values + +!----------------------------------------------------------------------- +! End of model_mod +!----------------------------------------------------------------------- +end module model_mod + diff --git a/models/aether_lat-lon/model_mod.nml b/models/aether_lat-lon/model_mod.nml new file mode 100644 index 0000000000..f68b1b0244 --- /dev/null +++ b/models/aether_lat-lon/model_mod.nml @@ -0,0 +1,9 @@ +&model_nml + template_file = 'if other than filter_input_0001.nc' + debug = 0 + variables = 'Temperature', 'QTY_TEMPERATURE', '1000.0', 'NA', 'UPDATE', + 'Opos', 'QTY_DENSITY_ION_OP', 'NA', 'NA', 'UPDATE' + time_step_days = 0 + time_step_seconds = 3600 + / + diff --git a/models/aether_lat-lon/readme.rst b/models/aether_lat-lon/readme.rst new file mode 100644 index 0000000000..4b146b9073 --- /dev/null +++ b/models/aether_lat-lon/readme.rst @@ -0,0 +1,230 @@ +Aether Rectangular Grid Interface +================================= + +Overview +-------- + +The Aether ("eether") space weather model can be implemented +on a logically rectangular grid "lat-lon", or on a cubed-sphere grid. +This is the interface to the lat-lon version. +The model code is available on +`GitHub `_ . + +Aether writes history and restart files. +The restart fields are divided among 2 types of files: neutrals and ions. +They are further divided into "blocks", which are subdomains of the globe. +The numbering of blocks starts in the southwest corner of the lat-lon grid +and goes east first, then to the west end of the next row north, +and ends in the northeast corner. +Each block has a halo around it filled with field values from neighboring blocks. +All of these need to be combined to make a single state vector for filter. +There's a unique set of these files for each member. +The restart file names reflect this information :: + + {neutrals,ions}_mMMMM_gBBBB.nc + MMMM = ensemble member (0-based) + BBBB = block number (0-based) + +The restart files do not have grid information in them. +Grid information must be read from :: + + grid_gBBBB.nc + +Programs ``aether_to_dart`` and ``dart_to_aether`` read the same namelist; +``transform_state_nml``. +The fields chosen to be part of the model state are specified in 'variables'. +``Aether_to_dart`` will read the specified fields, from all the restarts +for a member plus grid files, and repackage them into an ensemble state vector file +(filter_input.nc). Filter_input.nc has a single domain and no halos. +The field names will be transformed into CF-compliant names in filter_input.nc. + +``Filter`` will read a list of variables from ``model_nml`` (not ``transform_state_nml``), +then read the ensemble of filter_input.nc files, assimilate, +and write an ensemble of filter_output.nc files. + +``Dart_to_aether`` will convert the fields' names to the CF-compliant filter names, +find those names in filter_output.nc, extract the updated field data, +and overwrite those fields in the appropriate Aether restart files. + +Namelists +--------- + +- The namelists are read from the file ``input.nml``. +- Namelists start with an ampersand '&' and terminate with a slash '/'. +- Character strings that contain a '/' must be enclosed in quotes + to prevent them from prematurely terminating the namelist. + +transform_state_nml +................... + + aether_restart_dirname + The directory where the Aether restart files reside, + and will be transformed (the "run" directory). + + nblocks_lon, nblocks_lat, nblocks_lev + Number of Aether domain "blocks" in the longitudinal, latitudinal, + and vertical directions. Vertical is always 1 (2024-2). + The total number of blocks (nblocks_lon x nblocks_lat x nblocks_lev) + is defined by the number of processors used by Aether. + + variables + The Aether fields to be transformed into a model state are specified + in the 'variables' namelist variable in transform_state_nml. + The following information must be provided for each field + + 1) Aether field name + 2) which file contains the field ("neutrals" or "ions") + + Aether field names are not CF-compliant and are translated + to CF-compliant forms by aether_to_dart. + + In ``transform_state_nml`` there is no association of DART "quantities" + (QTY\_\*) with fields. + A subset of the transformed variables to be included in the model state + is specified in :ref:`model_nml:variables`, using the CF-compliant names. + That is where the associations with QTYs are made. + See the :ref:`QTY` section, below. + + The neutrals restart files contain the following fields. + The most important fields are **noted in bold text** + + | **Temperature**, **velocity_east**, **velocity_north**, + | velocity_up, N, O2, N2, NO, He, N_2D, N_2P, H, O_1D, CO2 + + Similarly for the ions restart files + + | **O+**, **O+_2D**, **O+_2P**, **O2+**, **N2+**, NO+, N+, He+, + | Temperature_bulk_ion, Temperature_electron + + In addition, there are 7 (independent) fields associated with *each* ion density + :: + + - Temperature\ \(O+\) + - velocity_parallel_east\ \(O+\) + - velocity_parallel_north\ \(O+\) + - velocity_parallel_up\ \(O+\) + - velocity_perp_east\ \(O+\) + - velocity_perp_north\ \(O+\) + - velocity_perp_up\ \(O+\) + +.. WARNING:: + As of this writing (2024-1-30) the electron density and solar radiation + parameter "f10.7" are not available through the restart files, + even though electron temperature is. + They may be available in the history files. + + +.. _model_nml: + +model_nml +......... + +template_file + = 'filter_input_0001.nc' is the default + +variables + Each field to be included in the state vector requires 5 descriptors: + + 1) field name (transformed to CF-compliant) + #) DART "quantity" to be associated with the field + #) min value + #) max value + #) update the field in the restart file? {UPDATE,NO_COPY_BACK} + + The field names listed in 'variables' must be the *transformed* names, + as found in the filter_input.nc files (see :ref:`Usage`). + In general the transformation does the following + + - Remove all '\\', '(', and ')' + - Replace blanks with underscores + - Replace '+' with 'pos' and '-' with 'neg' + - For ions, move the ion name from the end to the beginning. + + For example 'velocity_parallel_east\\ \\(O+_2D\\)' becomes 'Opos_2D_velocity_parallel_east'. + +.. _QTY: + + The DART QTY associated with each field is an open question, + depending on the forward operators required for the available observations + and on the scientific objective. The default choices are not necessarily correct + for your assimilation. For the fields identified as most important + in early Aether assimilation experiments, these are the defaults: + +============== ==================== +variables quantity (kind) +============== ==================== +Temperature QTY_TEMPERATURE +velocity_east QTY_U_WIND_COMPONENT +velocity_north QTY_V_WIND_COMPONENT +Opos QTY_DENSITY_ION_OP +O2pos QTY_DENSITY_ION_O2P +N2pos QTY_DENSITY_ION_N2P +O2pos_2D QTY_DENSITY_ION_O2DP +O2pos_2P QTY_DENSITY_ION_O2PP +============== ==================== + + Some fields could have one of several QTYs associated with them. + For example, the field 'Opos_velocity_parallel_up' + could potentially have these existing QTYs associated with it:: + + - QTY_VELOCITY_W + - QTY_VELOCITY_W_ION + - QTY_VERTICAL_VELOCITY + + It's possible that several fields could have the same QTY. + A third possibility is that the experiment may require the creation of a new QTY. + The example above may require something like QTY_VEL_PARALLEL_VERT_OP. + +.. WARNING:: + The size of these parameters may be limited to 31 characters (``types_mod.f90``) + +time_step_days, time_step_seconds + = 0, 3600 The hindcast period between assimilations. + +.. _Usage: + +Usage +----- + +The workflow and scripting for fully cycling assimilation +(ensemble hindcast, then assimilation, repeat as needed) +has not been defined yet for Aether (2024-2), +but we expect that all of the DART executables will be in a directory +which is defined in the script. +So the script will be able to run the programs using a full pathname. +In addition, all of the Aether restart files will be in a "run" directory, +which has plenty of space for the data. +The DART executables will be run in this directory using their full pathnames. + +To run a more limited test (no assimilation), +which is just the transformation of files for a member (0) +use the following csh commands, or equivalents in your preferred languange. +These build the ``aether_to_dart`` and ``dart_to_aether`` executables +in $DART/models/aether_lat-lon/work directory. +Also in that directory, edit input.nml to set ``transform_state_nml:`` ``aether_restart_dirname`` +to be the full pathname of the directory where the Aether restart and grid files are. + +:: + +> set exec_dir = $DART/models/aether_lat-lon/work +> cd $exec_dir +> ./quick_build.sh +> cd {aether_restart_dirname} +> mkdir Orig +> cp *m0000* Orig/ +> cp ${exec_dir}/input.nml . +> ${exec_dir}/aether_to_dart 0 +> cp filter_input_0001.nc filter_output_0001.nc +> ${exec_dir}/dart_to_aether 0 + +| Compare the modified Aether restart files with those in Orig. +| The filter\_ files will contain the CF-compliant field names + which must be used in ``model_nml:variables``. + +.. NOTE:: + Some halo parts may have no data in them because Aether currently (2024-2) + does not use those regions. +.. WARNING:: + The restart files have dimensions ordered such that common viewing tools + (e.g. ncview) may display the pictures transposed from what is expected. + diff --git a/models/aether_lat-lon/transform_state.nml b/models/aether_lat-lon/transform_state.nml new file mode 100644 index 0000000000..85275b0d88 --- /dev/null +++ b/models/aether_lat-lon/transform_state.nml @@ -0,0 +1,11 @@ +&transform_state_nml + aether_restart_dirname = + '/Users/raeder/DAI/Manhattan/models/aether_lat-lon/testdata4' + variables = + 'Temperature', 'neutrals', + 'O+', 'ions', + nblocks_lon = 2 + nblocks_lat = 2 + nblocks_lev = 1 + debug = 0 + / diff --git a/models/aether_lat-lon/transform_state_mod.f90 b/models/aether_lat-lon/transform_state_mod.f90 new file mode 100644 index 0000000000..382fc1599a --- /dev/null +++ b/models/aether_lat-lon/transform_state_mod.f90 @@ -0,0 +1,623 @@ +! DART software - Copyright 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 +! + +module transform_state_mod + +!----------------------------------------------------------------------- +! +! Routines used by aether_to_dart and dart_to_aether +! +!----------------------------------------------------------------------- + +use types_mod, only : & + r4, r8, MISSING_R4, MISSING_R8, vtablenamelength, MISSING_I, RAD2DEG + +use time_manager_mod, only : & + time_type, set_calendar_type, set_time, get_time, set_date, & + print_date, print_time + + +use utilities_mod, only : & + open_file, close_file, file_exist, & + error_handler, E_ERR, E_MSG, E_WARN, & + nmlfileunit, do_output, do_nml_file, do_nml_term, & + find_namelist_in_file, check_namelist_read + +use netcdf_utilities_mod, only : & + nc_open_file_readonly, nc_open_file_readwrite, nc_create_file, & + nc_get_dimension_size, nc_get_variable, & + nc_close_file + +implicit none +private + +public :: static_init_blocks, & + state_time, & + block_file_name, open_block_file, aether_name_to_dart, & + nblocks_lon, nblocks_lat, nblocks_lev, & + lons, lats, levs, & + nlon, nlat, nlev, & + nx_per_block, ny_per_block, nz_per_block, nghost, & + variables, VT_ORIGININDX, VT_VARNAMEINDX, & + nvar, nvar_neutral, nvar_ion, & + aether_restart_dirname, & + purge_chars, debug + +character(len=256), parameter :: source = 'aether_lat-lon/transform_state_mod.f90' + +logical :: module_initialized = .false. + +!----------------------------------------------------------------------- +! namelist parameters with default values. +!----------------------------------------------------------------------- + +character(len=256) :: aether_restart_dirname = '.' +! An ensemble of file names is created using this root and $member in it, + +integer, parameter :: MAX_STATE_VARIABLES = 100 +integer, parameter :: NUM_STATE_TABLE_COLUMNS = 2 +character(len=vtablenamelength) :: variables(NUM_STATE_TABLE_COLUMNS,MAX_STATE_VARIABLES) = ' ' + +! number of blocks along each dim +integer :: nblocks_lon=MISSING_I, nblocks_lat=MISSING_I, nblocks_lev=MISSING_I +! These are not used in DA, and lon_start is used only for 1D modeling +! real(r8) :: lat_start =MISSING_I, lat_end =MISSING_I, lon_start=MISSING_I + +integer :: debug = 0 + +namelist /transform_state_nml/ aether_restart_dirname, variables, debug, & + nblocks_lon, nblocks_lat, nblocks_lev + +!----------------------------------------------------------------------- +! Dimensions + +! To be assigned get_grid_from_blocks (aether_to_dart, dart_to_aether). +integer :: nlev, nlat, nlon +real(r8), allocatable :: levs(:), lats(:), lons(:) + +! Aether block parameters (nblocks_{lon,lat,lev} are read from a namelist) +integer :: nx_per_block, ny_per_block, nz_per_block + +integer, parameter :: nghost = 2 ! number of ghost cells on all edges + +!----------------------------------------------------------------------- +! Codes for interpreting the NUM_STATE_TABLE_COLUMNS of the variables table +! VT_ORIGININDX is used differently from the usual domains context. +integer, parameter :: VT_VARNAMEINDX = 1 ! ... variable name +integer, parameter :: VT_ORIGININDX = 2 ! file of origin + +!----------------------------------------------------------------------- +! Day 0 in Aether's calendar is (+/1 a day) -4710/11/24 0 UTC +! integer :: aether_ref_day = 2451545 ! cJULIAN2000 in Aether = day of date 2000/01/01. +character(len=32) :: calendar = 'GREGORIAN' + +! But what we care about is the ref time for the times in the files, which is 1965-1-1 00:00 +integer, dimension(:) :: aether_ref_date(5) = (/1965,1,1,0,0/) ! y,mo,d,h,m (secs assumed 0) +type(time_type) :: aether_ref_time, state_time +integer :: aether_ref_ndays, aether_ref_nsecs + +!----------------------------------------------------------------------- +! to be assigned in the verify_variables subroutine +integer :: nvar, nvar_neutral, nvar_ion + +!----------------------------------------------------------------------- +character(len=512) :: error_string_1, error_string_2 + +contains + +!----------------------------------------------------------------------- +! Like static_init_model, but for aether_to_dart and dart_to_aether. +! Read the namelist, +! parse the 'variables' table, +! get the Aether grid information +! convert the Aether time into a DART time. + +subroutine static_init_blocks(member) + +integer, intent(in) :: member + +character(len=128) :: aether_filename +integer :: iunit, io + +character(len=*), parameter :: routine = 'static_init_blocks' + +if (module_initialized) return ! only need to do this once + +! This prevents subroutines called from here from calling static_init_mod. +module_initialized = .true. + +!------------------ +! Read the namelist + +call find_namelist_in_file("input.nml", 'transform_state_nml', iunit) +read(iunit, nml = transform_state_nml, iostat = io) +! Record the namelist values used for the run +if (do_nml_file()) write(nmlfileunit, nml=transform_state_nml) +if (do_nml_term()) write( * , nml=transform_state_nml) +call check_namelist_read(iunit, io, 'transform_state_nml') ! closes, too. + + +! error-check, convert namelist input 'variables' to global variables. +call verify_variables(variables) + +! Aether uses Julian time internally, andor a Julian calendar +! (days from the start of the calendar), depending on the context) +call set_calendar_type( calendar ) + +!-------------------------------- +! 1) get grid dimensions +! 2) allocate space for the grids +! 3) read them from the block restart files, could be stretched ... +! Opens and closes the grid block file, but not the filter netcdf file. +call get_grid_from_blocks() + +if( debug > 0 ) then + write(error_string_1,'(A,3I5)') 'grid dims are ', nlon, nlat, nlev + call error_handler(E_MSG, routine, error_string_1, source) +endif + +! Convert the Aether reference date (not calendar day = 0 date) +! to the days and seconds of the calendar set in model_mod_nml. +aether_ref_time = set_date(aether_ref_date(1), aether_ref_date(2), aether_ref_date(3), & + aether_ref_date(4), aether_ref_date(5)) +call get_time(aether_ref_time, aether_ref_nsecs, aether_ref_ndays) + +! Get the model time from a restart file. +aether_filename = block_file_name(variables(VT_ORIGININDX,1), member, 0) +state_time = read_aether_time(trim(aether_restart_dirname)//'/'//trim(aether_filename)) + +if ( debug > 0 ) then + write(error_string_1,'("grid: nlon, nlat, nlev =",3(1x,i5))') nlon, nlat, nlev + call error_handler(E_MSG, routine, error_string_1, source) +endif + +end subroutine static_init_blocks + +!----------------------------------------------------------------------- +! Parse the table of variables' characteristics. + +subroutine verify_variables(variables) + +character(len=*), intent(in) :: variables(:,:) + +character(len=vtablenamelength) :: varname, rootstr +integer :: i + +character(len=*), parameter :: routine = 'verify_variables' + +nvar = 0 +MY_LOOP : do i = 1, size(variables,2) + + varname = variables(VT_VARNAMEINDX,i) + rootstr = variables(VT_ORIGININDX,i) + + if ( varname == ' ' .and. rootstr == ' ' ) exit MY_LOOP ! Found end of list. + + if ( varname == ' ' .or. rootstr == ' ' ) then + error_string_1 = 'variable list not fully specified' + call error_handler(E_ERR, routine, error_string_1, source) + endif + + if (i > 1) then + if (variables(VT_ORIGININDX,i-1) == 'ions' .and. rootstr /= 'ions' ) then + write(error_string_1,'(A,I1,A)') ' File type (',i, & + ') in transform_state_nml:variables is out of order or invalid.' + call error_handler(E_ERR, routine, error_string_1, source) + endif + endif + + ! The internal DART routines check if the variable name is valid. + + ! All good to here - fill the output variables + + nvar = nvar + 1 + if (variables(VT_ORIGININDX,i) == 'neutrals') nvar_neutral = nvar_neutral + 1 + if (variables(VT_ORIGININDX,i) == 'ions') nvar_ion = nvar_ion + 1 + + +enddo MY_LOOP + +if (nvar == MAX_STATE_VARIABLES) then + error_string_1 = 'WARNING: you may need to increase "MAX_STATE_VARIABLES"' + write(error_string_2,'(''you have specified at least '',i4,'' perhaps more.'')') nvar + call error_handler(E_MSG, routine, error_string_1, source, text2=error_string_2) +endif + +end subroutine verify_variables + +!----------------------------------------------------------------------- +! ? Will this need to open the grid_{below,corners,down,left} filetypes? +! This code can handle it; a longer filetype passed in, and no member. +! ? Aether output files? + +function block_file_name(filetype, memnum, blocknum) + +character(len=*), intent(in) :: filetype ! one of {grid,ions,neutrals} +integer, intent(in) :: blocknum +integer, intent(in) :: memnum +character(len=128) :: block_file_name + +character(len=*), parameter :: routine = 'block_file_name' + +block_file_name = trim(filetype) +if (memnum >= 0) write(block_file_name, '(A,A2,I0.4)') trim(block_file_name), '_m', memnum +if (blocknum >= 0) write(block_file_name, '(A,A2,I0.4)') trim(block_file_name), '_g', blocknum +block_file_name = trim(block_file_name)//'.nc' +if ( debug > 0 ) then + write(error_string_1,'("filename, memnum, blocknum = ",A,2(1x,i5))') & + trim(block_file_name), memnum, blocknum + call error_handler(E_MSG, routine, error_string_1, source) +endif + +end function block_file_name + +!----------------------------------------------------------------------- +! Read block grid values (2D arrays) from a grid NetCDF file. +! Allocate and fill the full-domain 1-D dimension arrays (lon, lat, levs) + +! This routine needs: +! +! 1. A base dirname for the restart files (aether_restart_dirname). +! The filenames have the format 'dirname/{neutrals,ions}_mMMMM_gBBBB.rst' +! where BBBB is the block number, MMMM is the member number, +! and they have leading 0s. Blocks start in the +! southwest corner of the lat/lon grid and go east first, +! then to the west end of the next row north and end in the northeast corner. +! +! In the process, the routine will find: +! +! 1. The number of blocks in Lon and Lat (nblocks_lon, nblocks_lat) +! +! 2. The number of lons and lats in a single grid block (nx_per_block, ny_per_block, nz_per_block) +! +! 3. The overall grid size, {nlon,nlat,nalt} when you've read in all the blocks. +! +! 4. The number of neutral species (and probably a mapping between +! the species number and the variable name) (nvar_neutral) +! +! 5. The number of ion species (ditto - numbers <-> names) (nvar_ion) +! +! In addition to reading in the state data, it fills Longitude, Latitude, and Altitude arrays. +! This grid is orthogonal and rectangular but can have irregular spacing along +! any of the three dimensions. + +subroutine get_grid_from_blocks() + +integer :: nb, offset, ncid, nboff +integer :: starts(3), ends(3), xcount, ycount, zcount +character(len=256) :: filename +real(r4), allocatable :: temp(:,:,:) + +character(len=*), parameter :: routine = 'get_grid_from_blocks' + +! Read the x,y,z from a NetCDF block file(s), +! in order to calculate the n[xyz]_per_block dimensions. +! grid_g0000.nc looks like a worthy candidate, but a restart could be used. +write (filename,'(2A)') trim(aether_restart_dirname),'/grid_g0000.nc' +ncid = nc_open_file_readonly(filename, routine) + +! The grid (and restart) file variables have halos, so strip them off +! to get the number of actual data values in each dimension of the block. +nx_per_block = nc_get_dimension_size(ncid, 'x', routine) - (2 * nghost) +ny_per_block = nc_get_dimension_size(ncid, 'y', routine) - (2 * nghost) +nz_per_block = nc_get_dimension_size(ncid, 'z', routine) + +nlon = nblocks_lon * nx_per_block +nlat = nblocks_lat * ny_per_block +nlev = nblocks_lev * nz_per_block + +write(error_string_1,'(3(A,I5))') 'nlon = ', nlon, ', nlat = ', nlat, ', nlev = ', nlev +call error_handler(E_MSG, routine, error_string_1, source) + +allocate( lons( nlon )) +allocate( lats( nlat )) +allocate( levs( nlev )) + +if (debug > 4) then + write(error_string_1,'(2A)') 'Successfully read Aether grid file:', trim(filename) + call error_handler(E_MSG, routine, error_string_1, source) + write(error_string_1,'(A,I5)') ' nx_per_block:', nx_per_block, & + ' ny_per_block:', ny_per_block, ' nz_per_block:', nz_per_block + call error_handler(E_MSG, routine, error_string_1, source) +endif + +! A temp array large enough to hold any of the 3D +! Lon, Lat or Alt arrays from a block plus ghost cells. +! The restart files have C-indexing (fastest changing dim is the last). +allocate(temp( 1:nz_per_block, & + 1-nghost:ny_per_block+nghost, & + 1-nghost:nx_per_block+nghost)) +temp = MISSING_R4 + +starts(1) = 1 - nghost +starts(2) = 1 - nghost +starts(3) = 1 +ends(1) = nx_per_block + nghost +ends(2) = ny_per_block + nghost +ends(3) = nz_per_block +xcount = nx_per_block + (2 * nghost) +ycount = ny_per_block + (2 * nghost) +zcount = nz_per_block +if ( debug > 0 ) then + write(error_string_1,'(2(A,3i5),A,3(1X,i5))') & + 'starts = ',starts, 'ends = ',ends, '[xyz]counts = ',xcount,ycount,zcount + call error_handler(E_MSG, routine, error_string_1, source) +endif + +! go across the south-most block row picking up all longitudes +do nb = 1, nblocks_lon + + ! filename is trimmed by passage to open_block_file + "len=*" there. + filename = trim(aether_restart_dirname)//'/'//block_file_name('grid', -1, nb-1) + ncid = open_block_file(filename, 'read') + + ! Read 3D array and extract the longitudes of the non-halo data of this block. + ! The restart files have C-indexing (fastest changing dim is the last), + ! So invert the dimension bounds. + call nc_get_variable(ncid, 'Longitude', & + temp(starts(3):ends(3), starts(2):ends(2), starts(1):ends(1)), & + context=routine, & + nc_count=(/ zcount,ycount,xcount /)) + + offset = (nx_per_block * (nb - 1)) + lons(offset+1:offset+nx_per_block) = temp(1,1,1:nx_per_block) + + call nc_close_file(ncid) +enddo + +! go up west-most block row picking up all latitudes +do nb = 1, nblocks_lat + + ! Aether's block name counter start with 0, but the lat values can come from + ! any lon=const column of blocks. + nboff = ((nb - 1) * nblocks_lon) + filename = trim(aether_restart_dirname)//'/'//block_file_name('grid', -1, nboff) + ncid = open_block_file(filename, 'read') + + call nc_get_variable(ncid, 'Latitude', & + temp(starts(3):ends(3), starts(2):ends(2), starts(1):ends(1)), & + context=routine, nc_count=(/zcount,ycount,xcount/)) + + + offset = (ny_per_block * (nb - 1)) + lats(offset+1:offset+ny_per_block) = temp(1,1:ny_per_block,1) + + call nc_close_file(ncid) +enddo + + +! this code assumes all columns share the same altitude array, +! so we can read it from the first block. +! if this is not the case, this code has to change. + +filename = trim(aether_restart_dirname)//'/'//block_file_name('grid', -1, 0) +ncid = open_block_file(filename, 'read') + +temp = MISSING_R8 +call nc_get_variable(ncid, 'Altitude', & + temp(starts(3):ends(3), starts(2):ends(2), starts(1):ends(1)), & + context=routine, nc_count=(/zcount,ycount,xcount/)) + +levs(1:nz_per_block) = temp(1:nz_per_block,1,1) + +call nc_close_file(ncid) + +deallocate(temp) + +! convert from radians into degrees +lons = lons * RAD2DEG +lats = lats * RAD2DEG + +if (debug > 4) then + print *, routine, 'All lons ', lons + print *, routine, 'All lats ', lats + print *, routine, 'All levs ', levs +endif + +if ( debug > 1 ) then ! Check dimension limits + write(error_string_1,'(A,2F15.4)') 'LON range ', minval(lons), maxval(lons) + call error_handler(E_MSG, routine, error_string_1, source) + write(error_string_1,'(A,2F15.4)') 'LAT range ', minval(lats), maxval(lats) + call error_handler(E_MSG, routine, error_string_1, source) + write(error_string_1,'(A,2F15.4)') 'ALT range ', minval(levs), maxval(levs) + call error_handler(E_MSG, routine, error_string_1, source) +endif + +end subroutine get_grid_from_blocks + +!----------------------------------------------------------------------- +! Read the Aether restart file time and convert to a DART time. + +function read_aether_time(filename) +type(time_type) :: read_aether_time +character(len=*), intent(in) :: filename + +integer :: ncid +integer :: tsimulation ! the time read from a restart file; seconds from aether_ref_date. +integer :: ndays, nsecs + +character(len=*), parameter :: routine = 'read_aether_time' + +tsimulation = MISSING_I + +ncid = open_block_file(filename, 'read') +call nc_get_variable(ncid, 'time', tsimulation, context=routine) +call nc_close_file(ncid, routine, filename) + +! Calculate the DART time of the file time. +ndays = tsimulation / 86400 +nsecs = tsimulation - (ndays * 86400) +! The ref day is not finished, but don't need to subtract 1 because +! that was accounted for in the integer calculation of ndays. +ndays = aether_ref_ndays + ndays +read_aether_time = set_time(nsecs, ndays) + +if (do_output()) & + call print_time(read_aether_time, routine//': time in restart file '//filename) +if (do_output()) & + call print_date(read_aether_time, routine//': date in restart file '//filename) + +if (debug > 8) then + write(error_string_1,'(A,I5)')'tsimulation ', tsimulation + call error_handler(E_MSG, routine, error_string_1, source) + write(error_string_1,'(A,I5)')'ndays ', ndays + call error_handler(E_MSG, routine, error_string_1, source) + write(error_string_1,'(A,I5)')'nsecs ', nsecs + call error_handler(E_MSG, routine, error_string_1, source) + + call print_date(aether_ref_time, routine//':model base date') + call print_time(aether_ref_time, routine//':model base time') +endif + +end function read_aether_time + +!----------------------------------------------------------------------- +! Convert Aether's non-CF-compliant names into CF-compliant names for filter. +! For the ions, it moves the name of the ion from the end of the variable names +! to the beginning. + +function aether_name_to_dart(varname) + +character(len=vtablenamelength), intent(in) :: varname + +character(len=vtablenamelength) :: aether_name_to_dart, aether +character(len=64) :: parts(8), var_root +integer :: char_num, first, i_parts, aether_len, end_str + +aether = trim(varname) +aether_len = len_trim(varname) +parts = '' + +! Look for the last ' '. The characters after that are the species. +! If there's no ' ', the whole string is the species. +char_num = 0 +char_num = scan(trim(aether),' ', back=.true.) +var_root = aether(char_num+1:aether_len) +! purge_chars removes unwanted [()\] +parts(1) = purge_chars( trim(var_root),')(\', plus_minus=.true.) +end_str = char_num + +! Tranform remaining pieces of varname into DART versions. +char_num = MISSING_I +first = 1 +i_parts = 2 +P_LOOP: do + ! This returns the position of the first blank *within the substring* passed in. + char_num = scan(aether(first:end_str),' ', back=.false.) + if (char_num > 0 .and. first < aether_len) then + parts(i_parts) = purge_chars(aether(first:first+char_num-1), '.)(\', plus_minus=.true.) + + first = first + char_num + i_parts = i_parts + 1 + else + exit P_LOOP + endif +enddo P_LOOP + +! Construct the DART field name from the parts +aether_name_to_dart = trim(parts(1)) +i_parts = 2 +Build : do + if (trim(parts(i_parts)) /= '') then + aether_name_to_dart = trim(aether_name_to_dart)//'_'//trim(parts(i_parts)) + i_parts = i_parts + 1 + else + exit Build + endif +enddo Build + +end function aether_name_to_dart + +!----------------------------------------------------------------------- +! Replace undesirable characters with better. + +function purge_chars(ugly_string, chars, plus_minus) + +character (len=*), intent(in) :: ugly_string, chars +logical, intent(in) :: plus_minus +character (len=64) :: purge_chars + +character (len=256) :: temp_str + +integer :: char_num, end_str, pm_num + +! Trim is not needed here +temp_str = ugly_string +end_str = len_trim(temp_str) +char_num = MISSING_I +Squeeze : do + ! Returns 0 if chars are not found + char_num = scan(temp_str, chars) + ! Need to change it to a char that won't be found by scan in the next iteration, + ! and can be easily removed. + if (char_num > 0) then + ! Squeeze out the character + temp_str(char_num:end_str-1) = temp_str(char_num+1:end_str) + temp_str(end_str:end_str) = '' +! temp_str(char_num:char_num) = ' ' + else + exit Squeeze + endif +enddo Squeeze + +! Replace + and - with pos and neg. Assume there's only 1. +temp_str = trim(adjustl(temp_str)) +end_str = len_trim(temp_str) +pm_num = scan(trim(temp_str),'+-', back=.false.) +if (pm_num == 0 .or. .not. plus_minus) then + purge_chars = trim(temp_str) +else + if (temp_str(pm_num:pm_num) == '+') then + purge_chars = temp_str(1:pm_num-1)//'pos' + else if (temp_str(pm_num:pm_num) == '-') then + purge_chars = temp_str(1:pm_num-1)//'neg' + endif + if (pm_num + 1 <= end_str) & + purge_chars = trim(purge_chars)//temp_str(pm_num+1:end_str) +endif + +end function purge_chars + +!----------------------------------------------------------------------- +! Open an Aether restart block file (neutral, ion, ...?) + +function open_block_file(filename, rw) + +! filename is trimmed by this definition +character(len=*), intent(in) :: filename +character(len=*), intent(in) :: rw ! 'read' or 'readwrite' +integer :: open_block_file + +character(len=*), parameter :: routine = 'open_block_file' + +if ( .not. file_exist(filename) ) then + write(error_string_1,'(4A)') 'cannot open file ', filename,' for ', rw + call error_handler(E_ERR, routine, error_string_1, source) +endif + +if (debug > 0) then + write(error_string_1,'(4A)') 'Opening file ', trim(filename), ' for ', rw + call error_handler(E_MSG, routine, error_string_1, source) +end if + + +if (rw == 'read') then + open_block_file = nc_open_file_readonly(filename, routine) +else if (rw == 'readwrite') then + open_block_file = nc_open_file_readwrite(filename, routine) +else + error_string_1 = ': must be called with rw={read,readwrite}, not '//rw + call error_handler(E_ERR, routine, error_string_1, source) +endif + + +if (debug > 80) then + write(error_string_1,'(4A)') 'Returned file descriptor is ', open_block_file + call error_handler(E_MSG, routine, error_string_1, source) +end if + +end function open_block_file + +end module transform_state_mod diff --git a/models/aether_lat-lon/work/filter_inputs.txt b/models/aether_lat-lon/work/filter_inputs.txt new file mode 100644 index 0000000000..c5891b8e29 --- /dev/null +++ b/models/aether_lat-lon/work/filter_inputs.txt @@ -0,0 +1,20 @@ +filter_input_0001.nc +filter_input_0002.nc +filter_input_0003.nc +filter_input_0004.nc +filter_input_0005.nc +filter_input_0006.nc +filter_input_0007.nc +filter_input_0008.nc +filter_input_0009.nc +filter_input_0010.nc +filter_input_0011.nc +filter_input_0012.nc +filter_input_0013.nc +filter_input_0014.nc +filter_input_0015.nc +filter_input_0016.nc +filter_input_0017.nc +filter_input_0018.nc +filter_input_0019.nc +filter_input_0020.nc diff --git a/models/aether_lat-lon/work/filter_outputs.txt b/models/aether_lat-lon/work/filter_outputs.txt new file mode 100644 index 0000000000..1b23ee7982 --- /dev/null +++ b/models/aether_lat-lon/work/filter_outputs.txt @@ -0,0 +1,20 @@ +filter_output_0001.nc +filter_output_0002.nc +filter_output_0003.nc +filter_output_0004.nc +filter_output_0005.nc +filter_output_0006.nc +filter_output_0007.nc +filter_output_0008.nc +filter_output_0009.nc +filter_output_0010.nc +filter_output_0011.nc +filter_output_0012.nc +filter_output_0013.nc +filter_output_0014.nc +filter_output_0015.nc +filter_output_0016.nc +filter_output_0017.nc +filter_output_0018.nc +filter_output_0019.nc +filter_output_0020.nc diff --git a/models/aether_lat-lon/work/input.nml b/models/aether_lat-lon/work/input.nml new file mode 100644 index 0000000000..c793d2cda5 --- /dev/null +++ b/models/aether_lat-lon/work/input.nml @@ -0,0 +1,327 @@ +&probit_transform_nml + / + +&algorithm_info_nml + qceff_table_filename = '' + / + +&quality_control_nml + input_qc_threshold = 3.0 + outlier_threshold = 3.0 +/ + +&state_vector_io_nml + / + +&perfect_model_obs_nml + read_input_state_from_file = .true. + single_file_in = .false. + input_state_files = 'pmo not ready for use' + init_time_days = -1 + init_time_seconds = -1 + + write_output_state_to_file = .false. + single_file_out = .false. + output_state_files = 'perfect_output_d01.nc' + output_interval = 1 + + obs_seq_in_file_name = "obs_seq.in" + obs_seq_out_file_name = "obs_seq.out" + first_obs_days = -1 + first_obs_seconds = -1 + last_obs_days = -1 + last_obs_seconds = -1 + + async = 0 + adv_ens_command = "../shell_scripts/advance_model.csh" + + trace_execution = .true. + output_timestamps = .false. + print_every_nth_obs = -1 + output_forward_op_errors = .true. + silence = .false. + / + +&filter_nml + single_file_in = .false., + input_state_files = '' + input_state_file_list = 'filter_inputs.txt' + init_time_days = 153131, + init_time_seconds = 0, + perturb_from_single_instance = .false., + perturbation_amplitude = 0.2, + + stages_to_write = 'preassim', 'analysis' + + single_file_out = .false., + output_state_files = '' + output_state_file_list = 'filter_outputs.txt' + output_interval = 1, + output_members = .true. + num_output_state_members = 20, + output_mean = .true. + output_sd = .true. + write_all_stages_at_end = .false. + compute_posterior = .true. + + ens_size = 20, + num_groups = 1, + distributed_state = .true. + + async = 4, + adv_ens_command = "./advance_model.csh", + tasks_per_model_advance = 1 + + obs_sequence_in_name = "obs_seq.out.1", + obs_sequence_out_name = "obs_seq.final", + num_output_obs_members = 20, + first_obs_days = -1, + first_obs_seconds = -1, + last_obs_days = -1, + last_obs_seconds = -1, + obs_window_days = -1, + obs_window_seconds = -1, + + inf_flavor = 5, 0, + inf_initial_from_restart = .false., .false., + inf_sd_initial_from_restart = .false., .false., + inf_deterministic = .true., .true., + inf_initial = 1.0, 1.0, + inf_lower_bound = 0.0, 1.0, + inf_upper_bound = 1000000.0, 1000000.0, + inf_damping = 1.0, 1.0, + inf_sd_initial = 0.6, 0.0, + inf_sd_lower_bound = 0.6, 0.0 + inf_sd_max_change = 1.05, 1.05, + + trace_execution = .false., + output_timestamps = .false., + output_forward_op_errors = .false., + write_obs_every_cycle = .false., + silence = .false., + + allow_missing_clm = .false. + / + + + +&ensemble_manager_nml + / + +&assim_tools_nml + cutoff = 0.2 + sort_obs_inc = .false. + spread_restoration = .false. + sampling_error_correction = .false. + adaptive_localization_threshold = -1 + output_localization_diagnostics = .false. + localization_diagnostics_file = 'localization_diagnostics' + print_every_nth_obs = 0 + / + +&transform_state_nml + aether_restart_dirname = + '/Users/raeder/DAI/Manhattan/models/aether_lat-lon/test4' + variables = + 'Temperature', 'neutrals', + 'O+', 'ions', + nblocks_lon = 2 + nblocks_lat = 2 + nblocks_lev = 1 + debug = 0 + / + +&model_nml + template_file = 'filter_input_0001.nc' + variables = 'Temperature', 'QTY_TEMPERATURE', '0.0', 'NA', 'UPDATE', + 'Opos', 'QTY_DENSITY_ION_OP', '0.0', 'NA', 'UPDATE' + time_step_days = 0 + time_step_seconds = 3600 + / + +&cov_cutoff_nml + select_localization = 1 + / + +®_factor_nml + select_regression = 1 + input_reg_file = "time_mean_reg" + save_reg_diagnostics = .false. + reg_diagnostics_file = "reg_diagnostics" + / + +&obs_sequence_nml + write_binary_obs_sequence = .false. + / + +&obs_kind_nml + assimilate_these_obs_types = 'AIRS_TEMPERATURE' + evaluate_these_obs_types = '' + / + +&location_nml + horiz_dist_only = .true. + vert_normalization_pressure = 100000.0 + vert_normalization_height = 10000.0 + vert_normalization_level = 20.0 + approximate_distance = .false. + nlon = 71 + nlat = 36 + output_box_info = .false. + / + +&preprocess_nml + overwrite_output = .true. + input_obs_qty_mod_file = '../../../assimilation_code/modules/observations/DEFAULT_obs_kind_mod.F90' + output_obs_qty_mod_file = '../../../assimilation_code/modules/observations/obs_kind_mod.f90' + input_obs_def_mod_file = '../../../observations/forward_operators/DEFAULT_obs_def_mod.F90' + output_obs_def_mod_file = '../../../observations/forward_operators/obs_def_mod.f90' + obs_type_files = '../../../observations/forward_operators/obs_def_upper_atm_mod.f90', + '../../../observations/forward_operators/obs_def_reanalysis_bufr_mod.f90', + '../../../observations/forward_operators/obs_def_altimeter_mod.f90', + '../../../observations/forward_operators/obs_def_metar_mod.f90', + '../../../observations/forward_operators/obs_def_dew_point_mod.f90', + '../../../observations/forward_operators/obs_def_rel_humidity_mod.f90', + '../../../observations/forward_operators/obs_def_gps_mod.f90', + '../../../observations/forward_operators/obs_def_vortex_mod.f90', + '../../../observations/forward_operators/obs_def_gts_mod.f90' + quantity_files = '../../../assimilation_code/modules/observations/atmosphere_quantities_mod.f90', + '../../../assimilation_code/modules/observations/space_quantities_mod.f90', + '../../../assimilation_code//modules/observations/chemistry_quantities_mod.f90' + / + +&utilities_nml + TERMLEVEL = 1 + module_details = .true. + logfilename = 'dart_log.out' + nmlfilename = 'dart_log.nml' + write_nml = 'file' + print_debug = .true. + / + +&mpi_utilities_nml + / + + +# The times in the namelist for the obs_diag program are vectors +# that follow the following sequence: +# year month day hour minute second +# max_num_bins can be used to specify a fixed number of bins +# in which case last_bin_center should be safely in the future. +# +# Acceptable latitudes range from [-90, 90] +# Acceptable longitudes range from [ 0, Inf] + +&obs_diag_nml + obs_sequence_name = 'obs_seq.final' + obs_sequence_list = '' + first_bin_center = 2005, 9, 9, 0, 0, 0 + last_bin_center = 2005, 9, 10, 0, 0, 0 + bin_separation = 0, 0, 0, 1, 0, 0 + bin_width = 0, 0, 0, 1, 0, 0 + time_to_skip = 0, 0, 0, 1, 0, 0 + max_num_bins = 1000 + trusted_obs = 'null' + Nregions = 4 + hlevel = 0, 100000, 200000, 300000, 400000, 500000, 600000, 700000, 800000, 900000, 1000000 + lonlim1 = 0.0, 0.0, 0.0, 235.0 + lonlim2 = 360.0, 360.0, 360.0, 295.0 + latlim1 = 20.0, -80.0, -20.0, 25.0 + latlim2 = 80.0, -20.0, 20.0, 55.0 + reg_names = 'Northern Hemisphere', 'Southern Hemisphere', 'Tropics', 'North America' + print_mismatched_locs = .false. + create_rank_histogram = .true. + outliers_in_histogram = .true. + use_zero_error_obs = .false. + verbose = .true. + / + +# obs_seq_to_netcdf also requires the schedule_nml. +# In this context, schedule_nml defines how many netcdf files get created. +# Each 'bin' results in an obs_epoch_xxxx.nc file. +# default is to put everything into one 'bin'. + +&obs_seq_to_netcdf_nml + obs_sequence_name = 'obs_seq.final' + obs_sequence_list = '' + append_to_netcdf = .false. + lonlim1 = 0.0 + lonlim2 = 360.0 + latlim1 = -90.0 + latlim2 = 90.0 + verbose = .false. + / + +&schedule_nml + calendar = 'Gregorian' + first_bin_start = 1601, 1, 1, 0, 0, 0 + first_bin_end = 2999, 1, 1, 0, 0, 0 + last_bin_end = 2999, 1, 1, 0, 0, 0 + bin_interval_days = 1000000 + bin_interval_seconds = 0 + max_num_bins = 1000 + print_table = .true. + / + +&obs_sequence_tool_nml + num_input_files = 1 + filename_seq = 'obs_seq.out' + filename_out = 'obs_seq.processed' + first_obs_days = -1 + first_obs_seconds = -1 + last_obs_days = -1 + last_obs_seconds = -1 + obs_types = '' + keep_types = .false. + print_only = .false. + min_lat = -90.0 + max_lat = 90.0 + min_lon = 0.0 + max_lon = 360.0 + / + + test 6 produces an exhaustive list of metadata for EVERY element in the DART state vector. + num_ens must = 1 + x_ind is for test3. The default (-1) will fail. + interp_test_dX are the model grid resolutions, + or numbers to use as such in the testing. + _d[xyz] is for cartesian grids, + _d{lon,lat,vert} is for spherical grids + interp_test_d[xyz] take precedence over d{lon,lat,vert} + all 3 must be specified. + aether (54 levels) dz ranges from ~1500 in the low levels to ~15,000 at the top. + interp_test_{lon,lat,vert}range; model domain limits (or a subdomain?) + Aether longitudes; in the filter_input_#.nc some are not whole numbers.; 75.00001 + Doc error: web page says run_tests uses entries from test1thru, + but that has test 0, which is not an option in model_mod_check. + tests_to_run is not dimensioned '(0:'. + +&model_mod_check_nml + num_ens = 1 + single_file = .FALSE. + input_state_files = 'filter_input_0001.nc' + output_state_files = 'filter_output_0001.nc' + quantity_of_interest = 'QTY_DENSITY_ION_OP' + all_metadata_file = 'test6_metadata.txt' + x_ind = 1234 + loc_of_interest = 15.0, -2.5, 100000. + interp_test_dlon = 10.0 + interp_test_dlat = 5.0 + interp_test_dvert = 1500.0 + interp_test_lonrange = 0, 360 + interp_test_latrange = -87.5, 87.5 + interp_test_vertrange = 96952.5625, 436360.25 + interp_test_dx = -888888.0 + interp_test_dy = -888888.0 + interp_test_dz = -888888.0 + interp_test_xrange = -888888.0, -888888.0 + interp_test_yrange = -888888.0, -888888.0 + interp_test_zrange = -888888.0, -888888.0 + interp_test_vertcoord = 'VERTISHEIGHT' + test1thru = -1 + run_tests = 1,2,3,4,5,7 + verbose = .FALSE. + / + +&quad_interpolate_nml +/ diff --git a/models/aether_lat-lon/work/obs_seq.out.1 b/models/aether_lat-lon/work/obs_seq.out.1 new file mode 100644 index 0000000000..1cd14ccb15 --- /dev/null +++ b/models/aether_lat-lon/work/obs_seq.out.1 @@ -0,0 +1,31 @@ + obs_sequence +obs_kind_definitions + 1 + 33 AIRS_TEMPERATURE + num_copies: 1 num_qc: 1 + num_obs: 2 max_num_obs: 2 +observation +Data QC + first: 1 last: 2 + OBS 1 + 271.330627441406 + 0.000000000000000E+000 + -1 2 -1 +obdef +loc3d + 3.406717740263719 0.5806184282903090 100000.0000000000 3 +kind + 33 +84601 153130 + 1.07229244766182 + OBS 2 + 27450.2966235645 + 0.000000000000000E+000 + 1 -1 -1 +obdef +loc3d + 3.484538383406885 0.5925166389933947 120000.0000000000 3 +kind + 33 +84601 153130 + 1.03153675838621 \ No newline at end of file diff --git a/models/aether_lat-lon/work/quickbuild.sh b/models/aether_lat-lon/work/quickbuild.sh new file mode 100755 index 0000000000..dd1b063219 --- /dev/null +++ b/models/aether_lat-lon/work/quickbuild.sh @@ -0,0 +1,51 @@ +#!/usr/bin/env bash + +# DART software - Copyright 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 + +main() { + +export DART=$(git rev-parse --show-toplevel) +source "$DART"/build_templates/buildfunctions.sh + +MODEL=aether_lat-lon +LOCATION=threed_sphere + +programs=( +filter +model_mod_check +perfect_model_obs +) + +serial_programs=( +create_fixed_network_seq +create_obs_sequence +obs_diag +obs_seq_to_netcdf +) + +model_serial_programs=( +aether_to_dart +dart_to_aether) + +arguments "$@" + +# clean the directory +\rm -f -- *.o *.mod Makefile .cppdefs + +# build any NetCDF files from .cdl files +cdl_to_netcdf + +# build and run preprocess before making any other DART executables +buildpreprocess + +# build DART +buildit + +# clean up +\rm -f -- *.o *.mod + +} + +main "$@"