mode_ekf.F90

Go to the documentation of this file.

!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
!SFX_LIC for details. version 1.
MODULE mode_ekf
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
INTERFACE get_file_name
  MODULE PROCEDURE get_file_name
END INTERFACE
INTERFACE b_big_loop
  MODULE PROCEDURE b_big_loop
END INTERFACE
INTERFACE inverse_matrix
  MODULE PROCEDURE inverse_matrix
END INTERFACE
INTERFACE choldc
  MODULE PROCEDURE choldc
END INTERFACE
INTERFACE cholsl
  MODULE PROCEDURE cholsl
END INTERFACE
INTERFACE cofswi
  MODULE PROCEDURE cofswi
END INTERFACE
INTERFACE set_filein
  MODULE PROCEDURE set_filein
END INTERFACE
!
 CONTAINS
!
SUBROUTINE get_file_name(KYEAR,KMONTH,KDAY,KHOUR,HMFILE)
!
USE modi_trans_chaine
!
IMPLICIT NONE
!
INTEGER, INTENT (IN) :: kyear, kmonth, kday, khour
 CHARACTER(LEN=50) :: ymonth,yday,yyear,yhour
 CHARACTER(LEN=200), INTENT(INOUT) :: hmfile
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('MODE_EKF:GET_FILE_NAME',0,zhook_handle)
!
 CALL trans_chaine(ymonth,kmonth,2)
 CALL trans_chaine(yday,kday,2)
 CALL trans_chaine(yyear,kyear,4)
 CALL trans_chaine(yhour,khour,2)
IF (khour.EQ.0 .OR. khour.EQ.24 .OR. khour.EQ.48) yhour='00'
!
hmfile = trim(hmfile)//yyear(3:4)
hmfile = trim(hmfile)//ymonth
hmfile = trim(hmfile)//yday
hmfile = trim(hmfile)//'H'
hmfile = trim(hmfile)//yhour
hmfile = trim(hmfile)
!
IF (lhook) CALL dr_hook('MODE_EKF:GET_FILE_NAME',1,zhook_handle)
!
END SUBROUTINE get_file_name
!
!
SUBROUTINE b_big_loop (I, &
                       haction,hfile,ptab,ptab_in)
!
!
USE modd_isba_n, ONLY : isba_t
!
USE modd_assim, ONLY : nvar
!
IMPLICIT NONE
!
!
TYPE(isba_t), INTENT(INOUT) :: i
!
 CHARACTER(LEN=4), INTENT(IN) :: haction
 CHARACTER(LEN=*), INTENT(IN) :: hfile
REAL, DIMENSION(:,:,:), INTENT(INOUT) :: ptab
REAL, DIMENSION(:), INTENT(IN), OPTIONAL :: ptab_in
INTEGER :: ji,jj,jk,jl,jjj,l1,k1
INTEGER :: iki
!RJ: missing declarations
INTEGER :: icpt,istat
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('MODE_EKF:B_BIG_LOOP',0,zhook_handle)
!
IF (haction=="READ" .OR. haction=="WRIT") THEN
  OPEN (unit=111,file=hfile,form='UNFORMATTED',access='SEQUENTIAL',iostat=istat)
ELSE
  icpt = 0
ENDIF
!
iki = SIZE(ptab,1)
!
DO jl = 1,nvar   ! control variable (x at previous time step)
  DO jk = 1,nvar
    DO ji = 1,iki
      DO jj = 1,i%NPATCH   
        DO jjj = 1,i%NPATCH
          !
          l1 = jj+i%NPATCH*(jl-1)
          k1 = jjj+i%NPATCH*(jk-1)
          !
          IF ( haction=="READ" ) THEN
            READ (111) ptab(ji,l1,k1)
          ELSEIF ( haction=="WRIT" ) THEN
            WRITE(111) ptab(ji,l1,k1)
          ELSE
            icpt = icpt+1
            ptab(ji,l1,k1) = ptab_in(icpt)
          ENDIF
          ! 
        ENDDO
      ENDDO
    ENDDO
  ENDDO
ENDDO
!
IF (haction=="READ" .OR. haction=="WRIT") CLOSE(111)
!
IF (lhook) CALL dr_hook('MODE_EKF:B_BIG_LOOP',1,zhook_handle)
!
END SUBROUTINE b_big_loop

SUBROUTINE inverse_matrix(KN,PA,PP)
!--------------------------------------------------------
!
! Explicit inversed matrix after Cholesky decomposition
!
!--------------------------------------------------------
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: kn
REAL, DIMENSION (KN,KN), INTENT(INOUT) :: pa
REAL, DIMENSION (KN),    INTENT(IN)    :: pp
REAL :: zsum
INTEGER :: ji, jj, jk
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('MODE_EKF:INVERSE_MATRIX',0,zhook_handle)
!
DO ji = 1,kn
  pa(ji,ji)=1./pp(ji)
  DO jj = ji+1,kn
    zsum = 0.
    DO jk = ji,jj-1
      zsum = zsum - pa(jj,jk)*pa(jk,ji)
    ENDDO
    pa(jj,ji) = zsum/pp(jj)
  ENDDO
