!---------------------------------------------------------------------- ! NEMO system team, System and Interface for oceanic RElocable Nesting !---------------------------------------------------------------------- ! ! ! PROGRAM: merge_bathy ! ! DESCRIPTION: !> @brief !> This program merge bathymetry file at boundaries. !> !> @details !> Coarse grid Bathymetry is interpolated on fine grid. !> Then fine Bathymetry and refined coarse bathymetry are merged at boundaries. !> !> BathyFine= weight * BathyCoarse + (1-weight)*BathyFine !> !> The weight function used is : 0.5 + 0.5*COS( (pi*dist) / width ) !> !> @author !> J.Paul ! REVISION HISTORY: !> @date Nov, 2013 - Initial Version ! !> @note Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) !> !> @todo !---------------------------------------------------------------------- !> @code PROGRAM merge_bathy USE netcdf ! nf90 library USE global ! global variable USE phycst ! physical constant USE kind ! F90 kind parameter USE logger ! log file manager USE fct ! basic useful function USE date ! date manager USE att ! attribute manager USE dim ! dimension manager USE var ! variable manager USE file ! file manager USE boundary ! boundary manager USE iom ! I/O manager USE dom ! domain manager USE grid ! grid manager USE extrap ! extrapolation manager USE interp ! interpolation manager USE filter ! filter manager USE mpp ! MPP manager USE iom_mpp ! MPP I/O manager IMPLICIT NONE ! local variable CHARACTER(LEN=lc) :: cl_namelist CHARACTER(LEN=lc) :: cl_date INTEGER(i4) :: il_narg INTEGER(i4) :: il_status INTEGER(i4) :: il_fileid INTEGER(i4) :: il_attid INTEGER(i4) :: il_imin0 INTEGER(i4) :: il_imax0 INTEGER(i4) :: il_jmin0 INTEGER(i4) :: il_jmax0 INTEGER(i4) , DIMENSION(ip_maxdim) :: il_rho INTEGER(i4) , DIMENSION(2,2) :: il_offset INTEGER(i4) , DIMENSION(2,2,2) :: il_ind ! INTEGER(i4) , DIMENSION(:,:,:,:), ALLOCATABLE :: il_value LOGICAL :: ll_exist REAL(dp) :: dl_fill REAL(dp) , DIMENSION(:,:,:,:), ALLOCATABLE :: dl_refined REAL(dp) , DIMENSION(:,:,:,:), ALLOCATABLE :: dl_weight TYPE(TFILE) :: tl_bathy0 TYPE(TFILE) :: tl_bathy1 TYPE(TFILE) :: tl_fileout TYPE(TATT) :: tl_att TYPE(TVAR) :: tl_var TYPE(TVAR) :: tl_lon TYPE(TVAR) :: tl_lat ! TYPE(TVAR) :: tl_depth ! TYPE(TVAR) :: tl_time TYPE(TDIM) , DIMENSION(ip_maxdim) :: tl_dim TYPE(TBDY) , DIMENSION(ip_ncard) :: tl_bdy ! loop indices INTEGER(i4) :: ji INTEGER(i4) :: jk INTEGER(i4) :: jl ! namelist variable CHARACTER(LEN=lc) :: cn_logfile = 'merge_bathy.log' CHARACTER(LEN=lc) :: cn_verbosity = 'warning' CHARACTER(LEN=lc) :: cn_bathy0 = '' INTEGER(i4) :: in_perio0 = -1 CHARACTER(LEN=lc) :: cn_bathy1 = '' INTEGER(i4) :: in_perio1 = -1 CHARACTER(LEN=lc) :: cn_varcfg = 'variable.cfg' CHARACTER(LEN=lc), DIMENSION(ig_maxvar) :: cn_varinfo = '' INTEGER(i4) :: in_imin0 = 0 INTEGER(i4) :: in_imax0 = 0 INTEGER(i4) :: in_jmin0 = 0 INTEGER(i4) :: in_jmax0 = 0 INTEGER(i4) :: in_rhoi = 0 INTEGER(i4) :: in_rhoj = 0 LOGICAL :: ln_north = .TRUE. LOGICAL :: ln_south = .TRUE. LOGICAL :: ln_east = .TRUE. LOGICAL :: ln_west = .TRUE. LOGICAL :: ln_oneseg= .TRUE. CHARACTER(LEN=lc) :: cn_north = '' CHARACTER(LEN=lc) :: cn_south = '' CHARACTER(LEN=lc) :: cn_east = '' CHARACTER(LEN=lc) :: cn_west = '' CHARACTER(LEN=lc) :: cn_fileout = 'bathy_merged.nc' !------------------------------------------------------------------- NAMELIST /namlog/ & !< logger namelist & cn_logfile, & !< log file & cn_verbosity !< log verbosity NAMELIST /namcfg/ & !< config namelist & cn_varcfg !< variable configuration file NAMELIST /namcrs/ & !< coarse grid namelist & cn_bathy0, & !< bathymetry file & in_perio0 !< periodicity index NAMELIST /namfin/ & !< fine grid namelist & cn_bathy1, & !< bathymetry file & in_perio1 !< periodicity index NAMELIST /namvar/ & !< variable namelist & cn_varinfo !< list of variable and interpolation !< method to be used. !< (ex: 'votemper|linear','vosaline|cubic' ) NAMELIST /namnst/ & !< nesting namelist & in_imin0, & !< i-direction lower left point indice on coarse grid & in_imax0, & !< i-direction upper right point indice on coarse grid & in_jmin0, & !< j-direction lower left point indice on coarse grid & in_jmax0, & !< j-direction upper right point indice on coarse grid & in_rhoi, & !< refinement factor in i-direction & in_rhoj !< refinement factor in j-direction NAMELIST /nambdy/ & !< boundary namelist & ln_north, & !< use north boundary & ln_south, & !< use south boundary & ln_east , & !< use east boundary & ln_west , & !< use west boundary & cn_north, & !< north boundary indices on fine grid & cn_south, & !< south boundary indices on fine grid & cn_east , & !< east boundary indices on fine grid & cn_west , & !< west boundary indices on fine grid & ln_oneseg !< use only one segment for each boundary or not NAMELIST /namout/ & !< output namelist & cn_fileout !< fine grid merged bathymetry file !------------------------------------------------------------------- !1- namelist !1-1 get namelist il_narg=COMMAND_ARGUMENT_COUNT() !f03 intrinsec IF( il_narg/=1 )THEN PRINT *,"MERGE BATHY: ERROR. need a namelist" STOP ELSE CALL GET_COMMAND_ARGUMENT(1,cl_namelist) !f03 intrinsec ENDIF !1-2 read namelist INQUIRE(FILE=TRIM(cl_namelist), EXIST=ll_exist) IF( ll_exist )THEN il_fileid=fct_getunit() OPEN( il_fileid, FILE=TRIM(cl_namelist), & & FORM='FORMATTED', & & ACCESS='SEQUENTIAL', & & STATUS='OLD', & & ACTION='READ', & & IOSTAT=il_status) CALL fct_err(il_status) IF( il_status /= 0 )THEN PRINT *,"MERGE BATHY: ERROR opening "//TRIM(cl_namelist) STOP ENDIF READ( il_fileid, NML = namlog ) !1-2-1 define log file CALL logger_open(TRIM(cn_logfile),TRIM(cn_verbosity)) CALL logger_header() READ( il_fileid, NML = namcfg ) !1-2-2 get variable extra information CALL var_def_extra(TRIM(cn_varcfg)) READ( il_fileid, NML = namcrs ) READ( il_fileid, NML = namfin ) READ( il_fileid, NML = namvar ) !1-2-3 add user change in extra information CALL var_chg_extra(cn_varinfo) READ( il_fileid, NML = namnst ) READ( il_fileid, NML = nambdy ) READ( il_fileid, NML = namout ) CLOSE( il_fileid, IOSTAT=il_status ) CALL fct_err(il_status) IF( il_status /= 0 )THEN CALL logger_error("MERGE BATHY: ERROR closing "//TRIM(cl_namelist)) ENDIF ELSE PRINT *,"MERGE BATHY: ERROR. can not find "//TRIM(cl_namelist) ENDIF !2- open files IF( TRIM(cn_bathy0) /= '' )THEN tl_bathy0=file_init(TRIM(cn_bathy0),id_perio=in_perio0) CALL iom_open(tl_bathy0) ELSE CALL logger_fatal("MERGE BATHY: can not find coarse grid bathymetry "//& & "file. check namelist") ENDIF IF( TRIM(cn_bathy1) /= '' )THEN tl_bathy1=file_init(TRIM(cn_bathy1),id_perio=in_perio1) CALL iom_open(tl_bathy1) ELSE CALL logger_fatal("MERGE BATHY: can not find fine grid bathymetry "//& & "file. check namelist") ENDIF !3- check !3-1 check output file do not already exist INQUIRE(FILE=TRIM(cn_fileout), EXIST=ll_exist) IF( ll_exist )THEN CALL logger_fatal("CREATE BATHY: output file "//TRIM(cn_fileout)//& & " already exist.") ENDIF !3-2 check namelist !3-2-1 check refinament factor il_rho(:)=1 IF( in_rhoi < 1 .OR. in_rhoj < 1 )THEN CALL logger_error("MERGE BATHY: invalid refinement factor."//& & " check namelist "//TRIM(cl_namelist)) ELSE il_rho(jp_I)=in_rhoi il_rho(jp_J)=in_rhoj ENDIF !3-2-2 check domain indices IF( in_imin0 < 1 .OR. in_imax0 < 1 .OR. in_jmin0 < 1 .OR. in_jmax0 < 1)THEN ! compute coarse grid indices around fine grid il_ind(:,:,:)=grid_get_coarse_index(tl_bathy0, tl_bathy1 ) il_imin0=il_ind(1,1,1) ; il_imax0=il_ind(1,2,1) il_jmin0=il_ind(2,1,1) ; il_jmax0=il_ind(2,2,1) il_offset(:,:)=il_ind(:,:,2) ELSE il_imin0=in_imin0 ; il_imax0=in_imax0 il_jmin0=in_jmin0 ; il_jmax0=in_jmax0 il_offset(1,:)=FLOOR(REAL(il_rho(jp_I)-1,dp)*0.5) il_offset(2,:)=FLOOR(REAL(il_rho(jp_J)-1,dp)*0.5) ENDIF !3-3 check domain validity CALL grid_check_dom(tl_bathy0, il_imin0, il_imax0, il_jmin0, il_jmax0) !3-4 check coincidence between coarse and fine grid CALL grid_check_coincidence( tl_bathy0, tl_bathy1, & & il_imin0, il_imax0, & & il_jmin0, il_jmax0, & & il_rho(:) ) !4- read or compute boundary tl_var=iom_read_var(tl_bathy1,'Bathymetry') tl_bdy(:)=boundary_init(tl_var, ln_north, ln_south, ln_east, ln_west, & & cn_north, cn_south, cn_east, cn_west, & & ln_oneseg ) !5- get boundary on coarse grid !5-1 define refined bathymetry table (for coarse grid) dl_fill=tl_var%d_fill ALLOCATE( dl_refined(tl_var%t_dim(1)%i_len, & & tl_var%t_dim(2)%i_len, & & tl_var%t_dim(3)%i_len, & & tl_var%t_dim(4)%i_len) ) dl_refined(:,:,:,:)=dl_fill !5-2 define weight table ALLOCATE( dl_weight(tl_var%t_dim(1)%i_len, & & tl_var%t_dim(2)%i_len, & & 1,1) ) dl_weight(:,:,:,:)=dl_fill !5-3 compute coarse grid refined bathymetry on boundary. DO jk=1,ip_ncard CALL merge_bathy_get_boundary(tl_bathy0, tl_bathy1, tl_bdy(jk), & & il_rho(:), & & dl_refined(:,:,:,:), dl_weight(:,:,:,:), & & dl_fill) ENDDO !6- merge bathy on boundary DO jl=1,tl_var%t_dim(4)%i_len DO jk=1,tl_var%t_dim(3)%i_len WHERE( dl_refined(:,:,jk,jl) /= dl_fill .AND. & & tl_var%d_value(:,:,jk,jl)/= tl_var%d_fill ) tl_var%d_value(:,:,jk,jl)= & & dl_weight(:,:,1,1) * dl_refined(:,:,jk,jl) + & & (1-dl_weight(:,:,1,1))*tl_var%d_value(:,:,jk,jl) ELSE WHERE( dl_refined(:,:,jk,jl)== dl_fill .AND. & & dl_weight(:,:,1,1) /= dl_fill ) ! to keep coarse grid mask on boundary tl_var%d_value(:,:,jk,jl)=dl_fill END WHERE ENDDO ENDDO DEALLOCATE(dl_refined) !7- create file tl_fileout=file_init(TRIM(cn_fileout),id_perio=in_perio1) !7-1 add dimension tl_dim(:)=tl_var%t_dim(:) DO ji=1,ip_maxdim IF( tl_dim(ji)%l_use ) CALL file_add_dim(tl_fileout, tl_dim(ji)) ENDDO !7-2 add variables IF( ALL( tl_dim(1:2)%l_use ) )THEN tl_lon=iom_read_var(tl_bathy1,'longitude') CALL file_add_var(tl_fileout, tl_lon) CALL var_clean(tl_lon) tl_lat=iom_read_var(tl_bathy1,'latitude') CALL file_add_var(tl_fileout, tl_lat) CALL var_clean(tl_lat) ENDIF CALL file_add_var(tl_fileout, tl_var) CALL var_clean(tl_var) ! only 2 first dimension to be used tl_dim(:)=tl_fileout%t_dim(:) tl_dim(3:4)%l_use=.FALSE. tl_var=var_init('weight',dl_weight(:,:,:,:),td_dim=tl_dim(:),dd_fill=dl_fill) CALL file_add_var(tl_fileout, tl_var) CALL var_clean(tl_var) !7-3 add some attribute tl_att=att_init("Created_by","SIREN merge_bathy") CALL file_add_att(tl_fileout, tl_att) cl_date=date_print(date_now()) tl_att=att_init("Creation_date",cl_date) CALL file_add_att(tl_fileout, tl_att) tl_att=att_init("coarse_grid_source_file",TRIM(fct_basename(tl_bathy0%c_name))) CALL file_add_att(tl_fileout, tl_att) tl_att=att_init("fine_grid_source_file",TRIM(fct_basename(tl_bathy1%c_name))) CALL file_add_att(tl_fileout, tl_att) ! a voir ! add attribute periodicity il_attid=0 IF( ASSOCIATED(tl_fileout%t_att) )THEN il_attid=att_get_id(tl_fileout%t_att(:),'periodicity') ENDIF IF( tl_bathy1%i_perio >= 0 .AND. il_attid == 0 )THEN tl_att=att_init('periodicity',tl_bathy1%i_perio) CALL file_add_att(tl_fileout,tl_att) ENDIF il_attid=0 IF( ASSOCIATED(tl_fileout%t_att) )THEN il_attid=att_get_id(tl_fileout%t_att(:),'ew_overlap') ENDIF IF( tl_bathy1%i_ew >= 0 .AND. il_attid == 0 )THEN tl_att=att_init('ew_overlap',tl_bathy1%i_ew) CALL file_add_att(tl_fileout,tl_att) ENDIF !7-4 create file CALL iom_create(tl_fileout) !7-5 write file CALL iom_write_file(tl_fileout) !7-6 close file CALL iom_close(tl_fileout) CALL iom_close(tl_bathy1) CALL iom_close(tl_bathy0) !8- clean CALL file_clean(tl_fileout) CALL file_clean(tl_bathy1) CALL file_clean(tl_bathy0) DEALLOCATE(dl_weight) ! close log file CALL logger_footer() CALL logger_close() !> @endcode CONTAINS !------------------------------------------------------------------- !> @brief !> This subroutine !> !> @details !> !> @author J.Paul !> - Nov, 2013- Initial Version !> !> @param[in] !> @todo !------------------------------------------------------------------- !> @code SUBROUTINE merge_bathy_get_boundary( td_bathy0, td_bathy1, td_bdy, & & id_rho, & & dd_refined, dd_weight, dd_fill ) IMPLICIT NONE ! Argument TYPE(TFILE) , INTENT(IN ) :: td_bathy0 TYPE(TFILE) , INTENT(IN ) :: td_bathy1 TYPE(TBDY) , INTENT(IN ) :: td_bdy INTEGER(i4), DIMENSION(:) , INTENT(IN ) :: id_rho REAL(dp) , DIMENSION(:,:,:,:), INTENT(INOUT) :: dd_refined REAL(dp) , DIMENSION(:,:,:,:), INTENT(INOUT) :: dd_weight REAL(dp) , INTENT(IN ) :: dd_fill ! local variable INTEGER(i4) :: il_imin1 INTEGER(i4) :: il_imax1 INTEGER(i4) :: il_jmin1 INTEGER(i4) :: il_jmax1 INTEGER(i4) :: il_imin0 INTEGER(i4) :: il_imax0 INTEGER(i4) :: il_jmin0 INTEGER(i4) :: il_jmax0 INTEGER(i4) :: il_imin INTEGER(i4) :: il_imax INTEGER(i4) :: il_jmin INTEGER(i4) :: il_jmax INTEGER(i4), DIMENSION(2,2) :: il_offset INTEGER(i4), DIMENSION(2,2,2) :: il_ind REAL(dp) , DIMENSION(:) , ALLOCATABLE :: dl_tmp1d REAL(dp) , DIMENSION(:,:) , ALLOCATABLE :: dl_tmp2d REAL(dp) , DIMENSION(:,:) , ALLOCATABLE :: dl_wseg TYPE(TVAR) :: tl_var0 TYPE(TVAR) :: tl_lon1 TYPE(TVAR) :: tl_lat1 TYPE(TFILE) :: tl_bathy1 TYPE(TFILE) :: tl_bathy0 TYPE(TMPP) :: tl_mppbathy1 TYPE(TMPP) :: tl_mppbathy0 TYPE(TDOM) :: tl_dom1 TYPE(TDOM) :: tl_dom0 ! loop indices INTEGER(i4) :: ji INTEGER(i4) :: jl !---------------------------------------------------------------- IF( td_bdy%l_use )THEN DO jl=1,td_bdy%i_nseg !1- get boundary definition SELECT CASE(TRIM(td_bdy%c_card)) CASE('north') il_imin1=td_bdy%t_seg(jl)%i_first il_imax1=td_bdy%t_seg(jl)%i_last il_jmin1=td_bdy%t_seg(jl)%i_index-(td_bdy%t_seg(jl)%i_width-1) il_jmax1=td_bdy%t_seg(jl)%i_index il_imin=1 il_imax=il_imax1-il_imin1+1 il_jmin=td_bdy%t_seg(jl)%i_width il_jmax=1 CASE('south') il_imin1=td_bdy%t_seg(jl)%i_first il_imax1=td_bdy%t_seg(jl)%i_last il_jmin1=td_bdy%t_seg(jl)%i_index il_jmax1=td_bdy%t_seg(jl)%i_index+(td_bdy%t_seg(jl)%i_width-1) il_imin=1 il_imax=il_imax1-il_imin1+1 il_jmin=1 il_jmax=td_bdy%t_seg(jl)%i_width CASE('east') il_imin1=td_bdy%t_seg(jl)%i_index-(td_bdy%t_seg(jl)%i_width-1) il_imax1=td_bdy%t_seg(jl)%i_index il_jmin1=td_bdy%t_seg(jl)%i_first il_jmax1=td_bdy%t_seg(jl)%i_last il_imin=td_bdy%t_seg(jl)%i_width il_imax=1 il_jmin=1 il_jmax=il_jmax1-il_jmin1+1 CASE('west') il_imin1=td_bdy%t_seg(jl)%i_index il_imax1=td_bdy%t_seg(jl)%i_index+(td_bdy%t_seg(jl)%i_width-1) il_jmin1=td_bdy%t_seg(jl)%i_first il_jmax1=td_bdy%t_seg(jl)%i_last il_imin=1 il_imax=td_bdy%t_seg(jl)%i_width il_jmin=1 il_jmax=il_jmax1-il_jmin1+1 END SELECT !2 -read fine grid domain tl_bathy1=td_bathy1 CALL iom_open(tl_bathy1) !2-1 compute domain tl_dom1=dom_init( tl_bathy1, & & il_imin1, il_imax1,& & il_jmin1, il_jmax1 ) !2-2 close file CALL iom_close(tl_bathy1) !2-3 read variables on domain (ugly way to do it, have to work on it) !2-3-1 init mpp structure tl_mppbathy1=mpp_init(tl_bathy1) CALL file_clean(tl_bathy1) !2-3-2 get processor to be used CALL mpp_get_use( tl_mppbathy1, tl_dom1 ) !2-3-3 open mpp files CALL iom_mpp_open(tl_mppbathy1) !2-3-4 read variable value on domain tl_lon1=iom_mpp_read_var(tl_mppbathy1,'longitude',td_dom=tl_dom1) tl_lat1=iom_mpp_read_var(tl_mppbathy1,'latitude' ,td_dom=tl_dom1) !2-3-5 close mpp files CALL iom_mpp_close(tl_mppbathy1) !2-3-6 clean structure CALL mpp_clean(tl_mppbathy1) !3- get coarse grid indices il_ind(:,:,:)=grid_get_coarse_index(td_bathy0, tl_lon1, tl_lat1, & & id_rho=id_rho(:)) il_imin0=il_ind(1,1,1) il_imax0=il_ind(1,2,1) il_jmin0=il_ind(2,1,1) il_jmax0=il_ind(2,2,1) il_offset(:,:)=il_ind(:,:,2) !4- read coarse grid bathymetry on domain tl_bathy0=td_bathy0 CALL iom_open(tl_bathy0) !4-1 compute domain tl_dom0=dom_init( tl_bathy0, & & il_imin0, il_imax0,& & il_jmin0, il_jmax0 ) !4-2 close file CALL iom_close(tl_bathy0) !4-3 add extra band (if possible) to compute interpolation CALL dom_add_extra(tl_dom0) !4-4 read variables on domain (ugly way to do it, have to work on it) !4-4-1 init mpp structure tl_mppbathy0=mpp_init(tl_bathy0) CALL file_clean(tl_bathy0) !4-4-2 get processor to be used CALL mpp_get_use( tl_mppbathy0, tl_dom0 ) !4-4-3 open mpp files CALL iom_mpp_open(tl_mppbathy0) !4-4-4 read variable value on domain tl_var0=iom_mpp_read_var(tl_mppbathy0,'Bathymetry',td_dom=tl_dom0) !4-4-5 close mpp files CALL iom_mpp_close(tl_mppbathy0) !4-4-6 clean structure CALL mpp_clean(tl_mppbathy0) !5- interpolate variable CALL merge_bathy_interp( tl_var0, & & id_rho(:), & & il_offset(:,:) ) !6- remove extraband added to domain CALL dom_del_extra( tl_var0, tl_dom0, id_rho(:) ) !6-1 remove extraband added to domain CALL dom_clean_extra( tl_dom0 ) !7- fill refined table !7-1 keep only useful point ! interpolation could create more point than necessary CALL boundary_clean_interp(tl_var0, td_bdy ) ! use add request ???? !7-2 fill refined table dd_refined( il_imin1:il_imax1, & & il_jmin1:il_jmax1, & & :,: )= tl_var0%d_value(:,:,:,:) !8- compute weight function ALLOCATE( dl_tmp1d(td_bdy%t_seg(jl)%i_width) ) SELECT CASE(TRIM(td_bdy%c_card)) CASE('north') ! compute "distance" dl_tmp1d(:)=(/(ji-1,ji=td_bdy%t_seg(jl)%i_width,1,-1)/) ! compute weight on segment dl_tmp1d(:)= 0.5 + 0.5*COS( (dg_pi*dl_tmp1d(:)) / & & (td_bdy%t_seg(jl)%i_width) ) ALLOCATE( dl_wseg(tl_dom1%t_dim(1)%i_len, & & tl_dom1%t_dim(2)%i_len) ) dl_wseg(:,:)=dd_fill dl_wseg(:,:)=SPREAD( dl_tmp1d(:), & & DIM=1, & & NCOPIES=tl_dom1%t_dim(1)%i_len ) CASE('south') ! compute "distance" dl_tmp1d(:)=(/(ji-1,ji=1,td_bdy%t_seg(jl)%i_width)/) ! compute weight on segment dl_tmp1d(:)= 0.5 + 0.5*COS( (dg_pi*dl_tmp1d(:)) / & & (td_bdy%t_seg(jl)%i_width) ) ALLOCATE( dl_wseg(tl_dom1%t_dim(1)%i_len, & & tl_dom1%t_dim(2)%i_len) ) dl_wseg(:,:)=dd_fill dl_wseg(:,:)=SPREAD( dl_tmp1d(:), & & DIM=1, & & NCOPIES=tl_dom1%t_dim(1)%i_len ) CASE('east') ! compute "distance" dl_tmp1d(:)=(/(ji-1,ji=td_bdy%t_seg(jl)%i_width,1,-1)/) ! compute weight on segment dl_tmp1d(:)= 0.5 + 0.5*COS( (dg_pi*dl_tmp1d(:)) / & & (td_bdy%t_seg(jl)%i_width) ) ALLOCATE( dl_wseg(tl_dom1%t_dim(1)%i_len, & & tl_dom1%t_dim(2)%i_len) ) dl_wseg(:,:)=dd_fill dl_wseg(:,:)=SPREAD( dl_tmp1d(:), & & DIM=2, & & NCOPIES=tl_dom1%t_dim(2)%i_len ) CASE('west') ! compute "distance" dl_tmp1d(:)=(/(ji-1,ji=1,td_bdy%t_seg(jl)%i_width)/) ! compute weight on segment dl_tmp1d(:)= 0.5 + 0.5*COS( (dg_pi*dl_tmp1d(:)) / & & (td_bdy%t_seg(jl)%i_width) ) ALLOCATE( dl_wseg(tl_dom1%t_dim(1)%i_len, & & tl_dom1%t_dim(2)%i_len) ) dl_wseg(:,:)=dd_fill dl_wseg(:,:)=SPREAD( dl_tmp1d(:), & & DIM=2, & & NCOPIES=tl_dom1%t_dim(2)%i_len ) END SELECT DEALLOCATE( dl_tmp1d ) !8-1 fill weight table ALLOCATE( dl_tmp2d( tl_dom1%t_dim(1)%i_len, & & tl_dom1%t_dim(2)%i_len) ) dl_tmp2d(:,:)=dd_weight( il_imin1:il_imax1, & & il_jmin1:il_jmax1, & & 1,1 ) WHERE( dl_tmp2d(:,:) == dd_fill ) dl_tmp2d(:,:)= dl_wseg(:,:) ELSE WHERE dl_tmp2d(:,:)= MAX( dl_tmp2d(:,:), dl_wseg(:,:) ) END WHERE dd_weight( il_imin1:il_imax1, & & il_jmin1:il_jmax1, & & 1,1 ) = dl_tmp2d(:,:) DEALLOCATE( dl_wseg ) DEALLOCATE( dl_tmp2d ) ENDDO ENDIF END SUBROUTINE merge_bathy_get_boundary !> @endcode !------------------------------------------------------------------- !> @brief !> This subroutine !> !> @details !> !> @author J.Paul !> - Nov, 2013- Initial Version !> !> @param[in] !> @todo !------------------------------------------------------------------- !> @code SUBROUTINE merge_bathy_interp( td_var, & & id_rho, & & id_offset, & & id_iext, id_jext) IMPLICIT NONE ! Argument TYPE(TVAR) , INTENT(INOUT) :: td_var INTEGER(i4), DIMENSION(:) , INTENT(IN ) :: id_rho INTEGER(i4), DIMENSION(:,:), INTENT(IN ) :: id_offset INTEGER(i4), INTENT(IN ), OPTIONAL :: id_iext INTEGER(i4), INTENT(IN ), OPTIONAL :: id_jext ! local variable TYPE(TVAR) :: tl_var TYPE(TVAR) :: tl_mask INTEGER(i1), DIMENSION(:,:,:,:), ALLOCATABLE :: bl_mask INTEGER(i4) :: il_iext INTEGER(i4) :: il_jext ! loop indices !---------------------------------------------------------------- ! copy variable tl_var=td_var !WARNING: two extrabands are required for cubic interpolation il_iext=3 IF( PRESENT(id_iext) ) il_iext=id_iext il_jext=3 IF( PRESENT(id_jext) ) il_jext=id_jext IF( il_iext < 2 .AND. td_var%c_interp(1) == 'cubic' )THEN CALL logger_warn("CREATE BATHY INTERP: at least extrapolation "//& & "on two points are required with cubic interpolation ") il_iext=2 ENDIF IF( il_jext < 2 .AND. td_var%c_interp(1) == 'cubic' )THEN CALL logger_warn("CREATE BATHY INTERP: at least extrapolation "//& & "on two points are required with cubic interpolation ") il_jext=2 ENDIF !1- work on mask !1-1 create mask ALLOCATE(bl_mask(tl_var%t_dim(1)%i_len, & & tl_var%t_dim(2)%i_len, & & tl_var%t_dim(3)%i_len, & & tl_var%t_dim(4)%i_len) ) bl_mask(:,:,:,:)=1 WHERE(tl_var%d_value(:,:,:,:)==tl_var%d_fill) bl_mask(:,:,:,:)=0 SELECT CASE(TRIM(tl_var%c_point)) CASE DEFAULT ! 'T' tl_mask=var_init('tmask', bl_mask(:,:,:,:), td_dim=tl_var%t_dim(:)) CASE('U') tl_mask=var_init('umask', bl_mask(:,:,:,:), td_dim=tl_var%t_dim(:)) CASE('V') tl_mask=var_init('vmask', bl_mask(:,:,:,:), td_dim=tl_var%t_dim(:)) CASE('F') tl_mask=var_init('fmask', bl_mask(:,:,:,:), td_dim=tl_var%t_dim(:)) END SELECT DEALLOCATE(bl_mask) !1-2 interpolate mask CALL interp_fill_value( tl_mask, id_rho(:), & & id_offset=id_offset(:,:) ) !2- work on variable !2-0 add extraband CALL extrap_add_extrabands(tl_var, il_iext, il_iext) !2-1 extrapolate variable CALL extrap_fill_value( tl_var, id_offset=id_offset(:,:), & & id_rho=id_rho(:), & & id_iext=il_iext, id_jext=il_jext ) !2-2 interpolate Bathymetry CALL interp_fill_value( tl_var, id_rho(:), & & id_offset=id_offset(:,:) ) !2-3 remove extraband CALL extrap_del_extrabands(tl_var, il_iext*id_rho(jp_I), il_jext*id_rho(jp_J)) !2-4 keep original mask WHERE( tl_mask%d_value(:,:,:,:) == 0 ) tl_var%d_value(:,:,:,:)=tl_var%d_fill END WHERE !3- save result td_var=tl_var ! clean variable structure CALL var_clean(tl_var) END SUBROUTINE merge_bathy_interp !> @endcode END PROGRAM merge_bathy