femwind/fortran/module_wrfout.f90

   1 module module_wrfout
   2 ! *** wrf specific netcdf i/o
   3
   4 use module_netcdf
   5 use module_msleep
   6
   7 contains
   8
   9
  10 subroutine netcdf_read_array_wrf(ncid,name,frame,sr,a1d,a2d,a3d)
  11     implicit none
  12
  13 !*** arguments
  14     integer, intent(in)::                ncid ! open netcdf file
  15     real, pointer, intent(out),optional:: a1d(:)  ! the array pointer; remember to deallocate when done with it
  16     real, pointer, intent(out),optional:: a2d(:,:)  ! the array pointer; remember to deallocate when done with it
  17     real, pointer, intent(out),optional:: a3d(:,:,:)  ! the array pointer; remember to deallocate when done with it
  18     character(LEN=*),intent(in):: name
  19     integer, intent(in), optional::                frame ! time step in the file
  20     integer, intent(in), dimension(2), optional::  sr ! strip to remove in i and j dimension
  21
  22 !*** local
  23     integer d1,d2,d3
  24     integer,dimension(:),allocatable::star,cnts,dims,ends
  25     integer::i,j,k,varid,ndims
  26     real,dimension(:,:),allocatable::a1
  27     real,dimension(:,:,:),allocatable::a2
  28     real,dimension(:,:,:,:),allocatable::a3
  29     integer:: istep=1, srf(2)
  30     character(len=256) msg
  31
  32 !*** executable
  33     if(present(frame))istep=frame
  34
  35     d1 = l2i(present(a1d))
  36     d2 = l2i(present(a2d))
  37     d3 = l2i(present(a3d))
  38
  39     if(d1+d2+d3.ne.1)call crash('netcdf_read_array_wrf: must have exactly one of a1d a2d a3d arguments')
  40     ndims = d1+2*d2+3*d3
  41
  42     if(ndims.eq.1)then
  43         srf = 0
  44     elseif(present(sr))then
  45         srf = sr
  46     else
  47         call get_wrf_dims(ncid,srf)
  48     endif
  49
  50     allocate(star(ndims+1))
  51     allocate(cnts(ndims+1))
  52     allocate(dims(ndims+1))
  53     allocate(ends(ndims+1))
  54
  55     write(msg,*)"netcdf_read_array_wrf reading variable ",trim(name), &
  56           " dimension ",ndims," time step ",istep
  57     call message(msg)
  58
  59     ! get idx
  60     call netcdf_var_info(ncid,name,dims,varid,netcdf_msglevel)
  61     write(msg,*)"got dimensions ",trim(name),dims
  62     call message(msg)
  63     if(dims(ndims+1).eq.0)then
  64        write(msg,*)'netcdf_read_array_wrf: wrong dimension ndims=',ndims, ' 1dims=',dims
  65        call crash(msg)
  66     endif
  67
  68     star   = 1
  69     star(ndims+1) = istep
  70     ends(1:ndims) = dims(1:ndims)
  71     ends(ndims+1) = istep
  72     cnts = ends - star + 1
  73
  74     write(msg,*)"reading ",trim(name),star," to ",ends
  75     call message(msg)
  76
  77     ! read from file
  78     select case (ndims)
  79     case(1)
  80         allocate(a1(star(1):ends(1),star(2):ends(2)))
  81         allocate(a1d(star(1):ends(1)))
  82         call check(nf90_get_var(ncid, varid, a1, start = star, count = cnts),"nf90_get_var:"//trim(name))
  83         a1d = a1(:,istep)
  84         write(msg,*)"returning 1D array length ",shape(a1d)
  85     case(2)
  86         allocate(a2(star(1):ends(1),star(2):ends(2),star(3):ends(3)))
  87         call check(nf90_get_var(ncid, varid, a2, start = star, count = cnts),"nf90_get_var:"//trim(name))
  88         ! postprocessing - strip
  89         write(msg,*)" stripping ",srf," at ij ends"
  90         call message(msg)
  91         allocate(a2d(star(1):ends(1)-srf(1),star(2):ends(2)-srf(2)))
  92         do j=star(2),ends(2)-srf(2)
  93             do i=star(1),ends(1)-srf(1)
  94                 a2d(i,j) = a2(i,j,istep)
  95             enddo
  96         enddo
  97         write(msg,*)"returning array ij shape ",shape(a2d)
  98     case(3)
  99         allocate(a3(star(1):ends(1),star(2):ends(2),star(3):ends(3),star(4):ends(4)))
 100         call check(nf90_get_var(ncid, varid, a3, start = star, count = cnts),"nf90_get_var:"//trim(name))
 101         ! transpose at -> a and strip
 102         write(msg,*)" stripping ",srf," at ij ends and transposing to ikj indexing order"
 103         call message(msg)
 104         allocate(a3d(star(1):ends(1)-srf(1),star(3):ends(3),star(2):ends(2)-srf(2)))
 105         do k=star(3),ends(3)
 106             do j=star(2),ends(2)-srf(2)
 107                 do i=star(1),ends(1)-srf(1)
 108                     a3d(i,k,j) = a3(i,j,k,istep)
 109                 enddo
 110             enddo
 111         enddo
 112         write(msg,*)"returning array ikj shape ",shape(a3d)
 113     case default
 114         call crash('wrong case ndims')
 115     end select
 116     call message(msg)
 117 end subroutine netcdf_read_array_wrf
 118
 119 subroutine get_wrf_dims(ncid,sr,dims3d)
 120     implicit none
 121 !*** arguments
 122     integer, intent(in)::ncid    ! open wrfout dataset
 123     integer, intent(out), dimension(2):: sr  ! fire "subgrid" refinement factors in x and y directions
 124     integer, intent(out), dimension(3), optional:: dims3d ! grid dimensions
 125 !*** local
 126     integer, dimension(3) :: dims_atm, dims_fire   ! 2D + allow time dimension
 127     integer, dimension(4) :: dims4d   ! 3D + allow time dimension
 128     integer::varid, prints
 129     character(len=256)::msg
 130 !*** executable
 131     prints = netcdf_msglevel
 132     dims_atm = 0
 133     dims_fire = 0
 134     call netcdf_var_info(ncid,"XLAT",dims_atm,varid,prints=prints)
 135     call netcdf_var_info(ncid,"FXLAT",dims_fire,varid,prints=prints)
 136     sr = dims_fire(1:2)/(dims_atm(1:2) + 1)
 137     write(msg,*)"get_wrf_dims: fire subgrid refinement factors are ",sr
 138     call message(msg)
 139     if(present(dims3d))then
 140         call netcdf_var_info(ncid,"PHB",dims4d,varid,prints=prints)
 141         write(msg,*)"get_wrf_dims: netcdf dimensions ",dims4d
 142         call message(msg)
 143         dims3d = dims4d(1:3)
 144     endif
 145 end subroutine get_wrf_dims
 146
 147 subroutine write_fire_wind(filename,frame0_fmw,uf,vf,u_fmw,v_fmw,w_fmw,frame)
 148 implicit none
 149 !*** purpose
 150 ! write solution
 151 !*** arguments
 152 character(len=*), intent(in)::filename  ! open file
 153 real, intent(in), dimension(:,:)::uf,vf
 154 real, intent(in), dimension(:,:,:),optional::u_fmw,v_fmw,w_fmw
 155 integer, intent(in)::frame0_fmw   ! frame number
 156 integer, intent(in), optional::frame ! the default frame in the file
 157 !*** local
 158 integer::istep,chsum,ncid
 159 character(len=256)::msg
 160
 161 !*** executable
 162 istep=1
 163 if(present(frame))istep=frame
 164 write(msg,*)"writing timestep ",frame0_fmw," to ", trim(filename), " frame ",istep
 165 call message(msg)
 166 call ncopen(filename,nf90_write,ncid)
 167 call netcdf_write_2d(ncid,uf,"UF",istep)
 168 call netcdf_write_2d(ncid,vf,"VF",istep)
 169 if(present(u_fmw))call netcdf_write_wrf_3d(ncid,u_fmw,"U_FMW",istep) ! diagnostics
 170 if(present(v_fmw))call netcdf_write_wrf_3d(ncid,v_fmw,"V_FMW",istep)
 171 if(present(w_fmw))call netcdf_write_wrf_3d(ncid,w_fmw,"W_FMW",istep)
 172 chsum=ieor(get_chsum_2d(uf),get_chsum_2d(vf))
 173 print *,'chsum=',chsum
 174 call netcdf_write_int(ncid,chsum,"CHSUM_FMW")
 175 call netcdf_write_int(ncid,frame0_fmw,"FRAME_FMW")
 176 call ncclose(ncid)
 177
 178 end subroutine write_fire_wind
 179
 180
 181 subroutine read_initial_wind(filename,u0_fmw,v0_fmw,w0_fmw,frame0_fmw,frame)
 182 implicit none
 183 !*** purpose
 184 ! read wrf data, cycle until frame_fmw matches and chsum is correct
 185 !*** arguments
 186 character(len=*), intent(in)::filename  ! open file
 187 real, pointer, intent(out), dimension(:,:,:)::u0_fmw,v0_fmw,w0_fmw
 188 integer, intent(in)::frame0_fmw   ! frame number to expect
 189 integer, intent(in), optional::frame ! the default frame in the file to read, default=1
 190 ! return: 0=OK, >0 timed out
 191 !*** local
 192 integer::istep=1,chsum0,sr(2),chsum0_fmw,ierr=0,maxtry=200,frame0_in,itry,ncid
 193 character(len=256)::msg
 194 real:: waitmax=200
 195 integer::msec=100
 196
 197 !*** executable
 198 if(present(frame))istep=frame
 199 write(msg,*)"reading from file frame ",istep
 200 call message(msg)
 201 write(msg,*)"expecting time step frame0_fmw=",frame0_fmw
 202 call message(msg)
 203 call ncopen(filename,nf90_nowrite,ncid)
 204 do itry=1,maxtry
 205   frame0_in = netcdf_read_int_wrf(ncid,"FRAME0_FMW",istep)
 206   write(msg,*)"try ",itry," got ",frame0_in," expecting",frame0_fmw
 207   call message(msg)
 208   if(frame0_fmw .eq. frame0_in)goto 1
 209   call ncclose(ncid)
 210   if(frame0_fmw .eq. -99)then
 211        call message('received stop request frame=-99')
 212        stop
 213   endif
 214   if(itry * msec * 1e-3 .gt. waitmax)then
 215      write(msg,*)'timed out after ',waitmax,' sec waiting for frame ',frame0_fmw,' got ',frame0_in
 216      call crash(trim(msg))
 217   endif
 218   write(msg,*)'try ',itry,' waiting ',msec,' msec'
 219   call message(msg)
 220   call call_msleep(msec)
 221   call ncopen(filename,nf90_nowrite,ncid)
 222 enddo
 223 1 continue
 224
 225 call get_wrf_dims(ncid,sr) ! submesh refinement factors
 226 do itry=1,maxtry
 227   chsum0_fmw = netcdf_read_int_wrf(ncid,"CHSUM0_FMW",istep)
 228   write(msg,*)"read CHSUM0_FMW=",chsum0_fmw
 229   call message(msg)
 230
 231   call netcdf_read_array_wrf(ncid,"U0_FMW",istep,sr,a3d=u0_fmw)
 232   call netcdf_read_array_wrf(ncid,"V0_FMW",istep,sr,a3d=v0_fmw)
 233   call netcdf_read_array_wrf(ncid,"W0_FMW",istep,sr,a3d=w0_fmw)
 234
 235   chsum0 = get_chsum(u0_fmw)
 236   chsum0 = ieor(chsum0,get_chsum(v0_fmw))
 237   chsum0 = ieor(chsum0,get_chsum(w0_fmw))
 238   write(msg,*)" computed chsum0 ", chsum0
 239   call message(msg)
 240   call ncclose(ncid)
 241   if (chsum0_fmw.eq.chsum0)goto 2
 242   write(msg,*)'try ',itry,' waiting ',msec, ' rmsec'
 243   call message(msg)
 244   call call_msleep(msec)
 245   if(itry * msec * 1e-3 .gt. waitmax)then
 246      write(msg,*)'timed out after ',waitmax,' sec waiting for correct check sum'
 247      call crash(trim(msg))
 248   endif
 249   call ncopen(filename,nf90_nowrite,ncid)
 250 enddo
 251 call ncclose(ncid)
 252 2 continue
 253 write(msg,*)'success check sum match for time step ',frame0_fmw
 254 call message(msg)
 255
 256 end subroutine read_initial_wind
 257
 258 end module module_wrfout