Wed Feb 13 09:29:55 CET 2019 Wed Jan 22 12:01:15 CET 2025 Updated for HPC2024 Precipitation patterns from CM2M_coarse_BLING are somewhat insufficient for coupling with LPJmL. Thus we try to use the higher-resolution CM2.1p1 together with LPJmL . The code is identical, we ``just'' have to merge the grid descriptions, INPUT/ and RESTART files. cd .../POEM/exp mkdir -p CM2.1p1_LPJ_05/INPUT cd CM2.1p1_LPJ_05/INPUT ln -s ../../CM2.1p1/INPUT/* . # use the 0.5 degree land grid rm grid_spec.nc atmos_mosaicXland_mosaic.nc land_hgrid.nc land_mosaic.nc land_mosaicXocean_mosaic.nc cp -av ../../CM2M_coarse_BLING_LPJ_05/INPUT/land_hgrid.nc . cp -av ../../CM2M_coarse_BLING_LPJ_05/INPUT/land_mosaic.nc . # the make_coupler_mosaic tool requires a dimension named ntiles # but this is not contained in the GFDL-provided topog.nc files # sigh. the name topog.nc is hardcoded in the ocean topography module ncap2 -o topog.nc -s 'defdim("ntiles",1)' ../../CM2.1p1/INPUT/topog.nc ../../../src/tools/make_coupler_mosaic/make_coupler_mosaic \ --ocean_mosaic ocean_mosaic.nc \ --land_mosaic land_mosaic.nc \ --atmos_mosaic atmos_mosaic.nc \ --ocean_topog topog.nc \ --check # created # ocean_mask.nc # land_mask.nc # atmos_mosaic_tile1Xland_mosaic_tile1.nc # atmos_mosaic_tile1Xocean_mosaic_tile1.nc # land_mosaic_tile1Xocean_mosaic_tile1.nc # mosaic.nc mv mosaic.nc grid_spec.nc #Now regrid the restart files for the land model. Keep fingers crossed that the coastlines are OK. #-> we get saturation vapor pressure table overflow at various points :-( #/p/projects/climber3/petri/POEM/src/tools/fregrid/fregrid \ # --input_mosaic /p/projects/climber3/petri/POEM/exp/CM2.1p1/INPUT/land_mosaic.nc \ # --input_dir /p/projects/climber3/petri/POEM/exp/CM2.1p1/INPUT/ \ # --input_file soil.res.nc \ # --output_mosaic /p/projects/climber3/petri/POEM/exp/CM2.1p1_LPJ_05/INPUT/land_mosaic.nc \ # --output_dir /p/projects/climber3/petri/POEM/exp/CM2.1p1_LPJ_05/INPUT/ \ # --output_file soil.res.nc \ # --scalar_field frac\,soil_temp_1\,soil_temp_2\,soil_temp_3\,soil_temp_4\,soil_temp_5\,soil_temp_6\,soil_temp_7\,soil_temp_8\,soil_temp_9\,soil_temp_10\,soil_temp_11\,soil_temp_12\,soil_temp_13\,soil_temp_14\,soil_temp_15\,soil_temp_16\,soil_temp_17\,soil_temp_18\,water\,snow\,groundwater\,fusion_1\,fusion_2\,fusion_3\,fusion_4\,fusion_5\,fusion_6\,fusion_7\,fusion_8\,fusion_9\,fusion_10\,fusion_11\,fusion_12\,fusion_13\,fusion_14\,fusion_15\,fusion_16\,fusion_17\,fusion_18 # # #/p/projects/climber3/petri/POEM/src/tools/fregrid/fregrid \ # --input_mosaic /p/projects/climber3/petri/POEM/exp/CM2.1p1/INPUT/land_mosaic.nc \ # --input_dir /p/projects/climber3/petri/POEM/exp/CM2.1p1/INPUT/ \ # --input_file vegetation.res.nc \ # --output_mosaic /p/projects/climber3/petri/POEM/exp/CM2.1p1_LPJ_05/INPUT/land_mosaic.nc \ # --output_dir /p/projects/climber3/petri/POEM/exp/CM2.1p1_LPJ_05/INPUT/ \ # --output_file vegetation.res.nc \ # --scalar_field q_ca Move the land restart files out of the way so that a cold start will be done: mv soil.res.nc soil.res.nc-merk mv vegetation.res.nc vegetation.res.nc-merk In field_table, add CO2 tracer for land and atmosphere --- CM2.1p1/INPUT/field_table 2009-12-16 21:09:59.000000000 +0100 +++ CM2.1p1_LPJ_05/field_table 2019-02-12 15:14:52.626007815 +0100 @@ -110,6 +110,24 @@ "TRACER", "atmos_mod", "cld_amt" "longname", "cloud fraction" "units", "none" / +# +# +# CO2 tracer [from ESM2M_pi-control_C2 setup] +## Note: atm transport and land model expect moist mmr. +## set co2 units: kg/kg (mass mixing ratio) for both atm and lnd field table entries. +## convert atm IC field table entry from dry vmr value to **dry** mmr +## Note: for atm cold start from constant co2 value, instead of using field table to initialize co2, +## create atmos_tracers.res.nc with co2 moist mmr values +## + "TRACER", "atmos_mod", "co2" + "longname", "carbon dioxide" + "units", "kg/kg" + "convection", "all" + "profile_type","fixed","surface_value = 434.5626E-06"/ + "TRACER", "land_mod", "co2" + "longname", "carbon dioxide" + "units", "kg/kg" / + # test tracer for radon # "TRACER", "atmos_mod", "radon" # "longname", "radon test tracer" ------------- Thu Feb 14 08:45:12 CET 2019 Wed Jan 22 12:01:15 CET 2025 Updated for HPC2024 Now prepare LPJmL input files. The self-made tool netcdf2grid.c reads a netcdf file like the land_mask.nc that is generated by make_coupler_mosaic and writes out an LPJ coordinates file like grid.bin, which contains one entry with lat/lon values for each land cell. ln -s ../CM2M_coarse_BLING_LPJ_05/Data-For-LPJ/netcdf2grid* . ./netcdf2grid INPUT/land_mask.nc grid.bin dimension 0 name ny len 360 dimension 1 name nx len 720 delta_lat 0.5 off_lat 0.25 delta_lon 0.5 off_lon 0.25 Number of grid cells: 92113 mkdir Soil cd Soil ln -s ../../CM2M_coarse_BLING_LPJ_05/Data-For-LPJ/Soil/HWSD* . ln -s ../../CM2M_coarse_BLING_LPJ_05/Data-For-LPJ/Soil/hwsd.* . ln -s ../../CM2M_coarse_BLING_LPJ_05/Data-For-LPJ/Soil/README.Soil . We must copy the scripts, because we need to modify e.g. the hard-coded number of grid cells, which was printed by the ./netcdf2grid invocation above cp -av ../../CM2M_coarse_BLING_LPJ_05/Data-For-LPJ/Soil/*.r . sed -i -e 's/^ *NCRUCELL* <- *[0-9]\+/NCRUCELL <- 92113/' create_soil_map.r sed -i -e 's/^ *NCRUCELL* <- *[0-9]\+/NCRUCELL <- 92113/' create_lake_frac.r This needs some memory: module load R/4.3.2 Rscript process_hwsd.r Fuer die Funktion TT.points.in.classes in create_soil_map.r muss das R Paket "soiltexture" geladen werden. Doku gibts hier: http://cran.r-project.org/web/packages/soiltexture/soiltexture.pdf R CMD 'install.packages(pkgs="/home/petri/tar/soiltexture_1.2.9.tar.gz",lib="/home/petri/R/library")' #export R_LIBS=/home/petri/R/library #wget https://cran.r-project.org/src/contrib/soiltexture_1.5.3.tar.gz #R CMD 'install.packages(pkgs="/home/petri/tar/soiltexture_1.5.3.tar.gz",lib="/home/petri/R/library")' Sigh. From 1.3.0 on, this package includes a tcltk-based GUI, which has additional dependencies. Thus we stick to v1.2.9 above. # the R scripts require library(fields). We already have that in the SacreX installation. export R_LIBS=/home/petri/R/library:/home/petri/wk/SacreX-Sciencehack/R/library-hpc24 ln -s ../grid.bin . #Rscript create_soil_map.r # #Sigh. #Error in frac.sea[ilon[i], ilat[i]] : subscript out of bounds ------------------- Thu Jan 30 10:54:04 CET 2020 Wed Jan 22 12:01:15 CET 2025 Updated for HPC2024 Werner has built a new tool for regridding the soil map. Also: since the R-scripts were built, the header format of LPJ CLM files was changed. Now implemented some code in create_soil_map.r and create_lake_frac.r to obeye the header version and adapt to it. Rscript create_soil_map.r script create_lake_frac.r -> successfully created soil_new.bin lake_frac.bin great_lakes_frac.bin ----------- Build a climatology for LPJ spin-up: Tue Jan 28 14:24:26 CET 2025 Next try: use the original CM2.1p1 code, i.e. with land_lad instead of land_atlantes, i.e. without Matteos snow age albedo (and without any trace of LPJ). Use the CM2.1p1 initial conditions on the above-created grid, where land_lad runs on 720x360 cells. Run it for 50 years. Then run for 30 years with daily output. lprec\,fprec\,lwnet\,swdown\,t_surf\,t_surf_max\,t_surf_min\,wind on land_lad grid. Then use the high-res RESTART files for soil and vegetation to start the coupled model. Georg suggested to compare not individual timesteps but climatology of the runs with different timesteps, to exclude comparison of local and temporal weather phenomena. That climatology comparison looks qite good. Thus we now take the 30-year daily output and prepare it as input for LPJ spinup: Add lprec and fprec to prec, and shift the years from 51..80 to 1901..1930 so that it matches the dates for the spinup configuration. mkdir history-for-LPJmL-spinup cdo -v shifttime\,1850year \ -merge \ -setattribute\,prec@long_name='total precipitation rate (input)' \ -chname\,fprec\,prec \ -add \ -selvar\,fprec 00810101.land_daily.nc \ -selvar\,lprec 00810101.land_daily.nc \ 00810101.land_daily.nc \ 19310101.land_daily+prec.nc #ncatted -a long_name\,prec\,m\,c\,'total precipitation rate (input)' 00820101.land_daily+prec.nc #Mon Feb 3 10:14:22 CET 2025 # #To create the spinup configuration, we first copy the original config file from the LPJmL source: # #( cd POEM/src/land_atlantes/LPJmL # git checkout poem_59_experimental # cp -av lpjml_fms.js \ # lpjml_fms_spinup.js \ # lpjml_netcdf.js \ # input_fms.js \ # input_fms_spinup.js \ # input_netcdf.js \ # param_fms.js \ # start_poem.sh \ # slurm_POEM.jcf \ # ../../../work/CM2.1p1_LPJ_05/ ----------------- Tue Feb 4 13:20:57 CET 2025 Lakes mkdir Lakes cd Lakes ln -s ../../CM2M_coarse_BLING_LPJ_05/Data-For-LPJ/Lakes/GLWD-level* . ln -s ../../CM2M_coarse_BLING_LPJ_05/Data-For-LPJ/Lakes/glwd_3.nc . cp -p ../../CM2M_coarse_BLING_LPJ_05/Data-For-LPJ/Lakes/README.Lakes . cp -p ../../CM2M_coarse_BLING_LPJ_05/Data-For-LPJ/Lakes/process_glwd.r . #edit process_hwsd.r to adapt the hard-coded number grid cells. For reading grid.bin, we need to replace the hard-coded skipping of 35 bytes of header information by the header-parsing code that Werner has implemented into create_soil_map.r. That also reads the actual number of grid cells from grid.bin . module load R/4.3.2 export R_LIBS=/home/petri/R/library:/home/petri/wk/SacreX-Sciencehack/R/library-hpc24 Rscript process_glwd.r Sigh. The generated glwd_lakes.png looks really strange. Most prominently it has the caspian sea on the southern hemisphere. At least, the latitude axis is reverted. od -A d -t u1 glwd_lakes.bin |less od -A d -t u1 -w 720 -v glwd_lakes.bin |less Werner says: the following command converts the glwd_lakes.bin into NetCDF: ../../../src/land_atlantes/LPJmL/bin/clm2cdf -raw -byte -notime -scale 0.01 glwd_lakes ./grid.bin glwd_lakes.bin glwd_lakes.nc Indeed, really strange. reversing the Y-axis in process_glwd.r script resulted in correct-looking output files. And also, running in CM2M_coarse_BLING_LPJ_05 the old process_glwd.r script with the old grid.bin, but with the new R interpreter, results now in poel-reverted output files. ----------- For a historical excursion about creation of river drainage maps from DDM30 see CM2M_coarse_BLING_LPJ_05/Data-For-LPJ/Rivers/README.Rivers For a quick start, we just re-use the drainage maps that Sibyll created, also for this configuration, and keep fingers crossed that they also work reasonably with the slightly change land-sea mask here. Sigh. We cant use Sibylls drainage map, because it contains the wrong number of cells for the new land-sea-mask. See below for creation of new drainage map. ------------- Wed Feb 5 14:02:31 CET 2025 Now we can come to a spin-up run with LPJ stand-alone. David remarks that the experimental branch with his improvements for the soil temperature computations gives bad results when used in stand-alone mode. Thus, for the spinup run, David and Werner recommend to use branch poem_master-5.8 . For the spinup we need some more files on the new land-sea-mask resp. grid.bin: no3deposition nh4deposition soilpH lightning human_ignition There is a tool to regrid input files. Values for newly-defined land cells are searched-for in the neighborhood. bin/regridclm -search -size4 \ /p/projects/climber3/petri/POEM/exp/CM2M_coarse_BLING_LPJ_05/Data-For-LPJ/grid.bin \ /p/projects/climber3/petri/POEM/exp/CM2.1p1_LPJ_05/grid.bin \ /p/projects/lpjml/input/POEM/nitrogen/no3_deposition_rcp8p5_fms.clm \ /p/projects/climber3/petri/POEM/exp/CM2.1p1_LPJ_05/no3_deposition_rcp8p5_fms.clm bin/regridclm -search -size4 \ /p/projects/climber3/petri/POEM/exp/CM2M_coarse_BLING_LPJ_05/Data-For-LPJ/grid.bin \ /p/projects/climber3/petri/POEM/exp/CM2.1p1_LPJ_05/grid.bin \ /p/projects/lpjml/input/POEM/nitrogen/nh4_deposition_rcp8p5_fms.clm \ /p/projects/climber3/petri/POEM/exp/CM2.1p1_LPJ_05/nh4_deposition_rcp8p5_fms.clm bin/regridclm -search \ /p/projects/climber3/petri/POEM/exp/CM2M_coarse_BLING_LPJ_05/Data-For-LPJ/grid.bin \ /p/projects/climber3/petri/POEM/exp/CM2.1p1_LPJ_05/grid.bin \ /p/projects/lpjml/input/POEM/nitrogen/soil_ph_fms.clm \ /p/projects/climber3/petri/POEM/exp/CM2.1p1_LPJ_05/soil_ph_fms.clm bin/regridclm -search -size4 \ /p/projects/climber3/petri/POEM/exp/CM2M_coarse_BLING_LPJ_05/Data-For-LPJ/grid.bin \ /p/projects/climber3/petri/POEM/exp/CM2.1p1_LPJ_05/grid.bin \ /p/projects/lpjml/input/POEM/dlightning_fms.clm \ /p/projects/climber3/petri/POEM/exp/CM2.1p1_LPJ_05/dlightning_fms.clm bin/regridclm -search \ /p/projects/climber3/petri/POEM/exp/CM2M_coarse_BLING_LPJ_05/Data-For-LPJ/grid.bin \ /p/projects/climber3/petri/POEM/exp/CM2.1p1_LPJ_05/grid.bin \ /p/projects/lpjml/input/POEM/human_ignition_fms.clm \ /p/projects/climber3/petri/POEM/exp/CM2.1p1_LPJ_05/human_ignition_fms.clm lpjcheck tells us: ERROR155: gridcells [0,90880] in '/p/projects/open/sibyll/POEM/exp/CM2M_coarse_BLING_LPJ_05/Data-For-LPJ/Rivers/drainage_DDM30_v1.04.clm' not in [0,92112]. --> must be created from scratch ! --> see below REMARK401: Longitudinal values>180 in '/p/projects/climber3/petri/POEM/exp/CM2.1p1_LPJ_05/history-for-LPJmL-spinup/00810101.land_daily+prec.nc', will be transformed. ERROR237: Last year=80 in '/p/projects/climber3/petri/POEM/exp/CM2.1p1_LPJ_05/history-for-LPJmL-spinup/00810101.land_daily+prec.nc' is less than last simulation year 1901. --> cdo shiftyear to 1901...1930 --> see above cdo showtimestamp /p/projects/climber3/petri/POEM/exp/CM2.1p1_LPJ_05/history-for-LPJmL-spinup/00810101.land_daily+prec.nc is really slow. Maybe we need again to convert the time axis from record axis to fixed axis? ----------------- Wed Feb 5 20:00:04 CET 2025 mkdir Rivers cd Rivers ln -s ../../CM2M_coarse_BLING_LPJ_05/Data-For-LPJ/Rivers/DDM30.asc . The map is generated with a cellsize of 30 arc minutes and is distributed as a zipped ASCII-grid that can be easily imported in most GIS-software that support rasters or grids. The direction of outflow from cell is according to the ArcInfo/ArcView convention as follows: ``0'' is assigned to internal sinks and to cells draining into the ocean while ``1'', ``2'', ``4'', ``8'', ``16'', ``32'', ``64'' and ``128'' standing for the outflow directions E, SE, S, SW, W, NW, N or NE. -9 is missing_value, for ocean cells. http://www.geo.uni-frankfurt.de/ipg/ag/dl/datensaetze/3_drainage_direction_map/index.html http://www.uni-frankfurt.de/45217896/3_drainage_direction_map? http://www.uni-frankfurt.de/45218101/DDM30 32 64 128 16 0 1 8 4 2 Convert the DDM30 convention to an ASCII art map that allows to get a graphical overview: ` ^ ' < o > / v \ sed -e "s;128;';g" -e 's;64;^;g' -e 's;32;`;g' -e 's;16;<;g' \ -e 's;8;/;g' -e 's;4;v;g' -e 's;2;\\;g' -e 's;1;>;g' \ -e 's;0;o;g' -e 's; ;;g' -e 's;-9;.;g' \ < DDM30.asc > DDM30.art.txt can mode upcase define axis/units=degrees_east/x=-179.75:179.75:0.5 x05deg define axis/y=89.75s:89.75n:0.5 y05deg define grid/x=x05deg/y=y05deg grid05deg set data/ez/variables=ddm_orig/grid=grid05deg/skip=6/columns=720 "DDM30.asc" ! this data set appears upside-down let ddm_rev=yreverse(ddm_orig) let/title="DDM30" ddm=ddm_rev[g=grid05deg@asn] set att ddm.missing_value=-9 ! sigh. for whatever reason ferret an attribute long_name_mod to some bogus value can att/output ddm.long_name_mod set var/outtype=short ddm save/clobber/file="DDM30.nc" ddm shade/levels=h ddm frame/file="DDM30.gif" ln -s ../grid.bin . ../../../src/land_atlantes/LPJmL-poem_master-5.8/bin/drainage grid.bin DDM30.asc drainage.bin This maps -9 -> -9 0 -> -1 direction 2^n -> index of destination cell I.e. LPJ cells with river destination -9 are ocean cells according to DDM30, LPJ cells with river destination -1 are inland sinks or coast cells. DDM30 may contain drainage directions from coast cells into adjacent ocean cells, but that information is discarded in the drainage program. Thu Feb 6 11:55:04 CET 2025 Sigh. This does not contain river data for Antarctica. To prepare that, I first added to src/land_lad/land_model.F90 a bit of code for diagnostic output of orography height, obtained from get_topog_mean, which regrids it to the land_lad grid, which in this case is also the LPJmL grid. (the code for writing out such a topography in map that I added earlier to lpj_driver.F90 is only invoked if LPJ is already coupled...). 00510102.lpj_horo.nc contains the resulting horo, together with the land-sea mask. Note: it has longitudes 0 to 360 See also ../Cretaceous/README_Cretaceous_MOM5 cd exp/CM2.1p1_LPJ_05/Rivers ../../Cretaceous/expand-map 00510102.lpj_horo.nc horo mask lpj_horo-expanded.nc Note: that also has longitudes 0 to 360 Now install grass GIS on your desktop machine, and copy the lpj_horo-expanded.nc to there. Here I have grass 7.8.7 invoke grass make sure that region is set correctly!!! g.region n=90 s=-90 e=359.875 w=-0.125 rows=1440 cols=720 r.in.gdal --v --o input=NETCDF:/home/stefan/grass-POEM/lpj_horo-expanded.nc:horo output=horo g.region raster=horo g.region -p r.watershed --v memory=30000 elevation=horo \ drainage=drainage_direction_land+sea threshold=1 r.out.gdal --v input=drainage_direction_land+sea \ format=netCDF \ output=/home/stefan/grass-POEM/lpj_horo-expanded-drainage_direction_land+sea.nc Note: this has longitudes running 0 to 360 ferret -> /home/stefan/grass-POEM/lpj_horo-expanded-drainage_direction_land+sea.gif copy the generated lpj_horo-expanded-drainage_direction_land+sea.nc and lpj_horo-expanded-drainage_direction_land+sea.gif back to the PIK cluster. ./watershed-to-DDM30 \ lpj_horo-expanded-drainage_direction_land+sea.nc Band1 \ lpj_horo-expanded.nc mask \ DDM-lpj_horo.asc Note: This runs in longitude from Greenwich to Greeenwich, while the original DDM30.asc runs can mode upcase !define axis/units=degrees_east/x=-179.75:179.75:0.5 x05deg define axis/units=degrees_east/x=0.25:359.75:0.5 x05deg define axis/y=89.75s:89.75n:0.5 y05deg define grid/x=x05deg/y=y05deg grid05deg set data/ez/variables=ddm_orig/grid=grid05deg/skip=6/columns=720 "DDM-lpj_horo.asc" ! this data set appears upside-down let ddm_rev=yreverse(ddm_orig) let/title="DDM_lpj_horo" ddm=ddm_rev[g=grid05deg@asn] set att ddm.missing_value=-9 ! sigh. for whatever reason ferret an attribute long_name_mod to some bogus value can att/output ddm.long_name_mod set var/outtype=short ddm save/clobber/file="DDM-lpj_horo.nc" ddm shade/levels=h ddm frame/file="DDM-lpj_horo.gif" CAUTION: this file is rotated in longitude, compared to the original DDM30.asc . Longitudes must run 0 to 360 here, the original DDM30.asc is stored on longitudes -180 to 180. A little awk script to swap right and left half of input lines awk '{ if (NF < 4) { print; next; } for (i = int((NF+1)/2+1); i <= NF; i++) printf("%s ", $i); if (NF%2 == 1) printf("%s ", $(NF/2+1)); for (i = 1; i <= int((NF)/2); i++) printf("%s ", $i); print(""); } ' < DDM-lpj_horo.asc > DDM-lpj_horo-shifted.asc sed -e 's;128;/;g' -e 's;64;^;g' -e 's;32;\\;g' -e 's;16;<;g' \ -e 's;8;/;g' -e 's;4;v;g' -e 's;2;\\;g' -e 's;1;>;g' \ -e 's;0;o;g' -e 's; ;;g' -e 's;-9; ;g' \ < DDM-lpj_horo-shifted.asc > DDM-lpj_horo-shifted.art.txt This now goes from date-line to date-line. can mode upcase define axis/units=degrees_east/x=-179.75:179.75:0.5 x05deg define axis/y=89.75s:89.75n:0.5 y05deg define grid/x=x05deg/y=y05deg grid05deg set data/ez/variables=ddm_orig/grid=grid05deg/skip=6/columns=720 "DDM-lpj_horo-shifted.asc" ! this data set appears upside-down let ddm_rev=yreverse(ddm_orig) let/title="DDM_lpj_horo-shifted" ddm=ddm_rev[g=grid05deg@asn] set att ddm.missing_value=-9 ! sigh. for whatever reason ferret an attribute long_name_mod to some bogus value can att/output ddm.long_name_mod set var/outtype=short ddm save/clobber/file="DDM-lpj_horo-shifted.nc" ddm shade/levels=h ddm frame/file="DDM-lpj_horo-shifted.gif" OK. Now we have to merge the Antarctica part of DDM-lpj_horo-shifted.asc with DDM30.asc . First lets find out which lines in the DDM-style maps have only -9 in them. awk '{ for (i=1; i <= NF; i++) { if ($i != "-9") { print(NR, $i); next; }} print(NR); }' < DDM30.asc In DDM30.asc , Southern-most information appears in line 198, lines 299-366 contain only -9. The above-created DDM-lpj_horo-shifted.asc contains -9-only in lines 297-312 Thus: (head -n 300 DDM30.asc ; tail -n +301 DDM-lpj_horo-shifted.asc ) > DDM30+antarctica-v0.0.asc ../LPJmL-poem_master-5.8/bin/drainage grid.bin DDM30+antarctica-v0.0.asc drainage+antarctica-v0.0.bin ----------------- Feb 10 10:32 clone LPJmL branch poem_master-5.8 into /p/projects/climber3/petri/POEM/exp/CM2.1p1_LPJ_05/LPJmL-poem_master-5.8/ Important: it must be compiled with option -DCOUPLING_WITH_FMS so that the resulting restart file can be used in couples runs. Change the input configuration to use all the above-generated files. ---------------- Mon Feb 10 14:57:48 CET 2025 With this drainage map, LPJmL standalone spinup starts well. On 10.02.25 15:46, Sibyll Schaphoff wrote: > wenn die Fluesse nicht im Ozean enden kuemmert das LPJ nicht, daher hat > das fuer das river routing in POEM wenig Bedeutung. But note that sea-ice and ocean need correct freshwater inflow for the currents. ln -s ../../CM2M_coarse_BLING_LPJ_05/Data-For-LPJ/drainage2cdfvec.c . make drainage2cdfvec ./drainage2cdfvec grid.bin drainage+antarctica-v0.0.bin drainage+antarctica-v0.0-vec.nc def sym greenland= X=-93:-17/Y=60:85 def sym northamerica= X=-168:-54/Y=22:80 def sym southamerica= X=-85:-35/Y=-60:15 def sym mideurope= X=-10:30/Y=36:60 def sym scandinavia= X=4:40/Y=55:72 def sym europe= X=-10:40/Y=36:72 def sym mediteran= X=-10:40/Y=30:45 def sym antarctica= X=0:360/Y=-90:-60 def sym africa= X=-20:65/Y=-40:35 def sym china= X=60:150/Y=20:70 def sym arabia= X=30:65/Y=10:35 def sym amazon= x=80w:40w/y=10s:12n def sym siberia= X=60:192/Y=60:80 use drainage+antarctica-v0.0-vec.nc shade/($amazon) dcount plot/vs/line/over/color=red/nolabel basinboundsx,basinboundsy frame/file=drainage+antarctica-v0.0-amazon.gif vector/over/xskip=1/yskip=1/color=red/aspect/len=2 dirx,diry etc. St Lawrence River does not really reach the Atlantic. Mississippi might be critical. Siberian River mouths do not look too good. Nile might be critical. Amazon looks OK. Tue Feb 11 13:51:35 CET 2025 Lets try to improve the runoff map. With above sed-script, convert the DDM30 format into a format that can be edited with emacs in picture-mode. tail -n +7 DDM30+antarctica-v0.0.asc | \ sed -e "s;128;';g" -e 's;64;^;g' -e 's;32;`;g' -e 's;16;<;g' \ -e 's;8;/;g' -e 's;4;v;g' -e 's;2;\\;g' -e 's;1;>;g' \ -e 's;0;o;g' -e 's; ;;g' -e 's;-9;.;g' \ > DDM30+antarctica-v0.1.art.txt Edit that file: - continue Lawrence river north-east-ward, convert some "0" cells to ocean. - continue MacKenzie river northwards until ocean - continue Nile northwards ( head -6 DDM30+antarctica-v0.0.asc sed -e "s;';128 ;g" \ -e 's;\^;64 ;g' \ -e 's;`;32 ;g' \ -e 's;<;16 ;g' \ -e 's;/;8 ;g' \ -e 's;v;4 ;g' \ -e 's;\\;2 ;g' \ -e 's;>;1 ;g' \ -e 's;o;0 ;g' \ -e 's;\.;-9 ;g' \ DDM30+antarctica-v0.1.art.txt ) > DDM30+antarctica-v0.1.asc ../LPJmL-poem_master-5.8/bin/drainage grid.bin DDM30+antarctica-v0.1.asc drainage+antarctica-v0.1.bin ./drainage2cdfvec grid.bin drainage+antarctica-v0.1.bin drainage+antarctica-v0.1-vec.nc Wed Feb 12 11:53:17 CET 2025 worked on Lena delta worked on Yenisei mouth worked on Ob mouth Fri Feb 14 15:37:32 CET 2025 worked on Mesopotamia worked on Greenland worked on South East Asia worked on India East worked on India West TODO: zoom into China coast regions X=110:125,Y=28:45 (west coast of yellow sea) Liao He Hai River (Beijing) Huang He (Yellow River) Yangtsekiang Mon Feb 17 12:10:30 CET 2025 Continued improvement of river mouths in East Siberia, South-East Asia, Western India, Africa, South America, Gulf of Mexico, Gulf of California, USA and Canada West coast -------------- Now we can try a coupled run. We use Davids latest experimental branch, which includes improved soil temperature calculation. ( cd POEM/src/land_atlantes/LPJmL git checkout poem_59_experimental cp -p \ input_fms.cjson \ lpjml_fms.cjson \ bin/start_poem.sh \ bin/createinput.sh \ bin/createlpjml.sh \ ../../../work/CM2.1p1_LPJ_05/ ) Note: branch poem_59_experimental must be compiled with -DRESTART_BELOW_5_9 so that it can read the restart files produced by above-used branch poem_master-5.8 . Note that createinput.sh , despite its name suffix, is really a template file for input.nml, in which the start_poem.sh script replaces several variables by the experiment-specific actually configured values. Thus: cd ../../../work/CM2.1p1_LPJ_05/ mv createinput.sh createinput.sh-CM2Mc cp -p input.nml createinput.sh-CM2.1p1 and then modify that new createinput.sh-CM2.1p1 template with the variable names, in the same fashion that you find in the original createinput.sh-CM2Mc On 19.02.25 12:10, Werner von Bloh wrote: > der soildepth Input wird von 5.9 nicht benutzt, ist fuer die variable > roots Version. grassland_lsuah wird nur fuer landuse benutzt, ist fuer den > online spinup nicht nowendig. Laesst sich aber per regriclm erzeugen. > setclm firstyear 0 /p/projects/climber3/petri/POEM/exp/CM2.1p1_LPJ_05/human_ignition_fms.clm Try to satisfy all the year-specifications in a consistent way. Since the CM2.1p1 climatology for the LPJmL spinup was created with 1990 forcings, we try that as starting year. cp -av cp -av INPUT-00510101 INPUT-19910101 sed -i -e 's/ 51/1991/' INPUT-19910101/couper.res rm INPUT ln -s INPUT-19910101 INPUT ----------------- Tue Feb 25 11:33:18 CET 2025 The first coupled runs crash with "Error from ocean_thickness_mod: Free surface penetrating rock! Model unstable." I guess that indicates problems with ocean water density. ocean_thickness.F90 update_tcell_thickness() type(ocean_grid_type), intent(in) :: Grid type(ocean_external_mode_type), intent(in) :: Ext_mode Ext_mode%eta_t < -Grid%ht i j task 54 28 200 depth 40 107.5E 89.5N Arctic sea (in the tripolar zone) task 55 60-62 155-157 depth 80-1068 139.5E 140.5E 141.5E 46.5N 45.5N 44.5N North-Japan sea task 58 189:217 168:174 depth 70-349 91.5W-63.5W 63.5N-57.5N Hudson Bay task 61 332-338 196-200 depth always 40 51.5E-57.5E 89.5N-85.5N Arctic Sea (in the tripolar zone) depth( 28 , 200 ) = 40.0000000000000 depth( 61 , 155 ) = 220.000000000000 depth( 62 , 155 ) = 70.0000000000000 depth( 60 , 156 ) = 1068.77367355403 depth( 62 , 156 ) = 80.0000000000000 depth( 61 , 157 ) = 506.413097905128 depth( 189 , 168 ) = 40.0000000000000 depth( 190 , 168 ) = 40.0000000000000 depth( 188 , 169 ) = 60.0000000000000 depth( 187 , 170 ) = 70.0000000000000 depth( 212 , 170 ) = 120.000000000000 depth( 214 , 170 ) = 130.000000000000 depth( 188 , 171 ) = 100.000000000000 depth( 212 , 171 ) = 120.000000000000 depth( 214 , 171 ) = 170.000000000000 depth( 215 , 171 ) = 170.000000000000 depth( 216 , 171 ) = 200.000000000000 depth( 217 , 171 ) = 349.182690022309 depth( 189 , 172 ) = 130.000000000000 depth( 213 , 172 ) = 130.000000000000 depth( 215 , 172 ) = 170.000000000000 depth( 216 , 172 ) = 200.000000000000 depth( 217 , 172 ) = 349.182690022309 depth( 190 , 173 ) = 110.000000000000 depth( 216 , 173 ) = 40.0000000000000 depth( 191 , 174 ) = 70.0000000000000 depth( 334 , 196 ) = 40.0000000000000 depth( 335 , 197 ) = 40.0000000000000 depth( 337 , 197 ) = 40.0000000000000 depth( 339 , 197 ) = 40.0000000000000 depth( 332 , 198 ) = 40.0000000000000 depth( 333 , 198 ) = 40.0000000000000 depth( 334 , 198 ) = 40.0000000000000 depth( 336 , 198 ) = 40.0000000000000 depth( 338 , 198 ) = 40.0000000000000 depth( 332 , 199 ) = 40.0000000000000 depth( 333 , 199 ) = 40.0000000000000 depth( 334 , 199 ) = 40.0000000000000 depth( 335 , 199 ) = 40.0000000000000 depth( 337 , 199 ) = 40.0000000000000 depth( 338 , 199 ) = 40.0000000000000 depth( 332 , 200 ) = 40.0000000000000 depth( 333 , 200 ) = 40.0000000000000 depth( 334 , 200 ) = 40.0000000000000 depth( 335 , 200 ) = 40.0000000000000 depth( 338 , 200 ) = 40.0000000000000 fms_MOM_ATLANTES_ 0000000002CDAFD5 mpp_mod_mp_mpp_er 58 mpp_util_mpi.inc fms_MOM_ATLANTES_ 00000000009989B9 ocean_thickness_m 3398 ocean_thickness.F90 fms_MOM_ATLANTES_ 0000000000476D24 ocean_model_mod_m 1821 ocean_model.F90 fms_MOM_ATLANTES_ 0000000000408B9B MAIN__ 708 coupler_main.F90 ocean_model.F90 call eta_and_pbot_update(Time, Ext_mode) call update_tcell_thickness(Time, Grid, Ext_mode, Dens, Thickness) ------------ Thu Mar 6 11:44:45 CET 2025 Try to run sequentially for debugging. That hangs during init of exchange grid :-( #0 0x00001555531817cd in MPID_Progress_wait (state=0x305e7f1c, req=0x7ffffffdc6e0) at ../../src/mpid/ch4/src/ch4_progress.c:192 #1 0x000015555370d6af in MPIR_Wait_impl (request_ptr=, status=) at ../../src/mpi/request/wait.c:36 #2 MPID_Wait (request_ptr=, status=) at ../../src/mpid/ch4/include/mpidpost.h:118 #3 MPIR_Wait (request=, status=) at ../../src/mpi/request/wait.c:93 #4 PMPI_Wait (request=0x305e7f1c, status=0x7ffffffdc6e0) at ../../src/mpi/request/wait.c:188 #5 0x0000155554a4a3ed in pmpi_wait_ (v1=0x305e7f1c, v2=0x7ffffffdc6e0, ierr=0x3060f720) at ../../src/binding/fortran/mpif_h/waitf.c:274 #6 0x0000000003c843f8 in mpp_mod::mpp_sync_self ( pelist=, check=, request=, msg_size=, msg_type=) at /p/projects/climber3/petri/POEM/src/shared/mpp/include/mpp_util_mpi.inc:226 #7 0x0000000003bb1fc4 in xgrid_mod::load_xgrid (xmap=..., grid=..., grid_file=..., grid1_id=..., grid_id=..., tile1=1, tile2=1, use_higher_order=4294967295, .tmp.GRID_FILE.len_V$2e87=256, .tmp.GRID1_ID.len_V$2e8a=3, .tmp.GRID_ID.len_V$2e8c=3) at /p/projects/climber3/petri/POEM/src/shared/exchange/xgrid.F90:983 #8 0x0000000003bcad5a in xgrid_mod::setup_xmap (xmap=..., grid_ids=..., grid_domains=..., grid_file=..., atm_grid=..., .tmp.GRID_FILE.len_V$3908=18) at /p/projects/climber3/petri/POEM/src/shared/exchange/xgrid.F90:1756 #9 0x0000000000431dcc in flux_exchange_mod::flux_exchange_init (time=..., atm=..., land=..., ice=..., ocean=..., ocean_state=0x96c9c260, atmos_ice_boundary=..., land_ice_atmos_boundary=..., land_ice_boundary=..., ice_ocean_boundary=..., ocean_ice_boundary=..., dt_atmos=1800, dt_cpld=7200) at /p/projects/climber3/petri/POEM/src/coupler/flux_exchange.F90:841 #10 0x00000000004141a9 in coupler_init () at /p/projects/climber3/petri/POEM/src/coupler/coupler_main.F90:1379 #11 0x00000000004078b3 in coupler_main () at /p/projects/climber3/petri/POEM/src/coupler/coupler_main.F90:442 load_xgrid: nxgrid > npes 88082 ---------------- Try with poem_master_58: fms.out-3043649-1-62-28 dt_fast 1800 swnet 0 lwnet 77.0934 sublimation 0.443694 latent_sublimation 2.834e+06 latent_fusion 334000 latent 2.264e+06 snow_melt 0 heat_snow 1.47822 evap 1.75069e-08 t_flux nan dhdt 86.7737 drdt 3.30272 dedt 1.93472e-13 soilheatcap(soil,0) 291624 forrtl: error (65): floating invalid Image PC Routine Line Source libpthread-2.28.s 0000146FF2CE5CE0 Unknown Unknown Unknown fms_MOM_ATLANTES_ 0000000002B5CCED soiltemp 157 soiltemp.c fms_MOM_ATLANTES_ 0000000002AF2E4D update_fast 112 update_fast.c fms_MOM_ATLANTES_ 0000000002AE1B18 lpj_update_ 1739 lpj_climber4.c fms_MOM_ATLANTES_ 0000000002A88A9B lpj_driver_mod_mp 1379 lpj_driver.F90 fms_MOM_ATLANTES_ 0000000002A13FCE land_model_mod_mp 841 land_model.F90 fms_MOM_ATLANTES_ 0000000000408149 MAIN__ 612 coupler_main.F90 fms_MOM_ATLANTES_ 00000000004075CD Unknown Unknown Unknown libc-2.28.so 0000146FF25C6CA3 __libc_start_main Unknown Unknown fms_MOM_ATLANTES_ 00000000004074EE Unknown Unknown Unknown lpj_climber4.c fast.t_flux = tmp_t_flux[cell]; lpj_driver.F90 subroutine lpj_driver_update_fast(..., 18=t_flux, ...) call lpj_update(... t_flux, ...) land_model.F90: subroutine update_land_model_fast ( atmos2land, bnd ) type(atmos_land_boundary_type), intent(in) :: atmos2land ! quantities exchanged between ! the atmosphere and the land call lpj_driver_update_fast(theLand%soil%time, & ..., atmos2land%t_flux(:,:,1), ... ! [W/m2] Sensible heat flux coupler_main.F90: type(atmos_land_boundary_type) :: Atmos_land_boundary do na = 1, num_atmos_calls call flux_down_from_atmos( Time_atmos, Atm, Land, Ice, & Land_ice_atmos_boundary, & Atmos_land_boundary, & Atmos_ice_boundary ) call update_land_model_fast( Atmos_land_boundary, Land ) flux_exchange.F90: subroutine sfc_boundary_layer call surface_flux (& ... ex_flux_t ...) subroutine flux_down_from_atmos (Time, Atm, Land, Ice, & Atmos_boundary, Land_boundary, Ice_boundary ) type(atmos_land_boundary_type), intent(inout):: Land_boundary where(ex_avail) ! temperature ex_gamma = 1./ (1.0 - ex_dtmass*(ex_dflux_t + ex_dhdt_atm*cp_inv)) ex_e_t_n = ex_dtmass*ex_dhdt_surf*cp_inv*ex_gamma ex_f_t_delt_n = (ex_delta_t + ex_dtmass * ex_flux_t*cp_inv) * ex_gamma ex_flux_t = ex_flux_t + ex_dhdt_atm * ex_f_t_delt_n ... end where call get_from_xgrid (Land_boundary%t_flux, 'LND', ex_flux_t, xmap_sfc) surface_flux.F90: subroutine surface_flux_1d (... flux_t ...) flux_t = rho_drag * (t_surf0 - th_atm) ! flux of sensible heat (W/m**2) The error occurs in Task 28: lpj_driver_update_fast 1991.01.01 05:00:00 goodmorning F goodnight F before lpj update number of NaNs in t_flux 0 cpl_fms_to_lpj(temp_mean) cpl_fms_to_lpj(prec) cpl_fms_to_lpj(wind) cpl_fms_to_lpj(swdown) cpl_fms_to_lpj(lwnet) cpl_fms_to_lpj(t_flux) cpl_fms_to_lpj: NaN output at cell i=1707 cpl_fms_to_lpj: NaN output at cell i=1709 cpl_fms_to_lpj: NaN output at cell i=1710 cpl_fms_to_lpj: NaN output at cell i=1711 cpl_fms_to_lpj: NaN output at cell i=1712 cpl_fms_to_lpj: NaN output at cell i=1713 cpl_fms_to_lpj: NaN output at cell i=1714 cpl_fms_to_lpj: NaN output at cell i=1715 cpl_fms_to_lpj: NaN output at cell i=1716 The NaN values come from Task 25: lpj_driver_update_fast 1991.01.01 05:00:00 goodmorning F goodnight F before lpj update number of NaNs in t_flux 9 cpl_fms_to_lpj(temp_mean) cpl_fms_to_lpj(prec) cpl_fms_to_lpj(wind) cpl_fms_to_lpj(swdown) cpl_fms_to_lpj(lwnet) cpl_fms_to_lpj(t_flux) cpl_fms_to_lpj: NaN at i 4371 xcoord 25 ycoord 37 cpl_fms_to_lpj: NaN at i 4373 xcoord 32 ycoord 37 cpl_fms_to_lpj: NaN at i 4374 xcoord 33 ycoord 37 cpl_fms_to_lpj: NaN at i 4375 xcoord 34 ycoord 37 cpl_fms_to_lpj: NaN at i 4376 xcoord 35 ycoord 37 cpl_fms_to_lpj: NaN at i 4377 xcoord 36 ycoord 37 cpl_fms_to_lpj: NaN at i 4378 xcoord 37 ycoord 37 cpl_fms_to_lpj: NaN at i 4379 xcoord 38 ycoord 37 cpl_fms_to_lpj: NaN at i 4380 xcoord 39 ycoord 37