Changed cova reading inequality constraint

constraint/readinput.f90

@@ -160,6 +160,11 @@ subroutine readinput()
   open(100,file=trim(dirinvert)//'/izwork.txt',action='read',status='old',iostat=ierr)
   if(ierr.ne.0) then
     write(logid,*) 'Error reading izwork.txt'
+    !CGZ: try reading cova.txt
+    open(100,file=trim(dirinvert)//'/cova.txt',action='read',status='old',iostat=ierr)
+    if(ierr.ne.0) then
+      write(logid,*) 'Error reading cova.txt'
+    endif
   endif
   do i=1,nvar
     read(100,*) cova(i,:)

prep_fluxes/convertgrid.f90

@@ -183,6 +183,9 @@ subroutine convertgrid(midpoints, lon_in, lat_in, flx_tmp, &
         endif
         mass = 0.
         aolaptot = 0.
+
+       ! print*, 'CGZ debug; ix_out', ix, ', jy_out', jy, 'ix_in', ix1, ix2, 'iy_in', jy1, jy2
+
         ! loop over grid cells in input grid which overlap with the grid cell in output grid
         do jj = jy1, jy2
           do ii = ix1, ix2
@@ -191,6 +194,8 @@ subroutine convertgrid(midpoints, lon_in, lat_in, flx_tmp, &
             xolap2 = min(wlon_in(ii)+dx_in(ii),wlon_out(ix)+dx_out)
             yolap1 = max(slat_in(jj),slat_out(jy))
             yolap2 = min(slat_in(jj)+dy_in(jj),slat_out(jy)+dy_out)
+            !print*, 'CGZ debug; xolap', xolap1, xolap2, ', yolap', yolap1, yolap2
+
             ! calculate overlapping area
             call gridarea((yolap1+yolap2)/2., abs(yolap2-yolap1), abs(xolap2-xolap1), aolap)
             mass = mass + aolap*flx_tmp(ii,jj,nt)

prep_fluxes/job_prep_flux.sh

 #!/bin/bash 
 #---------------------------------------------------
 
-file='./FLUXES_CH4_edgar5.txt'
+file='./FLUXES_CH4_VERIFY.txt'
 
 ./prep_flux ${file} 

prep_fluxes/read_ncdf.f90

@@ -137,7 +137,19 @@ subroutine read_ncdf(filename, varname, lonname, latname, timename,&
   endif
   if(trim(timename).ne.'none') then
     ierr=nf90_inq_varid(nid,timename,tid)
-    ierr=nf90_get_var(nid,tid,time(:),start=(/indextime/),count=(/ntime/))
+    if(ierr.ne.0) then
+      write(*,*) 'ERROR: reading variable time, timename not available'
+      stop
+    endif
+    ! ierr=nf90_inquire_dimension(nid,tid,len=nt)
+    ! if(ierr.ne.0) then
+    !   write(*,*) 'ERROR: inquire dimension variable time'
+    !   stop
+    ! else
+    !   write(*,*) 'Dimension of time variable;', nt, ntime
+    ! endif
+    !ierr=nf90_get_var(nid,tid,time(:),start=(/indextime/),count=(/ntime/))
+    ierr=nf90_get_var(nid,tid,time)
     if(ierr.ne.0) then
       write(*,*) 'ERROR: reading variable time'
       stop
@@ -153,7 +165,13 @@ subroutine read_ncdf(filename, varname, lonname, latname, timename,&
   ierr=nf90_get_var(nid,varid,flx(:,:,:),start=(/1,1,indextime/),count=(/nlon,nlat,ntime/))
   if(ierr.ne.0) then
     write(*,*) 'ERROR reading variable ',trim(varname)
-    stop
+    write(*,*) 'Test reading complete array'
+    ierr=nf90_get_var(nid,varid,flx)
+    if(ierr.ne.0) then
+      write(*,*) 'ERROR reading variable ',trim(varname)
+      stop
+    endif
+  
   endif
  
   ! read fill value attribute

prep_fluxes/sbatch_prep_flux.sh

 #!/bin/bash 
 #---------------------------------------------------
-file='./FLUXES_ghg'
+file='./FLUXES_CH4_VERIFY.txt
 #---------------------------------------------------
 
 cat <<EOF > run_job.sh

prep_regions/mod_strings.f90

@@ -117,11 +117,16 @@ module mod_strings
       replen = len_trim(repstr)
       ! find first position of substr in string_in 
       n = index( string_in, substr, back=.false. )
-      strtmp1 = string_in(1:n-1)
-      ! find last position of substr in string_in 
-      n = index( string_in, substr, back=.true. )      
-      strtmp2 = string_in(n+replen:strlen)
-      string_out = trim(strtmp1)//trim(repstr)//trim(strtmp2) 
+      if (n.gt.0) then
+        strtmp1 = string_in(1:n-1)
+        ! find last position of substr in string_in 
+        n = index( string_in, substr, back=.true. )      
+        strtmp2 = string_in(n+replen:strlen)
+        string_out = trim(strtmp1)//trim(repstr)//trim(strtmp2) 
+      else
+        write(*,*) 'YYYY in flux at 0 position?! Assume no yyyy in filename.'
+        string_out=string_in
+      endif
 
     end function str_replace
 

prep_regions/read_reclist.f90

@@ -65,11 +65,17 @@ subroutine read_reclist(filename)
   open(100,file=trim(filename),action='read',status='old',iostat=ierr)
   write(*,*) 'Receptors: '
 
-  do cnt = 1, nrec
-    read(100,fmt='(A4,1X,F6.2,1X,F6.2,1X,F7.2)') recname(cnt), reclat(cnt), reclon(cnt), recalt(cnt)
-    write(*,*) recname(cnt), reclat(cnt), reclon(cnt), recalt(cnt)
+!  do cnt = 1, nrec
+!    read(100,fmt='(A4,1X,F6.2,1X,F6.2,1X,F7.2)') recname(cnt), reclat(cnt), reclon(cnt), recalt(cnt)
+!    write(*,*) recname(cnt), reclat(cnt), reclon(cnt), recalt(cnt)
 !    read(100,*) recname(cnt)
 !    write(*,*) recname(cnt)
+!  end do
+
+  !CGZ; test old version
+  do cnt = 1, nrec
+    read(100,*) recname(cnt)
+    write(*,*) recname(cnt)
   end do
 
   close(100)

source/mod_analytic.f90

@@ -85,10 +85,7 @@ module mod_analytic
     if ( nvar.le.nobs ) then
       call invert_nvar(files, config, obs, states, covar)
     else
- !CGZ for testing of constraint need izwork; use invert_nvar
-      call invert_nvar(files, config, obs, states, covar)
-!      call invert_nobs(files, config, obs, states, covar)
-
+      call invert_nobs(files, config, obs, states, covar)
     endif
 
     call calc_cost(config, obs, states, cost_o)
@@ -560,6 +557,15 @@ module mod_analytic
     end do
     close(100)
 
+    filename = trim(files%path_output)//'cova.txt'
+    open(100,file=trim(filename),status='replace',action='write',iostat=ierr)
+    write(rowfmt,'(A,I6,A)') '(',nvar,'(E11.4,1X))'
+    do i = 1, nvar
+      write(100,rowfmt) cova(i,:)
+    end do
+    close(100)
+
+
     deallocate(gain)
     deallocate(imwork)
     deallocate(cova)

source/mod_flexpart.f90

@@ -55,7 +55,7 @@ module mod_flexpart
 
     implicit none
 
-    character(len=max_path_len),    intent(in)     :: filename
+    character(*),                   intent(in)     :: filename
     integer,                        intent(in out) :: numpoint
     integer,                        intent(in out) :: ibdate, ibtime   
     real(kind=8), dimension(maxpoint,2), intent(in out) :: releases           

source/read_ncdf.f90

@@ -97,9 +97,9 @@ subroutine read_ncdf(filename, varname, lonname, latname, timename,&
   xres = lon(2) - lon(1)
   yres = lat(2) - lat(1)
 
-  !write(logid,*) 'read_ncdf: lono, lato = ',lono,lato
-  !write(logid,*) 'read_ncdf: xres, yres = ',xres,yres 
-  !write(logid,*) 'read_ncdf: lon(1), lat(1) = ',lon(1),lat(1) 
+  write(logid,*) 'read_ncdf: lono, lato = ',lono,lato
+  write(logid,*) 'read_ncdf: xres, yres = ',xres,yres 
+  write(logid,*) 'read_ncdf: lon(1), lat(1) = ',lon(1),lat(1) 
   write(logid,*) 'cgz read_ncdf: varname = ',varname
 
   if ( ((lon(1)-xres/2.).gt.lono).or.((lat(1)-yres/2.).gt.lato) ) then

source/simulate.f90

@@ -247,27 +247,6 @@ subroutine simulate(iter, files, config, fluxes, obs, states, grad_o, cost_o)
         write(logid,*) 'WARNING: cannot find '//trim(filename)
         go to 10
       endif
-<<<<<<< HEAD
-      call system_clock(tread2, countrate)
-      
-      write(logid,*) 'simulate: time to read grid (s) = ',(tread2-tread1)/countrate
-      ! correction for dry air (only if obs are mixing ratios) 
-!      do n = 1, ngrid
-!        ind = minloc(ftime,dim=1,mask=abs(ftime-gtime(n)).eq.minval(abs(ftime-gtime(n))))
-!        grid(:,:,n) = grid(:,:,n)*factor(:,:,ind)
-!      end do
-      ! convert from ppt to ppmv or ppbv
-      grid = grid*config%coeff*mmair/config%molarmass
-      ! convert s.m3/kg to s.m2/kg
-      grid = grid/outheight(1)
-      ! apply numerical scaling
-      grid = grid/numscale
-    else
-      write(logid,*) 'WARNING: cannot find '//trim(filename)
-      go to 10
-    endif
-    
-=======
     endif
     call system_clock(tread2, countrate)
 !      write(logid,*) 'simulate: time to read grid (s) = ',(tread2-tread1)/countrate
@@ -285,7 +264,6 @@ subroutine simulate(iter, files, config, fluxes, obs, states, grad_o, cost_o)
     grid = grid/numscale
 
     ! read nested flexpart footprints
->>>>>>> master
     if ( config%nested ) then
       if ( .not.config%average_fp ) then
         ! footprint times match observation timestamps