From 33bd3e17719bca5e74bf33729bfd0b3d4ef1ed48 Mon Sep 17 00:00:00 2001 From: Alexander Robinson Date: Sat, 18 Dec 2021 20:53:50 +0100 Subject: [PATCH] Disabled reading of grid information, since current version is outdated and not consistent with input files. Removed unused parameter ydyn.beta_gl_sep. Reduced initial timestep for SIA solver for 0.1yr equilibrium, before moving forward. Also updated ncio/nml to latest versions (though not really necessary in the end). --- libs/ncio.f90 | 180 ++++++++++++++++++++++++++++++++++------ libs/nml.f90 | 29 +++++-- par/yelmo_initmip.nml | 1 + src/yelmo_defs.f90 | 1 - src/yelmo_dynamics.f90 | 1 - src/yelmo_grid.f90 | 6 +- tests/yelmo_initmip.f90 | 3 + 7 files changed, 185 insertions(+), 36 deletions(-) diff --git a/libs/ncio.f90 b/libs/ncio.f90 index 82083a7a..1e43e1c7 100644 --- a/libs/ncio.f90 +++ b/libs/ncio.f90 @@ -22,6 +22,7 @@ module ncio + use ieee_arithmetic use netcdf implicit none @@ -36,6 +37,8 @@ module ncio character(len=NC_STRLEN), parameter :: NC_CHARDIM = "strlen" + character(len=3), parameter :: GRID_MAPPING_NAME_DEFAULT = "crs" + type ncvar character (len=NC_STRLEN) :: name, long_name, standard_name, units character (len=NC_STRLEN) :: axis, calendar, grid_mapping @@ -381,6 +384,8 @@ subroutine nc4_read_internal_numeric(filename,name,dat,size_in,start,count,xtype integer, optional :: iostat integer :: i, status + double precision, parameter :: missing_value_default = -9999.0 + ! Open the file. call nc_check_open(filename, ncid, nf90_nowrite, nc_id) @@ -440,6 +445,15 @@ subroutine nc4_read_internal_numeric(filename,name,dat,size_in,start,count,xtype if (.not. present(ncid)) call nc_check( nf90_close(nc_id) ) if (v%missing_set) then + + ! === SPECIAL CASE: missing_value == NaN ==== + + ! Replace NaNs with internal missing value to avoid crashes + where ( ieee_is_nan(dat) ) dat = missing_value_default + v%missing_value = missing_value_default + + ! =========================================== + where( dabs(dat-v%missing_value) .gt. NC_TOL ) dat = dat*v%scale_factor + v%add_offset ! Fill with user-desired missing value @@ -543,12 +557,12 @@ subroutine nc_v_init(v,name,xtype,ndims_in,long_name,standard_name, & v%missing_set = .TRUE. v%missing_value = -9999d0 v%FillValue = v%missing_value - v%FillValue_set = .FALSE. + v%FillValue_set = .TRUE. v%xtype = "NF90_DOUBLE" v%coord = .FALSE. - v%grid_mapping = "" + v%grid_mapping = GRID_MAPPING_NAME_DEFAULT ! If args are present, reassign these variables if ( present(long_name) ) v%long_name = trim(long_name) @@ -880,14 +894,15 @@ subroutine nc_put_att(ncid, v) if (v%FillValue_set) then select case(trim(v%xtype)) case("NF90_INT") - call nc_check( nf90_put_att(ncid, v%varid, "FillValue", int(v%FillValue)) ) + call nc_check( nf90_put_att(ncid, v%varid, "_FillValue", int(v%FillValue)) ) case("NF90_FLOAT") - call nc_check( nf90_put_att(ncid, v%varid, "FillValue", real(v%FillValue)) ) + call nc_check( nf90_put_att(ncid, v%varid, "_FillValue", real(v%FillValue)) ) case("NF90_DOUBLE") - call nc_check( nf90_put_att(ncid, v%varid, "FillValue", v%FillValue) ) + call nc_check( nf90_put_att(ncid, v%varid, "_FillValue", v%FillValue) ) end select end if + end if ! Add additional variable attributes @@ -1014,10 +1029,17 @@ subroutine nc_get_att(ncid, v, readmeta) call nc_get_att_double(ncid,v%varid,"missing_value",v%missing_value,stat) if (stat .eq. noerr) v%missing_set = .TRUE. - stat = nc_check_att( nf90_get_att(ncid, v%varid, "FillValue", tmpi) ) + stat = nc_check_att( nf90_get_att(ncid, v%varid, "_FillValue", tmpi) ) if (stat .eq. noerr) then v%FillValue = dble(tmpi) v%FillValue_set = .TRUE. + + ! ajr, 2020-08-20 + ! Overwrite missing value with Fillvalue, since missing_value + ! attribute has been deprecated. + v%missing_value = v%FillValue + v%missing_set = .TRUE. + end if case DEFAULT @@ -1034,7 +1056,7 @@ subroutine nc_get_att(ncid, v, readmeta) call nc_get_att_double(ncid,v%varid,"missing_value",v%missing_value,stat) if (stat .eq. noerr) v%missing_set = .TRUE. - stat = nc_check_att( nf90_get_att(ncid, v%varid, "FillValue", tmp) ) + stat = nc_check_att( nf90_get_att(ncid, v%varid, "_FillValue", tmp) ) if (stat .eq. noerr) then v%FillValue = tmp v%FillValue_set = .TRUE. @@ -1418,61 +1440,166 @@ subroutine nc_write_attr_global_dp(filename, name, value) call nc_check( nf90_close(ncid) ) end subroutine - subroutine nc_write_map(filename,name,lambda,phi,alpha,x_e,y_n, ncid) + subroutine nc_write_map(filename,grid_mapping_name,lambda,phi,alpha,x_e,y_n, & + is_sphere,semi_major_axis,inverse_flattening,ncid) ! CF map conventions can be found here: ! http://cfconventions.org/Data/cf-conventions/cf-conventions-1.6/build/cf-conventions.html#appendix-grid-mappings implicit none - character(len=*) :: filename, name - + character(len=*), intent(IN) :: filename, grid_mapping_name + double precision, intent(IN), optional :: lambda, phi + double precision, intent(IN), optional :: alpha, x_e, y_n + logical, intent(IN), optional :: is_sphere + double precision, intent(IN), optional :: semi_major_axis + double precision, intent(IN), optional :: inverse_flattening + integer, intent(in), optional :: ncid + + ! Local variables integer :: nc_id, varid, stat - integer, intent(in), optional :: ncid - double precision, optional :: lambda, phi, alpha, x_e, y_n + character(len=56) :: crs_name + double precision :: phi_proj_orig integer, parameter :: noerr = NF90_NOERR + ! Define general name of integer variable to + ! hold coordinate reference system (CRS) information + ! Always, use crs by default to be consistent with others. + crs_name = "crs" + ! Open the file, set for redefinition call nc_check_open(filename, ncid, nf90_write, nc_id) call nc_check( nf90_redef(nc_id) ) ! Check if grid mapping has been defined in this file ! (if not, define it according to input arguments) - stat = nf90_inq_varid(nc_id, trim(name), varid) + stat = nf90_inq_varid(nc_id, trim(crs_name), varid) if ( stat .ne. noerr ) then ! Define the mapping variable as an integer with no dimensions, ! and include the grid mapping name - call nc_check( nf90_def_var(nc_id, trim(name), NF90_INT, varid) ) - call nc_check( nf90_put_att(nc_id,varid, "grid_mapping_name", trim(name)) ) - + + call nc_check( nf90_def_var(nc_id, trim(crs_name), NF90_INT, varid) ) + ! Add grid attributes depending on grid_mapping type - select case(trim(name)) + select case(trim(grid_mapping_name)) case("stereographic") + ! Including 'oblique_stereographic' + + if ( (.not. present(lambda)) .or. & + (.not. present(phi)) .or. & + (.not. present(alpha)) .or. & + (.not. present(x_e)) .or. & + (.not. present(y_n)) ) then + + write(*,"(a,a)") "ncio:: nc_write_map:: Error: ", & + "All grid_mapping arguments must be provided ", & + "(lambda, phi, alpha, x_e, y_n)." + stop + + end if + + ! Add grid mapping attributes + call nc_check( nf90_put_att(nc_id,varid, "grid_mapping_name", trim(grid_mapping_name)) ) call nc_check( nf90_put_att(nc_id,varid, "longitude_of_projection_origin", lambda) ) call nc_check( nf90_put_att(nc_id,varid, "latitude_of_projection_origin", phi) ) - if (present(alpha)) & call nc_check( nf90_put_att(nc_id,varid, "angle_of_oblique_tangent", alpha) ) call nc_check( nf90_put_att(nc_id,varid, "scale_factor_at_projection_origin", 1.d0) ) call nc_check( nf90_put_att(nc_id,varid, "false_easting", x_e) ) call nc_check( nf90_put_att(nc_id,varid, "false_northing", y_n) ) case("polar_stereographic") + + if ( (.not. present(lambda)) .or. & + (.not. present(phi)) .or. & + (.not. present(x_e)) .or. & + (.not. present(y_n)) ) then + + write(*,"(a,a)") "ncio:: nc_write_map:: Error: ", & + "All grid_mapping arguments must be provided ", & + "(lambda, phi, x_e, y_n)." + stop + + end if + + ! Determine latitude_of_projection_origin, since it must + ! be either -90 or +90 for a polar_stereographic projection: + if (phi .gt. 0.0d0) then + phi_proj_orig = 90.0d0 + else + phi_proj_orig = -90.0d0 + end if + + ! Add grid mapping attributes + call nc_check( nf90_put_att(nc_id,varid, "grid_mapping_name", trim(grid_mapping_name)) ) call nc_check( nf90_put_att(nc_id,varid, "straight_vertical_longitude_from_pole", lambda) ) - call nc_check( nf90_put_att(nc_id,varid, "latitude_of_projection_origin", phi) ) - if (present(alpha)) & - call nc_check( nf90_put_att(nc_id,varid, "angle_of_oblique_tangent", alpha) ) - call nc_check( nf90_put_att(nc_id,varid, "scale_factor_at_projection_origin", 1.d0) ) + call nc_check( nf90_put_att(nc_id,varid, "latitude_of_projection_origin", phi_proj_orig) ) + call nc_check( nf90_put_att(nc_id,varid, "standard_parallel", phi) ) + + call nc_check( nf90_put_att(nc_id,varid, "false_easting", x_e) ) call nc_check( nf90_put_att(nc_id,varid, "false_northing", y_n) ) + case("lambert_conformal_conic") + + call nc_check( nf90_put_att(nc_id,varid, "longitude_of_central_meridian", lambda) ) + call nc_check( nf90_put_att(nc_id,varid, "latitude_of_projection_origin", phi) ) + if (present(alpha)) & + call nc_check( nf90_put_att(nc_id,varid, "standard_parallel", alpha) ) + + case("latitude_longitude","latlon","gaussian") + + ! Pass - no grid mapping needed for a latitude_longitude grid + + ! Add grid mapping attributes + !call nc_check( nf90_put_att(nc_id,varid, "grid_mapping_name", "latitude_longitude") ) + ! == No additional parameters needed == + case DEFAULT ! Do nothing - end select + end select + + + select case(trim(grid_mapping_name)) + + case("latitude_longitude","latlon","gaussian") - write(*,"(a,a,a)") "ncio:: nc_write_map:: ",trim(filename)//" : ",trim(name) + ! Pass - do nothing for latitude_longitude grids + + case DEFAULT + ! For other projections, add planet information if available + + ! Add planet information if desired + if (present(is_sphere) .and. & + present(semi_major_axis) .and. & + present(inverse_flattening)) then + + if (is_sphere) then + call nc_check( nf90_put_att(nc_id,varid, "semi_major_axis", semi_major_axis) ) + call nc_check( nf90_put_att(nc_id,varid, "inverse_flattening", 0.0d0) ) + else + call nc_check( nf90_put_att(nc_id,varid, "semi_major_axis", semi_major_axis) ) + call nc_check( nf90_put_att(nc_id,varid, "inverse_flattening", inverse_flattening) ) + + end if + + else if (present(is_sphere) .or. & + present(semi_major_axis) .or. & + present(inverse_flattening)) then + + write(*,*) "ncio:: nc_write_map:: Error: to write planet information, & + &all three planet parameters must be provided: & + &is_sphere, semi_major_axis, inverse_flattening. & + &Try again." + + stop + end if + + end select + + write(*,"(a,a,a)") "ncio:: nc_write_map:: ",trim(filename)//" : ",trim(grid_mapping_name) end if @@ -1495,7 +1622,6 @@ end subroutine nc_write_map !! @param units NetCDF attribute of the units of the variable (optional) !! @param axis NetCDF attribute of the standard axis of the variable (optional) !! @param calendar NetCDF attribute of the calendar type to be used for time dimensions (optional) - subroutine nc_write_dim_int_pt(filename,name,x,dx,nx, & long_name,standard_name,units,axis,calendar,unlimited, ncid) @@ -1723,6 +1849,10 @@ subroutine nc_write_dim_internal(filename,name,xtype,x, & call nc_v_init(v,name=trim(name),xtype=xtype,coord=.TRUE.) + ! Ensure grid_mapping is set to blank since a dimension cannot + ! have a grid mapping. + v%grid_mapping = "" + !! Now fill in values of arguments that are present if ( present(long_name) ) v%long_name = trim(long_name) if ( present(standard_name) ) v%standard_name = trim(standard_name) diff --git a/libs/nml.f90 b/libs/nml.f90 index b1796b0e..d2106790 100644 --- a/libs/nml.f90 +++ b/libs/nml.f90 @@ -1,10 +1,14 @@ module nml + use, intrinsic :: iso_fortran_env, only : input_unit, output_unit, error_unit implicit none - logical :: VERBOSE = .FALSE. !! should freshly read namelist be printed to screen? + logical :: VERBOSE = .FALSE. !! should freshly read namelist be printed to screen? + logical :: ERROR_NO_PARAM = .TRUE. !! Should error be thrown if parameter isn't found? + + integer, parameter :: io_unit_err = error_unit interface nml_read module procedure nml_read_string, nml_read_double, nml_read_float @@ -46,7 +50,7 @@ subroutine nml_set_verbose(verb) return - end subroutine nml_set_verbose + end subroutine nml_set_verbose subroutine nml_replace(s,text,rep,outs) ! Adapted from FUNCTION Replace_Text: @@ -90,10 +94,12 @@ subroutine nml_read_internal(filename,group,name,value,comment) character(len=*), intent(IN) :: filename, group, name character(len=*), intent(INOUT), optional :: comment + ! Local variables integer :: io, file_unit integer :: iostat, l, ltype character(len=1000) :: line, name1, value1, comment1 logical :: ingroup + logical :: found_param ! Open the nml filename to be read, or get file units if already open inquire(file=trim(filename),NUMBER=file_unit) @@ -111,7 +117,8 @@ subroutine nml_read_internal(filename,group,name,value,comment) end if - ingroup = .FALSE. + ingroup = .FALSE. + found_param = .FALSE. do l = 1, 5000 read(io,"(a1000)",iostat=iostat) line @@ -121,8 +128,9 @@ subroutine nml_read_internal(filename,group,name,value,comment) ! Check if the parameter has been found if (ingroup .and. ltype == 3) then if (trim(name1) == trim(name)) then - value = trim(value1) - comment = trim(comment1) + value = trim(value1) + comment = trim(comment1) + found_param = .TRUE. exit end if end if @@ -141,6 +149,15 @@ subroutine nml_read_internal(filename,group,name,value,comment) close(io) endif + if (ERROR_NO_PARAM .and. (.not. found_param)) then + write(io_unit_err,*) "" + write(io_unit_err,*) "nml:: Error: parameter not found." + write(io_unit_err,*) "Filename: ", trim(filename) + write(io_unit_err,*) "Group: ", trim(group) + write(io_unit_err,*) "Parameter: ", trim(name) + stop "Program stopped." + end if + return end subroutine nml_read_internal @@ -170,7 +187,7 @@ subroutine nml_read_string(filename,group,name,value,comment,init) return - end subroutine nml_read_string + end subroutine nml_read_string subroutine nml_read_double(filename,group,name,value,comment,init) diff --git a/par/yelmo_initmip.nml b/par/yelmo_initmip.nml index ffb28192..9e398d95 100644 --- a/par/yelmo_initmip.nml +++ b/par/yelmo_initmip.nml @@ -94,6 +94,7 @@ beta_gl_stag = 1 ! 0: simple staggering, 1: Upstream beta at gl, 2: downstream, 3: f_grnd_ac scaling beta_gl_f = 1.0 ! [-] Scaling of beta at the grounding line (for beta_gl_scale=0) taud_gl_method = 0 ! 0: binary, no subgrid, 1: Two-sided gradient + calc_diffusivity = False H_grnd_lim = 500.0 ! [m] For beta_gl_scale=1, reduce beta linearly between H_grnd=0 and H_grnd_lim H_sed_sat = 250.0 ! [m] Sediment thickness at which sediment effect is saturated cb_method = 1 ! -1: set externally; 1: calculate cb online diff --git a/src/yelmo_defs.f90 b/src/yelmo_defs.f90 index 530746eb..e107bfa6 100644 --- a/src/yelmo_defs.f90 +++ b/src/yelmo_defs.f90 @@ -179,7 +179,6 @@ module yelmo_defs real(prec) :: beta_q ! Friction law exponent real(prec) :: beta_u0 ! [m/a] Friction law velocity threshold integer :: beta_gl_scale ! Beta grounding-line scaling method (beta => 0 at gl?) - integer :: beta_gl_sep ! Beta grounding-line sub-element (subgrid) parameterization integer :: beta_gl_stag ! Beta grounding-line staggering method real(prec) :: beta_gl_f ! Fraction of beta at gl integer :: taud_gl_method ! Driving stress grounding line treatment diff --git a/src/yelmo_dynamics.f90 b/src/yelmo_dynamics.f90 index 138e11d1..e5f35f89 100644 --- a/src/yelmo_dynamics.f90 +++ b/src/yelmo_dynamics.f90 @@ -961,7 +961,6 @@ subroutine ydyn_par_load(par,filename,zeta_aa,zeta_ac,nx,ny,dx,init) call nml_read(filename,"ydyn","beta_q", par%beta_q, init=init_pars) call nml_read(filename,"ydyn","beta_u0", par%beta_u0, init=init_pars) call nml_read(filename,"ydyn","beta_gl_scale", par%beta_gl_scale, init=init_pars) - call nml_read(filename,"ydyn","beta_gl_sep", par%beta_gl_sep, init=init_pars) call nml_read(filename,"ydyn","beta_gl_stag", par%beta_gl_stag, init=init_pars) call nml_read(filename,"ydyn","beta_gl_f", par%beta_gl_f, init=init_pars) call nml_read(filename,"ydyn","taud_gl_method", par%taud_gl_method, init=init_pars) diff --git a/src/yelmo_grid.f90 b/src/yelmo_grid.f90 index 6dfc4723..763032aa 100644 --- a/src/yelmo_grid.f90 +++ b/src/yelmo_grid.f90 @@ -42,7 +42,7 @@ subroutine yelmo_init_grid_fromfile(grd,filename,grid_name) ! Define the ygrid name grd%name = trim(grid_name) - + ! Determine grid axis sizes grd%nx = nc_size(filename,"xc") grd%ny = nc_size(filename,"yc") @@ -71,7 +71,7 @@ subroutine yelmo_init_grid_fromfile(grd,filename,grid_name) call nc_read(filename,"lon2D",grd%lon) call nc_read(filename,"lat2D",grd%lat) call nc_read(filename,"area", grd%area) - + ! Determine axis units call nc_read_attr(filename,"xc","units",units) @@ -118,7 +118,7 @@ subroutine yelmo_init_grid_fromfile(grd,filename,grid_name) if (grd%is_projection) then ! Load additional projection information if available - if ( nc_exists_attr(filename,"lon2D","grid_mapping") ) then + if (.FALSE. .and. nc_exists_attr(filename,"lon2D","grid_mapping") ) then ! Read grid map name (projection name, eg "polar_stereographic") call nc_read_attr(filename,"lon2D","grid_mapping",grd%mtype) diff --git a/tests/yelmo_initmip.f90 b/tests/yelmo_initmip.f90 index bd7b8e21..47945155 100644 --- a/tests/yelmo_initmip.f90 +++ b/tests/yelmo_initmip.f90 @@ -232,6 +232,9 @@ program yelmo_test ! Note: run this with SIA only dynamics for now yelmo1%dyn%par%solver = "sia" + call yelmo_update_equil(yelmo1,time,time_tot=0.1_prec,topo_fixed=.FALSE., & + dt=0.01_prec,ssa_vel_max=5000.0_prec) + call yelmo_update_equil(yelmo1,time,time_tot=10.0_prec,topo_fixed=.FALSE., & dt=0.2_prec,ssa_vel_max=5000.0_prec)