diff --git a/GEOSldas_App/ldas_setup b/GEOSldas_App/ldas_setup index 7ff6646..d5e7c35 100755 --- a/GEOSldas_App/ldas_setup +++ b/GEOSldas_App/ldas_setup @@ -752,8 +752,8 @@ class LDASsetup: if self.isZoomIn: newlocalTile = tile+'.domain' print ("\nCreating local tile file :"+ newlocalTile) - print ("\n by excluding land type MAPL_Land_ExcludeFromDomain=1100...\n") - cmd = self.bindir +'/preprocess_ldas.x c_localtile ' + tile + ' ' + newlocalTile + ' '+ tmp_f2g_file.name + print ("\nExcluding tiles not in the domain by adding 1000 to the type ...\n") + cmd = self.bindir +'/preprocess_ldas.x zoomin_tile ' + tile + ' ' + newlocalTile + ' '+ tmp_f2g_file.name print ("cmd: " + cmd) sp.call(shlex.split(cmd)) short_tile=short_tile +'.domain' @@ -782,7 +782,7 @@ class LDASsetup: print ("Creating the boundary files for the simulation domain...\n") bcs_tmp=[] for bcf in bcs : - cmd = self.bindir +'/preprocess_ldas.x c_localbc ' + bcf + ' '+ bcf+'.domain' + ' '+ tmp_f2g_file.name + cmd = self.bindir +'/preprocess_ldas.x zoomin_bc ' + bcf + ' '+ bcf+'.domain' + ' '+ tmp_f2g_file.name print ("cmd: " + cmd) sp.call(shlex.split(cmd)) bcs_tmp=bcs_tmp+[bcf+'.domain'] @@ -910,6 +910,7 @@ class LDASsetup: 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' myPertRst = myRstDir+'/'+ 'landpert' +ensid +'_internal_rst' @@ -938,7 +939,7 @@ class LDASsetup: catchLocal = self.rstdir+ensdir +'/'+ y4m2+'/'+self.rqdExeInp['EXP_ID']+'.'+self.catch+'_internal_rst.'+y4m2d2_h2m2 if self.isZoomIn : print( "Creating local catchment restart file... \n") - cmd=self.bindir +'/preprocess_ldas.x c_localcatchrst '+ catchRstFile +' ' + catchLocal + ' '+ tmp_f2g_file.name + cmd=self.bindir +'/preprocess_ldas.x zoomin_catchrst '+ catchRstFile +' ' + catchLocal + ' '+ tmp_f2g_file.name print ("cmd: "+cmd) sp.call(shlex.split(cmd)) else : @@ -956,7 +957,7 @@ class LDASsetup: vegdynLocal = self.rstdir+ensdir +'/'+self.rqdExeInp['EXP_ID']+'.vegdyn_internal_rst' if self.isZoomIn : print ("Creating the local veg restart file... \n") - cmd=self.bindir + '/preprocess_ldas.x c_localvegrst '+ vegdynRstFile +' ' + vegdynLocal + ' '+ tmp_f2g_file.name + cmd=self.bindir + '/preprocess_ldas.x zoomin_vegrst '+ vegdynRstFile +' ' + vegdynLocal + ' '+ tmp_f2g_file.name print ("cmd: " + cmd) sp.call(shlex.split(cmd)) else : @@ -969,6 +970,24 @@ class LDASsetup: else : vegdynRstFile = vegdynRstFile0 + landiceRstFile = rstpath+ensdir +'/'+ y4m2+'/'+self.rqdExeInp['RESTART_ID']+'.'+'landice_internal_rst.'+y4m2d2_h2m2 + if os.path.isfile(landiceRstFile) : + landiceLocal = self.rstdir+ensdir +'/'+self.rqdExeInp['EXP_ID']+'.landice_internal_rst' + if self.isZoomIn : + print ("Creating coom in landice restart file... \n") + cmd=self.bindir + '/preprocess_ldas.x zoomin_landicerst '+ landiceRstFile +' ' + landiceLocal + ' '+ tmp_f2g_file.name + print ("cmd: " + cmd) + sp.call(shlex.split(cmd)) + else : + shutil.copy(landiceRstFile,landiceLocal) + + landiceRstFile = landiceLocal + + if '0000' in ensdir : + landiceRstFile0 = landiceRstFile + else : + landiceRstFile = landiceRstFile0 + if (self.has_geos_pert and self.perturb == 1) : pertRstFile = rstpath+ensdir +'/'+ y4m2+'/'+self.rqdExeInp['RESTART_ID']+'.landpert_internal_rst.'+y4m2d2_h2m2 pertLocal = self.rstdir+ensdir +'/'+ y4m2+'/'+self.rqdExeInp['EXP_ID']+'.landpert_internal_rst.'+y4m2d2_h2m2 @@ -979,6 +998,7 @@ class LDASsetup: print ('vegdynRstFile: ' + vegdynRstFile) os.symlink(catchRstFile, myCatchRst) os.symlink(vegdynRstFile, myVegRst) + os.symlink(landiceRstFile, myLandiceRst) if ( self.has_geos_pert and self.perturb == 1 ): os.symlink(pertRstFile, myPertRst) @@ -991,7 +1011,7 @@ class LDASsetup: mwRTMLocal = self.bcsdir+'/'+ y4m2+'/'+self.rqdExeInp['EXP_ID']+'.ldas_mwRTMparam.'+y4m2d2_h2m2+'z.nc4' if self.isZoomIn : print ("Creating the local mwRTM restart file... \n") - cmd= self.bindir +'/preprocess_ldas.x c_localmwrtmrst '+ mwRTMRstFile +' ' + mwRTMLocal + ' '+ tmp_f2g_file.name + cmd= self.bindir +'/preprocess_ldas.x zoomin_mwrtmrst '+ mwRTMRstFile +' ' + mwRTMLocal + ' '+ tmp_f2g_file.name print ("cmd: " + cmd) sp.call(shlex.split(cmd)) diff --git a/GEOSldas_App/preprocess_ldas.F90 b/GEOSldas_App/preprocess_ldas.F90 index a2dbc6f..460b795 100644 --- a/GEOSldas_App/preprocess_ldas.F90 +++ b/GEOSldas_App/preprocess_ldas.F90 @@ -6,11 +6,12 @@ program main use preprocess_ldas_routines, ONLY: & create_mapping, & - createLocalTilefile, & - createLocalBC, & - createLocalVegRestart, & - createLocalmwRTMRestart, & - createLocalCatchRestart, & + createZoominTilefile, & + createZoominBC, & + createZoominVegRestart, & + createZoominmwRTMRestart, & + createZoominCatchRestart, & + createZoominLandiceRestart, & correctEase, & convert_pert_rst, & optimize_latlon @@ -46,6 +47,8 @@ program main character(len=512) :: f2g_file character(len=12 ) :: ymdhm character(len=12 ) :: SURFLAY + character(len=:), allocatable :: new_r, orig_r + call get_command_argument(1,option) call get_command_argument(2,arg1) @@ -73,44 +76,52 @@ program main call create_mapping(orig_tile,domain_def_file,trim(out_path),catch_def_file,trim(exp_id),ymdhm, SURFLAY, f2g_file) - else if (trim(option) == "c_localtile") then + else if (trim(option) == "zoomin_tile") then orig_tile = arg1 new_tile = arg2 f2g_file = arg3 - call createLocalTilefile(f2g_file, orig_tile,new_tile) + call createZoominTilefile(f2g_file, orig_tile,new_tile) - else if (trim(option) == "c_localbc" ) then + else if (trim(option) == "zoomin_bc" ) then orig_BC = arg1 new_BC = arg2 f2g_file = arg3 - call createLocalBC(f2g_file, orig_BC, new_BC) + call createZoominBC(f2g_file, orig_BC, new_BC) - else if (trim(option) == "c_localvegrst") then + else if (trim(option) == "zoomin_vegrst") then orig_veg = arg1 new_veg = arg2 f2g_file = arg3 - call createLocalVegRestart(f2g_file, orig_veg, new_veg) + call createZoominVegRestart(f2g_file, orig_veg, new_veg) - else if (trim(option) == "c_localmwrtmrst") then + else if (trim(option) == "zoomin_mwrtmrst") then orig_rtm = arg1 new_rtm = arg2 f2g_file = arg3 - call createLocalmwRTMRestart(f2g_file, orig_rtm, new_rtm) + call createZoominmwRTMRestart(f2g_file, orig_rtm, new_rtm) - else if (trim(option) == "c_localcatchrst") then + else if (trim(option) == "zoomin_catchrst") then orig_catch = arg1 new_catch = arg2 f2g_file = arg3 - call createLocalCatchRestart(f2g_file, orig_catch, new_catch) + call createZoominCatchRestart(f2g_file, orig_catch, new_catch) + + else if (trim(option) == "zoomin_landicerst") then + + orig_r = trim(arg1) + new_r = trim(arg2) + f2g_file = trim(arg3) + + call createZoominLandiceRestart(f2g_file, orig_r, new_r) else if (trim(option)=="correctease") then diff --git a/GEOSldas_App/preprocess_ldas_routines.F90 b/GEOSldas_App/preprocess_ldas_routines.F90 index 0e23705..7288f24 100644 --- a/GEOSldas_App/preprocess_ldas_routines.F90 +++ b/GEOSldas_App/preprocess_ldas_routines.F90 @@ -86,11 +86,12 @@ module preprocess_ldas_routines private public :: create_mapping - public :: createLocalTilefile - public :: createLocalBC - public :: createLocalCatchRestart - public :: createLocalVegRestart - public :: createLocalmwRTMRestart + public :: createZoominTilefile + public :: createZoominBC + public :: createZoominCatchRestart + public :: createZoominLandiceRestart + public :: createZoominVegRestart + public :: createZoominmwRTMRestart public :: correctEase public :: optimize_latlon public :: convert_pert_rst @@ -1535,7 +1536,7 @@ end subroutine read_mapping ! ******************************************************************** - subroutine createLocalTilefile(mapping_file, orig_tile, new_tile) + subroutine createZoominTilefile(mapping_file, orig_tile, new_tile) implicit none character(*), intent(in) :: mapping_file @@ -1554,12 +1555,12 @@ subroutine createLocalTilefile(mapping_file, orig_tile, new_tile) character(len=4) :: typ_str, typ_str_exclude - character(len=*), parameter :: Iam = 'createLocalTilefile' + character(len=*), parameter :: Iam = 'createZoominTilefile' inquire(file=trim(orig_tile),exist=file_exist) if( .not. file_exist) stop ("original tile file does not exist") - ! Set default local tile file name + ! Set default Zoom in tile file name call read_mapping( mapping_file, N_types, tile_types=tile_types, N_tiles_r=N_tiles_r, N_tiles_f=N_tiles_f, f2r=f2r, r2g=r2g) if( all(N_tiles_r == N_tiles_f)) then print*, "Domain is the same, no need to create tile file..." @@ -1602,7 +1603,7 @@ subroutine createLocalTilefile(mapping_file, orig_tile, new_tile) if( any( tile_types == ty)) then ! here g_id is (consecutive) id of the global *land* tiles if (f2g(f_id) /= g_id) then - ! if tile is not in local domain, replace ty in "line" with 1000+ty" + ! if tile is not in Zoom in domain, replace ty in "line" with 1000+ty" write(typ_str, '(I0)') ty typ_str = adjustr(typ_str) n=index(line, typ_str) @@ -1620,11 +1621,11 @@ subroutine createLocalTilefile(mapping_file, orig_tile, new_tile) close(40) close(50) - end subroutine createLocalTilefile + end subroutine createZoominTilefile ! ******************************************************************** - subroutine createLocalBC(mapping_file, orig_BC, new_BC) + subroutine createZoominBC(mapping_file, orig_BC, new_BC) implicit none character(*),intent(in) :: mapping_file @@ -1657,11 +1658,125 @@ subroutine createLocalBC(mapping_file, orig_BC, new_BC) close(10) close(20) deallocate(tmpvec) - end subroutine createLocalBC + end subroutine createZoominBC ! ******************************************************************** - - subroutine createLocalCatchRestart(mapping_file, orig_catch, new_catch) + + subroutine createZoominLandiceRestart(mapping_file, orig_landice, new_landice) + + implicit none + character(*),intent(in):: mapping_file + character(*),intent(in):: orig_landice + character(*),intent(in):: new_landice + integer,parameter :: subtile=4 + integer :: filetype, rc,i, j, ndims + real,allocatable :: tmp1(:) + type(Netcdf4_FileFormatter) :: InFmt,OutFmt + type(FileMetadata) :: OutCfg + type(FileMetadata) :: InCfg + integer :: dim1,dim2 + type(StringVariableMap), pointer :: variables + type(Variable), pointer :: var + type(StringVariableMapIterator) :: var_iter + type(StringVector), pointer :: var_dimensions + character(len=:), pointer :: vname,dname + integer :: n, N_landicer, N_landicef, N_types, r_starts, f_starts + + integer,dimension(:),allocatable :: f2r, r2g, f2r_landice, tile_types, N_tiles_r, N_tiles_f + logical :: have_landice + + call read_mapping( mapping_file, N_types, tile_types=tile_types, N_tiles_r=N_tiles_r, N_tiles_f=N_tiles_f, f2r=f2r, r2g=r2g) + + have_landice = .false. + do n = 1, N_types + if (tile_types(n) == MAPL_LANDICE) then + have_landice = .true. + N_landicef = N_tiles_f(n) + N_landicer = N_tiles_r(n) + exit + endif + enddo + if ( .not. have_landice) return + if (N_landicer == N_landicef) return + + r_starts = sum(N_tiles_r(1:n-1)) + f_starts = sum(N_tiles_f(1:n-1)) + + f2r_landice = f2r( f_starts+1: f_starts+N_landicef) + f2r_landice = f2r_landice - r_starts + + allocate(tmp1(N_landicer)) + + ! check file type + + call MAPL_NCIOGetFileType(orig_landice, filetype,rc=rc) + + if (filetype /= 0) then + print*, "Does not support binary landice restart" + else + + ! filetype = 0 : nc4 output file will also be nc4 + + call InFmt%open(trim(orig_landice), pFIO_READ,rc=rc) + InCfg = InFmt%read(rc=rc) + OutCfg = InCfg + + call OutCfg%modify_dimension('tile', size(f2r_landice), rc=rc) + + call OutFmt%create(trim(new_landice),rc=rc) + call OutFmt%write(OutCfg,rc=rc) + + variables => InCfg%get_variables() + var_iter = variables%begin() + do while (var_iter /= variables%end()) + + vname => var_iter%key() + var => var_iter%value() + var_dimensions => var%get_dimensions() + + ndims = var_dimensions%size() + + if (trim(vname) =='time') then + call var_iter%next() + cycle + endif + + if (ndims == 1) then + call MAPL_VarRead (InFmt,vname,tmp1) + call MAPL_VarWrite(OutFmt,vname,tmp1(f2r_landice)) + else if (ndims == 2) then + + dname => var%get_ith_dimension(2) + dim1=InCfg%get_dimension(dname) + do j=1,dim1 + call MAPL_VarRead ( InFmt,vname,tmp1 ,offset1=j) + call MAPL_VarWrite(OutFmt,vname,tmp1(f2r_landice),offset1=j) + enddo + else if (ndims == 3) then + + dname => var%get_ith_dimension(2) + dim1=InCfg%get_dimension(dname) + dname => var%get_ith_dimension(3) + dim2=InCfg%get_dimension(dname) + do i=1,dim2 + do j=1,dim1 + call MAPL_VarRead ( InFmt,vname,tmp1 ,offset1=j,offset2=i) + call MAPL_VarWrite(OutFmt,vname,tmp1(f2r_landice) ,offset1=j,offset2=i) + enddo + enddo + end if + call var_iter%next() + + enddo + call inFmt%close(rc=rc) + call OutFmt%close(rc=rc) + end if ! file type nc4 + print*, "done create Zoom in landice restart" + end subroutine createZoominLandIceRestart + +! *************************************** + + subroutine createZoominCatchRestart(mapping_file, orig_catch, new_catch) implicit none character(*),intent(in):: mapping_file @@ -1791,12 +1906,12 @@ subroutine createLocalCatchRestart(mapping_file, orig_catch, new_catch) call inFmt%close(rc=rc) call OutFmt%close(rc=rc) end if ! file type nc4 - print*, "done create local catchment restart" - end subroutine createLocalCatchRestart + print*, "done create Zoom in catchment restart" + end subroutine createZoominCatchRestart ! ******************************************************************** - subroutine createLocalmwRTMRestart(mapping_file, orig_mwrtm, new_mwrtm) + subroutine createZoominmwRTMRestart(mapping_file, orig_mwrtm, new_mwrtm) implicit none character(*),intent(in):: mapping_file @@ -1847,11 +1962,11 @@ subroutine createLocalmwRTMRestart(mapping_file, orig_mwrtm, new_mwrtm) call inFmt%close(rc=rc) call OutFmt%close(rc=rc) - end subroutine createLocalmwRTMRestart + end subroutine createZoominmwRTMRestart ! ******************************************************************** - subroutine createLocalVegRestart(mapping_file, orig_veg, new_veg) + subroutine createZoominVegRestart(mapping_file, orig_veg, new_veg) implicit none character(*),intent(in):: mapping_file @@ -1929,7 +2044,7 @@ subroutine createLocalVegRestart(mapping_file, orig_veg, new_veg) deallocate(tmp) endif - end subroutine createLocalVegRestart + end subroutine createZoominVegRestart ! ********************************************************************