ENDDO  
!
DO ji = 1,kn
  DO jj = ji+1,kn
    pa(ji,jj) = 0.0
  ENDDO
ENDDO
!
pa = matmul(transpose(pa),pa)
!
IF (lhook) CALL dr_hook('MODE_EKF:INVERSE_MATRIX',1,zhook_handle)
!
END SUBROUTINE inverse_matrix

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE choldc(KN,PA,PP)
!
USE modi_abor1_sfx        
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: kn
REAL, DIMENSION(KN,KN), INTENT(INOUT) :: pa
REAL, DIMENSION(KN), INTENT(OUT) :: pp
INTEGER :: ji
REAL :: zsum
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('MODE_EKF:CHOLDC',0,zhook_handle)
!
DO ji = 1,kn
  zsum = pa(ji,ji)- dot_product(pa(ji,1:ji-1),pa(ji,1:ji-1))
  IF ( zsum<=0.0 ) CALL abor1_sfx('MODE_EKF: ERROR IN CHOLDC')
  pp(ji) = sqrt(zsum) 
  pa(ji+1:kn,ji) = ( pa(ji,ji+1:kn) - matmul(pa(ji+1:kn,1:ji-1),pa(ji,1:ji-1)) )/pp(ji)
ENDDO
!
IF (lhook) CALL dr_hook('MODE_EKF:CHOLDC',1,zhook_handle)
!
END SUBROUTINE choldc 

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE cholsl(KN,PA,PP,PB,PX)
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: kn
REAL, DIMENSION (KN,KN), INTENT(IN) :: pa
REAL, DIMENSION (KN), INTENT(IN)    :: pp,pb
REAL, DIMENSION (KN), INTENT(INOUT) :: px
INTEGER :: ji
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('MODE_EKF:CHOLSL',0,zhook_handle)
!
DO ji = 1,kn
  px(ji) = (pb(ji) - dot_product(pa(ji,1:ji-1),px(1:ji-1)))/pp(ji)
ENDDO
DO ji = kn,1,-1
  px(ji) = (px(ji) - dot_product(pa(ji+1:kn,ji),px(ji+1:kn)))/pp(ji)
ENDDO
!
IF (lhook) CALL dr_hook('MODE_EKF:CHOLSL',1,zhook_handle)
!
END SUBROUTINE cholsl
!
SUBROUTINE cofswi(PCLAY,PCOFSWI)
!
IMPLICIT NONE
REAL, DIMENSION(:), INTENT(IN) :: pclay
REAL, DIMENSION(:), INTENT(OUT) :: pcofswi
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('MODE_EKF:COFSWI',0,zhook_handle)
!
pcofswi(:) = 0.001 * (89.0467 * ((100.*pclay(:))**0.3496) - 37.1342*((100.*pclay(:))**0.5))
!
IF (lhook) CALL dr_hook('MODE_EKF:COFSWI',1,zhook_handle)
!
END SUBROUTINE cofswi
!
!
SUBROUTINE set_filein(HFILE)
!
USE modn_io_offline, ONLY : csurf_filetype
!
#ifdef SFX_NC
USE modd_io_surf_nc,   ONLY : cfilein_nc, cfilein_nc_save
#endif
#ifdef SFX_ASC
USE modd_io_surf_asc,  ONLY : cfilein, cfilein_save
#endif
#ifdef SFX_FA
USE modd_io_surf_fa,   ONLY : cfilein_fa, cfilein_fa_save
#endif
#ifdef SFX_LFI
USE modd_io_surf_lfi,    ONLY : cfilein_lfi, cfilein_lfi_save, cfile_lfi
#endif
!
USE modi_abor1_sfx
!
IMPLICIT NONE
!
 CHARACTER(LEN=*), INTENT(IN) :: hfile
 CHARACTER(LEN=28) :: yfilein
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('MODE_EKF:SET_FILEIN',0,zhook_handle)
!
IF ( csurf_filetype == "LFI   " ) THEN
  yfilein = trim(hfile)
#ifdef SFX_LFI
  cfilein_lfi      = yfilein 
  cfile_lfi        = yfilein
  cfilein_lfi_save = yfilein
#endif
ELSEIF ( csurf_filetype == "FA    " ) THEN
  yfilein = trim(hfile)//'.fa'
#ifdef SFX_FA
  cfilein_fa      = yfilein
  cfilein_fa_save = yfilein
#endif
ELSEIF ( csurf_filetype == "ASCII " ) THEN
  yfilein = trim(hfile)//'.txt'
#ifdef SFX_ASC
  cfilein      = yfilein
  cfilein_save = yfilein
#endif
ELSEIF ( csurf_filetype == "NC    " ) THEN
  yfilein = trim(hfile)//'.nc'
#ifdef SFX_NC
  cfilein_nc      = yfilein
  cfilein_nc_save = yfilein
#endif
ELSE
  CALL abor1_sfx(trim(csurf_filetype)//" is not implemented!")
ENDIF
!
IF (lhook) CALL dr_hook('MODE_EKF:SET_FILEIN',1,zhook_handle)
!
END SUBROUTINE set_filein
!
END MODULE mode_ekf