p_fixhycom_eco.F90 16 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518
  1. ! File: p_fixhycom.F90
  2. !
  3. ! Created: ???
  4. !
  5. ! Last modified: 29/06/2010
  6. !
  7. ! Purpose: Fixes EnKF output.
  8. !
  9. ! Description:
  10. !
  11. ! Modifications:
  12. ! 1/03/2011 Ehouarn:
  13. ! -modification of the fixhycom subroutine: interpolation of
  14. ! biogeochemical tracers on the analysis grid according to
  15. ! hycom remapping in order to insure conservation.
  16. ! 29/06/2010 PS:
  17. ! - set the maximum ICEC to 0.995 to match the model
  18. ! ?/?/? KAL:
  19. ! - Modification of the "fixhycom" subroutine, into separate
  20. ! program, working on a file-by-file basis
  21. ! Prior history:
  22. ! Not documented.
  23. ! KAL -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
  24. ! KAL -- Input arguments:
  25. ! KAL -- template restart file
  26. ! KAL -- ensemble member
  27. ! KAL -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
  28. ! KAL -- TODO: Fix for ice fields
  29. program fixhycom_eco
  30. use mod_raw_io
  31. use m_parse_blkdat
  32. use m_put_mod_fld
  33. use m_get_mod_fld
  34. use m_get_mod_grid
  35. #if defined (ECO)
  36. use m_fixhycom_eco_metno
  37. #endif
  38. implicit none
  39. integer*4, external :: iargc
  40. real, parameter :: onem=9806.
  41. real, parameter :: PSTARB0 = 1000;
  42. integer imem ! ensemble member
  43. character(len=80) :: restart,icerestart ! restart template
  44. character(len=80) :: afile,newfile, char80
  45. integer :: fnd, rstind, tmpindx, iafile
  46. logical :: ex, allok, nomatch
  47. character(len=8) :: cfld, ctmp
  48. character(len=3) :: cproc,cmem
  49. integer :: tlevel, vlevel, nproc
  50. real :: bmin, bmax, rdummy
  51. integer :: idm,jdm,kdm
  52. real, allocatable:: fld(:,:)
  53. real*8, allocatable, dimension(:,:) :: &
  54. ficem,hicem,hsnwm,ticem,tsrfm
  55. real, allocatable, dimension(:,:) :: depths,modlon,modlat,saln
  56. real*4, allocatable:: fldr4(:,:)
  57. real*4 :: spval,amin,amax
  58. real, allocatable :: press(:)
  59. real, allocatable, dimension(:,:) :: dpsum
  60. real, allocatable, dimension(:,:,:) :: dp, dpold
  61. integer,parameter :: numfields=2
  62. integer :: ios,ios2, reclICE,ifld
  63. integer :: i,j,k,ktr
  64. real :: mindx,meandx
  65. #if defined (ECO)
  66. real,dimension(:,:,:), allocatable::dpfor,cfi
  67. character(len=80) :: restfor
  68. real,dimension(:,:,:,:), allocatable::tracerf
  69. real, dimension(:,:), allocatable::trcraij
  70. real, dimension(:), allocatable::prsf,dpf
  71. integer::ntracr,ktrcr
  72. real::dpthin
  73. character(2)::ctrcr
  74. logical, dimension(:), allocatable::lcm
  75. real, dimension(:,:,:), allocatable::temp,sal
  76. integer::kisop
  77. #endif
  78. icerestart=''
  79. if (iargc()==2 .or. iargc()==3) then
  80. call getarg(1,restart)
  81. call getarg(2,ctmp)
  82. read(ctmp,*) imem
  83. write(cmem,'(i3.3)') imem
  84. if (iargc()==3) call getarg(3,icerestart)
  85. else
  86. print *,'"fixhycom" -- A crude routine to correct restart files for obvious errors'
  87. print *
  88. print *,'usage: '
  89. print *,' fixhycom restart_file ensemble_member <ice_file>'
  90. print *,' "restart_file" the restart file you want to fix (.a-file)"'
  91. print *,' "ensemble_member" is the ensemble member - should corr. to that of restart file'
  92. print *,' "ice_file" is optional - it is the restart file for ice fields'
  93. call exit(1)
  94. endif
  95. ! Get dimensions from blkdat
  96. call parse_blkdat('idm ','integer',rdummy,idm)
  97. call parse_blkdat('jdm ','integer',rdummy,jdm)
  98. call parse_blkdat('kdm ','integer',rdummy,kdm)
  99. #if defined (ECO)
  100. call parse_blkdat('ntracr','integer',rdummy,ntracr)
  101. #endif
  102. if (idm>0 .and. idm < 1e4 .and. jdm>0 .and. jdm<1e4) then
  103. allocate(fld (idm,jdm))
  104. allocate(fldr4(idm,jdm))
  105. allocate(saln (idm,jdm))
  106. allocate(ficem(idm,jdm))
  107. allocate(hicem(idm,jdm))
  108. allocate(hsnwm(idm,jdm))
  109. allocate(ticem(idm,jdm))
  110. allocate(tsrfm(idm,jdm))
  111. allocate(dpsum(idm,jdm))
  112. allocate(depths(idm,jdm))
  113. allocate(modlon(idm,jdm))
  114. allocate(modlat(idm,jdm))
  115. allocate(dpold(idm,jdm,kdm))
  116. allocate(dp (idm,jdm,kdm))
  117. allocate(press(kdm+1))
  118. #if defined (ECO)
  119. allocate(dpfor(idm,jdm,kdm))
  120. allocate(cfi(kdm,ntracr,2))
  121. allocate(prsf(kdm+1))
  122. allocate(tracerf(idm,jdm,kdm,ntracr))
  123. allocate(trcraij(kdm,ntracr))
  124. allocate(lcm(kdm))
  125. allocate(dpf(kdm))
  126. allocate(sal(idm,jdm,kdm))
  127. allocate(temp(idm,jdm,kdm))
  128. #endif
  129. else
  130. print *,'fld allocate error'
  131. stop '(EnKF_postprocess)'
  132. end if
  133. ! Remove postfix of restart file
  134. fnd=max(index(restart,'.a'),index(restart,'.b'))
  135. ! Inquire for existence
  136. inquire(exist=ex,file=restart(1:fnd-1)//'.b')
  137. if (.not.ex) then
  138. write(*,*) 'Can not find '//restart(1:fnd-1)//'.b'
  139. stop '(EnKF_postprocess)'
  140. end if
  141. print *,restart(1:fnd-1)//'.b'
  142. newfile='fix'//restart(1:fnd-1)
  143. ! Get model grid
  144. call get_mod_grid(modlon,modlat,depths,mindx,meandx,idm,jdm)
  145. #if defined (ECO)
  146. !files where are stored the forecast fields!
  147. restfor='forecast'//cmem
  148. dpthin = onem*0.001
  149. #endif
  150. ! Get layer thickness
  151. dpsum=0.
  152. do k=1,kdm
  153. !call get_mod_fld(dp(:,:,k),1,'dp ',k,1)
  154. call get_mod_fld_new(restart(1:fnd-1),dp(:,:,k),imem,'dp ',k,1,idm,jdm)
  155. dpsum=dpsum+dp(:,:,k)
  156. #if defined (ECO)
  157. !reading of forecast fields: tracers, T and S
  158. call get_mod_fld_new(trim(restfor),dpfor(:,:,k),imem,'dp ',k,1,idm,jdm)
  159. do ktrcr=1,ntracr
  160. write(ctrcr,'(i2.2)') ktrcr
  161. cfld='tracer'//ctrcr
  162. call get_mod_fld_new(trim(restfor),tracerf(:,:,k,ktrcr),imem,cfld,k,1,idm,jdm)
  163. enddo
  164. call get_mod_fld_new(trim(restfor),temp(:,:,k),imem,'temp ',k,1,idm,jdm)
  165. call get_mod_fld_new(trim(restfor),sal(:,:,k),imem,'saln ',k,1,idm,jdm)
  166. #endif
  167. end do
  168. print *,maxval(dpsum-depths*onem)
  169. dpold=dp
  170. ! DP correction
  171. do j=1,jdm
  172. do i=1,idm
  173. !!! Move negative layers to neighbouring layers.
  174. do k = 1, kdm-1
  175. dp(i,j,k+1) = dp(i,j,k+1) + min(0.0,dp(i,j,k))
  176. dp(i,j,k ) = max(dp(i,j,k),0.0)
  177. !dp(i,j,k ) = max(dp(i,j,k),1.e-3*onem)
  178. end do
  179. !!! Go backwards to fix lowermost layer.
  180. do k = kdm, 3, -1
  181. dp(i,j,k-1) = dp(i,j,k-1) + min(0.0,dp(i,j,k))
  182. dp(i,j,k) = max(dp(i,j,k),0.0)
  183. !dp(i,j,k) = max(dp(i,j,k),1.e-3*onem)
  184. end do
  185. !!! No layers below the sea bed.
  186. press( 1) = 0.0
  187. #if defined (ECO)
  188. !computation of the forecast layer interfaces (prsf)!
  189. prsf(1)=0.
  190. #endif
  191. do k = 1, kdm-1
  192. press(k+1) = press(k) + dp(i,j,k)
  193. press(k+1) = min(depths(i,j)*onem,press(k+1))
  194. #if defined (ECO)
  195. prsf(k+1) = prsf(k) + dpfor(i,j,k)
  196. #endif
  197. end do
  198. press(kdm+1) = depths(i,j)*onem
  199. #if defined (ECO)
  200. prsf(kdm+1)=depths(i,j)*onem
  201. #endif
  202. do k = 1, kdm
  203. dp(i,j,k) = press(k+1) - press(k)
  204. #if defined (ECO)
  205. !definition of the isopycnal layers!
  206. dpf(k)=max(dpfor(i,j,k),dpthin)
  207. kisop=compute_kisop(temp(i,j,:),sal(i,j,:),kdm)
  208. if (k.le.max(2,kisop)) then
  209. lcm(k) = .false. !fixed layers are never PCM
  210. else
  211. ! --- thin and isopycnal layers remapped with PCM.
  212. lcm(k) = dpfor(i,j,k).le.dpthin
  213. endif
  214. #endif
  215. end do
  216. !eho 2/11/11
  217. if(depths(i,j)>10000. .or. depths(i,j)<1.)then
  218. dp(i,j,:)=dpold(i,j,:)
  219. endif
  220. #if defined (ECO)
  221. ! if(depths(i,j)==0.)then
  222. if(depths(i,j)>10000. .or. depths(i,j)<1.)then
  223. cycle
  224. else
  225. call hybgen_weno_coefs(tracerf(i,j,:,:),dpf,lcm,cfi,kdm,ntracr,dpthin)
  226. call hybgen_weno_remap(tracerf(i,j,:,:),prsf,dpfor(i,j,:),cfi,trcraij,&
  227. press,kdm,kdm,ntracr,dpthin)
  228. tracerf(i,j,1:kdm,1:ntracr)=trcraij(1:kdm,1:ntracr)
  229. endif
  230. #endif
  231. end do
  232. end do
  233. do k = 1, kdm
  234. print *,maxval(dpold(:,:,k)-dp(:,:,k))/onem
  235. end do
  236. #if defined (ECO)
  237. deallocate(trcraij,dpf,lcm,cfi,prsf,dpfor)
  238. #endif
  239. ! Loop over restart file
  240. rstind=1 ! Restart index
  241. allok=.true.
  242. do while ( allok)
  243. ! Get header info from restart
  244. call rst_header_from_index(restart(1:fnd-1)//'.b', &
  245. cfld,vlevel,tlevel,rstind,bmin,bmax,.true.)
  246. allok=tlevel/=-1 ! test to see if read was ok
  247. print *,cfld
  248. if (allok) then
  249. ! Get actual field - for now we use the READRAW routine (later on we
  250. ! should switch to m_get_mod_fld
  251. ! call READRAW(fldr4,amin,amax,idm,jdm,.false.,spval,restart(1:fnd-1)//'.a',rstind)
  252. ! fld=fldr4
  253. ! Sjekk p at vi har lest rett - samanlign max/min fr filene
  254. ! if (abs(amin-bmin).gt.abs(bmin)*1.e-4 .or. &
  255. ! abs(bmax-amax).gt.abs(bmax)*1.e-4 ) then
  256. ! print *,'Inconsistency between .a and .b files'
  257. ! print *,'.a : ',amin,amax
  258. ! print *,'.b : ',bmin,bmax
  259. ! print *,cfld,vlevel,tlevel
  260. ! call exit(1)
  261. ! end if
  262. call get_mod_fld_new(restart(1:fnd-1),fld(:,:),imem,cfld,vlevel,1,idm,jdm)
  263. if (trim(cfld)=='temp') then
  264. ! need salinity as well
  265. call get_mod_fld_new(restart(1:fnd-1),saln(:,:),imem,'saln ',vlevel,1,idm,jdm)
  266. !if (tlevel==-1) then
  267. ! print *,'Could not get salinity field'
  268. ! call exit(1)
  269. !end if
  270. ! keep water warmer than freezing point
  271. do j=1,jdm
  272. do i=1,idm
  273. fld(i,j)=max(-.057*saln(i,j),fld(i,j))
  274. fld(i,j)=min(35.,fld(i,j)) !FC: cut off values that are too warm
  275. end do
  276. end do
  277. else if (trim(cfld)=='saln') then
  278. do j=1,jdm
  279. do i=1,idm
  280. fld(i,j)=max(5.,fld(i,j)) ! LB :no water fresher than 5 psu (Baltic)
  281. fld(i,j)=min(41.,fld(i,j))! FC: no water saltier than 40 psu
  282. end do
  283. end do
  284. else if (trim(cfld)=='pstarb') then
  285. fld(:,:) = (sqrt(PSTARB0 ** 2 + fld(:,:) ** 2) + fld(:,:)) / 2.0d0;
  286. else if (trim(cfld)=='dp') then
  287. fld = dp(:,:,vlevel) ! NB, one time level
  288. #if defined (ECO)
  289. else if (cfld(1:6)=='tracer') then
  290. !updating the file!
  291. ktrcr=tracr_get_incr(cfld(7:8))
  292. if (ktrcr==-1)then
  293. print*,'alert tracer unknow'
  294. exit
  295. endif
  296. fld(:,:)= tracerf(:,:,vlevel,ktrcr)
  297. #endif
  298. end if ! No correction for other fields in the hycom restart file
  299. call put_mod_fld(trim(newfile),fld,imem,cfld,vlevel,tlevel,rstind,idm,jdm)
  300. rstind=rstind+1
  301. end if ! read of template was ok
  302. end do
  303. #if defined (ECO)
  304. deallocate(tracerf)
  305. #endif
  306. ! put_mod_fld does not include the header from the original file
  307. open(10,file=trim(newfile)//'.b',status='old')
  308. open(20,file=restart(1:fnd-1)//'.b',status='old')
  309. ! Supports parallel execution with different members
  310. open(11,file='tmp'//cmem//'.b',status='replace')
  311. ! Header from original file
  312. read(20,'(a80)') char80 ; write(11,'(a80)') char80
  313. read(20,'(a80)') char80 ; write(11,'(a80)') char80
  314. close(20)
  315. ! The rest from the newly created file
  316. ios=0
  317. do while(ios==0)
  318. read(10,'(a80)',iostat=ios) char80
  319. if (ios==0) write(11,'(a80)') char80
  320. end do
  321. close(10)
  322. close(11)
  323. close(20)
  324. ! Move the new .b file to "newfile"
  325. !SERIAL call system('mv tmp'//cmem//'.b '//trim(newfile)//'.b')
  326. !###################################################################
  327. !####################### FIX ICE MODEL #########################
  328. !###################################################################
  329. if (iargc()==3) then
  330. inquire(exist=ex,file=trim(icerestart))
  331. if (.not.ex) then
  332. print *,icerestart//' does not exist!'
  333. print *,'(fixhycom)'
  334. stop
  335. end if
  336. inquire(iolength=reclICE)ficem,hicem,hsnwm,ticem,tsrfm !,iceU,iceV
  337. open(10,file=icerestart,form='unformatted',access='direct',recl=reclICE)
  338. read(10,rec=imem)ficem,hicem,hsnwm,ticem,tsrfm !,iceU,iceV
  339. close(10)
  340. ! PS 25/06/2010 max(ficem) = 0.995 - max the model allows
  341. ficem=min(max(0.,ficem), 0.995)
  342. hicem=min(max(0.,hicem),15.)
  343. hsnwm=min(max(0.,hsnwm), 4.)
  344. open (10,file='fix'//trim(icerestart),form='unformatted',access='direct',recl=reclICE)
  345. write(10,rec=imem)ficem,hicem,hsnwm,ticem,tsrfm !,iceU,iceV
  346. close(10)
  347. end if
  348. !KAL
  349. !KAL ! This is to make it easier for supporting programs
  350. !KAL open(10,file='file.EnKF_postprocess',status='replace')
  351. !KAL write(10,'(a)') 'analysis'//cmem
  352. !KAL close(10)
  353. !KAL
  354. !KAL
  355. !KAL
  356. !KAL ! ice processing Loop over restart file
  357. !KAL print *
  358. !KAL print *
  359. !KAL print *,'processing ice restart file'
  360. !KAL rstind=1 ! Restart index
  361. !KAL allok=.true.
  362. !KAL call system("cp ensemble_TMP_ICE.uf ensemble_TMP_ICE_final.uf")
  363. !KAL do ifld=1,numfields
  364. !KAL cfld= fieldnames (ifld)
  365. !KAL vlevel=fieldlevels(ifld)
  366. !KAL if (trim(cfld)=='icec' .or. trim(cfld)=='hice') then
  367. !KAL nomatch=.true.
  368. !KAL do iafile=1,nproc ! List of procs used in analysis
  369. !KAL write(cproc,'(i3.3)') iafile-1
  370. !KAL
  371. !KAL ! NB - time level=1
  372. !KAL ! NB2 - the files dumped in the analysis lack a header (last argument
  373. !KAL ! is false)
  374. !KAL call rst_index_from_header(trim(afile)//'.b',cfld,vlevel,1, &
  375. !KAL tmpindx,bmin,bmax,.false.)
  376. !KAL
  377. !KAL
  378. !KAL if (tmpindx/=-1) then
  379. !KAL print '(a8," -- layer:",i4," match : record, file",i4," ",a)', cfld, vlevel,tmpindx,trim(afile)
  380. !KAL nomatch=.false.
  381. !KAL exit
  382. !KAL end if
  383. !KAL
  384. !KAL ! Read field from analysed file
  385. !KAL call READRAW(fldr4,amin,amax,idm,jdm,.false.,spval,trim(afile)//'.a',tmpindx)
  386. !KAL fld=fldr4
  387. !KAL
  388. !KAL ! Sjekk p at vi har lest rett - samanlign max/min fr filene
  389. !KAL if (abs(amin-bmin).gt.abs(bmin)*1.e-4 .or. &
  390. !KAL abs(bmax-amax).gt.abs(bmax)*1.e-4 ) then
  391. !KAL print *,'Inconsistency between .a and .b files'
  392. !KAL print *,'.a : ',amin,amax
  393. !KAL print *,'.b : ',bmin,bmax
  394. !KAL print *,cfld,vlevel,tlevel
  395. !KAL call exit(1)
  396. !KAL end if
  397. !KAL
  398. !KAL end do
  399. !KAL
  400. !KAL
  401. !KAL ! Check if we got ice concentration or ice thickness
  402. !KAL if (.not.nomatch) then
  403. !KAL
  404. !KAL
  405. !KAL inquire(iolength=reclICE)ficem,hicem,hsnwm,ticem,tsrfm !,iceU,iceV
  406. !KAL open(10,file=trim(icetemplate),form='unformatted',access='direct',recl=reclICE,status='old')
  407. !KAL open(11,file='ensemble_TMP_ICE_final.uf',form='unformatted',access='direct',recl=reclICE,status='unknown')
  408. !KAL read(10,rec=imem,iostat=ios)ficem,hicem,hsnwm,ticem,tsrfm !,iceU,iceV
  409. !KAL if (trim(cfld)=='icec') ficem=fld
  410. !KAL if (trim(cfld)=='hice') hicem=fld
  411. !KAL write(11,rec=imem,iostat=ios2)ficem,hicem,hsnwm,ticem,tsrfm !,iceU,iceV
  412. !KAL close(10)
  413. !KAL close(11)
  414. !KAL
  415. !KAL if (ios/=0 .or.ios2/=0) then
  416. !KAL print *,ios
  417. !KAL print *,ios2
  418. !KAL print *,'Error when writing to ice ens file'
  419. !KAL call exit(1)
  420. !KAL end if
  421. !KAL end if
  422. !KAL end if
  423. !KAL end do
  424. !KAL
  425. !KAL
  426. !KAL print *,'Normal exit of EnKF_postprocess'
  427. !KAL print *,'TODO: Process dp inconsistencies'
  428. !KAL call exit(0)
  429. !KAL
  430. !KAL
  431. !KAL
  432. end program