From bd57a5f2603617636e2dca5c400390461ad5197a Mon Sep 17 00:00:00 2001
From: Philippe WAUTELET <philippe.wautelet@aero.obs-mip.fr>
Date: Wed, 14 Oct 2020 14:48:04 +0200
Subject: [PATCH] Philippe 14/10/2020: IO: add new Les_diachro_* subroutines
---
src/MNH/mode_les_diachro.f90 | 413 ++++++++++++++++++++++++++++++++++-
1 file changed, 404 insertions(+), 9 deletions(-)
diff --git a/src/MNH/mode_les_diachro.f90 b/src/MNH/mode_les_diachro.f90
index 0b799994e..aa55d6390 100644
--- a/src/MNH/mode_les_diachro.f90
+++ b/src/MNH/mode_les_diachro.f90
@@ -24,9 +24,15 @@ implicit none
private
-public :: LES_DIACHRO, Les_diachro_2pt, LES_DIACHRO_MASKS, Les_diachro_spec, &
+public :: Les_diachro, Les_diachro_2pt, LES_DIACHRO_MASKS, Les_diachro_spec, &
LES_DIACHRO_SURF, LES_DIACHRO_SURF_SV, LES_DIACHRO_SV, LES_DIACHRO_SV_MASKS
+interface Les_diachro
+ module procedure LES_DIACHRO_old
+ module procedure Les_diachro_1D, Les_diachro_2D, Les_diachro_3D, Les_diachro_4D
+end interface
+
+
CONTAINS
!
!---------------------------------------------------------------------
@@ -799,7 +805,7 @@ else
end if
! Write the profile
-if ( iresp==0 .and. ( gsurf .or. any( zwork6 /= xundef ) ) ) then
+if ( iresp==0 .and. ( gsurf .or. any( zwork6 /= XUNDEF ) ) ) then
allocate( tzfields( Size( pfield, 3 ) ) )
tzfields(:)%cmnhname = ytitle(:)
@@ -814,11 +820,11 @@ if ( iresp==0 .and. ( gsurf .or. any( zwork6 /= xundef ) ) ) then
tzfields(jp)%ndimlist(:) = NMNHDIM_UNKNOWN
end do
- call Write_diachro( tpdiafile, tzfields, ygroup, "SSOL", tzdates, &
- zwork6, &
- oicp = .false., ojcp = .false., okcp = .false., &
- kil = iil, kih = iih, kjl = ijl, kjh = ijh, kkl = ikl, kkh = ikh, &
- ptrajx = ztrajx, ptrajy = ztrajy, ptrajz = ztrajz )
+ call Write_diachro( tpdiafile, tzfields, ygroup, "SSOL", tzdates, &
+ zwork6, &
+ oicp = .false., ojcp = .false., okcp = .false., &
+ kil = iil, kih = iih, kjl = ijl, kjh = ijh, kkl = ikl, kkh = ikh, &
+ ptrajx = ztrajx, ptrajy = ztrajy, ptrajz = ztrajz )
deallocate( tzfields )
end if
@@ -826,8 +832,396 @@ end if
!-------------------------------------------------------------------------------
end subroutine Les_diachro_gen
!------------------------------------------------------------------------------
+
+!#######################################################################
+subroutine Les_diachro_1D( tpdiafile, tpfield, odoavg, odonorm, pfield )
+!#######################################################################
+
+use modd_field, only: NMNHDIM_BUDGET_LES_TIME, NMNHDIM_UNUSED, tfield_metadata_base
+use modd_io, only: tfiledata
+
+type(tfiledata), intent(in) :: tpdiafile ! File to write
+type(tfield_metadata_base), intent(in) :: tpfield ! Metadata of field
+logical, intent(in) :: odoavg ! Compute and store time average
+logical, intent(in) :: odonorm ! Compute and store normalized field
+real, dimension(:), intent(in) :: pfield ! Data array
+
+type(tfield_metadata_base) :: tzfield
+
+tzfield = tpfield
+
+if ( tzfield%ndims /= 1 ) then
+ call Print_msg( NVERB_ERROR, 'IO', 'Les_diachro_1D', 'ndims /= 1 for ' // Trim( tzfield%cmnhname ) )
+ tzfield%ndims = 1
+end if
+
+if ( Any( tzfield%ndimlist(2:) /= NMNHDIM_UNUSED ) ) then
+ call Print_msg( NVERB_ERROR, 'IO', 'Les_diachro_1D', 'unexpected type for some dimensions of ' &
+ // Trim( tzfield%cmnhname ) )
+ tzfield%ndimlist(2:) = NMNHDIM_UNUSED
+end if
+
+if ( tzfield%ndimlist(1) == NMNHDIM_BUDGET_LES_TIME ) then
+ tzfield%ndimlist(2) = tzfield%ndimlist(1)
+ tzfield%ndimlist(1) = NMNHDIM_UNUSED
+ tzfield%ndimlist(3) = NMNHDIM_UNUSED
+ tzfield%ndimlist(4) = NMNHDIM_UNUSED
+ call Les_diachro_common( tpdiafile, tzfield, reshape( pfield, [ 1, size( pfield, 1 ), 1, 1 ] ), &
+ odoavg, odonorm )
+else
+ call Print_msg( NVERB_ERROR, 'IO', 'Les_diachro_1D', &
+ 'ndimlist configuration not yet implemented for ' // Trim( tzfield%cmnhname ) )
+end if
+
+end subroutine Les_diachro_1D
+
+!#######################################################################
+subroutine Les_diachro_2D( tpdiafile, tpfield, odoavg, odonorm, pfield )
+!#######################################################################
+
+use modd_field, only: NMNHDIM_BUDGET_LES_LEVEL, NMNHDIM_BUDGET_LES_SV, NMNHDIM_BUDGET_LES_TIME, NMNHDIM_UNUSED, &
+ tfield_metadata_base
+use modd_io, only: tfiledata
+
+type(tfiledata), intent(in) :: tpdiafile ! File to write
+type(tfield_metadata_base), intent(in) :: tpfield ! Metadata of field
+logical, intent(in) :: odoavg ! Compute and store time average
+logical, intent(in) :: odonorm ! Compute and store normalized field
+real, dimension(:,:), intent(in) :: pfield ! Data array
+
+type(tfield_metadata_base) :: tzfield
+
+tzfield = tpfield
+
+if ( tzfield%ndims /= 2 ) then
+ call Print_msg( NVERB_ERROR, 'IO', 'Les_diachro_2D', 'ndims /= 2 for ' // Trim( tzfield%cmnhname ) )
+ tzfield%ndims = 2
+end if
+
+if ( Any( tzfield%ndimlist(3:) /= NMNHDIM_UNUSED ) ) then
+ call Print_msg( NVERB_ERROR, 'IO', 'Les_diachro_2D', 'unexpected type for some dimensions of ' &
+ // Trim( tzfield%cmnhname ) )
+ tzfield%ndimlist(3:) = NMNHDIM_UNUSED
+end if
+
+if ( tzfield%ndimlist(1) == NMNHDIM_BUDGET_LES_LEVEL &
+ .and. tzfield%ndimlist(2) == NMNHDIM_BUDGET_LES_TIME ) then
+ tzfield%ndimlist(3) = NMNHDIM_UNUSED
+ tzfield%ndimlist(4) = NMNHDIM_UNUSED
+ call Les_diachro_common( tpdiafile, tzfield, &
+ reshape( pfield, [ size( pfield, 1 ), size( pfield, 2 ), 1, 1 ] ), &
+ odoavg, odonorm )
+else if ( tzfield%ndimlist(1) == NMNHDIM_BUDGET_LES_TIME &
+ .and. tzfield%ndimlist(2) == NMNHDIM_BUDGET_LES_SV ) then
+ tzfield%ndimlist(4) = tzfield%ndimlist(2)
+ tzfield%ndimlist(2) = tzfield%ndimlist(1)
+ tzfield%ndimlist(1) = NMNHDIM_UNUSED
+ tzfield%ndimlist(3) = NMNHDIM_UNUSED
+ call Les_diachro_common( tpdiafile, tzfield, &
+ reshape( pfield, [ 1, size( pfield, 1 ), 1, size( pfield, 2 ) ] ), &
+ odoavg, odonorm )
+else
+ call Print_msg( NVERB_ERROR, 'IO', 'Les_diachro_2D', &
+ 'ndimlist configuration not yet implemented for ' // Trim( tzfield%cmnhname ) )
+
+end if
+
+end subroutine Les_diachro_2D
+
+!##################################################################################
+subroutine Les_diachro_3D( tpdiafile, tpfield, odoavg, odonorm, pfield, hsuffixes )
+!##################################################################################
+
+use modd_field, only: NMNHDIM_BUDGET_LES_LEVEL, NMNHDIM_BUDGET_LES_MASK, NMNHDIM_BUDGET_LES_SV, &
+ NMNHDIM_BUDGET_LES_TIME, NMNHDIM_UNUSED, &
+ tfield_metadata_base
+use modd_io, only: tfiledata
+
+type(tfiledata), intent(in) :: tpdiafile ! File to write
+type(tfield_metadata_base), intent(in) :: tpfield ! Metadata of field
+logical, intent(in) :: odoavg ! Compute and store time average
+logical, intent(in) :: odonorm ! Compute and store normalized field
+real, dimension(:,:,:), intent(in) :: pfield ! Data array
+character(len=*), dimension(:), optional, intent(in) :: hsuffixes
+
+type(tfield_metadata_base) :: tzfield
+
+tzfield = tpfield
+
+if ( tzfield%ndims /= 3 ) then
+ call Print_msg( NVERB_ERROR, 'IO', 'Les_diachro_3D', 'ndims /= 3 for ' // Trim( tzfield%cmnhname ) )
+ tzfield%ndims = 3
+end if
+
+if ( Any( tzfield%ndimlist(4:) /= NMNHDIM_UNUSED ) ) then
+ call Print_msg( NVERB_ERROR, 'IO', 'Les_diachro_3D', 'unexpected type for some dimensions of ' &
+ // Trim( tzfield%cmnhname ) )
+ tzfield%ndimlist(4:) = NMNHDIM_UNUSED
+end if
+
+if ( tzfield%ndimlist(1) == NMNHDIM_BUDGET_LES_LEVEL &
+ .and. tzfield%ndimlist(2) == NMNHDIM_BUDGET_LES_TIME &
+ .and. tzfield%ndimlist(3) == NMNHDIM_BUDGET_LES_MASK ) then
+ if ( .not. Present( hsuffixes ) ) &
+ call Print_msg( NVERB_ERROR, 'IO', 'Les_diachro_3D', &
+ 'optional dummy argument hsuffixes is needed for tpfield (' // Trim( tzfield%cmnhname ) // ')' )
+
+ if ( Size( hsuffixes ) /= Size( pfield, 3) ) &
+ call Print_msg( NVERB_FATAL, 'IO', 'Les_diachro_3D', 'wrong size for hsuffixes (' // Trim( tzfield%cmnhname ) // ')' )
+
+ tzfield%ndimlist(4) = NMNHDIM_UNUSED
+ call Les_diachro_common( tpdiafile, tzfield, &
+ reshape( pfield, [ size( pfield, 1 ), size( pfield, 2 ), size( pfield, 3 ), 1 ] ), &
+ odoavg, odonorm, hsuffixes )
+else if ( tzfield%ndimlist(1) == NMNHDIM_BUDGET_LES_LEVEL &
+ .and. tzfield%ndimlist(2) == NMNHDIM_BUDGET_LES_TIME &
+ .and. tzfield%ndimlist(3) == NMNHDIM_BUDGET_LES_SV ) then
+ if ( Present( hsuffixes ) ) &
+ call Print_msg( NVERB_ERROR, 'IO', 'Les_diachro_3D', &
+ 'optional dummy argument hsuffixes is not needed for tpfield (' // Trim( tzfield%cmnhname ) // ')' )
+
+ tzfield%ndimlist(4) = tzfield%ndimlist(3)
+ tzfield%ndimlist(3) = NMNHDIM_UNUSED
+ call Les_diachro_common( tpdiafile, tzfield, &
+ reshape( pfield, [ size( pfield, 1 ), size( pfield, 2 ), 1, size( pfield, 3 ) ] ), &
+ odoavg, odonorm )
+else
+ call Print_msg( NVERB_ERROR, 'IO', 'Les_diachro_3D', &
+ 'ndimlist configuration not yet implemented for ' // Trim( tzfield%cmnhname ) )
+end if
+
+end subroutine Les_diachro_3D
+
+!##################################################################################
+subroutine Les_diachro_4D( tpdiafile, tpfield, odoavg, odonorm, pfield, hsuffixes )
+!##################################################################################
+
+use modd_field, only: NMNHDIM_BUDGET_LES_LEVEL, NMNHDIM_BUDGET_LES_MASK, NMNHDIM_BUDGET_LES_SV, &
+ NMNHDIM_BUDGET_LES_TIME, NMNHDIM_BUDGET_TERM, NMNHDIM_UNUSED, &
+ tfield_metadata_base
+use modd_io, only: tfiledata
+
+type(tfiledata), intent(in) :: tpdiafile ! File to write
+type(tfield_metadata_base), intent(in) :: tpfield ! Metadata of field
+logical, intent(in) :: odoavg ! Compute and store time average
+logical, intent(in) :: odonorm ! Compute and store normalized field
+real, dimension(:,:,:,:), intent(in) :: pfield ! Data array
+character(len=*), dimension(:), optional, intent(in) :: hsuffixes
+
+type(tfield_metadata_base) :: tzfield
+
+tzfield = tpfield
+
+if ( tzfield%ndims /= 4 ) then
+ call Print_msg( NVERB_ERROR, 'IO', 'Les_diachro_4D', 'ndims /= 4 for ' // Trim( tzfield%cmnhname ) )
+ tzfield%ndims = 4
+end if
+
+if ( Any( tzfield%ndimlist(5:) /= NMNHDIM_UNUSED ) ) then
+ call Print_msg( NVERB_ERROR, 'IO', 'Les_diachro_4D', 'unexpected type for some dimensions of ' &
+ // Trim( tzfield%cmnhname ) )
+ tzfield%ndimlist(5:) = NMNHDIM_UNUSED
+end if
+
+if ( tzfield%ndimlist(1) == NMNHDIM_BUDGET_LES_LEVEL &
+ .and. tzfield%ndimlist(2) == NMNHDIM_BUDGET_LES_TIME &
+ .and. ( tzfield%ndimlist(3) == NMNHDIM_BUDGET_LES_MASK .or.tzfield%ndimlist(3) == NMNHDIM_BUDGET_TERM ) &
+ .and. tzfield%ndimlist(4) == NMNHDIM_BUDGET_LES_SV ) then
+ if ( .not. Present( hsuffixes ) ) &
+ call Print_msg( NVERB_ERROR, 'IO', 'Les_diachro_4D', &
+ 'optional dummy argument hsuffixes is needed for tpfield (' // Trim( tzfield%cmnhname ) // ')' )
+
+ if ( Size( hsuffixes ) /= Size( pfield, 3) ) &
+ call Print_msg( NVERB_FATAL, 'IO', 'Les_diachro_4D', 'wrong size for hsuffixes (' // Trim( tzfield%cmnhname ) // ')' )
+
+ call Les_diachro_common( tpdiafile, tzfield, pfield, odoavg, odonorm, hsuffixes )
+else
+ call Print_msg( NVERB_ERROR, 'IO', 'Les_diachro_4D', &
+ 'ndimlist configuration not yet implemented for ' // Trim( tzfield%cmnhname ) )
+end if
+
+end subroutine Les_diachro_4D
+
+!#######################################################################################
+subroutine Les_diachro_common( tpdiafile, tpfield, pfield, odoavg, odonorm, hsuffixes )
+!#######################################################################################
+
+use modd_field, only: tfield_metadata_base
+use modd_io, only: tfiledata
+use modd_les, only: nles_current_iinf, nles_current_isup, nles_current_jinf, nles_current_jsup, &
+ nles_levels, xles_current_z
+use modd_parameters, only: XUNDEF
+use modd_type_date, only: date_time
+
+implicit none
+
+type(tfiledata), intent(in) :: tpdiafile ! File to write
+type(tfield_metadata_base), intent(in) :: tpfield
+real, dimension(:,:,:,:), intent(in) :: pfield ! Data array
+logical, intent(in) :: odoavg ! Compute and store time average
+logical, intent(in) :: odonorm ! Compute and store normalized field
+character(len=*), dimension(:), optional, intent(in) :: hsuffixes
+
+character(len=100), dimension(:), allocatable :: ycomment ! Comment string
+character(len=100), dimension(:), allocatable :: ytitle ! Title
+integer :: iavg
+integer :: iles_k ! Number of vertical levels
+integer :: iil, iih, ijl, ijh, ikl, ikh ! Cartesian area relatively to the
+ ! entire domain
+integer :: jk ! Vertical loop counter
+real, dimension(:,:,:), allocatable :: ztrajx ! Localization of the temporal
+real, dimension(:,:,:), allocatable :: ztrajy ! series in x,y and z. remark:
+real, dimension(:,:,:), allocatable :: ztrajz ! x and y are not used for LES
+type(tfield_metadata_base), dimension(:), allocatable :: tzfields
+!------------------------------------------------------------------------------
+
+iles_k = Size( pfield, 1 )
+
+! Initialization of diachro variables for les (z,t) profiles
+Allocate( ztrajx(1, 1, Size( pfield, 4 )) )
+Allocate( ztrajy(1, 1, Size( pfield, 4 )) )
+Allocate( ztrajz(iles_k, 1, Size( pfield, 4 )) )
+Allocate( ycomment(Size( pfield, 3 )) )
+Allocate( ytitle (Size( pfield, 3 )) )
+Allocate( tzfields(Size( pfield, 3 )) )
+
+iil = nles_current_iinf
+iih = nles_current_isup
+ijl = nles_current_jinf
+ijh = nles_current_jsup
+ikl = nles_levels(1)
+ikh = nles_levels(iles_k)
+
+ztrajx(:, :, :) = ( iil + iih ) / 2
+ztrajy(:, :, :) = ( ijl + ijh ) / 2
+do jk = 1, iles_k
+ ztrajz(jk, :, :) = xles_current_z(jk)
+end do
+
+if ( Present( hsuffixes ) ) then
+ if ( Size( hsuffixes ) /= Size( pfield, 3) ) &
+ call Print_msg( NVERB_FATAL, 'IO', 'Les_diachro_common', 'wrong size for hsuffixes (' // Trim( tpfield%cmnhname ) // ')' )
+ ycomment(:) = Trim( tpfield%ccomment(:) ) // hsuffixes(:)
+else
+ ycomment(:) = tpfield%ccomment(:)
+end if
+ call Les_diachro_common_intern( .false., .false. )
+if ( odoavg ) call Les_diachro_common_intern( .true., .false. )
+if ( odonorm ) call Les_diachro_common_intern( .false., .true. )
+if ( odoavg .and. odonorm ) call Les_diachro_common_intern( .true., .true. )
+
+
+contains
+
+!#######################################################################################
+subroutine Les_diachro_common_intern( oavg, onorm )
+!#######################################################################################
+
+use modd_field, only: NMNHDIM_BUDGET_LES_TIME, NMNHDIM_BUDGET_LES_AVG_TIME, NMNHDIM_BUDGET_LES_SV, NMNHDIM_UNUSED
+use modd_les, only: nles_current_times
+
+use mode_write_diachro, only: Write_diachro
+
+logical, intent(in) :: oavg
+logical, intent(in) :: onorm
+
+character(len=10) :: ygroup ! Group title
+integer :: iresp ! Return code
+integer :: ji
+integer :: jp ! Process loop counter
+integer :: jsv ! Scalar loop counter
+logical :: gsv
+real, dimension(:,:,:,:), allocatable :: zfield ! Normalized field
+real, dimension(:,:,:,:,:,:), allocatable :: zwork6 ! Contains physical field
+type(date_time), dimension(:), allocatable :: tzdates
+
+!Reallocate each time necessary because can be reallocated to an other size in Les_time_avg
+Allocate( zfield(Size( pfield, 1 ), Size( pfield, 2 ), Size( pfield, 3 ), Size( pfield, 4 )) )
+Allocate( tzdates( nles_current_times ) )
+
+tzdates(:) = xles_dates(:)
+
+!Copy all fields from tpfield
+tzfields(:) = tpfield
+
+if ( onorm ) then
+ gsv = tzfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_SV
+ call Les_norm_4d( tzfields(1)%cunits, pfield, zfield, gsv )
+
+ ! Normalization of vertical dimension
+ tzfields(:)%cunits = '1'
+ call Les_z_norm( oavg, ztrajz, zfield(:,:,:,:) )
+else
+ zfield(:, :, :, :) = pfield(:, :, :, :)
+end if
+
+! Time average
+iresp = 0
+if ( oavg ) call Les_time_avg_4d( zfield, tzdates, iresp )
+
+if ( oavg ) then
+ if ( onorm ) then
+ ygroup = 'H_' // tpfield%cmnhname
+ else
+ ygroup = 'A_' // tpfield%cmnhname
+ end if
+else
+ if ( onorm ) then
+ ygroup = 'E_' // tpfield%cmnhname
+ else
+ ygroup = tpfield%cmnhname
+ end if
+endif
+
+if ( Present( hsuffixes ) ) then
+ ytitle(:) = ygroup // hsuffixes(:)
+else
+ ytitle(:) = ygroup
+endif
+
+! Write the profile
+if ( iresp == 0 .and. any( zfield /= XUNDEF ) ) then
+ allocate(zwork6(1,1,size(zfield,1),size(zfield,2),size(zfield,4),size(zfield,3)))
+ do jsv = 1, size( zfield, 4 )
+ do jp = 1, size( zfield, 3 )
+ zwork6(1, 1, :, :, jsv, jp) = zfield (:, :, jp, jsv)
+ end do
+ end do
+
+ tzfields(:)%ndimlist(6) = tpfield%ndimlist(3)
+ tzfields(:)%ndimlist(5) = tpfield%ndimlist(4)
+ tzfields(:)%ndimlist(4) = tpfield%ndimlist(2)
+ tzfields(:)%ndimlist(3) = tpfield%ndimlist(1)
+ tzfields(:)%ndimlist(2) = NMNHDIM_UNUSED
+ tzfields(:)%ndimlist(1) = NMNHDIM_UNUSED
+
+ if ( oavg ) then
+ do ji = 1, 6
+ if ( tzfields(1)%ndimlist(ji) == NMNHDIM_BUDGET_LES_TIME ) tzfields(:)%ndimlist(ji) = NMNHDIM_BUDGET_LES_AVG_TIME
+ end do
+ end if
+
+ tzfields(:)%cmnhname = ytitle(:)
+ tzfields(:)%clongname = ytitle(:)
+ tzfields(:)%ccomment = ycomment(:)
+
+ call Write_diachro( tpdiafile, tzfields, ygroup, "SSOL", tzdates, &
+ zwork6, &
+ oicp = .false., ojcp = .false., okcp = .false., &
+ kil = iil, kih = iih, kjl = ijl, kjh = ijh, kkl = ikl, kkh = ikh, &
+ ptrajx = ztrajx, ptrajy = ztrajy, ptrajz = ztrajz )
+end if
+
+!-------------------------------------------------------------------------------
+end subroutine Les_diachro_common_intern
+!-------------------------------------------------------------------------------
+
+!-------------------------------------------------------------------------------
+end subroutine Les_diachro_common
+!------------------------------------------------------------------------------
+
!########################################################
-SUBROUTINE LES_DIACHRO(TPDIAFILE,HGROUP,HCOMMENT,HUNIT,PFIELD,HAVG)
+SUBROUTINE LES_DIACHRO_old(TPDIAFILE,HGROUP,HCOMMENT,HUNIT,PFIELD,HAVG)
!########################################################
IMPLICIT NONE
@@ -845,8 +1239,9 @@ call Les_diachro_gen( tpdiafile, hgroup, [ hcomment ], hunit, &
reshape( pfield, [ size( pfield, 1), size( pfield, 2), 1, 1 ] ), havg )
!-------------------------------------------------------------------------------
-END SUBROUTINE LES_DIACHRO
+END SUBROUTINE LES_DIACHRO_old
!-------------------------------------------------------------------------------
+
!###########################################################
SUBROUTINE LES_DIACHRO_SV(TPDIAFILE,HGROUP,HCOMMENT,HUNIT,PFIELD,HAVG)
!###########################################################
--
GitLab