|
XBeach
|
00001 module beachwizard_module 00002 implicit none 00003 save 00004 private 00005 public bwinit, assim, assim_update 00006 type beachwiz 00007 00008 real*8, dimension(1000,8) :: timesnew !time in minutes when remote sensed sources apply - eight sources max 00009 integer, dimension(8) :: ntimesnew = 0 !number of maps per source 00010 integer :: nflagsigpr = 0 !flags 00011 integer :: nflagsetobs = 0 00012 real*8, dimension(1000,8) :: timvalnew !duration in minutes of each map 00013 character(232), dimension(1000,8) :: fname !filename of the map 00014 character(232), dimension(1000,8) :: fnameerr !filename of the uncertainties in the map 00015 character(232), dimension(:), 00016 pointer :: fnamim 00017 character(232), dimension(:), 00018 pointer :: fnamzb 00019 character(232), dimension(:), 00020 pointer :: fnamcx 00021 character(232), dimension(:), 00022 pointer :: fnamsh 00023 integer, dimension(8) :: status = 0 !status of the map (1=read, 0=not read/does not exist) 00024 integer, dimension(8) :: statuserr = 0 !status of the error map 00025 integer, dimension(8) :: itimnew !counter of the time "timesnew" 00026 integer, dimension(8) :: itimrdnew !flag to indicate time has been read 00027 integer :: itim 00028 integer :: itimrd 00029 integer :: itim2 00030 integer :: itimrd2 00031 integer :: itim3 00032 integer :: itimrd3 00033 integer :: itim4 00034 integer :: itimrd4 00035 real*8, dimension(:,:), pointer :: fobs1 !value of the observed quantity in the coordinate system of the file 00036 real*8, dimension(:,:), pointer :: fobs1err !same, for error file 00037 real*8, dimension(:,:), pointer :: disobs 00038 real*8, dimension(:,:), pointer :: zbobs1 00039 real*8, dimension(:,:), pointer :: shobs1 00040 ! real*8,dimension(:,:),allocatable :: dobs !observed dissipation 00041 real*8,dimension(:,:),allocatable :: fobs 00042 real*8,dimension(:,:),allocatable :: cobs !observed celerity 00043 real*8,dimension(:,:),allocatable :: ccom !computed celerity 00044 real*8,dimension(:,:),allocatable :: kobs !observed wave number 00045 ! real*8,dimension(:,:),allocatable :: shobs !observed shoreline 00046 ! real*8,dimension(:,:),allocatable :: zbobs !observed bathymetry 00047 real*8,dimension(:,:),allocatable :: dcmdo !computed dissipation minus observed dissipation 00048 real*8,dimension(:,:),allocatable :: ccmco !computed celerity minus observed celerity 00049 real*8,dimension(:,:),allocatable :: scmso !computed shoreline minus observed shoreline 00050 real*8,dimension(:,:),allocatable :: dassim 00051 real*8,dimension(:,:),allocatable :: sigD !standard dev. of dissipation 00052 real*8,dimension(:,:),allocatable :: sigC !standard dev. of celerity 00053 real*8,dimension(:,:),allocatable :: sigK !standard dev. of wave number 00054 real*8,dimension(:,:),allocatable :: sigZ !standard dev. of bathymetry 00055 real*8,dimension(:,:),allocatable :: sigS !standard dev. of shoreline 00056 ! real*8,dimension(:,:),allocatable :: sig2prior 00057 real*8 :: xll !origin of the coordinate system used in the map 00058 real*8 :: yll !idem 00059 real*8 :: dx !dx in that c.s. 00060 real*8 :: dy !dy in that c.s. 00061 real*8 :: angle !angle in degrees 00062 real*8 :: tradar !period used in the wavenumber or celerity file 00063 integer :: nx !number of points in the x-direction in the c.s. of the map 00064 integer :: ny !same in y 00065 integer :: ix !counter 00066 integer :: iy !counter 00067 real*8 :: sigDdef !default value of the dissipation uncertainty 00068 real*8 :: sigCdef !default value of the celerity uncertainty 00069 real*8 :: sigZdef !default value of the bathymetry uncertainty 00070 real*8 :: sigSdef !default value of the shoreline uncertainty 00071 00072 00073 end type beachwiz 00074 00075 ! Type of beachwizard, stored local. Please store persistent info in s, par. 00076 type(beachwiz) :: bw 00077 00078 contains 00079 00080 subroutine bwinit(s) 00081 use spaceparamsdef, only: spacepars 00082 use xmpi_module, only: xmaster 00083 00084 00085 implicit none 00086 00087 type(spacepars), target :: s 00088 00089 if(.not. xmaster) return 00090 00091 if(.true.) then 00092 allocate(s%dobs(1:s%nx+1,1:s%ny+1)) 00093 allocate(s%zbobs(1:s%nx+1,1:s%ny+1)) 00094 allocate(s%sig2prior(1:s%nx+1,1:s%ny+1)) 00095 allocate(s%shobs(1:s%nx+1,1:s%ny+1)) 00096 allocate(s%bwalpha(1:s%nx+1,1:s%ny+1)) 00097 allocate(s%dcmdo(1:s%nx+1,1:s%ny+1)) 00098 allocate(s%dassim(1:s%nx+1,1:s%ny+1)) 00099 allocate(s%cobs(1:s%nx+1,1:s%ny+1)) 00100 endif 00101 00102 00103 end subroutine bwinit 00104 00105 subroutine assim(s,par) 00106 00107 use params, only: parameters 00108 use spaceparamsdef, only: spacepars 00109 00110 00111 00112 implicit none 00113 00114 type(spacepars), target :: s 00115 type(parameters) :: par 00116 00117 character*80 :: infofile 00118 integer :: im,in 00119 integer :: srcnr 00120 real*8, external :: readkey_dbl 00121 integer, external :: readkey_int 00122 logical :: exists 00123 real*8, allocatable, dimension(:,:) :: h1 00124 00125 integer :: jj 00126 real*8 :: om 00127 real*8, allocatable, dimension(:,:) :: ome2 00128 real*8, allocatable, dimension(:,:) :: num 00129 real*8, allocatable, dimension(:,:) :: den 00130 real*8, allocatable, dimension(:,:) :: rk 00131 real*8, allocatable, dimension(:,:) :: Hrms 00132 real*8, allocatable, dimension(:,:) :: E 00133 real*8, allocatable, dimension(:,:) :: hh 00134 00135 real*8, parameter :: a1 = 5.060219360721177D-01 00136 real*8, parameter :: a2 = 2.663457535068147D-01 00137 real*8, parameter :: a3 = 1.108728659243231D-01 00138 real*8, parameter :: a4 = 4.197392043833136D-02 00139 real*8, parameter :: a5 = 8.670877524768146D-03 00140 real*8, parameter :: a6 = 4.890806291366061D-03 00141 real*8, parameter :: b1 = 1.727544632667079D-01 00142 real*8, parameter :: b2 = 1.191224998569728D-01 00143 real*8, parameter :: b3 = 4.165097693766726D-02 00144 real*8, parameter :: b4 = 8.674993032204639D-03 00145 00146 00147 00148 if (.not. allocated(bw%fobs)) then 00149 ! allocate(s%dobs(1:s%nx+1,1:s%ny+1)) 00150 allocate(bw%fobs(1:s%nx+1,1:s%ny+1)) 00151 allocate(bw%cobs(1:s%nx+1,1:s%ny+1)) 00152 allocate(bw%ccom(1:s%nx+1,1:s%ny+1)) 00153 allocate(bw%kobs(1:s%nx+1,1:s%ny+1)) 00154 ! allocate(s%shobs(1:s%nx+1,1:s%ny+1)) 00155 ! allocate(s%zbobs(1:s%nx+1,1:s%ny+1)) 00156 allocate(bw%dcmdo(1:s%nx+1,1:s%ny+1)) 00157 allocate(bw%ccmco(1:s%nx+1,1:s%ny+1)) 00158 allocate(bw%scmso(1:s%nx+1,1:s%ny+1)) 00159 allocate(bw%dassim(1:s%nx+1,1:s%ny+1)) 00160 allocate(bw%sigD(1:s%nx+1,1:s%ny+1)) 00161 allocate(bw%sigC(1:s%nx+1,1:s%ny+1)) 00162 allocate(bw%sigK(1:s%nx+1,1:s%ny+1)) 00163 allocate(bw%sigZ(1:s%nx+1,1:s%ny+1)) 00164 allocate(bw%sigS(1:s%nx+1,1:s%ny+1)) 00165 ! allocate(s%sig2prior(1:s%nx+1,1:s%ny+1)) 00166 end if 00167 00168 allocate(h1(1:s%nx+1,1:s%ny+1)) 00169 allocate(hh(1:s%nx+1,1:s%ny+1)) 00170 allocate(ome2(1:s%nx+1,1:s%ny+1)) 00171 allocate(num(1:s%nx+1,1:s%ny+1)) 00172 allocate(rk(1:s%nx+1,1:s%ny+1)) 00173 allocate(E(1:s%nx+1,1:s%ny+1)) 00174 allocate(Hrms(1:s%nx+1,1:s%ny+1)) 00175 00176 ! wwvv somebody forgot: 00177 allocate(den(s%nx+1,s%ny+1)) 00178 00179 bw%dcmdo=0. 00180 bw%ccmco=0. 00181 bw%scmso=0. 00182 00183 if (bw%nflagsetobs==0) then !do this only once. 00184 00185 s%dobs=-999. 00186 bw%cobs=-999. 00187 s%shobs=-999. 00188 s%zbobs = -999. 00189 00190 bw%nflagsetobs = 1 00191 00192 end if 00193 00194 if (bw%nflagsigpr==0) then !read prior uncertainty in terms of variance 00195 00196 inquire (file='sigpr.dep', exist=exists) 00197 if (exists) then 00198 open(444,file='sigpr.dep',err=999) 00199 do im = 1,s%ny+1 00200 read(444,*)(s%sig2prior(in,im),in =1,s%nx+1) !Leo: nx and ny --> nx+1 ny+1 00201 enddo 00202 close(444) 00203 else 00204 s%sig2prior = 1**2 !if the file does not exist, set to hardcoded 1**2 00205 end if 00206 00207 inquire (file='bwdefaults.inp', exist=exists) !see if defaults file exists 00208 if (exists) then 00209 open(444,file='bwdefaults.inp',err=999) 00210 read(444,*) bw%sigDdef 00211 read(444,*) bw%sigCdef 00212 read(444,*) bw%sigZdef 00213 read(444,*) bw%sigSdef 00214 close(444) 00215 bw%sigD=bw%sigDdef 00216 bw%sigC=bw%sigCdef 00217 bw%sigZ=bw%sigZdef 00218 bw%sigS=bw%sigSdef 00219 else 00220 bw%sigCdef=1. !set to hardcoded defaults 00221 bw%sigDdef=20. 00222 bw%sigZdef=0.5 00223 bw%sigSdef=0.5 00224 end if 00225 bw%nflagsigpr=1 00226 00227 end if 00228 00229 999 continue 00230 00231 00232 inquire (file='imageinfoimap',exist=exists) !read file for dissipation 00233 if (exists) then 00234 infofile='imageinfoimap' 00235 srcnr = 1 00236 call assim_rd(bw, s, par, infofile, srcnr) 00237 if (bw%status(srcnr) ==1) then 00238 do in=1,s%nx 00239 do im=1,s%ny 00240 if ((s%wetu(in,im)>0).and.(s%dobs(in,im)>-5.)) then 00241 bw%dcmdo(in,im) = s%Dr(in,im) - s%dobs(in,im) 00242 end if 00243 end do 00244 end do 00245 end if 00246 if (bw%statuserr(srcnr) ==0) then 00247 h1 = max(s%hh+par%delta*s%H,par%hmin) 00248 !! bw%sigD=bw%sigDdef+(999.-bw%sigDdef)*(0.5+0.5*tanh((s%dobs/h1-bw%sigDdef)/10.)) ! change 50 --> bw%sigDdef 00249 bw%sigD = bw%sigDdef !!using only default dissipation uncertainty for now... 00250 00251 end if 00252 end if 00253 00254 inquire (file='imageinfozb',exist=exists) !read file for bathymetry 00255 if (exists) then 00256 infofile='imageinfozb' 00257 srcnr = 2 00258 call assim_rd(bw, s, par, infofile, srcnr) 00259 if (bw%statuserr(srcnr) ==0) then 00260 bw%sigZ=bw%sigZDef 00261 end if 00262 end if 00263 00264 inquire (file='imageinfocx',exist=exists) !read file for celerity 00265 if (exists) then 00266 infofile='imageinfocx' 00267 srcnr = 3 00268 call assim_rd(bw, s, par, infofile, srcnr) 00269 if (bw%status(srcnr) ==1) then 00270 do in=1,s%nx+1 00271 do im=1,s%ny+1 00272 om = 2.*3.1415/bw%tradar 00273 ome2(in,im) = om**2*s%hh(in,im)/par%g 00274 num(in,im) = 1.0 + & 00275 & ome2(in,im)*(a1 + ome2(in,im)*(a2 + ome2(in,im)*(a3 + ome2(in,im) & 00276 & *(a4 + ome2(in,im)*(a5 + ome2(in,im)*a6))))) 00277 den(in,im) = 1.0 + ome2(in,im)*(b1 + ome2(in,im)*(b2 + ome2(in,im)*(b3 + & 00278 & ome2(in,im)*(b4 + ome2(in,im)*a6)))) 00279 rk(in,im) = sqrt(ome2(in,im)*num(in,im)/den(in,im))/s%hh(in,im) 00280 Hrms(in,im) = max(0.01d0,sqrt(8*max(0.01d0,s%E(in,im))/par%rho/par%g)) 00281 bw%ccom(in,im) = sqrt(par%g/rk(in,im)*tanh(rk(in,im)*s%hh(in,im)+Hrms(in,im))) 00282 if ((s%wetu(in,im)>0).and.(bw%cobs(in,im)>-990.)) then 00283 bw%ccmco(in,im) = bw%ccom(in,im) - bw%cobs(in,im) 00284 end if 00285 enddo 00286 enddo 00287 end if 00288 if (bw%statuserr(srcnr) ==0) then 00289 bw%sigC=bw%sigCdef 00290 endif 00291 end if 00292 00293 inquire (file='imageinfokabs',exist=exists) !do the same as above for wavenumbers 00294 if (exists) then 00295 infofile='imageinfokabs' 00296 srcnr = 4 00297 call assim_rd(bw, s, par, infofile, srcnr) 00298 if (bw%status(srcnr) ==1) then 00299 do in = 1,s%nx+1 00300 do im = 1,s%ny+1 00301 om = 2.*3.1415/bw%tradar 00302 ome2(in,im) = om**2*s%hh(in,im)/par%g 00303 num(in,im) = 1.0 + & 00304 & ome2(in,im)*(a1 + ome2(in,im)*(a2 + ome2(in,im)*(a3 + ome2(in,im) & 00305 & *(a4 + ome2(in,im)*(a5 + ome2(in,im)*a6))))) 00306 den(in,im) = 1.0 + ome2(in,im)*(b1 + ome2(in,im)*(b2 + ome2(in,im)*(b3 + & 00307 & ome2(in,im)*(b4 + ome2(in,im)*a6)))) 00308 rk(in,im) = sqrt(ome2(in,im)*num(in,im)/den(in,im))/s%hh(in,im) 00309 Hrms(in,im) = max(0.01d0,sqrt(8*max(0.01d0,s%E(in,im))/par%rho/par%g)) 00310 bw%ccom(in,im) = sqrt(par%g/rk(in,im)*tanh(rk(in,im)*s%hh(in,im)+Hrms(in,im))) 00311 bw%cobs(in,im) = om/max(bw%kobs(in,im),0.01d0) ! here we revert back to celerities. 00312 if (bw%kobs(in,im)<-990.) then 00313 bw%cobs(in,im)=-999. 00314 end if 00315 if ((s%wetu(in,im)>0).and.(bw%kobs(in,im)>-990.)) then 00316 bw%ccmco(in,im) = bw%ccom(in,im) - bw%cobs(in,im) 00317 end if 00318 end do 00319 end do 00320 end if 00321 if (bw%statuserr(srcnr) ==1) then 00322 do in = 1,s%nx+1 00323 do im = 1,s%ny+1 00324 om = 2.*3.1415/bw%tradar 00325 ome2(in,im) = om**2*s%hh(in,im)/par%g 00326 num(in,im) = 1.0 + & 00327 & ome2(in,im)*(a1 + ome2(in,im)*(a2 + ome2(in,im)*(a3 + ome2(in,im) & 00328 & *(a4 + ome2(in,im)*(a5 + ome2(in,im)*a6))))) 00329 den(in,im) = 1.0 + ome2(in,im)*(b1 + ome2(in,im)*(b2 + ome2(in,im)*(b3 + & 00330 & ome2(in,im)*(b4 + ome2(in,im)*a6)))) 00331 rk(in,im) = sqrt(ome2(in,im)*num(in,im)/den(in,im))/s%hh(in,im) 00332 bw%sigC(in,im) = om/rk(in,im)*(bw%sigK(in,im)/(rk(in,im)+bw%sigK(in,im))) ! express sigK in terms of sigK, 00333 ! from c+sigC = om/(k+sigK) 00334 end do 00335 end do 00336 end if 00337 if (bw%statuserr(srcnr) ==0) then 00338 bw%sigC = bw%sigCdef 00339 end if 00340 end if 00341 00342 inquire (file='imageinfoibathy',exist=exists) !read file for shorelines 00343 if (exists) then 00344 infofile='imageinfoibathy' 00345 srcnr = 5 00346 call assim_rd(bw, s, par, infofile, srcnr) 00347 if (bw%status(srcnr) ==1) then 00348 00349 do im=1,s%ny 00350 do in = 1,s%nx 00351 if (abs(s%shobs(in,im)).le.5.) then 00352 ! bw%scmso(in,im) = s%hh(in,im) + s%shobs(in,im) !!plus? 00353 bw%scmso(in,im) = s%zb(in,im) - s%shobs(in,im) 00354 end if 00355 end do 00356 end do 00357 00358 00359 do jj=1,2 00360 do in = 2,s%nx 00361 do im = 2,s%ny 00362 bw%scmso(in,im) = (bw%scmso(in-1,im-1)+bw%scmso(in,im-1)+bw%scmso(in+1,im-1)+ & 00363 & bw%scmso(in-1,im)+ bw%scmso(in,im)+ bw%scmso(in+1,im)+ & 00364 & bw%scmso(in-1,im+1)+bw%scmso(in,im+1)+bw%scmso(in+1,im+1))/9. 00365 end do 00366 end do 00367 end do 00368 00369 end if !if status 00370 if (bw%statuserr(srcnr) ==0) then 00371 bw%sigS=bw%sigSdef 00372 end if 00373 end if !if exists 00374 00375 !expose variables to output module 00376 00377 s%dcmdo = bw%dcmdo 00378 00379 call comp_depchg(bw, s, par) 00380 00381 end subroutine assim 00382 00383 subroutine assim_rd(bw, s, par, infofile, srcnr) 00384 00385 00386 00387 use params, only: parameters 00388 use spaceparamsdef, only: spacepars 00389 implicit none 00390 ! Global variables 00391 type(spacepars), target :: s 00392 type(parameters) :: par 00393 type(beachwiz) :: bw 00394 character*80 :: infofile 00395 integer :: im,in,iy,ix 00396 integer :: srcnr 00397 00398 ! Local variables 00399 real*8 :: degrad 00400 real*8 :: a 00401 real*8 :: b 00402 real*8 :: w11 00403 real*8 :: w21 00404 real*8 :: w12 00405 real*8 :: w22 00406 real*8 :: x1 00407 real*8 :: y1 00408 real*8 :: xs 00409 real*8 :: ys 00410 real*8 :: cs 00411 real*8 :: sn 00412 real*8 :: x1max 00413 real*8 :: y1max 00414 integer :: ixp1 00415 integer :: iyp1 00416 integer :: i 00417 integer :: istat 00418 logical :: exists 00419 00420 ! Actual time based on timhr relative to itdate 00421 ! 00422 00423 00424 if (bw%ntimesnew(srcnr)==0) then 00425 inquire (file=infofile,exist=exists) 00426 00427 00428 if (exists) then 00429 00430 open(31,file=infofile,err=999) 00431 read(31,*)bw%ntimesnew(srcnr) 00432 00433 do i=1,bw%ntimesnew(srcnr) 00434 read(31,*)bw%timesnew(i,srcnr),bw%fname(i,srcnr),bw%fnameerr(i,srcnr), & 00435 & bw%timvalnew(i,srcnr) 00436 write(*,*) bw%timesnew(i,srcnr), bw%fname(i,srcnr), bw%fnameerr(i,srcnr), & 00437 & bw%timvalnew(i,srcnr) 00438 enddo 00439 !pause 00440 close(31) 00441 bw%timesnew(bw%ntimesnew(srcnr)+1:1000,srcnr)=1.e10 00442 00443 00444 else 00445 stop ' Argus assimilation active and no imageinfoimap file' 00446 endif 00447 bw%itimnew(srcnr)=1 00448 bw%itimrdnew(srcnr)=0 00449 00450 endif 00451 00452 bw%status(srcnr) = 0 00453 00454 00455 ! Check if time matches time of an image 00456 if (par%t<60.*bw%timesnew(1,srcnr)) then 00457 ! 00458 ! Time less than time first image; leave dcmdo at 0 00459 else 00460 10 if (par%t>=60.*bw%timesnew(bw%itimnew(srcnr)+1,srcnr)) then 00461 bw%itimnew(srcnr)=bw%itimnew(srcnr)+1 00462 goto 10 00463 endif 00464 ! 00465 00466 if (par%t-60.*bw%timesnew(bw%itimnew(srcnr),srcnr)<=60*bw%timvalnew(bw%itimnew(srcnr),srcnr)) then 00467 ! Read image if not read before 00468 ! write(*,*) 'smaller than timval',timmin,timesnew(itimnew(srcnr),srcnr), & 00469 !& timvalnew(itimnew(srcnr),srcnr) 00470 if (bw%itimnew(srcnr)/=bw%itimrdnew(srcnr)) then 00471 ! 00472 ! Read header file that defines argus grid 00473 inquire(file=bw%fname(bw%itimnew(srcnr),srcnr),exist=exists) 00474 if (exists) then 00475 open(31,file=bw%fname(bw%itimnew(srcnr),srcnr),err=999) 00476 read(31,*,err=999)bw%xll 00477 read(31,*,err=999)bw%yll 00478 read(31,*,err=999)bw%dx 00479 read(31,*,err=999)bw%dy 00480 read(31,*,err=999)bw%nx 00481 read(31,*,err=999)bw%ny 00482 read(31,*,err=999)bw%angle 00483 00484 if (associated(bw%fobs1)) deallocate (bw%fobs1, STAT = istat) 00485 allocate (bw%fobs1(bw%nx,bw%ny)) 00486 do iy=1,bw%ny 00487 read(31,*)(bw%fobs1(ix,iy),ix=1,bw%nx) 00488 enddo 00489 close(31) 00490 bw%itimrdnew(srcnr)=bw%itimnew(srcnr) 00491 else 00492 goto 999 00493 endif 00494 ! inquire(file=fnameerr(itimnew(srcnr),srcnr),exist=exists) 00495 ! write(*,*) fnameerr(itimnew(srcnr),srcnr) 00496 !pause 00497 ! if (exists) then 00498 ! 31=NEWLUN(GDP) 00499 ! open(31,file=fnameerr(itimnew(srcnr),srcnr),err=999) 00500 ! read(31,*,err=999)xll 00501 ! read(31,*,err=999)yll 00502 ! read(31,*,err=999)dx 00503 ! read(31,*,err=999)dy 00504 ! read(31,*,err=999)nx 00505 ! read(31,*,err=999)ny 00506 ! read(31,*,err=999)angle 00507 ! write(*,*) xll, yll, dx, dy, nx, ny, angle 00508 !pause 00509 ! if (associated(fobs1err)) deallocate (fobs1err, STAT = istat) 00510 ! allocate (fobs1err(nx,ny)) 00511 ! do iy=1,ny 00512 ! read(31,*)(fobs1err(ix,iy),ix=1,nx) 00513 ! enddo 00514 ! close(31) 00515 ! statuserr(srcnr)=1 00516 ! write(*,*) 'here read ', statuserr(srcnr), srcnr 00517 !pause 00518 ! else 00519 ! fobs1err=-999 00520 ! statuserr(srcnr)=0 00521 ! write(*,*) 'here not read ', statuserr(srcnr), srcnr 00522 ! !pause 00523 ! endif 00524 degrad=atan(1.)/45. 00525 cs=cos(bw%angle*degrad) 00526 sn=sin(bw%angle*degrad) 00527 x1max=(bw%nx-2)*bw%dx 00528 y1max=(bw%ny-2)*bw%dy 00529 00530 do in=1,s%nx+1 !Leo and Roberto changed nx --> nx+1 00531 do im=1,s%ny+1 !Leo and Roberto changed ny --> ny+1 00532 bw%fobs(in,im)=-999.d0 !Changed by Ap 31/5 00533 ! fobserr(n,m) = 999.d0 ! changed by Ap 31/5 00534 ! if (kfs(n,m)>0) then !Ap 00535 xs = s%xz(in,im) - bw%xll ! Jaap changed x --> xz 00536 ys = s%yz(in,im) - bw%yll 00537 x1 = xs*cs + ys*sn 00538 y1 =-xs*sn + ys*cs 00539 x1 = min(max(x1,bw%dx),x1max) 00540 y1 = min(max(y1,bw%dy),y1max) 00541 ix = int(x1/bw%dx)+1 00542 iy = int(y1/bw%dy)+1 00543 ixp1 = max(min(ix+1,bw%nx-1),2) 00544 iyp1 = max(min(iy+1,bw%ny-1),2) 00545 ix = max(min(ix,bw%nx-1),2) 00546 iy = max(min(iy,bw%ny-1),2) 00547 a = mod(x1,bw%dx)/bw%dx 00548 b = mod(y1,bw%dy)/bw%dy 00549 w11 = (1.-b)*(1.-a) 00550 w21 = (1.-b)* a 00551 w12 = b *(1.-a) 00552 w22 = b * a 00553 ! if (kfs(n,m)>0) then 00554 bw%fobs(in,im) = w11*bw%fobs1(ix ,iy ) + & 00555 & w21*bw%fobs1(ixp1,iy ) + & 00556 & w12*bw%fobs1(ix ,iyp1) + & 00557 & w22*bw%fobs1(ixp1,iyp1) 00558 00559 ! endif 00560 ! if (statuserr(srcnr)==1) then 00561 ! 00562 ! fobserr(n,m) = w11*fobs1err(ix ,iy ) + & 00563 ! & w21*fobs1err(ixp1,iy ) + & 00564 ! & w12*fobs1err(ix ,iyp1) + & 00565 ! & w22*fobs1err(ixp1,iyp1) 00566 ! end if 00567 00568 enddo 00569 enddo 00570 if (srcnr==1) s%dobs=bw%fobs 00571 if (srcnr==2) s%zbobs=bw%fobs 00572 if (srcnr==3) bw%cobs=bw%fobs 00573 if (srcnr==4) bw%cobs=bw%fobs 00574 if (srcnr==5) s%shobs=bw%fobs 00575 00576 endif 00577 00578 bw%status(srcnr)=1 00579 00580 else 00581 bw%fobs=0. 00582 endif 00583 endif 00584 00585 999 continue 00586 00587 end subroutine assim_rd 00588 00589 subroutine assim_update(s, par) 00590 use params, only: parameters 00591 use spaceparamsdef, only: spacepars 00592 type(spacepars), target :: s 00593 type(parameters) :: par 00594 00595 ! These are doing the same thing now 00596 if (par%bchwiz .eq. 1) then 00597 ! s%dzbdt(i,j) = s%dzbdt(i,j)-bw%dassim(i,j)/par%dt 00598 ! Times dt*morfac gives 00599 s%zb = s%zb - bw%dassim 00600 endif 00601 if (par%bchwiz .eq. 2) then 00602 s%zb = s%zb - bw%dassim !if we have only beachwizard 00603 end if 00604 end subroutine assim_update 00605 subroutine comp_depchg(bw, s, par) 00606 00607 use params, only: parameters 00608 use spaceparamsdef, only: spacepars 00609 implicit none 00610 ! Global variables 00611 type(spacepars), target :: s 00612 type(parameters) :: par 00613 type(beachwiz) :: bw 00614 00615 00616 ! Local variables 00617 integer :: in,i 00618 integer :: im,j 00619 integer :: ind 00620 real*8 :: Hrmsmax 00621 real*8 :: mxdDdh 00622 real*8 :: alp_as 00623 real*8 :: Nassim 00624 real*8, allocatable, dimension(:,:) :: errD 00625 real*8, allocatable, dimension(:,:) :: errC 00626 real*8, allocatable, dimension(:,:) :: errS 00627 real*8, allocatable, dimension(:,:) :: Hb 00628 real*8, allocatable, dimension(:,:) :: kh 00629 real*8, allocatable, dimension(:,:) :: gambal 00630 real*8, allocatable, dimension(:,:) :: Ga 00631 real*8, allocatable, dimension(:,:) :: dDdGa 00632 real*8, allocatable, dimension(:,:) :: dGadHb 00633 real*8, allocatable, dimension(:,:) :: dHbdh 00634 real*8, allocatable, dimension(:,:) :: h1 00635 real*8, allocatable, dimension(:,:) :: dDdh 00636 real*8, allocatable, dimension(:,:) :: sig2obs 00637 ! real*8, allocatable, dimension(:,:) :: alpha 00638 real*8, allocatable, dimension(:,:) :: Hrms 00639 real*8, allocatable, dimension(:,:) :: alphafac 00640 ! 00641 real*8 :: backdis,disfac 00642 integer :: index 00643 real*8, allocatable, dimension(:,:) :: h2 00644 if (.not.(allocated(errD))) then 00645 allocate (errD(1:s%nx+1,1:s%ny+1)) 00646 allocate (errC(1:s%nx+1,1:s%ny+1)) 00647 allocate (errS(1:s%nx+1,1:s%ny+1)) 00648 allocate (Hb(1:s%nx+1,1:s%ny+1)) 00649 allocate (kh(1:s%nx+1,1:s%ny+1)) 00650 allocate (gambal(1:s%nx+1,1:s%ny+1)) 00651 allocate (Ga(1:s%nx+1,1:s%ny+1)) 00652 allocate (dGadHb(1:s%nx+1,1:s%ny+1)) 00653 allocate (dHbdh(1:s%nx+1,1:s%ny+1)) 00654 allocate (dDdGa(1:s%nx+1,1:s%ny+1)) 00655 allocate (h1(1:s%nx+1,1:s%ny+1)) 00656 allocate (dDdh(1:s%nx+1,1:s%ny+1)) 00657 allocate (sig2obs(1:s%nx+1,1:s%ny+1)) 00658 ! allocate (alpha(1:s%nx+1,1:s%ny+1)) 00659 allocate (Hrms(1:s%nx+1,1:s%ny+1)) 00660 allocate (h2(1:s%nx+1,1:s%ny+1)) 00661 allocate (alphafac(1:s%nx+1,1:s%ny+1)) 00662 end if 00663 00664 00665 ! general 00666 Hrms = max(s%H,0.01d0); 00667 h2 = 0.d0 ! modified wave length, initially set to L1 00668 do j = 2,s%ny 00669 do i = 2,s%nx+1 00670 index = i ! start index 00671 backdis = 0.d0 ! relative distance backward 00672 do while (backdis<1.d0) 00673 ! disfac = s%dsc(index,j)/(par%facsd*s%L1(index,j)) 00674 ! use average wavelength over distance dsc 00675 disfac = s%dsc(index,j)/(0.5d0*(s%L1(index,j)+s%L1(max(index-1,1),j))) 00676 disfac = min(disfac,1.d0-backdis) 00677 00678 h2(i,j) = h2(i,j) + disfac*s%hh(index,j) 00679 backdis = backdis+disfac 00680 00681 index = max(index-1,1) 00682 enddo 00683 enddo 00684 enddo 00685 h2(:,1)=h2(:,2) 00686 h2(:,s%ny+1)=h2(:,s%ny) 00687 h1 = max(s%hh+par%delta*Hrms,par%hmin) !total water depth 00688 kh = min(10.0, s%k*h1) 00689 00690 ! ! The following equation numbers and variables follow the Beach Wizard Paper 00691 !! (Van Dongeren et al., 2008) 00692 00693 ! derivatives for dissipation assimilation 00694 gambal = 0.29+0.76*kh 00695 Hb = 0.88/s%k*tanh(gambal*kh/0.88) ! eq. (A4) 00696 Hrmsmax = maxval(Hrms) 00697 Ga = (Hb/Hrms)**2 00698 00699 dDdGa = -0.25*par%rho*par%g/par%Trep*(Hrms**2)*Ga*exp(-Ga) !eq. (A5) 00700 dGadHb = 2*Hb/(Hrms**2) 00701 00702 dHbdh = ((0.29+2*0.76*kh)*(1.-(kh/(sinh(kh)*cosh(kh)+kh)))) & ! eq. (A6) 00703 & /(cosh((0.29*kh+0.76*kh**2)/0.88)**2) + & 00704 & 0.88*tanh((0.29*kh+0.76*kh**2)/0.88)/(sinh(kh)*cosh(kh)+kh) 00705 00706 dDdh = dDdGa*dGadHb*dHbdh ! eq. (A3), derivative of dissipation with respect to water depth 00707 00708 00709 mxdDdh = maxval(dDdh**2) 00710 00711 00712 !! bw%sigD = bw%sigD+(100.d0-bw%sigD)*(1-tanh((3.7*Hrmsmax/h1)**20)) !measurement error for dissipatiopn, used in eq. (6) 00713 !bw%sigD above is commented out for now to use default value of measurement error (Roberto) 00714 00715 00716 bw%sigD = 0.15*maxval(s%dobs)!! dissipation measurement error equal to a percentage (15%) of the maximum observed Dissipation. 00717 00718 !! do in=1,s%nx+1 00719 !! do im=1,s%ny+1 !! here the measurement error is spatially varied, increasing for shallow waters where sandy areas or persistent foam may be. (Roberto) 00720 !! if (bw%sigD(in,im)/(h1(in,im)+1.e-16)>bw%sigD(in,im)) bw%sigD(in,im)=bw%sigD(in,im)/(h1(in,im)+1.e-16) 00721 !! end do 00722 !! end do 00723 00724 alp_as = 0.005 00725 errD = ((bw%dcmdo**2+bw%sigD**2)/(dDdh**2+alp_as*mxdDdh+1.e-16)) !eq.(6) uncertainty in observed data (Dissipation) 00726 00727 do in=1,s%nx+1 00728 do im=1,s%ny+1 00729 if (s%dobs(in,im)<-990) errD(in,im)=999 00730 end do 00731 end do 00732 00733 errC=999 !!! for now, if we want c as a source this needs the change to 00734 ! errC(nm) = ((ccmco(nm)**2+sigC(nm)**2)/(dcdh(nm)**2+1.e-16)) without indices 00735 do in=1,s%nx+1 00736 do im=1,s%ny+1 00737 if (bw%cobs(in,im)<-990) errC(in,im)=999 00738 end do 00739 end do 00740 00741 ! for output 00742 s%cobs = bw%cobs 00743 00744 errS = ((bw%scmso**2+bw%sigS**2)/(1.d0)) !eq.(6) uncertainty in observed data (ibathy) 00745 do in=1,s%nx+1 00746 do im=1,s%ny+1 00747 if (s%shobs(in,im)<-990) errS(in,im)=999 00748 end do 00749 end do 00750 00751 00752 00753 Nassim = par%tstop/par%dt !! changed par%tstop/par%dt -->par%tstop/par%wavint 00754 00755 sig2obs = 1.0d0/(1./errD+1./errS) 00756 s%bwalpha = s%sig2prior/(s%sig2prior+Nassim*sig2obs) ! eq.(2) optimal weighting of prior and observed estimates 00757 !ap2 s%bwalpha = s%bwalpha*tanh(h1**20) ! tanh(h1**20) included to reduce effect at shallow waters... 00758 00759 00760 00761 do in=1,s%nx+1 00762 do im=1,s%ny+1 00763 alphafac(in,im) = (cosh(s%sdist(in,im)/100.d0-0.65*par%t*s%sdist(s%nx+1,im)/par%tstop/100.d0-2.d0))**(-10.d0) 00764 enddo 00765 enddo 00766 ! max left hand to 1 00767 do im=1,s%ny+1 00768 ind = minval(maxloc(alphafac(:,im))) 00769 alphafac(1:ind,im)=1.d0 00770 enddo 00771 00772 ! s%bwalpha = s%bwalpha*alphafac 00773 00774 bw%dassim = -s%bwalpha*tanh((h1/0.85)**5)*((dDdh-sqrt(alp_as*mxdDdh))/(dDdh**2+alp_as*mxdDdh)*bw%dcmdo )*alphafac !eq.(5) depth change 00775 !ap2 00776 !& +dcdh(nm)/(dcdh(nm)**2+alp_as*dcdhmean**2)*ccmco(nm)) 00777 00778 do in=1,s%nx+1 00779 do im=1,s%ny+1 00780 if (bw%dassim(in,im)>0.) then 00781 bw%dassim(in,im) = min(bw%dassim(in,im),0.1*(h1(in,im))) 00782 elseif (bw%dassim(in,im)<0.) then 00783 bw%dassim(in,im) = max(bw%dassim(in,im),-0.1*(h1(in,im))) 00784 else 00785 bw%dassim(in,im) = 0. 00786 end if 00787 end do 00788 end do 00789 00790 00791 bw%dassim = bw%dassim + s%bwalpha*bw%scmso !!! in xbeach zb=zb-bw%dassim, so zb=zb-s%bwalpha*(zbcomp-zbobs) 00792 ! bw%dassim = bw%dassim - 0.0001*bw%scmso 00793 00794 !expose to output 00795 s%dassim = bw%dassim 00796 00797 do in=1,s%nx+1 00798 do im=1,s%ny+1 !!Roberto changed tanh term below 00799 if (s%dobs(in,im)>-990) s%sig2prior(in,im) = s%bwalpha(in,im)*tanh((h1(in,im)/0.85)**5)*Nassim*sig2obs(in,im) 00800 !! if (s%dobs(in,im)>-990) s%sig2prior(in,im) = s%bwalpha(in,im)*tanh((h2(in,im)/0.85)**5)*Nassim*sig2obs(in,im) 00801 end do 00802 end do 00803 00804 end subroutine comp_depchg 00805 00806 end module beachwizard_module
1.7.6.1
