123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707 |
- ! File: m_prep_4_EnKF.F90
- !
- ! Created: ???
- !
- ! Last modified: 29/06/2010
- !
- ! Purpose: Calculation of HA ("S")
- !
- ! Description: Calculates HA by going sequentially through each data type.
- !
- ! Modifications:
- ! 09/11/2012 Geir Arne Waagbo:
- ! - Added support for OSISAF ice drift obs
- ! 29/07/2010 PS:
- ! - merged insitu_QC() with generic obs_QC(). Moved
- ! insitu_writeforecast() to the point after QC.
- ! 29/06/2010 PS:
- ! - added generic observation QC: increase the observation
- ! error when observation and ensemble mean are much too far
- ! away than expected
- ! Prior history:
- ! Not documented.
- module m_prep_4_EnKF
- integer, parameter, private :: STRLEN = 512
- private read_mean_ssh
- contains
- ! This subroutine uses the observation and ensembles from the model
- ! to prepare the input to the EnKF analysis scheme.
- ! The output from this routine is used directly in the global analysis
- ! while the output has to be run through a "filter" to be used in the
- ! local analysis scheme.
- ! S = HA (ensemble observation anomalies)
- ! d = d - Hx (innovations)
- !
- ! S is calculated in two steps:
- ! 1. S = HE
- ! 2. S = S - repmat(s, 1, m),
- ! where s = mean(S')';
- ! Note that in reality (with HYCOM) H is different for each member...
- ! So that HX must be read "HX" rather than "H * X".
- !
- subroutine prep_4_EnKF(nrens, enslist, d, S, depths, meandx, nx, ny, nz)
- #if defined (QMPI)
- use qmpi, only : master, stop_mpi
- #else
- use qmpi_fake, only : master, stop_mpi
- #endif
- use mod_measurement
- use m_obs
- use m_Generate_element_Si
- use m_get_mod_fld
- use m_read_icemod
- use m_parameters
- implicit none
- integer, intent(in) :: nx, ny, nz ! Model size
- integer, intent(in) :: nrens ! Size of ensemble
- integer, dimension(:),intent(in) :: enslist ! [CKB,FM] List of existing ens members
- real, intent(in) :: depths(nx, ny)
- real, intent(in) :: meandx ! mean grid size
- real, intent(inout) :: d(nobs)
- real, intent(inout) :: S(nobs, nrens)
- real :: x(nobs)
- integer :: i, j, m, iens
- real*4, dimension(nx,ny) :: fldr4
- real :: readfld(nx, ny), ai1(nx,ny), ai2(nx,ny), ai3(nx,ny), ai4(nx,ny), ai5(nx,ny), uice(nx,ny), vice(nx,ny)
- real :: vi1(nx,ny), vi2(nx,ny), vi3(nx,ny), vi4(nx,ny), vi5(nx,ny)
- real :: vs1(nx,ny), vs2(nx,ny), vs3(nx,ny), vs4(nx,ny), vs5(nx,ny)
- ! hard-coded for now
- integer, parameter :: drnx = 152, drny = 132
- real*4, dimension(drnx, drny) :: modzon, modmer
- integer, parameter :: drnx_osisaf = 119, drny_osisaf = 177
- real*4, dimension(drnx_osisaf, drny_osisaf) :: dX, dY
- integer :: reclSLA, ios, reclDRIFT
- character*3 :: cmem
- character*2 :: day
- character*1 :: offset
- logical :: ex
- character(STRLEN) :: fname
- integer :: iuobs
- ! FANF: For track assim we launch m_Generate_Si for each day (t=1:Wd)
- ! which fills in S at the approriate indices.
- ! Wd is is the assimilation cycle (i.e. 7 days)
- !
- integer :: Wd, t
- integer :: tlevel
- real :: field2(nx, ny), field3(nx, ny) ! auxiliary fields (e.g. mean SSH,
- ! field bias estimate etc.)
- integer :: nthisobs, thisobs(nobs)
- if (any(obs(:) % id == 'TSLA ')) then
- Wd = 6
- else
- Wd = 0
- endif
- ! security check
- !
- if (any(obs(:) % id == 'SSH ') .or. any(obs(:) % id == 'SLA ')) then
- if (any(obs(:) % id == 'SLA ')) then
- inquire(exist = ex, file = 'model_SLA.uf')
- if (.not.ex) then
- if (master) print *,'model_SLA.uf does not exist'
- call stop_mpi()
- end if
- end if
- if (any(obs(:) % id == 'SSH ')) then
- inquire(exist = ex, file = 'model_SSH.uf')
- if (.not.ex) then
- if (master) print *,'model_SSH.uf does not exist'
- call stop_mpi()
- end if
- end if
- end if
- ! construct S=HA
- !
- do iuobs = 1, nuobs
- if (master) then
- print *, 'prep_4_EnKF: now preparing "', trim(unique_obs(iuobs)), '" observations'
- end if
- if (trim(unique_obs(iuobs)) == 'ICEC') then
- do iens = 1, nrens
- write(cmem,'(i3.3)') iens
- tlevel = 1
- call get_mod_fld_new(trim('forecast'//cmem), readfld, iens,&
- 'icec', 0, tlevel, nx, ny)
- if (tlevel == -1) then
- if (master) then
- print *, 'ERROR: get_mod_fld_new(): failed for "icec"'
- end if
- stop
- end if
- call Generate_element_Si(S(:, iens), unique_obs(iuobs),&
- readfld, depths, nx, ny, nz, 0)
- end do
- ! [FM, May 2013: for LIM3 sea ice model]
- elseif (trim(unique_obs(iuobs)) == 'AT_I') then
- do iens = 1, nrens
- write(cmem,'(i3.3)') iens
- tlevel = 1
- call io_mod_fld(ai1,iens,enslist, &
- 'a_i_htc1', 1, 0, 1, nx,ny, 'get',FLOAT(obs(1)%date))
- call io_mod_fld(ai2,iens,enslist, &
- 'a_i_htc2', 1, 0, 1, nx,ny, 'get',FLOAT(obs(1)%date))
- call io_mod_fld(ai3,iens,enslist, &
- 'a_i_htc3', 1, 0, 1, nx,ny, 'get',FLOAT(obs(1)%date))
- call io_mod_fld(ai4,iens,enslist, &
- 'a_i_htc4', 1, 0, 1, nx,ny, 'get',FLOAT(obs(1)%date))
- call io_mod_fld(ai5,iens,enslist, &
- 'a_i_htc5', 1, 0, 1, nx,ny, 'get',FLOAT(obs(1)%date))
- if (tlevel == -1) then
- if (master) then
- print *, 'ERROR: io_mod_fld_new(): failed for "at_i"'
- end if
- stop
- end if
- ! Multipply by 100 to match obs conventions
- readfld=(ai1+ai2+ai3+ai4+ai5) * 100.0
- call Generate_element_Si(S(:, iens), unique_obs(iuobs),&
- readfld, depths, nx, ny, nz, 0)
- end do
- ! freeboard
- elseif(trim(unique_obs(iuobs)) == 'VT_I') then
- do iens = 1, nrens
- write(cmem, '(i3.3)') iens
- tlevel = 1
- call io_mod_fld(ai1,iens,enslist, &
- 'a_i_htc1', 1, 0, 1, nx,ny, 'get',FLOAT(obs(1)%date))
- call io_mod_fld(ai2,iens,enslist, &
- 'a_i_htc2', 1, 0, 1, nx,ny, 'get',FLOAT(obs(1)%date))
- call io_mod_fld(ai3,iens,enslist, &
- 'a_i_htc3', 1, 0, 1, nx,ny, 'get',FLOAT(obs(1)%date))
- call io_mod_fld(ai4,iens,enslist, &
- 'a_i_htc4', 1, 0, 1, nx,ny, 'get',FLOAT(obs(1)%date))
- call io_mod_fld(ai5,iens,enslist, &
- 'a_i_htc5', 1, 0, 1, nx,ny, 'get',FLOAT(obs(1)%date))
- call io_mod_fld(vi1,iens,enslist, &
- 'v_i_htc1', 1, 0, 1, nx,ny, 'get',FLOAT(obs(1)%date))
- call io_mod_fld(vi2,iens,enslist, &
- 'v_i_htc2', 1, 0, 1, nx,ny, 'get',FLOAT(obs(1)%date))
- call io_mod_fld(vi3,iens,enslist, &
- 'v_i_htc3', 1, 0, 1, nx,ny, 'get',FLOAT(obs(1)%date))
- call io_mod_fld(vi4,iens,enslist, &
- 'v_i_htc4', 1, 0, 1, nx,ny, 'get',FLOAT(obs(1)%date))
- call io_mod_fld(vi5,iens,enslist, &
- 'v_i_htc5', 1, 0, 1, nx,ny, 'get',FLOAT(obs(1)%date))
- if (tlevel == -1) then
- if (master) then
- print *, 'ERROR: io_mod_fld_nex(): failed for "SIFB"'
- end if
- stop
- end if
- readfld=(vi1+vi2+vi3+vi4+vi5)
- call Generate_element_Si(S(:, iens), unique_obs(iuobs),&
- readfld, depths, nx, ny, nz, 0)
- end do
- ! freeboard
- elseif(trim(unique_obs(iuobs)) == 'RFB') then
- do iens = 1, nrens
- write(cmem, '(i3.3)') iens
- tlevel = 1
- call io_mod_fld(ai1,iens,enslist, &
- 'a_i_htc1', 1, 0, 1, nx,ny, 'get',FLOAT(obs(1)%date))
- call io_mod_fld(ai2,iens,enslist, &
- 'a_i_htc2', 1, 0, 1, nx,ny, 'get',FLOAT(obs(1)%date))
- call io_mod_fld(ai3,iens,enslist, &
- 'a_i_htc3', 1, 0, 1, nx,ny, 'get',FLOAT(obs(1)%date))
- call io_mod_fld(ai4,iens,enslist, &
- 'a_i_htc4', 1, 0, 1, nx,ny, 'get',FLOAT(obs(1)%date))
- call io_mod_fld(ai5,iens,enslist, &
- 'a_i_htc5', 1, 0, 1, nx,ny, 'get',FLOAT(obs(1)%date))
- call io_mod_fld(vi1,iens,enslist, &
- 'v_i_htc1', 1, 0, 1, nx,ny, 'get',FLOAT(obs(1)%date))
- call io_mod_fld(vi2,iens,enslist, &
- 'v_i_htc2', 1, 0, 1, nx,ny, 'get',FLOAT(obs(1)%date))
- call io_mod_fld(vi3,iens,enslist, &
- 'v_i_htc3', 1, 0, 1, nx,ny, 'get',FLOAT(obs(1)%date))
- call io_mod_fld(vi4,iens,enslist, &
- 'v_i_htc4', 1, 0, 1, nx,ny, 'get',FLOAT(obs(1)%date))
- call io_mod_fld(vi5,iens,enslist, &
- 'v_i_htc5', 1, 0, 1, nx,ny, 'get',FLOAT(obs(1)%date))
- call io_mod_fld(vs1,iens,enslist, &
- 'v_s_htc1', 1, 0, 1, nx,ny, 'get',FLOAT(obs(1)%date))
- call io_mod_fld(vs2,iens,enslist, &
- 'v_s_htc2', 1, 0, 1, nx,ny, 'get',FLOAT(obs(1)%date))
- call io_mod_fld(vs3,iens,enslist, &
- 'v_s_htc3', 1, 0, 1, nx,ny, 'get',FLOAT(obs(1)%date))
- call io_mod_fld(vs4,iens,enslist, &
- 'v_s_htc4', 1, 0, 1, nx,ny, 'get',FLOAT(obs(1)%date))
- call io_mod_fld(vs5,iens,enslist, &
- 'v_s_htc5', 1, 0, 1, nx,ny, 'get',FLOAT(obs(1)%date))
- if (tlevel == -1) then
- if (master) then
- print *, 'ERROR: io_mod_fld_nex(): failed for "SIFB"'
- end if
- stop
- end if
- readfld=(((vi1+vi2+vi3+vi4+vi5) * (1024.0 - 899.5) - 330 * (vs1+vs2+vs3+vs4+vs5)) / &
- 1024.0-0.25*(vs1 +vs2+vs3+vs4+vs5))
- !readfld=(((vi1+vi2+vi3+vi4+vi5) * (1024.0 - 899.5) - 330 * (vs1+vs2+vs3+vs4+vs5)) / 1024.0 - 0.25 * (vs1+vs2+vs3+vs4+vs5)) / (ai1+ai2+ai3+ai4+ai5)
-
- ! Conversion of models' sea ice thickness and snow thickness to
- ! model's freeboard using fixed densities for snow (330 kg/m3), ice
- ! (899.5 kg/m3 = average of MYI and FYI from Guerreiro et al. 2017
- ! and seawater (1024 kg/m3). The model freeboard is then lowered by
- ! 25% of the snow depth to account for the fact that the radar
- ! measurement underestimates the actual freeboard due to the lower
- ! propagation speed of the wave into the snow than in the air.
- ! Everything is converted from grid cell mean to in situ by
- ! dividing by concentration (if it is not zero). See exchanges
- ! e-mail with Sara Fleury 7 December 2020.
- call Generate_element_Si(S(:, iens), unique_obs(iuobs),&
- readfld, depths, nx, ny, nz, 0)
- end do
- elseif (trim(unique_obs(iuobs)) == 'SST') then
- do iens = 1, nrens
- write(cmem,'(i3.3)') iens
- tlevel = 1
- call get_mod_fld_new(trim('forecast'//cmem), readfld, iens,&
- 'tn', 1, tlevel, nx, ny)
- PRINT *, "FRANCOIS"
- if (tlevel == -1) then
- if (master) then
- print *, 'ERROR: get_mod_fld_new(): failed for "SST"'
- end if
- stop
- end if
- if (prm_prmestexists('sstb')) then
- tlevel = 1
- call get_mod_fld_new(trim('forecast'//cmem), field2, iens,&
- 'sstb', 0, tlevel, nx, ny)
- if (tlevel == -1) then
- if (master) then
- print *, 'ERROR: get_mod_fld_new(): failed for "sstb"'
- end if
- stop
- end if
- readfld = readfld - field2
- end if
- call Generate_element_Si(S(:, iens), unique_obs(iuobs),&
- readfld, depths, nx, ny, nz, 0)
- end do
- elseif (trim(unique_obs(iuobs)) == 'SLA' .or.&
- trim(unique_obs(iuobs)) == 'TSLA') then
- if (trim(unique_obs(iuobs)) == 'TSLA') then
- call read_mean_ssh(field2, nx, ny)
- end if
-
- inquire(iolength=reclSLA) fldr4
- ! FANF loop over each day of the week
- do t = 0, Wd
- if (trim(unique_obs(iuobs)) == 'TSLA') then
- write(day,'(i2.2)') t
- fname = trim('model_TSSH_'//day//'.uf')
- else
- fname = 'model_SLA.uf'
- endif
- if (master) then
- print *, 'TSLA, day', t, ': nobs = ',&
- count(obs(uobs_begin(iuobs) : uobs_end(iuobs)) % date == t)
- end if
- do iens = 1, nrens
- open(10, file = trim(fname), access = 'direct',&
- status = 'old', recl = reclSLA, action = 'read')
- read(10, rec = iens, iostat = ios) fldr4
- if (ios /= 0) then
- if (master) print *, 'Error reading ', trim(fname), ', member #', iens
- call stop_mpi()
- end if
- close(10)
- readfld = fldr4
-
- if (prm_prmestexists('msshb')) then
- write(cmem,'(i3.3)') iens
- tlevel = 1
- call get_mod_fld_new(trim('forecast'//cmem), field3, iens,&
- 'msshb', 0, tlevel, nx, ny)
- if (tlevel == -1) then
- if (master) then
- print *, 'ERROR: get_mod_fld_new(): failed for "msshb"'
- end if
- stop
- end if
- readfld = readfld - field3 ! mean SSH bias for this member
- end if
- if (trim(unique_obs(iuobs)) == 'TSLA') then
- readfld = readfld - field2 ! mean SSH
- end if
-
- call Generate_element_Si(S(:, iens), unique_obs(iuobs),&
- readfld, depths, nx, ny, nz, t)
- end do
- if (master) then
- print *, 'forming S, day', t
- print *, ' # of non-zero ens observations = ', count(S /= 0.0)
- print *, ' # of zero ens observations = ', count(S == 0.0)
- print *, ' # of non-zero observations for member 1 = ', count(S(:, 1) /= 0.0)
- end if
- end do
- elseif (trim(unique_obs(iuobs)) == 'SAL' .or.&
- trim(unique_obs(iuobs)) == 'TEM' .or.&
- trim(unique_obs(iuobs)) == 'GSAL' .or.&
- trim(unique_obs(iuobs)) == 'GTEM') then
- if (master) then
- print *, ' Interpolating ensemble vectors to the locations of "',&
- trim(unique_obs(iuobs)), '" observations'
- end if
- !
- ! for each ensemble member process all profiles "in parallel",
- ! reading the fields layer by layer
- !
- do iens = 1, nrens
- call get_S(S(:, iens), trim(unique_obs(iuobs)), nobs, obs, iens)
- end do
- if (master) then
- print *, ' Interpolation completed'
- end if
-
- elseif ((trim(unique_obs(iuobs)) == 'U_ICE') .or. trim(unique_obs(iuobs)) == 'V_ICE') then
- do iens = 1, nrens
- ! [FM] Read the file
- !inquire(iolength=reclDRIFT) modzon, modmer
- !open(10, file = 'model_ICEDRIFT.uf', access = 'direct',&
- ! status = 'old', recl = reclDRIFT, action = 'read')
- !read(10, rec = iens, iostat = ios) modzon, modmer
- !close(10)
- !if (ios /= 0) then
- ! if (master) then
- ! print *,'ERROR: could not read ensemble ice drift for member ', iens
- ! end if
- ! call stop_mpi()
- !end if
-
- call io_mod_fld(uice,iens,enslist, &
- 'u_ice', 1, 0, 1, nx,ny, 'get',FLOAT(obs(1)%date))
- call io_mod_fld(vice,iens,enslist, &
- 'v_ice', 1, 0, 1, nx,ny, 'get',FLOAT(obs(1)%date))
- do m = 1, nobs
- i = obs(m) % i_orig_grid
- j = obs(m) % j_orig_grid
- if (trim(obs(m) % id) == 'U_ICE') then
- S(m, iens) = uice(i, j)
- elseif (trim(obs(m) % id) == 'V_ICE') then
- S(m, iens) = vice(i, j)
- end if
- end do
- end do
-
- elseif ((trim(unique_obs(iuobs)) == 'U2D_I') .OR. trim(unique_obs(iuobs)) == 'V2D_I' ) THEN
- ! ADDED BY FM FRANCOIS MASSONNET. u_ice_2d or v_ice_2d is the sea ice u or v-velocity
- ! obtained as follows:
- ! 1) Rotate OSISAF Low resolution 2-day sea ice drift in a {east,north}
- ! reference frame
- ! 2) Interpolate to the ORCA grid
- ! 3) Rotate to align with the ORCA grid
- ! 4) Multiply by 1000 and divide by 2*86400 to convert the 2-day
- ! displacement from km to m/s
- DO iens=1,nrens
- CALL read_icemod(uice,iens,enslist,'iicevelu',nx,ny)
- CALL read_icemod(vice,iens,enslist,'iicevelv',nx,ny)
- DO m = 1, nobs
- i = obs(m) % i_orig_grid
- j = obs(m) % j_orig_grid
-
- IF (trim(obs(m) % id) == 'U2D_I') THEN
- S(m,iens) = uice(i,j)
- ELSEIF (trim(obs(m) % id) == 'V2D_I') THEN
- S(m,iens) = vice(i,j)
- END IF
- END DO ! nobs
- END DO ! iens
-
- elseif ((index(unique_obs(iuobs),'DX') > 0 ) .or. (index(unique_obs(iuobs),'DY') > 0)) then
- ! OSISAF Ice drift observations (d-2-offset -> d-offset)
- !print *, 'Ice drift observation type: ', unique_obs(iuobs)
- offset = unique_obs(iuobs)(3:3)
- ! Use offset (1,2,3,4 or 5) to open correct model drift file
- inquire(iolength=reclDRIFT) dX, dY
- open(10, file = 'model_ICEDRIFT_OSISAF'//offset//'.uf', access = 'direct',&
- status = 'old', recl = reclDRIFT, action = 'read')
- do iens = 1, nrens
- read(10, rec = iens, iostat = ios) dX, dY
- if (ios /= 0) then
- if (master) then
- print *,'ERROR: could not read ensemble ice drift for member ', iens
- end if
- call stop_mpi()
- end if
- do m = 1, nobs
- i = obs(m) % i_orig_grid
- j = obs(m) % j_orig_grid
- if (index(obs(m)%id,'DX') > 0) then
- S(m, iens) = dX(i, j)
- elseif (index(obs(m)%id,'DY') > 0) then
- S(m, iens) = dY(i, j)
- end if
- end do
- end do
- close(10)
- else
- if (master) then
- print *,'ERROR: unknown obs type ' // trim(unique_obs(iuobs))
- end if
- call stop_mpi()
- end if
- end do ! iuobs
- ! some generic QC - relax fitting if the model and obs are too far apart
- !
- call obs_QC(nrens, S)
- ! add calculated HA to to observations-<type>.nc files for each data type
- !
- do iuobs = 1, nuobs
- if (master) then
- nthisobs = 0
- do m = 1, nobs
- if (trim(unique_obs(iuobs)) == trim(obs(m) % id)) then
- nthisobs = nthisobs + 1
- thisobs(nthisobs) = m
- end if
- end do
- ! add forecast values to the observation-<TYPE>.nc files
- !
- call add_forecast(unique_obs(iuobs), S(thisobs(1 : nthisobs), :), obs(thisobs(1 : nthisobs)))
- ! append the superobed values (and modified observation error
- ! variances) to the file with pre-processed observations (SAL.nc,
- ! TEM.nc, GSAL.nc or GTEM.nc)
- !
- if (trim(unique_obs(iuobs)) == 'SAL' .or.&
- trim(unique_obs(iuobs)) == 'TEM' .or.&
- trim(unique_obs(iuobs)) == 'GSAL' .or.&
- trim(unique_obs(iuobs)) == 'GTEM') then
-
- call insitu_writeforecast(unique_obs(iuobs), nobs, nrens, S, obs)
- end if
- end if
- end do
- if (master) then
- print *, 'm_prep_4_EnKF: end calculating S = HA'
- end if
- x = sum(S, DIM = 2) / real(nrens) ! [ FM ] The mean forecast interpolated in the obs.space
- if (master) print*,'m_prep_4_EnKF: end calculating Hx'
- if (master) then
- print *, 'Hx range = ', minval(x), '-', maxval(x)
- print *, 'mean(Hx) = ', sum(x) / real(nobs)
- end if
- if (master) then
- print *, 'observation range = ', minval(obs % d), '-', maxval(obs % d)
- print *, 'mean(observation) = ', sum(obs % d) / nobs
- end if
- ! Compute HA = HE - mean(HE)
- !
- if (master) print*,'prep_4_EnKF(): calculating S = S - x'
- do j = 1, nrens
- S(:, j) = S(:, j) - x ! [ FM ] This is really where we switch from actual model data to anomalies
- enddo
- ! Compute innovation
- !
- d = obs % d - x ! [ FM ] This is exactly was is also done in add_forecast. This is the mean innovation.
- if (master) then
- print *, ' innovation range = ', minval(d), '-', maxval(d)
- if (minval(d) < -1000.0d0) then
- print *, 'm_prep_4_EnKF: error: innovation too small detected'
- call stop_mpi()
- end if
- if (maxval(d) > 1000.0d0) then
- print *, 'm_prep_4_EnKF: error: innovation too big detected'
- call stop_mpi()
- end if
- end if
- end subroutine prep_4_EnKF
- subroutine read_mean_ssh(mean_ssh, nx, ny)
- #if defined (QMPI)
- use qmpi
- #else
- use qmpi_fake
- #endif
- use m_parameters
- integer, intent(in) :: nx, ny
- real, intent(out):: mean_ssh(nx, ny)
- logical :: exists
- inquire(file = trim(MEANSSHFNAME), exist = exists)
- if (.not. exists) then
- if (master) then
- print *,'ERROR: read_mean_ssh(): file "', trim(MEANSSHFNAME), '" not found'
- end if
- stop
- end if
-
- open (10, file = trim(MEANSSHFNAME), status = 'unknown',form = 'unformatted', action = 'read')
- read (10) mean_ssh
- close (10)
- end subroutine read_mean_ssh
- ! This subroutine adds forecast observations (i.e Hx) to the NetCDF
- ! observation files for each data type.
- !
- subroutine add_forecast(obstag, S, obs)
- use mod_measurement
- use nfw_mod
- implicit none
-
- character(OBSTYPESTRLEN), intent(in) :: obstag
- real, dimension(:, :), intent(in) :: S
- type(measurement), dimension(:) :: obs
- character(STRLEN) :: fname
- logical :: exists
- integer :: ncid
- integer :: dids(2), dimlen
- logical :: addsobs
- integer :: for_id, inn_id, forvar_id, slon_id, slat_id,&
- sdepth_id, sipiv_id, sjpiv_id, sd_id, svar_id,&
- newvar_id
-
- real, allocatable, dimension(:) :: x, Svar, innovation
-
- integer :: m, p, o
- write(fname, '(a, a, a)') 'observations-', trim(obstag), '.nc'
- inquire(file = trim(fname), exist = exists)
- if (.not. exists) then
- print *, 'file "', trim(fname), 'not found, skip adding forecast'
- return
- else
- print *, 'dumping forecast to "', trim(fname), '"'
- end if
- p = size(S, DIM = 1)
- m = size(S, DIM = 2)
- allocate(x(p), Svar(p), innovation(p))
- x = sum(S, DIM = 2) / real(m); ! [ FM the mean of S=HA ]
- Svar = 0.0
- do o = 1, p
- Svar(o) = sum((S(o, :) - x(o))** 2) ! [ FM thus each row of Svar is the variance (see below) of the forecast]
- end do
- Svar = Svar / real(m - 1)
- innovation = obs % d - x ! [ FM ] the innovation for the mean forecast (or mean of the innovation forecasts)
-
- addsobs = .false.
- call nfw_open(fname, nf_write, ncid)
- call nfw_inq_dimid(fname, ncid, 'nobs', dids(1))
- call nfw_inq_dimlen(fname, ncid, dids(1), dimlen)
- call nfw_redef(fname, ncid)
- if (dimlen == p) then
- dids(2) = dids(1)
- elseif (.not. nfw_dim_exists(ncid, 'nsobs')) then
- addsobs = .true.
- call nfw_def_dim(fname, ncid, 'nsobs', p, dids(2))
- call nfw_def_var(fname, ncid, 'slon', nf_float, 1, dids(2), slon_id)
- call nfw_def_var(fname, ncid, 'slat', nf_float, 1, dids(2), slat_id)
- call nfw_def_var(fname, ncid, 'sdepth', nf_float, 1, dids(2), sdepth_id)
- call nfw_def_var(fname, ncid, 'sipiv', nf_int, 1, dids(2), sipiv_id)
- call nfw_def_var(fname, ncid, 'sjpiv', nf_int, 1, dids(2), sjpiv_id)
- call nfw_def_var(fname, ncid, 'sd', nf_float, 1, dids(2), sd_id)
- call nfw_def_var(fname, ncid, 'svar', nf_float, 1, dids(2), svar_id)
- end if
- if (.not. nfw_var_exists(ncid, 'innovation')) then
- call nfw_def_var(fname, ncid, 'innovation', nf_double, 1, dids(2), inn_id)
- else
- call nfw_inq_varid(fname, ncid, 'innovation', inn_id)
- end if
- if (.not. nfw_var_exists(ncid, 'forecast')) then
- call nfw_def_var(fname, ncid, 'forecast', nf_double, 1, dids(2), for_id)
- else
- call nfw_inq_varid(fname, ncid, 'forecast', for_id)
- end if
- if (.not. nfw_var_exists(ncid, 'forecast_variance')) then
- call nfw_def_var(fname, ncid, 'forecast_variance', nf_double, 1, dids(2), forvar_id)
- else
- call nfw_inq_varid(fname, ncid, 'forecast_variance', forvar_id)
- end if
- if (.not. addsobs) then
- if (dimlen == p) then
- if (.not. nfw_var_exists(ncid, 'new_var')) then
- call nfw_def_var(fname, ncid, 'new_var', nf_double, 1, dids(2), newvar_id)
- else
- call nfw_inq_varid(fname, ncid, 'new_var', newvar_id)
- end if
- else
- if (.not. nfw_var_exists(ncid, 'new_svar')) then
- call nfw_inq_dimid(fname, ncid, 'nsobs', dids(2))
- call nfw_def_var(fname, ncid, 'new_svar', nf_double, 1, dids(2), newvar_id)
- else
- call nfw_inq_varid(fname, ncid, 'new_svar', newvar_id)
- end if
- end if
- end if
- call nfw_enddef(fname, ncid)
- call nfw_put_var_double(fname, ncid, forvar_id, Svar)
- call nfw_put_var_double(fname, ncid, for_id, x)
- call nfw_put_var_double(fname, ncid, inn_id, innovation)
- if (addsobs) then
- call nfw_put_var_double(fname, ncid, slon_id, obs % lon)
- call nfw_put_var_double(fname, ncid, slat_id, obs % lat)
- call nfw_put_var_double(fname, ncid, sdepth_id, obs % depth)
- call nfw_put_var_int(fname, ncid, sipiv_id, obs % ipiv)
- call nfw_put_var_int(fname, ncid, sjpiv_id, obs % jpiv)
- call nfw_put_var_double(fname, ncid, sd_id, obs % d)
- call nfw_put_var_double(fname, ncid, svar_id, obs % var)
- else
- call nfw_put_var_double(fname, ncid, newvar_id, obs % var)
- end if
- call nfw_close(fname, ncid)
- deallocate(x)
- deallocate(Svar)
- deallocate(innovation)
- end subroutine add_forecast
- end module m_prep_4_EnKF
|