Commit 9de4289f authored by ymeur's avatar ymeur
Browse files

Improve nudging :

Vertical remapping can be done using XIOS axis interpolation or using inter dynamico remapping. Switch between this two choice using def parameter :
nudging_vertical_levels = 'pressure' / 'model'
'pressure' (default) : vertical level of forcing fields are in fixed pressure level, no need to read pressure field (or soil pressure) to make the guiding, so vertical remapping is done with DYNAMICO.
'model' : vertical level of forcing fields are in model level, so pressure forcing field is needed to make correct vertical remapping.
YM
parent 8a73a328
Pipeline #157230 failed with stages
in 29 seconds
......@@ -27,14 +27,15 @@ MODULE time_field_mod
& f_levels(:) ! temporary field on read level
LOGICAL :: is_ps ! surface pressure ps plays a special role when reading 3D input files
LOGICAL :: is_q !
INTEGER :: ndim
REAL(rstd),ALLOCATABLE :: levels(:) ! vertical axis avlue of the read field
INTEGER :: nlevels ! nb levels of read field
LOGICAL :: do_vertical_remap = .TRUE. ! do the vertical remapping (True) or let XIOS doing it (false)
END type t_time_field
TYPE(t_field), POINTER :: f_p(:) ! 3D pressure field reconstructed from ps and hybrid coordinate
TYPE(t_field), POINTER :: f_ps(:) ! ps
REAL(rstd),ALLOCATABLE :: levels(:)
INTEGER :: nlevels
PUBLIC :: t_time_interpolation, init_time_interpolation, update_time_interpolation, &
& t_time_field, init_time_field, allocate_time_field, &
& read_time_field, swap_time_field, interpolate_time_field
......@@ -51,14 +52,6 @@ CONTAINS
CALL allocate_field(f_p, field_t, type_real, llm, name='pressure_time_field')
CALL allocate_field(f_ps, field_t, type_real, name='pressure_surface_time_field')
IF (is_omp_master) CALL xios_get_axis_attr("axis_nudged",n_glo=nlevels)
CALL bcast_omp(nlevels)
ALLOCATE(levels(nlevels))
IF (is_omp_master) CALL xios_get_axis_attr("axis_nudged",value=levels)
CALL bcast_omp(levels)
! levels=levels*100 ! millibars -> Pa
END SUBROUTINE init_time_field
SUBROUTINE init_time_interpolation(interval, interp)
......@@ -67,13 +60,15 @@ CONTAINS
interp%interval = interval
END SUBROUTINE init_time_interpolation
SUBROUTINE allocate_time_field(field, field_type, data_type, dim3, dim4, name,suffix, pname)
SUBROUTINE allocate_time_field(do_vertical_remap, field, field_type, data_type, dim3, dim4, name,suffix, pname)
USE icosa
USE field_mod, ONLY : allocate_field
USE domain_mod, ONLY : ndomain, assigned_domain
USE xios_mod
USE omp_para
! value of argument pname must reflect the contents of XIOS XML configuration files
LOGICAL, INTENT(IN) :: do_vertical_remap
TYPE(t_time_field), INTENT(OUT) :: field
INTEGER, INTENT(IN) :: field_type
INTEGER, INTENT(IN) :: data_type
......@@ -81,7 +76,7 @@ CONTAINS
CHARACTER(*), OPTIONAL, INTENT(IN) :: name, suffix, pname
INTEGER :: ind
field%do_vertical_remap = do_vertical_remap
field%name = TRIM(name) ! note that name is in fact a compulsory argument
field%suffix =""
IF (PRESENT(suffix)) field%suffix = TRIM(suffix)
......@@ -89,11 +84,35 @@ CONTAINS
IF(field%is_ps) field%pname = TRIM(pname)
field%is_q = .FALSE.
IF (PRESENT(dim4) .AND. dim4 == nqtot) field%is_q = .TRUE.
CALL allocate_field(field%f_t, field_type, data_type, dim3, dim4, name=TRIM(name)//'_t')
CALL allocate_field(field%f_t0, field_type, data_type, dim3, dim4, name=TRIM(name)//'_t0')
CALL allocate_field(field%f_t1, field_type, data_type, dim3, dim4, name=TRIM(name)//'_t1')
CALL allocate_field(field%f_tmp, field_type, data_type, dim3, name=TRIM(name)//'_tmp')
CALL allocate_field(field%f_levels, field_type, data_type, nlevels, name=TRIM(name)//'_level')
IF (is_omp_master) CALL xios_get_axis_attr("axis_nudged",n_glo=field%nlevels)
CALL bcast_omp(field%nlevels)
ALLOCATE(field%levels(field%nlevels))
IF (is_omp_master) CALL xios_get_axis_attr("axis_nudged",value=field%levels)
CALL bcast_omp(field%levels)
IF (PRESENT(dim4)) THEN
field%ndim=4
ELSE IF (PRESENT(dim3)) THEN
field%ndim=3
ELSE
field%ndim=2
ENDIF
IF (field%do_vertical_remap .AND. PRESENT(dim3)) THEN
CALL allocate_field(field%f_t0, field_type, data_type, field%nlevels, dim4, name=TRIM(name)//'_t0')
CALL allocate_field(field%f_t1, field_type, data_type, field%nlevels, dim4, name=TRIM(name)//'_t1')
CALL allocate_field(field%f_tmp, field_type, data_type, field%nlevels, name=TRIM(name)//'_tmp')
CALL allocate_field(field%f_levels, field_type, data_type, field%nlevels, dim4, name=TRIM(name)//'_level')
CALL allocate_field(field%f_t, field_type, data_type, dim3, dim4, name=TRIM(name)//'_t')
ELSE
CALL allocate_field(field%f_t0, field_type, data_type, dim3, dim4, name=TRIM(name)//'_t0')
CALL allocate_field(field%f_t1, field_type, data_type, dim3, dim4, name=TRIM(name)//'_t1')
CALL allocate_field(field%f_tmp, field_type, data_type, dim3, name=TRIM(name)//'_tmp')
CALL allocate_field(field%f_levels, field_type, data_type, field%nlevels, dim4, name=TRIM(name)//'_level')
CALL allocate_field(field%f_t, field_type, data_type, dim3, dim4, name=TRIM(name)//'_t')
ENDIF
DO ind=1, ndomain
IF(.NOT. assigned_domain(ind)) CYCLE
......@@ -199,32 +218,35 @@ CONTAINS
IF(.NOT. assigned_domain(ind)) CYCLE
CALL swap_dimensions(ind)
CALL swap_geometry(ind)
DO l=ll_begin,ll_end
field%f_t0(ind)%rval4d(:,l,iq) = field%f_tmp(ind)%rval4d(:,l,1)
ENDDO
! DO l=ll_begin,ll_end
field%f_t0(ind)%rval4d(:,:,iq) = field%f_tmp(ind)%rval4d(:,:,1)
! ENDDO
ENDDO
ENDDO
ELSE
CALL xios_read_field(TRIM(field%name)//"0", field%f_t0)
! CALL writeField(TRIM(field%name),field%f_t0)
ENDIF
SELECT CASE(field%f_t0(1)%field_type)
CASE(field_T)
CALL transfert_request(field%f_t0, req_i0)
CALL transfert_request(field%f_t0, req_i1)
IF (.NOT. field%do_vertical_remap) THEN
CALL transfert_request(field%f_t0, req_i0)
CALL transfert_request(field%f_t0, req_i1)
ENDIF
CASE DEFAULT
STOP 'read_time_field : unsupported field_type'
END SELECT
IF(field%is_ps) THEN
DO ind=1,ndomain
f_ps(ind)%rval4d = field%f_t0(ind)%rval4d
ENDDO
CALL pression_mid(field%f_t0, f_p)
CALL xios_write_field(TRIM(field%pname)//"0",f_p)
! CALL writeField(TRIM(field%pname),f_p)
IF (.NOT. field%do_vertical_remap) THEN
IF(field%is_ps) THEN
DO ind=1,ndomain
f_ps(ind)%rval4d = field%f_t0(ind)%rval4d
ENDDO
CALL pression_mid(field%f_t0, f_p)
CALL xios_write_field(TRIM(field%pname)//"0",f_p)
ENDIF
END IF
END IF
......@@ -245,178 +267,70 @@ CONTAINS
IF(.NOT. assigned_domain(ind)) CYCLE
CALL swap_dimensions(ind)
CALL swap_geometry(ind)
DO l=ll_begin,ll_end
field%f_t1(ind)%rval4d(:,l,iq) = field%f_tmp(ind)%rval4d(:,l,1)
ENDDO
! DO l=ll_begin,ll_end
field%f_t1(ind)%rval4d(:,:,iq) = field%f_tmp(ind)%rval4d(:,:,1)
! ENDDO
ENDDO
ENDDO
ELSE
CALL xios_read_field(TRIM(field%name), field%f_t1)
! CALL writeField(TRIM(field%name),field%f_t1)
ENDIF
SELECT CASE(field%f_t1(1)%field_type)
CASE(field_T)
CALL transfert_request(field%f_t1, req_i0)
CALL transfert_request(field%f_t1, req_i1)
CASE DEFAULT
STOP 'read_time_field : unsupported field_type'
END SELECT
IF(field%is_ps) THEN
CALL pression_mid(field%f_t1, f_p)
CALL xios_write_field(field%pname, f_p)
! CALL writeField(TRIM(field%pname),f_p)
END IF
ENDIF
END SUBROUTINE read_time_field
SUBROUTINE read_time_field_remap(interp, field)
! surface pressure ps must always be read before 3D fields
! in the special case that we read ps, we compute the 3D pressure field
! and provide it to XIOS for later use in vertical interpolation of 3D fields
! this special role of ps is marked by the flag is_ps
USE layout_mod, ONLY : field_T, field_U, field_Z
USE transfert_mod, ONLY : transfert_request, req_i0, req_i1
USE pression_mod, ONLY : pression_mid
USE xios_mod, ONLY : xios_read_field, xios_write_field, xios_is_valid_field
USE tracer_icosa_mod, ONLY : tracers
USE icosa
USE write_field_mod
USE omp_para
USE vertical_remap_mod
TYPE(t_time_interpolation), INTENT(IN) :: interp
TYPE(t_time_field), INTENT(INOUT) :: field
REAL(rstd) :: t0,t1
INTEGER :: iq,l,ind
IF(interp%read_t0) THEN
! we get data for t0 only once, from XIOS field XXX0
IF (field%is_q) THEN ! for tracers
DO iq=1,nqtot
! CALL xios_read_field(TRIM(tracers(iq)%name)//TRIM(field%suffix)//"0", field%f_tmp)
IF (xios_is_valid_field(TRIM(tracers(iq)%name)//TRIM(field%suffix))) THEN
CALL xios_read_field(TRIM(tracers(iq)%name)//TRIM(field%suffix)//"0", field%f_levels)
ELSE
DO ind=1,ndomain
field%f_levels(ind)%rval4d(:,:,:)=0
ENDDO
ENDIF
CALL vertical_remap(levels,field%f_levels,f_ps,field%f_tmp)
DO ind=1,ndomain
IF(.NOT. assigned_domain(ind)) CYCLE
CALL swap_dimensions(ind)
CALL swap_geometry(ind)
DO l=ll_begin,ll_end
field%f_t0(ind)%rval4d(:,l,iq) = field%f_tmp(ind)%rval4d(:,l,1)
ENDDO
ENDDO
ENDDO
ELSE
IF(field%is_ps) THEN
CALL xios_read_field(TRIM(field%name)//"0", field%f_t0)
ELSE
CALL xios_read_field(TRIM(field%name)//"0", field%f_levels)
CALL vertical_remap(levels,field%f_levels,f_ps,field%f_t0)
IF (.NOT. field%do_vertical_remap) THEN
CALL transfert_request(field%f_t1, req_i0)
CALL transfert_request(field%f_t1, req_i1)
ENDIF
CALL writeField(TRIM(field%name),field%f_t0)
ENDIF
SELECT CASE(field%f_t0(1)%field_type)
CASE(field_T)
CALL transfert_request(field%f_t0, req_i0)
CALL transfert_request(field%f_t0, req_i1)
CASE DEFAULT
STOP 'read_time_field : unsupported field_type'
END SELECT
IF(field%is_ps) THEN
DO ind=1,ndomain
f_ps(ind)%rval4d = field%f_t0(ind)%rval4d
ENDDO
CALL pression_mid(field%f_t0, f_p)
CALL xios_write_field(TRIM(field%pname)//"0",f_p)
CALL writeField(TRIM(field%pname),f_p)
END IF
END IF
IF(interp%swap) CALL swap_time_field(field)
IF(interp%read_t1) THEN
IF (field%is_q) THEN ! for tracers
DO iq=1,nqtot
! CALL xios_read_field(TRIM(tracers(iq)%name)//TRIM(field%suffix), field%f_tmp)
IF (xios_is_valid_field(TRIM(tracers(iq)%name)//TRIM(field%suffix))) THEN
CALL xios_read_field(TRIM(tracers(iq)%name)//TRIM(field%suffix), field%f_levels)
ELSE
DO ind=1,ndomain
field%f_levels(ind)%rval4d(:,:,:)=0
ENDDO
ENDIF
CALL vertical_remap(levels,field%f_levels,f_ps,field%f_tmp)
DO ind=1,ndomain
IF(.NOT. assigned_domain(ind)) CYCLE
CALL swap_dimensions(ind)
CALL swap_geometry(ind)
DO l=ll_begin,ll_end
field%f_t1(ind)%rval4d(:,l,iq) = field%f_tmp(ind)%rval4d(:,l,1)
ENDDO
ENDDO
ENDDO
ELSE
IF (.NOT. field%do_vertical_remap) THEN
IF(field%is_ps) THEN
CALL xios_read_field(TRIM(field%name), field%f_t1)
ELSE
CALL xios_read_field(TRIM(field%name), field%f_levels)
CALL vertical_remap(levels,field%f_levels,f_ps,field%f_t1)
CALL pression_mid(field%f_t1, f_p)
CALL xios_write_field(field%pname, f_p)
ENDIF
! CALL writeField(TRIM(field%name),field%f_t1)
ENDIF
SELECT CASE(field%f_t1(1)%field_type)
CASE(field_T)
CALL transfert_request(field%f_t1, req_i0)
CALL transfert_request(field%f_t1, req_i1)
CASE DEFAULT
STOP 'read_time_field : unsupported field_type'
END SELECT
IF(field%is_ps) THEN
CALL pression_mid(field%f_t1, f_p)
CALL xios_write_field(field%pname, f_p)
! CALL writeField(TRIM(field%pname),f_p)
END IF
ENDIF
END SUBROUTINE read_time_field_remap
SUBROUTINE interpolate_time_field(interp, field)
END SUBROUTINE read_time_field
SUBROUTINE interpolate_time_field(interp, field, f_ps)
USE field_mod, ONLY : t_field
USE domain_mod, ONLY : ndomain, assigned_domain
TYPE(t_time_interpolation), INTENT(IN) :: interp
TYPE(t_time_field), INTENT(INOUT) :: field
USE vertical_remap_mod
IMPLICIT NONE
TYPE(t_time_interpolation), INTENT(IN) :: interp
TYPE(t_time_field), INTENT(INOUT) :: field
TYPE(t_field), POINTER,OPTIONAL,INTENT(IN) :: f_ps(:)
INTEGER :: ind, n
DO ind=1, ndomain
DO ind=1, ndomain
IF(.NOT. assigned_domain(ind)) CYCLE
n = SIZE(field%f_t(ind)%rval4d)
CALL compute_interpolate_time_field(interp%coef0, interp%coef1, n, &
field%f_t0(ind)%rval4d, field%f_t1(ind)%rval4d, field%f_t(ind)%rval4d)
END DO
IF (field%do_vertical_remap .AND. field%ndim>2) THEN
n = SIZE(field%f_levels(ind)%rval4d)
CALL compute_interpolate_time_field(interp%coef0, interp%coef1, n, &
field%f_t0(ind)%rval4d, field%f_t1(ind)%rval4d, field%f_levels(ind)%rval4d)
ELSE
n = SIZE(field%f_t(ind)%rval4d)
CALL compute_interpolate_time_field(interp%coef0, interp%coef1, n, &
field%f_t0(ind)%rval4d, field%f_t1(ind)%rval4d, field%f_t(ind)%rval4d)
ENDIF
END DO
IF (field%do_vertical_remap) THEN
IF (field%ndim==2) THEN
! no vertical interpolation ; do nothing
ELSE
CALL vertical_remap(field%levels,field%f_levels, f_ps, field%f_t)
ENDIF
ENDIF
! CALL writeField(TRIM(field%name)//"_time", field%f_t)
END SUBROUTINE interpolate_time_field
SUBROUTINE compute_interpolate_time_field(coef0, coef1, n, f0, f1, f)
......
......@@ -6,6 +6,30 @@ MODULE guided_mod
CONTAINS
SUBROUTINE pre_init_guided
USE icosa
USE nudging_mod, ONLY : pre_init_guided_nudging => pre_init_guided
IMPLICIT NONE
guided_type='none'
CALL getin("guided_type",guided_type)
SELECT CASE(TRIM(guided_type))
CASE ('none')
CASE ('dcmip1')
CASE ('nudging')
CALL pre_init_guided_nudging
CASE DEFAULT
PRINT*,"Bad selector for varaible guided_type >",TRIM(guided_type),"> option are <none>, <dcmip1>, <nudging>"
STOP
END SELECT
END SUBROUTINE pre_init_guided
SUBROUTINE init_guided
USE icosa
USE guided_ncar_mod, ONLY : init_guided_ncar => init_guided
......@@ -25,7 +49,7 @@ CONTAINS
CALL init_guided_nudging
CASE DEFAULT
PRINT*,"Bad selector for varaible guided_type >",TRIM(guided_type),"> option are <none>, <dcmip1>"
PRINT*,"Bad selector for varaible guided_type >",TRIM(guided_type),"> option are <none>, <dcmip1>, <nudging>"
STOP
END SELECT
......@@ -51,7 +75,7 @@ CONTAINS
CASE ('nudging')
CALL guided_nudging(tt, f_ps, f_mass, f_theta_rhodz, f_u, f_q)
CASE DEFAULT
PRINT*,"Bad selector for varaible guided_type >",TRIM(guided_type),"> option are <none>, <dcmip1>"
PRINT*,"Bad selector for varaible guided_type >",TRIM(guided_type),"> option are <none>, <dcmip1>, <nudging>"
STOP
END SELECT
......
......@@ -57,17 +57,21 @@ MODULE nudging_mod
INTEGER :: itau_nudging
LOGICAL,SAVE :: do_vertical_remap
!$OMP THREADPRIVATE(do_vertical_remap)
PUBLIC :: init_guided, guided, update_border
PUBLIC :: pre_init_guided, init_guided, guided, update_border
CONTAINS
SUBROUTINE init_guided
SUBROUTINE pre_init_guided
USE getin_mod, ONLY : getin
USE time_field_mod, ONLY : init_time_interpolation, init_time_field, allocate_time_field
USE icosa
USE omp_para
USE xios_mod
CHARACTER(LEN=255) :: nudging_zone_str
CHARACTER(LEN=255) :: nudging_vertical_levels
CALL init_time_field
......@@ -87,15 +91,6 @@ CONTAINS
PRINT*,"Bad selector for variable <nudging_zone> ",TRIM(nudging_zone_str), " option are <lonlat>, <circular>, <area>"
STOP
END SELECT
CALL init_time_interpolation(itau_nudging, interp)
CALL allocate_time_field(tf_ps, field_T, type_real, name='ps_nudged', pname='p_nudged')
CALL allocate_time_field(tf_ulon, field_T, type_real, llm, name='ulon_nudged')
CALL allocate_time_field(tf_ulat, field_T, type_real, llm, name='ulat_nudged')
CALL allocate_time_field(tf_u, field_U, type_real, llm, name='u_nudged')
CALL allocate_time_field(tf_temp, field_t, type_real, llm, name='temp_nudged')
CALL allocate_time_field(tf_q, field_T, type_real, llm, nqtot, name='q', suffix="_nudged")
CALL allocate_time_field(tf_theta_rhodz, field_T, type_real, llm, 1, name='theta_rhodz_nudged')
CALL getin("guide_T",guide_T)
T_relax_in=-1
......@@ -145,6 +140,55 @@ CONTAINS
ENDIF
q_relax_out = dt/q_relax_out
nudging_vertical_levels = 'pressure'
CALL getin("nudging_vertical_levels", nudging_vertical_levels)
IF (TRIM(nudging_vertical_levels)=="pressure") THEN
do_vertical_remap = .TRUE.
ELSE IF (TRIM(nudging_vertical_levels)=="model") THEN
do_vertical_remap = .FALSE.
ELSE
PRINT*, "bad option for 'nudging_vertical_levels' parameter : possible option are <pressure>, <model> and get : " &
,TRIM(nudging_vertical_levels)
ENDIF
! activate nudging in XIOS workflow
CALL xios_set_filegroup_attr("nudging_files",enabled=.TRUE.)
CALL xios_set_fieldgroup_attr("nudging_fields",read_access=.TRUE.)
CALL xios_set_fieldgroup_attr("nudging_fields0",read_access=.TRUE.)
! switch XIOS worklow to do or not vertical interpolation in the workflow
IF (do_vertical_remap) THEN
! vertical remapping is done with the model
CALL xios_set_axis_attr("from_nudged",axis_ref="axis_nudged")
CALL xios_set_axis_attr("from_nudged0",axis_ref="axis_nudged")
CALL xios_set_field_attr("p_nudged_horiz_interp",read_access=.FALSE.)
CALL xios_set_field_attr("p_nudged_horiz_interp0",read_access=.FALSE.)
ELSE
! vertical remapping is done in XIOS workflow
CALL xios_set_axis_attr("from_nudged",axis_ref="nudged_vert_interp_src_dest")
CALL xios_set_axis_attr("from_nudged0",axis_ref="nudged_vert_interp_src_dest")
CALL xios_set_field_attr("p_nudged_horiz_interp",read_access=.TRUE.)
CALL xios_set_field_attr("p_nudged_horiz_interp0",read_access=.TRUE.)
ENDIF
END SUBROUTINE pre_init_guided
SUBROUTINE init_guided
USE time_field_mod, ONLY : init_time_interpolation, init_time_field, allocate_time_field
USE icosa
IMPLICIT NONE
CALL init_time_interpolation(itau_nudging, interp)
CALL allocate_time_field(do_vertical_remap, tf_ps, field_T, type_real, name='ps_nudged', pname='p_nudged')
CALL allocate_time_field(do_vertical_remap, tf_ulon, field_T, type_real, llm, name='ulon_nudged')
CALL allocate_time_field(do_vertical_remap, tf_ulat, field_T, type_real, llm, name='ulat_nudged')
CALL allocate_time_field(do_vertical_remap, tf_u, field_U, type_real, llm, name='u_nudged')
CALL allocate_time_field(do_vertical_remap, tf_temp, field_t, type_real, llm, name='temp_nudged')
CALL allocate_time_field(do_vertical_remap, tf_q, field_T, type_real, llm, nqtot, name='q', suffix="_nudged")
CALL allocate_time_field(do_vertical_remap, tf_theta_rhodz, field_T, type_real, llm, 1, name='theta_rhodz_nudged')
SELECT CASE(nudging_zone)
CASE(nudging_lonlat)
......@@ -157,9 +201,9 @@ CONTAINS
CALL init_nudging_global
CASE DEFAULT
END SELECT
END SUBROUTINE init_guided
SUBROUTINE init_nudging_lonlat
USE icosa
......@@ -209,6 +253,7 @@ CONTAINS
stiffness = 8
CALL getin("nudging_stiffness", stiffness)
ncell = 0
cell_size = 0
......@@ -506,7 +551,8 @@ CONTAINS
USE transfert_omp_mod
USE output_field_mod
USE omp_para
USE write_field_mod
USE xios_mod
IMPLICIT NONE
TYPE(t_field),POINTER,SAVE :: f_neighbours(:)
REAL(rstd),POINTER :: neighbours(:)
......@@ -598,7 +644,7 @@ CONTAINS
coeff_i(ij)=1
border(ij) = 0
ELSE
coeff_i(ij) = 0.5*tanh(-2.5*(neighbours(ij)-lam_halo_scheme-stiffness/2)/stiffness)+0.5
coeff_i(ij) = 1-tanh(atanh(1-0.02)*(neighbours(ij)-lam_halo_scheme)/stiffness)
border(ij) = 1
ENDIF
ENDDO
......@@ -611,7 +657,7 @@ CONTAINS
CALL transfert_request(f_coeff_i,req_i1)
CALL transfert_request(f_border,req_i0)
CALL transfert_request(f_border,req_i1)
IF (is_omp_master) THEN
DO ind = 1 , ndomain
CALL swap_dimensions(ind)
......@@ -674,17 +720,20 @@ CONTAINS
END SUBROUTINE init_nudging_global
SUBROUTINE check_time_interp(time)
SUBROUTINE check_time_interp(time, f_ps)
USE output_field_mod, ONLY : output_field
USE xios_mod, ONLY : xios_write_field
USE time_field_mod, ONLY : update_time_interpolation, read_time_field, swap_time_field, interpolate_time_field
USE time_field_mod, ONLY : update_time_interpolation, read_time_field, swap_time_field,&
interpolate_time_field
USE pression_mod
USE wind_mod
USE theta2theta_rhodz_mod
USE write_field_mod
IMPLICIT NONE
REAL(rstd), INTENT(IN) :: time
TYPE(t_field), POINTER :: f_ps(:)
INTEGER :: itau
itau=NINT(time/dt)
IF (last_itau /= itau) THEN