module mod_gridp c --- TODO: Complete module info C --- Module contains tools for generating gridpoint data C --- C --- subroutine gridpinp - initialize gridpoint information arrays C --- subroutine gridpdef - Allocate arrays used use mod_xc implicit none c --- Shield everything by default private type grid_points integer*4 i integer*4 j real*4 z(kdm) real*4 u(kdm) real*4 v(kdm) real*4 t(kdm) real*4 s(kdm) real*4 sealevel real*4 uwind real*4 vwind #if defined (ICE) real*4 fice real*4 hice real*4 uice real*4 vice #endif end type grid_points type(grid_points), allocatable :: gp(:) ! Array of grid point information type(grid_points), allocatable :: gpstore(:,:) ! Keeps info for diff times and points integer, allocatable :: itime(:) ! holds time information for each record integer, allocatable :: irecnr(:) ! record number logical, allocatable :: lrecnr(:) ! Flag for record number which has been dumped integer gpdim ! Total nuber of grid points to be saved integer gptdim ! Total times in temporary array integer gpit c --- Public flag for activating gridp - set in m_limits logical, save, public :: lgridp c c --- These are made available for hycom public :: gridpinit, gridp_process contains c --- ----------------------------------------------------------------- c --- Initialization of gridp routines c --- ----------------------------------------------------------------- subroutine gridpinit() use mod_diagnostics ! contains iday1, dday implicit none if (lgridp) then call gridpinp() call gridpdef(iday1,dday) else if (mnproc==1) & write(lp,*) 'NO GRIDP information stored: lgridp=',lgridp endif end subroutine c --- ----------------------------------------------------------------- c --- gridp processing logic, moved from hycom in here c --- ----------------------------------------------------------------- subroutine gridp_process(rt,dtime,m,n) use mod_diagnostics use mod_year_info implicit none type(year_info) , intent(in) :: rt integer , intent(in) :: m,n real*8 , intent(in) :: dtime type(year_info) :: rttmp include 'common_blocks.h' c --- Grid points if (lgridp .and. rt%iss < nint(baclin)) then call xcstop('(NERSC gridp - fix times )') call gridpstore(n,m,rt) endif c --- Gridp dump -- not for initial dump if (nstep>nstep1+1) then call year_day(dtime+baclin/86400.d0,refyear,rttmp,yrflag) if (lgridp.and.loutput) then if (mnproc==1) write(lp,*)'Gridp output on diagnostic day' call gridpdump(n,rt) if (outday < iday2) then call gridpdef(outday,dday) end if c --- Always dump at end of month else if (lgridp .and.rttmp%imm /= rt%imm) then if (mnproc==1) write(lp,*)'Gridp output on month changeover' call gridpdump(n,rt) if (outday < iday2) then call gridpdef(outday,dday) end if endif endif end subroutine C --- This routine reads the input files, and sets a mask to determine C --- which points are dumped later on. The array "gp" (see above) gets C --- initialized with grid indices corresponding to lon/lat positions subroutine gridpinp() implicit none integer i,j,m,n,nrgps real mind,maxd,minlat,maxlat,minlon,maxlon integer igpskip, jgpskip, ibnd logical ex logical :: gpmask(itdm,jtdm) real, dimension(itdm,jtdm) :: gdepths,modlon,modlat character*3 cthread include 'common_blocks.h' if (mnproc==1) write(lp,*) 'GRIDP activated' call xcaget(gdepths,depths,0) call xcaget(modlon,plon,0) call xcaget(modlat,plat,0) C$OMP PARALLEL DO PRIVATE(i,j) do j=1,jtdm do i=1,itdm gpmask(i,j)=.false. end do end do C$OMP END PARALLEL DO c --- Read infile_gp.in inquire(file='infile_gp.in',exist=ex) if (.not.ex) then if (mnproc==1) & write(lp,*) 'infile_gp.in file does not exist' call xcstop('(gridpinp)') stop '(gridpinp)' end if open(11,file='infile_gp.in') read(11,*)mind ! Minimum depth read(11,*)maxd ! Max depth read(11,*)minlat ! Minimum lat (lonlat box) read(11,*)maxlat ! Maximum Lat (lonlat box) read(11,*)minlon ! Minimum lon (lonlat box) read(11,*)maxlon ! inximum lon (lonlat box) read(11,*)igpskip ! Index skip read(11,*)jgpskip ! Index skip c --- counting number of grid points for grid C$OMP PARALLEL DO PRIVATE(i,j) do j=2,jtdm-1,jgpskip do i=2,itdm-1,igpskip if (mind <= gdepths(i,j).and.gdepths(i,j)<=maxd & .and.minlon <= modlon(i,j) .and.modlon(i,j) <=maxlon & .and.minlat <= modlat(i,j) .and.modlat(i,j) <=maxlat ) then gpmask(i,j)=.true. end if enddo enddo C$OMP END PARALLEL DO if (mnproc==1)print *,'Number of regular points : ',count(gpmask) c --- From now on we read specified stations. Input is now i,j points c --- (NB: used to be lon lat values) do m=1,1000 read(11,*,end=100,err=100) i,j if (mnproc==1) print '(a,2i5)','gridpinp: ',i,j if (1 < i .and. i < itdm .and. 1 < j .and. j < jtdm) then gpmask(i ,j )=gdepths(i ,j ) > 0.0 gpmask(i+1,j )=gdepths(i+1,j ) > 0.0 gpmask(i ,j+1)=gdepths(i ,j+1) > 0.0 gpmask(i+1,j+1)=gdepths(i+1,j+1) > 0.0 else if (mnproc==1) write(lp,*) 'Station NOT on grid !!' call xcstop ('(gridpinp)') stop '(gridpinp)' endif enddo 100 close(11) c --- Stop if there are no points: if (all(.not.gpmask)) then if (mnproc==1) print *,'Error: There are no gp points!' call xcstop('(gridpinp)') stop '(gridpinp)' end if gpdim=0 C$OMP PARALLEL DO PRIVATE (i,j) REDUCTION(+:gpdim) do j=1,jtdm do i=1,itdm if (gpmask(i,j)) gpdim=gpdim+1 enddo enddo C$OMP END PARALLEL DO if (mnproc==1) then write(lp,*) 'gridpinp: Number of grid points is : ', & gpdim write(lp,*) 'gridpinp: Size of each file per year is (Mb): ', & (kdm*5+3)*24*4*365*1.0E-06 write(lp,*) 'gridpinp: Total archive size per year (Gb)) : ', & (kdm*5+3)*24*4*365*1.0E-09*gpdim end if c --- Now calculate for tiles: gpdim=0 C$OMP PARALLEL DO PRIVATE (i,j) REDUCTION(+:gpdim) do j=1,jj do i=1,ii if (gpmask(i0+i,j0+j)) then gpdim=gpdim+1 end if end do end do C$OMP END PARALLEL DO if (gpdim>0) allocate(gp(gpdim)) call xcsync(flush_lp) write(lp,'(a,2i5,a,i5)') &' gridpinp: Number of grid points for tile ', mproc,nproc, &' is ', gpdim call xcsync(flush_lp) c --- define local indexes for gpdata if (gpdim>0) then m=0 do j=1,jj do i=1,ii if (gpmask(i0+i,j0+j)) then m=m+1 gp(m)%i=i gp(m)%j=j endif enddo enddo c --- Save to tecplot file - one for each tile write(cthread,'(i3.3)') mnproc open(10,file='gridploc'//cthread//'.tec') do m=1,gpdim i=gp(m)%i j=gp(m)%j write(10,'(2f10.5,2i4)')plon(i,j),plat(i,j),i,j enddo close(10) end if c -- lgridp not switched off for now, but we might deactivate it later for "empty" tiles if (.not.lgridp) then write(lp,'(a,i5)') 'GRIDP DE-activated for tile ',ijqr end if end subroutine gridpinp C --- Routine Allocates global (in module) vars C --- Allocates: C --- gpstore C --- itime C --- irecnr C --- lrecnr subroutine gridpdef(iday,dday) use mod_diagnostics use mod_year_info implicit none type (diass), intent(in) :: dday(0:nrdday) integer, intent(in) :: iday ! Number of times in dump gptdim=nint((dday(iday+1)%day-dday(iday )%day)*24.0) if (mnproc==1) print *,'gridpdef: gptdim=',gptdim if (allocated(gpstore)) deallocate(gpstore) allocate(gpstore(gptdim,gpdim)) if (allocated(itime)) deallocate(itime) allocate(itime(gptdim)) itime=0 if (allocated(irecnr)) deallocate(irecnr) allocate(irecnr(gptdim)) irecnr=0 if (allocated(lrecnr)) deallocate(lrecnr) allocate(lrecnr(gptdim)) lrecnr(:)=.false. end subroutine gridpdef C --- Get gp data from hycom fields, once an hour C --- Put them into the gp variables subroutine gridpstore(n,m,rt) use mod_year_info use mod_forcing_nersc #if defined (ICE) use mod_common_ice #endif implicit none integer , intent(in) :: n ! New time indexes integer , intent(in) :: m ! old time indexes type(year_info) , intent(in) :: rt ! Time variable real, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,kdm+1) :: & pn,pm,tprs ! pressure at interfaces real, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: & tuwind,tvwind, tslp real, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,kdm) :: & tsaln,ttemp,tu,tv type(year_info) rtt character(len=17) tag17 real fli,flj real up,vp,ud,vd,theta_up,theta_vp real dlon,dlat integer ilon,ilat character(len=10)ctmp integer i,j,k,kn,igp,km,l real :: wn,wm include 'common_blocks.h' ! Put integer description of time into itime ctmp(1:4)=rt%cyy ctmp(5:6)=rt%cmm ctmp(7:8)=rt%cdm ctmp(9:10)=rt%chh read(ctmp,'(i10.10)')itime(gpit) !print *,ctmp,' ',rt%cyy,' ',rt%cmm,' ',rt%cdm,itime(gpit) if (gpit == gptdim .and. gptdim > 1) then irecnr(gpit)=irecnr(gpit-1)+1 else irecnr(gpit)=rt%idm*24+rt%ihh+1 endif if (mod(gpit,12)==0) print '(a,i5,a,i10,a,i5,a,4i4,a,i5)', & 'gridpstore: gpit=',gpit, & ' itime=',itime(gpit), & ' recnr=',irecnr(gpit), & ' rt= ',rt%iyy, rt%imm, rt%idm, rt%ihh, & ' tile= ',ijqr if (lrecnr(gpit)) then if (mnproc==1) print *, & 'gpstore: trying to rewrite a written record' call xcstop('gpstore') stop 'gpstore' endif lrecnr(gpit)=.true. ! This routine can now be called when rt%iss/=0 (see hycom.F90). We must do ! some sort of interpolation in order to get hourly values. Enter linear ! interpolation, which is good enough if (rt%iss>baclin) then if (mnproc==1) & print *,'(gridpstore should only be called on the hour or'// & ' immediately after the hour has passed)' call xcstop('gripstore') stop '(gripstore)' end if wm=rt%iss/baclin wn=1.-wm !print *,'gpstore wn wm iss baclin',wn,wm,rt%iss,baclin ! Get slp and winds !$OMP PARALLEL DO PRIVATE(j,i) !$OMP& SCHEDULE(STATIC,jblk) do j=1,jj do i=1,ii tuwind(i,j)=uwind(i,j,l0)*w0+uwind(i,j,l1)*w1 & +uwind(i,j,l2)*w2+uwind(i,j,l3)*w3 tvwind(i,j)=vwind(i,j,l0)*w0+vwind(i,j,l1)*w1 & +vwind(i,j,l2)*w2+vwind(i,j,l3)*w3 tslp (i,j)=slp (i,j,l0)*w0+slp (i,j,l1)*w1 & +slp (i,j,l2)*w2+slp (i,j,l3)*w3 end do end do c do k=1,kk+1 !$OMP PARALLEL DO PRIVATE(j,l,i) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do l=1,isp(j) do i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l)) tprs(i,j,k) = wn*pn(i,j,k) + wm*pm(i,j,k) end do end do end do C$OMP END PARALLEL DO end do ! Get U,V,saln and temp do k=1,kk !$OMP PARALLEL DO PRIVATE(j,l,i) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do l=1,isp(j) do i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l)) tsaln(i,j,k) = wn*saln (i,j,k,n) + wm*saln (i,j,k,m) ttemp(i,j,k) = wn*temp (i,j,k,n) + wm*temp (i,j,k,m) tu (i,j,k) = wn*u (i,j,k,n) + wm*u (i,j,k,m) & + wn*ubavg(i,j,n) + wm*ubavg(i,j,m) tv (i,j,k) = wn*v (i,j,k,n) + wm*v (i,j,k,m) & + wn*vbavg(i,j,n) + wm*vbavg(i,j,m) end do end do end do C$OMP END PARALLEL DO end do ! Interpolation time !$OMP PARALLEL DO PRIVATE(igp,i,j,dlon,dlat,theta_up,theta_vp, !$OMP& k,kn,ud,vd,km) SHARED(gpit,m,n) !$OMP& SCHEDULE(STATIC,1) do igp=1,gpdim !if (mod(igp,10) == 0) write(*,'(a1)',advance='no')'.' i=gp(igp)%i j=gp(igp)%j dlon=plon(i+1,j)-plon(i-1,j) dlat=plat(i+1,j)-plat(i-1,j) ! Grid directions -- "i" rel lon/lat if(dlon.LT.180.)dlon=360.0+dlon if(dlon.GT.180.)dlon=dlon-360.0 CKAL -- From HYCOM 2.1.34 - radian is pi/180 theta_up=atan2(dlat,dlon*cos(radian*.5* & (plat(i-1,j)+plat(i+1,j)))) ! Grid directions -- "j" rel lon/lat dlon=plon(i,j+1)-plon(i,j-1) dlat=plat(i,j+1)-plat(i,j-1) if(dlon.LT.180.)dlon=360.0+dlon if(dlon.GT.180.)dlon=dlon-360.0 CKAL -- From HYCOM 2.1.34 - radian is pi/180 theta_vp = atan2(dlat,dlon*cos(radian*.5* & (plat(i,j-1)+plat(i,j+1)))) do k=1,kdm gp(igp)%z(k)=tprs(i,j,k+1)/onem gp(igp)%u(k)=tu(i,j,k)*cos(theta_up)+tv(i,j,k)*cos(theta_vp) gp(igp)%v(k)=tu(i,j,k)*sin(theta_up)+tv(i,j,k)*sin(theta_vp) gp(igp)%t(k)=ttemp(i,j,k) gp(igp)%s(k)=tsaln(i,j,k) enddo ! Wind rotation introduced AFTER NWAG and WANE projects by GE gp(igp)%uwind =tuwind(i,j)*cos(theta_up) & +tvwind(i,j)*cos(theta_vp) gp(igp)%vwind =tuwind(i,j)*sin(theta_up) & +tvwind(i,j)*sin(theta_vp) gp(igp)%sealevel=srfhgt(i,j) #if defined (ICE) gp(igp)%hice = hicem(i,j) gp(igp)%fice = ficem(i,j) gp(igp)%uice = iceu(i,j)*cos(theta_up)+ & icev(i,j)*cos(theta_vp) gp(igp)%vice = iceu(i,j)*sin(theta_up)+ & icev(i,j)*sin(theta_vp) #endif gpstore(gpit,igp)=gp(igp) enddo C$OMP END PARALLEL DO ! open(10,file='gptest.uf') ! write(*,*)'gridpstore: gptest.dat' ! write(10,*)gp ! close(10) ! write(*,*)'done' end subroutine gridpstore C --- Put grid point data on file subroutine gridpdump(n,rt) use mod_xc use mod_year_info use mod_hycom_nersc implicit none integer, intent(in) :: n type(year_info), intent(in) :: rt type(year_info) rtt character(len=17) tag17 real fli,flj real up,vp,ud,vd,theta_up,theta_vp real dlon,dlat integer ilon,ilat,it integer ilen integer m,i,j character(len=10)ctmp character(len=4) yyyy character(len=12) gpdir character(len=31) gpfile character(len=10) tag10a,tag10b logical ex integer :: firstrec, lastrec real*4 :: r4lon, r4lat, r4deep type(grid_points) :: tgp ! temporary gp logical,save :: lfirst=.true. include 'common_blocks.h' firstrec=gptdim+1 lastrec=1 do i=1,gptdim if(irecnr(i)/=0) then firstrec=min(i,firstrec) lastrec =max(i,lastrec) end if end do c --- Put itime into tag10 write(tag10a,'(i10.10)')itime(lastrec) write(tag10b,'(i10.10)')itime(lastrec) ! This should work now if (mnproc==1) then print *,'itime(lastrec) = ',itime(lastrec) end if c --- Name of directory to store GP data gpdir(1:3)=rungen gpdir(4:5)='GP' gpdir(6:9)=tag10b(1:4) gpdir(10:10)='_' gpdir(11:12)=tag10b(5:6) inquire(file=gpdir,exist=ex) if (.not.ex.and.mnproc==1) then write(*,'(a,a)')'GP: creating the directory ./', gpdir call mkdir_wrap( trim(gpdir)//char(0) ) endif call xcsync(flush_lp) c --- dump legend in text file in GP dir if (mnproc==1) then open(10,file=gpdir//'/'//rungen//'gprecord.asc', & status='unknown') write(10,'(a12,a3)') 'var name ','dim' write(10,'(a12,i3)') 'itime ',1 write(10,'(a12,i3)') 'i ',1 write(10,'(a12,i3)') 'j ',1 write(10,'(a12,i3)') 'lon ',1 write(10,'(a12,i3)') 'lat ',1 write(10,'(a12,i3)') 'depth ',1 write(10,'(a12,i3)') 'intf ',kdm write(10,'(a12,i3)') 'ut ',kdm write(10,'(a12,i3)') 'vt ',kdm write(10,'(a12,i3)') 'temp ',kdm write(10,'(a12,i3)') 'saln ',kdm write(10,'(a12,i3)') 'ssh ',1 write(10,'(a12,i3)') 'uwind ',1 write(10,'(a12,i3)') 'vwind ',1 #if defined (ICE) write(10,'(a12,i3)') 'fice ',1 write(10,'(a12,i3)') 'hice ',1 write(10,'(a12,i3)') 'uice ',1 write(10,'(a12,i3)') 'vice ',1 #endif close (10) end if c --- Get IO-length for one record inquire(iolength=ilen)itime(1),tgp%i,tgp%j, & r4lon, r4lat, r4deep, ! lon/lat/depth & tgp%z,tgp%u,tgp%v,tgp%t,tgp%s, & tgp%sealevel, tgp%uwind, tgp%vwind #if defined (ICE) & ,tgp%fice,tgp%hice,tgp%uice,tgp%vice #endif c --- Go through the points and dump data do m=1,gpdim i=gp(m)%i j=gp(m)%j if (plon(i,j) < 0.0) then tag17(1:1)='-' else tag17(1:1)='+' endif c --- Write lon/lat into file name ilon=int(abs(plon(i,j))) ilat=int(abs(plat(i,j))) write(tag17(2:4),'(I3.3)')ilon tag17(5:5)='.' write(tag17(6:8),'(I3.3)') & nint((abs(plon(i,j))-float(ilon))*1000.0) tag17(9:9)='x' if (plat(i,j) < 0.0) then tag17(10:10)='-' else tag17(10:10)='+' endif write(tag17(11:13),'(I3.3)')ilat tag17(14:14)='.' write(tag17(15:17),'(I3.3)')nint((abs(plat(i,j))-float(ilat)) & *1000.0) c --- Create file name gpfile( 1: 3)='gp_' gpfile( 4: 7)=tag10b(1:4) gpfile( 8: 8)='_' gpfile( 9:10)=tag10b(5:6) gpfile(11:11)='_' gpfile(12:28)=tag17 gpfile(29:31)='.uf' c --- Diagnostics for each tile if (m==1) then print '(a,2i4,a,i6,a,i10)', & 'gpdump: tile ',mproc,nproc, & ' First record and time : ', & irecnr(firstrec),' ',itime(firstrec) print '(a,2i4,a,i6,a,i10)', & 'gpdump: tile ',mproc,nproc, & ' Last record and time : ', & irecnr(lastrec),' ', itime(lastrec) end if if (lfirst) then print '(a,2i3,a)', & 'gpdump: tile ',mproc,nproc, & ' Dumping gridp data to '//gpdir//'/'//gpfile end if c --- No files should be the same for different tiles c --- so this should be ok. open(10,file=gpdir//'/'//gpfile,form='unformatted', & access='direct',recl=ilen,status='unknown') do it=firstrec,lastrec if (irecnr(it) /= 0) then tgp=gpstore(it,m) r4lon =plon(tgp%i,tgp%j) r4lat =plat(tgp%i,tgp%j) r4deep=depths(tgp%i,tgp%j) write(10,rec=irecnr(it)) itime(it),tgp%i,tgp%j, & r4lon, r4lat, r4deep, & tgp%z,tgp%u,tgp%v,tgp%t,tgp%s, & tgp%sealevel, tgp%uwind, tgp%vwind #if defined (ICE) & ,tgp%fice,tgp%hice,tgp%uice,tgp%vice #endif else print *,'gridpdump: record for it= ',it, & ' not written, value= ',irecnr(it), & ' for tile ',ijqr endif enddo close(10) enddo if (mnproc == 1) then print *,'gpdump: done...' print * end if lfirst=.false. end subroutine gridpdump end module mod_gridp