diff --git a/GEOSldas_App/GEOSldas_HIST.rc b/GEOSldas_App/GEOSldas_HIST.rc index d9d09dd..ef1a629 100644 --- a/GEOSldas_App/GEOSldas_HIST.rc +++ b/GEOSldas_App/GEOSldas_HIST.rc @@ -146,6 +146,7 @@ COLLECTIONS: 'TP6' , 'GridComp' , 'TSOIL6' , # ... 'PRLAND' , 'GridComp' , 'PRECTOTCORRLAND' , # assume "corrected" precip 'SNOLAND' , 'GridComp' , 'PRECSNOCORRLAND' , # assume "corrected" precip +>>>HIST_IRRIG<<< 'IRRLAND' , 'GridComp' , 'TSLAND' , 'GridComp' , 'SNOMASLAND' , 'SNOWDP' , 'GridComp' , 'SNODPLAND' , 'EVPSOI' , 'GridComp' , 'LHLANDSOIL' , @@ -214,7 +215,10 @@ COLLECTIONS: >>>HIST_CATCHCN<<< 'WAT10CM' , 'GridComp' , >>>HIST_CATCHCN<<< 'WATSOI' , 'GridComp' , >>>HIST_CATCHCN<<< 'ICESOI' , 'GridComp' , ->>>HIST_IRRIG<<< 'IRRIGRATE' , 'GridComp' , +>>>HIST_IRRIG<<< 'SPRINKLERRATE' , 'IRRIGATION' , +>>>HIST_IRRIG<<< 'FURROWRATE' , 'IRRIGATION' , +>>>HIST_IRRIG<<< 'FLOODRATE' , 'IRRIGATION' , +>>>HIST_IRRIG<<< 'DRIPRATE' , 'IRRIGATION' , :: tavg24_2d_lnd_Nx.format: 'CFIO', @@ -251,6 +255,7 @@ COLLECTIONS: 'TP6' , 'GridComp' , 'TSOIL6' , # ... 'PRLAND' , 'GridComp' , 'PRECTOTCORRLAND' , # assume "corrected" precip 'SNOLAND' , 'GridComp' , 'PRECSNOCORRLAND' , # assume "corrected" precip +>>>HIST_IRRIG<<< 'IRRLAND' , 'GridComp' , 'TSLAND' , 'GridComp' , 'SNOMASLAND' , 'SNOWDP' , 'GridComp' , 'SNODPLAND' , 'EVPSOI' , 'GridComp' , 'LHLANDSOIL' , @@ -319,7 +324,10 @@ COLLECTIONS: >>>HIST_CATCHCN<<< 'WAT10CM' , 'GridComp' , >>>HIST_CATCHCN<<< 'WATSOI' , 'GridComp' , >>>HIST_CATCHCN<<< 'ICESOI' , 'GridComp' , ->>>HIST_IRRIG<<< 'IRRIGRATE' , 'GridComp' , +>>>HIST_IRRIG<<< 'SPRINKLERRATE' , 'IRRIGATION' , +>>>HIST_IRRIG<<< 'FURROWRATE' , 'IRRIGATION' , +>>>HIST_IRRIG<<< 'FLOODRATE' , 'IRRIGATION' , +>>>HIST_IRRIG<<< 'DRIPRATE' , 'IRRIGATION' , :: diff --git a/GEOSldas_App/ldas_setup b/GEOSldas_App/ldas_setup index 1e795ce..45510b3 100755 --- a/GEOSldas_App/ldas_setup +++ b/GEOSldas_App/ldas_setup @@ -534,7 +534,14 @@ class LDASsetup: vegdynRstFile=tmpRstDir+'/'+tmpFile if not os.path.isfile(vegdynRstFile): assert int(self.rqdExeInp['RST_FROM_GLOBAL']) == 1, 'restart from LDASsa should be global' - + + tmpFile=self.rqdExeInp['RESTART_ID']+'.irrigation_internal_rst' + tmpRstDir=self.rqdExeInp['RESTART_PATH']+'/'.join([self.rqdExeInp['RESTART_ID'],'output', + self.rqdExeInp['RESTART_DOMAIN'],'rs',self.ensdirs[0]]) + irrigationRstFile=tmpRstDir+'/'+tmpFile + if not os.path.isfile(irrigationRstFile): + assert int(self.rqdExeInp['RST_FROM_GLOBAL']) == 1, 'restart from LDASsa should be global' + tmpFile=self.rqdExeInp['RESTART_ID']+'.landpert_internal_rst.'+y4m2d2_h2m2 tmpRstDir=self.rqdExeInp['RESTART_PATH']+'/'.join([self.rqdExeInp['RESTART_ID'],'output', self.rqdExeInp['RESTART_DOMAIN'],'rs',self.ensdirs[0],y4m2]) @@ -1037,38 +1044,50 @@ class LDASsetup: landice_obj = lake_landice_saltwater(config_obj = config) landice_obj.remap() - #for ens in self.ensdirs : - catchRstFile0 = '' - vegdynRstFile0 = '' - landiceRstFile0 = '' + # for ens in self.ensdirs : + catchRstFile0 = '' + vegdynRstFile0 = '' + landiceRstFile0 = '' + irrigationRstFile0 = '' + for iens in range(self.nens) : ensdir = self.ensdirs[iens] ensid = self.ensids[iens] myCatchRst = myRstDir+'/'+self.catch +ensid +'_internal_rst' myLandiceRst = myRstDir+'/'+ 'landice' +ensid +'_internal_rst' myVegRst = myRstDir+'/'+'vegdyn'+ensid +'_internal_rst' + myIrrRst = myRstDir+'/'+'irrigation'+ensid +'_internal_rst' myPertRst = myRstDir+'/'+ 'landpert' +ensid +'_internal_rst' catchRstFile = '' vegdynRstFile = '' + irrigationRstFile = '' pertRstFile = '' print ("restart: " + self.rqdExeInp['RESTART']) - + self.present_irrfile = glob.glob(self.bcs__dir_land+'irrigation_*.dat') if self.rqdExeInp['RESTART'].isdigit() : if int(self.rqdExeInp['RESTART']) == 0 or int(self.rqdExeInp['RESTART']) == 2 : vegdynRstFile = glob.glob(self.bcs_dir_land + 'vegdyn_*.dat')[0] catchRstFile = glob.glob(self.exphome+'/'+exp_id+'/mk_restarts/*'+self.catch+'_internal_rst.'+YYYYMMDD+'*')[0] + if len(self.present_irrfile) > 0: + irrigationRstFile = glob.glob(self.bcs_dir_land+'irrigation_*.dat')[0] else : # RESTART == 1 catchRstFile = rstpath+ensdir +'/'+ y4m2+'/'+self.rqdExeInp['RESTART_ID']+'.'+self.catch+'_internal_rst.'+y4m2d2_h2m2 vegdynRstFile= rstpath+ensdir +'/'+self.rqdExeInp['RESTART_ID']+ '.vegdyn_internal_rst' if not os.path.isfile(vegdynRstFile): # no vegdyn restart from LDASsa - if not os.path.isfile(vegdynRstFile0): - vegdynRstFile = glob.glob(self.bcs_dir_land + 'vegdyn_*.dat')[0] + if not os.path.isfile(vegdynRstFile0): + vegdynRstFile = glob.glob(self.bcs_dir_land + 'vegdyn_*.dat')[0] + irrigationRstFile= rstpath+ensdir +'/'+ y4m2+'/'+self.rqdExeInp['RESTART_ID']+'.irrigation_internal_rst.'+y4m2d2_h2m2 + if not os.path.isfile(irrigationRstFile): # no irrigation restart from LDASsa + if len(self.present_irrfile) > 0: + irrigationRstFile = glob.glob(self.bcs_dir_land+'irrigation_*.dat')[0] else : vegdynRstFile = glob.glob(self.bcs_dir_land + 'vegdyn_*.dat')[0] if self.with_land: - catchRstFile = glob.glob(self.exphome+'/'+exp_id+'/mk_restarts/*'+self.catch+'_internal_rst.'+YYYYMMDD+'*')[0] + catchRstFile = glob.glob(self.exphome+'/'+exp_id+'/mk_restarts/*'+self.catch+'_internal_rst.'+YYYYMMDD+'*')[0] + if len(self.present_irrfile) > 0: + irrigationRstFile = glob.glob(self.bcs_dir_land+'irrigation_*.dat')[0] # catchment restart file if os.path.isfile(catchRstFile) and self.with_land : @@ -1105,6 +1124,24 @@ class LDASsetup: vegdynRstFile0 = vegdynRstFile else : vegdynRstFile = vegdynRstFile0 + # irrigation restart file + if os.path.isfile(irrigationRstFile) : + irrigationLocal = self.rstdir+ensdir +'/'+ y4m2 +'/'+self.rqdExeInp['EXP_ID']+'.irrigation_internal_rst.'+y4m2d2_h2m2 + if self.islocal : + print ("Creating the local irrigation restart file... \n") + cmd='./preprocess_ldas.x zoomin_irrrst '+ irrigationRstFile +' ' + irrigationLocal + ' '+ tmp_f2g_file.name + + print ("cmd: " + cmd) + sp.call(shlex.split(cmd)) + else : + shutil.copy(irrigationRstFile,irrigationLocal) + + irrigationRstFile = irrigationLocal + + if '0000' in ensdir : + irrigationRstFile0 = irrigationRstFile + else : + irrigationRstFile = irrigationRstFile0 landiceRstFile = '' if self.with_landice : @@ -1143,9 +1180,12 @@ class LDASsetup: print ('vegdynRstFile: ' + vegdynRstFile) os.symlink(catchRstFile, myCatchRst) os.symlink(vegdynRstFile, myVegRst) + if os.path.isfile(irrigationRstFile): + os.symlink(irrigationRstFile, myIrrRst) if self.with_landice : print("link landice restart: " + myLandiceRst) os.symlink(landiceRstFile, myLandiceRst) + if ( self.has_geos_pert and self.perturb == 1 ): os.symlink(pertRstFile, myPertRst) @@ -1367,9 +1407,8 @@ class LDASsetup: ldasrcInp.pop('CO2_YEAR',None) ldasrcInp.pop('PRESCRIBE_DVG',None) - # create restart item in RC + # create restart item in RC catch_ = self.catch.upper() - if catch_+'_INTERNAL_RESTART_TYPE' in ldasrcInp : # avoid duplicate del ldasrcInp[ catch_ +'_INTERNAL_RESTART_TYPE'] @@ -1379,6 +1418,9 @@ class LDASsetup: if 'VEGDYN_INTERNAL_RESTART_TYPE' in ldasrcInp : # avoid duplicate del ldasrcInp['VEGDYN_INTERNAL_RESTART_TYPE'] + if 'IRRIGATION_INTERNAL_RESTART_TYPE' in ldasrcInp : + # avoid duplicate + del ldasrcInp['IRRIGATION_INTERNAL_RESTART_TYPE'] rstkey.append(catch_) rstkey.append('VEGDYN') @@ -1412,12 +1454,24 @@ class LDASsetup: keyn='LANDASSIM_OBSPERTRSEED_CHECKPOINT_FILE' valn='landassim_obspertrseed'+tmpl_+'_checkpoint' ldasrcInp[keyn]= valn + + for key,val in zip(rstkey,rstval) : + keyn = key+ '_INTERNAL_RESTART_FILE' + valn = '../input/restart/'+val+tmpl_+'_internal_rst' + ldasrcInp[keyn]= valn + rstkey=[catch_,'IRRIGATION'] + rstval=[self.catch,'irrigation'] + for key,val in zip(rstkey,rstval) : keyn = key+ '_INTERNAL_RESTART_FILE' valn = '../input/restart/'+val+tmpl_+'_internal_rst' ldasrcInp[keyn]= valn + keyn = key+ '_INTERNAL_CHECKPOINT_FILE' + valn = val+'_internal_checkpoint' + ldasrcInp[keyn]= valn + # checkpoint file and its type if self.with_land : keyn = catch_ + '_INTERNAL_CHECKPOINT_FILE' diff --git a/GEOSldas_App/lenkf_j_template.py b/GEOSldas_App/lenkf_j_template.py index 2909009..f00892c 100644 --- a/GEOSldas_App/lenkf_j_template.py +++ b/GEOSldas_App/lenkf_j_template.py @@ -740,6 +740,15 @@ /bin/ln -rs $tmp_file $EXPDIR/input/restart/${{rstf}}${{ENSID}}_rst endif + set rstf = 'irrigation' + echo ${{rstf}}${{ENSID}}_internal_checkpoint + if (-f ${{rstf}}${{ENSID}}_internal_checkpoint ) then + set tmp_file = $EXPDIR/output/$EXPDOMAIN/rs/$ENSDIR/Y${{eYEAR}}/M${{eMON}}/${{EXPID}}.${{rstf}}_internal_rst.${{eYEAR}}${{eMON}}${{eDAY}}_${{eHour}}${{eMin}} + /bin/mv ${{rstf}}${{ENSID}}_internal_checkpoint $tmp_file + /bin/rm -f $EXPDIR/input/restart/${{rstf}}${{ENSID}}_internal_rst + /bin/ln -s $tmp_file $EXPDIR/input/restart/${{rstf}}${{ENSID}}_internal_rst + endif + set rstf = 'landpert' if (-f ${{rstf}}${{ENSID}}_internal_checkpoint ) then set tmp_file = $EXPDIR/output/$EXPDOMAIN/rs/$ENSDIR/Y${{eYEAR}}/M${{eMON}}/${{EXPID}}.${{rstf}}_internal_rst.${{eYEAR}}${{eMON}}${{eDAY}}_${{eHour}}${{eMin}} @@ -760,6 +769,7 @@ set rstfiles2 = `ls landpert${{ENSID}}_internal_checkpoint.*` set rstfiles3 = `ls landassim_obspertrseed${{ENSID}}_checkpoint.*` set rstfiles4 = `ls landice${{ENSID}}_internal_checkpoint.*` + set rstfiles5 = `ls irrigation${{ENSID}}_internal_checkpoint.*` foreach rfile ( $rstfiles1 $rstfiles4 ) set ThisTime = `echo $rfile | rev | cut -d'.' -f2 | rev` @@ -791,6 +801,16 @@ /bin/mv $rfile ${{THISDIR}}${{EXPID}}.landassim_obspertrseed_rst.${{ThisTime}}.nc4 end + foreach rfile ( $rstfiles5 ) + set ThisTime = `echo $rfile | rev | cut -d'.' -f2 | rev` + set TY = `echo $ThisTime | cut -c1-4` + set TM = `echo $ThisTime | cut -c5-6` + set THISDIR = $EXPDIR/output/$EXPDOMAIN/rs/$ENSDIR/Y${{TY}}/M${{TM}}/ + if (! -e $THISDIR ) mkdir -p $THISDIR + /bin/mv $rfile ${{THISDIR}}${{EXPID}}.irrigation_internal_rst.${{ThisTime}}.nc4 + /usr/bin/gzip ${{THISDIR}}${{EXPID}}.irrigation_internal_rst.${{ThisTime}}.nc4 & + end + @ inens ++ end ## end of while ($inens < $NENS) wait diff --git a/GEOSldas_App/preprocess_ldas.F90 b/GEOSldas_App/preprocess_ldas.F90 index 67e3e9e..fe45af4 100644 --- a/GEOSldas_App/preprocess_ldas.F90 +++ b/GEOSldas_App/preprocess_ldas.F90 @@ -41,6 +41,8 @@ program main character(len=512) :: new_BC character(len=512) :: orig_Veg character(len=512) :: new_veg + character(len=512) :: orig_irr + character(len=512) :: new_irr character(len=512) :: orig_ease character(len=512) :: new_ease character(len=512) :: f2g_file @@ -100,7 +102,15 @@ program main new_veg = arg2 f2g_file = arg3 - call createZoominVegRestart(f2g_file, orig_veg, new_veg) + call createZoominVegRestart(f2g_file, orig_veg, new_veg) + + else if (trim(option) == "zoomin_irrrst") then + + orig_irr = arg1 + new_irr = arg2 + f2g_file = arg3 + + call createZoominRestart(f2g_file, orig_irr, new_irr, 100) else if (trim(option) == "zoomin_mwrtmrst") then @@ -108,7 +118,7 @@ program main new_rtm = arg2 f2g_file = arg3 - call createZoominRestart(f2g_file, orig_rtm, new_rtm, 100) + call createZoominRestart(f2g_file, orig_rtm, new_rtm, 100) else if (trim(option) == "zoomin_catchrst") then diff --git a/GEOSldas_App/preprocess_ldas_routines.F90 b/GEOSldas_App/preprocess_ldas_routines.F90 index 382a549..5e76f4a 100644 --- a/GEOSldas_App/preprocess_ldas_routines.F90 +++ b/GEOSldas_App/preprocess_ldas_routines.F90 @@ -1773,10 +1773,24 @@ subroutine createZoominRestart(mapping_file, orig_rst, new_rst, tile_type) call var_iter%next() cycle endif + + if (trim(vname) =='CROPCLASSNAME') then + call var_iter%next() + cycle + endif if (ndims == 1) then + call MAPL_VarRead (InFmt,vname,tmp1) + + ! set irrigation variables to zero (i.e., no remapping) [????? reichle, 13 Jun 2025] + if (trim(vname) == 'SPRINKLERRATE') tmp1 = 0. + if (trim(vname) == 'DRIPRATE' ) tmp1 = 0. + if (trim(vname) == 'FURROWRATE' ) tmp1 = 0. + if (trim(vname) == 'FLOODRATE' ) tmp1 = 0. + call MAPL_VarWrite(OutFmt,vname,tmp1(f2r_)) + else if (ndims == 2) then dname => var%get_ith_dimension(2) diff --git a/GEOSldas_App/tile_bin2nc4.F90 b/GEOSldas_App/tile_bin2nc4.F90 index 890bfd8..fda64ac 100644 --- a/GEOSldas_App/tile_bin2nc4.F90 +++ b/GEOSldas_App/tile_bin2nc4.F90 @@ -390,7 +390,12 @@ FUNCTION getAttribute (SHORT_NAME, LNAME, UNT) result (str_atr) case ('RMELTBC002'); LONG_NAME = 'flushed_out_black_carbon_mass_flux_from_the_bottom_layer_bin_2'; UNITS = 'kg m-2 s-1' case ('RMELTOC001'); LONG_NAME = 'flushed_out_organic_carbon_mass_flux_from_the_bottom_layer_bin_1'; UNITS = 'kg m-2 s-1' case ('RMELTOC002'); LONG_NAME = 'flushed_out_organic_carbon_mass_flux_from_the_bottom_layer_bin_2'; UNITS = 'kg m-2 s-1' - + case ('SPRINKLERRATE'); LONG_NAME = 'sprinkler_irrigation_rate'; UNITS = 'kg m-2 s-1' + case ('DRIPRATE'); LONG_NAME = 'drip_irrigation_rate'; UNITS = 'kg m-2 s-1' + case ('FURROWRATE'); LONG_NAME = 'furrow_irrigation_rate'; UNITS = 'kg m-2 s-1' + case ('FLOODRATE'); LONG_NAME = 'flood_irrigation_rate'; UNITS = 'kg m-2 s-1' + case ('IRRLAND'); LONG_NAME = 'Total_irrigation_land'; UNITS = 'kg m-2 s-1' + ! land constants case ('CDCR2'); LONG_NAME = 'maximum soil water content above wilting point'; UNITS = 'kg m-2' @@ -462,10 +467,10 @@ FUNCTION getAttribute (SHORT_NAME, LNAME, UNT) result (str_atr) case ('PRMC_ANA_ENSSTD'); LONG_NAME = 'soil_moisture_profile_analysis_ensstd'; UNITS = 'm3 m-3' case ('TSURF_ANA_ENSSTD'); LONG_NAME = 'surface_temperature_of_land_incl_snow_ensstd'; UNITS = 'K' case ('TSOIL1_ANA_ENSSTD'); LONG_NAME = 'soil_temperature_layer_1_analysis_ensstd'; UNITS = 'K' - + ! other land assimilation fields - case ('MWRTM_VEGOPACITY'); LONG_NAME = 'Lband_microwave_vegopacity_normalized_with_cos_inc_angle'; UNITS = '1' + case ('MWRTM_VEGOPACITY'); LONG_NAME = 'Lband_microwave_vegopacity_normalized_with_cos_inc_angle'; UNITS = '1' ! land ice fields