!/===========================================================================/
! Copyright (c) 2007, The University of Massachusetts Dartmouth
! Produced at the School of Marine Science & Technology
! Marine Ecosystem Dynamics Modeling group
! All rights reserved.
!
! FVCOM has been developed by the joint UMASSD-WHOI research team. For
! details of authorship and attribution of credit please see the FVCOM
! technical manual or contact the MEDM group.
!
!
! This file is part of FVCOM. For details, see http://fvcom.smast.umassd.edu
! The full copyright notice is contained in the file COPYRIGHT located in the
! root directory of the FVCOM code. This original header must be maintained
! in all distributed versions.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
! THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED.
!
!/---------------------------------------------------------------------------/
! CVS VERSION INFORMATION
! $Id$
! $Name$
! $Revision$
!/===========================================================================/
MODULE ENKFVAL
# if defined (ENKF)
USE mod_ncdio, only : update_iodata
USE CONTROL
IMPLICIT NONE
CHARACTER(LEN=2) ENKF_NUM
CHARACTER(LEN=80) ENKF_INIT !!INITIAL PERTUBATION FIELD OPTION(usr/default)
CHARACTER(LEN=80) ENKF_TYPE !!ENSEMBLE UPDATE OPTION(default/square root)
INTEGER ENKF_NENS !!GLOBAL NUMBER OF ENSEMBLES
INTEGER DELTA_ASS !!ASSIMILATON TIME INTERVAL IN SECONDS
! INTEGER ENKF_INT !!ASSIMILATION TIME INTERVAL/FILE OUTPUT INTERVAL (>=1)
INTEGER ENKF_NOBSMAX !!MAXIMUM NUMBER OF THE OBSERVATION STATIONS
INTEGER ENKF_START !!ASSIMILATION START TIME
INTEGER ENKF_END !!ASSIMILATION END TIME
REAL(SP) :: ENKF_CINF !!MAX DISTANCE OF CORRELATIN
REAL(SP) :: OBSERR_EL !!EL OBSERVATION ERROR SPECIFIED
REAL(SP) :: OBSERR_UV !!UV OBSERVATION ERROR SPECIFIED
REAL(SP) :: OBSERR_T !!TEMPERATURE OBSERVATION ERROR SPECIFIED
REAL(SP) :: OBSERR_S !!SALINITY OBSERVATION ERROR SPECIFIED
LOGICAL :: EL_ASSIM !!OPTION FOR CHOSING ELEVATION AS ASSIMILATION VARIABLES
LOGICAL :: UV_ASSIM !!OPTION FOR CHOSING CURRENT AS ASSIMILATION VARIABLES
LOGICAL :: T_ASSIM !!OPTION FOR CHOSING TEMPERATURE AS ASSIMILATION VARIABLES
LOGICAL :: S_ASSIM !!OPTION FOR CHOSING SALINITY AS ASSIMILATION VARIABLES
LOGICAL :: ENKF_LOCALIZED !! only do enkf for certain region
LOGICAL :: EL_OBS !!OPTION FOR ELEVATION OBSERVATION DATA
LOGICAL :: UV_OBS !!OPTION FOR CURRENT OBSERVATION DATA
LOGICAL :: T_OBS !!OPTION FOR TEMPERATURE OBSERVATION DATA
LOGICAL :: S_OBS !!OPTION FOR SALINITY OBSERVATION DATA
INTEGER :: ENKF_METHOD !! DIFFERENENT ENKF_METHODS
INTEGER :: MODE !! EVENSEN_KF_OPTIONS
INTEGER :: INOOB !! 72 RESERVE FOR I/O INPUT OF OBSERVATION FILE
INTEGER :: INOKF !! 73 FILE I/O PIPE NUMBER
INTEGER :: IOBCKF !! 76 I/O PIPE ONLY FOR B.C. TREATMENT IN ENKF
INTEGER :: NLOC=0 !!Number of observation locations
INTEGER :: IENS !!ensemble number index
INTEGER :: NCYC !!Number of cycles for assimilation
INTEGER :: ICYC = 0 !!cycle number index
CHARACTER(LEN=4) FCYC
REAL(DP),PARAMETER :: ZEROD = 0.0_DP
INTEGER :: KF_TIMES =0
INTEGER :: ARCHIVE_TIMES=0
TYPE(TIME) :: enkf_out_Time
# endif
END MODULE ENKFVAL
!=============================================
!MODULE MOD_STARTUP_MIMIC
!# if defined (ENKF)
! USE MOD_UTILS
! USE MOD_NCTOOLS
! USE MOD_INPUT
! USE ALL_VARS
! USE EQS_OF_STATE
! USE MOD_WD
! USE SINTER
!
!
!# if defined (ICE)
! USE MOD_ICE
! USE MOD_ICE2D
!# endif
!
!# if defined (WATER_QUALITY)
! USE MOD_WQM
!# endif
!
!# if defined (BioGen)
! USE MOD_BIO_3D
!# endif
!
!# if defined (NH)
! USE NON_HYDRO, ONLY: W4ZT, NHQDRX, NHQDRY, NHQDRZ, NHQ2DX, NHQ2DY
!# endif
!
! IMPLICIT NONE
!
! PRIVATE
!
! PUBLIC :: READ_SSH1
! PUBLIC :: READ_UV1
! PUBLIC :: READ_TURB1
! PUBLIC :: READ_TS1
! PUBLIC :: READ_WETDRY1
!
!
!CONTAINS
!!==============================================================================!
! SUBROUTINE READ_SSH1(NC_START)
!# if defined (SEDIMENT)
! USE MOD_SED, only : morpho_model,sed_hot_start
!# endif
! IMPLICIT NONE
! TYPE(NCFILE),POINTER :: NC_START
! TYPE(NCVAR), POINTER :: VAR
! TYPE(NCDIM), POINTER :: DIM
! LOGICAL :: FOUND
! INTEGER :: STKCNT
!
! IF(DBG_SET(DBG_SBR)) write(ipt,*) "Start: READ_SSH"
!
! STKCNT = NC_START%FTIME%PREV_STKCNT
!
! ! LOAD EL
! VAR => FIND_VAR(NC_START,'zeta',FOUND)
! IF(.not. FOUND) CALL FATAL_ERROR("COULD NOT FIND VARIABLE 'zeta'&
! & IN THE HOTSTART FILE OBJECT")
! CALL NC_CONNECT_AVAR(VAR, EL)
! CALL NC_READ_VAR(VAR,STKCNT)
!
! ! LOAD ET
! VAR => FIND_VAR(NC_START,'et',FOUND)
! IF(.not. FOUND) CALL FATAL_ERROR("COULD NOT FIND VARIABLE 'et'&
! & IN THE HOTSTART FILE OBJECT")
! CALL NC_CONNECT_AVAR(VAR, ET)
! CALL NC_READ_VAR(VAR,STKCNT)
!
! !----------------------------------------------------------------
! ! Read the most recent bathymetry if Morphodynamics is Active
! !----------------------------------------------------------------
!# if defined(SEDIMENT)
! IF(MORPHO_MODEL .and. SED_HOT_START)THEN
! VAR => FIND_VAR(NC_START,'h',FOUND)
! IF(.not. FOUND) CALL FATAL_ERROR("COULD NOT FIND VARIABLE 'h'&
! & IN THE HOTSTART FILE OBJECT")
! CALL NC_CONNECT_AVAR(VAR, h)
! CALL NC_READ_VAR(VAR,STKCNT)
! ENDIF
!# endif
!
!
! !----------------------------------------------------------------
! ! Given SSH and Bathy, Update the Bathymetry
! !----------------------------------------------------------------
! D = H + EL
! DT = H + ET
!
!
! CALL N2E2D(H,H1)
! CALL N2E2D(EL,EL1)
! CALL N2E2D(D,D1)
! CALL N2E2D(DT,DT1)
!
! IF(DBG_SET(DBG_SBR)) write(ipt,*) "Start: READ_SSH"
!
! END SUBROUTINE READ_SSH1
!!==============================================================================!
!!==============================================================================!
! SUBROUTINE READ_WETDRY1(NC_START)
! USE MOD_WD
! IMPLICIT NONE
! TYPE(NCFILE),POINTER :: NC_START
! TYPE(NCVAR), POINTER :: VAR
! TYPE(NCDIM), POINTER :: DIM
! LOGICAL :: FOUND
! INTEGER :: STKCNT
!
! IF(DBG_SET(DBG_SBR)) write(ipt,*) "Start: READ_WETDRY"
!
! STKCNT = NC_START%FTIME%PREV_STKCNT
!
! ! LOAD ISWETN
! VAR => FIND_VAR(NC_START,'wet_nodes',FOUND)
! IF(.not. FOUND) CALL FATAL_ERROR("COULD NOT FIND VARIABLE 'wet_nodes'&
! & IN THE HOTSTART FILE OBJECT")
! CALL NC_CONNECT_AVAR(VAR, ISWETN)
! CALL NC_READ_VAR(VAR,STKCNT)
!
! ! LOAD ISWETC
! VAR => FIND_VAR(NC_START,'wet_cells',FOUND)
! IF(.not. FOUND) CALL FATAL_ERROR("COULD NOT FIND VARIABLE 'wet_cells'&
! & IN THE HOTSTART FILE OBJECT")
! CALL NC_CONNECT_AVAR(VAR, ISWETC)
! CALL NC_READ_VAR(VAR,STKCNT)
!
!
! ! LOAD ISWETNT
! VAR => FIND_VAR(NC_START,'wet_nodes_prev_int',FOUND)
! IF(.not. FOUND) CALL FATAL_ERROR("COULD NOT FIND VARIABLE 'wet_nodes_prev_int'&
! & IN THE HOTSTART FILE OBJECT")
! CALL NC_CONNECT_AVAR(VAR, ISWETNT)
! CALL NC_READ_VAR(VAR,STKCNT)
!
! ! LOAD ISWETCT
! VAR => FIND_VAR(NC_START,'wet_cells_prev_int',FOUND)
! IF(.not. FOUND) CALL FATAL_ERROR("COULD NOT FIND VARIABLE 'wet_cells_prev_int'&
! & IN THE HOTSTART FILE OBJECT")
! CALL NC_CONNECT_AVAR(VAR, ISWETCT)
! CALL NC_READ_VAR(VAR,STKCNT)
!
! ! LOAD ISWETCE
! VAR => FIND_VAR(NC_START,'wet_cells_prev_ext',FOUND)
! IF(.not. FOUND) CALL FATAL_ERROR("COULD NOT FIND VARIABLE 'wet_cells_prev_ext'&
! & IN THE HOTSTART FILE OBJECT")
! CALL NC_CONNECT_AVAR(VAR, ISWETC)
! CALL NC_READ_VAR(VAR,STKCNT)
!
! IF(DBG_SET(DBG_SBR)) write(ipt,*) "End: READ_WETDRY"
!
! END SUBROUTINE READ_WETDRY1
!!==============================================================================!
! SUBROUTINE READ_TS1(NC_START)
! IMPLICIT NONE
! TYPE(NCFILE),POINTER :: NC_START
! TYPE(NCVAR), POINTER :: VAR
! TYPE(NCDIM), POINTER :: DIM
! LOGICAL :: FOUND
! INTEGER :: STKCNT, K
! REAL(SP), DIMENSION(0:MT,KB) :: PRESSURE
!
! IF(DBG_SET(DBG_SBR)) write(ipt,*) "Start: READ_TS"
!
! STKCNT = NC_START%FTIME%PREV_STKCNT
!
!
! ! LOAD TEMPERATURE
! VAR => FIND_VAR(NC_START,'temp',FOUND)
! IF(.not. FOUND) CALL FATAL_ERROR("COULD NOT FIND VARIABLE 'temp'&
! & IN THE HOTSTART FILE OBJECT")
! CALL NC_CONNECT_AVAR(VAR, T1)
! CALL NC_READ_VAR(VAR,STKCNT)
!
!
! ! LOAD MEAN INITIAL TEMPERATURE
! VAR => FIND_VAR(NC_START,'tmean1',FOUND)
! IF(.not. FOUND) CALL FATAL_ERROR("COULD NOT FIND VARIABLE 'tmean1'&
! & IN THE HOTSTART FILE OBJECT")
! CALL NC_CONNECT_AVAR(VAR, tmean1)
! CALL NC_READ_VAR(VAR)
!
! ! LOAD SALINITY
! VAR => FIND_VAR(NC_START,'salinity',FOUND)
! IF(.not. FOUND) CALL FATAL_ERROR("COULD NOT FIND VARIABLE 'saltinity'&
! & IN THE HOTSTART FILE OBJECT")
! CALL NC_CONNECT_AVAR(VAR, S1)
! CALL NC_READ_VAR(VAR,STKCNT)
!
! ! LOAD MEAN INITIAL SALINITY
! VAR => FIND_VAR(NC_START,'smean1',FOUND)
! IF(.not. FOUND) CALL FATAL_ERROR("COULD NOT FIND VARIABLE 'smean1'&
! & IN THE HOTSTART FILE OBJECT")
! CALL NC_CONNECT_AVAR(VAR, smean1)
! CALL NC_READ_VAR(VAR)
!
!
! ! LOAD DENSITY
! VAR => FIND_VAR(NC_START,'rho1',FOUND)
! IF(.not. FOUND) CALL FATAL_ERROR("COULD NOT FIND VARIABLE 'rho1'&
! & IN THE HOTSTART FILE OBJECT")
! CALL NC_CONNECT_AVAR(VAR, RHO1)
! CALL NC_READ_VAR(VAR,STKCNT)
!
! ! LOAD MEAN DENSITY
! VAR => FIND_VAR(NC_START,'rmean1',FOUND)
! IF(.not. FOUND) CALL FATAL_ERROR("COULD NOT FIND VARIABLE 'rmean1'&
! & IN THE HOTSTART FILE OBJECT")
! CALL NC_CONNECT_AVAR(VAR, rmean1)
! CALL NC_READ_VAR(VAR)
!
! ! AVERAGE FROM CELLS TO FACE CENTERS
!
!
!!JQI SELECT CASE(SEA_WATER_DENSITY_FUNCTION)
!!JQI CASE(SW_DENS1)
!
!!JQI ! SET MEAN DENSITY
!!JQI DO K=1,KBM1
!!JQI PRESSURE(:,K) = -GRAV_N*1.025_SP*(ZZ(:,K)*D(:))*0.1_SP
!!JQI END DO
!!JQI CALL FOFONOFF_MILLARD(SMEAN1,TMEAN1,PRESSURE,0.0_SP,RMEAN1)
!!JQI RMEAN1(0,:)=0.0_SP
!!JQI RMEAN1(:,KB)=0.0_SP
!
!!JQI ! SET REAL DENSITY
!!JQI CALL DENS1 ! GENERIC CALL TO FOFONOFF_MILLARD FOR S1,T1...
!
!!JQI CASE(SW_DENS2)
!!JQI ! SET MEAN DENSITY
!!JQI CALL DENS2G(SMEAN1,TMEAN1,RMEAN1)
!!JQI RMEAN1(0,:)=0.0_SP
!!JQI RMEAN1(:,KB)=0.0_SP
!
!!JQI ! SET REAL DENSITY
!!JQI CALL DENS2 ! GENERIC CALL TO DENS2G FOR S1,T1...
!
!!JQI CASE(SW_DENS3)
!
!!JQI ! SET MEAN DENSITY
!!JQI DO K=1,KBM1
!!JQI PRESSURE(:,K) = -GRAV_N*1.025_SP*(ZZ(:,K)*D(:))*0.1_SP
!!JQI END DO
!!JQI CALL JACKET_MCDOUGALL(SMEAN1,TMEAN1,PRESSURE,RMEAN1)
!!JQI RMEAN1(0,:)=0.0_SP
!!JQI RMEAN1(:,KB)=0.0_SP
!
!!JQI ! SET REAL DENSITY
!!JQI CALL DENS3 ! GENERIC CALL TO JACKET_MCDOUGALL FOR S1,T1.
!
!!JQI CASE DEFAULT
!!JQI CALL FATAL_ERROR("INVALID DENSITY FUNCTION SELECTED:",&
!!JQI & " "//TRIM(SEA_WATER_DENSITY_FUNCTION) )
!!JQI END SELECT
!
! CALL N2E3D(T1,T)
! CALL N2E3D(S1,S)
! CALL N2E3D(Tmean1,Tmean)
! CALL N2E3D(Smean1,Smean)
! CALL N2E3D(Rmean1,Rmean)
!
!
! IF(DBG_SET(DBG_SBR)) write(ipt,*) "End: READ_TS"
!
! END SUBROUTINE READ_TS1
!!==============================================================================!
!
! SUBROUTINE READ_UV1(NC_START)
! IMPLICIT NONE
! TYPE(NCFILE),POINTER :: NC_START
! TYPE(NCVAR), POINTER :: VAR
! TYPE(NCDIM), POINTER :: DIM
! LOGICAL :: FOUND
! INTEGER :: STKCNT
!
! IF(DBG_SET(DBG_SBR)) write(ipt,*) "START: READ_UV"
!
!
! STKCNT = NC_START%FTIME%PREV_STKCNT
!
! VAR => FIND_VAR(NC_START,'u',FOUND)
! IF(.not. FOUND) CALL FATAL_ERROR("COULD NOT FIND VARIABLE 'u'&
! & IN THE HOTSTART FILE OBJECT")
! CALL NC_CONNECT_AVAR(VAR, U)
! CALL NC_READ_VAR(VAR,STKCNT)
!
!
! VAR => FIND_VAR(NC_START,'v',FOUND)
! IF(.not. FOUND) CALL FATAL_ERROR("COULD NOT FIND VARIABLE 'v'&
! & IN THE HOTSTART FILE OBJECT")
! CALL NC_CONNECT_AVAR(VAR, V)
! CALL NC_READ_VAR(VAR,STKCNT)
!
! VAR => FIND_VAR(NC_START,'omega',FOUND)
! IF(FOUND) THEN
! CALL NC_CONNECT_AVAR(VAR, WTS)
! CALL NC_READ_VAR(VAR,STKCNT)
!
! CALL N2E3D(WTS,W)
! ELSE
! VAR => FIND_VAR(NC_START,'w',FOUND)
! IF(.not. FOUND) CALL FATAL_ERROR("COULD NOT FIND VARIABLE 'w' &
! & or 'omega' IN THE HOTSTART FILE OBJECT")
! CALL NC_CONNECT_AVAR(VAR, W)
! CALL NC_READ_VAR(VAR,STKCNT)
! END IF
!
! VAR => FIND_VAR(NC_START,'ua',FOUND)
! IF(.not. FOUND) CALL FATAL_ERROR("COULD NOT FIND VARIABLE 'ua'&
! & IN THE HOTSTART FILE OBJECT")
! CALL NC_CONNECT_AVAR(VAR, UA)
! CALL NC_READ_VAR(VAR,STKCNT)
!
! VAR => FIND_VAR(NC_START,'va',FOUND)
! IF(.not. FOUND) CALL FATAL_ERROR("COULD NOT FIND VARIABLE 'va'&
! & IN THE HOTSTART FILE OBJECT")
! CALL NC_CONNECT_AVAR(VAR, VA)
! CALL NC_READ_VAR(VAR,STKCNT)
!
! IF(DBG_SET(DBG_SBR)) write(ipt,*) "END: READ_UV"
!
! END SUBROUTINE READ_UV1
!!==============================================================================!
! SUBROUTINE READ_TURB1(NC_START)
! IMPLICIT NONE
! TYPE(NCFILE),POINTER :: NC_START
! TYPE(NCVAR), POINTER :: VAR
! TYPE(NCDIM), POINTER :: DIM
! LOGICAL :: FOUND
! INTEGER :: STKCNT
!
! IF(DBG_SET(DBG_SBR)) write(ipt,*) "START: READ_TURB"
!
! STKCNT = NC_START%FTIME%PREV_STKCNT
!
!# if defined (GOTM)
!
! VAR => FIND_VAR(NC_START,'tke',FOUND)
! IF(.not. FOUND) CALL FATAL_ERROR("COULD NOT FIND VARIABLE 'tke'&
! & IN THE HOTSTART FILE OBJECT")
! CALL NC_CONNECT_AVAR(VAR, TKE)
! CALL NC_READ_VAR(VAR,STKCNT)
!
! VAR => FIND_VAR(NC_START,'teps',FOUND)
! IF(.not. FOUND) CALL FATAL_ERROR("COULD NOT FIND VARIABLE 'teps'&
! & IN THE HOTSTART FILE OBJECT")
! CALL NC_CONNECT_AVAR(VAR, TEPS)
! CALL NC_READ_VAR(VAR,STKCNT)
!
! L = .001
! L(1:MT,2:KBM1) = (.5544**3)*TKE(1:MT,2:KBM1)**1.5/TEPS(1:MT,2:KBM1)
!
!# else
!
! VAR => FIND_VAR(NC_START,'q2',FOUND)
! IF(.not. FOUND) CALL FATAL_ERROR("COULD NOT FIND VARIABLE 'q2'&
! & IN THE HOTSTART FILE OBJECT")
! CALL NC_CONNECT_AVAR(VAR, Q2)
! CALL NC_READ_VAR(VAR,STKCNT)
!
! VAR => FIND_VAR(NC_START,'q2l',FOUND)
! IF(.not. FOUND) CALL FATAL_ERROR("COULD NOT FIND VARIABLE 'q2l'&
! & IN THE HOTSTART FILE OBJECT")
! CALL NC_CONNECT_AVAR(VAR, Q2L)
! CALL NC_READ_VAR(VAR,STKCNT)
!
! VAR => FIND_VAR(NC_START,'l',FOUND)
! IF(.not. FOUND) CALL FATAL_ERROR("COULD NOT FIND VARIABLE 'l'&
! & IN THE HOTSTART FILE OBJECT")
! CALL NC_CONNECT_AVAR(VAR, L)
! CALL NC_READ_VAR(VAR,STKCNT)
!
!# endif
!
! VAR => FIND_VAR(NC_START,'km',FOUND)
! IF(.not. FOUND) CALL FATAL_ERROR("COULD NOT FIND VARIABLE 'km'&
! & IN THE HOTSTART FILE OBJECT")
! CALL NC_CONNECT_AVAR(VAR, km)
! CALL NC_READ_VAR(VAR,STKCNT)
!
! VAR => FIND_VAR(NC_START,'kq',FOUND)
! IF(.not. FOUND) CALL FATAL_ERROR("COULD NOT FIND VARIABLE 'kq'&
! & IN THE HOTSTART FILE OBJECT")
! CALL NC_CONNECT_AVAR(VAR, KQ)
! CALL NC_READ_VAR(VAR,STKCNT)
!
! VAR => FIND_VAR(NC_START,'kh',FOUND)
! IF(.not. FOUND) CALL FATAL_ERROR("COULD NOT FIND VARIABLE 'kh'&
! & IN THE HOTSTART FILE OBJECT")
! CALL NC_CONNECT_AVAR(VAR, KH)
! CALL NC_READ_VAR(VAR,STKCNT)
!
!
! CALL N2E3D(KM,KM1)
!
! IF(DBG_SET(DBG_SBR)) write(ipt,*) "END: READ_TURB"
!
! END SUBROUTINE READ_TURB1
!!==============================================================================|
!# endif
!END MODULE MOD_STARTUP_MIMIC
MODULE MOD_ENKF
USE ENKFVAL
USE mod_ncdio, only : update_iodata
USE MOD_INPUT, only : NC_START
USE CONTROL
USE MOD_UTILS
USE MOD_NCTOOLS
IMPLICIT NONE
SAVE
LOGICAL :: ENKF_ON
# if defined (ENKF)
INTEGER :: ENKF_memberCONTR !, ICYC
CHARACTER(LEN=80) :: ENKF_START_DATE
CHARACTER(LEN=80) :: ENKF_END_DATE
TYPE(TIME) :: ENKF_START_TIME
TYPE(TIME) :: ENKF_END_Time
TYPE(TIME) :: ENKF_cyc
CHARACTER(LEN=80) :: ENKF_ASSIM_INTERVAL
TYPE(TIME) :: ENKF_INTERVAL
! TYPE(NCFILE), POINTER :: NCF_RRKFDATA
TYPE(NCFILE), POINTER :: NCF_ENKFRE
TYPE(NCFILE), POINTER :: NCF_ENKFfct
TYPE(NCFILE), POINTER :: enkf_ncfout,enkf_mean_ncfout
LOGICAL :: LOCAL_DISK
CHARACTER(LEN=80), parameter :: Local_file ="/scratch/enkf_restart.nc"
CHARACTER(LEN=80), parameter :: group_file ="rrktemp/enkf_restart.nc"
CHARACTER(LEN=80), parameter :: Local_file1 ="/scratch/enkf_forecast.nc"
CHARACTER(LEN=80), parameter :: group_file1 ="rrktemp/enkf_forecast.nc"
INTEGER(ITIME) :: EKINT_COUNT, EKINT_START, EKINT_END
! TYPE(NCFILE), POINTER :: enkf_ncfout_enkf(20)
! INTERFACE READRESTART
! MODULE PROCEDURE READRESTART_TIME
! MODULE PROCEDURE READRESTART_STATE
! END INTERFACE
! PRIVATE READRESTART_TIME
! PRIVATE READRESTART_STATE
NAMELIST /NML_ENKF/ &
& ENKF_ON, &
& ENKF_START_DATE, &
& ENKF_END_DATE, &
& ENKF_ASSIM_INTERVAL,&
& ENKF_NOBSMAX, &
& ENKF_NENS, &
& ENKF_CINF, &
& EKINT_START, &
& EL_ASSIM, &
& EL_OBS, &
& UV_ASSIM, &
& UV_OBS, &
& T_ASSIM, &
& T_OBS, &
& S_ASSIM, &
& S_OBS, &
& ENKF_LOCALIZED, &
& ENKF_METHOD, &
& MODE, &
& OBSERR_EL, &
& OBSERR_UV, &
& OBSERR_T, &
& OBSERR_S, &
& LOCAL_DISK
! WRITE(IPT,*) '! # SPECIFICY INITIAL PERTUBATION FIELD :',ENKF_INIT
! WRITE(IPT,*) '! # GLOBAL NUMBER OF ENSEMBLES :',ENKF_NENS
! WRITE(IPT,*) '! # ASSIMILATON TIME INTERVAL IN SECONDS :',DELTA_ASS
! WRITE(IPT,*) '! # ASSIMILATION TIME INTERVAL/FILE OUTPUT INTERVAL :',ENKF_INT
! WRITE(IPT,*) '! # MAXIMUM NUMBER OF THE OBSERVATION STATIONS :',ENKF_NOBSMAX
! WRITE(IPT,*) '! # ASSIMILATION START TIME :',ENKF_START
! WRITE(IPT,*) '! # ASSIMILATION END TIME :',ENKF_END
! WRITE(IPT,*) '! # TIDAL AMPLITUDE ERROR RANGE SPECIFIED :',(BC_AMP_ERR(I),I=1,6)
! WRITE(IPT,*) '! # TIDAL PHASE ERROR RANGE SPECIFIED :',(BC_PHA_ERR(I),I=1,6)
! WRITE(IPT,*) '! # MAX DISTANCE OF CORRELATIN :',ENKF_CINF
! WRITE(IPT,*) '! # ELEVATION AS ASSIMILATION VARIABLES :',EL_ASSIM
! WRITE(IPT,*) '! # ELEVATION AS OBSERVATION DATA OPTION :',EL_OBS
! WRITE(IPT,*) '! # CURRENTS AS ASSIMILATION VARIABLES :',UV_ASSIM
! WRITE(IPT,*) '! # CURRENTS OBSERVATION DATA OPTION :',UV_OBS
! WRITE(IPT,*) '! # TEMPERATURE AS ASSIMILATION VARIABLES :',T_ASSIM
! WRITE(IPT,*) '! # TEMPERATURE OBSERVATION DATA OPTION :',T_OBS
! WRITE(IPT,*) '! # SALINITY AS ASSIMILATION VARIABLES :',S_ASSIM
! WRITE(IPT,*) '! # SALINITY OBSERVATION DATA OPTION :',S_OBS
! IF(EL_OBS) WRITE(IPT,*) '! # EL OBSERVATION ERROR SPECIFIED :',OBSERR_EL
! IF(UV_OBS) WRITE(IPT,*) '! # UV OBSERVATION ERROR SPECIFIED :',OBSERR_UV
! IF(T_OBS) WRITE(IPT,*) '! # T OBSERVATION ERROR SPECIFIED :',OBSERR_T
! IF(S_OBS) WRITE(IPT,*) '! # S OBSERVATION ERROR SPECIFIED :',OBSERR_S
! WRITE(IPT,*) '! !'
CONTAINS
SUBROUTINE NAME_LIST_INITIALIZE_ENKF
IMPLICIT NONE
ENKF_START_DATE = "RRK ASSIMILATION START AND END TIME"
ENKF_END_DATE = "For an idealized case specify 'seconds=(flt)','days=(flt)', or 'cycles=(int)'"
ENKF_ASSIM_INTERVAL= "A length of time: 'seconds= ','days= ', or 'cycles= '"
ENKF_NOBSMAX = -1
ENKF_NENS = -1
ENKF_CINF = 1E20
EKINT_START=0
ENKF_ON = .FALSE.
EL_ASSIM = .FALSE.
EL_OBS = .FALSE.
UV_ASSIM = .FALSE.
UV_OBS = .FALSE.
T_ASSIM = .FALSE.
T_OBS = .FALSE.
S_ASSIM = .FALSE.
S_OBS = .FALSE.
ENKF_LOCALIZED = .FALSE.
ENKF_METHOD = 1
MODE = 11
LOCAL_DISK = .FALSE.
OBSERR_EL = 0.001
OBSERR_UV = 0.001
OBSERR_T = 0.01
OBSERR_S = 0.01
END SUBROUTINE NAME_LIST_INITIALIZE_ENKF
SUBROUTINE NAME_LIST_PRINT_ENKF
IMPLICIT NONE
write(UNIT=IPT,NML=NML_ENKF)
RETURN
END SUBROUTINE NAME_LIST_PRINT_ENKF
SUBROUTINE NAME_LIST_READ_ENKF
IMPLICIT NONE
CHARACTER(LEN=120) :: FNAME
INTEGER :: IOS
FNAME = TRIM(INPUT_DIR)//"/"//trim(casename)//"_assim_enkf.nml"
CALL FOPEN(KFUNIT,trim(FNAME),'cfr')
! Read RRKF parameters
READ(UNIT=KFUNIT, NML=NML_ENKF,IOSTAT=ios)
if(ios .NE. 0 ) then
if(DBG_SET(dbg_log)) write(UNIT=IPT,NML=NML_ENKF)
Call Fatal_Error("Can Not Read NameList NML_ENKF from file: "//trim(FNAME))
end if
REWIND(KFUNIT)
if(DBG_SET(dbg_scl)) &
& write(IPT,*) "Read_Name_List: NML_ENKF"
if(DBG_SET(dbg_scl)) &
& write(UNIT=IPT,NML=NML_ENKF)
CLOSE(NMLUNIT)
! write(*,NML=NML_ENKF) ! should be marked
if (ENKF_ASSIM_INTERVAL .ne. NC_OUT_INTERVAL) then
print *, "NC_OUT_INTERVAL must be same as ENKF_ASSIM_INTERVAL"
print *, 'you can set the ENKF_ASSIM_INTERVAL smaller than your observation interal to fit your nc_output_interval even though you dont have observation data. '
print *, 'stop'
call pstop
end if
print *, 'finish reading namelist enkf'
RETURN
END SUBROUTINE NAME_LIST_READ_ENKF
SUBROUTINE ENKF_SET_TIME
USE MOD_SET_TIME
IMPLICIT NONE
CHARACTER(LEN=4) :: FLAG,BFLAG,EFLAG
INTEGER(ITIME):: begins, ends,steps
integer :: status
! GET THE RRK ASSIMILATION INTERVAL PERIOD
CALL IDEAL_TIME_STRING2TIME(ENKF_ASSIM_INTERVAL,FLAG,ENKF_INTERVAL,steps)
IF (FLAG == 'time') THEN ! IF START AND END TIME WERE SPECIFIED
IF(MOD(ENKF_INTERVAL,IMDTI) /= ZeroTime) THEN
CALL FATAL_ERROR("ENKF_ASSIM_INTERVAL must be an integer number of internal time steps!")
END IF
steps = seconds(ENKF_INTERVAL)/seconds(IMDTI)
ELSEIF(FLAG == 'step') THEN ! IF SPECIFIED AS A NUMBER OF STEPS
ENKF_INTERVAL = steps * IMDTI
END IF
! EKInt_start=0
EKint_end = steps
SELECT CASE(USE_REAL_WORLD_TIME)
CASE(.TRUE.)
! REF_START_TIME = READ_DATETIME(REF_START_DATE,DATE_FORMAT,TIMEZONE,status)
! if (.not. status) &
! & Call Fatal_Error("Could not read the date string REF_START_DATE: "//trim(REF_START_DATE))
! ! GET THE END TIME
! REF_END_Time = READ_DATETIME(REF_END_DATE,DATE_FORMAT,TIMEZONE,status)
! if (.not. status) &
! & Call Fatal_Error("Could not read the date string REF_END_DATE: "&
! &//trim(REF_END_DATE))
! ! SANITY CHECK
! if(REF_Start_Time .GE. REF_End_Time) &
! & Call Fatal_Error("RRKF REF_Start_Date exceeds or equal to REF_End_Date")
ENKF_START_TIME = READ_DATETIME(ENKF_START_DATE,DATE_FORMAT,TIMEZONE,status)
!!$ if (.not. status) &
if (status == 0) &
& Call Fatal_Error("Could not read the date string ENKF_START_DATE: "//trim(ENKF_START_DATE))
! GET THE END TIME
ENKF_END_Time = READ_DATETIME(ENKF_END_DATE,DATE_FORMAT,TIMEZONE,status)
!!$ if (.not. status) &
if (status == 0) &
& Call Fatal_Error("Could not read the date string ENKF_END_DATE: "&
&//trim(ENKF_END_DATE))
! SANITY CHECK
if(ENKF_START_TIME .GE. ENKF_END_Time) &
& Call Fatal_Error("ENKF_START_DATE exceeds or equal to ENKF_END_DATE")
CASE (.FALSE.)
! print *, 'REF_START_DATE',REF_START_DATE
! GET THE REFERENCE START AND END INFORMATION
! CALL IDEAL_TIME_STRING2TIME(REF_START_DATE,BFLAG,REF_START_TIME,begins)
! CALL IDEAL_TIME_STRING2TIME(REF_END_DATE,EFLAG,REF_End_TIME,ends)
! SANITY CHECK
! IF (BFLAG /= EFLAG) CALL FATAL_ERROR&
! ('IDEALIZED MODEL TIME SPECIFICATION IS INCORRENT',&
! &'RRK REF BEGIN AND END CAN BE IN EITHER CYCLES OR TIME BUT NOT MIXED',&
! & trim(REF_start_date),trim(ref_end_date) )
! IF (EFLAG == 'time') THEN ! IF START AND END TIME WERE SPECIFIED
! DO NOTHING
! ELSE IF(EFLAG == 'step') THEN ! IF START AND END IINT WERE SPECIFIED
! REF_START_TIME = StartTime + IMDTI*(begins - ISTART)
! REF_END_TIME = StartTime + IMDTI*(ends - ISTART)
! ELSE
! print *, 'EFLAG:' , eflag
! CALL FATAL_ERROR&
! &('IDEAL_TIME_STRING2TIME returned invalid flag for rrk ref date and time')
! END IF
! GET THE RRK START AND END INFORMATION
CALL IDEAL_TIME_STRING2TIME(ENKF_START_DATE,BFLAG,ENKF_START_TIME,begins)
CALL IDEAL_TIME_STRING2TIME(ENKF_END_DATE,EFLAG,ENKF_END_Time,ends)
! SANITY CHECK
IF (BFLAG /= EFLAG) CALL FATAL_ERROR&
('IDEALIZED MODEL TIME SPECIFICATION IS INCORRENT',&
&'ENKF BEGIN AND END CAN BE IN EITHER CYCLES OR TIME BUT NOT MIXED',&
& trim(enkf_start_date),trim(enkf_end_date) )
IF (EFLAG == 'time') THEN ! IF START AND END TIME WERE SPECIFIED
! DO NOTHING
ELSE IF(EFLAG == 'step') THEN ! IF START AND END IINT WERE SPECIFIED
! ENKF_START_TIME = StartTime + IMDTI*(begins - ISTART)
ENKF_START_TIME = StartTime + IMDTI*(begins - ISTART+1)
! ENKF_END_Time = StartTime + IMDTI*(ends - ISTART)
ENKF_END_Time = StartTime + IMDTI*(ends - ISTART+1)
ELSE
CALL FATAL_ERROR&
&('IDEAL_TIME_STRING2TIME returned invalid flag for rrk ref date and time')
END IF
END SELECT
! IF(MOD(REF_START_TIME-StartTime,IMDTI) /= ZeroTime) THEN
! CALL FATAL_ERROR("REF_START_TIME must be an integer number of internal time steps from the start time!")
! END IF
! IF(MOD(REF_END_TIME-StartTime,IMDTI) /= ZeroTime) THEN
! CALL FATAL_ERROR("REF_END_TIME must be an integer number of internal time steps from the start time!")
! END IF
IF(MOD(ENKF_START_TIME-StartTime,IMDTI) /= ZeroTime) THEN
CALL FATAL_ERROR("ENKF_START_TIME must be an integer number of internal time steps from the start time!")
END IF
IF(MOD(ENKF_END_Time-StartTime,IMDTI) /= ZeroTime) THEN
CALL FATAL_ERROR("ENKF_END_Time must be an integer number of internal time steps from the start time!")
END IF
END SUBROUTINE ENKF_SET_TIME
! SUBROUTINE RRK_SET_STARTUP
! USE MOD_SET_TIME
! USE MOD_INPUT
! IMPLICIT NONE
! character(len=160) :: pathnfile
! TYPE(NCFILE), POINTER :: NCF
!
! OPEN(75,FILE='./err.dat')
!
! !RESTART FILE
! RESTART_FILE_NAME = trim(casename)//"_restart_0001.nc"
! pathnfile= trim(OUTPUT_DIR)//TRIM(RESTART_FILE_NAME)
!
! CALL SEARCH_FOR_LAST_MATCHING_NAME(PATHNFILE)
!
! ! OPEN THE FILE AND LOAD FOR STARTUP
! NCF => NEW_FILE()
! NCF%FNAME=trim(pathnfile)
! Call NC_OPEN(NCF)
! CALL NC_LOAD(NCF)
! NC_START => NCF
!
! Nullify(NCF)
! CALL SET_STARTUP_FILE_STACK(REF_START_TIME)
!
!! STARTUP_TYPE = STARTUP_TYPE_CRASHRESTART
!
!! STARTUP_UV_TYPE = STARTUP_TYPE_SETVALUES
!
!! STARTUP_TURB_TYPE = STARTUP_TYPE_SETVALUES
!
!! STARTUP_TS_TYPE = STARTUP_TYPE_SETVALUES
!
! END SUBROUTINE RRK_SET_STARTUP
!!=====================================================================================
!!
!! READ IN THE SIMULATION DATA AND CALCULATE THE EOF
!!
!!=====================================================================================
SUBROUTINE ENKF_REF
! USE CONTROL
! USE ENKFVAL
! USE MOD_NCTOOLS
IMPLICIT NONE
! TYPE(NCFILE), POINTER :: NCF
! TYPE(TIME) :: tTEST
! LOGICAL :: FOUND
! INTEGER :: NTIMES
! INTEGER IDUMMY, I, K, J
! INTEGER SS_DIM
! INTEGER STDIM
! ! REAL(DP) :: ELSD, USD, TSD, SSD
! INTEGER LWORK4, LDVT
! REAL(DP),ALLOCATABLE :: WORK4(:)
! REAL(DP) :: VT
! REAL(DP),ALLOCATABLE :: RKSF(:,:)
! REAL(DP),ALLOCATABLE :: RKSF1(:,:)
! ! REAL(DP),ALLOCATABLE :: SEOF(:,:)
! REAL(DP),ALLOCATABLE :: SFSF(:,:)
! REAL(DP),ALLOCATABLE :: SFD(:)
! REAL(DP),ALLOCATABLE :: SFU(:,:)
! REAL(DP),ALLOCATABLE :: STVAR(:)
! ! REAL(DP),ALLOCATABLE :: STSD(:)
! REAL(DP) :: SUM0
! INTEGER RCODE
! INTEGER IDUMMY1
! REAL(DP),ALLOCATABLE :: RRKU2(:,:,:)
! REAL(DP),ALLOCATABLE :: RRKV2(:,:,:)
! REAL(DP),ALLOCATABLE :: RRKEL2(:,:)
! REAL(DP),ALLOCATABLE :: RRKT2(:,:,:)
! REAL(DP),ALLOCATABLE :: RRKS2(:,:,:)
! REAL(DP),ALLOCATABLE :: STMEAN(:)
! TYPE(NCDIM), POINTER :: DIM
TYPE(NCVAR), POINTER :: VAR
! INTEGER :: IRANGE(2),CNT, STATUS
! character(len=160) :: pathnfile
! pathnfile = trim(INPUT_DIR)//trim(casename)//"_true_0001.nc"
!
! NCF => NEW_FILE()
! NCF%FNAME=trim(pathnfile)
!
! Call NC_OPEN(NCF)
! CALL NC_LOAD(NCF)
!
! ! GET THE TIME DIMENSION
! DIM => FIND_UNLIMITED(NCF,FOUND)
! IF(.not. FOUND) CALL FATAL_ERROR &
! & ("IN RRK_REF FILE",&
! & "FILE NAME: "//TRIM(NCF%FNAME),&
! &"COULD NOT FIND THE UNLIMITED DIMENSION")
! NTIMES = DIM%DIM
!
!
! ! TEST THE FILE START AND END
! TTEST = get_file_time(NCF,1)
!
! IF(REF_START_TIME < TTEST) CALL FATAL_ERROR &
! &("RRK_REF: REF_START_TIME IS BEFORE FIRST TIME IN FILE.",&
! & "FILE NAME:"//TRIM(NCF%FNAME))
!
!
! TTEST = get_file_time(NCF,NTIMES)
! IF(REF_END_TIME > TTEST) CALL FATAL_ERROR &
! &("RRK_REF: REF_START_TIME IS BEFORE FIRST TIME IN FILE.",&
! & "FILE NAME:"//TRIM(NCF%FNAME))
!
!
! ! GET THE START AND END INDEX
!
! IRANGE = 0
!
! CNT = 0
! DO WHILE (product(irange)==0)
! CNT = CNT +1
!
! IF (CNT > NTIMES) THEN
! CALL PRINT_TIME(REF_START_TIME,IPT,"RRK_REF START TIME")
! CALL PRINT_TIME(REF_END_TIME,IPT,"RRK_REF END TIME")
!
!
! CALL FATAL_ERROR&
! &("RRK_REF: CAN NOT FIND SPECIFIED START AND END TIME IN THE REFERENCE FILE!", &
! & "FILE NAME:"//TRIM(NCF%FNAME))
! END IF
!
!
! TTEST = get_file_time(NCF,CNT)
! IF (REF_START_TIME == TTEST) THEN
! IRANGE(1) = CNT
! END IF
!
! IF (REF_END_TIME == TTEST) THEN
! IRANGE(2) = CNT
! END IF
!
! END DO
!
!
! nrange = irange(2)-irange(1)+1
! IF(EL_ASSIM) THEN
!
! VAR => FIND_VAR(NCF,"zeta",FOUND)
! IF(.not. FOUND) CALL FATAL_ERROR &
! & ("IN RRK_REF FILE",&
! & "FILE NAME: "//TRIM(NCF%FNAME),&
! &"COULD NOT FIND THE VARIABLE 'zeta'")
!
! allocate(RRKEL_REF(MGL,nrange),stat=status)
! IF (STATUS /=0 ) CALL FATAL_ERROR("COULD NOT ALLOCATE MEMORY FOR RRKEL_REF")
! RRKEL_REF = 0.0_SP
!
! CALL nc_connect_avar(VAR,RRKEL_REF)
!
! CALL NC_READ_VAR(var,stkrng=irange,parallel=.false.)
!
!
! ENDIF
!
! IF(UV_ASSIM) THEN
!
! VAR => FIND_VAR(NCF,"u",FOUND)
! IF(.not. FOUND) CALL FATAL_ERROR &
! & ("IN RRK_REF FILE",&
! & "FILE NAME: "//TRIM(NCF%FNAME),&
! &"COULD NOT FIND THE VARIABLE 'u'")
!
! allocate(RRKU_REF(NGL,kbm1,nrange),stat=status)
! IF (STATUS /=0 ) CALL FATAL_ERROR("COULD NOT ALLOCATE MEMORY FOR RRKU_REF")
! RRKU_REF = 0.0_SP
! CALL nc_connect_avar(VAR,RRKU_REF)
!
! CALL NC_READ_VAR(var,stkrng=irange,parallel=.false.)
!
! VAR => FIND_VAR(NCF,"v",FOUND)
! IF(.not. FOUND) CALL FATAL_ERROR &
! & ("IN RRK_REF FILE",&
! & "FILE NAME: "//TRIM(NCF%FNAME),&
! &"COULD NOT FIND THE VARIABLE 'v'")
!
! allocate(RRKV_REF(NGL,kbm1,nrange),stat=status)
! IF (STATUS /=0 ) CALL FATAL_ERROR("COULD NOT ALLOCATE MEMORY FOR RRKV_REF")
! RRKV_REF = 0.0_SP
! CALL nc_connect_avar(VAR,RRKV_REF)
!
! CALL NC_READ_VAR(var,stkrng=irange,parallel=.false.)
!
!
! ENDIF
!
! IF(T_ASSIM) THEN
! VAR => FIND_VAR(NCF,"temp",FOUND)
! IF(.not. FOUND) CALL FATAL_ERROR &
! & ("IN RRK_REF FILE",&
! & "FILE NAME: "//TRIM(NCF%FNAME),&
! &"COULD NOT FIND THE VARIABLE 'temp'")
!
! allocate(RRKT_REF(MGL,kbm1,nrange),stat=status)
! IF (STATUS /=0 ) CALL FATAL_ERROR("COULD NOT ALLOCATE MEMORY FOR RRKT_REF")
! RRKT_REF = 0.0_SP
! CALL nc_connect_avar(VAR,RRKT_REF)
!
! CALL NC_READ_VAR(var,stkrng=irange,parallel=.false.)
!
! ENDIF
!
! IF(S_ASSIM) THEN
! VAR => FIND_VAR(NCF,"salinity",FOUND)
! IF(.not. FOUND) CALL FATAL_ERROR &
! & ("IN RRK_REF FILE",&
! & "FILE NAME: "//TRIM(NCF%FNAME),&
! &"COULD NOT FIND THE VARIABLE 'salinity'")
!
! allocate(RRKS_REF(MGL,kbm1,nrange),stat=status)
! IF (STATUS /=0 ) CALL FATAL_ERROR("COULD NOT ALLOCATE MEMORY FOR RRKS_REF")
! RRKS_REF = 0.0_SP
! CALL nc_connect_avar(VAR,RRKS_REF)
!
! CALL NC_READ_VAR(var,stkrng=irange,parallel=.false.)
!
! ENDIF
! CALL KILL_FILE(NCF)
!
! WRITE(IPT,*) 'Calculate the EOFs from the control run......'
!
! STDIM = 0
! IF(EL_ASSIM) STDIM = STDIM + MGL
! IF(UV_ASSIM) STDIM = STDIM + 2*NGL*KBM1
! IF(T_ASSIM) STDIM = STDIM + MGL*KBM1
! IF(S_ASSIM) STDIM = STDIM + MGL*KBM1
! SS_DIM = nrange
! LWORK4 = 5*SS_DIM
! LDVT = 1
!
! ALLOCATE(WORK4(LWORK4)) ;WORK4 = ZERO
! ALLOCATE(RKSF(STDIM,SS_DIM)) ;RKSF = ZERO
! ALLOCATE(RKSF1(STDIM,SS_DIM)) ;RKSF1 = ZERO
! ALLOCATE(SEOF(STDIM,SS_DIM)) ;SEOF = ZERO
! ALLOCATE(SFSF(SS_DIM,SS_DIM)) ;SFSF = ZERO
! ALLOCATE(SFD(SS_DIM)) ;SFD = ZERO
! ALLOCATE(SFU(SS_DIM,SS_DIM)) ;SFU = ZERO
! ALLOCATE(STVAR(STDIM)) ;STVAR = ZERO
! ALLOCATE(STSD(STDIM)) ;STSD = ZERO
! ALLOCATE(RRKU(NGL,KBM1)) ;RRKU = ZERO
! ALLOCATE(RRKV(NGL,KBM1)) ;RRKV = ZERO
! ALLOCATE(RRKEL(MGL)) ;RRKEL = ZERO
! ALLOCATE(RRKT(MGL,KBM1)) ;RRKT = ZERO
! ALLOCATE(RRKS(MGL,KBM1)) ;RRKS = ZERO
! ALLOCATE(RRKU2(NGL,KBM1,SS_DIM)) ;RRKU2 = ZERO
! ALLOCATE(RRKV2(NGL,KBM1,SS_DIM)) ;RRKV2 = ZERO
! ALLOCATE(RRKEL2(MGL,SS_DIM)) ;RRKEL2 = ZERO
! ALLOCATE(RRKT2(MGL,KBM1,SS_DIM)) ;RRKT2 = ZERO
! ALLOCATE(RRKS2(MGL,KBM1,SS_DIM)) ;RRKS2 = ZERO
! ALLOCATE(STMEAN(STDIM)) ;STMEAN = ZERO
! STORE IN RKSF
! DO I = 1, SS_DIM
!
! IDUMMY = 0
!
! IF(EL_ASSIM) THEN
! DO J = 1, MGL
! IDUMMY = IDUMMY + 1
! RKSF(IDUMMY,I) = RRKEL_REF(J,I)
! ENDDO
! ENDIF
!
! IF(UV_ASSIM) THEN
! DO K = 1, KBM1
! DO J = 1, NGL
! IDUMMY = IDUMMY + 1
! RKSF(IDUMMY,I) = RRKU_REF(J,K,I)
! ENDDO
! ENDDO
!
! DO K = 1, KBM1
! DO J = 1, NGL
! IDUMMY = IDUMMY + 1
! RKSF(IDUMMY,I) = RRKV_REF(J,K,I)
! ENDDO
! ENDDO
! ENDIF
!
! IF(T_ASSIM) THEN
! DO K = 1, KBM1
! DO J = 1, MGL
! IDUMMY = IDUMMY + 1
! RKSF(IDUMMY,I) = RRKT_REF(J,K,I)
! ENDDO
! ENDDO
! ENDIF
!
! IF(S_ASSIM) THEN
! DO K = 1, KBM1
! DO J = 1, MGL
! IDUMMY = IDUMMY + 1
! RKSF(IDUMMY,I) = RRKS_REF(J,K,I)
! ENDDO
! ENDDO
! ENDIF
!
! ENDDO
!
!
! ! CALCULATE THE MEAN AND THE STANDARD DEVIATION
!
! DO I=1, STDIM
! SUM0=0.0_DP
! DO J=1, SS_DIM
! SUM0=SUM0+RKSF(I,J)
! ENDDO
!
! STMEAN(I) = SUM0/DBLE(SS_DIM)
! SUM0=0.0_DP
! DO J=1, SS_DIM
! RKSF1(I,J)=(RKSF(I,J)-STMEAN(I))
! SUM0=SUM0+RKSF1(I,J)**2.0_DP
! ENDDO
!
! STVAR(I)=SUM0
! STSD(I)=DSQRT(SUM0/DBLE(SS_DIM-1))
! ENDDO
!
!
! IDUMMY = 0
! IDUMMY1 = 0
! ELSD = 0.0_DP
! USD = 0.0_DP
! TSD = 0.0_DP
! SSD = 0.0_DP
!
! IF(EL_ASSIM) THEN
! SUM0=0.0_DP
! DO I=1, MGL
! IDUMMY = IDUMMY + 1
! SUM0 = SUM0 + STVAR(IDUMMY)
! ENDDO
!
! ELSD=DSQRT(SUM0/DBLE(MGL)/DBLE(SS_DIM-1))
!
! DO I=1, MGL
! DO J=1, SS_DIM
! RKSF1(I,J)=RKSF1(I,J)/ELSD/DSQRT(DBLE(SS_DIM-1))
! ENDDO
! ENDDO
!
! IDUMMY1 = IDUMMY
! ENDIF
!
! IF(UV_ASSIM) THEN
! SUM0=0.0_DP
! DO K=1, 2*KBM1
! DO I=1, NGL
! IDUMMY = IDUMMY + 1
! SUM0 = SUM0 +STVAR(IDUMMY)
! ENDDO
! ENDDO
! USD=DSQRT(SUM0/DBLE(NGL*KBM1)/DBLE(SS_DIM-1))
!
! DO I=IDUMMY1+1, IDUMMY1+2*NGL*KBM1
! DO J=1, SS_DIM
! RKSF1(I,J)=RKSF1(I,J)/USD/DSQRT(DBLE(KBM1))/DSQRT(DBLE(SS_DIM-1))
! ENDDO
! ENDDO
!
! IDUMMY1 = IDUMMY
! ENDIF
!
! IF(T_ASSIM) THEN
! SUM0=0.0_DP
! DO K=1, KBM1
! DO I=1, MGL
! IDUMMY = IDUMMY + 1
! SUM0 = SUM0 +STVAR(IDUMMY)
! ENDDO
! ENDDO
! TSD=DSQRT(SUM0/DBLE(MGL*KBM1)/DBLE(SS_DIM-1))
!
! DO I=IDUMMY1+1, IDUMMY1+MGL*KBM1
! DO J=1, SS_DIM
! RKSF1(I,J)=RKSF1(I,J)/TSD/DSQRT(DBLE(SS_DIM-1))
! ENDDO
! ENDDO
!
! IDUMMY1 = IDUMMY
! ENDIF
!
! IF(S_ASSIM) THEN
! SUM0=0.0_DP
! DO K=1, KBM1
! DO I=1, MGL
! IDUMMY = IDUMMY + 1
! SUM0 = SUM0 +STVAR(IDUMMY)
! ENDDO
! ENDDO
! SSD=DSQRT(SUM0/DBLE(MGL*KBM1)/DBLE(SS_DIM-1))
!
! DO I=IDUMMY1+1, IDUMMY1+MGL*KBM1
! DO J=1, SS_DIM
! RKSF1(I,J)=RKSF1(I,J)/SSD/DSQRT(DBLE(SS_DIM-1))
! ENDDO
! ENDDO
!
! IDUMMY1 = IDUMMY
! ENDIF
!
!
! ! CALCULATE THE EOFs BY SVD OF RKSF1' RKSF1 = C LAMBDA C', RKSF1 RKSF1' = E LAMBDA E', E = RKSF1 C LAMBDA ^-1/2
!
! DO I=1, SS_DIM
! DO J=1, SS_DIM
! SUM0=0.0_DP
! DO K=1, STDIM
! SUM0 = SUM0 + RKSF1(K,I)*RKSF1(K,J)
! ENDDO
! SFSF(I,J) =SUM0
! ENDDO
! ENDDO
!
! CALL DGESVD('A','N',SS_DIM,SS_DIM,SFSF,SS_DIM,SFD,SFU,SS_DIM,VT,LDVT,WORK4,LWORK4,RCODE)
!
! OPEN(72,FILE=TRIM(OUTPUT_DIR)//'/rrktemp/'//'eigenvalue.dat')
!
! DO I=1, SS_DIM
! WRITE(72,'(I5,E15.7)') I, SFD(I)
! ENDDO
! DO I=1, SS_DIM
! DO J=1, SS_DIM
! SFU(I,J)=SFU(I,J)/DSQRT(SFD(J))
! ENDDO
! ENDDO
!
! DO I=1, STDIM
! DO J=1, SS_DIM
! SUM0=0.0_DP
! DO K=1, SS_DIM
! SUM0=SUM0+RKSF1(I,K)*SFU(K,J)
! ENDDO
! SEOF(I,J)=SUM0
! ENDDO
! ENDDO
! write(1000,'(871f12.4)') (seof(i,1),i=1,stdim)
! DEALLOCATE(WORK4,RKSF,RKSF1,SFSF,SFD,SFU,STVAR)
! DEALLOCATE(RRKU2,RRKV2,RRKEL2,RRKT2,RRKS2)
END SUBROUTINE ENKF_REF
!------------------------------------------------------------------------------|
! Main program to calculate the reduced rank kalman gain matrix |
!------------------------------------------------------------------------------|
! SUBROUTINE RRK_RRK(CHOICE)
!
! USE LIMS
! USE CONTROL
! USE ENKFVAL
! IMPLICIT NONE
!
! INTEGER SS_DIM
! INTEGER STDIM
! INTEGER I_EOF
! INTEGER I,J,K
! INTEGER CHOICE
! INTEGER NLOC
! INTEGER RCODE
! INTEGER IT
! INTEGER ILAST
! INTEGER IDUMMY
! INTEGER SS_TOT
! INTEGER,ALLOCATABLE :: STLOC(:)
! REAL(DP) ERR1
! REAL(DP),ALLOCATABLE :: KAL(:,:)
! REAL(DP),ALLOCATABLE :: HEOF(:,:)
! REAL(DP),ALLOCATABLE :: R(:,:)
! REAL(DP),ALLOCATABLE :: HTR(:,:)
! REAL(DP),ALLOCATABLE :: R2(:,:)
! REAL(DP),ALLOCATABLE :: PHT(:,:)
! REAL(DP),ALLOCATABLE :: RALPHA(:,:)
! REAL(DP),ALLOCATABLE :: BETA2(:,:)
! REAL(DP),ALLOCATABLE :: GAMMA(:,:)
! REAL(DP),ALLOCATABLE :: GAMMA2(:,:)
! REAL(DP),ALLOCATABLE :: TRANS2(:,:)
!
!! WORK ARRAYS FOR LAPACK SUBROUTINES
! INTEGER LWORK, LWORK2, LDVT
! INTEGER,ALLOCATABLE :: IPIV(:)
! INTEGER,ALLOCATABLE :: IPIV2(:)
! REAL(DP),ALLOCATABLE :: WORK(:)
! REAL(DP),ALLOCATABLE :: WORK2(:)
! REAL(DP) :: VT
!
!
!! CHARACTER STRINGS
! CHARACTER(LEN=80) KALNAME
! CHARACTER(LEN=80) FILENAME
! CHARACTER(LEN=4) STRCYC
! CHARACTER(LEN=8) CH8
! CHARACTER(LEN=80) INAME
! CHARACTER(LEN=80) INAME2
!
!! REAL(DP),ALLOCATABLE :: SEOF(:,:)
!! REAL(DP),ALLOCATABLE :: STSD(:)
!! REAL(DP) :: ELSD, USD, TSD, SSD
! REAL(DP),ALLOCATABLE :: HU(:,:)
! REAL(DP),ALLOCATABLE :: HUL(:,:)
! REAL(DP),ALLOCATABLE :: EVAL(:)
! REAL(DP) :: RRKSUM
!
! STDIM = 0
! IF(EL_ASSIM) STDIM = STDIM + MGL
! IF(UV_ASSIM) STDIM = STDIM + 2*NGL*KBM1
! IF(T_ASSIM) STDIM = STDIM + MGL*KBM1
! IF(S_ASSIM) STDIM = STDIM + MGL*KBM1
!
! SS_DIM = ENKF_NENS
! LWORK = 4*SS_DIM
! LWORK2 = 4*ENKF_NOBSMAX
! SS_TOT = nrange
! LDVT = 1
! ILAST = 10 ! THE NUMBER OF ITERATION IN DOUBLING ALGORITHM
!
!! TEMPORARILY ALLOCATE ARRY TO ARRYS
!
!
! ALLOCATE(STLOC(ENKF_NOBSMAX)) ;STLOC = 0
! ALLOCATE(IPIV(SS_DIM)) ;IPIV = 0
! ALLOCATE(IPIV2(ENKF_NOBSMAX)) ;IPIV2 = 0
! ALLOCATE(KAL(SS_DIM,ENKF_NOBSMAX)) ;KAL = ZERO
! ALLOCATE(HEOF(ENKF_NOBSMAX,SS_DIM)) ;HEOF = ZERO
! ALLOCATE(R(ENKF_NOBSMAX,ENKF_NOBSMAX)) ;R = ZERO
! ALLOCATE(HTR(SS_DIM,ENKF_NOBSMAX)) ;HTR = ZERO
! ALLOCATE(R2(ENKF_NOBSMAX,ENKF_NOBSMAX)) ;R2 = ZERO
! ALLOCATE(PHT(SS_DIM,ENKF_NOBSMAX)) ;PHT = ZERO
! ALLOCATE(RALPHA(SS_DIM,SS_DIM)) ;RALPHA = ZERO
! ALLOCATE(BETA2(SS_DIM,SS_DIM)) ;BETA2 = ZERO
! ALLOCATE(GAMMA(SS_DIM,SS_DIM)) ;GAMMA = ZERO
! ALLOCATE(GAMMA2(SS_DIM,SS_DIM)) ;GAMMA2 = ZERO
! ALLOCATE(TRANS2(SS_DIM,SS_DIM)) ;TRANS2 = ZERO
! ALLOCATE(WORK(LWORK)) ;WORK = ZERO
! ALLOCATE(WORK2(LWORK2)) ;WORK2 = ZERO
!! ALLOCATE(SEOF(STDIM,SS_DIM)) ;SEOF = ZERO
!! ALLOCATE(STSD(STDIM)) ;STSD = ZERO
! ALLOCATE(HU(ENKF_NOBSMAX,SS_TOT)) ;HU = ZERO
! ALLOCATE(HUL(ENKF_NOBSMAX,SS_TOT)) ;HUL = ZERO
! ALLOCATE(EVAL(SS_TOT)) ;EVAL = ZERO
!
!
!! END ALLOCATION
!
! IF(CHOICE == 1) THEN
! IF(MSR) WRITE(IPT,*) 'Starting perturb the base state in the eigenvector direction......'
!
! CALL PERTURB
!
! DEALLOCATE(STLOC,IPIV,IPIV2,KAL,HEOF,R,HTR,R2,PHT,RALPHA,TRANS2)
! DEALLOCATE(BETA2,GAMMA,GAMMA2,WORK,WORK2,HU,HUL,EVAL)
!
! RETURN
! ENDIF
!
! IF(CHOICE == 2) THEN
! IF(MSR) WRITE(IPT,*) 'Starting calculate the linearized model matrix in the reduced space......'
!! CALCULATE THE LINEARIZED MODEL MATRIX IN THE REDUCED SPACE
! CALL MREDUCED
!
! IF(MSR) THEN
!! OUTPUT LINEAR TRANSITION MATRIX
! FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/Amatr.dat'
! OPEN(INORRK,FILE=TRIM(FILENAME),FORM='UNFORMATTED')
! DO J=1, SS_DIM
! WRITE(INORRK) (TRANS(I,J), I=1, SS_DIM)
! ENDDO
! CLOSE(INORRK)
!
! ! ALSO IN ASCII
! FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/Amatr1.txt'
! OPEN(INORRK, FILE=TRIM(FILENAME),FORM='FORMATTED')
! DO I=1, SS_DIM
! DO J=1, SS_DIM
! WRITE(INORRK,*) TRANS(J,I)
! ENDDO
! ENDDO
! CLOSE(INORRK)
! ENDIF
!
! DEALLOCATE(STLOC,IPIV,IPIV2,KAL,HEOF,R,HTR,R2,PHT,RALPHA,TRANS2)
! DEALLOCATE(BETA2,GAMMA,GAMMA2,WORK,WORK2,HU,HUL,EVAL)
!! DEALLOCATE(BETA2,GAMMA,GAMMA2,WORK,WORK2,SEOF,STSD,HU,HUL,EVAL)
!! DEALLOCATE(STTEMP0,STTEMP1,STEOF,SDEOF,TRANS,RRKEL,RRKU,RRKV,RRKT,RRKS)
! RETURN
! ENDIF
!
! IF(CHOICE == 4) THEN
! CALL RRK_ALLOC_VAR
! WRITE(IPT,*) 'Calculate the Kalman gain matrix by doubling algorith......'
! KALNAME = TRIM(OUTPUT_DIR)//'/rrktemp/rrK.dat'
!
!! READ THE EIGENVALUES AND SET THE INITIAL ERROR COVARIANCE
! INAME = TRIM(OUTPUT_DIR)//'/rrktemp/eigenvalue.dat'
! OPEN(INORRK,FILE=INAME,STATUS='OLD')
!
! WRITE(IPT,*) 'ss_tot:', ss_tot, ss_dim
!
! DO I=1, SS_TOT
! READ(INORRK,*) J, EVAL(I)
! ENDDO
!
! DO I=1, SS_DIM
! DO J=1, SS_DIM
! GAMMA2(I,J) = 0._SP
! ENDDO
! GAMMA2(I,I) = EVAL(I)*RRK_PSCALE
! ENDDO
! CLOSE(INORRK)
!
!! READ THE LINEAR TRANSITION MATRIX
! FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/Amatr.dat'
! OPEN(INORRK,FILE=TRIM(FILENAME),FORM='UNFORMATTED')
! DO J=1, SS_DIM
! READ(INORRK) (TRANS(I,J), I=1, SS_DIM)
! ENDDO
! CLOSE(INORRK)
!
!! READ THE OBSERVATION NUMBER AND LOCATION
! CALL READOBS(STLOC,NLOC)
!
! print *, 'obs=', STLOC(1), NLOC
!
!! SET RALPHA
! DO I=1, SS_DIM
! DO J=1, SS_DIM
! TRANS2(I,J)=TRANS(J,I)
! ENDDO
!! WRITE(IPT,*) 'Amatr:', I, TRANS(I,I)
! ENDDO
!
! print *, 'debug - 4'
!
!! CALCULATE OBSERVATION MATRIX (EOF TO OBS TRANSFORMATION) H*D*EOFs
! DO J=1, SS_DIM
! DO I=1, STDIM
! STTEMP1(I) = SEOF(I,J)
! ENDDO
!
! DO I=1, NLOC
! HEOF(I,J)=STTEMP1(STLOC(I))*STSD(STLOC(I))
! WRITE(IPT,*) 'Heof:', I,J,HEOF(I,J)
! ENDDO
! ENDDO
!
!
! FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/ObsOp.dat'
! OPEN(INORRK,FILE=FILENAME, FORM='UNFORMATTED')
! DO J=1, SS_DIM
! WRITE(INORRK) (HEOF(I,J), I=1, NLOC)
! ENDDO
! CLOSE(INORRK)
! FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/Heof.txt'
! OPEN(INORRK,FILE=FILENAME, FORM='FORMATTED')
! DO J=1, SS_DIM
! DO I=1, NLOC
! WRITE(INORRK,*) I,J,HEOF(I,J)
! ENDDO
! ENDDO
! CLOSE(INORRK)
!
!! OUTPUT RALPHA MATRIX, WHICH IS TRANSPOSE OF TRANS
! FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/Alpha.dat'
! OPEN(INORRK,FILE=FILENAME, FORM='UNFORMATTED')
! DO J=1, SS_DIM
! WRITE(INORRK) (TRANS2(I,J), I=1, SS_DIM)
! ENDDO
! CLOSE(INORRK)
!
!! OUTPUT GAMMA MATRIX
! FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/Gamma.dat'
! OPEN(INORRK,FILE=FILENAME, FORM='UNFORMATTED')
! DO J=1, SS_DIM
! WRITE(INORRK) (GAMMA2(I,J), I=1, SS_DIM)
! ENDDO
! CLOSE(INORRK)
!
!! CALCULATE REPRESENTATIVE OBSERVATION ERROR DIRECTLY USING EOFs IN UNRESOLVED SUBSPACE
! DO I=SS_DIM+1, SS_TOT
!! CALL READEOF(I,2)
! DO J=1, STDIM
! STTEMP1(J) = SEOF(J,I)
! ENDDO
!
! DO J=1, NLOC
! HU(J,I) = STTEMP1(STLOC(J))*STSD(STLOC(J))
! ENDDO
! ENDDO
! DO I=1, NLOC
! RRKSUM=0.0_DP
! DO J=SS_DIM+1, SS_TOT
! HUL(I,J)=RRKSUM+HU(I,J)*EVAL(J)
! ENDDO
! ENDDO
!
! DO I=1, NLOC
! DO J=1, NLOC
! RRKSUM=0.0_DP
! DO K=SS_DIM+1,SS_TOT
! RRKSUM = RRKSUM + HUL(I,K)*HU(J,K)
! ENDDO
! R(I,J) = RRKSUM
! ENDDO
!
!! ASSUME THAT THE MEASUREMENT ERROR IS 1% OF THE ONE STANDARD DEVIATION OF THE CONTROL RUN
! IF(RRK_OPTION == 0) THEN
! R(I,I) = R(I,I) + (STSD(STLOC(I))*0.01_DP)**2.0_DP
! ELSE
! R(I,I) = R(I,I) + (STSD(STLOC(I))/DSQRT(DBLE(KBM1))*0.01_DP)**2.0_DP
! ENDIF
! ENDDO
!
! DO I=1, NLOC
! DO J=1, I-1
! R(I,J) = 0.5_DP*(R(I,J)+R(J,I))
! R(J,I) = R(I,J)
! ENDDO
! ENDDO
! FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/Robs.dat'
! OPEN(INORRK,FILE=FILENAME, FORM='UNFORMATTED')
! DO J=1, NLOC
! WRITE(INORRK) (R(I,J), I=1, NLOC)
! ENDDO
! CLOSE(INORRK)
! FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/Robs.txt'
! OPEN(INORRK,FILE=FILENAME, FORM='FORMATTED')
! DO J=1, NLOC
! DO I=1, NLOC
! WRITE(INORRK,*) I,J, R(I,J)
! ENDDO
! ENDDO
! CLOSE(INORRK)
!
!! INVERT OBSERVATION ERROR COVARIANCE
! CALL DGETRF(NLOC,NLOC,R,ENKF_NOBSMAX,IPIV2,RCODE)
! IF(RCODE/=0) WRITE(IPT,*) 'error in computing LU factorization'
! call DGETRI(NLOC,R,ENKF_NOBSMAX,IPIV2,WORK2,LWORK2,RCODE)
! IF(RCODE/=0) WRITE(IPT,*) 'error in inverting the matrix'
!
!! FORM BETA MATRIX = H_T*R**(-1)*H, WHERE H IS A NORMALIZED MATRIX H = H_ORIG*S
! CALL DGEMM('t','n',SS_DIM,NLOC,NLOC,1.0d0,HEOF,ENKF_NOBSMAX,R,ENKF_NOBSMAX,0.0d0,HTR,SS_DIM)
! CALL DGEMM('n','n',SS_DIM,SS_DIM,NLOC,1.0d0,HTR,SS_DIM,HEOF,ENKF_NOBSMAX,0.0d0,BETA2,SS_DIM)
! FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/Beta.dat'
! OPEN(INORRK,FILE=FILENAME, FORM='UNFORMATTED')
! DO J=1, SS_DIM
! WRITE(INORRK) (BETA2(I,J), I=1, SS_DIM)
! ENDDO
! CLOSE(INORRK)
!
! WRITE(IPT,*) 'Running the doubling algorithm......'
!
! DO IT =1, ILAST
!
!! READ GAMMA FROM FILE
! FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/Gamma.dat'
! OPEN(INORRK,FILE=TRIM(FILENAME),FORM='UNFORMATTED')
! DO J=1, SS_DIM
! READ(INORRK) (GAMMA(I,J), I=1, SS_DIM)
! ENDDO
! CLOSE(INORRK)
!
!! READ BETA FROM FILE
! FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/Beta.dat'
! OPEN(INORRK,FILE=TRIM(FILENAME),FORM='UNFORMATTED')
! DO J=1, SS_DIM
! READ(INORRK) (BETA2(I,J), I=1, SS_DIM)
! ENDDO
! CLOSE(INORRK)
!
!! COMPUTE EYE + BETA*GAMMA = STORE TEMPORARILY IN RALPHA (RALPHA = EYE + BETA*GAMMA)
!
! DO I=1, SS_DIM
! DO J=1, SS_DIM
! RALPHA(I,J)=0.0_DP
! ENDDO
! RALPHA(I,I)=1.0_DP
! ENDDO
! CALL DGEMM('n','n',SS_DIM,SS_DIM,SS_DIM,1.0d0,BETA2,SS_DIM,GAMMA,SS_DIM,1.0d0,RALPHA,SS_DIM)
!
!! COMPUTE INVERSE OF ABOBE RALPHA (=EYE+BETA*GAMMA) BY LAPACK ROUTINES AND STORE IT IN BETA (BETA=(EYE+BETA*GAMMA)**(-1))
! CALL DGETRF(SS_DIM,SS_DIM,RALPHA,SS_DIM,IPIV,RCODE)
!
! IF(RCODE/=0) WRITE(IPT,*) 'error in computing LU factorization'
! CALL DGETRI(SS_DIM,RALPHA,SS_DIM,IPIV,WORK,LWORK,RCODE)
! IF(RCODE/=0) WRITE(IPT,*) 'error in inverting the matrix'
! DO I=1, SS_DIM
! DO J=1, SS_DIM
! BETA2(I,J) = RALPHA(I,J)
! ENDDO
! ENDDO
!
!! READ RALPHA FROM FILE
! FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/Alpha.dat'
! OPEN(INORRK,FILE=TRIM(FILENAME),FORM='UNFORMATTED')
! DO J=1, SS_DIM
! READ(INORRK) (RALPHA(I,J), I=1, SS_DIM)
! ENDDO
! CLOSE(INORRK)
!
!! COMPUTE PRODUCT (EYE+BETA*GAMMA)**(-1)*RALPHA (GAMMA = BETA'*RALPHA)
! CALL DGEMM('n','n',SS_DIM,SS_DIM,SS_DIM,1.0d0,BETA2,SS_DIM,RALPHA,SS_DIM,0.0d0,GAMMA,SS_DIM)
!
!! OUTPUT THIS PRODUCT TO TEMPORARY FILE
! FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/temp.02'
! OPEN(INORRK,FILE=TRIM(FILENAME),FORM='UNFORMATTED')
! DO J=1, SS_DIM
! WRITE(INORRK) (GAMMA(I,J), I=1, SS_DIM)
! ENDDO
! CLOSE(INORRK)
!
!! COMPUTE PRODUCT RALPHA*(EYE+BETA*GAMMA)**(-1) AND STORE IN GAMMA (GAMMA = RALPHA*BETA')
! CALL DGEMM('n','n',SS_DIM,SS_DIM,SS_DIM,1.0d0,RALPHA,SS_DIM,BETA2,SS_DIM,0.0d0,GAMMA,SS_DIM)
!
!! READ BACK OLD FILE
! FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/Beta.dat'
! OPEN(INORRK,FILE=TRIM(FILENAME),FORM='UNFORMATTED')
! DO J=1, SS_DIM
! READ(INORRK) (Beta2(I,J), I=1, SS_DIM)
! ENDDO
! CLOSE(INORRK)
!
!! COMPUTE PRODUCT RALPHA*(EYE+BETA*GAMMA)**(-1)*BETA (RALPHA' = GAMMA'*BETA')
! CALL DGEMM('n','n',SS_DIM,SS_DIM,SS_DIM,1.0d0,GAMMA,SS_DIM,BETA2,SS_DIM,0.0d0,RALPHA,SS_DIM)
!
!! READ BACK OLD RALPHA AND PUT IN GAMMA
! FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/Alpha.dat'
! OPEN(INORRK,FILE=TRIM(FILENAME),FORM='UNFORMATTED')
! DO J=1, SS_DIM
! READ(INORRK) (GAMMA(I,J), I=1, SS_DIM)
! ENDDO
! CLOSE(INORRK)
!
!! COMPUTE BETA+RALPHA*(EYE+BETA*GAMMA)**(-1)*BETA*RALPHA_T (BETA = BETA+RALPHA'*GAMMA_T)
! CALL DGEMM('n','t',SS_DIM,SS_DIM,SS_DIM,1.0d0,RALPHA,SS_DIM,GAMMA,SS_DIM,1.0d0,BETA2,SS_DIM)
!
!! MAKE SURE BETA IS SYMMETRIC
! DO I=1, SS_DIM
! DO J=1, I-1
! BETA2(I,J) = 0.5_DP*(BETA2(I,J)+BETA2(J,I))
! BETA2(J,I) = BETA2(I,J)
! ENDDO
! ENDDO
!
!! SAVE THIS NEW BETA TO FILE
! FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/Beta.dat'
! OPEN(INORRK,FILE=FILENAME, FORM='UNFORMATTED')
! DO J=1, SS_DIM
! WRITE(INORRK) (BETA2(I,J), I=1, SS_DIM)
! ENDDO
! CLOSE(INORRK)
!
!! READ BACK OLD GAMMA
! FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/Gamma.dat'
! OPEN(INORRK,FILE=FILENAME, FORM='UNFORMATTED')
! DO J=1, SS_DIM
! READ(INORRK) (GAMMA(I,J), I=1, SS_DIM)
! ENDDO
! CLOSE(INORRK)
!
!! READ TEMPORARY FILE INTO RALPHA
! FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/temp.02'
! OPEN(INORRK,FILE=FILENAME, FORM='UNFORMATTED')
! DO J=1, SS_DIM
! READ(INORRK) (RALPHA(I,J), I=1, SS_DIM)
! ENDDO
! CLOSE(INORRK)
!
!! COMPUTE GAMMA*(EYE+BETA*GAMMA)**(-1)*RALPHA (BETA'=GAMMA*RALPHA')
! CALL DGEMM('n','n',SS_DIM,SS_DIM,SS_DIM,1.0d0,GAMMA,SS_DIM,RALPHA,SS_DIM,0.0d0,BETA2,SS_DIM)
!
!! READ BACK OLD RALPHA
! FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/Alpha.dat'
! OPEN(INORRK,FILE=FILENAME, FORM='UNFORMATTED')
! DO J=1, SS_DIM
! READ(INORRK) (RALPHA(I,J), I=1, SS_DIM)
! ENDDO
! CLOSE(INORRK)
!
!! COMPUTE GAMMA + RALPHA_T*GAMMA*(EYE+BETA*GAMMA)**(-1)*RALPHA (GAMMA = GAMMA+RALPHA_T*BETA')
! CALL DGEMM('t','n',SS_DIM,SS_DIM,SS_DIM,1.0d0,RALPHA,SS_DIM,BETA2,SS_DIM,1.0d0,GAMMA,SS_DIM)
!
!! ENSURE GAMMA IS SYMMETRIC
! DO I=1, SS_DIM
! DO J=1, I-1
! GAMMA(I,J) = 0.5_DP*(GAMMA(I,J)+GAMMA(J,I))
! GAMMA(J,I) = GAMMA(I,J)
! ENDDO
! ENDDO
!
!! SAVE THIS NEW GAMMA TO FILE
! FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/Gamma.dat'
! OPEN(INORRK,FILE=FILENAME, FORM='UNFORMATTED')
! DO J=1, SS_DIM
! WRITE(INORRK) (GAMMA(I,J), I=1, SS_DIM)
! ENDDO
! CLOSE(INORRK)
! FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/Pfctr.txt'
! OPEN(INORRK,FILE=FILENAME, FORM='FORMATTED')
! DO J=1, SS_DIM
! DO I=1, SS_DIM
! WRITE(INORRK,*) I,J,GAMMA(I,J)
! ENDDO
! ENDDO
! CLOSE(INORRK)
!
!! READ TEMPORARY FILE INTO GAMMA
! FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/temp.02'
! OPEN(INORRK,FILE=FILENAME, FORM='UNFORMATTED')
! DO J=1, SS_DIM
! READ(INORRK) (GAMMA(I,J), I=1, SS_DIM)
! ENDDO
! CLOSE(INORRK)
!
!! COMPUTE PRODUCT RALPHA*(EYE+BETA*GAMMA)**(-1)*RALPHA (BETA = RALPHA*GAMMA)
! CALL DGEMM('n','n',SS_DIM,SS_DIM,SS_DIM,1.0d0,RALPHA,SS_DIM,GAMMA,SS_DIM,0.0d0,BETA2,SS_DIM)
!
!! SAVE THIS NEW RALPHA TO FILE
! FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/Alpha.dat'
! OPEN(INORRK,FILE=FILENAME, FORM='UNFORMATTED')
! DO J=1, SS_DIM
! WRITE(INORRK) (BETA2(I,J), I=1, SS_DIM)
! ENDDO
! CLOSE(INORRK)
!
!!===================================================================================================
!! END OF ITERATION FOR THE DOUBLING ALGORITHM
! ENDDO
!!===================================================================================================
!
! WRITE(IPT,*) 'Setting up the Kalman gain matrix in the reduced space......'
!
!! READ IN PF FROM DOUBLING ALGORITHM
! FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/Gamma.dat'
! OPEN(INORRK,FILE=FILENAME, FORM='UNFORMATTED')
! DO J=1, SS_DIM
! READ(INORRK) (GAMMA(I,J), I=1, SS_DIM)
! ENDDO
! CLOSE(INORRK)
!
!! READ IN OBSERVATION OPERATOR
! FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/ObsOp.dat'
! OPEN(INORRK,FILE=FILENAME, FORM='UNFORMATTED')
! DO J=1, SS_DIM
! READ(INORRK) (HEOF(I,J), I=1, NLOC)
! ENDDO
! CLOSE(INORRK)
!
!! READ OBSERVATION ERROR COVARIANCE
! FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/Robs.dat'
! OPEN(INORRK,FILE=FILENAME, FORM='UNFORMATTED')
! DO J=1, NLOC
! READ(INORRK) (R(I,J), I=1, NLOC)
! ENDDO
! CLOSE(INORRK)
!
!! COMPUTE P*H^T (PHT = GAMMA * HEOF^T)
! CALL DGEMM('n','t',SS_DIM,NLOC,SS_DIM,1.0d0,GAMMA,SS_DIM,HEOF,ENKF_NOBSMAX,0.0d0,PHT,SS_DIM)
!
!! COMPUTE H*P*H^T (R2 = H*P*HT)
! CALL DGEMM('n','n',NLOC,NLOC,SS_DIM,1.0d0,HEOF,ENKF_NOBSMAX,PHT,SS_DIM,0.0d0,R2,ENKF_NOBSMAX)
!
!! OUPUT BOTH FORCAST AND OBSERVATION ERROR STANDARD DEVIATION IN ASCII
! FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/Fctobserr.txt'
! OPEN(INORRK,FILE=FILENAME, FORM='FORMATTED')
! DO I=1, NLOC
! WRITE(INORRK,*) STLOC(I),DSQRT(R2(I,I)),DSQRT(R(I,I))
! ENDDO
! CLOSE(INORRK)
!
!! ADD R + (H P H^t) ---> R
! DO I = 1, NLOC
! DO J=1, NLOC
! R(J,I) = R(J,I) + R2(J,I)
! ENDDO
! ENDDO
!
!! INVERT R (H*P*H^T+R)
! CALL DGETRF(NLOC,NLOC,R,ENKF_NOBSMAX,IPIV2,RCODE)
! IF(RCODE/=0) WRITE(IPT,*) 'error in computing LU factorization'
! CALL DGETRI(NLOC,R,ENKF_NOBSMAX,IPIV2,WORK2,LWORK2,RCODE)
! IF(RCODE/=0) WRITE(IPT,*) 'error in inverting the matrix'
!
!! COMPUTE KALMAN GAIN: K = P*HT*(H*P*H^T+R)^(-1)
! CALL DGEMM('n','n',SS_DIM,NLOC,NLOC,1.0d0,PHT,SS_DIM,R,ENKF_NOBSMAX,0.0d0,KAL,SS_DIM)
!
!! OUTPUT RESULT
! OPEN(INORRK,FILE=KALNAME, FORM='UNFORMATTED')
! DO J=1, NLOC
! WRITE(INORRK) (KAL(I,J), I=1, SS_DIM)
! ENDDO
! CLOSE(INORRK)
! FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/rrK.txt'
! OPEN(INORRK,FILE=FILENAME, FORM='FORMATTED')
! DO J=1, NLOC
! DO I=1, SS_DIM
! WRITE(INORRK,*) I,J,KAL(I,J)
! ENDDO
! ENDDO
! CLOSE(INORRK)
!
!! OUTPUT KALMAN GAIN MATRIX IN THE FULL SPACE: D*ER*KAL
!! DO J=1, SS_DIM
!! CALL READEOF(J,2)
!! DO I=1, STDIM
!! SEOF(I,J)=STTEMP1(I)
!! ENDDO
!! ENDDO
!
! DO J=1, NLOC
! WRITE(STRCYC,'(I4.4)') J
! DO I=1, STDIM
! RRKSUM = 0.0_DP
! DO K=1, SS_DIM
! RRKSUM = RRKSUM + STSD(I)*SEOF(I,K)*KAL(K,J)
! ENDDO
! STTEMP1(I) = RRKSUM
! ENDDO
! FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/'//'BigK'//STRCYC//'.txt'
! OPEN(INORRK,FILE=FILENAME,FORM='FORMATTED')
! DO I=1, STDIM
! WRITE(INORRK,*) STTEMP1(I)
! ENDDO
! CLOSE(INORRK)
! ENDDO
!
!! COMPUTE ANALYSIS ERROR COVARIANCE: P_AN1 = (I-K*H)*P (JUST AS DIAGNOSIS)
! CALL DGEMM('n','t',SS_DIM,SS_DIM,NLOC,-1.0d0,KAL,SS_DIM,PHT,SS_DIM,1.0d0,GAMMA,SS_DIM)
! FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/Panlr.txt'
! OPEN(INORRK,FILE=FILENAME, FORM='FORMATTED')
! DO J=1, SS_DIM
! DO I=1, SS_DIM
! WRITE(INORRK,*) I,J,GAMMA(I,J)
! ENDDO
! ENDDO
! CLOSE(INORRK)
!
! CALL DGEMM('n','n',SS_DIM,SS_DIM,NLOC,1.0d0,KAL,SS_DIM,HEOF,ENKF_NOBSMAX,0.0d0,RALPHA,SS_DIM)
!
!! COMPUTE (I-KH)M AND ITS SINGULAR VALUE
! DO I=1, SS_DIM
! DO J=1, SS_DIM
! IF(I/=J) THEN
! RALPHA(I,J) = (-1.0_DP)*RALPHA(I,J)
! ELSE
! RALPHA(I,I) = 1.0_DP - RALPHA(I,I)
! ENDIF
! ENDDO
! ENDDO
!
! FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/I_KH.txt'
! OPEN(INORRK,FILE=FILENAME, FORM='FORMATTED')
! DO J=1, SS_DIM
! DO I=1, SS_DIM
! WRITE(INORRK,*) I,J,RALPHA(I,J)
! ENDDO
! ENDDO
! CLOSE(INORRK)
!
! FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/Amatr.dat'
! OPEN(INORRK,FILE=FILENAME, FORM='UNFORMATTED')
! DO J=1, SS_DIM
! READ(INORRK) (TRANS(I,J), I=1, SS_DIM)
! ENDDO
! CLOSE(INORRK)
! CALL DGEMM('n','n',SS_DIM,SS_DIM,SS_DIM,1.0d0,RALPHA,SS_DIM,TRANS,SS_DIM,0.0d0,RALPHA,SS_DIM)
!
! FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/M_KHM.txt'
! OPEN(INORRK,FILE=FILENAME, FORM='FORMATTED')
! DO J=1, SS_DIM
! DO I=1, SS_DIM
! WRITE(INORRK,*) I,J,RALPHA(I,J)
! ENDDO
! ENDDO
! CLOSE(INORRK)
!
! DEALLOCATE(STLOC,IPIV,IPIV2,KAL,HEOF,R,HTR,R2,PHT,RALPHA,TRANS2)
!! DEALLOCATE(BETA2,GAMMA,GAMMA2,WORK,WORK2,SEOF,STSD,HU,HUL,EVAL)
! DEALLOCATE(BETA2,GAMMA,GAMMA2,WORK,WORK2,HU,HUL,EVAL)
!! DEALLOCATE(STTEMP0,STTEMP1,STEOF,SDEOF,TRANS,RRKU,RRKV,RRKEL,RRKT,RRKS)
! ENDIF
!
! RETURN
! END SUBROUTINE RRK_RRK
! UTILITIES PROGRAMS
! SUBROUTINE RRK_ALLOC_VAR
!
! USE CONTROL
! USE ENKFVAL
! IMPLICIT NONE
!
! INTEGER STDIM
! INTEGER SS_DIM
!
! STDIM = 0
! IF(EL_ASSIM) STDIM = STDIM + MGL
! IF(UV_ASSIM) STDIM = STDIM + 2*NGL*KBM1
! IF(T_ASSIM) STDIM = STDIM + MGL*KBM1
! IF(S_ASSIM) STDIM = STDIM + MGL*KBM1
!
! SS_DIM = ENKF_NENS
!
!! ALLOCATE ARRYS
!
! IF(.not.ALLOCATED(STTEMP0)) ALLOCATE(STTEMP0(STDIM)) ;STTEMP0 = ZERO
! IF(.not.ALLOCATED(STTEMP1)) ALLOCATE(STTEMP1(STDIM)) ;STTEMP1 = ZERO
! IF(.not.ALLOCATED(STEOF)) ALLOCATE(STEOF(STDIM)) ;STEOF = ZERO
! IF(.not.ALLOCATED(SDEOF)) ALLOCATE(SDEOF(STDIM)) ;SDEOF = ZERO
! IF(.not.ALLOCATED(TRANS)) ALLOCATE(TRANS(SS_DIM,SS_DIM)) ;TRANS = ZERO
! IF(.not.ALLOCATED(RRKEL)) ALLOCATE(RRKEL(0:MGL)) ;RRKEL = ZERO
! IF(.not.ALLOCATED(RRKU)) ALLOCATE(RRKU(0:NGL,KB)) ;RRKU = ZERO
! IF(.not.ALLOCATED(RRKV)) ALLOCATE(RRKV(0:NGL,KB)) ;RRKV = ZERO
! IF(.not.ALLOCATED(RRKT)) ALLOCATE(RRKT(0:MGL,KB)) ;RRKT = ZERO
! IF(.not.ALLOCATED(RRKS)) ALLOCATE(RRKS(0:MGL,KB)) ;RRKS = ZERO
!
!! END ALLOCATION
!
! RETURN
! END SUBROUTINE RRK_ALLOC_VAR
!------------------------------------------------------------------------------|
! PERTURB THE BASE STATE IN THE EIGENVECTOR DIRECTIONS |
!------------------------------------------------------------------------------|
! SUBROUTINE PERTURB
!
!# if defined(WET_DRY)
! USE MOD_WD
!# endif
! USE LIMS
! USE ALL_VARS
! USE ENKFVAL
! USE CONTROL
!# if defined (MULTIPROCESSOR)
! USE MOD_PAR
!# endif
! IMPLICIT NONE
!
! INTEGER SS_DIM
! INTEGER STDIM
! INTEGER I_EOF
! INTEGER II,I,J
! INTEGER IDUMMY
! INTEGER I_START
! CHARACTER(LEN=80) IFILE
! CHARACTER(LEN=80) IFILE2
! CHARACTER(LEN=80) TEMP1FILE
! CHARACTER(LEN=80) TEMP2FILE
! CHARACTER(LEN=4) FEOF
!
! STDIM = 0
! IF(EL_ASSIM) STDIM = STDIM + MGL
! IF(UV_ASSIM) STDIM = STDIM + 2*NGL*KBM1
! IF(T_ASSIM) STDIM = STDIM + MGL*KBM1
! IF(S_ASSIM) STDIM = STDIM + MGL*KBM1
!
! SS_DIM = ENKF_NENS
!
! ALLOCATE(STTEMP0(STDIM)) ;STTEMP0 = ZERO
! ALLOCATE(STTEMP1(STDIM)) ;STTEMP1 = ZERO
! ALLOCATE(STEOF(STDIM)) ;STEOF = ZERO
! ALLOCATE(SDEOF(STDIM)) ;SDEOF = ZERO
! ALLOCATE(TRANS(SS_DIM,SS_DIM)) ;TRANS = ZERO
!
! IDUMMY = 0
! IF(EL_ASSIM) THEN
! DO I=1, MGL
! IDUMMY = IDUMMY + 1
! SDEOF(IDUMMY) = ELSD
! ENDDO
! ENDIF
! IF(UV_ASSIM) THEN
! DO I=1, 2*KBM1
! DO J=1, NGL
! IDUMMY = IDUMMY + 1
! SDEOF(IDUMMY)=USD*DSQRT(DBLE(KBM1))
! ENDDO
! ENDDO
! ENDIF
! IF(T_ASSIM) THEN
! DO I=1, KBM1
! DO J=1, MGL
! IDUMMY = IDUMMY + 1
! SDEOF(IDUMMY)=TSD
! ENDDO
! ENDDO
! ENDIF
! IF(S_ASSIM) THEN
! DO I=1, KBM1
! DO J=1, MGL
! IDUMMY = IDUMMY + 1
! SDEOF(IDUMMY)=SSD
! ENDDO
! ENDDO
! ENDIF
!
! print *, 'before rrk_rrk read restart, file name:',nc_start%fname ! should be marked
! CALL READRESTART(NC_START,REF_START_TIME)
! CALL GR2ST(0)
!
!
!! READ THE EOF AND PERTURB THE BASE STATE IN THIS DIRECTION
!
! DO I_EOF=1, SS_DIM
!! CALL READEOF(I_EOF,1)
! DO I=1, STDIM
! STEOF(I) = SEOF(I,I_EOF)
! ENDDO
!
!! PERTURB THE BASE STATE IN THE DIRECTION OF THE I_EOF'TH EOF AND STORE IT IN THE MP1FILE (FOR RESTART)
! DO I=1, STDIM
! STTEMP1(I) = STTEMP0(I) + STEOF(I)*SDEOF(I)*DBLE(RRK_PSIZE)
! ENDDO
! CALL ST2GR
!
!
! IF(SERIAL) THEN
! EL = RRKEL
! U = RRKU
! V = RRKV
! T1 = RRKT
! S1 = RRKS
! ELSE
!# if defined (MULTIPROCESSOR)
! DO I=1, M
! EL(I)=RRKEL(NGID(I))
! ENDDO
!
! DO I=1, N
! DO J=1, KB
! U(I,J)=RRKU(EGID(I),J)
! V(I,J)=RRKV(EGID(I),J)
! ENDDO
! ENDDO
!
! DO I=1, M
! DO J=1, KB
! T1(I,J)=RRKT(NGID(I),J)
! S1(I,J)=RRKS(NGID(I),J)
! ENDDO
! ENDDO
!# endif
! ENDIF
!
!
!! WRITE THE PERTURBED STATE IN TEMP1FILE
! print *, 'where to read to perturb:', nc_start%fname
! print *, 'where to write perturb:', NCF_ENKFRE%fname
!
! CALL WRITERESTART(NCF_ENKFRE,I_EOF)
!
! ENDDO
!
!
! RETURN
! END SUBROUTINE PERTURB
!------------------------------------------------------------------------------|
! CALCULATE THE LINEARIZED MODEL IN THE REDUCED SPACE |
!------------------------------------------------------------------------------|
! SUBROUTINE MREDUCED
!
! USE MOD_WD
! USE LIMS
! USE CONTROL
! USE ENKFVAL
! IMPLICIT NONE
!
! INTEGER SS_DIM
! INTEGER STDIM
! INTEGER I_EOF
! INTEGER I,J
! INTEGER IDUMMY
! INTEGER I_START
! CHARACTER(LEN=80) IFILE
! CHARACTER(LEN=80) IFILE2
! CHARACTER(LEN=80) TEMPFILE
! CHARACTER(LEN=4) FEOF
! REAL(DP) :: SUM0
!
! STDIM = 0
! IF(EL_ASSIM) STDIM = STDIM + MGL
! IF(UV_ASSIM) STDIM = STDIM + 2*NGL*KBM1
! IF(T_ASSIM) STDIM = STDIM + MGL*KBM1
! IF(S_ASSIM) STDIM = STDIM + MGL*KBM1
!
! SS_DIM = ENKF_NENS
!
! IDUMMY = 0
! IF(EL_ASSIM) THEN
! DO I=1, MGL
! IDUMMY = IDUMMY + 1
! SDEOF(IDUMMY) = ELSD
! ENDDO
! ENDIF
! IF(UV_ASSIM) THEN
! DO I =1, 2*KBM1
! DO J=1, NGL
! IDUMMY = IDUMMY + 1
! IF(RRK_OPTION == 0) THEN
! SDEOF(IDUMMY) = USD
! ELSE
! SDEOF(IDUMMY) = USD*DSQRT(DBLE(KBM1))
! ENDIF
! ENDDO
! ENDDO
! ENDIF
! IF(T_ASSIM) THEN
! DO I =1, KBM1
! DO J=1, MGL
! IDUMMY = IDUMMY + 1
! SDEOF(IDUMMY) = TSD
! ENDDO
! ENDDO
! ENDIF
! IF(S_ASSIM) THEN
! DO I =1, KBM1
! DO J=1, MGL
! IDUMMY = IDUMMY + 1
! SDEOF(IDUMMY) = SSD
! ENDDO
! ENDDO
! ENDIF
!
! CALL READRESTART(NC_START,REF_START_TIME)
! CALL GR2ST(0)
!
!
!! READ THE FORECAST WHICH IS PROPAGATED FROM THE PERTURBED STATE
!
! CALL NC_OPEN(NCF_ENKFfct)
!
! DO I_EOF =1, SS_DIM
!! READ THE EACH PERTURBED STATE AT THE END OF LINEARUZATION TIME STEP
!
! CALL READRESTART(NCF_ENKFfct,I_EOF)
!
! CALL GR2ST(1)
!
!! PROJECT EVOLVED PERTURBATION ONTO EOFs
! DO I=1, SS_DIM
! DO J=1, STDIM
! STEOF(J) = SEOF(J,I)
! ENDDO
! SUM0 = 0.0_DP
! DO J=1, STDIM
! SUM0 = SUM0 + STEOF(J)/SDEOF(J)*(STTEMP1(J)-STTEMP0(J))
! ENDDO
! TRANS(I,I_EOF) = SUM0/DBLE(RRK_PSIZE)
!
! ENDDO
!! EDN OF LOOP OVER EACH EOF
! ENDDO
!
! END SUBROUTINE MREDUCED
!=====================================================================================/
! CONVERT THE STATE FROM THE GRID SPACE TO THE STATE VECTOR SPACE /
!=====================================================================================/
! SUBROUTINE GR2ST(OPT)
!
! USE LIMS
! USE ALL_VARS
! USE ENKFVAL
! IMPLICIT NONE
!
! INTEGER IDUMMY, I, J
! INTEGER OPT
!
! IDUMMY=0
!
! IF(EL_ASSIM) THEN
! DO I=1, MGL
! IDUMMY = IDUMMY + 1
! IF (OPT == 0) THEN
! STTEMP0(IDUMMY) = EL(I)
! ELSE
! STTEMP1(IDUMMY) = EL(I)
! ENDIF
! ENDDO
! ENDIF
!
! IF(UV_ASSIM) THEN
! DO I=1, KBM1
! DO J=1, NGL
! IDUMMY = IDUMMY + 1
! IF (OPT == 0) THEN
! STTEMP0(IDUMMY) = U(J,I)
!
! ELSE
! STTEMP1(IDUMMY) = U(J,I)
! ENDIF
! ENDDO
! ENDDO
! DO I=1, KBM1
! DO J=1, NGL
! IDUMMY = IDUMMY + 1
! IF (OPT == 0) THEN
! STTEMP0(IDUMMY) = V(J,I)
! ELSE
! STTEMP1(IDUMMY) = V(J,I)
! ENDIF
! ENDDO
! ENDDO
! ENDIF
!
! IF(T_ASSIM) THEN
! DO I=1, KBM1
! DO J=1, MGL
! IDUMMY = IDUMMY + 1
! IF (OPT == 0) THEN
! STTEMP0(IDUMMY) = T(J,I)
! ELSE
! STTEMP1(IDUMMY) = T(J,I)
! ENDIF
! ENDDO
! ENDDO
! ENDIF
!
! IF(S_ASSIM) THEN
! DO I=1, KBM1
! DO J=1, MGL
! IDUMMY = IDUMMY + 1
! IF (OPT == 0) THEN
! STTEMP0(IDUMMY) = S(J,I)
! ELSE
! STTEMP1(IDUMMY) = S(J,I)
! ENDIF
! ENDDO
! ENDDO
! ENDIF
!
! RETURN
! END SUBROUTINE GR2ST
!=====================================================================================/
! CONVERT THE STATE FROM THE STATE VECTOR SPACE TO THE GRID SPACE /
!=====================================================================================/
! SUBROUTINE ST2GR
!
! USE LIMS
! USE ALL_VARS
! USE ENKFVAL
! IMPLICIT NONE
!
! INTEGER IDUMMY, I, J
!
! IDUMMY=0
! IF(EL_ASSIM) THEN
! DO I=1, MGL
! IDUMMY = IDUMMY + 1
! RRKEL(I) = STTEMP1(IDUMMY)
! ENDDO
! ENDIF
! IF(UV_ASSIM) THEN
! DO I=1, KBM1
! DO J=1, NGL
! IDUMMY = IDUMMY + 1
! RRKU(J,I) = STTEMP1(IDUMMY)
! ENDDO
! ENDDO
! DO I=1, KBM1
! DO J=1, NGL
! IDUMMY = IDUMMY + 1
! RRKV(J,I) = STTEMP1(IDUMMY)
! ENDDO
! ENDDO
! ENDIF
! IF(T_ASSIM) THEN
! DO I=1, KBM1
! DO J=1, MGL
! IDUMMY = IDUMMY + 1
! RRKT(J,I) = STTEMP1(IDUMMY)
! ENDDO
! ENDDO
! ENDIF
! IF(S_ASSIM) THEN
! DO I=1, KBM1
! DO J=1, MGL
! IDUMMY = IDUMMY + 1
! RRKS(J,I) = STTEMP1(IDUMMY)
! ENDDO
! ENDDO
! ENDIF
!
! RETURN
! END SUBROUTINE ST2GR
!=====================================================================================/
! DETERMINE IF NODES/ELEMENTS ARE WET OR DRY /
!=====================================================================================/
! SUBROUTINE WET_JUDGE_EL
!
! USE MOD_PREC
! USE ALL_VARS
!# if defined (MULTIPROCESSOR)
! USE MOD_PAR
!# endif
!# if defined(WET_DRY)
! USE MOD_WD
!# endif
! IMPLICIT NONE
! REAL(SP) :: DTMP
! INTEGER :: ITA_TEMP
! INTEGER :: I,IL,IA,IB,K1,K2,K3,K4,K5,K6
!
!# if defined(WET_DRY)
!
!!
!!--Determine If Node Points Are Wet/Dry Based on Depth Threshold---------------!
!!
! ISWETN = 1
! DO I = 1, M
! DTMP = H(I) + EL(I)
! IF((DTMP - MIN_DEPTH) < 1.0E-5_SP) ISWETN(I) = 0
! END DO
!
!!
!!--Determine if Cells are Wet/Dry Based on Depth Threshold---------------------!
!!
! ISWETC = 1
! DO I = 1, N
! DTMP = MAX(EL(NV(I,1)),EL(NV(I,2)),EL(NV(I,3))) + &
! MIN( H(NV(I,1)), H(NV(I,2)), H(NV(I,3)))
! IF((DTMP - MIN_DEPTH) < 1.0E-5_SP) ISWETC(I) = 0
! END DO
!
!!
!!--A Secondary Condition for Nodal Dryness-(All Elements Around Node Are Dry)--!
!!
! DO I = 1, M
! IF(SUM(ISWETC(NBVE(I,1:NTVE(I)))) == 0) ISWETN(I) = 0
! END DO
!
!!
!!--Adjust Water Surface So It Does Not Go Below Minimum Depth------------------!
!!
! EL = MAX(EL,-H + MIN_DEPTH)
!
!!
!!--Recompute Element Based Depths----------------------------------------------!
!!
! DO I = 1, N
! EL1(I) = ONE_THIRD*(EL(NV(I,1))+EL(NV(I,2))+EL(NV(I,3)))
! END DO
!
!!
!!--Extend Element/Node Based Wet/Dry Flags to Domain Halo----------------------!
!!
!# if defined (MULTIPROCESSOR)
! IF(PAR)THEN
!
! CALL AEXCHANGE(EC,MYID,NPROCS,ISWETC)
! CALL AEXCHANGE(NC,MYID,NPROCS,ISWETN)
!
! END IF
!
!# endif
!# endif
!
! RETURN
!
! END SUBROUTINE WET_JUDGE_EL
!!$!==============================================================================|
!!$! DUMP WET/DRY FLAG DATA FOR RESTART |
!!$!==============================================================================|
!!$
!!$ SUBROUTINE WD_DUMP_EL(INF,I_START)
!!$
!!$!------------------------------------------------------------------------------|
!!$
!!$ USE ALL_VARS
!!$# if defined (MULTIPROCESSOR)
!!$ USE MOD_PAR
!!$# endif
!!$# if defined(WET_DRY)
!!$ USE MOD_WD
!!$# endif
!!$
!!$ IMPLICIT NONE
!!$ INTEGER, ALLOCATABLE,DIMENSION(:) :: NTEMP1,NTEMP2
!!$ INTEGER I, INF
!!$ INTEGER I_START
!!$!==============================================================================|
!!$
!!$# if defined(WET_DRY)
!!$
!!$ IF(MSR)THEN
!!$ REWIND(INF)
!!$ WRITE(INF,*) I_START
!!$ WRITE(INF,*) NGL,MGL
!!$ END IF
!!$
!!$ IF(SERIAL)THEN
!!$ WRITE(INF,*) (ISWETC(I), I=1,N)
!!$ WRITE(INF,*) (ISWETN(I), I=1,M)
!!$ ELSE
!!$ ALLOCATE(NTEMP1(NGL),NTEMP2(MGL))
!!$# if defined (MULTIPROCESSOR)
!!$ CALL IGATHER(LBOUND(ISWETC,1),UBOUND(ISWETC,1),N,NGL,1,MYID,NPROCS,EMAP,ISWETC,NTEMP1)
!!$ CALL IGATHER(LBOUND(ISWETN,1),UBOUND(ISWETN,1),M,MGL,1,MYID,NPROCS,NMAP,ISWETN,NTEMP2)
!!$ IF(MSR)THEN
!!$ WRITE(INF,*) (NTEMP1(I), I=1,NGL)
!!$ WRITE(INF,*) (NTEMP2(I), I=1,MGL)
!!$ END IF
!!$ DEALLOCATE(NTEMP1,NTEMP2)
!!$# endif
!!$ END IF
!!$
!!$ CLOSE(INF)
!!$
!!$# endif
!!$
!!$ RETURN
!!$ END SUBROUTINE WD_DUMP_EL
SUBROUTINE READOBS(STLOC,NLOC)
USE LIMS
USE CONTROL
IMPLICIT NONE
INTEGER :: NUM = 0
INTEGER :: NLOC
INTEGER :: SWITCH = 0
INTEGER :: J,K
INTEGER :: IDUMMY
INTEGER :: TMP
CHARACTER(LEN=80) FILENAME
CHARACTER(LEN=24) HEADINFO
INTEGER STLOC(ENKF_NOBSMAX)
INTEGER LAY(ENKF_NOBSMAX)
INTEGER IERR
FILENAME = TRIM(INPUT_DIR)//"/"//trim(casename)//"_assim_enkf.nml"
OPEN(73,FILE=TRIM(FILENAME),FORM='FORMATTED')
REWIND(73)
NLOC = 0
IDUMMY = 0
STLOC = ZERO
DO WHILE ( .TRUE. )
READ(73,'(A24)',IOSTAT=IERR) HEADINFO
IF(IERR /= 0) exit
IF(SWITCH/=1) THEN
IF(HEADINFO=='!===READ IN OBSERVATION') SWITCH = 1
CYCLE
ENDIF
IF(TRIM(HEADINFO)=='!EL') THEN
IF(EL_OBS) THEN
READ(73,*) NUM
NLOC = NLOC + NUM
IF(NLOC>ENKF_NOBSMAX) THEN
WRITE(IPT,*) 'not enough storage for observations:', 'Nloc=', Nloc, 'Nobsmax=', ENKF_NOBSMAX
CALL PSTOP
ENDIF
READ(73,*) (STLOC(K), K=1,NLOC)
ENDIF
IF(EL_ASSIM) THEN
IDUMMY = IDUMMY + MGL
ENDIF
ENDIF
IF(TRIM(HEADINFO)=='!UV') THEN
IF(UV_OBS) THEN
READ(73,*) NUM
NLOC = NLOC + NUM
IF(NLOC+NUM>ENKF_NOBSMAX) THEN
WRITE(IPT,*) 'not enough storage for observations:', 'Nloc=', Nloc+num, 'Nobsmax=', ENKF_NOBSMAX
CALL PSTOP
ENDIF
READ(73,*) (STLOC(K), K=NLOC-NUM+1,NLOC)
READ(73,*) (LAY(K), K=NLOC-NUM+1,NLOC)
DO K=NLOC-NUM+1, NLOC
STLOC(K)=STLOC(K)+IDUMMY+NGL*(LAY(K)-1)
ENDDO
IDUMMY = IDUMMY + NGL*KBM1
NLOC = NLOC + NUM
DO K=NLOC-NUM+1, NLOC
STLOC(K)=STLOC(K-NUM)+NGL*KBM1+NGL*(LAY(K-NUM)-1)
ENDDO
ENDIF
IF(UV_ASSIM) THEN
IDUMMY = IDUMMY + NGL*KBM1
ENDIF
ENDIF
IF(TRIM(HEADINFO)=='!T') THEN
IF(T_OBS) THEN
READ(73,*) NUM
NLOC = NLOC + NUM
IF(NLOC>ENKF_NOBSMAX) THEN
WRITE(IPT,*) 'not enough storage for observations:', 'Nloc=', Nloc, 'Nobsmax=', ENKF_NOBSMAX
CALL PSTOP
ENDIF
READ(73,*) (STLOC(K), K=NLOC-NUM+1,NLOC)
READ(73,*) (LAY(K), K=NLOC-NUM+1,NLOC)
DO K=NLOC-NUM+1, NLOC
STLOC(K)=STLOC(K)+IDUMMY+MGL*(LAY(K)-1)
ENDDO
ENDIF
IF(T_ASSIM) THEN
IDUMMY = IDUMMY + MGL*KBM1
ENDIF
ENDIF
IF(TRIM(HEADINFO)=='!S') THEN
IF(S_OBS) THEN
READ(73,*) NUM
PRINT *, 'SALINITY READING OBSERVATION',NUM
NLOC = NLOC + NUM
IF(NLOC>ENKF_NOBSMAX) THEN
WRITE(IPT,*) 'not enough storage for observations:', 'Nloc=', Nloc, 'Nobsmax=', ENKF_NOBSMAX
CALL PSTOP
ENDIF
READ(73,*) (STLOC(K),K=NLOC-NUM+1,NLOC)
READ(73,*) (LAY(K), K=NLOC-NUM+1,NLOC)
DO K=NLOC-NUM+1, NLOC
STLOC(K)=STLOC(K)+IDUMMY+MGL*(LAY(K)-1)
ENDDO
ENDIF
IF(S_ASSIM) THEN
IDUMMY = IDUMMY + MGL*KBM1
ENDIF
ENDIF
ENDDO
CLOSE(73)
RETURN
END SUBROUTINE READOBS
!-------------old--------------------------------
SUBROUTINE ENKF_READRESTART_first(NCF)
USE MOD_INPUT
USE MOD_STARTUP
IMPLICIT NONE
TYPE(NCFILE), POINTER :: NCF, NCF_tmP
INTEGER :: STATE
NCF_tmP=>NC_START
IF(.not. associated(NCF)) CALL FATAL_ERROR&
& ("MOD_ENKF: READRESTART: THE FILE POINTER PASSED IS NOT ASSOCIATED?")
NCF%FTIME%PREV_STKCNT = 1
NC_START=>NCF
print *, 'read ssh enkf'
CALL READ_SSH
print *, 'finish read ssh enkf'
IF(WETTING_DRYING_ON) CALL READ_WETDRY
CALL READ_UV
CALL READ_TURB
CALL READ_TS
NC_START=>NCF_tmP
END SUBROUTINE ENKF_READRESTART_first
!-------------------new-------------------------------------------------------
!-------------old--------------------------------
SUBROUTINE ENKF_READRESTART_last(NCF)
USE MOD_INPUT
USE MOD_STARTUP
IMPLICIT NONE
TYPE(NCFILE), POINTER :: NCF, NCF_tmP
INTEGER :: STATE
TYPE(NCDIM), POINTER :: DIM
TYPE(NCVAR), POINTER :: VAR
LOGICAL :: FOUND
INTEGER :: NTIMES
! GET THE TIME DIMENSION
DIM => FIND_UNLIMITED(NCF,FOUND)
IF(.not. FOUND) CALL FATAL_ERROR &
& ("IN ENKF_READ_FIELD FILE",&
& "FILE NAME: "//TRIM(NCF%FNAME),&
&"COULD NOT FIND THE UNLIMITED DIMENSION")
NTIMES = DIM%DIM
NCF_tmP=>NC_START
IF(.not. associated(NCF)) CALL FATAL_ERROR&
& ("MOD_ENKF: READRESTART: THE FILE POINTER PASSED IS NOT ASSOCIATED?")
NCF%FTIME%PREV_STKCNT = NTIMES
NC_START=>NCF
print *, 'read ssh enkf'
CALL READ_SSH
print *, 'finish read ssh enkf'
IF(WETTING_DRYING_ON) CALL READ_WETDRY
CALL READ_UV
CALL READ_TURB
CALL READ_TS
NC_START=>NCF_tmP
END SUBROUTINE ENKF_READRESTART_last
!-------------------new-------------------------------------------------------
SUBROUTINE ENKF_READRESTART(NCF)
USE CONTROL
USE ENKFVAL
USE MOD_NCTOOLS
USE MOD_INPUT
USE MOD_STARTUP
IMPLICIT NONE
TYPE(NCFILE), POINTER :: NCF, NCF_tmP
INTEGER :: STATE
TYPE(TIME) :: tTEST
LOGICAL :: FOUND
INTEGER :: NTIMES
TYPE(NCDIM), POINTER :: DIM
TYPE(NCVAR), POINTER :: VAR
INTEGER :: nrange,IRANGE(2),CNT, STATUS
integer :: istkcnt
character(len=120) :: pathnfile
! GET THE TIME DIMENSION
DIM => FIND_UNLIMITED(NCF,FOUND)
IF(.not. FOUND) CALL FATAL_ERROR &
& ("IN ENKF_READ_FIELD FILE",&
& "FILE NAME: "//TRIM(NCF%FNAME),&
&"COULD NOT FIND THE UNLIMITED DIMENSION")
NTIMES = DIM%DIM
! TEST THE FILE START AND END
TTEST = get_file_time(NCF,1)
IF(INTTIME < TTEST) CALL FATAL_ERROR &
&("ENKF_READ_FIELD: INTTIME IS BEFORE FIRST TIME IN FILE.",&
& "FILE NAME:"//TRIM(NCF%FNAME))
TTEST = get_file_time(NCF,NTIMES)
IF(INTTIME > TTEST) CALL FATAL_ERROR &
&("ENKD_READ_FIELD: INTTIME IS AFTER LAST TIME IN FILE.",&
& "FILE NAME:"//TRIM(NCF%FNAME))
! GET THE START AND END INDEX
IRANGE = 0
istkcnt=0
CNT = 0
DO WHILE (product(irange)==0)
CNT = CNT +1
IF (CNT > NTIMES) THEN
CALL PRINT_TIME(INTTIME,IPT,"INTTIME")
CALL FATAL_ERROR&
&("ENKF_READ_FIELD: CAN NOT FIND INTTIME IN THE THIS FILE!", &
& "FILE NAME:"//TRIM(NCF%FNAME))
END IF
TTEST = get_file_time(NCF,CNT)
IF (INTTIME == TTEST) THEN
IRANGE(1) = CNT
istkcnt=cnt
END IF
IF (inttime == TTEST) THEN
IRANGE(2) = CNT
END IF
END DO
nrange = irange(2)-irange(1)+1
!------------------FINISH FIND CNT---------------------
NCF_tmP=>NC_START
IF(.not. associated(NCF)) CALL FATAL_ERROR&
& ("MOD_ENKF: READRESTART: THE FILE POINTER PASSED IS NOT ASSOCIATED?")
NCF%FTIME%PREV_STKCNT = CNT
NC_START=>NCF
print *, 'read ssh enkf'
CALL READ_SSH
print *, 'finish read ssh enkf'
IF(WETTING_DRYING_ON) CALL READ_WETDRY
CALL READ_UV
CALL READ_TURB
CALL READ_TS
NC_START=>NCF_tmP
END SUBROUTINE ENKF_READRESTART
! SUBROUTINE ENKF_READRESTART(NCF)
! USE MOD_INPUT
! USE MOD_STARTUP_mimic
! IMPLICIT NONE
! TYPE(NCFILE), POINTER :: NCF
! INTEGER :: STATE
! IF(.not. associated(NCF)) CALL FATAL_ERROR&
! & ("MOD_ENKF: READRESTART: THE FILE POINTER PASSED IS NOT ASSOCIATED?")
! NCF%FTIME%PREV_STKCNT = 1
! CALL READ_SSH1(NCF)
! IF(WETTING_DRYING_ON) CALL READ_WETDRY1(NCF)
! CALL READ_UV1(NCF)
! CALL READ_TURB1(NCF)
! CALL READ_TS1(NCF)
! END SUBROUTINE ENKF_READRESTART
! SUBROUTINE READRESTART_STATE(NCF,STATE)
! USE MOD_INPUT
! USE MOD_STARTUP_mimic
! IMPLICIT NONE
! TYPE(NCFILE), POINTER :: NCF
! INTEGER :: STATE
!
! IF(.not. associated(NCF)) CALL FATAL_ERROR&
! & ("MOD_ENKF: READRESTART: THE FILE POINTER PASSED IS NOT ASSOCIATED?")
!
! NCF%FTIME%PREV_STKCNT = STATE
!
!
!
! CALL READ_SSH1(NCF)
!
! IF(WETTING_DRYING_ON) CALL READ_WETDRY1(NCF)
!
! CALL READ_UV1(NCF)
!
! CALL READ_TURB1(NCF)
!
! CALL READ_TS1(NCF)
!
! END SUBROUTINE READRESTART_STATE
!
! SUBROUTINE READRESTART_TIME(NCF,NOW)
! USE MOD_INPUT
! USE MOD_STARTUP_MIMIC
! IMPLICIT NONE
! TYPE(NCFILE), POINTER :: NCF
! TYPE(TIME) :: NOW,TTEST
! INTEGER :: STATE,NSTATE
! TYPE(NCDIM), POINTER :: DIM
! LOGICAL FOUND
!
! IF(.not. associated(NCF)) CALL FATAL_ERROR&
! & ("MOD_ENKF: READRESTART: THE FILE POINTER PASSED IS NOT ASSOCIATED?")
!
! ! GET THE TIME DIMENSION
! DIM => FIND_UNLIMITED(NCF,FOUND)
! IF(.not. FOUND) CALL FATAL_ERROR &
! & ("IN RESTART FILE",&
! & "FILE NAME: "//TRIM(NCF%FNAME),&
! &"COULD NOT FIND THE UNLIMITED DIMENSION")
! NSTATE = DIM%DIM
! FOUND = .false.
! DO STATE=1,NSTATE
!
! TTEST = get_file_time(NCF,STATE)
! IF(TTEST == NOW) THEN
! FOUND = .true.
! exit
! END IF
!
! END DO
!
! IF(FOUND)THEN
! CALL READRESTART_STATE(NCF,STATE)
! ELSE
! CALL PRINT_TIME(NOW,IPT,"TIME REQUSTED")
! CALL FATAL_ERROR&
! &("COULD NOT FIND CORRECT TIME IN RESTART FILE",&
! &"FILE NAME: "//TRIM(NCF%FNAME))
!
! END IF
!
! END SUBROUTINE READRESTART_TIME
! SUBROUTINE WRITERESTART(NCF,STATE)
! IMPLICIT NONE
! TYPE(NCFILE), POINTER :: NCF
! INTEGER :: STATE
!
! IF(.not. associated(NCF)) CALL FATAL_ERROR&
! & ("MOD_ENKF: WRITERESTART: THE FILE POINTER PASSED IS NOT ASSOCIATED?")
!
! NCF%FTIME%NEXT_STKCNT = state
! CALL UPDATE_IODATA(NCF,IntTime)
! CALL NC_WRITE_FILE(NCF,LOCAL_DISK)
!
! END SUBROUTINE WRITERESTART
!-----------------old-----------------------
! SUBROUTINE ENKF_WRITERESTART(NCF)
! IMPLICIT NONE
! TYPE(NCFILE), POINTER :: NCF
! INTEGER :: STATE
!
! IF(.not. associated(NCF)) CALL FATAL_ERROR&
! & ("MOD_ENKF: WRITERESTART: THE FILE POINTER PASSED IS NOT ASSOCIATED?")
!
!! NCF%FTIME%NEXT_STKCNT = state
! NCF%FTIME%NEXT_STKCNT = 1
! CALL UPDATE_IODATA(NCF,IntTime)
! CALL NC_WRITE_FILE(NCF,LOCAL_DISK)
!
! END SUBROUTINE ENKF_WRITERESTART
!-------------new----------------------------------
SUBROUTINE ENKF_WRITERESTART(NCF)
IMPLICIT NONE
TYPE(NCFILE), POINTER :: NCF
INTEGER :: STATE
TYPE(NCFTIME), POINTER :: FTM
IF(.not. associated(NCF)) CALL FATAL_ERROR&
& ("MOD_ENKF: WRITERESTART: THE FILE POINTER PASSED IS NOT ASSOCIATED?")
! NCF%FTIME%NEXT_STKCNT = state
! NCF%FTIME%NEXT_STKCNT = 1
! CALL UPDATE_IODATA(NCF,IntTime)
! CALL NC_WRITE_FILE(NCF,LOCAL_DISK)
FTM => NCF%FTIME
CALL UPDATE_IODATA(NCF,IntTime)
CALL NC_WRITE_FILE(NCF)
IF (ENKF_memberCONTR == ENKF_NENS) then
! ONCE THE FILE IS WRITEN INCRIMENT THE FILE OBJECT TIME
FTM%PREV_IO = IntTime
FTM%NEXT_IO = FTM%NEXT_IO + ENKF_INTERVAL
! INCRIMENT THE STACK COUNT
FTM%PREV_STKCNT = FTM%NEXT_STKCNT
FTM%NEXT_STKCNT = FTM%NEXT_STKCNT + 1
END IF
END SUBROUTINE ENKF_WRITERESTART
!=======================================================|
SUBROUTINE ENKF_SETUP_MEAN_NC(NCF)
IMPLICIT NONE
INTEGER :: STKCNT
TYPE(NCFILE), POINTER :: NCF
TYPE(NCFTIME), POINTER :: FTM
IF(.not. associated(NCF)) CALL FATAL_ERROR&
& ("MOD_ENKF: WRITERESTART: THE FILE POINTER PASSED IS NOT ASSOCIATED?")
if (ARCHIVE_TIMES==0) then
enkf_MEAN_ncfout => COPY_FILE(NCF)
enkf_MEAN_ncfout%FNAME=TRIM(OUTPUT_DIR)//"enkftemp/enkf_forecast_mean.nc"
enkf_MEAN_ncfout%FTIME%MAX_STKCNT =0
FTM => enkf_mean_ncfout%FTIME
CALL NC_WRITE_FILE(enkf_mean_ncfout,stkcnt=0) ! because this program is called after the call archive , it FTM%NEXT_IO already became 1, then I don't need below
! IF (ENKF_memberCONTR == ENKF_NENS) then
! ! ONCE THE FILE IS WRITEN INCRIMENT THE FILE OBJECT TIME
!
!
! FTM%PREV_IO = IntTime
! FTM%NEXT_IO = FTM%NEXT_IO + FTM%INTERVAL
!
! ! INCRIMENT THE STACK COUNT
! FTM%PREV_STKCNT = FTM%NEXT_STKCNT
! FTM%NEXT_STKCNT = FTM%NEXT_STKCNT + 1
! END IF
WRITE(IPT,*) "END DUMP mean DATA kf_times =0"
enkf_MEAN_ncfout => COPY_FILE(NCF)
enkf_MEAN_ncfout%FNAME=TRIM(OUTPUT_DIR)//"enkftemp/enkf_analysis_mean.nc"
enkf_MEAN_ncfout%FTIME%MAX_STKCNT =0
FTM => enkf_mean_ncfout%FTIME
CALL NC_WRITE_FILE(enkf_mean_ncfout,stkcnt=0) ! because this program is called after the call archive , it FTM%NEXT_IO already became 1, then I don't need below
! IF (ENKF_memberCONTR == ENKF_NENS) then
! ! ONCE THE FILE IS WRITEN INCRIMENT THE FILE OBJECT TIME
!
!
! FTM%PREV_IO = IntTime
! FTM%NEXT_IO = FTM%NEXT_IO + FTM%INTERVAL
!
! ! INCRIMENT THE STACK COUNT
! FTM%PREV_STKCNT = FTM%NEXT_STKCNT
! FTM%NEXT_STKCNT = FTM%NEXT_STKCNT + 1
! END IF
WRITE(IPT,*) "END DUMP mean DATA kf_times =0"
! acturall I didn't use call below
! else
! FTM => enkf_mean_ncfout%FTIME
! CALL UPDATE_IODATA(enkf_mean_ncfout,IntTime)
!
! CALL NC_WRITE_FILE(enkf_mean_ncfout)
! IF (ENKF_memberCONTR == ENKF_NENS) then
! ! ONCE THE FILE IS WRITEN INCRIMENT THE FILE OBJECT TIME
! FTM%PREV_IO = IntTime
! FTM%NEXT_IO = FTM%NEXT_IO + FTM%INTERVAL
!
! ! INCRIMENT THE STACK COUNT
! FTM%PREV_STKCNT = FTM%NEXT_STKCNT
! FTM%NEXT_STKCNT = FTM%NEXT_STKCNT + 1
! END IF
! WRITE(IPT,*) "END DUMP mean DATA kftimes > 0, FTM%NEXT_STKCNT",FTM%NEXT_STKCNT
end if
END SUBROUTINE ENKF_SETUP_MEAN_NC
!=======================================================|
!=======================================================|
SUBROUTINE ENKF_SETUP_NC(NCF)
IMPLICIT NONE
INTEGER :: STKCNT
TYPE(NCFILE), POINTER :: NCF
TYPE(NCFTIME), POINTER :: FTM
IF(.not. associated(NCF)) CALL FATAL_ERROR&
& ("MOD_ENKF: WRITERESTART: THE FILE POINTER PASSED IS NOT ASSOCIATED?")
if (ARCHIVE_TIMES==0) then
enkf_ncfout => COPY_FILE(NCF)
enkf_ncfout%FNAME=TRIM(OUTPUT_DIR)//"enkftemp/enkf_output_"//enkf_num//".nc"
enkf_ncfout%FTIME%MAX_STKCNT =0
FTM => enkf_ncfout%FTIME
! CALL UPDATE_IODATA(enkf_ncfout,IntTime)
CALL NC_WRITE_FILE(enkf_ncfout,stkcnt=0) ! because this program is called after the call archive , it FTM%NEXT_IO already became 1, then I don't need below
IF (ENKF_memberCONTR == ENKF_NENS) then
! ! ONCE THE FILE IS WRITEN INCRIMENT THE FILE OBJECT TIME
!
!
! FTM%PREV_IO = IntTime
! FTM%NEXT_IO = FTM%NEXT_IO + FTM%INTERVAL
!
! ! INCRIMENT THE STACK COUNT
! FTM%PREV_STKCNT = FTM%NEXT_STKCNT
! FTM%NEXT_STKCNT = FTM%NEXT_STKCNT + 1
END IF
WRITE(IPT,*) "END DUMP_DATA kf_times =0"
else
enkf_ncfout%FNAME=TRIM(OUTPUT_DIR)//"enkftemp/enkf_output_"//enkf_num//".nc"
FTM => enkf_ncfout%FTIME
CALL UPDATE_IODATA(enkf_ncfout,IntTime)
CALL NC_WRITE_FILE(enkf_ncfout)
IF (ENKF_memberCONTR == ENKF_NENS) then
! ONCE THE FILE IS WRITEN INCRIMENT THE FILE OBJECT TIME
FTM%PREV_IO = IntTime
FTM%NEXT_IO = FTM%NEXT_IO + FTM%INTERVAL
! INCRIMENT THE STACK COUNT
FTM%PREV_STKCNT = FTM%NEXT_STKCNT
FTM%NEXT_STKCNT = FTM%NEXT_STKCNT + 1
END IF
WRITE(IPT,*) "END DUMP_DATA kftimes > 0, FTM%NEXT_STKCNT",FTM%NEXT_STKCNT
! call NC_WRITE_FILE(enkf_ncfout,stkcnt=STKCNT)
end if
END SUBROUTINE ENKF_SETUP_NC
!=======================================================|
SUBROUTINE ENKF_RESTART
USE MOD_SET_TIME
IMPLICIT NONE
if (kf_times>1) then
write(enkf_num,'(i2.2)') ENKF_memberCONTR
NCF_ENKFre => NEW_FILE()
NCF_ENKFre%FNAME=TRIM(OUTPUT_DIR)//"enkftemp/enkf_forecast_"//enkf_num//".nc"
Call NC_OPEN(NCF_ENKFre)
CALL NC_LOAD(NCF_ENKFre)
CALL ENKF_READRESTART_last(NCF_ENKFre)
CALL SET_LAST_FILE_TIME(NCF_ENKFre,StartTime, IINT)
IntTime = starttime
ExtTime = starttime
CALL KILL_FILE(NCF_ENKFre)
else
write(enkf_num,'(i2.2)') ENKF_memberCONTR
NCF_ENKFre => NEW_FILE()
NCF_ENKFre%FNAME=TRIM(INPUT_DIR)//"enkf_initial_"//enkf_num//".nc"
Call NC_OPEN(NCF_ENKFre)
CALL NC_LOAD(NCF_ENKFre)
CALL ENKF_READRESTART_last(NCF_ENKFre)
CALL SET_LAST_FILE_TIME(NCF_ENKFre,StartTime, IINT)
starttime=enkf_start_time
IntTime = enkf_start_time
ExtTime = enkf_start_time
iint=ekint_start
CALL KILL_FILE(NCF_ENKFre)
end if
END SUBROUTINE ENKF_RESTART
!======================================================|
FUNCTION SETUP_RESTART(FNAME) RESULT(NCF)
USE MOD_NCDIO
IMPLICIT NONE
CHARACTER(LEN=*) :: FNAME
TYPE(GRID), SAVE :: MYGRID
TYPE(NCFILE), POINTER :: NCF
TYPE(NCFILE), POINTER :: NCF2
IF(DBG_SET(DBG_SBR)) WRITE(IPT,*) "START SETUP_RSTFILE"
IF(DBG_SET(DBG_LOG)) THEN
write(IPT,*)"!--------------------------------------------------"
write(IPT,*)"! SETTING UP RESTART FILE OUTPUT..."
END IF
NCF => NEW_FILE()
NCF%FNAME=trim(fname)
print *, 'running setup_restart:',NCF%FNAME ! should be mark
NCF%FTIME=>NEW_FTIME()
CALL SET_FVCOM_GRID(MYGRID)
CALL DEFINE_DIMENSIONS(MYGRID)
NCF2 => GRID_FILE_OBJECT(MYGRID)
NCF => ADD(NCF,NCF2)
!!$ NCF => ADD(NCF,GRID_FILE_OBJECT(MYGRID) )
NCF2 => TIME_FILE_OBJECT()
NCF => ADD(NCF,NCF2)
!!$ NCF => ADD(NCF,TIME_FILE_OBJECT() )
NCF2 => ZETA_FILE_OBJECT()
NCF => ADD(NCF,NCF2)
!!$ NCF => ADD(NCF,ZETA_FILE_OBJECT() )
NCF2 => FILE_DATE_OBJECT()
NCF => ADD(NCF,NCF2)
!!$ NCF => ADD(NCF,FILE_DATE_OBJECT() )
NCF2 => VELOCITY_FILE_OBJECT()
NCF => ADD(NCF,NCF2)
!!$ NCF => ADD(NCF,VELOCITY_FILE_OBJECT() )
NCF2 => AVERAGE_VEL_FILE_OBJECT()
NCF => ADD(NCF,NCF2)
!!$ NCF => ADD(NCF,AVERAGE_VEL_FILE_OBJECT() )
NCF2 => VERTICAL_VEL_FILE_OBJECT()
NCF => ADD(NCF,NCF2)
!!$ NCF => ADD(NCF,VERTICAL_VEL_FILE_OBJECT() )
NCF2 => TURBULENCE_FILE_OBJECT()
NCF => ADD(NCF,NCF2)
!!$ NCF => ADD(NCF,TURBULENCE_FILE_OBJECT() )
NCF2 => SALT_TEMP_FILE_OBJECT()
NCF => ADD(NCF,NCF2)
!!$ NCF => ADD(NCF,SALT_TEMP_FILE_OBJECT() )
NCF2 => DENSITY_FILE_OBJECT()
NCF => ADD(NCF,NCF2)
!!$ NCF => ADD(NCF,DENSITY_FILE_OBJECT() )
NCF2 => RESTART_EXTRAS_FILE_OBJECT()
NCF => ADD(NCF,NCF2)
!!$ NCF => ADD(NCF,RESTART_EXTRAS_FILE_OBJECT() )
IF(WETTING_DRYING_ON) THEN
NCF2 => WET_DRY_FILE_OBJECT()
NCF => ADD(NCF, NCF2)
!!$ NCF => ADD(NCF, WET_DRY_FILE_OBJECT() )
END IF
# if defined (BioGen)
NCF2 => BIO_FILE_OBJECT()
NCF => ADD(NCF,NCF2)
!!$ NCF => ADD(NCF,BIO_FILE_OBJECT() )
# endif
# if defined (WATER_QUALITY)
NCF2 => WQM_FILE_OBJECT()
NCF => ADD(NCF,NCF2)
!!$ NCF => ADD(NCF,WQM_FILE_OBJECT() )
# endif
# if defined (SEDIMENT)
IF(SEDIMENT_MODEL)THEN
NCF2 => SED_FILE_OBJECT()
NCF => ADD(NCF,NCF2)
!!$ NCF => ADD(NCF,SED_FILE_OBJECT() )
ENDIF
# endif
# if defined (WAVE_CURRENT_INTERACTION)
NCF2 => WAVE_PARA_FILE_OBJECT()
NCF => ADD(NCF,NCF2)
!!$ NCF => ADD(NCF,WAVE_PARA_FILE_OBJECT() )
NCF2 => WAVE_STRESS_FILE_OBJECT()
NCF => ADD(NCF,NCF2)
!!$ NCF => ADD(NCF,WAVE_STRESS_FILE_OBJECT() )
# endif
# if defined (NH)
NCF2 => NH_RST_FILE_OBJECT()
NCF => ADD(NCF,NCF2)
!!$ NCF => ADD(NCF,NH_RST_FILE_OBJECT() )
# endif
# if defined (ICE)
!-----------------------------------------------------------------
! state variables
NCF2 => ICE_RESTART_STATE_FILE_OBJECT()
NCF => ADD(NCF,NCF2)
!!$ NCF => ADD(NCF,ICE_RESTART_STATE_FILE_OBJECT() )
!-----------------------------------------------------------------
! velocity
!-----------------------------------------------------------------
NCF2 => ICE_VEL_FILE_OBJECT()
NCF => ADD(NCF,NCF2)
!!$ NCF => ADD(NCF,ICE_VEL_FILE_OBJECT() )
!-----------------------------------------------------------------
! fresh water, salt, and heat flux
!-----------------------------------------------------------------
NCF2 => ICE_EXTRA_FILE_OBJECT()
NCF => ADD(NCF,NCF2)
!!$ NCF => ADD(NCF,ICE_EXTRA_FILE_OBJECT() )
# endif
IF (STARTUP_TYPE /= "crashrestart") THEN
CALL NC_WRITE_FILE(NCF)
NCF%FTIME%NEXT_STKCNT = 1
ELSE
NCF%CONNECTED = .TRUE.
END IF
CALL KILL_DIMENSIONS
IF(DBG_SET(DBG_SBR)) WRITE(IPT,*) "END SETUP_RSTFILE"
END FUNCTION SETUP_RESTART
! FUNCTION SETUP_RESTART(FNAME) RESULT(NCF)
! USE MOD_NCDIO
! IMPLICIT NONE
! CHARACTER(LEN=*) :: FNAME
! TYPE(NCFILE), POINTER :: NCF
!
! IF(DBG_SET(DBG_SBR)) WRITE(IPT,*) "START: SETUP_RESTART"
!
!
! NCF => NEW_FILE()
! NCF%FNAME=trim(fname)
! print *, 'running setup_restart:',NCF%FNAME ! should be mark
! NCF%FTIME=>NEW_FTIME()
!
! NCF => ADD(NCF,GRID_FILE_OBJECT() )
!
! NCF => ADD(NCF,RESTART_EXTRAS_FILE_OBJECT() )
!
! NCF => ADD(NCF,VELOCITY_FILE_OBJECT() )
!
! NCF => ADD(NCF,AVERAGE_VEL_FILE_OBJECT() )
!
! NCF => ADD(NCF,VERTICAL_VEL_FILE_OBJECT() )
!
! NCF => ADD(NCF,TURBULENCE_FILE_OBJECT() )
!
! NCF => ADD(NCF,SALT_TEMP_FILE_OBJECT() )
!
! IF(WETTING_DRYING_ON) THEN
! NCF => ADD(NCF, WET_DRY_FILE_OBJECT() )
! END IF
!
! CALL NC_WRITE_FILE(NCF,LOCAL_DISK)
! NCF%FTIME%NEXT_STKCNT = 1
!
!
!
! IF(DBG_SET(DBG_SBR)) WRITE(IPT,*) "END SETUP_RESTART"
!
! END FUNCTION SETUP_RESTART
!
# endif
END MODULE MOD_ENKF