The basic idea is to derive the Cretaceous setups from the GFDL-supplied EBM example case. For each of the cretaceous periods, a topography file is supplied, which contains the complete land and sea map, at resolution 128x64 cells. @Article{cp-3-647-2007, AUTHOR = {Sewall, J. O. and van de Wal, R. S. W. and van der Zwan, K. and van Oosterhout, C. and Dijkstra, H. A. and Scotese, C. R.}, TITLE = {Climate model boundary conditions for four Cretaceous time slices}, JOURNAL = {Climate of the Past}, VOLUME = {3}, YEAR = {2007}, NUMBER = {4}, PAGES = {647--657}, URL = {http://www.clim-past.net/3/647/2007/}, DOI = {10.5194/cp-3-647-2007} } Dr. Jacob O. Sewall Assistant Professor of Geology Department of Physical Sciences Kutztown University of Pennsylvania Kutztown, PA 19530 http://faculty.kutztown.edu/sewall/ sewall@kutztown.edu 484-646-5864 Show the land elevation map in ferret: use maast_70Ma_vtop.nc shade/levels=h topo shade/over if oro eq 0 then oro We need: - ocean topography, with land areas cut off above zero make_coupler_mosaic requires: The topography data is positive down. - land topography, with oceans cut off below zero - grid_spec with land, ocean and atmosphere all set to 128x64 resolution. We can try to use tripolar grid for the ocean, however, the land masses in the archaean data sets are distributed differently from today, thus we might actually get even three instead of on open-sea singularity. 0. lets start with 070Ma cd 070Ma mkdir INPUT cd INPUT 1a. Generate ocean horizontal grid: (nlon and nlat are on supergrid size, i.e. doubled values) ('c_cell' should be used for the grid used in MOM.) In the GFDL-delivered examples, the poles are located at 280W (Siberia) and 100W (Canada). But for the Cretacious world maps, this must be reconsidered. They must be moved a little to the east. ../../../../src/tools/make_hgrid/make_hgrid \ --grid_type tripolar_grid \ --nxbnd 2 --xbnd -270\,90 \ --nybnd 2 --ybnd -90\,90 \ --nlon 256 --nlat 128 \ --grid_name tripolar_128x64 \ --lat_join 65 \ --center c_cell To visualize generated output files with tripolar grid in ferret, see: http://www.gfdl.noaa.gov/~atw/ferret/tripolar/ Instead of geolon_vertex and geolat_vertex use geolon_c and geolat_c 1b. Generate ocean mosaic: ../../../../src/tools/make_solo_mosaic/make_solo_mosaic \ --num_tiles 1 --dir ./ --mosaic_name ocean_mosaic \ --tile_file tripolar_128x64.nc --periodx 360 1c. Generate ocean vertical grid ../../../../src/tools/make_vgrid/make_vgrid \ --center c_cell \ --grid_name ocean_vgrid \ --nbnds 6 --bnds 0\,70\,80\,200\,500\,5900 --nz 14\,2\,10\,10\,20 (Oder ein ocean_vgrid.nc kopieren aus dem EBM-Beispiel, oder aus dem om3_core1 Beispiel) #cp -p ../../ebm_mosaic/ocean_vgrid.nc . 1d. Generate ocean topography ocean_topog_init: Regions with 0 < kmt < 2 have been found. Recommend kmt > 2 for wet ocean regions. Thus, use parameter --kmt_min 3 (after fixing the Bug in MOM 5.0.2 topog.c) ncatted -a missing_value\,topo\,c\,f\,-99999.f ../maast_70Ma_vtop.nc #../../../../src/tools/make_topog/make_topog \ # --verbose \ # --mosaic ocean_mosaic.nc \ # --topog_type realistic \ # --vgrid_file ocean_vgrid.nc \ # --topog_file ../maast_70Ma_vtop.nc --topog_field topo --scale_factor -1 \ # --bottom_depth 5900 \ # --min_depth 10 \ # --kmt_min 3 \ # --output topog.nc Sigh deeply. Command line parameter --min_depth is not used at all for the case --topog_type realistic. Parameter --kmt_min is useless, because the MOM5 model computes that on its own during ocean init. Thus adjust ocean depth via ncap2, and also make sure that all land has positive height, so that the land mask gets calculated correctly. ncap2 -O -v -c -s ' oro=oro; where(oro>=0.5 && topo < 10) topo=1; where(oro< 0.5 && topo>-50) topo=-50; ' ../maast_70Ma_vtop.nc ../maast_70Ma_vtop-min_depth50.nc ../../../../src/tools/make_topog/make_topog \ --verbose \ --mosaic ocean_mosaic.nc \ --topog_type realistic \ --vgrid_file ocean_vgrid.nc \ --topog_file ../maast_70Ma_vtop-min_depth50.nc --topog_field topo --scale_factor -1 \ --bottom_depth 5900 \ --output topog.nc Ferret: use maast_70Ma_vtop-min_depth50.nc let ht=if oro eq 0 then topo shade if ht gt -50 then ht Question: is this now on tripolar grid already? Or do we need to use regrid onto the tripolar grid first? Answer: it seems that Yes Question: can we convert the archaean topography input file to a grid with higher resolution? just using make_topog? or do we need to apply fregrid first? Answer: the code in topog.c looks like it does conservative interpolation from source grid onto above-generated model ocean grid. 2a. Generate atmos horizontal grid. First try was with --grid_type regular_lonlat_grid, but the coupler does not like that. Next try: (nlon and nlat are for supergrid resolution) ../../../../src/tools/make_hgrid/make_hgrid \ --grid_type spectral_grid \ --nxbnd 2 --nybnd 2 --xbnd 0\,360 --ybnd -90\,90 \ --nlon 256 --nlat 128 \ --verbose \ --grid_name spectral_128x64 2b. Generate atmos mosaic. ../../../../src/tools/make_solo_mosaic/make_solo_mosaic \ --num_tiles 1 --dir ./ --mosaic_name atmos_mosaic \ --tile_file spectral_128x64.nc --periodx 360 3a. Generate land horizontal grid. (nlon and nlat are for supergrid resolution) ../../../../src/tools/make_hgrid/make_hgrid \ --grid_type regular_lonlat_grid \ --nxbnd 2 --nybnd 2 --xbnd 0\,360 --ybnd -90\,90 \ --nlon 256 --nlat 128 \ --verbose \ --grid_name regular_128x64 3b. Generate land mosaic. #../../../../src/tools/make_solo_mosaic/make_solo_mosaic \ # --num_tiles 1 --dir ./ --mosaic_name land_mosaic \ # --tile_file regular_128x64.nc --periodx 360 ../../../../src/tools/make_solo_mosaic/make_solo_mosaic \ --num_tiles 1 --dir ./ --mosaic_name land_mosaic \ --tile_file regular_128x64.nc --periody 1 4. Generate grid_spec mosaic file This also generates land_mask.nc and ocean_mask.nc ../../../../src/tools/make_coupler_mosaic/make_coupler_mosaic \ --verbose \ --atmos_mosaic atmos_mosaic.nc \ --land_mosaic land_mosaic.nc \ --ocean_mosaic ocean_mosaic.nc \ --ocean_topog topog.nc \ --check \ --mosaic grid_spec NOTE: tiling error is -0.000000% ***** Congratulation! You have successfully run make_coupler_mosaic 5. Generate land topography netcdf navy_topography.data { dimensions: x = 2160 ; y = 1080 ; x1 = 1 ; x2 = 2161 ; x3 = 1081 ; variables: float zdat(y, x) ; float xdat(x2) ; float ydat(x3) ; float ipts(x1) ; float jpts(x1) ; data: ipts = 2160 ; jpts = 1080 ; xdat = 0, 0.002908882, 0.005817764, 0.008726647, 0.01163553, 0.01454441, ... 6.268529, 6.271438, 6.274347, 6.277256, 6.280164, 6.283185 ; ydat = -1.570796, -1.567888, -1.564979, -1.56207, -1.559161, -1.556252, ... 1.564959, 1.567868, 1.570796 ; } We can do that with NCO. However, on the iplex login nodes, ncap2 segfaults with this. On the visualization server, it works fine: ncap2 -O -v -c -s 'defdim("x",128);defdim("y",64);defdim("x1",1);defdim("x2",$x.size+1);defdim("x3",$y.size+1); ipts[$x1]=$x.size; jpts[$x1]=$y.size; xdat[$x2]=0.0; ydat[$x3]=0.0; zdat[$y,$x]=0.0; where (topo>0) zdat=topo; xdat(0:$x2.size-2)=lon(:); xdat($x2.size-1)=360; ydat(0)=-90; ydat($x3.size-1)=90; ydat(1:$y.size-1) = (lat(0:$y.size-2) + lat(1:$y.size-1)) / 2; *pi=4*atan(1.0); xdat=xdat*pi/180; ydat=ydat*pi/180;' \ ../maast_70Ma_vtop-min_depth50.nc \ maast_70ma_topography.nc #for(*i=1;i<$y.size;i++){ydat(i)=(lat(i)+lat(i-1))/2; }' \ In the file input.nml, this file must be specified in the section topography_nml: &topography_nml topog_file = 'maast_70ma_topography.nc' / 6. Generate cover type map alias vegetation map * Cretaceous: The 12 plant functional types, number code, and associated characteristics are: 0 = ocean 1 = land ice 3 = High altitude/latitude evergreen conifer closed canopy forest. [arctic,antarctic] 6 = High altitude/latitude mixed forest with equal percentage broad v. needle leaf and evergreen v. deciduous. [arctic,antarctic] 7 = Low altitude/latitude evergreen conifer closed canopy forest. [not in 070Ma] 10 = Closed canopy, broad leaved, moist evergreen forest. [rainforest?] 11 = Closed canopy, broad leaved, dry deciduous forest. [?] 12 = Savanna: dry, low understory with sparse, broad leaved overstory. 15 = High altitude/latitude moist, open canopy evergreen forest with shrub understory. 16 = Low altitude/latitude moist, open canopy, mixed evergreen/deciduous forest with shrub understory. 20 = Wet or cool shrubland (evergreen). [Kalte Steppe] 21 = Dry or warm shrubland (deciduous). [Warme Steppe] * land_lad: ! veg_to_use choice of method for defining vegetation [character] ! 'variable' - use input map to index vegetation parameter vectors ! 'constant' - use veg_index_constant to index veg parameter vectors ! 'cons_tun' - use tuned global constant vegetation ! 'cons_ssl' - use Manabe Climate Model-like vegetation (forces use_climap to be true) ! !---- Cover (vegetation) type. There are 12 possible cases: veg_to_use can ! take 4 values, and glaciers can be (1) absent, (2) based on input cover ! field, or (3) based on climap albedo field. (use of climap albedo forces ! use of climap to locate glaciers, if glaciers are used.) the code ! ensures consistency between cover_type field and glacier mask, with ! the latter taking precedence: where glacier based on input cover field ! must be removed in deference to climap or if no glaciers are used, ! such points are assigned tundra (or global constant) type. ! ! the following namelist vectors contain properties indexed by cover (veg) type: ! 0 ocean ! 1 (BE) broadleaf evergreen trees [rainforest] ! 2 (BD) broadleaf deciduous trees [middle europe / E USA] ! 3 (BN) broadleaf/needleleaf trees [middle europe] ! 4 (NE) needleleaf evergreen trees [N Russia / S Canada] ! 5 (ND) needleleaf deciduous trees [taiga] ! 6 (G) grassland ! 7 (D) desert ! 8 (T) tundra ! 9 (A) agriculture ! 10 (I) ice ! 11 (L) lake ! 12 (TV) tuned global vegetation ! 13 (SV) Manabe Climate Model-like vegetation ! 14 (SI) Manabe Climate Model-like ice/glacier Cret lad GF SP 0 -> 0 0 ocean 1 -> 10 10 ice 3 -> 4 4 6 -> 3 3 (apparently it was much warmer in that time) 7 -> 4 4 (not in 070Ma) 10 -> 1 1 rainforest 11 -> 2 5 ? Cretaceous grows around the equator, land_lad taiga around polar circle 12 -> 6 6 15 -> 1o.3 4 ? 16 -> 3 6 20 -> 6 6 21 -> 6 6 Note that in maast_70Ma_vtop.nc longitude goes 0..360, but in cover_type_field.nc it goes -180..180. Thus we also rotate the longitudes, and also shift coordinates by 0.5 grid cells to get consistent behaviour with the GFDL-supplied cover_type_field.nc file. ncap2 -O -v -c -s 'defdim("time",1);defdim("blon",$lon.size+1);defdim("blat",$lat.size+1); *c[$time,$lat,$lon]=0L; where(pft==1) c=10; where(pft==3) c=4; where(pft==6) c=3; where(pft==7) c=4; where(pft==10) c=1; where(pft==11) c=5; where(pft==12) c=6; where(pft==15) c=4; where(pft==16) c=6; where(pft==20) c=6; where(pft==21) c=6; cover[$time,$lat,$lon]=0L; cover@long_name="cover type"; cover@units="unitless"; cover@missing_value=0; cover(:,:,0:$lon.size/2-1) = c(:,:,$lon.size/2:$lon.size-1); cover(:,:,$lon.size/2:$lon.size-1) = c(:,:,0:$lon.size/2-1); *l = lon; *dlon_2=(lon(1)-lon(0))/2; lon(0:$lon.size/2-1) = l($lon.size/2:$lon.size-1)-360+dlon_2; lon($lon.size/2:$lon.size-1) = l(0:$lon.size/2-1)+dlon_2; blon[$blon] = 0.0; blon(0) = -180; blon($blon.size-1) = 180; blon(1:$blon.size-2) = (lon(0:$lon.size-2) + lon(1:$lon.size-1)) / 2; blat[$blat] = 0.0; blat(0)=-90; blat($lat.size)=90; blat(1:$lat.size-1)=(lat(0:$lat.size-2) + lat(1:$lat.size-1)) / 2; lon@edges="blon"; lat@edges="blat"; blon@long_name="longitude edges"; blon@units="degrees_east"; blat@long_name="latitude edges"; blat@units="degrees_north"; ' \ ../maast_70Ma_vtop.nc \ cover_type_field.nc Das File von GFDL enthaelt uebrigens 442 Zeitscheiben, in years since 1860, also fuer den Zeitraum 1860-2301. Keine Ahnung wie sie die erzeugt haben. Bevor man mit Climber-4 was mit present-day produziert, sollte man das wohl nochmal versuchen in deren Papers nachzulesen. Naechstes TODO: mal pruefen, wie das cover_type_field gegenueber de beim Gitter-Bauen erzeugten Land-See-Maske aussieht. ------------------------ * river routing : ===== land_lad ==== The file river_destination_field contains two times 96x80 values (2*96*80=2*7680). Apparently, there is one pair of values for each surface cell. The first value denotes column, the second value denotes row index of the destination cell. Cells which point to themselfes are ocean cells. The input data is regridded dynamically to the actually used model grid resolution. After regridding, for each land cell it is made sure that its runoff goes into an ocean cell. Apparently, runoff from each cell goes directly into its corresponding ocean cell (???) The regridded river routing map is sent to static diag fields rivers|dest|destination points|none|2|F| 0.000000000000000E+000|||lon,lat rivers|basin|river basins|none|2|F| 0.000000000000000E+000|||lon,lat subroutine init_routing ( glonb, glatb, gfrac, is,ie, js,je, i_dest,j_dest ) call mpp_open (unit, 'INPUT/river_destination_field', action = MPP_RDONLY) read (unit,*) n_lon_in, n_lat_in, wb_degrees read (unit,*) input_format do j=1,n_lat_in read(unit,input_format) (input_1(i,j),i=1,n_lon_in) end do do j=1,n_lat_in read(unit,input_format) (input_2(i,j),i=1,n_lon_in) end do route_input = input_1 + n_lon_in*(input_2-1) deallocate ( input_1, input_2 ) subroutine calc_discharge ( rivers, runoff, discharge ) ! obtain the global field call mpp_sum ( g_runoff, rivers%gnlon*rivers%gnlat ) ! do i = lbound(g_runoff,1),ubound(g_runoff,1) do j = lbound(g_runoff,2),ubound(g_runoff,2) if ( rivers%is<=rivers%i_dest(i,j).and.rivers%i_dest(i,j)<=rivers%ie.and.& rivers%js<=rivers%j_dest(i,j).and.rivers%j_dest(i,j)<=rivers%je ) & discharge(rivers%i_dest(i,j),rivers%j_dest(i,j)) = & discharge(rivers%i_dest(i,j),rivers%j_dest(i,j)) + g_runoff(i,j)*rivers%larea(i,j) enddo enddo define axis/x=179.75w:179.75e:0.5 x05deg define axis/y=89.75s:89.75n:0.5 y05deg set data/ez/variables=ddm/grid=grid05deg/skip=6/columns=720 "/data/biosx/LPJ/input_new/DDM30.asc" ! this data set appears upside-down fill/vlimits=90:-90 ddm ln -s river_destination_field river_destination_field2 define axis/x=-178:178:`360/96` x96 define axis/y=-89:89:`180/80` y80 define grid/x=x96/y=y80 g96x80 set data/ez/variables=rdf1/grid=g96x80/skip=2/columns=96 "river_destination_field" set data/ez/variables=rdf2/grid=g96x80/skip=402/columns=96 "river_destination_field2" fill rdf1[d=1]+96*(rdf2[d=2]-1) read (unit,*) n_lon_in, n_lat_in, wb_degrees read (unit,*) input_format do j=1,n_lat_in read(unit,input_format) (input_1(i,j),i=1,n_lon_in) end do do j=1,n_lat_in read(unit,input_format) (input_2(i,j),i=1,n_lon_in) end do route_input = input_1 + n_lon_in*(input_2-1) 85 85 85 85 85 85 85 85 85 85 84 84 84 84 84 84 84 84 84 84 84 84 84 84 84 84 84 84 84 84 84 84 84 84 84 84 84 84 84 84 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 84 84 84 84 84 84 84 84 84 84 84 84 84 83 83 84 84 84 84 84 84 84 84 84 84 84 84 84 84 84 84 84 84 85 85 85 85 85 85 85 85 85 85 85 85 85 85 85 85 86 86 86 86 86 86 86 86 86 86 49 49 49 49 49 49 49 49 49 49 49 49 49 49 49 49 49 49 49 49 49 49 49 49 49 50 50 50 51 51 51 51 51 51 51 51 51 51 51 51 51 51 51 51 51 51 83 83 83 83 83 83 83 83 83 83 83 84 84 84 84 84 84 85 85 85 85 85 85 85 85 85 85 ===== land_lad2 ==== !--- read the data from the source file call read_data(river_src_file, 'tocell', River%tocell, domain) where (River%tocell(:,:).eq. 4) River%tocell(:,:)=3 where (River%tocell(:,:).eq. 8) River%tocell(:,:)=4 where (River%tocell(:,:).eq. 16) River%tocell(:,:)=5 where (River%tocell(:,:).eq. 32) River%tocell(:,:)=6 where (River%tocell(:,:).eq. 64) River%tocell(:,:)=7 where (River%tocell(:,:).eq.128) River%tocell(:,:)=8 do j = jsc, jec do i = isc, iec if(River%tocell(i,j) > 0) then River%i_tocell(i,j) = i + di(River%tocell(i,j)) River%j_tocell(i,j) = j + dj(River%tocell(i,j)) end if end do end do 1 2 3 4 5 6 7 8 integer, parameter, dimension(8) :: di=(/1, 1, 0,-1,-1,-1,0,1/) integer, parameter, dimension(8) :: dj=(/0,-1,-1,-1, 0, 1,1,1/) tocell: 0 1 1 E 2 2 SE 4 3 S 8 4 SW 16 5 W 32 6 NW 64 7 N 128 8 NE ====================== On 04/09/13 09:16, Fred Hattermann wrote: > Das macht so ziemlich jedes GIS, z.b. GRASS mit dem Befehl > d.watershed. Wir haben auch ein Interface fuer SWIM, wo dies > errechnet wird. #ncap2 -s 'landheight=oro*topo;' maast_70Ma_vtop.nc /home/petri/wk/grassdata/maast_70Ma_vtop+landheight.nc Caution: when creating a new ``Mapset'' in GRASS, it is important to define apropriate region and resolution. The values are stored in the File WIND. For the Cretaceous case, I use proj: 0 zone: 0 north: 90 south: -90 east: 361.40625 west: -1.40625 cols: 129 rows: 64 e-w resol: 2.8125 n-s resol: 2.8125 ## import NetCDF into GRASS #r.in.gdal -o -e --v --o #input=NETCDF:"/home/petri/wk/grassdata/maast_70Ma_vtop+landheight.nc":landheight output=landheight # Compute drainage direction. # Once for land only, with ocean set to height zero. # And once for entire world, land+sea. Drainage direction for ocean # areas does not make sense, but this improve drainage direction # calculation at the land/ocean boundary. Further processing must then # obey the land/sea mask. # I would expect, that inside the land regions both invocations yield # identical results. # I could not use the given land/sea mask in NetCDF variable oro as # GRASS mask, because GRASS insists on importing that variable as # floating point values, while GRASS masks must be integer-valued # maps. Sigh. #r.watershed --v memory=30000 elevation=landheight \ # drainage=drainage_direction_landonly threshold=1 #r.watershed --v memory=30000 elevation=topo \ # drainage=drainage_direction_land+sea threshold=1 #drainage #Output map: drainage direction. Provides the "aspect" for each cell #measured CCW from East. Multiplying positive values by 45 will give #the direction in degrees that the surface runoff will travel from that #cell. The value 0 (zero) indicates that the cell is a depression area #(defined by the depression input map). Negative values indicate that #surface runoff is leaving the boundaries of the current geographic #region. The absolute value of these negative cells indicates the #direction of flow. # that gives: # 0 1 2 3 4 5 6 7 8 # - NE N NW W SW S SE E # r.out.gdal names the variables always BAND1. # There are hints in the documentation that it should be possible to # set option for writing metadata (metaopt=NAME=VALUE), but I could not get that to work. #r.out.gdal --v input=drainage_direction_landonly \ # format=netCDF \ # output=/home/petri/wk/grassdata/maast_70Ma-drainage_direction_landonly.nc #r.out.gdal --v input=drainage_direction_land+sea \ # format=netCDF \ # output=/home/petri/wk/grassdata/maast_70Ma-drainage_direction_land+sea.nc # the generated NetCDF files have 129 columns, longitudes 0, 2.8125, # ..., 357.1875, 360 . Apparently the first column is duplicated into # the last one. This now needs to be converted into something that land_lad/soil/rivers.F90 can read. For that we need to write a small program: Read the GRASS-generated drainage map and the original land/sea mask. For each land cell c, trace its outflow, until it reaches an ocean cell o. Record that o as discharge destination cell for c. Possible problems: along the map edges, r.watershed places negative drainage direction values to indicate that the drainage leaves the region. Apparently, that means that r.watershed does not ``wrap'' the world around the longitude? Yes. And, more important, not across the poles. That hits us with a bone crushing sound in Antarctica. icc -g -Wall \ -L/home/petri/netcdf-4.2.1.1-intel-nohdf5/lib \ -I/home/petri/netcdf-4.2.1.1-intel-nohdf5/include \ watershed-to-lad.c -lnetcdf -o watershed-to-lad Above suspicion turns out to be true. Moreover, tracing drainage across the map edges results in circular flow in several places, especially in antarctica, but also across the east/west boundary in two places. 6,0 +E -> 7,0 -S -> 71,0 +W -> 70,0 -S -> 6,0 ! sigh. 14,0 -S -> 78,0 -S -> 14,0 82,0 +E -> 83,0 -S -> 19,0 +W -> 18,0 +W -> 17,0 -S -> 81,0 +E -> 82,0 88,0 +E -> 89,0 +E -> 90,0 -S -> 26,0 +W -> 25,0 +W -> 24,0 -S -> 88,0 127,35 -E -> 0,35 +S -> 0,34 -W -> 127,34 +N -> 127,35 127,56 -E -> 0,56 +N -> 0,57 -W -> 127,57 +S -> 127,56 One possible solution might be to try and rotate the map , such that the poles fall into the ocean, before applying r.watershed. That should work with CDO remap operator. But how to rotate the result back onto standard grid? Interpolation of the drainage directions will not give the desired results. Another possibility might be to extend r.watershed in a way that it also obeyes the topography across the map borders. Sigh. Modification of r.watershed seems , lets say, difficult. Another approach: Extend the original map, by appending the western half to the eastern edge, the eastern half to the western edge, south half rotated and mirrored to southern edge, northern half rotated and mirrored to northern edge. Ignore lat / lon axis labels. Then import it into GRASS, run r.watershed on it, export, and run the inverse procedure, which cuts out the original area. Rotation and mirroring: lon_r = (lon + nlons/2) mod nlons lat_m = nlats/2 - lat (???) icc -g -Wall -L/home/petri/netcdf-4.2.1.1-intel-nohdf5/lib -I/home/petri/netcdf-4.2.1.1-intel-nohdf5/include expand-map.c -lnetcdf -o expand-map ./expand-map 070Ma/maast_70Ma_vtop.nc topo oro 070Ma/maast_70Ma_vtop-expanded.nc invoke grass. make sure that region is set correctly!!! r.in.gdal --v --o input=NETCDF:/iplex/01/climber3/petri/mom5.0.2/exp/Cretaceous/070Ma/maast_070Ma_vtop-expanded.nc:topo output=topo g.region -p projection: 3 (Latitude-Longitude) zone: 0 datum: ** unknown (default: WGS84) ** ellipsoid: sphere north: 90N south: 90S west: 0 east: 0 nsres: 1:24:22.5 ewres: 1:24:22.5 rows: 128 cols: 256 cells: 32768 r.watershed --v memory=30000 elevation=topo \ drainage=drainage_direction_land+sea threshold=1 r.out.gdal --v input=drainage_direction_land+sea \ format=netCDF \ output=/iplex/01/climber3/petri/mom5.0.2/exp/Cretaceous/070Ma/maast_70Ma-expanded-drainage_direction_land+sea.nc For some reasons, I still get loops when I first cut out the middle region of the drainage map and then convert it to land_lad format. However, it works fine, when I do the conversion on the expanded map, then cut out the middle region of the converted map. ./watershed-to-lad 070Ma/maast_70Ma-expanded-drainage_direction_land+sea.nc Band1 \ 070Ma/maast_70Ma_vtop-expanded.nc oro \ 070Ma/maast_70Ma_vtop.nc topo \ 070Ma/INPUT/river_destination_field #### Das Astronomie-Module kennt verschiedene Parameter, mit denen sich vermutlich auch die Verhaeltnisse zur Kreidezeit einstellen lassen. &ASTRONOMY_NML ECC = 1.671000000000000E-002, OBLIQ = 23.4390000000000 , PER = 102.932000000000 , PERIOD = 0, YEAR_AE = 1998, MONTH_AE = 9, DAY_AE = 23, HOUR_AE = 5, MINUTE_AE = 37, SECOND_AE = 0, NUM_ANGLES = 3600 / Dazu kommen die Parameter des EBM-Atmosphaeren-Moduls &ATMOSPHERE_NML T_BOT_ATM = 30.0000000000000 , SOLAR_CONSTANT = 1360.00000000000 , SEASONAL_SOLAR = T, DIFF_IS_UNIFORM = T, TG_INIT = 250.000000000000 , [..] / Parameter des land-properties-Moduls. Da wirds leicht unuebersichtlich. Ich habe jetzt mal folgende Parameter eingesetzt: ------------- &land_properties_nml veg_to_use = 'variable' !veg_index_constant = 3 dynamic_cover_type = .false. ! read one static time slice from cover_type_field.nc read_old_ascii_cover = .false. cover_dataset_init_year = 1860 cover_dataset_entry = 1860, 1, 1, 0, 0, 0, soil_to_use = 'constant' soil_index_constant = 2 use_glaciers = .true. use_climap = .false. ! dont use albedo_data.nc use_desert_albedo_map = .false. ! dont use albedo_data.nc use_single_geo = .true. ! dont use groundwater_residence_time_field but constant geo_res_time geo_res_time = 5184000.0 ! default value use_topo_rough = .true., ! use INPUT/maast_70ma_topography.nc as specified in topography_nml. max_topo_rough = 100.0, ! ... But is the Cretacious map resolution sufficient for this??? topo_rough_factor = 0.01, use_single_basin = .false. / -------------- Vergleich mit unseren beiden Lieblings-GFDL-Beispielen: EBM CM2M_coarse_BLING 070Ma &LAND_PROPERTIES_NML DO_ALL_MCM = F, F F VEG_TO_USE = constant, variable variable SOIL_TO_USE = constant, variable constant USE_GLACIERS = T, T T USE_CLIMAP = T, F F USE_DESERT_ALBEDO_MAP = T, T F DO_MCM_MASKING = F, F USE_SINGLE_BASIN = F, F F USE_SINGLE_GEO = T, F T I_DEST0 = 1, 1 J_DEST0 = 1, 1 VEG_INDEX_CONSTANT = 3, 3 SOIL_INDEX_CONSTANT = 2, 2 2 SOIL_INDEX_ICE_SUBSTITUTE = 1, 1 GEO_RES_TIME = 5184000.0, 5184000.0 5184000.0 T_RANGE = 10.0, 10.0 USE_TOPO_ROUGH = F, T T MAX_TOPO_ROUGH = 100.0, 100.0 100.0 TOPO_ROUGH_FACTOR = 1.0, 0.01 0.01 SFC_HEAT_FACTOR = 1.0, 0.25 NUM_SFC_LAYERS = 0, 6 0 DYNAMIC_COVER_TYPE = T, F F READ_OLD_ASCII_COVER = F, T F COVER_DATASET_INIT_YEAR = 1860, 1860 1860 COVER_DATASET_ENTRY = 1860,1,1,0,0,0,1860,1,1,0,0,0 1860,1,1,0,0,0 RECONCILE_LAKES = F, F F SOIL_INDEX_LAKE_SUBSTITUTE = 1 1 1 / ----------------- Damit laeuft jetzt die Initialisierung von EBM-Atmosphaere und land_lad durch, im Ozean bleibts haengen bei den Tracern: ---------------- ==>Note from ocean_tracer_mod(ocean_prog_tracer_init): Reading prognostic tracer initial conditions or restarts Initializing tracer number 1 at time level tau. This tracer is called temp FATAL: temp is not initialized as an ideal passive tracer, it is not initialized as constant, and it does not exist in the file INPUT/ocean_temp_salt.res.ncAll tracers must have initialization specified. There is no default. --- This is specified via the tracer name lists in field_table . We do not have a suitable restart file for temp and salt. To solve this, add options "prog_tracers","ocean_mod","temp" ... const_init_tracer = t const_init_value = 274.0 / "prog_tracers","ocean_mod","salt" ... const_init_tracer = t const_init_value = 10.0 / Those constant initial conditions can be modified subsequently: MOM5 has options for initializing with idealized conditions, subroutine ocean_tempsalt_ideal_reinit() in ocean_tracers/ocean_tempsalt.F90 &ocean_tempsalt_nml ! added for Cretaceous experiment reinit_ts_with_ideal = .true. / Also, comment out restoring in ocean_sbc_nml: &ocean_sbc_nml use_waterflux=.true. temp_restore_tscale=-10. !salt_restore_tscale=60. salt_restore_tscale=-10. ! do not use salt_sfc_restore.nc ... / ! ocean_shortwave_gfdl_mod wants either to read Chlorophyll data from file, ! or to compute it as tracer. For now, do neither. &ocean_shortwave_gfdl_nml use_this_module=.false. debug_this_module=.false. !read_chl=.true. read_chl=.false. ! avoid reading chl.nc zmax_pen=100.0 enforce_sw_frac=.true. optics_morel_antoine=.false. optics_manizza=.true. / &ocean_shortwave_nml use_this_module=.true. use_shortwave_gfdl=.false. ! see above use_shortwave_csiro=.false. / Next problem: ocean thickness. &ocean_thickness_nml debug_this_module=.false. debug_this_module_detail=.false. thickness_dzt_min_init=2.0 thickness_dzt_min=1.0 thickness_method='energetic' initialize_zero_eta=.false. read_rescale_rho0_mask=.false. ! in maast_070Ma map no isolated small basins like black sea !read_rescale_rho0_mask=.true. ! read basin_mask.nc !rescale_rho0_mask_gfdl=.true. !rescale_rho0_basin_label=7.0 ! black sea !rescale_rho0_value=.75 / After this modification, we get lots of NAN&Co in ocean_thickness.F90 subroutine update_tcell_thickness(), in variable Thickness%dzt_dst(:,:,:), and subsequently in Thickness%dzwt(:,:,:) Apparently, things are computed correctly for the indices, where the tracer grid mask is 1, but not elsewhere. But afterwards, the entire fields, regardless of tracer mask, are used for sunsequent calculations. Notably in ocean_nphysics_util.F90 subroutine tracer_derivs(), where the inverse value of Thickness%dzwt(:,:,:) is used to compute vertical derivatives. There, the tracer mask value is multiplied with 1/Thickness%dzwt(:,:,:), which then also operates on all the NANs, which the mask should normally mask out... Solution attempt: try to set up a basin mask, even though we (apparently) dont have basins. Or just copy the ocean_mask.nc file as dummy-basin-mask. ferret: use ocean_mask.nc can mode upcase let basin_mask=if mask ge 0.00001 then 1 else 0 save/file="basin_mask.nc" basin_mask Sigh. Does not help. What is Thickness%pbot0r ? Apparently that one is causing the trouble. The first ocean field, for which the Grid%tmask is 1, is at i=67,j=2,k=1 There, Thickness%dzt_dst is -0.00011, Dens%rho is 892.5, etc. pbot0r is 1/Thickness%pbot0 . It is set on thos cells, where Grid%kmt > 1. But Grid%kmt(67,2) is just 1.0, thus pbot0r remains at initial 0.0 there. Apparently that means that there is only one T-cell layer at that location, i.e. very shallow water. meine aktuelle These ist, dass MOM5 Probleme hat mit flachen Kuestengewaeesern, wo die Wassersaeule nur aus einer einzelnen Tracer-Zelle besteht. Kommt sowas in ``normalen'' Modellen nicht vor? In ocean_thickness.F90, subroutine thickness_initialize() wird der Bodendruck berechnet in Thickness%pbot0 . Spaeter wird oefter mal der Kehrwert gebraucht, in Thickness%pbot0r . Der wird aber nur dort berechnet, wo Grid%kmt > 1 , also mehr als eine T-Zelle in der Wassersaeule uebereinander liegt (wenn ich das richtig verstehe). ! tentative setting for reference bottom pressure (Pa) ! on T-cells, as well as its inverse. Will recompute ! this field in ocean_thickness_init_adjust [..] do j=jsd,jed do i=isd,ied if(Grid%kmt(i,j) > 1) then Thickness%pbot0r(i,j) = 1.0/Thickness%pbot0(i,j) endif enddo enddo In ocean_thickness_init_adjust() sieht die Berechnung dann genauso aus. Eigentlich wuerd ich doch erwarten, dass die Maskierungsbedingung heisst if(Grid%kmt(i,j) > 0) - oder hab ich was falsch verstanden? MH> ... in MOM-3 ... Der Grid-Generator hatte das gewaehrleistet. Ich MH> wuerde das Gitter Deines MOM-5 Modells manuell so veraendern, dass MH> keine "Ein-Zellen Wassersaeulen" vorkommen. HK> Am besten waere es also aus meiner Sicht, einfach alle Stellen, die nur HK> eine Box tief sind, auf 2 Tracer-Zellen zu setzen. Aarrgghhh. Bug in make_topog utility. File topog.c line 498 invocation of process_topog(): parameters 18 and 19 must be re-ordered (kmt_min and dont_change_landmask). -> Created issue in MOM5 bug tracker. Ich habe das nun gefixt, und dann eine vertikale Einteilung mit 26 Stufen, aehnlich wie beim CM2M_coarse_BLING, gebaut, mit linear 5m-Teilung in den oberen Schichten, und darunter dann groeber werdend. Dann erzaehlt mir ocean_topog_init ganz oft: ==>ocean_topog_init: WARNING: kmt( 121, 55) = 2 is generally not permitted. kmt>2 recommended for wet ocean domain. WARNING: ==>ocean_topog_init: Regions with 0 < kmt < 2 have been found. Recommend kmt > 2 for wet ocean regions. Also habe ich dem make_topog noch mitgegeben den Parameter --kmt_min 3. Damit jammert ocean_topog_init nicht mehr rum, aber dann: ==>Error from ocean_thickness: Surface undulations too negative; model unstable eta_t( 27 , 55 ) = -90.8861974383728 depth( 27 , 55 ) = 66.2495615481052 [.. noch ein paar aehnlich Meldungen fuer andere Positionen..] FATAL: ==>Error from ocean_thickness_mod: Free surface penetrating rock! Model unstable. Da geht wohl der Meeresspiegel unter dem Meeresboden durch. Seufz. MH> 5m fuer die oberen Boxen ist sehr sehr wenig! Bei 26 Schichten wuerde ich oben mindestens 20m nehmen. Ach ja, die 5m sind im ``Supergrid''. Der Ozean rechnet mit der halben Supergrid-Aufloesung, in dem Fall also 10 m. mmhh. At the location (27,55) shown above, there are 7 T-cells in the water column. Orographic depths is 66.24 meter. When eta_t reaches the -90.886 shown above, the corresponding values of Thickness%dz(25,55,1:7) show negative values, which is strange. Even stranger is that in the following time step, Thickness%dzt even is set to NaN in those 7 cells. update_tcell_thickness() in ocean_thickness.F90, around line 3204 elseif(vert_coordinate==PSTAR) then do k=1,nk do j=jsd,jed do i=isd,ied ! leave values over land untouched if(Grid%tmask(i,j,k) > 0.0) then Thickness%dzt_dst(i,j,k) = -grav_r/(Dens%rho(i,j,k,tau)+epsln) & *(Ext_mode%pbot_t(i,j,taup1)-Ext_mode%patm_t(i,j,taup1)) & *Thickness%pbot0r(i,j) Thickness%dztlo(i,j,k) = Thickness%dstlo(i,j,k)*Thickness%dzt_dst(i,j,k) Thickness%dztup(i,j,k) = Thickness%dstup(i,j,k)*Thickness%dzt_dst(i,j,k) Thickness%dzt(i,j,k) = Thickness%dztlo(i,j,k) + Thickness%dztup(i,j,k) Thickness%rho_dzt(i,j,k,taup1) = Dens%rho(i,j,k,tau)*Thickness%dzt(i,j,k) Thickness%rho_dztr(i,j,k) = 1.0/(Thickness%rho_dzt(i,j,k,taup1)+epsln) endif wrk1(i,j,k) = 1.0/Thickness%dzt_dst(i,j,k) enddo enddo enddo maybe this warning is related to the problem? WARNING: diag_manager_mod::send_data_3d: A value for ocean_model in field dzt (Min: -1.80333E+001, Max: 8.19243E+002) is outside the range [ -1.00000E+001, 1.00000E+005], and not equal to the missing value. WARNING: diag_manager_mod::send_data_3d: A value for ocean_model in field dzt (Min: -1.80333E+001, Max: 8.19243E+002) is outside the range [ -1.00000E+001, 1.00000E+005], and not equal to the missing value. Solution: file field_table section "prog_tracers","ocean_mod","temp" const_init_value = 3.65 The const_init_value is in degC, not in K . First symptom is that the initial average ocean density is similar to strong alcoholic drink. (Apparently the density calculation does not take into account boiling of ocean water). > ==>Note: From ocean_density_mod: Boussinesq reference density rho0(kg/m3) = > 0.103500000000E+04 > Initial rho_average(kg/m3) = 0.896007384138E+03 > Since rho0 .ne. rho_average, consider changing rho0 in > ocean_parameters.F90 to be equal to rho_average for better accuracy. The global average for salt in mom4p1_ebm1/INPUT/ocean_temp_salt.res.nc is 34.37 . The global average for temp in mom4p1_ebm1/INPUT/ocean_temp_salt.res.nc is 3.654 degC With those settings, ocean water mass in the initial stocks report looks realistic. But energy and salt content of the ocean still look bad. To fix that, the above explained option reinit_ts_with_ideal must be switched off in &ocean_tempsalt_nml. It causes the const values to be re-calculated with some exponentionl dependency on depth. Maybe we can utilize this feature lateron, after we better understand it. If the initial stock report for the land module shows strange numbers or even NaN, then the most likely reason is that the land cell area is read incorrectly from the grid spec. As distributed, MOM5.0.2 has a bug that makes the EBM model incompatible with mosaic-style grid_spec. The reason is confusion between F90 and C about the storage size of real*8 . my fix is to put an unconditional #undef OVERLOAD_R4 at the top of src/shared/moasic/read_mosaic.c . I am not usre if that is a general fix for everyone. At PIK we use Intel Fortran with option -r8 to make real*8 the default precision for type real. But the code for the EBM model is compiled with option -DOVERLOAD_R4, which caused read_mosaic.c to assume that the caller wants the data in single precision floats. ------------- Fri Apr 9 17:38:10 CEST 2021 The problem with read_mosaic.c still persists! NOTE: the EBM atmosphere does not consider CO2 nor any other form of radiative gases.