Changeset 6225 for branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/TOOLS/SIREN/src/merge_bathy.f90
- Timestamp:
- 2016-01-08T10:35:19+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/TOOLS/SIREN/src/merge_bathy.f90
r4213 r6225 7 7 ! 8 8 ! DESCRIPTION: 9 !> @file 9 10 !> @brief 10 11 !> This program merge bathymetry file at boundaries. 11 12 !> 12 13 !> @details 14 !> @section sec1 method 13 15 !> Coarse grid Bathymetry is interpolated on fine grid. 14 !> Then fine Bathymetry and refined coarse bathymetry are merged at boundaries. 15 !> 16 !> BathyFine= weight * BathyCoarse + (1-weight)*BathyFine 17 !> 18 !> The weight function used is : 0.5 + 0.5*COS( (pi*dist) / width ) 19 !> 20 !> @author 21 !> J.Paul 16 !> Then fine Bathymetry and refined coarse bathymetry are merged at boundaries.<br/> 17 !> @f[BathyFine= Weight * BathyCoarse + (1-Weight)*BathyFine@f] 18 !> The weight function used is :<br/> 19 !> @f[Weight = 0.5 + 0.5*COS( \frac{\pi*dist}{width} )@f]<br/> 20 !> with 21 !> - dist : number of point to border 22 !> - width : boundary size 23 !> 24 !> @section sec2 how to 25 !> to merge bathymetry file:<br/> 26 !> @code{.sh} 27 !> ./SIREN/bin/merge_bathy merge_bathy.nam 28 !> @endcode 29 !> 30 !> @note 31 !> you could find a template of the namelist in templates directory. 32 !> 33 !> merge_bathy.nam comprise 8 namelists: 34 !> - logger namelist (namlog) 35 !> - config namelist (namcfg) 36 !> - coarse grid namelist (namcrs) 37 !> - fine grid namelist (namfin) 38 !> - variable namelist (namvar) 39 !> - nesting namelist (namnst) 40 !> - boundary namelist (nambdy) 41 !> - output namelist (namout) 42 !> 43 !> @note 44 !> All namelists have to be in file merge_bathy.nam, 45 !> however variables of those namelists are all optional. 46 !> 47 !> * _logger namelist (namlog)_: 48 !> - cn_logfile : logger filename 49 !> - cn_verbosity : verbosity ('trace','debug','info', 50 !> 'warning','error','fatal','none') 51 !> - in_maxerror : maximum number of error allowed 52 !> 53 !> * _config namelist (namcfg)_: 54 !> - cn_varcfg : variable configuration file (see ./SIREN/cfg/variable.cfg) 55 !> 56 !> * _coarse grid namelist (namcrs)_: 57 !> - cn_bathy0 : bathymetry file 58 !> - in_perio0 : NEMO periodicity index (see Model Boundary Condition in 59 !> [NEMO documentation](http://www.nemo-ocean.eu/About-NEMO/Reference-manuals)) 60 !> 61 !> * _fine grid namelist (namfin)_: 62 !> - cn_bathy1 : bathymetry file 63 !> - in_perio1 : NEMO periodicity index 64 !> 65 !> * _variable namelist (namvar)_: 66 !> - cn_varinfo : list of variable and extra information about request(s) 67 !> to be used (separated by ',').<br/> 68 !> each elements of *cn_varinfo* is a string character.<br/> 69 !> it is composed of the variable name follow by ':', 70 !> then request(s) to be used on this variable.<br/> 71 !> request could be: 72 !> - int = interpolation method 73 !> 74 !> requests must be separated by ';'.<br/> 75 !> order of requests does not matter.<br/> 76 !> 77 !> informations about available method could be find in 78 !> @ref interp modules.<br/> 79 !> Example: 'bathymetry: int=cubic' 80 !> @note 81 !> If you do not specify a method which is required, 82 !> default one is apply. 83 !> @warning 84 !> variable name must be __Bathymetry__ here. 85 !> 86 !> * _nesting namelist (namnst)_: 87 !> - in_rhoi : refinement factor in i-direction 88 !> - in_rhoj : refinement factor in j-direction 89 !> 90 !> * _boundary namelist (nambdy)_: 91 !> - ln_north : use north boundary or not 92 !> - ln_south : use south boundary or not 93 !> - ln_east : use east boundary or not 94 !> - ln_west : use west boundary or not 95 !> - cn_north : north boundary indices on fine grid<br/> 96 !> *cn_north* is a string character defining boundary 97 !> segmentation.<br/> 98 !> segments are separated by '|'.<br/> 99 !> each segments of the boundary is composed of: 100 !> - indice of velocity (orthogonal to boundary .ie. 101 !> for north boundary, J-indice). 102 !> - indice of segment start (I-indice for north boundary) 103 !> - indice of segment end (I-indice for north boundary)<br/> 104 !> indices must be separated by ':' .<br/> 105 !> - optionally, boundary size could be added between '(' and ')' 106 !> in the first segment defined. 107 !> @note 108 !> boundary size is the same for all segments of one boundary. 109 !> 110 !> Examples: 111 !> - cn_north='index1,first1:last1(width)' 112 !> - cn_north='index1(width),first1:last1|index2,first2:last2' 113 !> 114 !> - cn_south : south boundary indices on fine grid<br/> 115 !> - cn_east : east boundary indices on fine grid<br/> 116 !> - cn_west : west boundary indices on fine grid<br/> 117 !> - ln_oneseg: use only one segment for each boundary or not 118 !> 119 !> * _output namelist (namout)_: 120 !> - cn_fileout : merged bathymetry file 121 !> 122 !> @author J.Paul 22 123 ! REVISION HISTORY: 23 !> @date Nov, 2013 - Initial Version 24 ! 124 !> @date November, 2013 - Initial Version 125 !> @date Sepember, 2014 126 !> - add header for user 127 !> @date July, 2015 128 !> - extrapolate all land points 129 !> - add attributes with boundary string character (as in namelist) 130 !> 25 131 !> @note Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 26 !>27 !> @todo28 132 !---------------------------------------------------------------------- 29 !> @code30 133 PROGRAM merge_bathy 31 134 … … 43 146 USE boundary ! boundary manager 44 147 USE iom ! I/O manager 45 USE dom ! domain manager46 148 USE grid ! grid manager 47 149 USE extrap ! extrapolation manager … … 49 151 USE filter ! filter manager 50 152 USE mpp ! MPP manager 153 USE dom ! domain manager 51 154 USE iom_mpp ! MPP I/O manager 155 USE iom_dom ! DOM I/O manager 52 156 53 157 IMPLICIT NONE … … 56 160 CHARACTER(LEN=lc) :: cl_namelist 57 161 CHARACTER(LEN=lc) :: cl_date 162 CHARACTER(LEN=lc) :: cl_tmp 58 163 59 164 INTEGER(i4) :: il_narg 60 165 INTEGER(i4) :: il_status 61 166 INTEGER(i4) :: il_fileid 62 INTEGER(i4) :: il_atti d167 INTEGER(i4) :: il_attind 63 168 INTEGER(i4) :: il_imin0 64 169 INTEGER(i4) :: il_imax0 65 170 INTEGER(i4) :: il_jmin0 66 171 INTEGER(i4) :: il_jmax0 172 INTEGER(i4) :: il_shift 67 173 INTEGER(i4) , DIMENSION(ip_maxdim) :: il_rho 68 INTEGER(i4) , DIMENSION(2,2) :: il_offset 69 INTEGER(i4) , DIMENSION(2,2,2) :: il_ind 70 ! INTEGER(i4) , DIMENSION(:,:,:,:), ALLOCATABLE :: il_value 174 INTEGER(i4) , DIMENSION(2,2) :: il_ind 71 175 72 176 LOGICAL :: ll_exist … … 76 180 REAL(dp) , DIMENSION(:,:,:,:), ALLOCATABLE :: dl_weight 77 181 78 TYPE(T FILE):: tl_bathy079 TYPE(T FILE):: tl_bathy1182 TYPE(TMPP) :: tl_bathy0 183 TYPE(TMPP) :: tl_bathy1 80 184 TYPE(TFILE) :: tl_fileout 81 185 … … 85 189 TYPE(TVAR) :: tl_lon 86 190 TYPE(TVAR) :: tl_lat 87 ! TYPE(TVAR) :: tl_depth88 ! TYPE(TVAR) :: tl_time89 191 90 192 TYPE(TDIM) , DIMENSION(ip_maxdim) :: tl_dim … … 98 200 99 201 ! namelist variable 202 ! namlog 100 203 CHARACTER(LEN=lc) :: cn_logfile = 'merge_bathy.log' 101 204 CHARACTER(LEN=lc) :: cn_verbosity = 'warning' 102 205 INTEGER(i4) :: in_maxerror = 5 206 207 ! namcfg 208 CHARACTER(LEN=lc) :: cn_varcfg = 'variable.cfg' 209 210 ! namcrs 103 211 CHARACTER(LEN=lc) :: cn_bathy0 = '' 104 212 INTEGER(i4) :: in_perio0 = -1 105 213 214 ! namfin 106 215 CHARACTER(LEN=lc) :: cn_bathy1 = '' 107 216 INTEGER(i4) :: in_perio1 = -1 108 217 109 CHARACTER(LEN=lc) :: cn_varcfg = 'variable.cfg' 110 111 CHARACTER(LEN=lc), DIMENSION(ig_maxvar) :: cn_varinfo = '' 112 113 INTEGER(i4) :: in_imin0 = 0 114 INTEGER(i4) :: in_imax0 = 0 115 INTEGER(i4) :: in_jmin0 = 0 116 INTEGER(i4) :: in_jmax0 = 0 218 ! namvar 219 CHARACTER(LEN=lc), DIMENSION(ip_maxvar) :: cn_varinfo = '' 220 221 ! namnst 117 222 INTEGER(i4) :: in_rhoi = 0 118 223 INTEGER(i4) :: in_rhoj = 0 119 224 225 ! nambdy 120 226 LOGICAL :: ln_north = .TRUE. 121 227 LOGICAL :: ln_south = .TRUE. … … 128 234 CHARACTER(LEN=lc) :: cn_west = '' 129 235 236 ! namout 130 237 CHARACTER(LEN=lc) :: cn_fileout = 'bathy_merged.nc' 131 238 !------------------------------------------------------------------- … … 133 240 NAMELIST /namlog/ & !< logger namelist 134 241 & cn_logfile, & !< log file 135 & cn_verbosity !< log verbosity 242 & cn_verbosity, & !< log verbosity 243 & in_maxerror !< logger maximum error 136 244 137 245 NAMELIST /namcfg/ & !< config namelist … … 152 260 153 261 NAMELIST /namnst/ & !< nesting namelist 154 & in_imin0, & !< i-direction lower left point indice on coarse grid155 & in_imax0, & !< i-direction upper right point indice on coarse grid156 & in_jmin0, & !< j-direction lower left point indice on coarse grid157 & in_jmax0, & !< j-direction upper right point indice on coarse grid158 262 & in_rhoi, & !< refinement factor in i-direction 159 263 & in_rhoj !< refinement factor in j-direction … … 171 275 172 276 NAMELIST /namout/ & !< output namelist 173 & cn_fileout !< fine grid merged bathymetry file277 & cn_fileout !< fine grid merged bathymetry file 174 278 !------------------------------------------------------------------- 175 279 176 ! 1-namelist177 ! 1-1get namelist280 ! namelist 281 ! get namelist 178 282 il_narg=COMMAND_ARGUMENT_COUNT() !f03 intrinsec 179 283 IF( il_narg/=1 )THEN … … 184 288 ENDIF 185 289 186 ! 1-2read namelist290 ! read namelist 187 291 INQUIRE(FILE=TRIM(cl_namelist), EXIST=ll_exist) 188 292 IF( ll_exist )THEN … … 203 307 204 308 READ( il_fileid, NML = namlog ) 205 ! 1-2-1define log file206 CALL logger_open(TRIM(cn_logfile),TRIM(cn_verbosity) )309 ! define log file 310 CALL logger_open(TRIM(cn_logfile),TRIM(cn_verbosity),in_maxerror) 207 311 CALL logger_header() 208 312 209 313 READ( il_fileid, NML = namcfg ) 210 ! 1-2-2get variable extra information314 ! get variable extra information 211 315 CALL var_def_extra(TRIM(cn_varcfg)) 212 316 … … 214 318 READ( il_fileid, NML = namfin ) 215 319 READ( il_fileid, NML = namvar ) 216 ! 1-2-3add user change in extra information320 ! add user change in extra information 217 321 CALL var_chg_extra(cn_varinfo) 218 322 … … 234 338 ENDIF 235 339 236 ! 2-open files340 ! open files 237 341 IF( TRIM(cn_bathy0) /= '' )THEN 238 tl_bathy0= file_init(TRIM(cn_bathy0),id_perio=in_perio0)239 CALL iom_open(tl_bathy0)342 tl_bathy0=mpp_init( file_init(TRIM(cn_bathy0)), id_perio=in_perio0) 343 CALL grid_get_info(tl_bathy0) 240 344 ELSE 241 345 CALL logger_fatal("MERGE BATHY: can not find coarse grid bathymetry "//& … … 244 348 245 349 IF( TRIM(cn_bathy1) /= '' )THEN 246 tl_bathy1= file_init(TRIM(cn_bathy1),id_perio=in_perio1)247 CALL iom_open(tl_bathy1)350 tl_bathy1=mpp_init( file_init(TRIM(cn_bathy1)), id_perio=in_perio1) 351 CALL grid_get_info(tl_bathy1) 248 352 ELSE 249 353 CALL logger_fatal("MERGE BATHY: can not find fine grid bathymetry "//& … … 251 355 ENDIF 252 356 253 ! 3-check254 ! 3-1check output file do not already exist357 ! check 358 ! check output file do not already exist 255 359 INQUIRE(FILE=TRIM(cn_fileout), EXIST=ll_exist) 256 360 IF( ll_exist )THEN … … 259 363 ENDIF 260 364 261 ! 3-2check namelist262 ! 3-2-1check refinament factor365 ! check namelist 366 ! check refinament factor 263 367 il_rho(:)=1 264 368 IF( in_rhoi < 1 .OR. in_rhoj < 1 )THEN … … 270 374 ENDIF 271 375 272 !3-2-2 check domain indices 273 IF( in_imin0 < 1 .OR. in_imax0 < 1 .OR. in_jmin0 < 1 .OR. in_jmax0 < 1)THEN 274 ! compute coarse grid indices around fine grid 275 il_ind(:,:,:)=grid_get_coarse_index(tl_bathy0, tl_bathy1 ) 276 277 il_imin0=il_ind(1,1,1) ; il_imax0=il_ind(1,2,1) 278 il_jmin0=il_ind(2,1,1) ; il_jmax0=il_ind(2,2,1) 279 280 il_offset(:,:)=il_ind(:,:,2) 281 ELSE 282 il_imin0=in_imin0 ; il_imax0=in_imax0 283 il_jmin0=in_jmin0 ; il_jmax0=in_jmax0 284 285 il_offset(1,:)=FLOOR(REAL(il_rho(jp_I)-1,dp)*0.5) 286 il_offset(2,:)=FLOOR(REAL(il_rho(jp_J)-1,dp)*0.5) 287 ENDIF 288 289 !3-3 check domain validity 376 ! check domain indices 377 ! compute coarse grid indices around fine grid 378 il_ind(:,:)=grid_get_coarse_index(tl_bathy0, tl_bathy1, & 379 & id_rho=il_rho(:) ) 380 381 il_imin0=il_ind(1,1) ; il_imax0=il_ind(1,2) 382 il_jmin0=il_ind(2,1) ; il_jmax0=il_ind(2,2) 383 384 ! check domain validity 290 385 CALL grid_check_dom(tl_bathy0, il_imin0, il_imax0, il_jmin0, il_jmax0) 291 386 292 ! 3-4check coincidence between coarse and fine grid387 ! check coincidence between coarse and fine grid 293 388 CALL grid_check_coincidence( tl_bathy0, tl_bathy1, & 294 389 & il_imin0, il_imax0, & … … 296 391 & il_rho(:) ) 297 392 298 !4- read or compute boundary 299 tl_var=iom_read_var(tl_bathy1,'Bathymetry') 393 ! open mpp files 394 CALL iom_mpp_open(tl_bathy1) 395 396 ! read or compute boundary 397 tl_var=iom_mpp_read_var(tl_bathy1,'Bathymetry') 398 399 ! close mpp files 400 CALL iom_mpp_close(tl_bathy1) 300 401 301 402 tl_bdy(:)=boundary_init(tl_var, ln_north, ln_south, ln_east, ln_west, & … … 303 404 & ln_oneseg ) 304 405 305 ! 5-get boundary on coarse grid306 ! 5-1 define refined bathymetry table(for coarse grid)406 ! get boundary on coarse grid 407 ! define refined bathymetry array (for coarse grid) 307 408 dl_fill=tl_var%d_fill 308 409 ALLOCATE( dl_refined(tl_var%t_dim(1)%i_len, & … … 313 414 dl_refined(:,:,:,:)=dl_fill 314 415 315 ! 5-2 define weight table416 ! define weight array 316 417 ALLOCATE( dl_weight(tl_var%t_dim(1)%i_len, & 317 418 & tl_var%t_dim(2)%i_len, & … … 320 421 dl_weight(:,:,:,:)=dl_fill 321 422 322 ! 5-3compute coarse grid refined bathymetry on boundary.423 ! compute coarse grid refined bathymetry on boundary. 323 424 DO jk=1,ip_ncard 324 325 425 CALL merge_bathy_get_boundary(tl_bathy0, tl_bathy1, tl_bdy(jk), & 326 426 & il_rho(:), & … … 330 430 ENDDO 331 431 332 ! 6-merge bathy on boundary432 ! merge bathy on boundary 333 433 DO jl=1,tl_var%t_dim(4)%i_len 334 434 DO jk=1,tl_var%t_dim(3)%i_len … … 348 448 DEALLOCATE(dl_refined) 349 449 350 ! 7-create file450 ! create file 351 451 tl_fileout=file_init(TRIM(cn_fileout),id_perio=in_perio1) 352 452 353 ! 7-1add dimension354 tl_dim(:)= tl_var%t_dim(:)453 ! add dimension 454 tl_dim(:)=dim_copy(tl_var%t_dim(:)) 355 455 356 456 DO ji=1,ip_maxdim … … 358 458 ENDDO 359 459 360 ! 7-2add variables460 ! add variables 361 461 IF( ALL( tl_dim(1:2)%l_use ) )THEN 362 363 tl_lon=iom_read_var(tl_bathy1,'longitude') 462 ! open mpp files 463 CALL iom_mpp_open(tl_bathy1) 464 465 ! add longitude 466 tl_lon=iom_mpp_read_var(tl_bathy1,'longitude') 364 467 CALL file_add_var(tl_fileout, tl_lon) 365 468 CALL var_clean(tl_lon) 366 469 367 tl_lat=iom_read_var(tl_bathy1,'latitude') 470 ! add latitude 471 tl_lat=iom_mpp_read_var(tl_bathy1,'latitude') 368 472 CALL file_add_var(tl_fileout, tl_lat) 369 473 CALL var_clean(tl_lat) 370 474 475 ! close mpp files 476 CALL iom_mpp_close(tl_bathy1) 371 477 ENDIF 372 478 … … 375 481 376 482 ! only 2 first dimension to be used 377 tl_dim(:)= tl_fileout%t_dim(:)483 tl_dim(:)=dim_copy(tl_fileout%t_dim(:)) 378 484 tl_dim(3:4)%l_use=.FALSE. 379 485 tl_var=var_init('weight',dl_weight(:,:,:,:),td_dim=tl_dim(:),dd_fill=dl_fill) … … 381 487 CALL var_clean(tl_var) 382 488 383 ! 7-3add some attribute489 ! add some attribute 384 490 tl_att=att_init("Created_by","SIREN merge_bathy") 385 491 CALL file_add_att(tl_fileout, tl_att) … … 395 501 CALL file_add_att(tl_fileout, tl_att) 396 502 397 ! a voir398 503 ! add attribute periodicity 399 il_atti d=0504 il_attind=0 400 505 IF( ASSOCIATED(tl_fileout%t_att) )THEN 401 il_atti d=att_get_id(tl_fileout%t_att(:),'periodicity')402 ENDIF 403 IF( tl_bathy1%i_perio >= 0 .AND. il_atti d == 0 )THEN506 il_attind=att_get_index(tl_fileout%t_att(:),'periodicity') 507 ENDIF 508 IF( tl_bathy1%i_perio >= 0 .AND. il_attind == 0 )THEN 404 509 tl_att=att_init('periodicity',tl_bathy1%i_perio) 405 510 CALL file_add_att(tl_fileout,tl_att) 406 511 ENDIF 407 512 408 il_atti d=0513 il_attind=0 409 514 IF( ASSOCIATED(tl_fileout%t_att) )THEN 410 il_atti d=att_get_id(tl_fileout%t_att(:),'ew_overlap')411 ENDIF 412 IF( tl_bathy1%i_ew >= 0 .AND. il_atti d == 0 )THEN515 il_attind=att_get_index(tl_fileout%t_att(:),'ew_overlap') 516 ENDIF 517 IF( tl_bathy1%i_ew >= 0 .AND. il_attind == 0 )THEN 413 518 tl_att=att_init('ew_overlap',tl_bathy1%i_ew) 414 519 CALL file_add_att(tl_fileout,tl_att) 415 520 ENDIF 416 521 417 !7-4 create file 522 523 IF( tl_bdy(jp_north)%l_use )THEN 524 ! add shift on north boundary 525 ! boundary compute on T point but express on U or V point 526 il_shift=1 527 528 cl_tmp=TRIM(fct_str(tl_bdy(jp_north)%t_seg(1)%i_index-il_shift))//','//& 529 & TRIM(fct_str(tl_bdy(jp_north)%t_seg(1)%i_first))//':'//& 530 & TRIM(fct_str(tl_bdy(jp_north)%t_seg(1)%i_last))//& 531 & '('//TRIM(fct_str(tl_bdy(jp_north)%t_seg(1)%i_width))//')' 532 DO ji=2,tl_bdy(jp_north)%i_nseg 533 cl_tmp=TRIM(cl_tmp)//'|'//& 534 & TRIM(fct_str(tl_bdy(jp_north)%t_seg(ji)%i_index-il_shift))//','//& 535 & TRIM(fct_str(tl_bdy(jp_north)%t_seg(ji)%i_first))//':'//& 536 & TRIM(fct_str(tl_bdy(jp_north)%t_seg(ji)%i_last)) 537 ENDDO 538 tl_att=att_init("bdy_north",TRIM(cl_tmp)) 539 CALL file_add_att(tl_fileout, tl_att) 540 ENDIF 541 542 IF( tl_bdy(jp_south)%l_use )THEN 543 544 cl_tmp=TRIM(fct_str(tl_bdy(jp_south)%t_seg(1)%i_index))//','//& 545 & TRIM(fct_str(tl_bdy(jp_south)%t_seg(1)%i_first))//':'//& 546 & TRIM(fct_str(tl_bdy(jp_south)%t_seg(1)%i_last))//& 547 & '('//TRIM(fct_str(tl_bdy(jp_south)%t_seg(1)%i_width))//')' 548 DO ji=2,tl_bdy(jp_south)%i_nseg 549 cl_tmp=TRIM(cl_tmp)//'|'//& 550 & TRIM(fct_str(tl_bdy(jp_south)%t_seg(ji)%i_index))//','//& 551 & TRIM(fct_str(tl_bdy(jp_south)%t_seg(ji)%i_first))//':'//& 552 & TRIM(fct_str(tl_bdy(jp_south)%t_seg(ji)%i_last)) 553 ENDDO 554 555 tl_att=att_init("bdy_south",TRIM(cl_tmp)) 556 CALL file_add_att(tl_fileout, tl_att) 557 ENDIF 558 559 IF( tl_bdy(jp_east)%l_use )THEN 560 ! add shift on east boundary 561 ! boundary compute on T point but express on U or V point 562 il_shift=1 563 564 cl_tmp=TRIM(fct_str(tl_bdy(jp_east)%t_seg(1)%i_index-il_shift))//','//& 565 & TRIM(fct_str(tl_bdy(jp_east)%t_seg(1)%i_first))//':'//& 566 & TRIM(fct_str(tl_bdy(jp_east)%t_seg(1)%i_last))//& 567 & '('//TRIM(fct_str(tl_bdy(jp_east)%t_seg(1)%i_width))//')' 568 DO ji=2,tl_bdy(jp_east)%i_nseg 569 cl_tmp=TRIM(cl_tmp)//'|'//& 570 & TRIM(fct_str(tl_bdy(jp_east)%t_seg(ji)%i_index-il_shift))//','//& 571 & TRIM(fct_str(tl_bdy(jp_east)%t_seg(ji)%i_first))//':'//& 572 & TRIM(fct_str(tl_bdy(jp_east)%t_seg(ji)%i_last)) 573 ENDDO 574 575 tl_att=att_init("bdy_east",TRIM(cl_tmp)) 576 CALL file_add_att(tl_fileout, tl_att) 577 ENDIF 578 579 IF( tl_bdy(jp_west)%l_use )THEN 580 581 cl_tmp=TRIM(fct_str(tl_bdy(jp_west)%t_seg(1)%i_index))//','//& 582 & TRIM(fct_str(tl_bdy(jp_west)%t_seg(1)%i_first))//':'//& 583 & TRIM(fct_str(tl_bdy(jp_west)%t_seg(1)%i_last))//& 584 & '('//TRIM(fct_str(tl_bdy(jp_west)%t_seg(1)%i_width))//')' 585 DO ji=2,tl_bdy(jp_west)%i_nseg 586 cl_tmp=TRIM(cl_tmp)//'|'//& 587 & TRIM(fct_str(tl_bdy(jp_west)%t_seg(ji)%i_index))//','//& 588 & TRIM(fct_str(tl_bdy(jp_west)%t_seg(ji)%i_first))//':'//& 589 & TRIM(fct_str(tl_bdy(jp_west)%t_seg(ji)%i_last)) 590 ENDDO 591 592 tl_att=att_init("bdy_west",TRIM(cl_tmp)) 593 CALL file_add_att(tl_fileout, tl_att) 594 ENDIF 595 596 ! create file 418 597 CALL iom_create(tl_fileout) 419 598 420 ! 7-5write file599 ! write file 421 600 CALL iom_write_file(tl_fileout) 422 601 423 ! 7-6close file602 ! close file 424 603 CALL iom_close(tl_fileout) 425 604 426 CALL iom_close(tl_bathy1) 427 CALL iom_close(tl_bathy0) 428 429 !8- clean 605 ! clean 606 CALL att_clean(tl_att) 430 607 CALL file_clean(tl_fileout) 431 CALL file_clean(tl_bathy1)432 CALL file_clean(tl_bathy0)608 CALL mpp_clean(tl_bathy1) 609 CALL mpp_clean(tl_bathy0) 433 610 DEALLOCATE(dl_weight) 611 CALL boundary_clean(tl_bdy(:)) 434 612 435 613 ! close log file … … 437 615 CALL logger_close() 438 616 439 !> @endcode440 617 CONTAINS 441 618 !------------------------------------------------------------------- 442 619 !> @brief 443 !> This subroutine 620 !> This subroutine compute refined bathymetry on boundary from coarse grid. 444 621 !> 445 !> @details 622 !> @author J.Paul 623 !> @date November, 2013 - Initial Version 446 624 !> 447 !> @author J.Paul 448 !> - Nov, 2013- Initial Version 449 !> 450 !> @param[in] 451 !> @todo 625 !> @param[in] td_bathy0 coarse grid bathymetry file structure 626 !> @param[in] td_bathy1 fine grid bathymetry file structure 627 !> @param[in] td_bdy boundary structure 628 !> @param[in] id_rho array of refinement factor 629 !> @param[inout] dd_refined array of refined bathymetry 630 !> @param[inout] dd_weight array of weight 631 !> @param[in] dd_fill fillValue 452 632 !------------------------------------------------------------------- 453 !> @code454 633 SUBROUTINE merge_bathy_get_boundary( td_bathy0, td_bathy1, td_bdy, & 455 634 & id_rho, & … … 459 638 460 639 ! Argument 461 TYPE(T FILE), INTENT(IN ) :: td_bathy0462 TYPE(T FILE), INTENT(IN ) :: td_bathy1640 TYPE(TMPP) , INTENT(IN ) :: td_bathy0 641 TYPE(TMPP) , INTENT(IN ) :: td_bathy1 463 642 TYPE(TBDY) , INTENT(IN ) :: td_bdy 464 643 INTEGER(i4), DIMENSION(:) , INTENT(IN ) :: id_rho … … 478 657 INTEGER(i4) :: il_jmax0 479 658 480 INTEGER(i4) :: il_imin481 INTEGER(i4) :: il_imax482 INTEGER(i4) :: il_jmin483 INTEGER(i4) :: il_jmax484 485 659 INTEGER(i4), DIMENSION(2,2) :: il_offset 486 INTEGER(i4), DIMENSION(2,2 ,2):: il_ind660 INTEGER(i4), DIMENSION(2,2) :: il_ind 487 661 488 662 REAL(dp) , DIMENSION(:) , ALLOCATABLE :: dl_tmp1d … … 494 668 TYPE(TVAR) :: tl_lat1 495 669 496 TYPE(TFILE) :: tl_bathy1 497 TYPE(TFILE) :: tl_bathy0 498 499 TYPE(TMPP) :: tl_mppbathy1 500 TYPE(TMPP) :: tl_mppbathy0 670 TYPE(TMPP) :: tl_bathy1 671 TYPE(TMPP) :: tl_bathy0 501 672 502 673 TYPE(TDOM) :: tl_dom1 … … 510 681 IF( td_bdy%l_use )THEN 511 682 DO jl=1,td_bdy%i_nseg 512 513 !1- get boundary definition 683 ! get boundary definition 514 684 SELECT CASE(TRIM(td_bdy%c_card)) 515 685 CASE('north') … … 520 690 il_jmax1=td_bdy%t_seg(jl)%i_index 521 691 522 il_imin=1523 il_imax=il_imax1-il_imin1+1524 il_jmin=td_bdy%t_seg(jl)%i_width525 il_jmax=1526 527 692 CASE('south') 528 693 … … 532 697 il_jmax1=td_bdy%t_seg(jl)%i_index+(td_bdy%t_seg(jl)%i_width-1) 533 698 534 il_imin=1535 il_imax=il_imax1-il_imin1+1536 il_jmin=1537 il_jmax=td_bdy%t_seg(jl)%i_width538 539 699 CASE('east') 540 700 … … 544 704 il_jmax1=td_bdy%t_seg(jl)%i_last 545 705 546 il_imin=td_bdy%t_seg(jl)%i_width547 il_imax=1548 il_jmin=1549 il_jmax=il_jmax1-il_jmin1+1550 551 706 CASE('west') 552 707 … … 556 711 il_jmax1=td_bdy%t_seg(jl)%i_last 557 712 558 il_imin=1559 il_imax=td_bdy%t_seg(jl)%i_width560 il_jmin=1561 il_jmax=il_jmax1-il_jmin1+1562 563 713 END SELECT 564 714 565 !2 -read fine grid domain 566 tl_bathy1=td_bathy1 567 CALL iom_open(tl_bathy1) 568 569 !2-1 compute domain 715 ! -read fine grid domain 716 tl_bathy1=mpp_copy(td_bathy1) 717 718 ! compute domain 570 719 tl_dom1=dom_init( tl_bathy1, & 571 720 & il_imin1, il_imax1,& 572 & il_jmin1, il_jmax1 ) 573 574 !2-2 close file 575 CALL iom_close(tl_bathy1) 576 577 !2-3 read variables on domain (ugly way to do it, have to work on it) 578 !2-3-1 init mpp structure 579 tl_mppbathy1=mpp_init(tl_bathy1) 580 581 CALL file_clean(tl_bathy1) 582 583 !2-3-2 get processor to be used 584 CALL mpp_get_use( tl_mppbathy1, tl_dom1 ) 585 586 !2-3-3 open mpp files 587 CALL iom_mpp_open(tl_mppbathy1) 588 589 !2-3-4 read variable value on domain 590 tl_lon1=iom_mpp_read_var(tl_mppbathy1,'longitude',td_dom=tl_dom1) 591 tl_lat1=iom_mpp_read_var(tl_mppbathy1,'latitude' ,td_dom=tl_dom1) 592 593 !2-3-5 close mpp files 594 CALL iom_mpp_close(tl_mppbathy1) 595 596 !2-3-6 clean structure 597 CALL mpp_clean(tl_mppbathy1) 598 599 !3- get coarse grid indices 600 il_ind(:,:,:)=grid_get_coarse_index(td_bathy0, tl_lon1, tl_lat1, & 601 & id_rho=id_rho(:)) 602 603 il_imin0=il_ind(1,1,1) 604 il_imax0=il_ind(1,2,1) 605 606 il_jmin0=il_ind(2,1,1) 607 il_jmax0=il_ind(2,2,1) 608 609 il_offset(:,:)=il_ind(:,:,2) 610 611 !4- read coarse grid bathymetry on domain 612 tl_bathy0=td_bathy0 613 CALL iom_open(tl_bathy0) 614 615 !4-1 compute domain 721 & il_jmin1, il_jmax1,& 722 & TRIM(td_bdy%c_card)) 723 724 ! add extra band to fine grid domain (if possible) 725 ! to avoid dimension of one and so be able to compute offset 726 CALL dom_add_extra(tl_dom1, id_rho(jp_I), id_rho(jp_J)) 727 728 ! open mpp files over domain 729 CALL iom_dom_open(tl_bathy1, tl_dom1) 730 731 ! read variable value on domain 732 tl_lon1=iom_dom_read_var(tl_bathy1,'longitude',tl_dom1) 733 tl_lat1=iom_dom_read_var(tl_bathy1,'latitude' ,tl_dom1) 734 735 ! close mpp files 736 CALL iom_dom_close(tl_bathy1) 737 738 ! clean structure 739 CALL mpp_clean(tl_bathy1) 740 741 ! get coarse grid indices 742 il_ind(:,:)=grid_get_coarse_index(td_bathy0, tl_lon1, tl_lat1, & 743 & id_rho=id_rho(:)) 744 745 il_imin0=il_ind(1,1) 746 il_imax0=il_ind(1,2) 747 748 il_jmin0=il_ind(2,1) 749 il_jmax0=il_ind(2,2) 750 751 ! read coarse grid bathymetry on domain 752 tl_bathy0=mpp_copy(td_bathy0) 753 754 ! compute domain 616 755 tl_dom0=dom_init( tl_bathy0, & 617 756 & il_imin0, il_imax0,& 618 757 & il_jmin0, il_jmax0 ) 619 758 620 !4-2 close file 621 CALL iom_close(tl_bathy0) 622 623 !4-3 add extra band (if possible) to compute interpolation 759 il_offset(:,:)= grid_get_fine_offset(tl_bathy0, & 760 & il_imin0, il_jmin0,& 761 & il_imax0, il_jmax0,& 762 & tl_lon1%d_value(:,:,1,1), & 763 & tl_lat1%d_value(:,:,1,1), & 764 & id_rho=id_rho(:)) 765 766 ! clean 767 CALL var_clean(tl_lon1) 768 CALL var_clean(tl_lat1) 769 770 ! add extra band (if possible) to compute interpolation 624 771 CALL dom_add_extra(tl_dom0) 625 772 626 !4-4 read variables on domain (ugly way to do it, have to work on it) 627 !4-4-1 init mpp structure 628 tl_mppbathy0=mpp_init(tl_bathy0) 629 630 CALL file_clean(tl_bathy0) 631 632 !4-4-2 get processor to be used 633 CALL mpp_get_use( tl_mppbathy0, tl_dom0 ) 634 635 !4-4-3 open mpp files 636 CALL iom_mpp_open(tl_mppbathy0) 637 638 !4-4-4 read variable value on domain 639 tl_var0=iom_mpp_read_var(tl_mppbathy0,'Bathymetry',td_dom=tl_dom0) 640 641 !4-4-5 close mpp files 642 CALL iom_mpp_close(tl_mppbathy0) 643 644 !4-4-6 clean structure 645 CALL mpp_clean(tl_mppbathy0) 646 647 !5- interpolate variable 773 ! open mpp files over domain 774 CALL iom_dom_open(tl_bathy0, tl_dom0) 775 776 ! read variable value on domain 777 tl_var0=iom_dom_read_var(tl_bathy0,'Bathymetry',tl_dom0) 778 779 ! close mpp files 780 CALL iom_dom_close(tl_bathy0) 781 782 ! clean structure 783 CALL mpp_clean(tl_bathy0) 784 785 ! interpolate variable 648 786 CALL merge_bathy_interp( tl_var0, & 649 787 & id_rho(:), & 650 788 & il_offset(:,:) ) 651 789 652 ! 6-remove extraband added to domain790 ! remove extraband added to domain 653 791 CALL dom_del_extra( tl_var0, tl_dom0, id_rho(:) ) 654 792 655 ! 6-1remove extraband added to domain793 ! remove extraband added to domain 656 794 CALL dom_clean_extra( tl_dom0 ) 657 795 658 !7- fill refined table 659 !7-1 keep only useful point 660 ! interpolation could create more point than necessary 661 CALL boundary_clean_interp(tl_var0, td_bdy ) 662 663 ! use add request ???? 664 665 !7-2 fill refined table 796 ! remove extraband added to fine grid domain 797 CALL dom_del_extra( tl_var0, tl_dom1 ) 798 799 ! remove extraband added to fine grid domain 800 CALL dom_clean_extra( tl_dom1 ) 801 802 ! fill refined array 666 803 dd_refined( il_imin1:il_imax1, & 667 804 & il_jmin1:il_jmax1, & 668 805 & :,: )= tl_var0%d_value(:,:,:,:) 669 806 670 !8- compute weight function 807 ! clean 808 CALL var_clean(tl_var0) 809 810 ! compute weight function 671 811 ALLOCATE( dl_tmp1d(td_bdy%t_seg(jl)%i_width) ) 672 812 … … 678 818 679 819 ! compute weight on segment 680 dl_tmp1d(:)= 0.5 + 0.5*COS( (d g_pi*dl_tmp1d(:)) / &820 dl_tmp1d(:)= 0.5 + 0.5*COS( (dp_pi*dl_tmp1d(:)) / & 681 821 & (td_bdy%t_seg(jl)%i_width) ) 682 822 … … 694 834 695 835 ! compute weight on segment 696 dl_tmp1d(:)= 0.5 + 0.5*COS( (d g_pi*dl_tmp1d(:)) / &836 dl_tmp1d(:)= 0.5 + 0.5*COS( (dp_pi*dl_tmp1d(:)) / & 697 837 & (td_bdy%t_seg(jl)%i_width) ) 698 838 … … 710 850 711 851 ! compute weight on segment 712 dl_tmp1d(:)= 0.5 + 0.5*COS( (d g_pi*dl_tmp1d(:)) / &852 dl_tmp1d(:)= 0.5 + 0.5*COS( (dp_pi*dl_tmp1d(:)) / & 713 853 & (td_bdy%t_seg(jl)%i_width) ) 714 854 … … 726 866 727 867 ! compute weight on segment 728 dl_tmp1d(:)= 0.5 + 0.5*COS( (d g_pi*dl_tmp1d(:)) / &868 dl_tmp1d(:)= 0.5 + 0.5*COS( (dp_pi*dl_tmp1d(:)) / & 729 869 & (td_bdy%t_seg(jl)%i_width) ) 730 870 … … 740 880 DEALLOCATE( dl_tmp1d ) 741 881 742 ! 8-1 fill weight table882 ! fill weight array 743 883 ALLOCATE( dl_tmp2d( tl_dom1%t_dim(1)%i_len, & 744 884 & tl_dom1%t_dim(2)%i_len) ) … … 764 904 ENDIF 765 905 END SUBROUTINE merge_bathy_get_boundary 766 !> @endcode767 906 !------------------------------------------------------------------- 768 907 !> @brief 769 !> This subroutine 908 !> This subroutine interpolate variable. 770 909 !> 771 !> @details 910 !> @author J.Paul 911 !> @date November, 2013 - Initial Version 772 912 !> 773 !> @ author J.Paul774 !> - Nov, 2013- Initial Version775 !> 776 !> @param[in] 777 !> @ todo913 !> @param[inout] td_var variable structure 914 !> @param[in] id_rho array of refinment factor 915 !> @param[in] id_offset array of offset between fine and coarse grid 916 !> @param[in] id_iext i-direction size of extra bands (default=im_minext) 917 !> @param[in] id_jext j-direction size of extra bands (default=im_minext) 778 918 !------------------------------------------------------------------- 779 !> @code780 919 SUBROUTINE merge_bathy_interp( td_var, & 781 920 & id_rho, & … … 793 932 794 933 ! local variable 795 TYPE(TVAR) :: tl_var796 934 TYPE(TVAR) :: tl_mask 797 935 … … 803 941 ! loop indices 804 942 !---------------------------------------------------------------- 805 806 ! copy variable807 tl_var=td_var808 943 809 944 !WARNING: two extrabands are required for cubic interpolation … … 815 950 816 951 IF( il_iext < 2 .AND. td_var%c_interp(1) == 'cubic' )THEN 817 CALL logger_warn(" CREATE BATHY INTERP: at least extrapolation "//&952 CALL logger_warn("MERGE BATHY INTERP: at least extrapolation "//& 818 953 & "on two points are required with cubic interpolation ") 819 954 il_iext=2 … … 821 956 822 957 IF( il_jext < 2 .AND. td_var%c_interp(1) == 'cubic' )THEN 823 CALL logger_warn(" CREATE BATHY INTERP: at least extrapolation "//&958 CALL logger_warn("MERGE BATHY INTERP: at least extrapolation "//& 824 959 & "on two points are required with cubic interpolation ") 825 960 il_jext=2 826 961 ENDIF 827 962 828 ! 1-work on mask829 ! 1-1create mask830 ALLOCATE(bl_mask(t l_var%t_dim(1)%i_len, &831 & t l_var%t_dim(2)%i_len, &832 & t l_var%t_dim(3)%i_len, &833 & t l_var%t_dim(4)%i_len) )963 ! work on mask 964 ! create mask 965 ALLOCATE(bl_mask(td_var%t_dim(1)%i_len, & 966 & td_var%t_dim(2)%i_len, & 967 & td_var%t_dim(3)%i_len, & 968 & td_var%t_dim(4)%i_len) ) 834 969 835 970 bl_mask(:,:,:,:)=1 836 WHERE(t l_var%d_value(:,:,:,:)==tl_var%d_fill) bl_mask(:,:,:,:)=0837 838 SELECT CASE(TRIM(t l_var%c_point))971 WHERE(td_var%d_value(:,:,:,:)==td_var%d_fill) bl_mask(:,:,:,:)=0 972 973 SELECT CASE(TRIM(td_var%c_point)) 839 974 CASE DEFAULT ! 'T' 840 tl_mask=var_init('tmask', bl_mask(:,:,:,:), td_dim=tl_var%t_dim(:))841 CASE('U')842 tl_mask=var_init('umask', bl_mask(:,:,:,:), td_dim=tl_var%t_dim(:))843 CASE('V')844 tl_mask=var_init('vmask', bl_mask(:,:,:,:), td_dim=tl_var%t_dim(:))845 CASE('F')846 tl_mask=var_init('fmask', bl_mask(:,:,:,:), td_dim=tl_var%t_dim(:))975 tl_mask=var_init('tmask',bl_mask(:,:,:,:),td_dim=td_var%t_dim(:),& 976 & id_ew=td_var%i_ew ) 977 CASE('U','V','F') 978 CALL logger_fatal("MERGE BATHY INTERP: can not computed "//& 979 & "interpolation on "//TRIM(td_var%c_point)//& 980 & " grid point (variable "//TRIM(td_var%c_name)//& 981 & "). check namelist.") 847 982 END SELECT 848 983 849 984 DEALLOCATE(bl_mask) 850 985 851 ! 1-2interpolate mask986 ! interpolate mask 852 987 CALL interp_fill_value( tl_mask, id_rho(:), & 853 988 & id_offset=id_offset(:,:) ) 854 989 855 !2- work on variable 856 !2-0 add extraband 857 CALL extrap_add_extrabands(tl_var, il_iext, il_iext) 858 859 !2-1 extrapolate variable 860 CALL extrap_fill_value( tl_var, id_offset=id_offset(:,:), & 861 & id_rho=id_rho(:), & 862 & id_iext=il_iext, id_jext=il_jext ) 863 864 !2-2 interpolate Bathymetry 865 CALL interp_fill_value( tl_var, id_rho(:), & 990 ! work on variable 991 ! add extraband 992 CALL extrap_add_extrabands(td_var, il_iext, il_iext) 993 994 ! extrapolate variable 995 CALL extrap_fill_value( td_var ) 996 997 ! interpolate Bathymetry 998 CALL interp_fill_value( td_var, id_rho(:), & 866 999 & id_offset=id_offset(:,:) ) 867 1000 868 ! 2-3remove extraband869 CALL extrap_del_extrabands(t l_var, il_iext*id_rho(jp_I), il_jext*id_rho(jp_J))870 871 ! 2-4keep original mask1001 ! remove extraband 1002 CALL extrap_del_extrabands(td_var, il_iext*id_rho(jp_I), il_jext*id_rho(jp_J)) 1003 1004 ! keep original mask 872 1005 WHERE( tl_mask%d_value(:,:,:,:) == 0 ) 873 t l_var%d_value(:,:,:,:)=tl_var%d_fill1006 td_var%d_value(:,:,:,:)=td_var%d_fill 874 1007 END WHERE 875 1008 876 !3- save result877 td_var=tl_var878 879 ! clean variable structure880 CALL var_clean(tl_var)881 882 1009 END SUBROUTINE merge_bathy_interp 883 !> @endcode884 1010 END PROGRAM merge_bathy
Note: See TracChangeset
for help on using the changeset viewer.