!---------------------------------------------------------------------- ! NEMO system team, System and Interface for oceanic RElocable Nesting !---------------------------------------------------------------------- ! ! ! PROGRAM: create_boundary ! ! DESCRIPTION: !> @brief !> This program create boundary files. !> !> @details !> Variables are read from standard output. !> Then theses variables are interpolated on fine grid boundaries. !> !> @author !> J.Paul ! REVISION HISTORY: !> @date Nov, 2013 - Initial Version ! !> @note Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) !> !> @todo !---------------------------------------------------------------------- !> @code PROGRAM create_boundary 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 multi ! multi file manager USE boundary ! boundary manager USE iom ! I/O manager USE dom ! domain manager USE grid ! grid manager USE vgrid ! vartical 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 CHARACTER(LEN=lc) :: cl_name CHARACTER(LEN=lc) :: cl_bdyout CHARACTER(LEN=lc) :: cl_data INTEGER(i4) :: il_narg INTEGER(i4) :: il_status INTEGER(i4) :: il_fileid INTEGER(i4) :: il_attid INTEGER(i4) :: il_dim 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 LOGICAL :: ll_exist TYPE(TFILE) :: tl_coord0 TYPE(TFILE) :: tl_bathy0 TYPE(TFILE) :: tl_coord1 TYPE(TFILE) :: tl_bathy1 TYPE(TFILE) :: tl_file TYPE(TFILE) :: tl_fileout TYPE(TMPP) :: tl_mpp TYPE(TMULTI) :: tl_multi TYPE(TATT) :: tl_att TYPE(TVAR) , DIMENSION(:) , ALLOCATABLE :: tl_level TYPE(TVAR) , DIMENSION(:,:,:) , ALLOCATABLE :: tl_seglvl1 TYPE(TVAR) :: tl_var1 TYPE(TVAR) , DIMENSION(:,:,:) , ALLOCATABLE :: tl_segvar1 TYPE(TVAR) , DIMENSION(:,:) , ALLOCATABLE :: tl_seglon1 TYPE(TVAR) , DIMENSION(:,:) , ALLOCATABLE :: tl_seglat1 TYPE(TVAR) , DIMENSION(:,:) , ALLOCATABLE :: tl_var TYPE(TVAR) , DIMENSION(:) , ALLOCATABLE :: tl_lon1 TYPE(TVAR) , DIMENSION(:) , ALLOCATABLE :: tl_lat1 TYPE(TVAR) :: tl_depth TYPE(TVAR) :: tl_time TYPE(TVAR) :: tl_tmp TYPE(TDIM) , DIMENSION(ip_maxdim) :: tl_dim TYPE(TBDY) , DIMENSION(ip_ncard) :: tl_bdy TYPE(TDOM) :: tl_dom0 TYPE(TDOM) , DIMENSION(:,:) , ALLOCATABLE :: tl_segdom1 ! loop indices INTEGER(i4) :: jvar INTEGER(i4) :: ji INTEGER(i4) :: jj INTEGER(i4) :: jk INTEGER(i4) :: jl ! namelist variable ! namlog CHARACTER(LEN=lc) :: cn_logfile = 'create_boundary.log' CHARACTER(LEN=lc) :: cn_verbosity = 'warning' ! namcrs CHARACTER(LEN=lc) :: cn_coord0 = '' INTEGER(i4) :: in_perio0 = -1 ! namfin CHARACTER(LEN=lc) :: cn_coord1 = '' CHARACTER(LEN=lc) :: cn_bathy1 = '' INTEGER(i4) :: in_perio1 = -1 ! namcfg CHARACTER(LEN=lc) :: cn_varcfg = 'variable.cfg' ! namvar CHARACTER(LEN=lc), DIMENSION(ig_maxvar) :: cn_varinfo = '' CHARACTER(LEN=lc), DIMENSION(ig_maxvar) :: cn_varfile = '' ! namnst INTEGER(i4) :: in_imin0 = 0 INTEGER(i4) :: in_imax0 = 0 INTEGER(i4) :: in_jmin0 = 0 INTEGER(i4) :: in_jmax0 = 0 INTEGER(i4) :: in_rhoi = 1 INTEGER(i4) :: in_rhoj = 1 ! nambdy LOGICAL :: ln_north = .TRUE. LOGICAL :: ln_south = .TRUE. LOGICAL :: ln_east = .TRUE. LOGICAL :: ln_west = .TRUE. CHARACTER(LEN=lc) :: cn_north = '' CHARACTER(LEN=lc) :: cn_south = '' CHARACTER(LEN=lc) :: cn_east = '' CHARACTER(LEN=lc) :: cn_west = '' LOGICAL :: ln_oneseg = .TRUE. INTEGER(i4) :: in_extrap = 0 ! namout CHARACTER(LEN=lc) :: cn_fileout = 'boundary.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_coord0, & !< coordinate file & in_perio0 !< periodicity index NAMELIST /namfin/ & !< fine grid namelist & cn_coord1, & !< coordinate file & cn_bathy1, & !< bathymetry file & in_perio1 !< periodicity index NAMELIST /namvar/ & !< variable namelist & cn_varinfo, & !< list of variable and method to apply on. (ex: 'votemper:linear','vosaline:cubic' ) & cn_varfile !< list of variable and file where find it. (ex: 'votemper:GLORYS_gridT.nc' ) 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 & in_extrap !< number of mask point to extrapolate NAMELIST /namout/ & !< output namelist & cn_fileout !< fine grid boundary file basename !------------------------------------------------------------------- !1- namelist !1-1 get namelist il_narg=COMMAND_ARGUMENT_COUNT() !f03 intrinsec IF( il_narg/=1 )THEN PRINT *,"CREATE BOUNDARY: 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 *,"CREATE BOUNDARY: 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) !1-2-4 match variable with file tl_multi=multi_init(cn_varfile) 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("CREATE BOUNDARY: ERROR closing "//TRIM(cl_namelist)) ENDIF ELSE PRINT *,"CREATE BOUNDARY: ERROR. can not find "//TRIM(cl_namelist) ENDIF !2- open files IF( TRIM(cn_coord0) /= '' )THEN tl_coord0=file_init(TRIM(cn_coord0),id_perio=in_perio0) CALL iom_open(tl_coord0) ELSE CALL logger_fatal("CREATE BOUNDARY: can not find coarse grid "//& & "coordinate file. check namelist") ENDIF IF( TRIM(cn_coord1) /= '' )THEN tl_coord1=file_init(TRIM(cn_coord1),id_perio=in_perio1) CALL iom_open(tl_coord1) ELSE CALL logger_fatal("CREATE BOUNDARY: can not find fine grid coordinate "//& & "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("CREATE BOUNDARY: can not find fine grid bathymetry "//& & "file. check namelist") ENDIF !3- check !3-1 check output file do not already exist DO jk=1,ip_ncard cl_bdyout=boundary_set_filename( TRIM(cn_fileout), & & TRIM(ip_card(jk)) ) INQUIRE(FILE=TRIM(cl_bdyout), EXIST=ll_exist) IF( ll_exist )THEN CALL logger_fatal("CREATE BOUNDARY: output file "//TRIM(cl_bdyout)//& & " already exist.") ENDIF ENDDO !3-1 check namelist !3-1-1 check refinement factor il_rho(:)=1 IF( in_rhoi < 1 .OR. in_rhoj < 1 )THEN CALL logger_error("CREATE BOUNDARY: invalid refinement factor."//& & " check namelist "//TRIM(cl_namelist)) ELSE il_rho(jp_I)=in_rhoi il_rho(jp_J)=in_rhoj ENDIF !3-1-2 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) ELSE il_imin0=in_imin0 ; il_imax0=in_imax0 il_jmin0=in_jmin0 ; il_jmax0=in_jmax0 ENDIF !3-2 check domain validity CALL grid_check_dom(tl_coord0, il_imin0, il_imax0, il_jmin0, il_jmax0) !3-3 check coordinate file CALL grid_check_coincidence( tl_coord0, tl_coord1, & & il_imin0, il_imax0, & & il_jmin0, il_jmax0, & & il_rho(:) ) !4- read or compute boundary tl_var1=iom_read_var(tl_bathy1,'Bathymetry') tl_bdy(:)=boundary_init(tl_var1, ln_north, ln_south, ln_east, ln_west, & & cn_north, cn_south, cn_east, cn_west, & & ln_oneseg ) CALL var_clean(tl_var1) !5- compute level ALLOCATE(tl_level(ig_npoint)) tl_level(:)=vgrid_get_level(tl_bathy1, cl_namelist ) !6- get coordinate on each segment of each boundary ALLOCATE( tl_seglon1(ip_ncard,ig_maxseg) ) ALLOCATE( tl_seglat1(ip_ncard,ig_maxseg) ) ALLOCATE( tl_segdom1(ip_ncard,ig_maxseg) ) ALLOCATE( tl_segvar1(tl_multi%i_nvar,ip_ncard,ig_maxseg) ) ALLOCATE( tl_seglvl1(ip_ncard,ig_maxseg,ig_npoint) ) DO jk=1,ip_ncard IF( tl_bdy(jk)%l_use )THEN DO jl=1,tl_bdy(jk)%i_nseg !6-1 get fine grid segment domain tl_segdom1(jk,jl)=create_boundary_get_dom( tl_bathy1, tl_bdy(jk), jl ) !6-2 get fine grid segment coordinate CALL create_boundary_get_coord( tl_bathy1, tl_segdom1(jk,jl), & & tl_seglon1(jk,jl), tl_seglat1(jk,jl) ) !6-2 get fine grid segment coordinate tl_seglvl1(jk,jl,:)=create_bdy_get_level(tl_level(:), tl_segdom1(jk,jl)) ENDDO ENDIF ENDDO DEALLOCATE(tl_level) !7- compute boundary for variable to be used (see namelist) IF( .NOT. ASSOCIATED(tl_multi%t_file) )THEN CALL logger_error("CREATE BOUNDARY: no file to work on. "//& & "check cn_varfile in namelist.") ELSE jvar=0 ! for each file DO ji=1,tl_multi%i_nfile WRITE(cl_data,'(a,i2.2)') 'data_',jvar+1 IF( .NOT. ASSOCIATED(tl_multi%t_file(ji)%t_var) )THEN CALL logger_error("CREATE BOUNDARY: no variable to work on for "//& & "file "//TRIM(tl_multi%t_file(ji)%c_name)//& & ". check cn_varfile in namelist.") ELSE IF( .NOT. ASSOCIATED(tl_multi%t_file(ji)%t_var) )THEN CALL logger_error("CREATE BOUNDARY: no variable to work on for "//& & "file "//TRIM(tl_multi%t_file(ji)%c_name)//& & ". check cn_varfile in namelist.") ELSEIF( TRIM(tl_multi%t_file(ji)%c_name) == TRIM(cl_data) )THEN !- use input matrix to fill variable ! for each variable initialise from matrix DO jj=1,tl_multi%t_file(ji)%i_nvar jvar=jvar+1 tl_tmp=tl_multi%t_file(ji)%t_var(jj) DO jk=1,ip_ncard IF( tl_bdy(jk)%l_use )THEN DO jl=1,tl_bdy(jk)%i_nseg !7-1 fill value with matrix data ! pb voir comment gerer nb de dimension tl_segvar1(jvar,jk,jl)=create_bdy_matrix(tl_tmp, tl_segdom1(jk,jl), tl_coord1) !7-2 use mask CALL create_bdy_use_mask(tl_segvar1(jvar,jk,jl), tl_seglvl1(jk,jl,:)) ENDDO ENDIF ENDDO ENDDO ELSE !- use file to fill variable ! open file tl_file=file_init(TRIM(tl_multi%t_file(ji)%c_name)) CALL iom_open(tl_file) ! get or check depth value IF( tl_file%i_depthid /= 0 )THEN IF( ASSOCIATED(tl_depth%d_value) )THEN IF( ANY( tl_depth%d_value(:,:,:,:) /= & & tl_tmp%d_value(:,:,:,:) ) )THEN CALL logger_fatal("CREATE BOUNDARY: depth value from "//& & TRIM(tl_multi%t_file(ji)%c_name)//" not conform "//& & " to those from former file(s).") ENDIF ELSE tl_depth=iom_read_var(tl_file,tl_file%i_depthid) ENDIF ENDIF ! get or check time value IF( tl_file%i_timeid /= 0 )THEN IF( ASSOCIATED(tl_time%d_value) )THEN IF( ANY( tl_time%d_value(:,:,:,:) /= & & tl_tmp%d_value(:,:,:,:) ) )THEN CALL logger_fatal("CREATE BOUNDARY: time value from "//& & TRIM(tl_multi%t_file(ji)%c_name)//" not conform "//& & " to those from former file(s).") ENDIF ELSE tl_time=iom_read_var(tl_file,tl_file%i_timeid) ENDIF ENDIF IF( ANY( tl_file%t_dim(1:2)%i_len /= & & tl_coord0%t_dim(1:2)%i_len) )THEN !- extract value from fine grid DO jk=1,ip_ncard IF( tl_bdy(jk)%l_use )THEN DO jl=1,tl_bdy(jk)%i_nseg !7-1 compute domain on fine grid ! open mpp file on domain !7-2 init mpp structure tl_mpp=mpp_init(tl_file) !7-3 get processor to be used CALL mpp_get_use( tl_mpp, tl_segdom1(jk,jl) ) !7-4 open mpp files CALL iom_mpp_open(tl_mpp) ! for each variable of this file DO jj=1,tl_multi%t_file(ji)%i_nvar cl_name=tl_multi%t_file(ji)%t_var(jj)%c_name !7-5 read variable over domain tl_segvar1(jvar+jj,jk,jl)=iom_mpp_read_var( tl_mpp, TRIM(cl_name), & & td_dom=tl_segdom1(jk,jl) ) !7-6 add attribute to variable tl_att=att_init('src_file',TRIM(fct_basename(tl_mpp%c_name))) CALL var_move_att(tl_segvar1(jvar+jj,jk,jl), tl_att) tl_att=att_init('src_i-indices',(/tl_segdom1(jk,jl)%i_imin, tl_segdom1(jk,jl)%i_imax/)) CALL var_move_att(tl_segvar1(jvar+jj,jk,jl), tl_att) tl_att=att_init('src_j-indices',(/tl_segdom1(jk,jl)%i_jmin, tl_segdom1(jk,jl)%i_jmax/)) CALL var_move_att(tl_segvar1(jvar+jj,jk,jl), tl_att) ! clean structure CALL att_clean(tl_att) !7-7 use mask CALL create_bdy_use_mask(tl_segvar1(jvar+jj,jk,jl), tl_seglvl1(jk,jl,:)) ENDDO !7-8 close mpp files CALL iom_mpp_close(tl_mpp) CALL mpp_clean(tl_mpp) ENDDO ENDIF ENDDO jvar=jvar+tl_multi%t_file(ji)%i_nvar ELSE !- interpolate value from coarse grid DO jk=1,ip_ncard IF( tl_bdy(jk)%l_use )THEN DO jl=1,tl_bdy(jk)%i_nseg !7-1 get coarse grid indices of this segment il_ind(:,:,:)=grid_get_coarse_index(tl_coord0, & & tl_seglon1(jk,jl), tl_seglat1(jk,jl), & & id_rho=il_rho(:) ) IF( ANY(il_ind(:,:,:)==0) )THEN CALL logger_error("CREATE BOUNDARY: error computing "//& & " coarse grid indices") ELSE 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) ENDIF !7-2 compute coarse grid segment domain tl_dom0=dom_init( tl_coord0, & & il_imin0, il_imax0,& & il_jmin0, il_jmax0 ) !7-3 add extra band (if possible) to compute interpolation CALL dom_add_extra(tl_dom0) !7-4 read variables on domain (ugly way to do it, have to work on it) !7-4-1 init mpp structure tl_mpp=mpp_init(tl_file) !7-4-2 get processor to be used CALL mpp_get_use( tl_mpp, tl_dom0 ) !7-4-3 open mpp files CALL iom_mpp_open(tl_mpp) ! check file dimension IF( ANY(tl_mpp%t_dim(1:2)%i_len /= tl_coord0%t_dim(1:2)%i_len) )THEN CALL logger_error("CREATE BOUNDARY INTERP: dimension of file "//& & TRIM(tl_mpp%c_name)//" not conform to those of "//& & TRIM(tl_coord0%c_name)) ELSE ! for each variable of this file DO jj=1,tl_multi%t_file(ji)%i_nvar cl_name=tl_multi%t_file(ji)%t_var(jj)%c_name !7-4-4 read variable value on domain tl_segvar1(jvar+jj,jk,jl)=iom_mpp_read_var( tl_mpp, TRIM(cl_name), & & td_dom=tl_dom0 ) !7-4-5 work on variable CALL create_boundary_interp(tl_segvar1(jvar+jj,jk,jl), & & il_rho(:), & & il_offset(:,:) ) !7-4-6 remove extraband added to domain CALL dom_del_extra( tl_segvar1(jvar+jj,jk,jl), tl_dom0, il_rho(:) ) !7-4-7 keep only useful point (width) ! interpolation could create more point than necessary CALL boundary_clean_interp(tl_segvar1(jvar+jj,jk,jl), tl_bdy(jk) ) !7-4-8 add attribute to variable tl_att=att_init('src_file',TRIM(fct_basename(tl_mpp%c_name))) CALL var_move_att(tl_segvar1(jvar+jj,jk,jl), tl_att) tl_att=att_init('src_i-indices',(/tl_dom0%i_imin, tl_dom0%i_imax/)) CALL var_move_att(tl_segvar1(jvar+jj,jk,jl), tl_att) tl_att=att_init('src_j-indices',(/tl_dom0%i_jmin, tl_dom0%i_jmax/)) CALL var_move_att(tl_segvar1(jvar+jj,jk,jl), tl_att) ! clean structure CALL att_clean(tl_att) !7-4-9 use mask CALL create_bdy_use_mask(tl_segvar1(jvar+jj,jk,jl), tl_seglvl1(jk,jl,:)) ENDDO ENDIF CALL dom_clean(tl_dom0) !7-5 close mpp files CALL iom_mpp_close(tl_mpp) !7-6 clean structure CALL mpp_clean(tl_mpp) ENDDO ENDIF ENDDO jvar=jvar+tl_multi%t_file(ji)%i_nvar ENDIF CALL file_clean(tl_file) ENDIF ENDIF ENDDO ENDIF IF( jvar /= tl_multi%i_nvar )THEN CALL logger_error("CREATE BOUNDARY: it seems some variable can not be read") ENDIF !8- concatenate file ALLOCATE( tl_lon1(ip_ncard) ) ALLOCATE( tl_lat1(ip_ncard) ) ALLOCATE( tl_var(tl_multi%i_nvar,ip_ncard) ) DO jk=1,ip_ncard IF( tl_bdy(jk)%l_use )THEN SELECT CASE(TRIM(tl_bdy(jk)%c_card)) CASE('north','south') il_dim=1 CASE('east','west') il_dim=2 END SELECT DO jl=1,tl_bdy(jk)%i_nseg !- concatenate variable IF( jl == 1 )THEN tl_lon1(jk)=tl_seglon1(jk,jl) tl_lat1(jk)=tl_seglat1(jk,jl) DO jvar=1,tl_multi%i_nvar tl_var(jvar,jk)=tl_segvar1(jvar,jk,jl) ENDDO ELSE tl_lon1(jk)=var_concat(tl_lon1(jk),tl_seglon1(jk,jl),DIM=il_dim) tl_lat1(jk)=var_concat(tl_lat1(jk),tl_seglat1(jk,jl),DIM=il_dim) DO jvar=1,tl_multi%i_nvar tl_var(jvar,jk)=var_concat(tl_var(jvar,jk),tl_segvar1(jvar,jk,jl),DIM=il_dim) ENDDO ENDIF ENDDO ! swap array CALL boundary_swap(tl_lon1(jk), tl_bdy(jk)) CALL boundary_swap(tl_lat1(jk), tl_bdy(jk)) DO jvar=1,tl_multi%i_nvar CALL boundary_swap(tl_var(jvar,jk), tl_bdy(jk)) !9- use additional request !9-1 forced min and max value CALL var_limit_value(tl_var(jvar,jk)) !9-2 filter CALL filter_fill_value(tl_var(jvar,jk)) !9-3 extrapolate CALL extrap_fill_value(tl_var(jvar,jk), id_iext=in_extrap, & & id_jext=in_extrap, & & id_kext=in_extrap) ENDDO ENDIF ENDDO DEALLOCATE( tl_seglon1 ) DEALLOCATE( tl_seglat1 ) DEALLOCATE( tl_segdom1 ) DEALLOCATE( tl_segvar1 ) DEALLOCATE( tl_seglvl1 ) DO jk=1,ip_ncard IF( tl_bdy(jk)%l_use )THEN !10 create file !10-1 create file structure !10-1-1 set file name cl_bdyout=boundary_set_filename( TRIM(cn_fileout), & & TRIM(tl_bdy(jk)%c_card) ) !10-1-2 tl_fileout=file_init(TRIM(cl_bdyout),id_perio=in_perio1) !10-2 add dimension tl_dim(:)=var_max_dim(tl_var(:,jk)) DO ji=1,ip_maxdim IF( tl_dim(ji)%l_use ) CALL file_add_dim(tl_fileout, tl_dim(ji)) ENDDO !10-3 add variables IF( ALL( tl_dim(1:2)%l_use ) )THEN ! add longitude CALL file_add_var(tl_fileout, tl_lon1(jk)) CALL var_clean(tl_lon1(jk)) ! add latitude CALL file_add_var(tl_fileout, tl_lat1(jk)) CALL var_clean(tl_lat1(jk)) ENDIF IF( tl_dim(3)%l_use )THEN ! add depth CALL file_add_var(tl_fileout, tl_depth) ENDIF IF( tl_dim(4)%l_use )THEN ! add time CALL file_add_var(tl_fileout, tl_time) ENDIF ! add other variable DO jvar=1,tl_multi%i_nvar !IF( TRIM(tl_var(jvar,jk)%c_name) /= 'X' .AND. & !& TRIM(tl_var(jvar,jk)%c_name) /= 'Y' )THEN CALL file_add_var(tl_fileout, tl_var(jvar,jk)) !ENDIF CALL var_clean(tl_var(jvar,jk)) ENDDO !10-4 add some attribute tl_att=att_init("Created_by","SIREN create_boundary") 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) ! 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_coord1%i_perio >= 0 .AND. il_attid == 0 )THEN tl_att=att_init('periodicity',tl_coord1%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_coord1%i_ew >= 0 .AND. il_attid == 0 )THEN tl_att=att_init('ew_overlap',tl_coord1%i_ew) CALL file_add_att(tl_fileout,tl_att) ENDIF !10-5 create file CALL iom_create(tl_fileout) !10-6 write file CALL iom_write_file(tl_fileout) !10-7 close file CALL iom_close(tl_fileout) CALL file_clean(tl_fileout) ENDIF ENDDO DEALLOCATE( tl_lon1 ) DEALLOCATE( tl_lat1 ) DEALLOCATE( tl_var ) !11- close file CALL iom_close(tl_bathy1) CALL iom_close(tl_coord1) CALL iom_close(tl_coord0) !12- clean CALL var_clean(tl_depth) CALL var_clean(tl_time) CALL file_clean(tl_fileout) CALL file_clean(tl_bathy1) CALL file_clean(tl_coord1) CALL file_clean(tl_coord0) ! close log file CALL logger_footer() CALL logger_close() !> @endcode CONTAINS !------------------------------------------------------------------- !> @brief !> This subroutine !> !> @details !> !> @author J.Paul !> - 2013- Initial Version !> !> @param[in] !> @todo !------------------------------------------------------------------- !> @code FUNCTION create_boundary_get_dom( td_bathy1, td_bdy, id_seg ) IMPLICIT NONE ! Argument TYPE(TFILE), INTENT(IN ) :: td_bathy1 TYPE(TBDY) , INTENT(IN ) :: td_bdy INTEGER(i4), INTENT(IN ) :: id_seg ! function TYPE(TDOM) :: create_boundary_get_dom ! local variable INTEGER(i4) :: il_imin1 INTEGER(i4) :: il_imax1 INTEGER(i4) :: il_jmin1 INTEGER(i4) :: il_jmax1 TYPE(TFILE) :: tl_bathy1 ! loop indices INTEGER(i4) :: jl !---------------------------------------------------------------- jl=id_seg !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 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) 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 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 END SELECT !2 -read fine grid domain tl_bathy1=td_bathy1 CALL iom_open(tl_bathy1) !2-1 compute domain create_boundary_get_dom=dom_init( tl_bathy1, & & il_imin1, il_imax1,& & il_jmin1, il_jmax1 ) !2-2 close file CALL iom_close(tl_bathy1) END FUNCTION create_boundary_get_dom !> @endcode !------------------------------------------------------------------- !> @brief !> This subroutine !> !> @details !> !> @author J.Paul !> - Nov, 2013- Initial Version !> !> @param[in] !> @todo !------------------------------------------------------------------- !> @code SUBROUTINE create_boundary_get_coord( td_bathy1, td_dom1, & & td_lon1, td_lat1 ) IMPLICIT NONE ! Argument TYPE(TFILE), INTENT(IN ) :: td_bathy1 TYPE(TDOM) , INTENT(IN ) :: td_dom1 TYPE(TVAR) , INTENT( OUT) :: td_lon1 TYPE(TVAR) , INTENT( OUT) :: td_lat1 ! local variable TYPE(TFILE) :: tl_bathy1 TYPE(TMPP) :: tl_mppbathy1 ! loop indices !---------------------------------------------------------------- !read variables on domain (ugly way to do it, have to work on it) !1 init mpp structure tl_bathy1=td_bathy1 tl_mppbathy1=mpp_init(tl_bathy1) CALL file_clean(tl_bathy1) !2 get processor to be used CALL mpp_get_use( tl_mppbathy1, td_dom1 ) !3 open mpp files CALL iom_mpp_open(tl_mppbathy1) !4 read variable value on domain td_lon1=iom_mpp_read_var(tl_mppbathy1,'longitude',td_dom=td_dom1) td_lat1=iom_mpp_read_var(tl_mppbathy1,'latitude' ,td_dom=td_dom1) !5 close mpp files CALL iom_mpp_close(tl_mppbathy1) !6 clean structure CALL mpp_clean(tl_mppbathy1) END SUBROUTINE create_boundary_get_coord !> @endcode !------------------------------------------------------------------- !> @brief !> This subroutine !> !> @details !> !> @author J.Paul !> - Nov, 2013- Initial Version !> !> @param[in] !> @todo !------------------------------------------------------------------- !> @code SUBROUTINE create_boundary_get_mask( td_level1, td_dom1, & & td_var, td_mask ) IMPLICIT NONE ! Argument TYPE(TFILE), INTENT(IN ) :: td_level1 TYPE(TDOM) , INTENT(IN ) :: td_dom1 TYPE(TVAR) , INTENT(IN ) :: td_var TYPE(TVAR) , INTENT( OUT) :: td_mask ! local variable TYPE(TFILE) :: tl_level1 TYPE(TMPP) :: tl_mpplevel1 ! loop indices !---------------------------------------------------------------- !read variables on domain (ugly way to do it, have to work on it) !1 init mpp structure tl_level1=td_level1 tl_mpplevel1=mpp_init(tl_level1) CALL file_clean(tl_level1) !2 get processor to be used CALL mpp_get_use( tl_mpplevel1, td_dom1 ) !3 open mpp files CALL iom_mpp_open(tl_mpplevel1) !4 read variable value on domain SELECT CASE(TRIM(td_var%c_point)) CASE('T') td_mask=iom_mpp_read_var(tl_mpplevel1,'tlevel',td_dom=td_dom1) CASE('U') td_mask=iom_mpp_read_var(tl_mpplevel1,'ulevel',td_dom=td_dom1) CASE('V') td_mask=iom_mpp_read_var(tl_mpplevel1,'vlevel',td_dom=td_dom1) CASE('F') td_mask=iom_mpp_read_var(tl_mpplevel1,'flevel',td_dom=td_dom1) END SELECT !5 close mpp files CALL iom_mpp_close(tl_mpplevel1) !6 clean structure CALL mpp_clean(tl_mpplevel1) END SUBROUTINE create_boundary_get_mask !> @endcode ! !------------------------------------------------------------------- ! !> @brief ! !> This subroutine ! !> ! !> @details ! !> ! !> @author J.Paul ! !> - Nov, 2013- Initial Version ! !> ! !> @param[in] ! !> @todo ! !------------------------------------------------------------------- ! !> @code ! SUBROUTINE create_boundary_get_var( td_var, td_bdy, & ! & td_coord0, td_dom0, & ! & td_mask, & ! & id_rhoi, id_rhoj ) ! ! IMPLICIT NONE ! ! ! Argument ! TYPE(TVAR) , INTENT(INOUT) :: td_var ! TYPE(TBDY) , INTENT(IN ) :: td_bdy ! TYPE(TFILE), INTENT(IN ) :: td_coord0 ! TYPE(TDOM) , INTENT(IN ) :: td_dom0 ! TYPE(TVAR) , INTENT(IN ) :: td_mask ! INTEGER(I4), INTENT(IN ) :: id_rhoi ! INTEGER(I4), INTENT(IN ) :: id_rhoj ! ! ! local variable ! TYPE(TVAR) :: tl_var0 ! ! TYPE(TDOM) :: tl_dom0 ! ! TYPE(TFILE) :: tl_file0 ! ! TYPE(TMPP) :: tl_mppfile0 ! ! ! loop indices ! INTEGER(i4) :: jk ! INTEGER(i4) :: jl ! !---------------------------------------------------------------- ! ! CALL logger_debug("CREATE BOUNDARY INTERP: read coarse grid"// TRIM(td_var%c_file) ) ! !1- read coarse grid variable on domain ! tl_file0=file_init( TRIM(td_var%c_file) ) ! ! !2- init ! tl_dom0=td_dom0 ! ! !3- add extra band (if possible) to compute interpolation ! CALL dom_add_extra(tl_dom0) ! ! !4- read variables on domain (ugly way to do it, have to work on it) ! !4-1 init mpp structure ! tl_mppfile0=mpp_init(tl_file0) ! ! CALL file_clean(tl_file0) ! ! !4-2 get processor to be used ! CALL mpp_get_use( tl_mppfile0, tl_dom0 ) ! ! !4-3 open mpp files ! CALL iom_mpp_open(tl_mppfile0) ! ! ! check file dimension ! IF( ANY(tl_mppfile0%t_dim(1:2)%i_len /= td_coord0%t_dim(1:2)%i_len) )THEN ! CALL logger_error("CREATE BOUNDARY INTERP: dimension of file "//& ! & TRIM(tl_mppfile0%c_name)//" not conform to those of "//& ! & TRIM(td_coord0%c_name)) ! ELSE ! ! !4-4 read variable value on domain ! tl_var0=iom_mpp_read_var( tl_mppfile0, TRIM(td_var%c_name), & ! & td_dom=tl_dom0 ) ! ! !5- work on variable ! CALL create_boundary_interp(tl_var0, id_rhoi, id_rhoj ) ! ! !6- remove extraband added to domain ! CALL dom_del_extra( tl_var0, tl_dom0, id_rhoi, id_rhoj ) ! ! !6-1 remove extraband added to domain ! CALL dom_clean_extra( tl_dom0 ) ! ! !7- keep only useful point (width) ! ! interpolation could create more point than necessary ! CALL boundary_clean_interp(tl_var0, td_bdy ) ! ! !8- forced min and max value ! CALL var_limit_value(tl_var0) ! ! !9- filter ! CALL filter_fill_value(tl_var0) ! ! td_var=tl_var0 ! ! CALL var_clean(tl_var0) ! ENDIF ! ! !4-5 close mpp files ! CALL iom_mpp_close(tl_mppfile0) ! ! !4-6 clean structure ! CALL mpp_clean(tl_mppfile0) ! ! !5- apply mask ! DO jl=1,td_var%t_dim(4)%i_len ! DO jk=1,td_var%t_dim(3)%i_len ! WHERE( td_mask%d_value(:,:,1,1) < jk ) ! td_var%d_value(:,:,jk,jl)=td_var%d_fill ! END WHERE ! ENDDO ! ENDDO ! ! END SUBROUTINE create_boundary_get_var ! !> @endcode !------------------------------------------------------------------- !> @brief !> This subroutine !> !> @details !> !> @author J.Paul !> - Nov, 2013- Initial Version !> !> @param[in] !> @todo !------------------------------------------------------------------- !> @code SUBROUTINE create_boundary_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 INTEGER(i4) :: il_iext INTEGER(i4) :: il_jext ! loop indices !---------------------------------------------------------------- ! copy variable tl_var=td_var !WARNING: at least two extrabands are required for cubic interpolation il_iext=2 IF( PRESENT(id_iext) ) il_iext=id_iext il_jext=2 IF( PRESENT(id_jext) ) il_jext=id_jext IF( il_iext < 2 .AND. td_var%c_interp(1) == 'cubic' )THEN CALL logger_warn("CREATE BOUNDARY 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 BOUNDARY INTERP: at least extrapolation "//& & "on two points are required with cubic interpolation ") il_jext=2 ENDIF !2- work on variable !2-0 add extraband CALL extrap_add_extrabands(tl_var, il_iext, il_jext) !2-1 extrapolate variable CALL extrap_fill_value( tl_var, 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)) !3- save result td_var=tl_var ! clean variable structure CALL var_clean(tl_var) END SUBROUTINE create_boundary_interp !> @endcode !------------------------------------------------------------------- !> @brief !> This function create variable, filled with matrix value !> !> @details !> A variable is create with the same name that the input variable, !> and with dimension of the coordinate file. !> Then the variable table of value is split into equal subdomain. !> Each subdomain is fill with the linked value of the matrix. !> !> @author J.Paul !> - Nov, 2013- Initial Version !> !> @param[in] td_var : variable structure !> @param[in] td_dom : domain structure !> @param[in] td_coord : coordinate !> @return variable structure !------------------------------------------------------------------- !> @code FUNCTION create_bdy_matrix(td_var, td_dom, td_coord) IMPLICIT NONE ! Argument TYPE(TVAR) , INTENT(IN) :: td_var TYPE(TDOM) , INTENT(IN) :: td_dom TYPE(TFILE), INTENT(IN) :: td_coord ! function TYPE(TVAR) :: create_bdy_matrix ! local variable INTEGER(i4) :: il_ighost INTEGER(i4) :: il_jghost INTEGER(i4) , DIMENSION(2) :: il_xghost INTEGER(i4) , DIMENSION(3) :: il_dim INTEGER(i4) , DIMENSION(3) :: il_size INTEGER(i4) , DIMENSION(3) :: il_rest INTEGER(i4) , DIMENSION(:) , ALLOCATABLE :: il_ishape INTEGER(i4) , DIMENSION(:) , ALLOCATABLE :: il_jshape INTEGER(i4) , DIMENSION(:) , ALLOCATABLE :: il_kshape REAL(dp) , DIMENSION(:,:,:,:), ALLOCATABLE :: dl_value TYPE(TVAR) :: tl_lon TYPE(TVAR) :: tl_lat TYPE(TVAR) :: tl_var TYPE(TDIM) , DIMENSION(ip_maxdim) :: tl_dim ! loop indices INTEGER(i4) :: ji INTEGER(i4) :: jj INTEGER(i4) :: jk !---------------------------------------------------------------- !1- read output grid tl_lon=iom_read_var(td_coord,'longitude') tl_lat=iom_read_var(td_coord,'latitude') !2- look for ghost cell il_xghost(:)=grid_get_ghost( tl_lon, tl_lat ) il_ighost=il_xghost(1)*ig_ghost il_jghost=il_xghost(2)*ig_ghost !3- write value on grid !3-1 get matrix dimension il_dim(:)=td_var%t_dim(1:3)%i_len !3-2 output dimension tl_dim(:)=tl_lon%t_dim(:) ! remove ghost cell tl_dim(1)%i_len=tl_dim(1)%i_len - 2*il_xghost(1)*ig_ghost tl_dim(2)%i_len=tl_dim(2)%i_len - 2*il_xghost(2)*ig_ghost !3-3 split output domain in N subdomain depending of matrix dimension il_size(:) = tl_dim(1:3)%i_len / il_dim(:) il_rest(:) = MOD(tl_dim(1:3)%i_len, il_dim(:)) ALLOCATE( il_ishape(il_dim(1)+1) ) il_ishape(:)=0 DO ji=2,il_dim(1)+1 il_ishape(ji)=il_ishape(ji-1)+il_size(1) ENDDO ! add rest to last cell il_ishape(il_dim(1)+1)=il_ishape(il_dim(1)+1)+il_rest(1) ALLOCATE( il_jshape(il_dim(2)+1) ) il_jshape(:)=0 DO jj=2,il_dim(2)+1 il_jshape(jj)=il_jshape(jj-1)+il_size(2) ENDDO ! add rest to last cell il_jshape(il_dim(2)+1)=il_jshape(il_dim(2)+1)+il_rest(2) ALLOCATE( il_kshape(il_dim(3)+1) ) il_kshape(:)=0 DO jk=2,il_dim(3)+1 il_kshape(jk)=il_kshape(jk-1)+il_size(3) ENDDO ! add rest to last cell il_kshape(il_dim(3)+1)=il_kshape(il_dim(3)+1)+il_rest(3) !3-3 write ouput table of value ALLOCATE(dl_value( tl_dim(1)%i_len, & & tl_dim(2)%i_len, & & tl_dim(3)%i_len, & & tl_dim(4)%i_len) ) dl_value(:,:,:,:)=0 DO jk=2,il_dim(3)+1 DO jj=2,il_dim(2)+1 DO ji=2,il_dim(1)+1 dl_value( 1+il_ishape(ji-1):il_ishape(ji), & & 1+il_jshape(jj-1):il_jshape(jj), & & 1+il_kshape(jk-1):il_kshape(jk), & & 1 ) = td_var%d_value(ji-1,jj-1,jk-1,1) ENDDO ENDDO ENDDO !3-4 initialise variable with value tl_var=var_init(TRIM(td_var%c_name),dl_value(:,:,:,:)) DEALLOCATE(dl_value) !4- add ghost cell CALL grid_add_ghost(tl_var,il_ighost,il_jghost) !5- save result create_bdy_matrix=tl_var END FUNCTION create_bdy_matrix !> @endcode !------------------------------------------------------------------- !> @brief !> This subroutine !> !> @details !> !> @author J.Paul !> - Nov, 2013- Initial Version !> !> @param[in] !> @todo !------------------------------------------------------------------- !> @code SUBROUTINE create_bdy_use_mask( td_var, td_mask ) IMPLICIT NONE ! Argument TYPE(TVAR) , INTENT(INOUT) :: td_var TYPE(TVAR), DIMENSION(:), INTENT(IN ) :: td_mask ! local variable INTEGER(i4), DIMENSION(:,:), ALLOCATABLE :: il_mask ! loop indices INTEGER(i4) :: jk INTEGER(i4) :: jl !---------------------------------------------------------------- ALLOCATE( il_mask(td_var%t_dim(1)%i_len, & & td_var%t_dim(2)%i_len) ) SELECT CASE(TRIM(td_var%c_point)) CASE('T') il_mask(:,:)=INT(td_mask(jp_T)%d_value(:,:,1,1)) CASE('U') il_mask(:,:)=INT(td_mask(jp_U)%d_value(:,:,1,1)) CASE('V') il_mask(:,:)=INT(td_mask(jp_V)%d_value(:,:,1,1)) CASE('F') il_mask(:,:)=INT(td_mask(jp_F)%d_value(:,:,1,1)) END SELECT DO jl=1,td_var%t_dim(4)%i_len DO jk=1,td_var%t_dim(3)%i_len WHERE( il_mask(:,:) < jk ) td_var%d_value(:,:,jk,jl)=td_var%d_fill ENDDO ENDDO DEALLOCATE( il_mask ) END SUBROUTINE create_bdy_use_mask !> @endcode !------------------------------------------------------------------- !> @brief !> !> @details !> !> @author J.Paul !> - 2013- Initial Version !> !------------------------------------------------------------------- !> @code FUNCTION create_bdy_get_level(td_level, td_dom) IMPLICIT NONE ! Argument TYPE(TVAR), DIMENSION(:), INTENT(IN) :: td_level TYPE(TDOM) , INTENT(IN) :: td_dom ! function TYPE(TVAR), DIMENSION(ig_npoint) :: create_bdy_get_level ! local variable TYPE(TVAR), DIMENSION(ig_npoint) :: tl_var ! loop indices INTEGER(i4) :: ji !---------------------------------------------------------------- IF( SIZE(td_level(:)) /= ig_npoint )THEN CALL logger_error("CREATE BDY GET LEVEL: invalid dimension. "//& & "check input table of level.") ELSE !tl_var(1:ig_npoint)=td_level(1:ig_npoint) create_bdy_get_level(:)=tl_var(:) DO ji=1,ig_npoint tl_var(ji)=td_level(ji) IF( ASSOCIATED(tl_var(ji)%d_value) ) DEALLOCATE( tl_var(ji)%d_value ) tl_var(ji)%t_dim(1)%i_len=td_dom%t_dim(1)%i_len tl_var(ji)%t_dim(2)%i_len=td_dom%t_dim(2)%i_len ALLOCATE(tl_var(ji)%d_value(tl_var(ji)%t_dim(1)%i_len, & & tl_var(ji)%t_dim(2)%i_len, & & tl_var(ji)%t_dim(3)%i_len, & & tl_var(ji)%t_dim(4)%i_len) ) tl_var(ji)%d_value(:,:,:,:) = & & td_level(ji)%d_value( td_dom%i_imin:td_dom%i_imax, & & td_dom%i_jmin:td_dom%i_jmax, :, : ) ENDDO !4 save result create_bdy_get_level(:)=tl_var(:) ENDIF END FUNCTION create_bdy_get_level !> @endcode END PROGRAM create_boundary