var/da/da_obs_io/da_write_obs_chem_sfc.inc

   1 subroutine da_write_obs_chem_sfc(it,ob, iv, re)
   2
   3    !-------------------------------------------------------------------------
   4    ! Purpose: Writes out components of iv=O-B structure.
   5    !-------------------------------------------------------------------------
   6
   7    implicit none
   8
   9    integer,        intent(in)    :: it
  10    type (y_type),  intent(in)    :: ob      ! Observation structure.
  11    type (iv_type), intent(in)    :: iv      ! O-B structure.
  12    type (y_type),  intent(inout) :: re      ! residual vector.
  13
  14    integer                     :: n, k, num_obs, ios, ios2
  15    integer                     :: ounit, ounit2     ! Output unit
  16    character(len=filename_len) :: filename, filename2
  17
  18    character(len=120) :: fmt_chem1 = '(i8,2x,a6,2f11.6,2f11.2,i8,2f11.2)'     ! pm2.5 or pm10
  19    character(len=120) :: fmt_chem2 = '(i8,2x,a6,2f11.6,2(2f11.2,i8,2f11.2))'  ! pm2.5 and pm10
  20    character(len=120) :: fmt_chem4 = '(i8,2x,a6,2f11.6,4(2f12.3,i8,2f12.3))'  ! 4 gas species
  21    integer :: itime, ifgat
  22    integer :: ipcc, jpcc, ic
  23    integer, parameter       :: ngas = 4
  24    integer, dimension(ngas) :: ipc
  25    ipc = (/p_chemsi_so2, p_chemsi_no2, p_chemsi_o3, p_chemsi_co/)
  26
  27    if (trace_use) call da_trace_entry("da_write_obs_chem_sfc")
  28
  29    !-------------------------------------------------------------------------
  30    ! Fix output unit
  31    !-------------------------------------------------------------------------
  32
  33    if (chemicda_opt == 5) then
  34      call da_get_unit(ounit)
  35      call da_get_unit(ounit2)
  36    else
  37      call da_get_unit(ounit)
  38    end if
  39
  40 #ifdef DM_PARALLEL
  41    if (chemicda_opt == 5) then
  42       write(unit=filename, fmt='(a,i2.2,a,i4.4)') 'chem_omb_oma_',it,'.', myproc
  43       write(unit=filename2, fmt='(a,i2.2,a,i4.4)') 'gas_omb_oma_',it,'.', myproc
  44    else if (chemicda_opt == 4) then
  45       write(unit=filename, fmt='(a,i2.2,a,i4.4)') 'gas_omb_oma_',it,'.', myproc
  46    else
  47       write(unit=filename, fmt='(a,i2.2,a,i4.4)') 'chem_omb_oma_',it,'.', myproc
  48    end if
  49 #else
  50    if (chemicda_opt == 5) then
  51       write(unit=filename, fmt='(a,i2.2,a)') 'chem_omb_oma_',it,'.0000'
  52       write(unit=filename2, fmt='(a,i2.2,a)') 'gas_omb_oma_',it,'.0000'
  53    else if (chemicda_opt == 4) then
  54       write(unit=filename, fmt='(a,i2.2,a)') 'gas_omb_oma_',it,'.0000'
  55    else
  56       write(unit=filename, fmt='(a,i2.2,a)') 'chem_omb_oma_',it,'.0000'
  57    end if
  58 #endif
  59
  60    open (unit=ounit,file=trim(filename),form='formatted',status='replace', &
  61       iostat=ios)
  62    if (ios /= 0) then
  63       call da_error(__FILE__,__LINE__, &
  64          (/"Cannot open chemical observation omb and oma file"//filename/))
  65    end if
  66
  67    if (chemicda_opt == 5) then
  68    open (unit=ounit2,file=trim(filename2),form='formatted',status='replace', &
  69       iostat=ios2)
  70    if (ios2 /= 0) then
  71       call da_error(__FILE__,__LINE__, &
  72          (/"Cannot open gas observation omb and oma file"//filename2/))
  73    end if
  74    end if
  75
  76    num_obs = 0
  77    do n = 1, iv%info(chemic_surf)%nlocal
  78       if (iv%info(chemic_surf)%proc_domain(1,n)) num_obs = num_obs + 1
  79    end do
  80    if (num_obs > 0) then
  81       if (chemicda_opt == 4) then
  82           write(ounit,'(a20,i8)') 'gas', num_obs
  83       else
  84           write(ounit,'(a20,i8)') 'chem', num_obs
  85       endif
  86       if (chemicda_opt == 5) then
  87           write(ounit2,'(a20,i8)') 'gas', num_obs
  88       end if
  89       num_obs = 0
  90       do n = 1, iv%info(chemic_surf)%nlocal
  91          do itime = 1, num_fgat_time
  92             if ( n >= iv%info(chemic_surf)%plocal(itime-1)+1 .and. &
  93                  n <= iv%info(chemic_surf)%plocal(itime) ) then
  94                ifgat = itime
  95                exit
  96             end if
  97          end do
  98
  99          if (iv%info(chemic_surf)%proc_domain(1,n)) then
 100             num_obs = num_obs + 1
 101              write(ounit,'(2i8)') 1, ifgat
 102              if (chemicda_opt == 5) then
 103                write(ounit2,'(2i8)') 1, ifgat
 104              end if
 105
 106             ipcc = p_chemsi_pm25
 107             jpcc = p_chemsi_pm10
 108             if (chemicda_opt <= 2) then
 109                if(chemicda_opt == 2) ipcc = p_chemsi_pm10
 110                write(ounit, fmt = fmt_chem1)         &
 111                num_obs , iv%info(chemic_surf)%id(n), &  ! Station
 112                iv%info(chemic_surf)%lat(1,n),        &  ! Latitude
 113                iv%info(chemic_surf)%lon(1,n),        &  ! Longitude
 114                ob%chemic_surf(n)%chem(ipcc),         &  ! observation
 115                iv%chemic_surf(n)%chem(ipcc)%inv,     &  ! o-b
 116                iv%chemic_surf(n)%chem(ipcc)%qc,      &  ! qc
 117                iv%chemic_surf(n)%chem(ipcc)%error,   &  ! obs error
 118                re%chemic_surf(n)%chem(ipcc)             ! o-a
 119             else if (chemicda_opt == 4) then
 120                write(ounit, fmt = fmt_chem4)         &
 121                num_obs , iv%info(chemic_surf)%id(n), &  ! Station
 122                iv%info(chemic_surf)%lat(1,n),        &  ! Latitude
 123                iv%info(chemic_surf)%lon(1,n),        &  ! Longitude
 124               (ob%chemic_surf(n)%chem(ipc(ic)),      &  ! observation
 125                iv%chemic_surf(n)%chem(ipc(ic))%inv,  &  ! o-b
 126                iv%chemic_surf(n)%chem(ipc(ic))%qc,   &  ! qc
 127                iv%chemic_surf(n)%chem(ipc(ic))%error,&  ! obs error
 128                re%chemic_surf(n)%chem(ipc(ic)), ic=1,ngas) ! o-a
 129             else  ! chemicda_opt == 3 or 5
 130                write(ounit, fmt = fmt_chem2)         &
 131                num_obs , iv%info(chemic_surf)%id(n), &  ! Station
 132                iv%info(chemic_surf)%lat(1,n),        &  ! Latitude
 133                iv%info(chemic_surf)%lon(1,n),        &  ! Longitude
 134                ob%chemic_surf(n)%chem(ipcc),         &  ! observation
 135                iv%chemic_surf(n)%chem(ipcc)%inv,     &  ! o-b
 136                iv%chemic_surf(n)%chem(ipcc)%qc,      &  ! qc
 137                iv%chemic_surf(n)%chem(ipcc)%error,   &  ! obs error
 138                re%chemic_surf(n)%chem(ipcc),         &  ! o-a
 139                ob%chemic_surf(n)%chem(jpcc),         &  ! observation
 140                iv%chemic_surf(n)%chem(jpcc)%inv,     &  ! o-b
 141                iv%chemic_surf(n)%chem(jpcc)%qc,      &  ! qc
 142                iv%chemic_surf(n)%chem(jpcc)%error,   &  ! obs error
 143                re%chemic_surf(n)%chem(jpcc)             ! o-a
 144               if(chemicda_opt == 5) then
 145                write(ounit2, fmt = fmt_chem4)        &
 146                num_obs , iv%info(chemic_surf)%id(n), &  ! Station
 147                iv%info(chemic_surf)%lat(1,n),        &  ! Latitude
 148                iv%info(chemic_surf)%lon(1,n),        &  ! Longitude
 149               (ob%chemic_surf(n)%chem(ipc(ic)),      &  ! observation
 150                iv%chemic_surf(n)%chem(ipc(ic))%inv,  &  ! o-b
 151                iv%chemic_surf(n)%chem(ipc(ic))%qc,   &  ! qc
 152                iv%chemic_surf(n)%chem(ipc(ic))%error,&  ! obs error
 153                re%chemic_surf(n)%chem(ipc(ic)), ic=1,ngas) ! o-a
 154               endif  ! (chemicda_opt == 5)
 155             endif
 156
 157          end if      ! if (iv%info(chemic_surf)%proc_domain(1,n)) then
 158       end do         ! do n = 1, iv%info(chemic_surf)%nlocal
 159    end if            ! if (num_obs > 0) then
 160
 161    close (ounit)
 162    call da_free_unit(ounit)
 163
 164    if (chemicda_opt == 5) then
 165      close (ounit2)
 166      call da_free_unit(ounit2)
 167    end if
 168
 169    if (trace_use) call da_trace_exit("da_write_obs_chem_sfc")
 170
 171 end subroutine da_write_obs_chem_sfc