From 5d826f73ef7ea07aab50bb9c02ca667a06117165 Mon Sep 17 00:00:00 2001
From: Philippe WAUTELET <philippe.wautelet@aero.obs-mip.fr>
Date: Fri, 25 Nov 2022 15:11:32 +0100
Subject: [PATCH] Philippe 25/11/2022: rewrite of aircraft_balloon_evol First
version: some cleaning, deduplication and optimisation not yet done
---
src/MNH/aircraft_balloon.f90 | 38 +-
src/MNH/aircraft_balloon_evol.f90 | 2816 +++++++++++++++++------------
src/MNH/diag.f90 | 20 +-
src/MNH/ini_aircraft.f90 | 44 +-
src/MNH/ini_aircraft_balloon.f90 | 121 +-
src/MNH/ini_balloon.f90 | 6 +-
src/MNH/modd_aircraft_balloon.f90 | 44 +-
src/MNH/modeln.f90 | 14 +-
8 files changed, 1807 insertions(+), 1296 deletions(-)
diff --git a/src/MNH/aircraft_balloon.f90 b/src/MNH/aircraft_balloon.f90
index 02b18e25d..a43f04c53 100644
--- a/src/MNH/aircraft_balloon.f90
+++ b/src/MNH/aircraft_balloon.f90
@@ -9,15 +9,12 @@ MODULE MODI_AIRCRAFT_BALLOON
!
INTERFACE
!
- SUBROUTINE AIRCRAFT_BALLOON(PTSTEP, &
- PXHAT, PYHAT, PZ, &
- PMAP, PLONOR, PLATOR, &
- PU, PV, PW, PP, PTH, PR, PSV, PTKE, &
- PTS, PRHODREF, PCIT, PSEA)
+ SUBROUTINE AIRCRAFT_BALLOON(PTSTEP, PZ, &
+ PMAP, PLONOR, PLATOR, &
+ PU, PV, PW, PP, PTH, PR, PSV, PTKE, &
+ PTS, PRHODREF, PCIT, PSEA )
!
REAL, INTENT(IN) :: PTSTEP ! time step
-REAL, DIMENSION(:), INTENT(IN) :: PXHAT ! x coordinate
-REAL, DIMENSION(:), INTENT(IN) :: PYHAT ! y coordinate
REAL, DIMENSION(:,:,:), INTENT(IN) :: PZ ! z
REAL, DIMENSION(:,:), INTENT(IN) :: PMAP ! map factor
REAL, INTENT(IN) :: PLONOR ! origine longitude
@@ -31,11 +28,9 @@ REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PR ! water mixing ratios
REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PSV ! Scalar variables
REAL, DIMENSION(:,:,:), INTENT(IN) :: PTKE ! turbulent kinetic energy
REAL, DIMENSION(:,:), INTENT(IN) :: PTS ! surface temperature
-! ++ OC
REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODREF ! dry air density of the reference state
REAL, DIMENSION(:,:,:), INTENT(IN) :: PCIT ! pristine ice concentration
-REAL, DIMENSION(:,:),OPTIONAL,INTENT(IN) :: PSEA
-! -- OC
+REAL, DIMENSION(:,:), OPTIONAL, INTENT(IN) :: PSEA
!
!-------------------------------------------------------------------------------
!
@@ -52,13 +47,12 @@ END INTERFACE
!
END MODULE MODI_AIRCRAFT_BALLOON
!
-! ###################################################################
- SUBROUTINE AIRCRAFT_BALLOON(PTSTEP, &
- PXHAT, PYHAT, PZ, &
- PMAP, PLONOR, PLATOR, &
- PU, PV, PW, PP, PTH, PR, PSV, PTKE, &
- PTS, PRHODREF, PCIT, PSEA)
-! ###################################################################
+! #################################################################
+ SUBROUTINE AIRCRAFT_BALLOON(PTSTEP, PZ, &
+ PMAP, PLONOR, PLATOR, &
+ PU, PV, PW, PP, PTH, PR, PSV, PTKE, &
+ PTS, PRHODREF, PCIT, PSEA )
+! #################################################################
!
!
!!**** *AIRCRAFT_BALLOON* - monitor for balloons and aircrafts
@@ -101,8 +95,8 @@ END MODULE MODI_AIRCRAFT_BALLOON
USE MODD_AIRCRAFT_BALLOON
!
USE MODD_TURB_FLUX_AIRCRAFT_BALLOON
-USE MODI_AIRCRAFT_BALLOON_EVOL
!
+USE MODE_AIRCRAFT_BALLOON_EVOL, ONLY: AIRCRAFT_BALLOON_EVOL
USE MODE_ll
!
!
@@ -113,8 +107,6 @@ IMPLICIT NONE
!
!
REAL, INTENT(IN) :: PTSTEP ! time step
-REAL, DIMENSION(:), INTENT(IN) :: PXHAT ! x coordinate
-REAL, DIMENSION(:), INTENT(IN) :: PYHAT ! y coordinate
REAL, DIMENSION(:,:,:), INTENT(IN) :: PZ ! z array
REAL, DIMENSION(:,:), INTENT(IN) :: PMAP ! map factor
REAL, INTENT(IN) :: PLONOR ! origine longitude
@@ -130,7 +122,7 @@ REAL, DIMENSION(:,:,:), INTENT(IN) :: PTKE ! turbulent kinetic energy
REAL, DIMENSION(:,:), INTENT(IN) :: PTS ! surface temperature
REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODREF ! dry air density of the reference state
REAL, DIMENSION(:,:,:), INTENT(IN) :: PCIT ! pristine ice concentration
-REAL, DIMENSION(:,:), INTENT(IN) :: PSEA
+REAL, DIMENSION(:,:), OPTIONAL, INTENT(IN) :: PSEA
!
!-------------------------------------------------------------------------------
!
@@ -146,13 +138,13 @@ IF(.NOT. ALLOCATED(XSVW_FLUX)) &
ALLOCATE(XSVW_FLUX(SIZE(PSV,1),SIZE(PSV,2),SIZE(PSV,3),SIZE(PSV,4)))
!
DO JI = 1, NBALLOONS
- CALL AIRCRAFT_BALLOON_EVOL( PTSTEP, PXHAT, PYHAT, PZ, PMAP, PLONOR, PLATOR, &
+ CALL AIRCRAFT_BALLOON_EVOL( PTSTEP, PZ, PMAP, PLONOR, PLATOR, &
PU, PV, PW, PP, PTH, PR, PSV, PTKE, PTS, PRHODREF, PCIT, &
TBALLOONS(JI), PSEA )
END DO
!
DO JI = 1, NAIRCRAFTS
- CALL AIRCRAFT_BALLOON_EVOL( PTSTEP, PXHAT, PYHAT, PZ, PMAP, PLONOR, PLATOR, &
+ CALL AIRCRAFT_BALLOON_EVOL( PTSTEP, PZ, PMAP, PLONOR, PLATOR, &
PU, PV, PW, PP, PTH, PR, PSV, PTKE, PTS, PRHODREF, PCIT, &
TAIRCRAFTS(JI), PSEA )
END DO
diff --git a/src/MNH/aircraft_balloon_evol.f90 b/src/MNH/aircraft_balloon_evol.f90
index 6a81850ba..5b52ea22f 100644
--- a/src/MNH/aircraft_balloon_evol.f90
+++ b/src/MNH/aircraft_balloon_evol.f90
@@ -4,59 +4,27 @@
!MNH_LIC for details. version 1.
!-----------------------------------------------------------------
! ##########################
-MODULE MODI_AIRCRAFT_BALLOON_EVOL
+MODULE MODE_AIRCRAFT_BALLOON_EVOL
! ##########################
-!
-INTERFACE
-!
- SUBROUTINE AIRCRAFT_BALLOON_EVOL(PTSTEP, &
- PXHAT, PYHAT, PZ, &
- PMAP, PLONOR, PLATOR, &
- PU, PV, PW, PP, PTH, PR, PSV, PTKE, &
- PTS, PRHODREF, PCIT,TPFLYER, PSEA )
-!
-USE MODD_AIRCRAFT_BALLOON
-!
-REAL, INTENT(IN) :: PTSTEP ! time step
-REAL, DIMENSION(:), INTENT(IN) :: PXHAT ! x coordinate
-REAL, DIMENSION(:), INTENT(IN) :: PYHAT ! y coordinate
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PZ ! z array
-REAL, DIMENSION(:,:), INTENT(IN) :: PMAP ! map factor
-REAL, INTENT(IN) :: PLONOR ! origine longitude
-REAL, INTENT(IN) :: PLATOR ! origine latitude
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PU ! horizontal wind X component
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PV ! horizontal wind Y component
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PW ! vertical wind
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PP ! pressure
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PTH ! potential temperature
-REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PR ! water mixing ratios
-REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PSV ! Scalar variables
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PTKE ! turbulent kinetic energy
-REAL, DIMENSION(:,:), INTENT(IN) :: PTS ! surface temperature
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODREF ! dry air density of the reference state
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PCIT ! pristine ice concentration
-!
-CLASS(TFLYERDATA), INTENT(INOUT) :: TPFLYER! balloon/aircraft
-REAL, DIMENSION(:,:), OPTIONAL, INTENT(IN) :: PSEA
-!
-!-------------------------------------------------------------------------------
-!
-END SUBROUTINE AIRCRAFT_BALLOON_EVOL
-!
-END INTERFACE
-!
-END MODULE MODI_AIRCRAFT_BALLOON_EVOL
-!
+
+USE MODE_MSG
+
+IMPLICIT NONE
+
+PRIVATE
+
+PUBLIC :: AIRCRAFT_BALLOON_EVOL
+
+CONTAINS
! ########################################################
SUBROUTINE AIRCRAFT_BALLOON_EVOL(PTSTEP, &
- PXHAT, PYHAT, PZ, &
- PMAP, PLONOR, PLATOR, &
+ PZ, PMAP, PLONOR, PLATOR, &
PU, PV, PW, PP, PTH, PR, PSV, PTKE, &
PTS, PRHODREF, PCIT,TPFLYER, PSEA )
! ########################################################
!
!
-!!**** *AIRCRAFT_BALLOON_EVOL* - (advects and) stores
+!!**** *AIRCRAFT_BALLOON_EVOL* - (advects and) stores
!! balloons/aircrafts in the model
!!
!! PURPOSE
@@ -65,7 +33,7 @@ END MODULE MODI_AIRCRAFT_BALLOON_EVOL
!
!!** METHOD
!! ------
-!!
+!!
!! 1) All the balloons are tested. If the current balloon is
!! a) in the current model
!! b) not crashed
@@ -86,10 +54,10 @@ END MODULE MODI_AIRCRAFT_BALLOON_EVOL
!! a) iso-density balloons are advected following horizontal wind.
!! the slope of the iso-density surfaces is neglected.
!! b) radio-sounding balloons are advected according to all wind velocities.
-!! the vertical ascent speed is added to the vertical wind speed.
+!! the vertical ascent speed is added to the vertical wind speed.
!! c) Constant Volume balloons are advected according to all wind velocities.
!! the vertical ascent speed is computed using the balloon equation
-!!
+!!
!!
!! EXTERNAL
!! --------
@@ -115,7 +83,7 @@ END MODULE MODI_AIRCRAFT_BALLOON_EVOL
!! April, 2014 (C.Lac) allow RARE calculation only if CCLOUD=ICE3
!! May, 2014 (O.Caumont) modify RARE for hydrometeors containing ice
!! add bright band calculation for RARE
-!! Feb, 2015 (C.Lac) Correction to prevent aircraft crash
+!! Feb, 2015 (C.Lac) Correction to prevent aircraft crash
!! July, 2015 (O.Nuissier/F.Duffourg) Add microphysics diagnostic for
!! aircraft, ballon and profiler
!! October, 2016 (G.DELAUTIER) LIMA
@@ -129,7 +97,7 @@ END MODULE MODI_AIRCRAFT_BALLOON_EVOL
! -do not use PMAP if cartesian domain
! P. Wautelet 06/2022: reorganize flyers
!! --------------------------------------------------------------------------
-!
+!
!* 0. DECLARATIONS
! ------------
!
@@ -138,7 +106,8 @@ USE MODD_CONF
USE MODD_CST
USE MODD_DIAG_IN_RUN
USE MODD_GRID
-USE MODD_GRID_n, ONLY: XXHATM, XYHATM
+USE MODD_GRID_n
+USE MODD_IO, ONLY: ISP
USE MODD_LUNIT_n, ONLY: TLUOUT
USE MODD_NESTING
USE MODD_NSV, ONLY : NSV_LIMA_NI,NSV_LIMA_NR,NSV_LIMA_NC
@@ -167,7 +136,7 @@ USE MODD_RAIN_ICE_DESCR, ONLY: XALPHAR_I=>XALPHAR,XNUR_I=>XNUR,XLBEXR_I=>XLBEX
XLBI_I=>XLBI,XAI_I=>XAI,XBI_I=>XBI,XC_I_I=>XC_I,&
XRTMIN_I=>XRTMIN,XCONC_LAND,XCONC_SEA
USE MODD_REF_n, ONLY: XRHODREF
-USE MODD_TIME, only: tdtexp
+USE MODD_TIME, only: TDTSEG
USE MODD_TIME_n, only: tdtcur
USE MODD_TURB_FLUX_AIRCRAFT_BALLOON
!
@@ -176,7 +145,7 @@ USE MODE_FGAU, ONLY: GAULAG
USE MODE_FSCATTER, ONLY: QEPSW,QEPSI,BHMIE,MOMG,MG
USE MODE_GRIDPROJ
USE MODE_ll
-USE MODE_MSG
+USE MODE_POSITION_TOOLS, ONLY: FIND_PROCESS_AND_MODEL_FROM_XY_POS
USE MODE_STATPROF_TOOLS, ONLY: STATPROF_INSTANT
!
USE MODI_GAMMA, ONLY: GAMMA
@@ -189,8 +158,6 @@ IMPLICIT NONE
!
!
REAL, INTENT(IN) :: PTSTEP ! time step
-REAL, DIMENSION(:), INTENT(IN) :: PXHAT ! x coordinate
-REAL, DIMENSION(:), INTENT(IN) :: PYHAT ! y coordinate
REAL, DIMENSION(:,:,:), INTENT(IN) :: PZ ! z array
REAL, DIMENSION(:,:), INTENT(IN) :: PMAP ! map factor
REAL, INTENT(IN) :: PLONOR ! origine longitude
@@ -208,7 +175,7 @@ REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODREF ! dry air density of the re
REAL, DIMENSION(:,:,:), INTENT(IN) :: PCIT ! pristine ice concentration
!
CLASS(TFLYERDATA), INTENT(INOUT) :: TPFLYER! balloon/aircraft
-REAL, DIMENSION(:,:), INTENT(IN) :: PSEA
+REAL, DIMENSION(:,:), OPTIONAL, INTENT(IN) :: PSEA
!
!-------------------------------------------------------------------------------
!
@@ -218,7 +185,7 @@ REAL, DIMENSION(:,:), INTENT(IN) :: PSEA
INTEGER :: IMI ! model index
REAL :: ZTHIS_PROC ! 1 if balloon is currently treated by this proc., else 0
!
-INTEGER :: IIB ! current processor domain sizes
+INTEGER :: IIB ! current process domain sizes
INTEGER :: IJB
INTEGER :: IIE
INTEGER :: IJE
@@ -241,7 +208,7 @@ REAL, DIMENSION(2,2,SIZE(PTH,3)) :: ZTEMP ! temperature
REAL, DIMENSION(2,2,SIZE(PTH,3)) :: ZEXN ! Exner function
REAL, DIMENSION(2,2,SIZE(PTH,3)) :: ZRHO ! air density
REAL :: ZFLYER_EXN ! balloon/aircraft Exner func.
-REAL, DIMENSION(2,2,SIZE(PTH,3)) :: ZTHW_FLUX !
+REAL, DIMENSION(2,2,SIZE(PTH,3)) :: ZTHW_FLUX !
REAL, DIMENSION(2,2,SIZE(PTH,3)) :: ZRCW_FLUX !
REAL, DIMENSION(2,2,SIZE(PSV,3),SIZE(PSV,4)) :: ZSVW_FLUX
!
@@ -249,45 +216,45 @@ REAL :: ZTDIST ! time until launch (sec)
LOGICAL :: GLAUNCH ! launch/takeoff is effective at this time-step (if true)
LOGICAL :: GSTORE ! storage occurs at this time step
!
-INTEGER :: II ! mass balloon position (x index)
-INTEGER :: IJ ! mass balloon position (y index)
-INTEGER :: IU ! U flux point balloon position (x index)
-INTEGER :: IV ! V flux point balloon position (y index)
-INTEGER :: IDU ! difference between IU and II
-INTEGER :: IDV ! difference between IV and IJ
-!
-INTEGER :: IK00 ! balloon position for II , IJ
-INTEGER :: IK01 ! balloon position for II , IJ+1
-INTEGER :: IK10 ! balloon position for II+1, IJ
-INTEGER :: IK11 ! balloon position for II+1, IJ+1
-INTEGER :: IU00 ! balloon position for IU , IJ
-INTEGER :: IU01 ! balloon position for IU , IJ+1
-INTEGER :: IU10 ! balloon position for IU+1, IJ
-INTEGER :: IU11 ! balloon position for IU+1, IJ+1
-INTEGER :: IV00 ! balloon position for II , IV
-INTEGER :: IV01 ! balloon position for II , IV+1
-INTEGER :: IV10 ! balloon position for II+1, IV
-INTEGER :: IV11 ! balloon position for II+1, IV+1
+INTEGER :: II_M ! mass balloon position (x index)
+INTEGER :: IJ_M ! mass balloon position (y index)
+INTEGER :: II_U ! U flux point balloon position (x index)
+INTEGER :: IJ_V ! V flux point balloon position (y index)
+INTEGER :: IDU ! difference between II_U and II_M
+INTEGER :: IDV ! difference between IJ_V and IJ_M
+!
+INTEGER :: IK00 ! balloon position for II_M , IJ_M
+INTEGER :: IK01 ! balloon position for II_M , IJ_M+1
+INTEGER :: IK10 ! balloon position for II_M+1, IJ_M
+INTEGER :: IK11 ! balloon position for II_M+1, IJ_M+1
+INTEGER :: IU00 ! balloon position for II_U , IJ_M
+INTEGER :: IU01 ! balloon position for II_U , IJ_M+1
+INTEGER :: IU10 ! balloon position for II_U+1, IJ_M
+INTEGER :: IU11 ! balloon position for II_U+1, IJ_M+1
+INTEGER :: IV00 ! balloon position for II_M , IJ_V
+INTEGER :: IV01 ! balloon position for II_M , IJ_V+1
+INTEGER :: IV10 ! balloon position for II_M+1, IJ_V
+INTEGER :: IV11 ! balloon position for II_M+1, IJ_V+1
!
REAL :: ZXCOEF ! X direction interpolation coefficient
REAL :: ZUCOEF ! X direction interpolation coefficient (for U)
REAL :: ZYCOEF ! Y direction interpolation coefficient
REAL :: ZVCOEF ! Y direction interpolation coefficient (for V)
!
-REAL :: ZZCOEF00 ! Z direction interpolation coefficient for II , IJ
-REAL :: ZZCOEF01 ! Z direction interpolation coefficient for II , IJ+1
-REAL :: ZZCOEF10 ! Z direction interpolation coefficient for II+1, IJ
-REAL :: ZZCOEF11 ! Z direction interpolation coefficient for II+1, IJ+1
-REAL :: ZUCOEF00 ! Z direction interpolation coefficient for IU , IJ
-REAL :: ZUCOEF01 ! Z direction interpolation coefficient for IU , IJ+1
-REAL :: ZUCOEF10 ! Z direction interpolation coefficient for IU+1, IJ
-REAL :: ZUCOEF11 ! Z direction interpolation coefficient for IU+1, IJ+1
-REAL :: ZVCOEF00 ! Z direction interpolation coefficient for II , IV
-REAL :: ZVCOEF01 ! Z direction interpolation coefficient for II , IV+1
-REAL :: ZVCOEF10 ! Z direction interpolation coefficient for II+1, IV
-REAL :: ZVCOEF11 ! Z direction interpolation coefficient for II+1, IV+1
-!
-INTEGER :: IN ! time index
+REAL :: ZZCOEF00 ! Z direction interpolation coefficient for II_M , IJ_M
+REAL :: ZZCOEF01 ! Z direction interpolation coefficient for II_M , IJ_M+1
+REAL :: ZZCOEF10 ! Z direction interpolation coefficient for II_M+1, IJ_M
+REAL :: ZZCOEF11 ! Z direction interpolation coefficient for II_M+1, IJ_M+1
+REAL :: ZUCOEF00 ! Z direction interpolation coefficient for II_U , IJ_M
+REAL :: ZUCOEF01 ! Z direction interpolation coefficient for II_U , IJ_M+1
+REAL :: ZUCOEF10 ! Z direction interpolation coefficient for II_U+1, IJ_M
+REAL :: ZUCOEF11 ! Z direction interpolation coefficient for II_U+1, IJ_M+1
+REAL :: ZVCOEF00 ! Z direction interpolation coefficient for II_M , IJ_V
+REAL :: ZVCOEF01 ! Z direction interpolation coefficient for II_M , IJ_V+1
+REAL :: ZVCOEF10 ! Z direction interpolation coefficient for II_M+1, IJ_V
+REAL :: ZVCOEF11 ! Z direction interpolation coefficient for II_M+1, IJ_V+1
+!
+INTEGER :: ISTORE ! time index for storage
INTEGER :: JLOOP,JLOOP2 ! loop counter
!
REAL :: ZU_BAL ! horizontal wind speed at balloon location (along x)
@@ -308,7 +275,7 @@ REAL, DIMENSION(SIZE(PR,3)) :: ZTEMPZ! vertical profile of temperature
REAL, DIMENSION(SIZE(PR,3)) :: ZRHODREFZ ! vertical profile of dry air density of the reference state
REAL, DIMENSION(SIZE(PR,3)) :: ZCIT ! pristine ice concentration
REAL, DIMENSION(SIZE(PR,3)) :: ZCCI,ZCCR,ZCCC ! ICE,RAIN CLOUD concentration (LIMA)
-REAL, DIMENSION(SIZE(PR,1),SIZE(PR,2),SIZE(PR,3)) :: ZR
+REAL, DIMENSION(SIZE(PR,1),SIZE(PR,2),SIZE(PR,3)) :: ZR
REAL, DIMENSION(SIZE(PR,3),SIZE(PR,4)+1) :: ZRZ ! vertical profile of hydrometeor mixing ratios
REAL :: ZA,ZB,ZCC,ZCX,ZALPHA,ZNU,ZLB,ZLBEX,ZRHOHYD ! generic microphysical parameters
INTEGER :: JJ ! loop counter for quadrature
@@ -324,318 +291,908 @@ REAL :: ZFPW ! weight for mixed-phase reflectivity
REAL,DIMENSION(:),ALLOCATABLE :: ZX,ZW ! Gauss-Laguerre points and weights
REAL,DIMENSION(:),ALLOCATABLE :: ZRTMIN ! local values for XRTMIN
LOGICAL :: GCALC
-!----------------------------------------------------------------------------
-!
-!* 1. PRELIMINARIES
-! -------------
-!
-IF(.NOT. ALLOCATED(XTHW_FLUX)) &
-ALLOCATE(XTHW_FLUX(SIZE(PTH,1),SIZE(PTH,2),SIZE(PTH,3)))
-IF(.NOT. ALLOCATED(XRCW_FLUX)) &
-ALLOCATE(XRCW_FLUX(SIZE(PTH,1),SIZE(PTH,2),SIZE(PTH,3)))
-IF(.NOT. ALLOCATED(XSVW_FLUX)) &
-ALLOCATE(XSVW_FLUX(SIZE(PSV,1),SIZE(PSV,2),SIZE(PSV,3),SIZE(PSV,4)))
-ILUOUT = TLUOUT%NLU
-!
-ZR = 0.
-GSTORE = .FALSE.
-!
-!* 1.0 initialization of processor test
-! --------------------------------
-!
-ZTHIS_PROC=0.
-!
-!
-!* 1.1 test on the model
-! -----------------
-!
-CALL GET_MODEL_NUMBER_ll (IMI)
-!
-!
-IF ( TPFLYER%CMODEL /= 'FIX' .AND. COUNT( NDAD(:) == IMI ) /= 0 &
- .AND. ( TPFLYER%NMODEL == IMI .OR. NDAD(TPFLYER%NMODEL) == IMI ) &
- .AND. TPFLYER%XX_CUR /= XUNDEF .AND. TPFLYER%XY_CUR /= XUNDEF &
- .AND. TPFLYER%LFLY .AND. .NOT. TPFLYER%LCRASH &
- .AND. CPROGRAM == 'MESONH' ) THEN
- CALL FLYER_CHANGE_MODEL( IMI )
-ENDIF
-!
-IF ( TPFLYER%NMODEL /= IMI ) RETURN
!
+INTEGER :: IRANK
+INTEGER :: IMODEL
+INTEGER :: IMODEL_OLD
+REAL :: ZX_OLD, ZY_OLD
+REAL :: ZDELTATIME
+REAL :: ZTSTEP
+REAL :: ZDIVTMP
+
!----------------------------------------------------------------------------
-!
-!* 2. PRELIMINARIES-2
-! -------------
-!
-!* 2.1 Indices
-! -------
-!
-CALL GET_INDICE_ll (IIB,IJB,IIE,IJE)
-IKB = 1 + JPVEXT
-IKE = SIZE(PZ,3) - JPVEXT
IKU = SIZE(PZ,3)
-!
-!
-!* 2.2 Interpolations of model variables to mass points
-! ------------------------------------------------
-!
-IIU=SIZE(PXHAT)
-IJU=SIZE(PYHAT)
-!
-!----------------------------------------------------------------------------
-!
-!* 2.3 Compute time until launch by comparison of dates and times
-! ----------------------------------------------------------
-!
-CALL DATETIME_DISTANCE( TPFLYER%TLAUNCH, TDTCUR, ZTDIST )
-!
-!* 3. LAUNCH
-! ------
-!
-GLAUNCH = .FALSE.
-!
-!
-IF ( .NOT. TPFLYER%LFLY ) THEN
-!
-!* 3.2 launch/takeoff is effective
-! ---------------------------
-!
- IF (ZTDIST >= - PTSTEP ) THEN
- SELECT TYPE ( TPFLYER )
- CLASS IS ( TAIRCRAFTDATA)
-!
-!* 3.2.1 Determination of flight segment
-! -------------------------------
-!
- TPFLYER%NSEGCURN = 1
- IL = TPFLYER%NSEGCURN
- !
- TPFLYER%XSEGCURT = ZTDIST
- !
- DO WHILE (TPFLYER%XSEGCURT>TPFLYER%XSEGTIME(IL) .AND. IL <= TPFLYER%NSEG)
- TPFLYER%NSEGCURN = TPFLYER%NSEGCURN + 1
- IL = TPFLYER%NSEGCURN
- TPFLYER%XSEGCURT = TPFLYER%XSEGCURT - TPFLYER%XSEGTIME(IL-1)
- IF (IL>TPFLYER%NSEG) EXIT
+
+CALL GET_MODEL_NUMBER_ll(IMI)
+
+SELECT TYPE ( TPFLYER )
+ CLASS IS ( TAIRCRAFTDATA)
+
+ !Do the positioning only if model 1 (data will be available to others after)
+ MODEL1: IF ( IMI == 1 ) THEN
+ !Do we have to store aircraft data?
+ CALL STATPROF_INSTANT( TPFLYER%TFLYER_TIME, ISTORE )
+ IF ( ISTORE < 1 ) THEN
+ !No profiler storage at this time step
+ TPFLYER%LSTORE = .FALSE.
+ RETURN
+ ELSE
+ TPFLYER%LSTORE = .TRUE.
+ END IF
+
+ ! Is the aircraft in flight ?
+ IF ( TDTCUR >= TPFLYER%TLAUNCH .AND. TDTCUR <= TPFLYER%TLAND ) THEN
+ TPFLYER%LFLY = .TRUE.
+
+ ! Find the flight segment
+ ZTDIST = TDTCUR - TPFLYER%TLAUNCH
+ IL = TPFLYER%NPOSCUR
+ DO WHILE ( ZTDIST > TPFLYER%XPOSTIME(IL+1) )
+ IL = IL + 1
+ IF ( IL > TPFLYER%NPOS-1 ) THEN
+ !Security (should not happen)
+ IL = TPFLYER%NPOS-1
+ EXIT
+ END IF
END DO
- !
- !* end of flight
- !
- IF (IL > TPFLYER%NSEG) THEN
- TPFLYER%LFLY = .FALSE.
+ TPFLYER%NPOSCUR = IL
+
+ ! Compute the current position
+ ZSEG_FRAC = ( ZTDIST - TPFLYER%XPOSTIME(IL) ) / ( TPFLYER%XPOSTIME(IL+1) - TPFLYER%XPOSTIME(IL) )
+
+ TPFLYER%XX_CUR = (1.-ZSEG_FRAC) * TPFLYER%XPOSX(IL ) &
+ + ZSEG_FRAC * TPFLYER%XPOSX(IL+1)
+ TPFLYER%XY_CUR = (1.-ZSEG_FRAC) * TPFLYER%XPOSY(IL ) &
+ + ZSEG_FRAC * TPFLYER%XPOSY(IL+1)
+ IF (TPFLYER%LALTDEF) THEN
+ TPFLYER%XP_CUR = (1.-ZSEG_FRAC) * TPFLYER%XPOSP(IL ) &
+ + ZSEG_FRAC * TPFLYER%XPOSP(IL+1)
ELSE
- TPFLYER%LFLY = .TRUE.
- GLAUNCH = .TRUE.
- TPFLYER%LCRASH =.FALSE.
- IF (ZTDIST <= PTSTEP ) THEN
- WRITE(ILUOUT,*) '-------------------------------------------------------------------'
- WRITE(ILUOUT,*) 'Aircraft ',TPFLYER%CTITLE,' takes off the ', &
- TDTCUR%nday,'/',TDTCUR%nmonth,'/', &
- TDTCUR%nyear,' at ',NINT(TDTCUR%xtime),' sec.'
- WRITE(ILUOUT,*) '-------------------------------------------------------------------'
- ENDIF
- ENDIF
-
- CLASS IS ( TBALLOONDATA)
- IF (ZTDIST <= PTSTEP ) THEN
- TPFLYER%LFLY = .TRUE.
- GLAUNCH = .TRUE.
- WRITE(ILUOUT,*) '-------------------------------------------------------------------'
- WRITE(ILUOUT,*) 'Balloon ',TPFLYER%CTITLE,' is launched the ', &
- TDTCUR%nday,'/',TDTCUR%nmonth,'/', &
- TDTCUR%nyear,' at ',NINT(TDTCUR%xtime),' sec.'
- WRITE(ILUOUT,*) '-------------------------------------------------------------------'
+ TPFLYER%XZ_CUR = (1.-ZSEG_FRAC) * TPFLYER%XPOSZ(IL ) &
+ + ZSEG_FRAC * TPFLYER%XPOSZ(IL +1 )
END IF
- CLASS DEFAULT
- CALL PRINT_MSG( NVERB_ERROR, 'GEN', 'AIRCRAFT_BALLOON_EVOL', 'unknown type for TPFLYER', OLOCAL = .TRUE. )
-
- END SELECT
-!
-!* 3.3 Initial horizontal positions
-! ----------------------------
-!
- SELECT TYPE ( TPFLYER )
- CLASS IS ( TBALLOONDATA)
- IF ( TPFLYER%CTYPE == 'RADIOS' .OR. TPFLYER%CTYPE == 'ISODEN' .OR. TPFLYER%CTYPE == 'CVBALL' ) THEN
- TPFLYER%XX_CUR = TPFLYER%XXLAUNCH
- TPFLYER%XY_CUR = TPFLYER%XYLAUNCH
+ ! Get rank of the process where the aircraft is and the model number
+ IF ( TPFLYER%CMODEL == 'FIX' ) THEN
+ IMODEL = TPFLYER%NMODEL
+ ELSE
+ IMODEL = 0
+ END IF
+ CALL FIND_PROCESS_AND_MODEL_FROM_XY_POS( TPFLYER%XX_CUR, TPFLYER%XY_CUR, IRANK, IMODEL )
+ IF ( IRANK < 1 ) THEN
+ TPFLYER%NMODEL = 1 !Set to 1 because it is always a valid model (to prevent crash if NDAD(TPFLYER%NMODEL) )
+ TPFLYER%LCRASH = .TRUE.
+ TPFLYER%NCRASH = NCRASH_OUT_HORIZ
+ TPFLYER%LFLY = .FALSE.
+ ELSE
+ TPFLYER%NMODEL = IMODEL
+ TPFLYER%LCRASH = .FALSE.
+ TPFLYER%NCRASH = NCRASH_NO
+ TPFLYER%NRANK_CUR = IRANK
END IF
+ ELSE
+ TPFLYER%LFLY = .FALSE.
- CLASS IS ( TAIRCRAFTDATA)
-!
-!
-!* 3.3.2 Determination of initial position
-! -----------------------------
-!
- IF (TPFLYER%LFLY) THEN
- ZSEG_FRAC = TPFLYER%XSEGCURT / TPFLYER%XSEGTIME(IL)
+ END IF
+ END IF MODEL1
+
+ IF ( IMI == TPFLYER%NMODEL .AND. TPFLYER%LFLY ) THEN
+ ! Check if it is the right moment to store data
+ IF ( TPFLYER%LSTORE ) THEN
+ ISTORE = TPFLYER%TFLYER_TIME%N_CUR
+ IF ( ABS( TDTCUR - TPFLYER%TFLYER_TIME%TPDATES(ISTORE) ) < 1e-10 ) THEN
+
+ ISOWNER: IF ( TPFLYER%NRANK_CUR == ISP ) THEN
+ ZTHIS_PROC = 1.
+ !
+ !* 2. PRELIMINARIES-2
+ ! -------------
+ !
+ !* 2.1 Indices
+ ! -------
+ !
+ CALL GET_INDICE_ll (IIB,IJB,IIE,IJE)
+
+ IKB = 1 + JPVEXT
+ IKE = SIZE(PZ,3) - JPVEXT
+ !
+ !
+ !* 2.2 Interpolations of model variables to mass points
+ ! ------------------------------------------------
+ !
+ IIU=SIZE(XXHAT)
+ IJU=SIZE(XYHAT)
+ ! X position
+ II_U = COUNT( XXHAT (:) <= TPFLYER%XX_CUR )
+ II_M = COUNT( XXHATM(:) <= TPFLYER%XX_CUR )
+
+ ! Y position
+ IJ_V=COUNT( XYHAT (:)<=TPFLYER%XY_CUR )
+ IJ_M=COUNT( XYHATM(:)<=TPFLYER%XY_CUR )
+ !
+ !* 4.5 Interpolations of model variables to mass points
+ ! ------------------------------------------------
+ !
+ ZZM(:,:,1:IKU-1)=0.5 *PZ(II_M :II_M+1,IJ_M :IJ_M+1,1:IKU-1)+0.5 *PZ(II_M :II_M+1,IJ_M :IJ_M+1,2:IKU )
+ ZZM(:,:, IKU )=1.5 *PZ(II_M :II_M+1,IJ_M :IJ_M+1, IKU-1)-0.5 *PZ(II_M :II_M+1,IJ_M :IJ_M+1, IKU-2)
+
+ IDU = II_U - II_M
+ ZZU(:,:,1:IKU-1)=0.25*PZ(IDU+II_M-1:IDU+II_M, IJ_M :IJ_M+1,1:IKU-1)+0.25*PZ(IDU+II_M-1:IDU+II_M ,IJ_M :IJ_M+1,2:IKU ) &
+ +0.25*PZ(IDU+II_M :IDU+II_M+1,IJ_M :IJ_M+1,1:IKU-1)+0.25*PZ(IDU+II_M :IDU+II_M+1,IJ_M :IJ_M+1,2:IKU )
+ ZZU(:,:, IKU )=0.75*PZ(IDU+II_M-1:IDU+II_M ,IJ_M :IJ_M+1, IKU-1)-0.25*PZ(IDU+II_M-1:IDU+II_M ,IJ_M :IJ_M+1, IKU-2) &
+ +0.75*PZ(IDU+II_M :IDU+II_M+1,IJ_M :IJ_M+1, IKU-1)-0.25*PZ(IDU+II_M :IDU+II_M+1,IJ_M :IJ_M+1, IKU-2)
+
+ IDV = IJ_V - IJ_M
+ ZZV(:,:,1:IKU-1)=0.25*PZ(II_M :II_M+1,IDV+IJ_M-1:IDV+IJ_M ,1:IKU-1)+0.25*PZ(II_M :II_M+1,IDV+IJ_M-1:IDV+IJ_M ,2:IKU ) &
+ +0.25*PZ(II_M :II_M+1,IDV+IJ_M :IDV+IJ_M+1,1:IKU-1)+0.25*PZ(II_M :II_M+1,IDV+IJ_M :IDV+IJ_M+1,2:IKU )
+ ZZV(:,:, IKU )=0.75*PZ(II_M :II_M+1,IDV+IJ_M-1:IDV+IJ_M , IKU-1)-0.25*PZ(II_M :II_M+1,IDV+IJ_M-1:IDV+IJ_M , IKU-2) &
+ +0.75*PZ(II_M :II_M+1,IDV+IJ_M :IDV+IJ_M+1, IKU-1)-0.25*PZ(II_M :II_M+1,IDV+IJ_M :IDV+IJ_M+1, IKU-2)
+
+ ZWM(:,:,1:IKU-1)=0.5*PW(II_M:II_M+1,IJ_M:IJ_M+1,1:IKU-1)+0.5*PW(II_M:II_M+1,IJ_M:IJ_M+1,2:IKU )
+ ZWM(:,:, IKU )=1.5*PW(II_M:II_M+1,IJ_M:IJ_M+1, IKU-1)-0.5*PW(II_M:II_M+1,IJ_M:IJ_M+1, IKU-2)
+ !
+ !----------------------------------------------------------------------------
+ !
+ !* 5. BALLOON/AIRCRAFT VERTICAL POSITION
+ ! ----------------------------------
+ !
+ !
+ !* 5.1 Density
+ ! -------
+ !
+ ZEXN(:,:,: ) = (PP(II_M:II_M+1,IJ_M:IJ_M+1,:)/XP00)**(XRD/XCPD)
+ DO JK=IKB-1,1,-1
+ ZEXN(:,:,JK) = 1.5 * ZEXN(:,:,JK+1) - 0.5 * ZEXN(:,:,JK+2)
+ END DO
+ DO JK=IKE+1,IKU
+ ZEXN(:,:,JK) = 1.5 * ZEXN(:,:,JK-1) - 0.5 * ZEXN(:,:,JK-2)
+ END DO
+ !
+ ! IF ( TPFLYER%CTYPE == 'ISODEN' .OR. TPFLYER%CTYPE == 'CVBALL' .OR. TPFLYER%CTYPE == 'AIRCRA' ) THEN
+ ZTHV(:,:,:) = PTH(II_M:II_M+1,IJ_M:IJ_M+1,:)
+ IF (SIZE(PR,4)>0) &
+ ZTHV(:,:,:) = ZTHV(:,:,:) * ( 1. + XRV/XRD*PR(II_M:II_M+1,IJ_M:IJ_M+1,:,1) ) &
+ / ( 1. + WATER_SUM(PR(II_M:II_M+1,IJ_M:IJ_M+1,:,:)) )
+ !
+ ZTV (:,:,:) = ZTHV(:,:,:) * ZEXN(:,:,:)
+ ZRHO(:,:,:) = PP(II_M:II_M+1,IJ_M:IJ_M+1,:) / (XRD*ZTV(:,:,:))
+ DO JK=IKB-1,1,-1
+ ZRHO(:,:,JK) = 1.5 * ZRHO(:,:,JK+1) - 0.5 * ZRHO(:,:,JK+2)
+ END DO
+ DO JK=IKE+1,IKU
+ ZRHO(:,:,JK) = 1.5 * ZRHO(:,:,JK-1) - 0.5 * ZRHO(:,:,JK-2)
+ END DO
+ ZTHW_FLUX(:,:,:) = ZRHO(:,:,:)*XCPD *XTHW_FLUX(II_M:II_M+1,IJ_M:IJ_M+1,:)
+ ZRCW_FLUX(:,:,:) = ZRHO(:,:,:)*XLVTT*XRCW_FLUX(II_M:II_M+1,IJ_M:IJ_M+1,:)
+ ZSVW_FLUX(:,:,:,:) = XSVW_FLUX(II_M:II_M+1,IJ_M:IJ_M+1,:,:)
+ ! END IF
+
+ ! Vertical position
+ IF ( TPFLYER%LALTDEF ) THEN
+ ZFLYER_EXN = (TPFLYER%XP_CUR/XP00)**(XRD/XCPD)
+ IK00 = MAX ( COUNT (ZFLYER_EXN <= ZEXN(1,1,:)), 1)
+ IK01 = MAX ( COUNT (ZFLYER_EXN <= ZEXN(1,2,:)), 1)
+ IK10 = MAX ( COUNT (ZFLYER_EXN <= ZEXN(2,1,:)), 1)
+ IK11 = MAX ( COUNT (ZFLYER_EXN <= ZEXN(2,2,:)), 1)
+ ELSE
+ IK00 = MAX ( COUNT (TPFLYER%XZ_CUR >= ZZM(1,1,:)), 1)
+ IK01 = MAX ( COUNT (TPFLYER%XZ_CUR >= ZZM(1,2,:)), 1)
+ IK10 = MAX ( COUNT (TPFLYER%XZ_CUR >= ZZM(2,1,:)), 1)
+ IK11 = MAX ( COUNT (TPFLYER%XZ_CUR >= ZZM(2,2,:)), 1)
+ END IF
+
+ IF ( ANY( [ IK00, IK01, IK10, IK11 ] < IKB ) ) THEN
+ !Minimum altitude is on the ground at IKB (no crash if too low)
+ IK00 = MAX ( IK00, IKB )
+ IK01 = MAX ( IK01, IKB )
+ IK10 = MAX ( IK10, IKB )
+ IK11 = MAX ( IK11, IKB )
+
+ CMNHMSG(1) = 'flyer ' // TRIM( TPFLYER%CTITLE ) // ' is near the ground'
+ WRITE( CMNHMSG(2), "( 'at ', I2, '/', I2, '/', I4, ' ', F18.12, 's' )" ) &
+ TDTCUR%NDAY, TDTCUR%NMONTH, TDTCUR%NYEAR, TDTCUR%XTIME
+ CALL PRINT_MSG( NVERB_INFO, 'GEN', 'AIRCRAFT_BALLOON_EVOL', OLOCAL = .TRUE. )
+ END IF
+
+ IF (IK00 >= IKE .OR. IK01 >= IKE .OR. IK10 >= IKE .OR. IK11 >= IKE ) THEN
+ TPFLYER%LCRASH = .TRUE.
+ TPFLYER%NCRASH = NCRASH_OUT_HIGH
+ END IF
+ !
+ !* 6. INITIALIZATIONS FOR INTERPOLATIONS
+ ! ----------------------------------
+ !
+ !* 6.1 Interpolation coefficient for X
+ ! -------------------------------
+ !
+ ZXCOEF = (TPFLYER%XX_CUR - XXHATM(II_M)) / (XXHATM(II_M+1) - XXHATM(II_M))
+ ZXCOEF = MAX (0.,MIN(ZXCOEF,1.))
+ !
+ !
+ !* 6.2 Interpolation coefficient for y
+ ! -------------------------------
+ !
+ ZYCOEF = (TPFLYER%XY_CUR - XYHATM(IJ_M)) / (XYHATM(IJ_M+1) - XYHATM(IJ_M))
+ ZYCOEF = MAX (0.,MIN(ZYCOEF,1.))
+ !
+ !
+ !* 6.3 Interpolation coefficients for the 4 suroundings verticals
+ ! ----------------------------------------------------------
+ !
+ ! SELECT TYPE ( TPFLYER )
+ ! CLASS IS ( TBALLOONDATA)
+ ! IF ( TPFLYER%CTYPE == 'ISODEN' ) THEN
+ ! ZZCOEF00 = (TPFLYER%XRHO - ZRHO(1,1,IK00)) / ( ZRHO(1,1,IK00+1) - ZRHO(1,1,IK00) )
+ ! ZZCOEF01 = (TPFLYER%XRHO - ZRHO(1,2,IK01)) / ( ZRHO(1,2,IK01+1) - ZRHO(1,2,IK01) )
+ ! ZZCOEF10 = (TPFLYER%XRHO - ZRHO(2,1,IK10)) / ( ZRHO(2,1,IK10+1) - ZRHO(2,1,IK10) )
+ ! ZZCOEF11 = (TPFLYER%XRHO - ZRHO(2,2,IK11)) / ( ZRHO(2,2,IK11+1) - ZRHO(2,2,IK11) )
+ ! TPFLYER%XZ_CUR = FLYER_INTERP(ZZM)
+ ! ELSE IF ( TPFLYER%CTYPE == 'RADIOS' .OR. TPFLYER%CTYPE == 'CVBALL' ) THEN
+ ! ZZCOEF00 = (TPFLYER%XZ_CUR - ZZM(1,1,IK00)) / ( ZZM(1,1,IK00+1) - ZZM(1,1,IK00) )
+ ! ZZCOEF01 = (TPFLYER%XZ_CUR - ZZM(1,2,IK01)) / ( ZZM(1,2,IK01+1) - ZZM(1,2,IK01) )
+ ! ZZCOEF10 = (TPFLYER%XZ_CUR - ZZM(2,1,IK10)) / ( ZZM(2,1,IK10+1) - ZZM(2,1,IK10) )
+ ! ZZCOEF11 = (TPFLYER%XZ_CUR - ZZM(2,2,IK11)) / ( ZZM(2,2,IK11+1) - ZZM(2,2,IK11) )
+ ! END IF
+ !
+ ! CLASS IS ( TAIRCRAFTDATA)
+ IF ( TPFLYER%LALTDEF ) THEN
+ ZZCOEF00 = (ZFLYER_EXN - ZEXN(1,1,IK00)) / ( ZEXN(1,1,IK00+1) - ZEXN(1,1,IK00) )
+ ZZCOEF01 = (ZFLYER_EXN - ZEXN(1,2,IK01)) / ( ZEXN(1,2,IK01+1) - ZEXN(1,2,IK01) )
+ ZZCOEF10 = (ZFLYER_EXN - ZEXN(2,1,IK10)) / ( ZEXN(2,1,IK10+1) - ZEXN(2,1,IK10) )
+ ZZCOEF11 = (ZFLYER_EXN - ZEXN(2,2,IK11)) / ( ZEXN(2,2,IK11+1) - ZEXN(2,2,IK11) )
+ TPFLYER%XZ_CUR = FLYER_INTERP(ZZM)
+ ELSE
+ ZZCOEF00 = (TPFLYER%XZ_CUR - ZZM(1,1,IK00)) / ( ZZM(1,1,IK00+1) - ZZM(1,1,IK00) )
+ ZZCOEF01 = (TPFLYER%XZ_CUR - ZZM(1,2,IK01)) / ( ZZM(1,2,IK01+1) - ZZM(1,2,IK01) )
+ ZZCOEF10 = (TPFLYER%XZ_CUR - ZZM(2,1,IK10)) / ( ZZM(2,1,IK10+1) - ZZM(2,1,IK10) )
+ ZZCOEF11 = (TPFLYER%XZ_CUR - ZZM(2,2,IK11)) / ( ZZM(2,2,IK11+1) - ZZM(2,2,IK11) )
+ TPFLYER%XP_CUR = FLYER_INTERP(PP)
+ END IF
+ ! END SELECT
+ !
+ !----------------------------------------------------------------------------
+ !
+ !* 7. INITIALIZATIONS FOR INTERPOLATIONS OF U AND V
+ ! ---------------------------------------------
+ !
+ !* 7.1 Interpolation coefficient for X (for U)
+ ! -------------------------------
+ !
+ ZUCOEF = (TPFLYER%XX_CUR - XXHAT(II_U)) / (XXHAT(II_U+1) - XXHAT(II_U))
+ ZUCOEF = MAX(0.,MIN(ZUCOEF,1.))
+ !
+ !
+ !* 7.2 Interpolation coefficient for y (for V)
+ ! -------------------------------
+ !
+ ZVCOEF = (TPFLYER%XY_CUR - XYHAT(IJ_V)) / (XYHAT(IJ_V+1) - XYHAT(IJ_V))
+ ZVCOEF = MAX(0.,MIN(ZVCOEF,1.))
+ !
+ !
+ !* 7.3 Interpolation coefficients for the 4 suroundings verticals (for U)
+ ! ----------------------------------------------------------
+ !
+ IU00 = MAX( COUNT (TPFLYER%XZ_CUR >= ZZU(1,1,:)), 1)
+ IU01 = MAX( COUNT (TPFLYER%XZ_CUR >= ZZU(1,2,:)), 1)
+ IU10 = MAX( COUNT (TPFLYER%XZ_CUR >= ZZU(2,1,:)), 1)
+ IU11 = MAX( COUNT (TPFLYER%XZ_CUR >= ZZU(2,2,:)), 1)
+ ZUCOEF00 = (TPFLYER%XZ_CUR - ZZU(1,1,IU00)) / ( ZZU(1,1,IU00+1) - ZZU(1,1,IU00) )
+ ZUCOEF01 = (TPFLYER%XZ_CUR - ZZU(1,2,IU01)) / ( ZZU(1,2,IU01+1) - ZZU(1,2,IU01) )
+ ZUCOEF10 = (TPFLYER%XZ_CUR - ZZU(2,1,IU10)) / ( ZZU(2,1,IU10+1) - ZZU(2,1,IU10) )
+ ZUCOEF11 = (TPFLYER%XZ_CUR - ZZU(2,2,IU11)) / ( ZZU(2,2,IU11+1) - ZZU(2,2,IU11) )
+ !
+ !
+ !* 7.4 Interpolation coefficients for the 4 suroundings verticals (for V)
+ ! ----------------------------------------------------------
+ !
+ IV00 = MAX ( COUNT (TPFLYER%XZ_CUR >= ZZV(1,1,:)), 1)
+ IV01 = MAX ( COUNT (TPFLYER%XZ_CUR >= ZZV(1,2,:)), 1)
+ IV10 = MAX ( COUNT (TPFLYER%XZ_CUR >= ZZV(2,1,:)), 1)
+ IV11 = MAX ( COUNT (TPFLYER%XZ_CUR >= ZZV(2,2,:)), 1)
+ ZVCOEF00 = (TPFLYER%XZ_CUR - ZZV(1,1,IV00)) / ( ZZV(1,1,IV00+1) - ZZV(1,1,IV00) )
+ ZVCOEF01 = (TPFLYER%XZ_CUR - ZZV(1,2,IV01)) / ( ZZV(1,2,IV01+1) - ZZV(1,2,IV01) )
+ ZVCOEF10 = (TPFLYER%XZ_CUR - ZZV(2,1,IV10)) / ( ZZV(2,1,IV10+1) - ZZV(2,1,IV10) )
+ ZVCOEF11 = (TPFLYER%XZ_CUR - ZZV(2,2,IV11)) / ( ZZV(2,2,IV11+1) - ZZV(2,2,IV11) )
+ !
+ !
+ !* 8. DATA RECORDING
+ ! --------------
+ !
+ TPFLYER%NMODELHIST(ISTORE) = TPFLYER%NMODEL
+ TPFLYER%XX (ISTORE) = TPFLYER%XX_CUR
+ TPFLYER%XY (ISTORE) = TPFLYER%XY_CUR
+ TPFLYER%XZ (ISTORE) = TPFLYER%XZ_CUR
+ !
+ CALL SM_LATLON(PLATOR,PLONOR, &
+ TPFLYER%XX_CUR, TPFLYER%XY_CUR, &
+ TPFLYER%XLAT(ISTORE), TPFLYER%XLON(ISTORE) )
+ !
+ ZU_BAL = FLYER_INTERP_U(PU)
+ ZV_BAL = FLYER_INTERP_V(PV)
+ ZGAM = (XRPK * (TPFLYER%XLON(ISTORE) - XLON0) - XBETA)*(XPI/180.)
+ TPFLYER%XZON (ISTORE) = ZU_BAL * COS(ZGAM) + ZV_BAL * SIN(ZGAM)
+ TPFLYER%XMER (ISTORE) = - ZU_BAL * SIN(ZGAM) + ZV_BAL * COS(ZGAM)
+ !
+ TPFLYER%XW (ISTORE) = FLYER_INTERP(ZWM)
+ TPFLYER%XTH (ISTORE) = FLYER_INTERP(PTH)
+ !
+ ZFLYER_EXN = FLYER_INTERP(ZEXN)
+ TPFLYER%XP (ISTORE) = XP00 * ZFLYER_EXN**(XCPD/XRD)
+
+ ZR(:,:,:) = 0.
+ DO JLOOP=1,SIZE(PR,4)
+ TPFLYER%XR (ISTORE,JLOOP) = FLYER_INTERP(PR(:,:,:,JLOOP))
+ IF (JLOOP>=2) ZR(:,:,:) = ZR(:,:,:) + PR(:,:,:,JLOOP)
+ END DO
+ DO JLOOP=1,SIZE(PSV,4)
+ TPFLYER%XSV (ISTORE,JLOOP) = FLYER_INTERP(PSV(:,:,:,JLOOP))
+ END DO
+ TPFLYER%XRTZ (ISTORE,:) = FLYER_INTERPZ(ZR(:,:,:))
+ DO JLOOP=1,SIZE(PR,4)
+ TPFLYER%XRZ (ISTORE,:,JLOOP) = FLYER_INTERPZ(PR(:,:,:,JLOOP))
+ END DO
+ ! Fin Modifs ON
+ TPFLYER%XFFZ (ISTORE,:) = FLYER_INTERPZ(SQRT(PU**2+PV**2))
+ IF (CCLOUD=="LIMA") THEN
+ TPFLYER%XCIZ (ISTORE,:) = FLYER_INTERPZ(PSV(:,:,:,NSV_LIMA_NI))
+ TPFLYER%XCCZ (ISTORE,:) = FLYER_INTERPZ(PSV(:,:,:,NSV_LIMA_NC))
+ TPFLYER%XCRZ (ISTORE,:) = FLYER_INTERPZ(PSV(:,:,:,NSV_LIMA_NR))
+ ELSE IF ( CCLOUD=="ICE3" .OR. CCLOUD=="ICE4" ) THEN
+ TPFLYER%XCIZ (ISTORE,:) = FLYER_INTERPZ(PCIT(:,:,:))
+ ENDIF
+ ! initialization CRARE and CRARE_ATT + LWC and IWC
+ TPFLYER%XCRARE(ISTORE,:) = 0.
+ TPFLYER%XCRARE_ATT(ISTORE,:) = 0.
+ TPFLYER%XLWCZ (ISTORE,:) = 0.
+ TPFLYER%XIWCZ (ISTORE,:) = 0.
+ IF (CCLOUD=="LIMA" .OR. CCLOUD=="ICE3" ) THEN ! only for ICE3 and LIMA
+ TPFLYER%XLWCZ (ISTORE,:) = FLYER_INTERPZ((PR(:,:,:,2)+PR(:,:,:,3))*PRHODREF(:,:,:))
+ TPFLYER%XIWCZ (ISTORE,:) = FLYER_INTERPZ((PR(:,:,:,4)+PR(:,:,:,5)+PR(:,:,:,6))*PRHODREF(:,:,:))
+ ZTEMPZ(:)=FLYER_INTERPZ(PTH(II_M:II_M+1,IJ_M:IJ_M+1,:) * ZEXN(:,:,:))
+ ZRHODREFZ(:)=FLYER_INTERPZ(PRHODREF(:,:,:))
+ IF (CCLOUD=="LIMA") THEN
+ ZCCI(:)=FLYER_INTERPZ(PSV(:,:,:,NSV_LIMA_NI))
+ ZCCR(:)=FLYER_INTERPZ(PSV(:,:,:,NSV_LIMA_NR))
+ ZCCC(:)=FLYER_INTERPZ(PSV(:,:,:,NSV_LIMA_NC))
+ ELSE
+ ZCIT(:)=FLYER_INTERPZ(PCIT(:,:,:))
+ ENDIF
+ DO JLOOP=3,6
+ ZRZ(:,JLOOP)=FLYER_INTERPZ(PR(:,:,:,JLOOP))
+ END DO
+ if ( csurf == 'EXTE' ) then
+ DO JK=1,IKU
+ ZRZ(JK,2)=FLYER_INTERP_2D(PR(:,:,JK,2)*PSEA(:,:)) ! becomes cloud mixing ratio over sea
+ ZRZ(JK,7)=FLYER_INTERP_2D(PR(:,:,JK,2)*(1.-PSEA(:,:))) ! becomes cloud mixing ratio over land
+ END DO
+ else
+ !if csurf/='EXTE', psea is not allocated
+ DO JK=1,IKU
+ ZRZ(JK,2)=FLYER_INTERP_2D(PR(:,:,JK,2))
+ ZRZ(JK,7) = 0.
+ END DO
+ end if
+ ALLOCATE(ZAELOC(IKU))
+ !
+ ZAELOC(:)=0.
+ ! initialization of quadrature points and weights
+ ALLOCATE(ZX(JPTS_GAULAG),ZW(JPTS_GAULAG))
+ CALL GAULAG(JPTS_GAULAG,ZX,ZW) ! for integration over diameters
+ ! initialize minimum values
+ ALLOCATE(ZRTMIN(SIZE(PR,4)+1))
+ IF (CCLOUD == 'LIMA') THEN
+ ZRTMIN(2)=XRTMIN_L(2) ! cloud water over sea
+ ZRTMIN(3)=XRTMIN_L(3)
+ ZRTMIN(4)=XRTMIN_L(4)
+ ZRTMIN(5)=1E-10
+ ZRTMIN(6)=XRTMIN_L(6)
+ ZRTMIN(7)=XRTMIN_L(2) ! cloud water over land
+ ELSE
+ ZRTMIN(2)=XRTMIN_I(2) ! cloud water over sea
+ ZRTMIN(3)=XRTMIN_I(3)
+ ZRTMIN(4)=XRTMIN_I(4)
+ ZRTMIN(5)=1E-10
+ ZRTMIN(6)=XRTMIN_I(6)
+ ZRTMIN(7)=XRTMIN_I(2) ! cloud water over land
+ ENDIF
+ ! compute cloud radar reflectivity from vertical profiles of temperature and mixing ratios
+ DO JK=1,IKU
+ QMW=SQRT(QEPSW(ZTEMPZ(JK),XLIGHTSPEED/XLAM_CRAD))
+ QMI=SQRT(QEPSI(ZTEMPZ(JK),XLIGHTSPEED/XLAM_CRAD))
+ DO JLOOP=2,7
+ IF (CCLOUD == 'LIMA') THEN
+ GCALC=(ZRZ(JK,JLOOP)>ZRTMIN(JLOOP).AND.(JLOOP.NE.4.OR.ZCCI(JK)>0.).AND.&
+ (JLOOP.NE.3.OR.ZCCR(JK)>0.).AND.((JLOOP.NE.2.AND. JLOOP.NE.7).OR.ZCCC(JK)>0.))
+ ELSE
+ GCALC=(ZRZ(JK,JLOOP)>ZRTMIN(JLOOP).AND.(JLOOP.NE.4.OR.ZCIT(JK)>0.))
+ ENDIF
+ IF(GCALC) THEN
+ SELECT CASE(JLOOP)
+ CASE(2) ! cloud water over sea
+ IF (CCLOUD == 'LIMA') THEN
+ ZA=XAC_L
+ ZB=XBC_L
+ ZCC=ZCCC(JK)*ZRHODREFZ(JK)
+ ZCX=0.
+ ZALPHA=XALPHAC_L
+ ZNU=XNUC_L
+ ZLBEX=1.0/(ZCX-ZB)
+ ZLB=( ZA*ZCC*MOMG(ZALPHA,ZNU,ZB) )**(-ZLBEX)
+ ELSE
+ ZA=XAC_I
+ ZB=XBC_I
+ ZCC=XCONC_SEA
+ ZCX=0.
+ ZALPHA=XALPHAC2_I
+ ZNU=XNUC2_I
+ ZLBEX=1.0/(ZCX-ZB)
+ ZLB=( ZA*ZCC*MOMG(ZALPHA,ZNU,ZB) )**(-ZLBEX)
+ ENDIF
+ CASE(3) ! rain water
+ IF (CCLOUD == 'LIMA') THEN
+ ZA=XAR_L
+ ZB=XBR_L
+ ZCC=ZCCR(JK)*ZRHODREFZ(JK)
+ ZCX=0.
+ ZALPHA=XALPHAR_L
+ ZNU=XNUR_L
+ ZLBEX=1.0/(ZCX-ZB)
+ ZLB=( ZA*ZCC*MOMG(ZALPHA,ZNU,ZB) )**(-ZLBEX)
+ ELSE
+ ZA=XAR_I
+ ZB=XBR_I
+ ZCC=XCCR_I
+ ZCX=-1.
+ ZALPHA=XALPHAR_I
+ ZNU=XNUR_I
+ ZLB=XLBR_I
+ ZLBEX=XLBEXR_I
+ ENDIF
+ CASE(4) ! pristine ice
+ IF (CCLOUD == 'LIMA') THEN
+ ZA=XAI_L
+ ZB=XBI_L
+ ZCC=ZCCI(JK)*ZRHODREFZ(JK)
+ ZCX=0.
+ ZALPHA=XALPHAI_L
+ ZNU=XNUI_L
+ ZLBEX=1.0/(ZCX-ZB)
+ ZLB=( ZA*ZCC*MOMG(ZALPHA,ZNU,ZB) )**(-ZLBEX) ! because ZCC not included in XLBI
+ ZFW=0
+ ELSE
+ ZA=XAI_I
+ ZB=XBI_I
+ ZCC=ZCIT(JK)
+ ZCX=0.
+ ZALPHA=XALPHAI_I
+ ZNU=XNUI_I
+ ZLBEX=XLBEXI_I
+ ZLB=XLBI_I*ZCC**(-ZLBEX) ! because ZCC not included in XLBI
+ ZFW=0
+ ENDIF
+ CASE(5) ! snow
+ IF (CCLOUD == 'LIMA') THEN
+ ZA=XAS_L
+ ZB=XBS_L
+ ZCC=XCCS_L
+ ZCX=XCXS_L
+ ZALPHA=XALPHAS_L
+ ZNU=XNUS_L
+ ZLB=XLBS_L
+ ZLBEX=XLBEXS_L
+ ZFW=0
+ ELSE
+ ZA=XAS_I
+ ZB=XBS_I
+ ZCC=XCCS_I
+ ZCX=XCXS_I
+ ZALPHA=XALPHAS_I
+ ZNU=XNUS_I
+ ZLB=XLBS_I
+ ZLBEX=XLBEXS_I
+ ZFW=0
+ ENDIF
+ CASE(6) ! graupel
+ !If temperature between -10 and 10°C and Mr and Mg over min threshold: melting graupel
+ ! with liquid water fraction Fw=Mr/(Mr+Mg) else dry graupel (Fw=0)
+ IF( ZTEMPZ(JK) > XTT-10 .AND. ZTEMPZ(JK) < XTT+10 &
+ .AND. ZRZ(JK,3) > ZRTMIN(3) ) THEN
+ ZFW=ZRZ(JK,3)/(ZRZ(JK,3)+ZRZ(JK,JLOOP))
+ ELSE
+ ZFW=0
+ ENDIF
+ IF (CCLOUD == 'LIMA') THEN
+ ZA=XAG_L
+ ZB=XBG_L
+ ZCC=XCCG_L
+ ZCX=XCXG_L
+ ZALPHA=XALPHAG_L
+ ZNU=XNUG_L
+ ZLB=XLBG_L
+ ZLBEX=XLBEXG_L
+ ELSE
+ ZA=XAG_I
+ ZB=XBG_I
+ ZCC=XCCG_I
+ ZCX=XCXG_I
+ ZALPHA=XALPHAG_I
+ ZNU=XNUG_I
+ ZLB=XLBG_I
+ ZLBEX=XLBEXG_I
+ ENDIF
+ CASE(7) ! cloud water over land
+ IF (CCLOUD == 'LIMA') THEN
+ ZA=XAC_L
+ ZB=XBC_L
+ ZCC=ZCCC(JK)*ZRHODREFZ(JK)
+ ZCX=0.
+ ZALPHA=XALPHAC_L
+ ZNU=XNUC_L
+ ZLBEX=1.0/(ZCX-ZB)
+ ZLB=( ZA*ZCC*MOMG(ZALPHA,ZNU,ZB) )**(-ZLBEX)
+ ELSE
+ ZA=XAC_I
+ ZB=XBC_I
+ ZCC=XCONC_LAND
+ ZCX=0.
+ ZALPHA=XALPHAC_I
+ ZNU=XNUC_I
+ ZLBEX=1.0/(ZCX-ZB)
+ ZLB=( ZA*ZCC*MOMG(ZALPHA,ZNU,ZB) )**(-ZLBEX)
+ ENDIF
+ END SELECT
+ ZLBDA=ZLB*(ZRHODREFZ(JK)*ZRZ(JK,JLOOP))**ZLBEX
+ ZREFLOC=0.
+ ZAETMP=0.
+ DO JJ=1,JPTS_GAULAG ! Gauss-Laguerre quadrature
+ ZDELTA_EQUIV=ZX(JJ)**(1./ZALPHA)/ZLBDA
+ SELECT CASE(JLOOP)
+ CASE(2,3,7)
+ QM=QMW
+ CASE(4,5,6)
+ ! pristine ice, snow, dry graupel
+ ZRHOHYD=MIN(6.*ZA*ZDELTA_EQUIV**(ZB-3.)/XPI,.92*XRHOLW)
+ QM=sqrt(MG(QMI**2,CMPLX(1,0),ZRHOHYD/.92/XRHOLW))
+
+ ! water inclusions in ice in air
+ QEPSWI=MG(QMW**2,QM**2,ZFW)
+ ! ice in air inclusions in water
+ QEPSIW=MG(QM**2,QMW**2,1.-ZFW)
+
+ !MG weighted rule (Matrosov 2008)
+ IF(ZFW .LT. 0.37) THEN
+ ZFPW=0
+ ELSE IF(ZFW .GT. 0.63) THEN
+ ZFPW=1
+ ELSE
+ ZFPW=(ZFW-0.37)/(0.63-0.37)
+ ENDIF
+ QM=sqrt(QEPSWI*(1.-ZFPW)+QEPSIW*ZFPW)
+ END SELECT
+ CALL BHMIE(XPI/XLAM_CRAD*ZDELTA_EQUIV,QM,ZQEXT,ZQSCA,ZQBACK)
+ ZREFLOC=ZREFLOC+ZQBACK*ZX(JJ)**(ZNU-1)*ZDELTA_EQUIV**2*ZW(JJ)
+ ZAETMP =ZAETMP +ZQEXT *ZX(JJ)**(ZNU-1)*ZDELTA_EQUIV**2*ZW(JJ)
+ END DO
+ ZREFLOC=ZREFLOC*(XLAM_CRAD/XPI)**4*ZCC*ZLBDA**ZCX/(4.*GAMMA(ZNU)*.93)
+ ZAETMP=ZAETMP * XPI *ZCC*ZLBDA**ZCX/(4.*GAMMA(ZNU))
+ TPFLYER%XCRARE(ISTORE,JK)=TPFLYER%XCRARE(ISTORE,JK)+ZREFLOC
+ ZAELOC(JK)=ZAELOC(JK)+ZAETMP
+ END IF
+ END DO
+ END DO
+
+ ! apply attenuation
+ ALLOCATE(ZZMZ(IKU))
+ ZZMZ(:)=FLYER_INTERPZ(ZZM(:,:,:))
+ ! nadir
+ ZAETOT=1.
+ DO JK=COUNT(TPFLYER%XZ_CUR >= ZZMZ(:)),1,-1
+ IF(JK.EQ.COUNT(TPFLYER%XZ_CUR >= ZZMZ(:))) THEN
+ IF(TPFLYER%XZ_CUR<=ZZMZ(JK)+.5*(ZZMZ(JK+1)-ZZMZ(JK))) THEN
+ ! only attenuation from ZAELOC(JK)
+ ZAETOT=ZAETOT*EXP(-2.*(ZAELOC(JK)*(TPFLYER%XZ_CUR-ZZMZ(JK))))
+ ELSE
+ ! attenuation from ZAELOC(JK) and ZAELOC(JK+1)
+ ZAETOT=ZAETOT*EXP(-2.*(ZAELOC(JK+1)*(TPFLYER%XZ_CUR-.5*(ZZMZ(JK+1)+ZZMZ(JK))) &
+ +ZAELOC(JK)*.5*(ZZMZ(JK+1)-ZZMZ(JK))))
+ END IF
+ ELSE
+ ! attenuation from ZAELOC(JK) and ZAELOC(JK+1)
+ ZAETOT=ZAETOT*EXP(-(ZAELOC(JK+1)+ZAELOC(JK))*(ZZMZ(JK+1)-ZZMZ(JK)))
+ END IF
+ TPFLYER%XCRARE_ATT(ISTORE,JK)=TPFLYER%XCRARE(ISTORE,JK)*ZAETOT
+ END DO
+ ! zenith
+ ZAETOT=1.
+ DO JK = MAX(COUNT(TPFLYER%XZ_CUR >= ZZMZ(:)),1)+1,IKU
+ IF ( JK .EQ. (MAX(COUNT(TPFLYER%XZ_CUR >= ZZMZ(:)),1)+1) ) THEN
+ IF(TPFLYER%XZ_CUR>=ZZMZ(JK)-.5*(ZZMZ(JK)-ZZMZ(JK-1))) THEN
+ ! only attenuation from ZAELOC(JK)
+ ZAETOT=ZAETOT*EXP(-2.*(ZAELOC(JK)*(ZZMZ(JK)-TPFLYER%XZ_CUR)))
+ ELSE
+ ! attenuation from ZAELOC(JK) and ZAELOC(JK-1)
+ ZAETOT=ZAETOT*EXP(-2.*(ZAELOC(JK-1)*(.5*(ZZMZ(JK)+ZZMZ(JK-1))-TPFLYER%XZ_CUR) &
+ +ZAELOC(JK)*.5*(ZZMZ(JK)-ZZMZ(JK-1))))
+ END IF
+ ELSE
+ ! attenuation from ZAELOC(JK) and ZAELOC(JK-1)
+ ZAETOT=ZAETOT*EXP(-(ZAELOC(JK-1)+ZAELOC(JK))*(ZZMZ(JK)-ZZMZ(JK-1)))
+ END IF
+ TPFLYER%XCRARE_ATT(ISTORE,JK)=TPFLYER%XCRARE(ISTORE,JK)*ZAETOT
+ END DO
+ TPFLYER%XZZ (ISTORE,:) = ZZMZ(:)
+ DEALLOCATE(ZZMZ,ZAELOC)
+ ! m^3 → mm^6/m^3 → dBZ
+ WHERE(TPFLYER%XCRARE(ISTORE,:)>0)
+ TPFLYER%XCRARE(ISTORE,:)=10.*LOG10(1.E18*TPFLYER%XCRARE(ISTORE,:))
+ ELSEWHERE
+ TPFLYER%XCRARE(ISTORE,:)=XUNDEF
+ END WHERE
+ WHERE(TPFLYER%XCRARE_ATT(ISTORE,:)>0)
+ TPFLYER%XCRARE_ATT(ISTORE,:)=10.*LOG10(1.E18*TPFLYER%XCRARE_ATT(ISTORE,:))
+ ELSEWHERE
+ TPFLYER%XCRARE_ATT(ISTORE,:)=XUNDEF
+ END WHERE
+ DEALLOCATE(ZX,ZW,ZRTMIN)
+ END IF ! end LOOP ICE3
+ ! vertical wind
+ TPFLYER%XWZ (ISTORE,:) = FLYER_INTERPZ(ZWM(:,:,:))
+ IF (SIZE(PTKE)>0) TPFLYER%XTKE (ISTORE) = FLYER_INTERP(PTKE)
+ IF (SIZE(PTS) >0) TPFLYER%XTSRAD(ISTORE) = FLYER_INTERP_2D(PTS)
+ IF (LDIAG_IN_RUN) TPFLYER%XTKE_DISS(ISTORE) = FLYER_INTERP(XCURRENT_TKE_DISS)
+ TPFLYER%XZS(ISTORE) = FLYER_INTERP_2D(PZ(:,:,1+JPVEXT))
+ TPFLYER%XTHW_FLUX(ISTORE) = FLYER_INTERP(ZTHW_FLUX)
+ TPFLYER%XRCW_FLUX(ISTORE) = FLYER_INTERP(ZRCW_FLUX)
+ DO JLOOP=1,SIZE(PSV,4)
+ TPFLYER%XSVW_FLUX(ISTORE,JLOOP) = FLYER_INTERP(ZSVW_FLUX(:,:,:,JLOOP))
+ END DO
+
+ ELSE ISOWNER
+
+ ZTHIS_PROC = 0.
+
+ END IF ISOWNER
+
+!----------------------------------------------------------------------------
+ !
+ !* 11. EXCHANGE OF INFORMATION BETWEEN PROCESSES
+ ! -----------------------------------------
+ !
+ !* 11.1 current position
+ ! ----------------
+ !
!
- TPFLYER%XX_CUR = (1.-ZSEG_FRAC) * TPFLYER%XSEGX(IL ) &
- + ZSEG_FRAC * TPFLYER%XSEGX(IL+1)
- TPFLYER%XY_CUR = (1.-ZSEG_FRAC) * TPFLYER%XSEGY(IL ) &
- + ZSEG_FRAC * TPFLYER%XSEGY(IL+1)
+ !* 11.2 data stored
+ ! -----------
+ !
+ ! IF ( GSTORE ) THEN
+ CALL DISTRIBUTE_FLYER_N(TPFLYER%NMODELHIST(ISTORE))
+ CALL DISTRIBUTE_FLYER(TPFLYER%XX (ISTORE))
+ CALL DISTRIBUTE_FLYER(TPFLYER%XY (ISTORE))
+ CALL DISTRIBUTE_FLYER(TPFLYER%XZ (ISTORE))
+ CALL DISTRIBUTE_FLYER(TPFLYER%XLON(ISTORE))
+ CALL DISTRIBUTE_FLYER(TPFLYER%XLAT(ISTORE))
+ CALL DISTRIBUTE_FLYER(TPFLYER%XZON(ISTORE))
+ CALL DISTRIBUTE_FLYER(TPFLYER%XMER(ISTORE))
+ CALL DISTRIBUTE_FLYER(TPFLYER%XW (ISTORE))
+ CALL DISTRIBUTE_FLYER(TPFLYER%XP (ISTORE))
+ CALL DISTRIBUTE_FLYER(TPFLYER%XTH (ISTORE))
+ DO JLOOP=1,SIZE(PR,4)
+ CALL DISTRIBUTE_FLYER(TPFLYER%XR (ISTORE,JLOOP))
+ END DO
+ DO JLOOP=1,SIZE(PSV,4)
+ CALL DISTRIBUTE_FLYER(TPFLYER%XSV (ISTORE,JLOOP))
+ END DO
+ DO JLOOP=1,IKU
+ CALL DISTRIBUTE_FLYER(TPFLYER%XRTZ (ISTORE,JLOOP))
+ DO JLOOP2=1,SIZE(PR,4)
+ CALL DISTRIBUTE_FLYER(TPFLYER%XRZ (ISTORE,JLOOP,JLOOP2))
+ ENDDO
+ CALL DISTRIBUTE_FLYER(TPFLYER%XFFZ (ISTORE,JLOOP))
+ CALL DISTRIBUTE_FLYER(TPFLYER%XCIZ (ISTORE,JLOOP))
+ IF (CCLOUD== 'LIMA' ) THEN
+ CALL DISTRIBUTE_FLYER(TPFLYER%XCRZ (ISTORE,JLOOP))
+ CALL DISTRIBUTE_FLYER(TPFLYER%XCCZ (ISTORE,JLOOP))
+ ENDIF
+ CALL DISTRIBUTE_FLYER(TPFLYER%XIWCZ (ISTORE,JLOOP))
+ CALL DISTRIBUTE_FLYER(TPFLYER%XLWCZ (ISTORE,JLOOP))
+ CALL DISTRIBUTE_FLYER(TPFLYER%XCRARE (ISTORE,JLOOP))
+ CALL DISTRIBUTE_FLYER(TPFLYER%XCRARE_ATT (ISTORE,JLOOP))
+ CALL DISTRIBUTE_FLYER(TPFLYER%XWZ (ISTORE,JLOOP))
+ CALL DISTRIBUTE_FLYER(TPFLYER%XZZ (ISTORE,JLOOP))
+ END DO
+ IF (SIZE(PTKE)>0) CALL DISTRIBUTE_FLYER(TPFLYER%XTKE (ISTORE))
+ IF (SIZE(PTS) >0) CALL DISTRIBUTE_FLYER(TPFLYER%XTSRAD(ISTORE))
+ IF (LDIAG_IN_RUN) CALL DISTRIBUTE_FLYER(TPFLYER%XTKE_DISS(ISTORE))
+ CALL DISTRIBUTE_FLYER(TPFLYER%XZS (ISTORE))
+ CALL DISTRIBUTE_FLYER(TPFLYER%XTHW_FLUX(ISTORE))
+ CALL DISTRIBUTE_FLYER(TPFLYER%XRCW_FLUX(ISTORE))
+ DO JLOOP=1,SIZE(PSV,4)
+ CALL DISTRIBUTE_FLYER(TPFLYER%XSVW_FLUX(ISTORE,JLOOP))
+ END DO
END IF
- END SELECT
- END IF
-END IF
-!
-!* 3.4 instant of storage
-! ------------------
-!
-CALL STATPROF_INSTANT( TPFLYER%TFLYER_TIME, IN )
-IF ( IN > 0 ) GSTORE = .TRUE. ! else no profiler storage at this time step
-!
-IF ( TPFLYER%LFLY ) THEN
-!
-!----------------------------------------------------------------------------
-!
-!* 4. FLYER POSITION
-! --------------
-!
-!* 4.1 X position
-! ----------
-!
- IU=COUNT( PXHAT (:)<=TPFLYER%XX_CUR )
- II=COUNT( XXHATM(:)<=TPFLYER%XX_CUR )
-!
- IF ( IU < IIB .AND. LWEST_ll() ) THEN
- IF ( TPFLYER%CMODEL == 'FIX' .OR. TPFLYER%NMODEL == 1 ) THEN
- TPFLYER%LCRASH = .TRUE.
- ELSE
- II = IIB
- IU = IIB
- END IF
- END IF
- IF ( IU > IIE .AND. LEAST_ll() ) THEN
- IF ( TPFLYER%CMODEL == 'FIX' .OR. TPFLYER%NMODEL == 1 ) THEN
- TPFLYER%LCRASH = .TRUE.
- ELSE
- II = IIE
- IU = IIE
- END IF
- END IF
-!
-!
-!* 4.2 Y position
-! ----------
-!
- IV=COUNT( PYHAT (:)<=TPFLYER%XY_CUR )
- IJ=COUNT( XYHATM(:)<=TPFLYER%XY_CUR )
-!
- IF ( IV < IJB .AND. LSOUTH_ll() ) THEN
- IF ( TPFLYER%CMODEL == 'FIX' .OR. TPFLYER%NMODEL == 1 ) THEN
- TPFLYER%LCRASH = .TRUE.
- ELSE
- IJ = IJB
- IV = IJB
- END IF
- END IF
- IF (IV > IJE .AND. LNORTH_ll() ) THEN
- IF ( TPFLYER%CMODEL == 'FIX' .OR. TPFLYER%NMODEL == 1 ) THEN
- TPFLYER%LCRASH = .TRUE.
- ELSE
- IJ = IJE
- IV = IJE
+ END IF
END IF
- END IF
-!
-!
-!* 4.3 Position of balloon according to processors
-! -------------------------------------------
-!
- IF (IU>=IIB .AND. IU<=IIE .AND. IV>=IJB .AND. IV<=IJE) ZTHIS_PROC=1.
-!
-!
-!* 4.4 Computations only on correct processor
-! --------------------------------------
-!
-!----------------------------------------------------------------------------
- IF ( ZTHIS_PROC > 0. .AND. .NOT. TPFLYER%LCRASH ) THEN
-!----------------------------------------------------------------------------
-!
-!* 4.5 Interpolations of model variables to mass points
-! ------------------------------------------------
-!
- ZZM(:,:,1:IKU-1)=0.5 *PZ(II :II+1,IJ :IJ+1,1:IKU-1)+0.5 *PZ(II :II+1,IJ :IJ+1,2:IKU )
- ZZM(:,:, IKU )=1.5 *PZ(II :II+1,IJ :IJ+1, IKU-1)-0.5 *PZ(II :II+1,IJ :IJ+1, IKU-2)
-!
- IDU = IU - II
- ZZU(:,:,1:IKU-1)=0.25*PZ(IDU+II-1:IDU+II, IJ :IJ+1,1:IKU-1)+0.25*PZ(IDU+II-1:IDU+II ,IJ :IJ+1,2:IKU ) &
- +0.25*PZ(IDU+II :IDU+II+1,IJ :IJ+1,1:IKU-1)+0.25*PZ(IDU+II :IDU+II+1,IJ :IJ+1,2:IKU )
- ZZU(:,:, IKU )=0.75*PZ(IDU+II-1:IDU+II ,IJ :IJ+1, IKU-1)-0.25*PZ(IDU+II-1:IDU+II ,IJ :IJ+1, IKU-2) &
- +0.75*PZ(IDU+II :IDU+II+1,IJ :IJ+1, IKU-1)-0.25*PZ(IDU+II :IDU+II+1,IJ :IJ+1, IKU-2)
- IDV = IV - IJ
- ZZV(:,:,1:IKU-1)=0.25*PZ(II :II+1,IDV+IJ-1:IDV+IJ ,1:IKU-1)+0.25*PZ(II :II+1,IDV+IJ-1:IDV+IJ ,2:IKU ) &
- +0.25*PZ(II :II+1,IDV+IJ :IDV+IJ+1,1:IKU-1)+0.25*PZ(II :II+1,IDV+IJ :IDV+IJ+1,2:IKU )
- ZZV(:,:, IKU )=0.75*PZ(II :II+1,IDV+IJ-1:IDV+IJ , IKU-1)-0.25*PZ(II :II+1,IDV+IJ-1:IDV+IJ , IKU-2) &
- +0.75*PZ(II :II+1,IDV+IJ :IDV+IJ+1, IKU-1)-0.25*PZ(II :II+1,IDV+IJ :IDV+IJ+1, IKU-2)
-!
-!
- ZWM(:,:,1:IKU-1)=0.5*PW(II:II+1,IJ:IJ+1,1:IKU-1)+0.5*PW(II:II+1,IJ:IJ+1,2:IKU )
- ZWM(:,:, IKU )=1.5*PW(II:II+1,IJ:IJ+1, IKU-1)-0.5*PW(II:II+1,IJ:IJ+1, IKU-2)
-!
-!----------------------------------------------------------------------------
-!
-!* 5. BALLOON/AIRCRAFT VERTICAL POSITION
-! ----------------------------------
-!
-!
-!* 5.1 Density
-! -------
-!
- ZEXN(:,:,: ) = (PP(II:II+1,IJ:IJ+1,:)/XP00)**(XRD/XCPD)
- DO JK=IKB-1,1,-1
- ZEXN(:,:,JK) = 1.5 * ZEXN(:,:,JK+1) - 0.5 * ZEXN(:,:,JK+2)
- END DO
- DO JK=IKE+1,IKU
- ZEXN(:,:,JK) = 1.5 * ZEXN(:,:,JK-1) - 0.5 * ZEXN(:,:,JK-2)
- END DO
- !
- IF ( TPFLYER%CTYPE == 'ISODEN' .OR. TPFLYER%CTYPE == 'CVBALL' .OR. TPFLYER%CTYPE == 'AIRCRA' ) THEN
- ZTHV(:,:,:) = PTH(II:II+1,IJ:IJ+1,:)
- IF (SIZE(PR,4)>0) &
- ZTHV(:,:,:) = ZTHV(:,:,:) * ( 1. + XRV/XRD*PR(II:II+1,IJ:IJ+1,:,1) ) &
- / ( 1. + WATER_SUM(PR(II:II+1,IJ:IJ+1,:,:)) )
- !
- ZTV (:,:,:) = ZTHV(:,:,:) * ZEXN(:,:,:)
- ZRHO(:,:,:) = PP(II:II+1,IJ:IJ+1,:) / (XRD*ZTV(:,:,:))
- DO JK=IKB-1,1,-1
- ZRHO(:,:,JK) = 1.5 * ZRHO(:,:,JK+1) - 0.5 * ZRHO(:,:,JK+2)
- END DO
- DO JK=IKE+1,IKU
- ZRHO(:,:,JK) = 1.5 * ZRHO(:,:,JK-1) - 0.5 * ZRHO(:,:,JK-2)
- END DO
- ZTHW_FLUX(:,:,:) = ZRHO(:,:,:)*XCPD *XTHW_FLUX(II:II+1,IJ:IJ+1,:)
- ZRCW_FLUX(:,:,:) = ZRHO(:,:,:)*XLVTT*XRCW_FLUX(II:II+1,IJ:IJ+1,:)
- ZSVW_FLUX(:,:,:,:) = XSVW_FLUX(II:II+1,IJ:IJ+1,:,:)
+ CLASS IS ( TBALLOONDATA)
+ GLAUNCH = .FALSE. !Set to true only at the launch instant (set to false in flight after launch)
+
+ ! Initialize model number (and rank)
+ ! This is not done in initialisation phase because some data is not yet available at this early stage
+ ! (XXHAT_ll of all models are needed by FIND_PROCESS_AND_MODEL_FROM_XY_POS)
+ IF ( .NOT. TPFLYER%LFLY .AND. .NOT. TPFLYER%LCRASH .AND. TPFLYER%NRANK_CUR < 0 ) THEN
+ ! Get rank of the process where the balloon is and the model number
+ IF ( TPFLYER%CMODEL == 'FIX' ) THEN
+ IMODEL = TPFLYER%NMODEL
+ ELSE
+ IMODEL = 0
+ END IF
+ CALL FIND_PROCESS_AND_MODEL_FROM_XY_POS( TPFLYER%XXLAUNCH, TPFLYER%XYLAUNCH, IRANK, IMODEL )
+
+ IF ( IRANK < 1 ) THEN
+ TPFLYER%NMODEL = 1 !Set to 1 because it is always a valid model (to prevent crash if NDAD(TPFLYER%NMODEL) )
+ TPFLYER%LCRASH = .TRUE.
+ TPFLYER%NCRASH = NCRASH_OUT_HORIZ
+ TPFLYER%LFLY = .FALSE.
+ CALL PRINT_MSG( NVERB_WARNING, 'GEN', 'AIRCRAFT_BALLOON_EVOL', 'balloon ' // TRIM( TPFLYER%CTITLE ) &
+ // ': launch coordinates are outside of horizontal physical domain' )
+ ELSE
+ TPFLYER%NMODEL = IMODEL
+ TPFLYER%LCRASH = .FALSE.
+ TPFLYER%NCRASH = NCRASH_NO
+ TPFLYER%NRANK_CUR = IRANK
+ END IF
END IF
-!
-!* 5.2 Initial vertical positions
-! --------------------------
-!
- IF (GLAUNCH) THEN
- SELECT TYPE ( TPFLYER )
- CLASS IS ( TBALLOONDATA)
- SELECT CASE ( TPFLYER%CTYPE )
-!
-!* 5.2.1 Iso-density balloon
-!
- CASE ( 'ISODEN' )
- ZXCOEF = (TPFLYER%XX_CUR - XXHATM(II)) / (XXHATM(II+1) - XXHATM(II))
- ZXCOEF = MAX (0.,MIN(ZXCOEF,1.))
- ZYCOEF = (TPFLYER%XY_CUR - XYHATM(IJ)) / (XYHATM(IJ+1) - XYHATM(IJ))
- ZYCOEF = MAX (0.,MIN(ZYCOEF,1.))
- IF ( TPFLYER%XALTLAUNCH /= XNEGUNDEF ) THEN
- IK00 = MAX ( COUNT (TPFLYER%XALTLAUNCH >= ZZM(1,1,:)), 1)
- IK01 = MAX ( COUNT (TPFLYER%XALTLAUNCH >= ZZM(1,2,:)), 1)
- IK10 = MAX ( COUNT (TPFLYER%XALTLAUNCH >= ZZM(2,1,:)), 1)
- IK11 = MAX ( COUNT (TPFLYER%XALTLAUNCH >= ZZM(2,2,:)), 1)
- ZZCOEF00 = (TPFLYER%XALTLAUNCH - ZZM(1,1,IK00)) / ( ZZM(1,1,IK00+1) - ZZM(1,1,IK00))
- ZZCOEF01 = (TPFLYER%XALTLAUNCH - ZZM(1,2,IK01)) / ( ZZM(1,2,IK01+1) - ZZM(1,2,IK01))
+ ! Launch?
+ LAUNCH: IF ( .NOT. TPFLYER%LFLY .AND. .NOT. TPFLYER%LCRASH .AND. TPFLYER%NMODEL == IMI ) THEN
+ LAUNCHTIME: IF ( ( TDTCUR - TPFLYER%TLAUNCH ) >= -1.e-10 ) THEN
+ TPFLYER%LFLY = .TRUE.
+ GLAUNCH = .TRUE.
+
+ TPFLYER%XX_CUR = TPFLYER%XXLAUNCH
+ TPFLYER%XY_CUR = TPFLYER%XYLAUNCH
+ TPFLYER%TPOS_CUR = TDTCUR
+ ! Get rank of the process where the balloon is and the model number
+ IF ( TPFLYER%CMODEL == 'FIX' ) THEN
+ IMODEL = TPFLYER%NMODEL
+ ELSE
+ IMODEL = 0
+ END IF
+ CALL FIND_PROCESS_AND_MODEL_FROM_XY_POS( TPFLYER%XX_CUR, TPFLYER%XY_CUR, IRANK, IMODEL )
+ IF ( IRANK < 1 ) THEN
+ TPFLYER%NMODEL = 1 !Set to 1 because it is always a valid model (to prevent crash if NDAD(TPFLYER%NMODEL) )
+ TPFLYER%LCRASH = .TRUE.
+ TPFLYER%NCRASH = NCRASH_OUT_HORIZ
+ TPFLYER%LFLY = .FALSE.
+ WRITE( CMNHMSG(1), "( 'Balloon ', A, ' crashed the ', I2, '/', I2, '/', I4, ' at ', F18.12, &
+ 's (out of the horizontal boundaries)' )" ) &
+ TRIM( TPFLYER%CTITLE ), TDTCUR%NDAY, TDTCUR%NMONTH, TDTCUR%NYEAR, TDTCUR%XTIME
+ CALL PRINT_MSG( NVERB_WARNING, 'GEN', 'AIRCRAFT_BALLOON_EVOL' )
+ ELSE
+ TPFLYER%NMODEL = IMODEL
+ TPFLYER%LCRASH = .FALSE.
+ TPFLYER%NCRASH = NCRASH_NO
+ TPFLYER%NRANK_CUR = IRANK
+ END IF
+ END IF LAUNCHTIME
+ END IF LAUNCH
+
+ MODEL1BAL: IF ( TPFLYER%NMODEL == IMI .AND. &
+ !Do we have to store balloon data?
+ CALL STATPROF_INSTANT( TPFLYER%TFLYER_TIME, ISTORE )
+ IF ( ISTORE < 1 ) THEN
+ !No profiler storage at this time step
+ TPFLYER%LSTORE = .FALSE.
+ ELSE
+ TPFLYER%LSTORE = .TRUE.
+ END IF
+ END IF MODEL1BAL
+
+ ! In flight
+ INFLIGHTONMODEL: IF ( TPFLYER%LFLY .AND. .NOT. TPFLYER%LCRASH .AND. TPFLYER%NMODEL == IMI &
+ .AND. ABS( TPFLYER%TPOS_CUR - TDTCUR ) < 1.e-8 ) THEN
+
+ ISOWNERBAL: IF ( TPFLYER%NRANK_CUR == ISP ) THEN
+ ZTHIS_PROC = 1.
+ !
+ !* 2. PRELIMINARIES-2
+ ! -------------
+ !
+ !* 2.1 Indices
+ ! -------
+ !
+ CALL GET_INDICE_ll (IIB,IJB,IIE,IJE)
+ IKB = 1 + JPVEXT
+ IKE = SIZE(PZ,3) - JPVEXT
+ !
+ !
+ !* 2.2 Interpolations of model variables to mass points
+ ! ------------------------------------------------
+ !
+ IIU=SIZE(XXHAT)
+ IJU=SIZE(XYHAT)
+ ! X position
+ II_U = COUNT( XXHAT (:) <= TPFLYER%XX_CUR )
+ II_M = COUNT( XXHATM(:) <= TPFLYER%XX_CUR )
+
+ ! Y position
+ IJ_V=COUNT( XYHAT (:)<=TPFLYER%XY_CUR )
+ IJ_M=COUNT( XYHATM(:)<=TPFLYER%XY_CUR )
+ !
+ !* 4.5 Interpolations of model variables to mass points
+ ! ------------------------------------------------
+ !
+ ZZM(:,:,1:IKU-1)=0.5 *PZ(II_M :II_M+1,IJ_M :IJ_M+1,1:IKU-1)+0.5 *PZ(II_M :II_M+1,IJ_M :IJ_M+1,2:IKU )
+ ZZM(:,:, IKU )=1.5 *PZ(II_M :II_M+1,IJ_M :IJ_M+1, IKU-1)-0.5 *PZ(II_M :II_M+1,IJ_M :IJ_M+1, IKU-2)
+
+ IDU = II_U - II_M
+ ZZU(:,:,1:IKU-1)=0.25*PZ(IDU+II_M-1:IDU+II_M, IJ_M :IJ_M+1,1:IKU-1)+0.25*PZ(IDU+II_M-1:IDU+II_M ,IJ_M :IJ_M+1,2:IKU ) &
+ +0.25*PZ(IDU+II_M :IDU+II_M+1,IJ_M :IJ_M+1,1:IKU-1)+0.25*PZ(IDU+II_M :IDU+II_M+1,IJ_M :IJ_M+1,2:IKU )
+ ZZU(:,:, IKU )=0.75*PZ(IDU+II_M-1:IDU+II_M ,IJ_M :IJ_M+1, IKU-1)-0.25*PZ(IDU+II_M-1:IDU+II_M ,IJ_M :IJ_M+1, IKU-2) &
+ +0.75*PZ(IDU+II_M :IDU+II_M+1,IJ_M :IJ_M+1, IKU-1)-0.25*PZ(IDU+II_M :IDU+II_M+1,IJ_M :IJ_M+1, IKU-2)
+
+ IDV = IJ_V - IJ_M
+ ZZV(:,:,1:IKU-1)=0.25*PZ(II_M :II_M+1,IDV+IJ_M-1:IDV+IJ_M ,1:IKU-1)+0.25*PZ(II_M :II_M+1,IDV+IJ_M-1:IDV+IJ_M ,2:IKU ) &
+ +0.25*PZ(II_M :II_M+1,IDV+IJ_M :IDV+IJ_M+1,1:IKU-1)+0.25*PZ(II_M :II_M+1,IDV+IJ_M :IDV+IJ_M+1,2:IKU )
+ ZZV(:,:, IKU )=0.75*PZ(II_M :II_M+1,IDV+IJ_M-1:IDV+IJ_M , IKU-1)-0.25*PZ(II_M :II_M+1,IDV+IJ_M-1:IDV+IJ_M , IKU-2) &
+ +0.75*PZ(II_M :II_M+1,IDV+IJ_M :IDV+IJ_M+1, IKU-1)-0.25*PZ(II_M :II_M+1,IDV+IJ_M :IDV+IJ_M+1, IKU-2)
+
+ ZWM(:,:,1:IKU-1)=0.5*PW(II_M:II_M+1,IJ_M:IJ_M+1,1:IKU-1)+0.5*PW(II_M:II_M+1,IJ_M:IJ_M+1,2:IKU )
+ ZWM(:,:, IKU )=1.5*PW(II_M:II_M+1,IJ_M:IJ_M+1, IKU-1)-0.5*PW(II_M:II_M+1,IJ_M:IJ_M+1, IKU-2)
+ !
+ !----------------------------------------------------------------------------
+ !
+ !* 5. BALLOON/AIRCRAFT VERTICAL POSITION
+ ! ----------------------------------
+ !
+ !
+ !* 5.1 Density
+ ! -------
+ !
+ ZEXN(:,:,: ) = (PP(II_M:II_M+1,IJ_M:IJ_M+1,:)/XP00)**(XRD/XCPD)
+ DO JK=IKB-1,1,-1
+ ZEXN(:,:,JK) = 1.5 * ZEXN(:,:,JK+1) - 0.5 * ZEXN(:,:,JK+2)
+ END DO
+ DO JK=IKE+1,IKU
+ ZEXN(:,:,JK) = 1.5 * ZEXN(:,:,JK-1) - 0.5 * ZEXN(:,:,JK-2)
+ END DO
+ !
+ IF ( TPFLYER%CTYPE == 'ISODEN' .OR. TPFLYER%CTYPE == 'CVBALL' .OR. TPFLYER%CTYPE == 'AIRCRA' ) THEN
+ ZTHV(:,:,:) = PTH(II_M:II_M+1,IJ_M:IJ_M+1,:)
+ IF (SIZE(PR,4)>0) &
+ ZTHV(:,:,:) = ZTHV(:,:,:) * ( 1. + XRV/XRD*PR(II_M:II_M+1,IJ_M:IJ_M+1,:,1) ) &
+ / ( 1. + WATER_SUM(PR(II_M:II_M+1,IJ_M:IJ_M+1,:,:)) )
+ !
+ ZTV (:,:,:) = ZTHV(:,:,:) * ZEXN(:,:,:)
+ ZRHO(:,:,:) = PP(II_M:II_M+1,IJ_M:IJ_M+1,:) / (XRD*ZTV(:,:,:))
+ DO JK=IKB-1,1,-1
+ ZRHO(:,:,JK) = 1.5 * ZRHO(:,:,JK+1) - 0.5 * ZRHO(:,:,JK+2)
+ END DO
+ DO JK=IKE+1,IKU
+ ZRHO(:,:,JK) = 1.5 * ZRHO(:,:,JK-1) - 0.5 * ZRHO(:,:,JK-2)
+ END DO
+ ZTHW_FLUX(:,:,:) = ZRHO(:,:,:)*XCPD *XTHW_FLUX(II_M:II_M+1,IJ_M:IJ_M+1,:)
+ ZRCW_FLUX(:,:,:) = ZRHO(:,:,:)*XLVTT*XRCW_FLUX(II_M:II_M+1,IJ_M:IJ_M+1,:)
+ ZSVW_FLUX(:,:,:,:) = XSVW_FLUX(II_M:II_M+1,IJ_M:IJ_M+1,:,:)
+ END IF
+ SELECT CASE ( TPFLYER%CTYPE )
+ !
+ !* 5.2.1 Iso-density balloon
+ !
+ CASE ( 'ISODEN' )
+ ZXCOEF = (TPFLYER%XX_CUR - XXHATM(II_M)) / (XXHATM(II_M+1) - XXHATM(II_M))
+ ZXCOEF = MAX (0.,MIN(ZXCOEF,1.))
+ ZYCOEF = (TPFLYER%XY_CUR - XYHATM(IJ_M)) / (XYHATM(IJ_M+1) - XYHATM(IJ_M))
+ ZYCOEF = MAX (0.,MIN(ZYCOEF,1.))
+ IF ( TPFLYER%XALTLAUNCH /= XNEGUNDEF ) THEN
+ IK00 = MAX ( COUNT (TPFLYER%XALTLAUNCH >= ZZM(1,1,:)), 1)
+ IK01 = MAX ( COUNT (TPFLYER%XALTLAUNCH >= ZZM(1,2,:)), 1)
+ IK10 = MAX ( COUNT (TPFLYER%XALTLAUNCH >= ZZM(2,1,:)), 1)
+ IK11 = MAX ( COUNT (TPFLYER%XALTLAUNCH >= ZZM(2,2,:)), 1)
+ ZZCOEF00 = (TPFLYER%XALTLAUNCH - ZZM(1,1,IK00)) / ( ZZM(1,1,IK00+1) - ZZM(1,1,IK00))
+ ZZCOEF01 = (TPFLYER%XALTLAUNCH - ZZM(1,2,IK01)) / ( ZZM(1,2,IK01+1) - ZZM(1,2,IK01))
ZZCOEF10 = (TPFLYER%XALTLAUNCH - ZZM(2,1,IK10)) / ( ZZM(2,1,IK10+1) - ZZM(2,1,IK10))
ZZCOEF11 = (TPFLYER%XALTLAUNCH - ZZM(2,2,IK11)) / ( ZZM(2,2,IK11+1) - ZZM(2,2,IK11))
TPFLYER%XRHO = FLYER_INTERP(ZRHO)
@@ -656,22 +1213,22 @@ IF ( TPFLYER%LFLY ) THEN
CMNHMSG(3) = 'Check your INI_BALLOON routine'
CALL PRINT_MSG( NVERB_FATAL, 'GEN', 'AIRCRAFT_BALLOON_EVOL', OLOCAL = .TRUE. )
END IF
-!
-!* 5.2.2 Radiosounding balloon
-!
+ !
+ !* 5.2.2 Radiosounding balloon
+ !
CASE ( 'RADIOS' )
TPFLYER%XZ_CUR = TPFLYER%XALTLAUNCH
TPFLYER%XZ_CUR = MAX ( TPFLYER%XZ_CUR , ZZM(1,1,IKB) )
TPFLYER%XZ_CUR = MAX ( TPFLYER%XZ_CUR , ZZM(2,1,IKB) )
TPFLYER%XZ_CUR = MAX ( TPFLYER%XZ_CUR , ZZM(1,2,IKB) )
TPFLYER%XZ_CUR = MAX ( TPFLYER%XZ_CUR , ZZM(2,2,IKB) )
-!
-!* 5.2.4 Constant Volume Balloon
-!
+ !
+ !* 5.2.4 Constant Volume Balloon
+ !
CASE ( 'CVBALL' )
- ZXCOEF = (TPFLYER%XX_CUR - XXHATM(II)) / (XXHATM(II+1) - XXHATM(II))
+ ZXCOEF = (TPFLYER%XX_CUR - XXHATM(II_M)) / (XXHATM(II_M+1) - XXHATM(II_M))
ZXCOEF = MAX (0.,MIN(ZXCOEF,1.))
- ZYCOEF = (TPFLYER%XY_CUR - XYHATM(IJ)) / (XYHATM(IJ+1) - XYHATM(IJ))
+ ZYCOEF = (TPFLYER%XY_CUR - XYHATM(IJ_M)) / (XYHATM(IJ_M+1) - XYHATM(IJ_M))
ZYCOEF = MAX (0.,MIN(ZYCOEF,1.))
IF ( TPFLYER%XALTLAUNCH /= XNEGUNDEF ) THEN
IK00 = MAX ( COUNT (TPFLYER%XALTLAUNCH >= ZZM(1,1,:)), 1)
@@ -731,27 +1288,13 @@ IF ( TPFLYER%LFLY ) THEN
END IF
END IF
END SELECT
+ END IF LAUNCHVERTPOS
+ !
+ !
+ !
+ !* 5.3 Vertical position
+ ! -----------------
!
-!* 5.2.3 Aircraft
-!
- CLASS IS ( TAIRCRAFTDATA)
- IF (TPFLYER%LALTDEF) THEN
- TPFLYER%XP_CUR = (1.-ZSEG_FRAC) * TPFLYER%XSEGP(IL ) &
- + ZSEG_FRAC * TPFLYER%XSEGP(IL+1)
- ELSE
- TPFLYER%XZ_CUR = (1.-ZSEG_FRAC) * TPFLYER%XSEGZ(IL ) &
- + ZSEG_FRAC * TPFLYER%XSEGZ(IL +1 )
- END IF
- END SELECT
- END IF
-!
-!
-!
-!* 5.3 Vertical position
-! -----------------
-!
- SELECT TYPE ( TPFLYER )
- CLASS IS ( TBALLOONDATA)
IF ( TPFLYER%CTYPE == 'ISODEN' ) THEN
IK00 = MAX ( COUNT (TPFLYER%XRHO <= ZRHO(1,1,:)), 1)
IK01 = MAX ( COUNT (TPFLYER%XRHO <= ZRHO(1,2,:)), 1)
@@ -764,88 +1307,63 @@ IF ( TPFLYER%LFLY ) THEN
IK11 = MAX ( COUNT (TPFLYER%XZ_CUR >= ZZM(2,2,:)), 1)
END IF
- CLASS IS ( TAIRCRAFTDATA)
- IF ( TPFLYER%LALTDEF ) THEN
- ZFLYER_EXN = (TPFLYER%XP_CUR/XP00)**(XRD/XCPD)
- IK00 = MAX ( COUNT (ZFLYER_EXN <= ZEXN(1,1,:)), 1)
- IK01 = MAX ( COUNT (ZFLYER_EXN <= ZEXN(1,2,:)), 1)
- IK10 = MAX ( COUNT (ZFLYER_EXN <= ZEXN(2,1,:)), 1)
- IK11 = MAX ( COUNT (ZFLYER_EXN <= ZEXN(2,2,:)), 1)
- ELSE
- IK00 = MAX ( COUNT (TPFLYER%XZ_CUR >= ZZM(1,1,:)), 1)
- IK01 = MAX ( COUNT (TPFLYER%XZ_CUR >= ZZM(1,2,:)), 1)
- IK10 = MAX ( COUNT (TPFLYER%XZ_CUR >= ZZM(2,1,:)), 1)
- IK11 = MAX ( COUNT (TPFLYER%XZ_CUR >= ZZM(2,2,:)), 1)
+ IF ( ANY( [ IK00, IK01, IK10, IK11 ] < IKB ) ) THEN
+ !Minimum altitude is on the ground at IKB (no crash if too low)
+ IK00 = MAX ( IK00, IKB )
+ IK01 = MAX ( IK01, IKB )
+ IK10 = MAX ( IK10, IKB )
+ IK11 = MAX ( IK11, IKB )
+
+ CMNHMSG(1) = 'flyer ' // TRIM( TPFLYER%CTITLE ) // ' is near the ground'
+ WRITE( CMNHMSG(2), "( 'at ', I2, '/', I2, '/', I4, ' ', F18.12, 's' )" ) &
+ TDTCUR%NDAY, TDTCUR%NMONTH, TDTCUR%NYEAR, TDTCUR%XTIME
+ CALL PRINT_MSG( NVERB_INFO, 'GEN', 'AIRCRAFT_BALLOON_EVOL' , OLOCAL = .TRUE.)
END IF
- END SELECT
- IK00 = MAX ( IK00, IKB )
- IK01 = MAX ( IK01, IKB )
- IK10 = MAX ( IK10, IKB )
- IK11 = MAX ( IK11, IKB )
-!
-!
-!* 5.4 Crash of the balloon
-! --------------------
-!
-!
- IF (IK00 < IKB .OR. IK01 < IKB .OR. IK10 < IKB .OR. IK11 < IKB .OR. &
- IK00 >= IKE .OR. IK01 >= IKE .OR. IK10 >= IKE .OR. IK11 >= IKE ) THEN
- TPFLYER%LCRASH = .TRUE.
- END IF
-!
- END IF
-!
-!
- IF ( TPFLYER%LCRASH ) THEN
- TPFLYER%LFLY = .FALSE.
- IF ( TPFLYER%CTYPE == 'AIRCRA' .AND. .NOT. GLAUNCH ) THEN
- WRITE(ILUOUT,*) 'Aircraft ',TPFLYER%CTITLE,' flew out of the domain the ', &
- TDTCUR%nday,'/',TDTCUR%nmonth,'/', &
- TDTCUR%nyear,' at ',TDTCUR%xtime,' sec.'
- ELSE IF (TPFLYER%CTYPE /= 'AIRCRA') THEN
- WRITE(ILUOUT,*) 'Balloon ',TPFLYER%CTITLE,' crashed the ', &
- TDTCUR%nday,'/',TDTCUR%nmonth,'/', &
- TDTCUR%nyear,' at ',TDTCUR%xtime,' sec.'
- END IF
- ELSE
- SELECT TYPE ( TPFLYER )
- CLASS IS ( TAIRCRAFTDATA)
- IF ( .NOT. GLAUNCH .AND. ZTDIST > PTSTEP ) THEN
- WRITE(ILUOUT,*) '-------------------------------------------------------------------'
- WRITE(ILUOUT,*) 'Aircraft ',TPFLYER%CTITLE,' flies in leg',TPFLYER%NSEGCURN ,' the ', &
- TDTCUR%nday,'/',TDTCUR%nmonth,'/', &
- TDTCUR%nyear,' at ',NINT(TDTCUR%xtime),' sec.'
- WRITE(ILUOUT,*) '-------------------------------------------------------------------'
- ENDIF
- END SELECT
-!
-!----------------------------------------------------------------------------
- IF (ZTHIS_PROC>0.) THEN
-!----------------------------------------------------------------------------
-!
-!* 6. INITIALIZATIONS FOR INTERPOLATIONS
-! ----------------------------------
-!
-!* 6.1 Interpolation coefficient for X
-! -------------------------------
-!
- ZXCOEF = (TPFLYER%XX_CUR - XXHATM(II)) / (XXHATM(II+1) - XXHATM(II))
- ZXCOEF = MAX (0.,MIN(ZXCOEF,1.))
-!
-!
-!* 6.2 Interpolation coefficient for y
-! -------------------------------
-!
- ZYCOEF = (TPFLYER%XY_CUR - XYHATM(IJ)) / (XYHATM(IJ+1) - XYHATM(IJ))
- ZYCOEF = MAX (0.,MIN(ZYCOEF,1.))
-!
-!
-!* 6.3 Interpolation coefficients for the 4 suroundings verticals
-! ----------------------------------------------------------
-!
- SELECT TYPE ( TPFLYER )
- CLASS IS ( TBALLOONDATA)
+ !
+ !
+ !* 5.4 Crash of the balloon
+ ! --------------------
+ !
+ !
+ IF (IK00 < IKB .OR. IK01 < IKB .OR. IK10 < IKB .OR. IK11 < IKB ) THEN
+ TPFLYER%LCRASH = .TRUE.
+ TPFLYER%NCRASH = NCRASH_OUT_LOW
+ END IF
+ IF (IK00 >= IKE .OR. IK01 >= IKE .OR. IK10 >= IKE .OR. IK11 >= IKE ) THEN
+ TPFLYER%LCRASH = .TRUE.
+ TPFLYER%NCRASH = NCRASH_OUT_HIGH
+ END IF
+ CRASH_VERT: IF ( TPFLYER%LCRASH ) THEN
+ TPFLYER%LFLY = .FALSE.
+ WRITE( CMNHMSG(1), "( 'Balloon ', A, ' crashed the ', I2, '/', I2, '/', I4, ' at ', F18.12, 's (too low or too high)' )" ) &
+ TRIM( TPFLYER%CTITLE ), TDTCUR%NDAY, TDTCUR%NMONTH, TDTCUR%NYEAR, TDTCUR%XTIME
+ CALL PRINT_MSG( NVERB_WARNING, 'GEN', 'AIRCRAFT_BALLOON_EVOL', OLOCAL = .TRUE. )
+ ELSE CRASH_VERT
+ !No vertical crash
+ !
+ !* 6. INITIALIZATIONS FOR INTERPOLATIONS
+ ! ----------------------------------
+ !
+ !* 6.1 Interpolation coefficient for X
+ ! -------------------------------
+ !
+ ZXCOEF = (TPFLYER%XX_CUR - XXHATM(II_M)) / (XXHATM(II_M+1) - XXHATM(II_M))
+ ZXCOEF = MAX (0.,MIN(ZXCOEF,1.))
+ !
+ !
+ !* 6.2 Interpolation coefficient for y
+ ! -------------------------------
+ !
+ ZYCOEF = (TPFLYER%XY_CUR - XYHATM(IJ_M)) / (XYHATM(IJ_M+1) - XYHATM(IJ_M))
+ ZYCOEF = MAX (0.,MIN(ZYCOEF,1.))
+ !
+ !
+ !* 6.3 Interpolation coefficients for the 4 suroundings verticals
+ ! ----------------------------------------------------------
+ !
+ ! SELECT TYPE ( TPFLYER )
+ ! CLASS IS ( TBALLOONDATA)
IF ( TPFLYER%CTYPE == 'ISODEN' ) THEN
ZZCOEF00 = (TPFLYER%XRHO - ZRHO(1,1,IK00)) / ( ZRHO(1,1,IK00+1) - ZRHO(1,1,IK00) )
ZZCOEF01 = (TPFLYER%XRHO - ZRHO(1,2,IK01)) / ( ZRHO(1,2,IK01+1) - ZRHO(1,2,IK01) )
@@ -858,431 +1376,433 @@ IF ( TPFLYER%LFLY ) THEN
ZZCOEF10 = (TPFLYER%XZ_CUR - ZZM(2,1,IK10)) / ( ZZM(2,1,IK10+1) - ZZM(2,1,IK10) )
ZZCOEF11 = (TPFLYER%XZ_CUR - ZZM(2,2,IK11)) / ( ZZM(2,2,IK11+1) - ZZM(2,2,IK11) )
END IF
+ !----------------------------------------------------------------------------
+ !
+ !* 7. INITIALIZATIONS FOR INTERPOLATIONS OF U AND V
+ ! ---------------------------------------------
+ !
+ !* 7.1 Interpolation coefficient for X (for U)
+ ! -------------------------------
+ !
+ ZUCOEF = (TPFLYER%XX_CUR - XXHAT(II_U)) / (XXHAT(II_U+1) - XXHAT(II_U))
+ ZUCOEF = MAX(0.,MIN(ZUCOEF,1.))
+ !
+ !
+ !* 7.2 Interpolation coefficient for y (for V)
+ ! -------------------------------
+ !
+ ZVCOEF = (TPFLYER%XY_CUR - XYHAT(IJ_V)) / (XYHAT(IJ_V+1) - XYHAT(IJ_V))
+ ZVCOEF = MAX(0.,MIN(ZVCOEF,1.))
+ !
+ !
+ !* 7.3 Interpolation coefficients for the 4 suroundings verticals (for U)
+ ! ----------------------------------------------------------
+ !
+ IU00 = MAX( COUNT (TPFLYER%XZ_CUR >= ZZU(1,1,:)), 1)
+ IU01 = MAX( COUNT (TPFLYER%XZ_CUR >= ZZU(1,2,:)), 1)
+ IU10 = MAX( COUNT (TPFLYER%XZ_CUR >= ZZU(2,1,:)), 1)
+ IU11 = MAX( COUNT (TPFLYER%XZ_CUR >= ZZU(2,2,:)), 1)
+ ZUCOEF00 = (TPFLYER%XZ_CUR - ZZU(1,1,IU00)) / ( ZZU(1,1,IU00+1) - ZZU(1,1,IU00) )
+ ZUCOEF01 = (TPFLYER%XZ_CUR - ZZU(1,2,IU01)) / ( ZZU(1,2,IU01+1) - ZZU(1,2,IU01) )
+ ZUCOEF10 = (TPFLYER%XZ_CUR - ZZU(2,1,IU10)) / ( ZZU(2,1,IU10+1) - ZZU(2,1,IU10) )
+ ZUCOEF11 = (TPFLYER%XZ_CUR - ZZU(2,2,IU11)) / ( ZZU(2,2,IU11+1) - ZZU(2,2,IU11) )
+ !
+ !
+ !* 7.4 Interpolation coefficients for the 4 suroundings verticals (for V)
+ ! ----------------------------------------------------------
+ !
+ IV00 = MAX ( COUNT (TPFLYER%XZ_CUR >= ZZV(1,1,:)), 1)
+ IV01 = MAX ( COUNT (TPFLYER%XZ_CUR >= ZZV(1,2,:)), 1)
+ IV10 = MAX ( COUNT (TPFLYER%XZ_CUR >= ZZV(2,1,:)), 1)
+ IV11 = MAX ( COUNT (TPFLYER%XZ_CUR >= ZZV(2,2,:)), 1)
+ ZVCOEF00 = (TPFLYER%XZ_CUR - ZZV(1,1,IV00)) / ( ZZV(1,1,IV00+1) - ZZV(1,1,IV00) )
+ ZVCOEF01 = (TPFLYER%XZ_CUR - ZZV(1,2,IV01)) / ( ZZV(1,2,IV01+1) - ZZV(1,2,IV01) )
+ ZVCOEF10 = (TPFLYER%XZ_CUR - ZZV(2,1,IV10)) / ( ZZV(2,1,IV10+1) - ZZV(2,1,IV10) )
+ ZVCOEF11 = (TPFLYER%XZ_CUR - ZZV(2,2,IV11)) / ( ZZV(2,2,IV11+1) - ZZV(2,2,IV11) )
- CLASS IS ( TAIRCRAFTDATA)
- IF ( TPFLYER%LALTDEF ) THEN
- ZZCOEF00 = (ZFLYER_EXN - ZEXN(1,1,IK00)) / ( ZEXN(1,1,IK00+1) - ZEXN(1,1,IK00) )
- ZZCOEF01 = (ZFLYER_EXN - ZEXN(1,2,IK01)) / ( ZEXN(1,2,IK01+1) - ZEXN(1,2,IK01) )
- ZZCOEF10 = (ZFLYER_EXN - ZEXN(2,1,IK10)) / ( ZEXN(2,1,IK10+1) - ZEXN(2,1,IK10) )
- ZZCOEF11 = (ZFLYER_EXN - ZEXN(2,2,IK11)) / ( ZEXN(2,2,IK11+1) - ZEXN(2,2,IK11) )
- TPFLYER%XZ_CUR = FLYER_INTERP(ZZM)
- ELSE
- ZZCOEF00 = (TPFLYER%XZ_CUR - ZZM(1,1,IK00)) / ( ZZM(1,1,IK00+1) - ZZM(1,1,IK00) )
- ZZCOEF01 = (TPFLYER%XZ_CUR - ZZM(1,2,IK01)) / ( ZZM(1,2,IK01+1) - ZZM(1,2,IK01) )
- ZZCOEF10 = (TPFLYER%XZ_CUR - ZZM(2,1,IK10)) / ( ZZM(2,1,IK10+1) - ZZM(2,1,IK10) )
- ZZCOEF11 = (TPFLYER%XZ_CUR - ZZM(2,2,IK11)) / ( ZZM(2,2,IK11+1) - ZZM(2,2,IK11) )
- TPFLYER%XP_CUR = FLYER_INTERP(PP)
- END IF
- END SELECT
-!
-!----------------------------------------------------------------------------
-!
-!* 7. INITIALIZATIONS FOR INTERPOLATIONS OF U AND V
-! ---------------------------------------------
-!
-!* 7.1 Interpolation coefficient for X (for U)
-! -------------------------------
-!
- ZUCOEF = (TPFLYER%XX_CUR - PXHAT(IU)) / (PXHAT(IU+1) - PXHAT(IU))
- ZUCOEF = MAX(0.,MIN(ZUCOEF,1.))
-!
-!
-!* 7.2 Interpolation coefficient for y (for V)
-! -------------------------------
-!
- ZVCOEF = (TPFLYER%XY_CUR - PYHAT(IV)) / (PYHAT(IV+1) - PYHAT(IV))
- ZVCOEF = MAX(0.,MIN(ZVCOEF,1.))
-!
-!
-!* 7.3 Interpolation coefficients for the 4 suroundings verticals (for U)
-! ----------------------------------------------------------
-!
- IU00 = MAX( COUNT (TPFLYER%XZ_CUR >= ZZU(1,1,:)), 1)
- IU01 = MAX( COUNT (TPFLYER%XZ_CUR >= ZZU(1,2,:)), 1)
- IU10 = MAX( COUNT (TPFLYER%XZ_CUR >= ZZU(2,1,:)), 1)
- IU11 = MAX( COUNT (TPFLYER%XZ_CUR >= ZZU(2,2,:)), 1)
- ZUCOEF00 = (TPFLYER%XZ_CUR - ZZU(1,1,IU00)) / ( ZZU(1,1,IU00+1) - ZZU(1,1,IU00) )
- ZUCOEF01 = (TPFLYER%XZ_CUR - ZZU(1,2,IU01)) / ( ZZU(1,2,IU01+1) - ZZU(1,2,IU01) )
- ZUCOEF10 = (TPFLYER%XZ_CUR - ZZU(2,1,IU10)) / ( ZZU(2,1,IU10+1) - ZZU(2,1,IU10) )
- ZUCOEF11 = (TPFLYER%XZ_CUR - ZZU(2,2,IU11)) / ( ZZU(2,2,IU11+1) - ZZU(2,2,IU11) )
-!
-!
-!* 7.4 Interpolation coefficients for the 4 suroundings verticals (for V)
-! ----------------------------------------------------------
-!
+ ! Check if it is the right moment to store data
+ IF ( TPFLYER%LSTORE ) THEN
+ ISTORE = TPFLYER%TFLYER_TIME%N_CUR
+ IF ( ABS( TDTCUR - TPFLYER%TFLYER_TIME%TPDATES(ISTORE) ) < 1e-10 ) THEN
+ !
+ !* 2. PRELIMINARIES-2
+ ! -------------
+ !
+ !* 2.1 Indices
+ ! -------
+ !
+ CALL GET_INDICE_ll (IIB,IJB,IIE,IJE)
- IV00 = MAX ( COUNT (TPFLYER%XZ_CUR >= ZZV(1,1,:)), 1)
- IV01 = MAX ( COUNT (TPFLYER%XZ_CUR >= ZZV(1,2,:)), 1)
- IV10 = MAX ( COUNT (TPFLYER%XZ_CUR >= ZZV(2,1,:)), 1)
- IV11 = MAX ( COUNT (TPFLYER%XZ_CUR >= ZZV(2,2,:)), 1)
- ZVCOEF00 = (TPFLYER%XZ_CUR - ZZV(1,1,IV00)) / ( ZZV(1,1,IV00+1) - ZZV(1,1,IV00) )
- ZVCOEF01 = (TPFLYER%XZ_CUR - ZZV(1,2,IV01)) / ( ZZV(1,2,IV01+1) - ZZV(1,2,IV01) )
- ZVCOEF10 = (TPFLYER%XZ_CUR - ZZV(2,1,IV10)) / ( ZZV(2,1,IV10+1) - ZZV(2,1,IV10) )
- ZVCOEF11 = (TPFLYER%XZ_CUR - ZZV(2,2,IV11)) / ( ZZV(2,2,IV11+1) - ZZV(2,2,IV11) )
-!
-!----------------------------------------------------------------------------
-!
-!* 8. DATA RECORDING
-! --------------
-!
- IF ( GSTORE ) THEN
- TPFLYER%XX (IN) = TPFLYER%XX_CUR
- TPFLYER%XY (IN) = TPFLYER%XY_CUR
- TPFLYER%XZ (IN) = TPFLYER%XZ_CUR
- !
- CALL SM_LATLON(PLATOR,PLONOR, &
- TPFLYER%XX_CUR, TPFLYER%XY_CUR, &
- TPFLYER%XLAT(IN), TPFLYER%XLON(IN) )
- !
- ZU_BAL = FLYER_INTERP_U(PU)
- ZV_BAL = FLYER_INTERP_V(PV)
- ZGAM = (XRPK * (TPFLYER%XLON(IN) - XLON0) - XBETA)*(XPI/180.)
- TPFLYER%XZON (IN) = ZU_BAL * COS(ZGAM) + ZV_BAL * SIN(ZGAM)
- TPFLYER%XMER (IN) = - ZU_BAL * SIN(ZGAM) + ZV_BAL * COS(ZGAM)
- !
- TPFLYER%XW (IN) = FLYER_INTERP(ZWM)
- TPFLYER%XTH (IN) = FLYER_INTERP(PTH)
- !
- ZFLYER_EXN = FLYER_INTERP(ZEXN)
- TPFLYER%XP (IN) = XP00 * ZFLYER_EXN**(XCPD/XRD)
- !
- DO JLOOP=1,SIZE(PR,4)
- TPFLYER%XR (IN,JLOOP) = FLYER_INTERP(PR(:,:,:,JLOOP))
- IF (JLOOP>=2) ZR(:,:,:) = ZR(:,:,:) + PR(:,:,:,JLOOP)
- END DO
- DO JLOOP=1,SIZE(PSV,4)
- TPFLYER%XSV (IN,JLOOP) = FLYER_INTERP(PSV(:,:,:,JLOOP))
- END DO
- TPFLYER%XRTZ (IN,:) = FLYER_INTERPZ(ZR(:,:,:))
- DO JLOOP=1,SIZE(PR,4)
- TPFLYER%XRZ (IN,:,JLOOP) = FLYER_INTERPZ(PR(:,:,:,JLOOP))
- END DO
- ! Fin Modifs ON
- TPFLYER%XFFZ (IN,:) = FLYER_INTERPZ(SQRT(PU**2+PV**2))
- IF (CCLOUD=="LIMA") THEN
- TPFLYER%XCIZ (IN,:) = FLYER_INTERPZ(PSV(:,:,:,NSV_LIMA_NI))
- TPFLYER%XCCZ (IN,:) = FLYER_INTERPZ(PSV(:,:,:,NSV_LIMA_NC))
- TPFLYER%XCRZ (IN,:) = FLYER_INTERPZ(PSV(:,:,:,NSV_LIMA_NR))
- ELSE IF ( CCLOUD=="ICE3" .OR. CCLOUD=="ICE4" ) THEN
- TPFLYER%XCIZ (IN,:) = FLYER_INTERPZ(PCIT(:,:,:))
- ENDIF
- ! initialization CRARE and CRARE_ATT + LWC and IWC
- TPFLYER%XCRARE(IN,:) = 0.
- TPFLYER%XCRARE_ATT(IN,:) = 0.
- TPFLYER%XLWCZ (IN,:) = 0.
- TPFLYER%XIWCZ (IN,:) = 0.
- IF (CCLOUD=="LIMA" .OR. CCLOUD=="ICE3" ) THEN ! only for ICE3 and LIMA
- TPFLYER%XLWCZ (IN,:) = FLYER_INTERPZ((PR(:,:,:,2)+PR(:,:,:,3))*PRHODREF(:,:,:))
- TPFLYER%XIWCZ (IN,:) = FLYER_INTERPZ((PR(:,:,:,4)+PR(:,:,:,5)+PR(:,:,:,6))*PRHODREF(:,:,:))
- ZTEMPZ(:)=FLYER_INTERPZ(PTH(II:II+1,IJ:IJ+1,:) * ZEXN(:,:,:))
- ZRHODREFZ(:)=FLYER_INTERPZ(PRHODREF(:,:,:))
- IF (CCLOUD=="LIMA") THEN
- ZCCI(:)=FLYER_INTERPZ(PSV(:,:,:,NSV_LIMA_NI))
- ZCCR(:)=FLYER_INTERPZ(PSV(:,:,:,NSV_LIMA_NR))
- ZCCC(:)=FLYER_INTERPZ(PSV(:,:,:,NSV_LIMA_NC))
- ELSE
- ZCIT(:)=FLYER_INTERPZ(PCIT(:,:,:))
- ENDIF
- DO JLOOP=3,6
- ZRZ(:,JLOOP)=FLYER_INTERPZ(PR(:,:,:,JLOOP))
- END DO
- if ( csurf == 'EXTE' ) then
- DO JK=1,IKU
- ZRZ(JK,2)=FLYER_INTERP_2D(PR(:,:,JK,2)*PSEA(:,:)) ! becomes cloud mixing ratio over sea
- ZRZ(JK,7)=FLYER_INTERP_2D(PR(:,:,JK,2)*(1.-PSEA(:,:))) ! becomes cloud mixing ratio over land
- END DO
- else
- !if csurf/='EXTE', psea is not allocated
- DO JK=1,IKU
- ZRZ(JK,2)=FLYER_INTERP_2D(PR(:,:,JK,2))
- ZRZ(JK,7) = 0.
- END DO
- end if
- ALLOCATE(ZAELOC(IKU))
- !
- ZAELOC(:)=0.
- ! initialization of quadrature points and weights
- ALLOCATE(ZX(JPTS_GAULAG),ZW(JPTS_GAULAG))
- CALL GAULAG(JPTS_GAULAG,ZX,ZW) ! for integration over diameters
- ! initialize minimum values
- ALLOCATE(ZRTMIN(SIZE(PR,4)+1))
- IF (CCLOUD == 'LIMA') THEN
- ZRTMIN(2)=XRTMIN_L(2) ! cloud water over sea
- ZRTMIN(3)=XRTMIN_L(3)
- ZRTMIN(4)=XRTMIN_L(4)
- ZRTMIN(5)=1E-10
- ZRTMIN(6)=XRTMIN_L(6)
- ZRTMIN(7)=XRTMIN_L(2) ! cloud water over land
- ELSE
- ZRTMIN(2)=XRTMIN_I(2) ! cloud water over sea
- ZRTMIN(3)=XRTMIN_I(3)
- ZRTMIN(4)=XRTMIN_I(4)
- ZRTMIN(5)=1E-10
- ZRTMIN(6)=XRTMIN_I(6)
- ZRTMIN(7)=XRTMIN_I(2) ! cloud water over land
- ENDIF
- ! compute cloud radar reflectivity from vertical profiles of temperature and mixing ratios
- DO JK=1,IKU
- QMW=SQRT(QEPSW(ZTEMPZ(JK),XLIGHTSPEED/XLAM_CRAD))
- QMI=SQRT(QEPSI(ZTEMPZ(JK),XLIGHTSPEED/XLAM_CRAD))
- DO JLOOP=2,7
- IF (CCLOUD == 'LIMA') THEN
- GCALC=(ZRZ(JK,JLOOP)>ZRTMIN(JLOOP).AND.(JLOOP.NE.4.OR.ZCCI(JK)>0.).AND.&
- (JLOOP.NE.3.OR.ZCCR(JK)>0.).AND.((JLOOP.NE.2.AND. JLOOP.NE.7).OR.ZCCC(JK)>0.))
- ELSE
- GCALC=(ZRZ(JK,JLOOP)>ZRTMIN(JLOOP).AND.(JLOOP.NE.4.OR.ZCIT(JK)>0.))
- ENDIF
- IF(GCALC) THEN
- SELECT CASE(JLOOP)
- CASE(2) ! cloud water over sea
- IF (CCLOUD == 'LIMA') THEN
- ZA=XAC_L
- ZB=XBC_L
- ZCC=ZCCC(JK)*ZRHODREFZ(JK)
- ZCX=0.
- ZALPHA=XALPHAC_L
- ZNU=XNUC_L
- ZLBEX=1.0/(ZCX-ZB)
- ZLB=( ZA*ZCC*MOMG(ZALPHA,ZNU,ZB) )**(-ZLBEX)
- ELSE
- ZA=XAC_I
- ZB=XBC_I
- ZCC=XCONC_SEA
- ZCX=0.
- ZALPHA=XALPHAC2_I
- ZNU=XNUC2_I
- ZLBEX=1.0/(ZCX-ZB)
- ZLB=( ZA*ZCC*MOMG(ZALPHA,ZNU,ZB) )**(-ZLBEX)
- ENDIF
- CASE(3) ! rain water
- IF (CCLOUD == 'LIMA') THEN
- ZA=XAR_L
- ZB=XBR_L
- ZCC=ZCCR(JK)*ZRHODREFZ(JK)
- ZCX=0.
- ZALPHA=XALPHAR_L
- ZNU=XNUR_L
- ZLBEX=1.0/(ZCX-ZB)
- ZLB=( ZA*ZCC*MOMG(ZALPHA,ZNU,ZB) )**(-ZLBEX)
- ELSE
- ZA=XAR_I
- ZB=XBR_I
- ZCC=XCCR_I
- ZCX=-1.
- ZALPHA=XALPHAR_I
- ZNU=XNUR_I
- ZLB=XLBR_I
- ZLBEX=XLBEXR_I
- ENDIF
- CASE(4) ! pristine ice
- IF (CCLOUD == 'LIMA') THEN
- ZA=XAI_L
- ZB=XBI_L
- ZCC=ZCCI(JK)*ZRHODREFZ(JK)
- ZCX=0.
- ZALPHA=XALPHAI_L
- ZNU=XNUI_L
- ZLBEX=1.0/(ZCX-ZB)
- ZLB=( ZA*ZCC*MOMG(ZALPHA,ZNU,ZB) )**(-ZLBEX) ! because ZCC not included in XLBI
- ZFW=0
- ELSE
- ZA=XAI_I
- ZB=XBI_I
- ZCC=ZCIT(JK)
- ZCX=0.
- ZALPHA=XALPHAI_I
- ZNU=XNUI_I
- ZLBEX=XLBEXI_I
- ZLB=XLBI_I*ZCC**(-ZLBEX) ! because ZCC not included in XLBI
- ZFW=0
- ENDIF
- CASE(5) ! snow
- IF (CCLOUD == 'LIMA') THEN
- ZA=XAS_L
- ZB=XBS_L
- ZCC=XCCS_L
- ZCX=XCXS_L
- ZALPHA=XALPHAS_L
- ZNU=XNUS_L
- ZLB=XLBS_L
- ZLBEX=XLBEXS_L
- ZFW=0
- ELSE
- ZA=XAS_I
- ZB=XBS_I
- ZCC=XCCS_I
- ZCX=XCXS_I
- ZALPHA=XALPHAS_I
- ZNU=XNUS_I
- ZLB=XLBS_I
- ZLBEX=XLBEXS_I
- ZFW=0
- ENDIF
- CASE(6) ! graupel
- !If temperature between -10 and 10°C and Mr and Mg over min threshold: melting graupel
- ! with liquid water fraction Fw=Mr/(Mr+Mg) else dry graupel (Fw=0)
- IF( ZTEMPZ(JK) > XTT-10 .AND. ZTEMPZ(JK) < XTT+10 &
- .AND. ZRZ(JK,3) > ZRTMIN(3) ) THEN
- ZFW=ZRZ(JK,3)/(ZRZ(JK,3)+ZRZ(JK,JLOOP))
- ELSE
- ZFW=0
- ENDIF
- IF (CCLOUD == 'LIMA') THEN
- ZA=XAG_L
- ZB=XBG_L
- ZCC=XCCG_L
- ZCX=XCXG_L
- ZALPHA=XALPHAG_L
- ZNU=XNUG_L
- ZLB=XLBG_L
- ZLBEX=XLBEXG_L
- ELSE
- ZA=XAG_I
- ZB=XBG_I
- ZCC=XCCG_I
- ZCX=XCXG_I
- ZALPHA=XALPHAG_I
- ZNU=XNUG_I
- ZLB=XLBG_I
- ZLBEX=XLBEXG_I
- ENDIF
- CASE(7) ! cloud water over land
- IF (CCLOUD == 'LIMA') THEN
- ZA=XAC_L
- ZB=XBC_L
- ZCC=ZCCC(JK)*ZRHODREFZ(JK)
- ZCX=0.
- ZALPHA=XALPHAC_L
- ZNU=XNUC_L
- ZLBEX=1.0/(ZCX-ZB)
- ZLB=( ZA*ZCC*MOMG(ZALPHA,ZNU,ZB) )**(-ZLBEX)
- ELSE
- ZA=XAC_I
- ZB=XBC_I
- ZCC=XCONC_LAND
- ZCX=0.
- ZALPHA=XALPHAC_I
- ZNU=XNUC_I
- ZLBEX=1.0/(ZCX-ZB)
- ZLB=( ZA*ZCC*MOMG(ZALPHA,ZNU,ZB) )**(-ZLBEX)
- ENDIF
- END SELECT
- ZLBDA=ZLB*(ZRHODREFZ(JK)*ZRZ(JK,JLOOP))**ZLBEX
- ZREFLOC=0.
- ZAETMP=0.
- DO JJ=1,JPTS_GAULAG ! Gauss-Laguerre quadrature
- ZDELTA_EQUIV=ZX(JJ)**(1./ZALPHA)/ZLBDA
- SELECT CASE(JLOOP)
- CASE(2,3,7)
- QM=QMW
- CASE(4,5,6)
- ! pristine ice, snow, dry graupel
- ZRHOHYD=MIN(6.*ZA*ZDELTA_EQUIV**(ZB-3.)/XPI,.92*XRHOLW)
- QM=sqrt(MG(QMI**2,CMPLX(1,0),ZRHOHYD/.92/XRHOLW))
-
- ! water inclusions in ice in air
- QEPSWI=MG(QMW**2,QM**2,ZFW)
- ! ice in air inclusions in water
- QEPSIW=MG(QM**2,QMW**2,1.-ZFW)
-
- !MG weighted rule (Matrosov 2008)
- IF(ZFW .LT. 0.37) THEN
- ZFPW=0
- ELSE IF(ZFW .GT. 0.63) THEN
- ZFPW=1
- ELSE
- ZFPW=(ZFW-0.37)/(0.63-0.37)
- ENDIF
- QM=sqrt(QEPSWI*(1.-ZFPW)+QEPSIW*ZFPW)
- END SELECT
- CALL BHMIE(XPI/XLAM_CRAD*ZDELTA_EQUIV,QM,ZQEXT,ZQSCA,ZQBACK)
- ZREFLOC=ZREFLOC+ZQBACK*ZX(JJ)**(ZNU-1)*ZDELTA_EQUIV**2*ZW(JJ)
- ZAETMP =ZAETMP +ZQEXT *ZX(JJ)**(ZNU-1)*ZDELTA_EQUIV**2*ZW(JJ)
+ IKB = 1 + JPVEXT
+ IKE = SIZE(PZ,3) - JPVEXT
+ !
+ !
+ !* 2.2 Interpolations of model variables to mass points
+ ! ------------------------------------------------
+ !
+ IIU=SIZE(XXHAT)
+ IJU=SIZE(XYHAT)
+
+ TPFLYER%NMODELHIST(ISTORE) = TPFLYER%NMODEL
+ TPFLYER%XX (ISTORE) = TPFLYER%XX_CUR
+ TPFLYER%XY (ISTORE) = TPFLYER%XY_CUR
+ TPFLYER%XZ (ISTORE) = TPFLYER%XZ_CUR
+ !
+ CALL SM_LATLON(PLATOR,PLONOR, &
+ TPFLYER%XX_CUR, TPFLYER%XY_CUR, &
+ TPFLYER%XLAT(ISTORE), TPFLYER%XLON(ISTORE) )
+ !
+ ZU_BAL = FLYER_INTERP_U(PU)
+ ZV_BAL = FLYER_INTERP_V(PV)
+ ZGAM = (XRPK * (TPFLYER%XLON(ISTORE) - XLON0) - XBETA)*(XPI/180.)
+ TPFLYER%XZON (ISTORE) = ZU_BAL * COS(ZGAM) + ZV_BAL * SIN(ZGAM)
+ TPFLYER%XMER (ISTORE) = - ZU_BAL * SIN(ZGAM) + ZV_BAL * COS(ZGAM)
+ !
+ TPFLYER%XW (ISTORE) = FLYER_INTERP(ZWM)
+ TPFLYER%XTH (ISTORE) = FLYER_INTERP(PTH)
+ !
+ ZFLYER_EXN = FLYER_INTERP(ZEXN)
+ TPFLYER%XP (ISTORE) = XP00 * ZFLYER_EXN**(XCPD/XRD)
+ !
+ ZR(:,:,:) = 0.
+ DO JLOOP=1,SIZE(PR,4)
+ TPFLYER%XR (ISTORE,JLOOP) = FLYER_INTERP(PR(:,:,:,JLOOP))
+ IF (JLOOP>=2) ZR(:,:,:) = ZR(:,:,:) + PR(:,:,:,JLOOP)
END DO
- ZREFLOC=ZREFLOC*(XLAM_CRAD/XPI)**4*ZCC*ZLBDA**ZCX/(4.*GAMMA(ZNU)*.93)
- ZAETMP=ZAETMP * XPI *ZCC*ZLBDA**ZCX/(4.*GAMMA(ZNU))
- TPFLYER%XCRARE(IN,JK)=TPFLYER%XCRARE(IN,JK)+ZREFLOC
- ZAELOC(JK)=ZAELOC(JK)+ZAETMP
- END IF
+ DO JLOOP=1,SIZE(PSV,4)
+ TPFLYER%XSV (ISTORE,JLOOP) = FLYER_INTERP(PSV(:,:,:,JLOOP))
+ END DO
+ TPFLYER%XRTZ (ISTORE,:) = FLYER_INTERPZ(ZR(:,:,:))
+ DO JLOOP=1,SIZE(PR,4)
+ TPFLYER%XRZ (ISTORE,:,JLOOP) = FLYER_INTERPZ(PR(:,:,:,JLOOP))
+ END DO
+ ! Fin Modifs ON
+ TPFLYER%XFFZ (ISTORE,:) = FLYER_INTERPZ(SQRT(PU**2+PV**2))
+ IF (CCLOUD=="LIMA") THEN
+ TPFLYER%XCIZ (ISTORE,:) = FLYER_INTERPZ(PSV(:,:,:,NSV_LIMA_NI))
+ TPFLYER%XCCZ (ISTORE,:) = FLYER_INTERPZ(PSV(:,:,:,NSV_LIMA_NC))
+ TPFLYER%XCRZ (ISTORE,:) = FLYER_INTERPZ(PSV(:,:,:,NSV_LIMA_NR))
+ ELSE IF ( CCLOUD=="ICE3" .OR. CCLOUD=="ICE4" ) THEN
+ TPFLYER%XCIZ (ISTORE,:) = FLYER_INTERPZ(PCIT(:,:,:))
+ ENDIF
+ ! initialization CRARE and CRARE_ATT + LWC and IWC
+ TPFLYER%XCRARE(ISTORE,:) = 0.
+ TPFLYER%XCRARE_ATT(ISTORE,:) = 0.
+ TPFLYER%XLWCZ (ISTORE,:) = 0.
+ TPFLYER%XIWCZ (ISTORE,:) = 0.
+ IF (CCLOUD=="LIMA" .OR. CCLOUD=="ICE3" ) THEN ! only for ICE3 and LIMA
+ TPFLYER%XLWCZ (ISTORE,:) = FLYER_INTERPZ((PR(:,:,:,2)+PR(:,:,:,3))*PRHODREF(:,:,:))
+ TPFLYER%XIWCZ (ISTORE,:) = FLYER_INTERPZ((PR(:,:,:,4)+PR(:,:,:,5)+PR(:,:,:,6))*PRHODREF(:,:,:))
+ ZTEMPZ(:)=FLYER_INTERPZ(PTH(II_M:II_M+1,IJ_M:IJ_M+1,:) * ZEXN(:,:,:))
+ ZRHODREFZ(:)=FLYER_INTERPZ(PRHODREF(:,:,:))
+ IF (CCLOUD=="LIMA") THEN
+ ZCCI(:)=FLYER_INTERPZ(PSV(:,:,:,NSV_LIMA_NI))
+ ZCCR(:)=FLYER_INTERPZ(PSV(:,:,:,NSV_LIMA_NR))
+ ZCCC(:)=FLYER_INTERPZ(PSV(:,:,:,NSV_LIMA_NC))
+ ELSE
+ ZCIT(:)=FLYER_INTERPZ(PCIT(:,:,:))
+ ENDIF
+ DO JLOOP=3,6
+ ZRZ(:,JLOOP)=FLYER_INTERPZ(PR(:,:,:,JLOOP))
+ END DO
+ if ( csurf == 'EXTE' ) then
+ DO JK=1,IKU
+ ZRZ(JK,2)=FLYER_INTERP_2D(PR(:,:,JK,2)*PSEA(:,:)) ! becomes cloud mixing ratio over sea
+ ZRZ(JK,7)=FLYER_INTERP_2D(PR(:,:,JK,2)*(1.-PSEA(:,:))) ! becomes cloud mixing ratio over land
+ END DO
+ else
+ !if csurf/='EXTE', psea is not allocated
+ DO JK=1,IKU
+ ZRZ(JK,2)=FLYER_INTERP_2D(PR(:,:,JK,2))
+ ZRZ(JK,7) = 0.
+ END DO
+ end if
+ ALLOCATE(ZAELOC(IKU))
+ !
+ ZAELOC(:)=0.
+ ! initialization of quadrature points and weights
+ ALLOCATE(ZX(JPTS_GAULAG),ZW(JPTS_GAULAG))
+ CALL GAULAG(JPTS_GAULAG,ZX,ZW) ! for integration over diameters
+ ! initialize minimum values
+ ALLOCATE(ZRTMIN(SIZE(PR,4)+1))
+ IF (CCLOUD == 'LIMA') THEN
+ ZRTMIN(2)=XRTMIN_L(2) ! cloud water over sea
+ ZRTMIN(3)=XRTMIN_L(3)
+ ZRTMIN(4)=XRTMIN_L(4)
+ ZRTMIN(5)=1E-10
+ ZRTMIN(6)=XRTMIN_L(6)
+ ZRTMIN(7)=XRTMIN_L(2) ! cloud water over land
+ ELSE
+ ZRTMIN(2)=XRTMIN_I(2) ! cloud water over sea
+ ZRTMIN(3)=XRTMIN_I(3)
+ ZRTMIN(4)=XRTMIN_I(4)
+ ZRTMIN(5)=1E-10
+ ZRTMIN(6)=XRTMIN_I(6)
+ ZRTMIN(7)=XRTMIN_I(2) ! cloud water over land
+ ENDIF
+ ! compute cloud radar reflectivity from vertical profiles of temperature and mixing ratios
+ DO JK=1,IKU
+ QMW=SQRT(QEPSW(ZTEMPZ(JK),XLIGHTSPEED/XLAM_CRAD))
+ QMI=SQRT(QEPSI(ZTEMPZ(JK),XLIGHTSPEED/XLAM_CRAD))
+ DO JLOOP=2,7
+ IF (CCLOUD == 'LIMA') THEN
+ GCALC=(ZRZ(JK,JLOOP)>ZRTMIN(JLOOP).AND.(JLOOP.NE.4.OR.ZCCI(JK)>0.).AND.&
+ (JLOOP.NE.3.OR.ZCCR(JK)>0.).AND.((JLOOP.NE.2.AND. JLOOP.NE.7).OR.ZCCC(JK)>0.))
+ ELSE
+ GCALC=(ZRZ(JK,JLOOP)>ZRTMIN(JLOOP).AND.(JLOOP.NE.4.OR.ZCIT(JK)>0.))
+ ENDIF
+ IF(GCALC) THEN
+ SELECT CASE(JLOOP)
+ CASE(2) ! cloud water over sea
+ IF (CCLOUD == 'LIMA') THEN
+ ZA=XAC_L
+ ZB=XBC_L
+ ZCC=ZCCC(JK)*ZRHODREFZ(JK)
+ ZCX=0.
+ ZALPHA=XALPHAC_L
+ ZNU=XNUC_L
+ ZLBEX=1.0/(ZCX-ZB)
+ ZLB=( ZA*ZCC*MOMG(ZALPHA,ZNU,ZB) )**(-ZLBEX)
+ ELSE
+ ZA=XAC_I
+ ZB=XBC_I
+ ZCC=XCONC_SEA
+ ZCX=0.
+ ZALPHA=XALPHAC2_I
+ ZNU=XNUC2_I
+ ZLBEX=1.0/(ZCX-ZB)
+ ZLB=( ZA*ZCC*MOMG(ZALPHA,ZNU,ZB) )**(-ZLBEX)
+ ENDIF
+ CASE(3) ! rain water
+ IF (CCLOUD == 'LIMA') THEN
+ ZA=XAR_L
+ ZB=XBR_L
+ ZCC=ZCCR(JK)*ZRHODREFZ(JK)
+ ZCX=0.
+ ZALPHA=XALPHAR_L
+ ZNU=XNUR_L
+ ZLBEX=1.0/(ZCX-ZB)
+ ZLB=( ZA*ZCC*MOMG(ZALPHA,ZNU,ZB) )**(-ZLBEX)
+ ELSE
+ ZA=XAR_I
+ ZB=XBR_I
+ ZCC=XCCR_I
+ ZCX=-1.
+ ZALPHA=XALPHAR_I
+ ZNU=XNUR_I
+ ZLB=XLBR_I
+ ZLBEX=XLBEXR_I
+ ENDIF
+ CASE(4) ! pristine ice
+ IF (CCLOUD == 'LIMA') THEN
+ ZA=XAI_L
+ ZB=XBI_L
+ ZCC=ZCCI(JK)*ZRHODREFZ(JK)
+ ZCX=0.
+ ZALPHA=XALPHAI_L
+ ZNU=XNUI_L
+ ZLBEX=1.0/(ZCX-ZB)
+ ZLB=( ZA*ZCC*MOMG(ZALPHA,ZNU,ZB) )**(-ZLBEX) ! because ZCC not included in XLBI
+ ZFW=0
+ ELSE
+ ZA=XAI_I
+ ZB=XBI_I
+ ZCC=ZCIT(JK)
+ ZCX=0.
+ ZALPHA=XALPHAI_I
+ ZNU=XNUI_I
+ ZLBEX=XLBEXI_I
+ ZLB=XLBI_I*ZCC**(-ZLBEX) ! because ZCC not included in XLBI
+ ZFW=0
+ ENDIF
+ CASE(5) ! snow
+ IF (CCLOUD == 'LIMA') THEN
+ ZA=XAS_L
+ ZB=XBS_L
+ ZCC=XCCS_L
+ ZCX=XCXS_L
+ ZALPHA=XALPHAS_L
+ ZNU=XNUS_L
+ ZLB=XLBS_L
+ ZLBEX=XLBEXS_L
+ ZFW=0
+ ELSE
+ ZA=XAS_I
+ ZB=XBS_I
+ ZCC=XCCS_I
+ ZCX=XCXS_I
+ ZALPHA=XALPHAS_I
+ ZNU=XNUS_I
+ ZLB=XLBS_I
+ ZLBEX=XLBEXS_I
+ ZFW=0
+ ENDIF
+ CASE(6) ! graupel
+ !If temperature between -10 and 10°C and Mr and Mg over min threshold: melting graupel
+ ! with liquid water fraction Fw=Mr/(Mr+Mg) else dry graupel (Fw=0)
+ IF( ZTEMPZ(JK) > XTT-10 .AND. ZTEMPZ(JK) < XTT+10 &
+ .AND. ZRZ(JK,3) > ZRTMIN(3) ) THEN
+ ZFW=ZRZ(JK,3)/(ZRZ(JK,3)+ZRZ(JK,JLOOP))
+ ELSE
+ ZFW=0
+ ENDIF
+ IF (CCLOUD == 'LIMA') THEN
+ ZA=XAG_L
+ ZB=XBG_L
+ ZCC=XCCG_L
+ ZCX=XCXG_L
+ ZALPHA=XALPHAG_L
+ ZNU=XNUG_L
+ ZLB=XLBG_L
+ ZLBEX=XLBEXG_L
+ ELSE
+ ZA=XAG_I
+ ZB=XBG_I
+ ZCC=XCCG_I
+ ZCX=XCXG_I
+ ZALPHA=XALPHAG_I
+ ZNU=XNUG_I
+ ZLB=XLBG_I
+ ZLBEX=XLBEXG_I
+ ENDIF
+ CASE(7) ! cloud water over land
+ IF (CCLOUD == 'LIMA') THEN
+ ZA=XAC_L
+ ZB=XBC_L
+ ZCC=ZCCC(JK)*ZRHODREFZ(JK)
+ ZCX=0.
+ ZALPHA=XALPHAC_L
+ ZNU=XNUC_L
+ ZLBEX=1.0/(ZCX-ZB)
+ ZLB=( ZA*ZCC*MOMG(ZALPHA,ZNU,ZB) )**(-ZLBEX)
+ ELSE
+ ZA=XAC_I
+ ZB=XBC_I
+ ZCC=XCONC_LAND
+ ZCX=0.
+ ZALPHA=XALPHAC_I
+ ZNU=XNUC_I
+ ZLBEX=1.0/(ZCX-ZB)
+ ZLB=( ZA*ZCC*MOMG(ZALPHA,ZNU,ZB) )**(-ZLBEX)
+ ENDIF
+ END SELECT
+ ZLBDA=ZLB*(ZRHODREFZ(JK)*ZRZ(JK,JLOOP))**ZLBEX
+ ZREFLOC=0.
+ ZAETMP=0.
+ DO JJ=1,JPTS_GAULAG ! Gauss-Laguerre quadrature
+ ZDELTA_EQUIV=ZX(JJ)**(1./ZALPHA)/ZLBDA
+ SELECT CASE(JLOOP)
+ CASE(2,3,7)
+ QM=QMW
+ CASE(4,5,6)
+ ! pristine ice, snow, dry graupel
+ ZRHOHYD=MIN(6.*ZA*ZDELTA_EQUIV**(ZB-3.)/XPI,.92*XRHOLW)
+ QM=sqrt(MG(QMI**2,CMPLX(1,0),ZRHOHYD/.92/XRHOLW))
- END DO
+ ! water inclusions in ice in air
+ QEPSWI=MG(QMW**2,QM**2,ZFW)
+ ! ice in air inclusions in water
+ QEPSIW=MG(QM**2,QMW**2,1.-ZFW)
- END DO
+ !MG weighted rule (Matrosov 2008)
+ IF(ZFW .LT. 0.37) THEN
+ ZFPW=0
+ ELSE IF(ZFW .GT. 0.63) THEN
+ ZFPW=1
+ ELSE
+ ZFPW=(ZFW-0.37)/(0.63-0.37)
+ ENDIF
+ QM=sqrt(QEPSWI*(1.-ZFPW)+QEPSIW*ZFPW)
+ END SELECT
+ CALL BHMIE(XPI/XLAM_CRAD*ZDELTA_EQUIV,QM,ZQEXT,ZQSCA,ZQBACK)
+ ZREFLOC=ZREFLOC+ZQBACK*ZX(JJ)**(ZNU-1)*ZDELTA_EQUIV**2*ZW(JJ)
+ ZAETMP =ZAETMP +ZQEXT *ZX(JJ)**(ZNU-1)*ZDELTA_EQUIV**2*ZW(JJ)
+ END DO
+ ZREFLOC=ZREFLOC*(XLAM_CRAD/XPI)**4*ZCC*ZLBDA**ZCX/(4.*GAMMA(ZNU)*.93)
+ ZAETMP=ZAETMP * XPI *ZCC*ZLBDA**ZCX/(4.*GAMMA(ZNU))
+ TPFLYER%XCRARE(ISTORE,JK)=TPFLYER%XCRARE(ISTORE,JK)+ZREFLOC
+ ZAELOC(JK)=ZAELOC(JK)+ZAETMP
+ END IF
+ END DO
+ END DO
- ! apply attenuation
- ALLOCATE(ZZMZ(IKU))
- ZZMZ(:)=FLYER_INTERPZ(ZZM(:,:,:))
- ! nadir
- ZAETOT=1.
- DO JK=COUNT(TPFLYER%XZ_CUR >= ZZMZ(:)),1,-1
- IF(JK.EQ.COUNT(TPFLYER%XZ_CUR >= ZZMZ(:))) THEN
- IF(TPFLYER%XZ_CUR<=ZZMZ(JK)+.5*(ZZMZ(JK+1)-ZZMZ(JK))) THEN
- ! only attenuation from ZAELOC(JK)
- ZAETOT=ZAETOT*EXP(-2.*(ZAELOC(JK)*(TPFLYER%XZ_CUR-ZZMZ(JK))))
- ELSE
- ! attenuation from ZAELOC(JK) and ZAELOC(JK+1)
- ZAETOT=ZAETOT*EXP(-2.*(ZAELOC(JK+1)*(TPFLYER%XZ_CUR-.5*(ZZMZ(JK+1)+ZZMZ(JK))) &
- +ZAELOC(JK)*.5*(ZZMZ(JK+1)-ZZMZ(JK))))
- END IF
- ELSE
- ! attenuation from ZAELOC(JK) and ZAELOC(JK+1)
- ZAETOT=ZAETOT*EXP(-(ZAELOC(JK+1)+ZAELOC(JK))*(ZZMZ(JK+1)-ZZMZ(JK)))
- END IF
- TPFLYER%XCRARE_ATT(IN,JK)=TPFLYER%XCRARE(IN,JK)*ZAETOT
- END DO
- ! zenith
- ZAETOT=1.
- DO JK = MAX(COUNT(TPFLYER%XZ_CUR >= ZZMZ(:)),1)+1,IKU
- IF ( JK .EQ. (MAX(COUNT(TPFLYER%XZ_CUR >= ZZMZ(:)),1)+1) ) THEN
- IF(TPFLYER%XZ_CUR>=ZZMZ(JK)-.5*(ZZMZ(JK)-ZZMZ(JK-1))) THEN
- ! only attenuation from ZAELOC(JK)
- ZAETOT=ZAETOT*EXP(-2.*(ZAELOC(JK)*(ZZMZ(JK)-TPFLYER%XZ_CUR)))
- ELSE
- ! attenuation from ZAELOC(JK) and ZAELOC(JK-1)
- ZAETOT=ZAETOT*EXP(-2.*(ZAELOC(JK-1)*(.5*(ZZMZ(JK)+ZZMZ(JK-1))-TPFLYER%XZ_CUR) &
- +ZAELOC(JK)*.5*(ZZMZ(JK)-ZZMZ(JK-1))))
+ ! apply attenuation
+ ALLOCATE(ZZMZ(IKU))
+ ZZMZ(:)=FLYER_INTERPZ(ZZM(:,:,:))
+ ! nadir
+ ZAETOT=1.
+ DO JK=COUNT(TPFLYER%XZ_CUR >= ZZMZ(:)),1,-1
+ IF(JK.EQ.COUNT(TPFLYER%XZ_CUR >= ZZMZ(:))) THEN
+ IF(TPFLYER%XZ_CUR<=ZZMZ(JK)+.5*(ZZMZ(JK+1)-ZZMZ(JK))) THEN
+ ! only attenuation from ZAELOC(JK)
+ ZAETOT=ZAETOT*EXP(-2.*(ZAELOC(JK)*(TPFLYER%XZ_CUR-ZZMZ(JK))))
+ ELSE
+ ! attenuation from ZAELOC(JK) and ZAELOC(JK+1)
+ ZAETOT=ZAETOT*EXP(-2.*(ZAELOC(JK+1)*(TPFLYER%XZ_CUR-.5*(ZZMZ(JK+1)+ZZMZ(JK))) &
+ +ZAELOC(JK)*.5*(ZZMZ(JK+1)-ZZMZ(JK))))
+ END IF
+ ELSE
+ ! attenuation from ZAELOC(JK) and ZAELOC(JK+1)
+ ZAETOT=ZAETOT*EXP(-(ZAELOC(JK+1)+ZAELOC(JK))*(ZZMZ(JK+1)-ZZMZ(JK)))
+ END IF
+ TPFLYER%XCRARE_ATT(ISTORE,JK)=TPFLYER%XCRARE(ISTORE,JK)*ZAETOT
+ END DO
+ ! zenith
+ ZAETOT=1.
+ DO JK = MAX(COUNT(TPFLYER%XZ_CUR >= ZZMZ(:)),1)+1,IKU
+ IF ( JK .EQ. (MAX(COUNT(TPFLYER%XZ_CUR >= ZZMZ(:)),1)+1) ) THEN
+ IF(TPFLYER%XZ_CUR>=ZZMZ(JK)-.5*(ZZMZ(JK)-ZZMZ(JK-1))) THEN
+ ! only attenuation from ZAELOC(JK)
+ ZAETOT=ZAETOT*EXP(-2.*(ZAELOC(JK)*(ZZMZ(JK)-TPFLYER%XZ_CUR)))
+ ELSE
+ ! attenuation from ZAELOC(JK) and ZAELOC(JK-1)
+ ZAETOT=ZAETOT*EXP(-2.*(ZAELOC(JK-1)*(.5*(ZZMZ(JK)+ZZMZ(JK-1))-TPFLYER%XZ_CUR) &
+ +ZAELOC(JK)*.5*(ZZMZ(JK)-ZZMZ(JK-1))))
+ END IF
+ ELSE
+ ! attenuation from ZAELOC(JK) and ZAELOC(JK-1)
+ ZAETOT=ZAETOT*EXP(-(ZAELOC(JK-1)+ZAELOC(JK))*(ZZMZ(JK)-ZZMZ(JK-1)))
+ END IF
+ TPFLYER%XCRARE_ATT(ISTORE,JK)=TPFLYER%XCRARE(ISTORE,JK)*ZAETOT
+ END DO
+ TPFLYER%XZZ (ISTORE,:) = ZZMZ(:)
+ DEALLOCATE(ZZMZ,ZAELOC)
+ ! m^3 → mm^6/m^3 → dBZ
+ WHERE(TPFLYER%XCRARE(ISTORE,:)>0)
+ TPFLYER%XCRARE(ISTORE,:)=10.*LOG10(1.E18*TPFLYER%XCRARE(ISTORE,:))
+ ELSEWHERE
+ TPFLYER%XCRARE(ISTORE,:)=XUNDEF
+ END WHERE
+ WHERE(TPFLYER%XCRARE_ATT(ISTORE,:)>0)
+ TPFLYER%XCRARE_ATT(ISTORE,:)=10.*LOG10(1.E18*TPFLYER%XCRARE_ATT(ISTORE,:))
+ ELSEWHERE
+ TPFLYER%XCRARE_ATT(ISTORE,:)=XUNDEF
+ END WHERE
+ DEALLOCATE(ZX,ZW,ZRTMIN)
+ END IF ! end LOOP ICE3
+ ! vertical wind
+ TPFLYER%XWZ (ISTORE,:) = FLYER_INTERPZ(ZWM(:,:,:))
+ IF (SIZE(PTKE)>0) TPFLYER%XTKE (ISTORE) = FLYER_INTERP(PTKE)
+ IF (SIZE(PTS) >0) TPFLYER%XTSRAD(ISTORE) = FLYER_INTERP_2D(PTS)
+ IF (LDIAG_IN_RUN) TPFLYER%XTKE_DISS(ISTORE) = FLYER_INTERP(XCURRENT_TKE_DISS)
+ TPFLYER%XZS(ISTORE) = FLYER_INTERP_2D(PZ(:,:,1+JPVEXT))
+ TPFLYER%XTHW_FLUX(ISTORE) = FLYER_INTERP(ZTHW_FLUX)
+ TPFLYER%XRCW_FLUX(ISTORE) = FLYER_INTERP(ZRCW_FLUX)
+ DO JLOOP=1,SIZE(PSV,4)
+ TPFLYER%XSVW_FLUX(ISTORE,JLOOP) = FLYER_INTERP(ZSVW_FLUX(:,:,:,JLOOP))
+ END DO
END IF
- ELSE
- ! attenuation from ZAELOC(JK) and ZAELOC(JK-1)
- ZAETOT=ZAETOT*EXP(-(ZAELOC(JK-1)+ZAELOC(JK))*(ZZMZ(JK)-ZZMZ(JK-1)))
END IF
- TPFLYER%XCRARE_ATT(IN,JK)=TPFLYER%XCRARE(IN,JK)*ZAETOT
- END DO
- TPFLYER%XZZ (IN,:) = ZZMZ(:)
- DEALLOCATE(ZZMZ,ZAELOC)
- ! m^3 → mm^6/m^3 → dBZ
- WHERE(TPFLYER%XCRARE(IN,:)>0)
- TPFLYER%XCRARE(IN,:)=10.*LOG10(1.E18*TPFLYER%XCRARE(IN,:))
- ELSEWHERE
- TPFLYER%XCRARE(IN,:)=XUNDEF
- END WHERE
- WHERE(TPFLYER%XCRARE_ATT(IN,:)>0)
- TPFLYER%XCRARE_ATT(IN,:)=10.*LOG10(1.E18*TPFLYER%XCRARE_ATT(IN,:))
- ELSEWHERE
- TPFLYER%XCRARE_ATT(IN,:)=XUNDEF
- END WHERE
- DEALLOCATE(ZX,ZW,ZRTMIN)
- END IF ! end LOOP ICE3
- ! vertical wind
- TPFLYER%XWZ (IN,:) = FLYER_INTERPZ(ZWM(:,:,:))
- IF (SIZE(PTKE)>0) TPFLYER%XTKE (IN) = FLYER_INTERP(PTKE)
- IF (SIZE(PTS) >0) TPFLYER%XTSRAD(IN) = FLYER_INTERP_2D(PTS)
- IF (LDIAG_IN_RUN) TPFLYER%XTKE_DISS(IN) = FLYER_INTERP(XCURRENT_TKE_DISS)
- TPFLYER%XZS(IN) = FLYER_INTERP_2D(PZ(:,:,1+JPVEXT))
- TPFLYER%XTHW_FLUX(IN) = FLYER_INTERP(ZTHW_FLUX)
- TPFLYER%XRCW_FLUX(IN) = FLYER_INTERP(ZRCW_FLUX)
- DO JLOOP=1,SIZE(PSV,4)
- TPFLYER%XSVW_FLUX(IN,JLOOP) = FLYER_INTERP(ZSVW_FLUX(:,:,:,JLOOP))
- END DO
- END IF
-!
-!----------------------------------------------------------------------------
-!
-!* 9. BALLOON ADVECTION
-! -----------------
-!
- SELECT TYPE ( TPFLYER )
- CLASS IS ( TBALLOONDATA)
+ !
+ !----------------------------------------------------------------------------
+ !
+ !* 9. BALLOON ADVECTION
+ ! -----------------
+ !
+ ! SELECT TYPE ( TPFLYER )
+ ! CLASS IS ( TBALLOONDATA)
+ ZTSTEP = PTSTEP
+
IF ( TPFLYER%CTYPE == 'RADIOS' .OR. TPFLYER%CTYPE == 'ISODEN' .OR. TPFLYER%CTYPE == 'CVBALL' ) THEN
ZU_BAL = FLYER_INTERP_U(PU)
ZV_BAL = FLYER_INTERP_V(PV)
@@ -1292,21 +1812,113 @@ IF ( TPFLYER%LFLY ) THEN
ZMAP = 1.
end if
!
- TPFLYER%XX_CUR = TPFLYER%XX_CUR + ZU_BAL * PTSTEP * ZMAP
- TPFLYER%XY_CUR = TPFLYER%XY_CUR + ZV_BAL * PTSTEP * ZMAP
+ ZX_OLD = TPFLYER%XX_CUR
+ ZY_OLD = TPFLYER%XY_CUR
+
+ TPFLYER%XX_CUR = TPFLYER%XX_CUR + ZU_BAL * ZTSTEP * ZMAP
+ TPFLYER%XY_CUR = TPFLYER%XY_CUR + ZV_BAL * ZTSTEP * ZMAP
END IF
!
+ !Compute rank and model for next position
+ !This is done here because we need to check if there is a change of model (for 'MOB' balloons)
+ !because position has to be adapted to the timestep of a coarser model (if necessary)
+ IMODEL_OLD = TPFLYER%NMODEL
+ ! Get rank of the process where the balloon is and the model number
+ IF ( TPFLYER%CMODEL == 'FIX' ) THEN
+ IMODEL = TPFLYER%NMODEL
+ ELSE
+ IMODEL = 0
+ END IF
+ CALL FIND_PROCESS_AND_MODEL_FROM_XY_POS( TPFLYER%XX_CUR, TPFLYER%XY_CUR, IRANK, IMODEL )
+ IF ( IRANK < 1 ) THEN
+ TPFLYER%NMODEL = 1 !Set to 1 because it is always a valid model (to prevent crash if NDAD(TPFLYER%NMODEL) )
+ TPFLYER%LCRASH = .TRUE.
+ TPFLYER%NCRASH = NCRASH_OUT_HORIZ
+ TPFLYER%LFLY = .FALSE.
+ WRITE( CMNHMSG(1), "( 'Balloon ', A, ' crashed the ', I2, '/', I2, '/', I4, ' at ', F18.12, &
+ 's (out of the horizontal boundaries)' )" ) &
+ TRIM( TPFLYER%CTITLE ), TDTCUR%NDAY, TDTCUR%NMONTH, TDTCUR%NYEAR, TDTCUR%XTIME
+ CALL PRINT_MSG( NVERB_WARNING, 'GEN', 'AIRCRAFT_BALLOON_EVOL', OLOCAL = .TRUE. )
+ ELSE
+ TPFLYER%NMODEL = IMODEL
+ TPFLYER%LCRASH = .FALSE.
+ TPFLYER%NCRASH = NCRASH_NO
+ TPFLYER%NRANK_CUR = IRANK
+ END IF
+ IF ( TPFLYER%NMODEL /= IMODEL_OLD .AND. .NOT. TPFLYER%LCRASH ) THEN
+ IF ( NDAD(TPFLYER%NMODEL ) == IMODEL_OLD ) THEN
+ !Nothing special to do
+ ELSE IF ( TPFLYER%NMODEL == NDAD(IMODEL_OLD) ) THEN
+ !Recompute position to be compatible with parent timestep
+
+ !Determine step compatible with parent model at next parent timestep
+ ZDELTATIME = TDTCUR - TDTSEG
+ ZDIVTMP = ZDELTATIME / ( PTSTEP * NDTRATIO(IMODEL_OLD) )
+ IF ( ABS( ZDIVTMP - NINT( ZDIVTMP ) ) < 1E-6 * PTSTEP * NDTRATIO(IMODEL_OLD) ) THEN
+ ZTSTEP = ZTSTEP * NDTRATIO(IMODEL_OLD)
+ ELSE
+ ZTSTEP = ZTSTEP * ( 1 + NDTRATIO(IMODEL_OLD) )
+ ! Detect if we need to skip a store (if time of next position is after time of next store)
+ ! This can happen when a ballon goes to its parent model
+ IF ( TDTCUR + ZTSTEP > &
+ TPFLYER%TFLYER_TIME%TPDATES(TPFLYER%TFLYER_TIME%N_CUR) + TPFLYER%TFLYER_TIME%XTSTEP + 1e-6 ) THEN
+ TPFLYER%TFLYER_TIME%N_CUR = TPFLYER%TFLYER_TIME%N_CUR + 1
+ ISTORE = TPFLYER%TFLYER_TIME%N_CUR
+ !Remark: by construction here, ISTORE is always > 1 => no risk with ISTORE-1 value
+ TPFLYER%TFLYER_TIME%TPDATES(ISTORE) = TPFLYER%TFLYER_TIME%TPDATES(ISTORE-1) + TPFLYER%TFLYER_TIME%XTSTEP
+
+ !Force a dummy store (nothing is computed, therefore default/initial values will be stored)
+ TPFLYER%LSTORE = .TRUE.
+
+ WRITE( CMNHMSG(1), "( 'Balloon ', A, ': store skipped at ', I2, '/', I2, '/', I4, ' at ', F18.12, 's' )" ) &
+ TRIM( TPFLYER%CTITLE ), &
+ TPFLYER%TFLYER_TIME%TPDATES(ISTORE)%NDAY, TPFLYER%TFLYER_TIME%TPDATES(ISTORE)%NMONTH, &
+ TPFLYER%TFLYER_TIME%TPDATES(ISTORE)%NYEAR, TPFLYER%TFLYER_TIME%TPDATES(ISTORE)%XTIME
+ CMNHMSG(2) = 'due to change of model (child to its parent)'
+ CALL PRINT_MSG( NVERB_WARNING, 'GEN', 'AIRCRAFT_BALLOON_EVOL', OLOCAL = .TRUE. )
+ END IF
+ END IF
+
+ TPFLYER%XX_CUR = TPFLYER%XX_CUR + ZU_BAL * ZTSTEP * ZMAP
+ TPFLYER%XY_CUR = TPFLYER%XY_CUR + ZV_BAL * ZTSTEP * ZMAP
+
+ ! Model number is now imposed
+ IMODEL = TPFLYER%NMODEL
+ CALL FIND_PROCESS_AND_MODEL_FROM_XY_POS( TPFLYER%XX_CUR, TPFLYER%XY_CUR, IRANK, IMODEL )
+ IF ( IRANK < 1 ) THEN
+ TPFLYER%NMODEL = 1 !Set to 1 because it is always a valid model (to prevent crash if NDAD(TPFLYER%NMODEL) )
+ TPFLYER%LCRASH = .TRUE.
+ TPFLYER%NCRASH = NCRASH_OUT_HORIZ
+ TPFLYER%LFLY = .FALSE.
+ WRITE( CMNHMSG(1), "( 'Balloon ', A, ' crashed the ', I2, '/', I2, '/', I4, ' at ', F18.12, &
+ 's (out of the horizontal boundaries)' )" ) &
+ TRIM( TPFLYER%CTITLE ), TDTCUR%NDAY, TDTCUR%NMONTH, TDTCUR%NYEAR, TDTCUR%XTIME
+ CALL PRINT_MSG( NVERB_WARNING, 'GEN', 'AIRCRAFT_BALLOON_EVOL', OLOCAL = .TRUE. )
+ ELSE
+ !TPFLYER%NMODEL = IMODEL !Do not change model because we are in transition to parent
+ TPFLYER%LCRASH = .FALSE.
+ TPFLYER%NCRASH = NCRASH_NO
+ TPFLYER%NRANK_CUR = IRANK
+ END IF
+ ELSE
+ !Special case not-managed (different dads, change of several models in 1 step (going to grand parent)...)
+ CMNHMSG(1) = 'unmanaged change of model for ballon ' // TPFLYER%CTITLE
+ CMNHMSG(2) = 'its trajectory might be wrong'
+ CALL PRINT_MSG( NVERB_WARNING, 'GEN', 'AIRCRAFT_BALLOON_EVOL', OLOCAL = .TRUE. )
+ END IF
+ END IF
+
IF ( TPFLYER%CTYPE == 'RADIOS' ) THEN
ZW_BAL = FLYER_INTERP(ZWM)
- TPFLYER%XZ_CUR = TPFLYER%XZ_CUR + ( ZW_BAL + TPFLYER%XWASCENT ) * PTSTEP
+ TPFLYER%XZ_CUR = TPFLYER%XZ_CUR + ( ZW_BAL + TPFLYER%XWASCENT ) * ZTSTEP
END IF
!
IF ( TPFLYER%CTYPE == 'CVBALL' ) THEN
ZW_BAL = FLYER_INTERP(ZWM)
ZRO_BAL = FLYER_INTERP(ZRHO)
! calculation with a time step of 1 second or less
- IF (INT(PTSTEP) .GT. 1 ) THEN
- DO JK=1,INT(PTSTEP)
+ IF (INT(ZTSTEP) .GT. 1 ) THEN
+ DO JK=1,INT(ZTSTEP)
TPFLYER%XWASCENT = TPFLYER%XWASCENT &
- ( 1. / (1. + TPFLYER%XINDDRAG ) ) * 1. * &
( XG * ( ( TPFLYER%XMASS / TPFLYER%XVOLUME ) - ZRO_BAL ) / ( TPFLYER%XMASS / TPFLYER%XVOLUME ) &
@@ -1316,161 +1928,114 @@ IF ( TPFLYER%LFLY ) THEN
TPFLYER%XZ_CUR = TPFLYER%XZ_CUR + ( ZW_BAL + TPFLYER%XWASCENT ) * 1.
END DO
END IF
- IF (PTSTEP .GT. INT(PTSTEP)) THEN
- TPFLYER%XWASCENT = TPFLYER%XWASCENT &
- - ( 1. / (1. + TPFLYER%XINDDRAG ) ) * (PTSTEP-INT(PTSTEP)) * &
- ( XG * ( ( TPFLYER%XMASS / TPFLYER%XVOLUME ) - ZRO_BAL ) / ( TPFLYER%XMASS / TPFLYER%XVOLUME ) &
- + TPFLYER%XWASCENT * ABS ( TPFLYER%XWASCENT ) * &
- TPFLYER%XDIAMETER * TPFLYER%XAERODRAG / ( 2. * TPFLYER%XVOLUME ) &
- )
- TPFLYER%XZ_CUR = TPFLYER%XZ_CUR + ( ZW_BAL + TPFLYER%XWASCENT ) * (PTSTEP-INT(PTSTEP))
+ IF (ZTSTEP .GT. INT(ZTSTEP)) THEN
+ TPFLYER%XWASCENT = TPFLYER%XWASCENT &
+ - ( 1. / (1. + TPFLYER%XINDDRAG ) ) * (ZTSTEP-INT(ZTSTEP)) * &
+ ( XG * ( ( TPFLYER%XMASS / TPFLYER%XVOLUME ) - ZRO_BAL ) / ( TPFLYER%XMASS / TPFLYER%XVOLUME ) &
+ + TPFLYER%XWASCENT * ABS ( TPFLYER%XWASCENT ) * &
+ TPFLYER%XDIAMETER * TPFLYER%XAERODRAG / ( 2. * TPFLYER%XVOLUME ) &
+ )
+ TPFLYER%XZ_CUR = TPFLYER%XZ_CUR + ( ZW_BAL + TPFLYER%XWASCENT ) * (ZTSTEP-INT(ZTSTEP))
END IF
END IF
- END SELECT
-!
-!----------------------------------------------------------------------------
- END IF
-!----------------------------------------------------------------------------
-!
-!* 10. AIRCRAFT MOVE (computations done on all processors, to limit exchanges)
-! -------------
-!
- SELECT TYPE ( TPFLYER )
- CLASS IS ( TAIRCRAFTDATA )
-!
-!
-!* 10.1 Determination of flight segment
-! -------------------------------
-!
- IL = TPFLYER%NSEGCURN
- !
- TPFLYER%XSEGCURT = TPFLYER%XSEGCURT + PTSTEP
- !
- DO WHILE (TPFLYER%XSEGCURT>TPFLYER%XSEGTIME(IL))
- TPFLYER%NSEGCURN = TPFLYER%NSEGCURN + 1
- IL = TPFLYER%NSEGCURN
- TPFLYER%XSEGCURT = TPFLYER%XSEGCURT - TPFLYER%XSEGTIME(IL-1)
- IF (IL>TPFLYER%NSEG) EXIT
- END DO
-! DO WHILE (TPFLYER%XSEGCURT>TPFLYER%XSEGTIME(IL) .AND. IL <= TPFLYER%NSEG)
-! TPFLYER%NSEGCURN = TPFLYER%NSEGCURN + 1
-! IL = TPFLYER%NSEGCURN
-! TPFLYER%XSEGCURT = TPFLYER%XSEGCURT - TPFLYER%XSEGTIME(IL-1)
-! END DO
- !
- !* end of flight
- !
- IF (IL > TPFLYER%NSEG) TPFLYER%LFLY = .FALSE.
-!
-!
-!* 10.2 Determination of new position
-! -----------------------------
-!
- IF (TPFLYER%LFLY) THEN
- ZSEG_FRAC = TPFLYER%XSEGCURT / TPFLYER%XSEGTIME(IL)
- !
- TPFLYER%XX_CUR = (1.-ZSEG_FRAC) * TPFLYER%XSEGX(IL ) &
- + ZSEG_FRAC * TPFLYER%XSEGX(IL+1)
- TPFLYER%XY_CUR = (1.-ZSEG_FRAC) * TPFLYER%XSEGY(IL ) &
- + ZSEG_FRAC * TPFLYER%XSEGY(IL+1)
- IF (TPFLYER%LALTDEF) THEN
- TPFLYER%XP_CUR = (1.-ZSEG_FRAC) * TPFLYER%XSEGP(IL ) &
- + ZSEG_FRAC * TPFLYER%XSEGP(IL+1)
- ELSE
- TPFLYER%XZ_CUR = (1.-ZSEG_FRAC) * TPFLYER%XSEGZ(IL ) &
- + ZSEG_FRAC * TPFLYER%XSEGZ(IL+1)
- END IF
- END IF
- END SELECT
- !
- END IF
-!
-END IF
-!
-!----------------------------------------------------------------------------
-!
-!* 11. EXCHANGE OF INFORMATION BETWEEN PROCESSORS
-! ------------------------------------------
-!
-!* 11.1 current position
-! ----------------
-!
-CALL DISTRIBUTE_FLYER_L(TPFLYER%LFLY)
-CALL DISTRIBUTE_FLYER_L(TPFLYER%LCRASH)
-CALL DISTRIBUTE_FLYER(TPFLYER%XX_CUR)
-CALL DISTRIBUTE_FLYER(TPFLYER%XY_CUR)
-SELECT TYPE ( TPFLYER )
- CLASS IS ( TBALLOONDATA )
- IF ( TPFLYER%CTYPE == 'CVBALL' ) THEN
- CALL DISTRIBUTE_FLYER(TPFLYER%XZ_CUR)
- CALL DISTRIBUTE_FLYER(TPFLYER%XWASCENT)
- ELSE IF ( TPFLYER%CTYPE == 'RADIOS' ) THEN
- CALL DISTRIBUTE_FLYER(TPFLYER%XZ_CUR)
- ELSE IF ( TPFLYER%CTYPE == 'ISODEN' ) THEN
- CALL DISTRIBUTE_FLYER(TPFLYER%XRHO)
- END IF
+ TPFLYER%TPOS_CUR = TDTCUR + ZTSTEP
+ END IF CRASH_VERT !end of no vertical crash branch
+ ELSE ISOWNERBAL
+ !The balloon is not present on this MPI process
+ ZTHIS_PROC = 0.
+ END IF ISOWNERBAL
+ !----------------------------------------------------------------------------
+ !
+ !* 11. EXCHANGE OF INFORMATION BETWEEN PROCESSES
+ ! -----------------------------------------
+ !
+ !* 11.1 current position
+ ! ----------------
+ !
+ IF ( TPFLYER%CMODEL == 'MOB' ) THEN
+ CALL DISTRIBUTE_FLYER_N(TPFLYER%NMODEL)
+ END IF
+ CALL DISTRIBUTE_FLYER_N(TPFLYER%NRANK_CUR)
+ CALL DISTRIBUTE_FLYER_L(TPFLYER%LFLY)
+ CALL DISTRIBUTE_FLYER_L(TPFLYER%LCRASH)
+ CALL DISTRIBUTE_FLYER_L(TPFLYER%LSTORE)
+ CALL DISTRIBUTE_FLYER(TPFLYER%XX_CUR)
+ CALL DISTRIBUTE_FLYER(TPFLYER%XY_CUR)
+ CALL DISTRIBUTE_FLYER_N(TPFLYER%TPOS_CUR%NYEAR)
+ CALL DISTRIBUTE_FLYER_N(TPFLYER%TPOS_CUR%NMONTH)
+ CALL DISTRIBUTE_FLYER_N(TPFLYER%TPOS_CUR%NDAY)
+ CALL DISTRIBUTE_FLYER (TPFLYER%TPOS_CUR%XTIME)
+
+ CALL DISTRIBUTE_FLYER_N(TPFLYER%TFLYER_TIME%N_CUR)
+
+ ! SELECT TYPE ( TPFLYER )
+ ! CLASS IS ( TBALLOONDATA )
+ IF ( TPFLYER%CTYPE == 'CVBALL' ) THEN
+ CALL DISTRIBUTE_FLYER(TPFLYER%XZ_CUR)
+ CALL DISTRIBUTE_FLYER(TPFLYER%XWASCENT)
+ ELSE IF ( TPFLYER%CTYPE == 'RADIOS' ) THEN
+ CALL DISTRIBUTE_FLYER(TPFLYER%XZ_CUR)
+ ELSE IF ( TPFLYER%CTYPE == 'ISODEN' ) THEN
+ CALL DISTRIBUTE_FLYER(TPFLYER%XRHO)
+ END IF
+
+ IF ( TPFLYER%LSTORE ) THEN
+ CALL DISTRIBUTE_FLYER_N(TPFLYER%TFLYER_TIME%TPDATES(ISTORE)%NYEAR)
+ CALL DISTRIBUTE_FLYER_N(TPFLYER%TFLYER_TIME%TPDATES(ISTORE)%NMONTH)
+ CALL DISTRIBUTE_FLYER_N(TPFLYER%TFLYER_TIME%TPDATES(ISTORE)%NDAY)
+ CALL DISTRIBUTE_FLYER (TPFLYER%TFLYER_TIME%TPDATES(ISTORE)%XTIME)
+
+ CALL DISTRIBUTE_FLYER_N(TPFLYER%NMODELHIST(ISTORE))
+ CALL DISTRIBUTE_FLYER(TPFLYER%XX (ISTORE))
+ CALL DISTRIBUTE_FLYER(TPFLYER%XY (ISTORE))
+ CALL DISTRIBUTE_FLYER(TPFLYER%XZ (ISTORE))
+ CALL DISTRIBUTE_FLYER(TPFLYER%XLON(ISTORE))
+ CALL DISTRIBUTE_FLYER(TPFLYER%XLAT(ISTORE))
+ CALL DISTRIBUTE_FLYER(TPFLYER%XZON(ISTORE))
+ CALL DISTRIBUTE_FLYER(TPFLYER%XMER(ISTORE))
+ CALL DISTRIBUTE_FLYER(TPFLYER%XW (ISTORE))
+ CALL DISTRIBUTE_FLYER(TPFLYER%XP (ISTORE))
+ CALL DISTRIBUTE_FLYER(TPFLYER%XTH (ISTORE))
+ DO JLOOP=1,SIZE(PR,4)
+ CALL DISTRIBUTE_FLYER(TPFLYER%XR (ISTORE,JLOOP))
+ END DO
+ DO JLOOP=1,SIZE(PSV,4)
+ CALL DISTRIBUTE_FLYER(TPFLYER%XSV (ISTORE,JLOOP))
+ END DO
+ DO JLOOP=1,IKU
+ CALL DISTRIBUTE_FLYER(TPFLYER%XRTZ (ISTORE,JLOOP))
+ DO JLOOP2=1,SIZE(PR,4)
+ CALL DISTRIBUTE_FLYER(TPFLYER%XRZ (ISTORE,JLOOP,JLOOP2))
+ ENDDO
+ CALL DISTRIBUTE_FLYER(TPFLYER%XFFZ (ISTORE,JLOOP))
+ CALL DISTRIBUTE_FLYER(TPFLYER%XCIZ (ISTORE,JLOOP))
+ IF (CCLOUD== 'LIMA' ) THEN
+ CALL DISTRIBUTE_FLYER(TPFLYER%XCRZ (ISTORE,JLOOP))
+ CALL DISTRIBUTE_FLYER(TPFLYER%XCCZ (ISTORE,JLOOP))
+ ENDIF
+ CALL DISTRIBUTE_FLYER(TPFLYER%XIWCZ (ISTORE,JLOOP))
+ CALL DISTRIBUTE_FLYER(TPFLYER%XLWCZ (ISTORE,JLOOP))
+ CALL DISTRIBUTE_FLYER(TPFLYER%XCRARE (ISTORE,JLOOP))
+ CALL DISTRIBUTE_FLYER(TPFLYER%XCRARE_ATT (ISTORE,JLOOP))
+ CALL DISTRIBUTE_FLYER(TPFLYER%XWZ (ISTORE,JLOOP))
+ CALL DISTRIBUTE_FLYER(TPFLYER%XZZ (ISTORE,JLOOP))
+ END DO
+ IF (SIZE(PTKE)>0) CALL DISTRIBUTE_FLYER(TPFLYER%XTKE (ISTORE))
+ IF (SIZE(PTS) >0) CALL DISTRIBUTE_FLYER(TPFLYER%XTSRAD(ISTORE))
+ IF (LDIAG_IN_RUN) CALL DISTRIBUTE_FLYER(TPFLYER%XTKE_DISS(ISTORE))
+ CALL DISTRIBUTE_FLYER(TPFLYER%XZS (ISTORE))
+ CALL DISTRIBUTE_FLYER(TPFLYER%XTHW_FLUX(ISTORE))
+ CALL DISTRIBUTE_FLYER(TPFLYER%XRCW_FLUX(ISTORE))
+ DO JLOOP=1,SIZE(PSV,4)
+ CALL DISTRIBUTE_FLYER(TPFLYER%XSVW_FLUX(ISTORE,JLOOP))
+ END DO
+ END IF
+ END IF INFLIGHTONMODEL
- CLASS IS ( TAIRCRAFTDATA )
- IF (TPFLYER%LALTDEF) THEN
- CALL DISTRIBUTE_FLYER(TPFLYER%XP_CUR)
- ELSE
- CALL DISTRIBUTE_FLYER(TPFLYER%XZ_CUR)
- ENDIF
END SELECT
-!
-!* 11.2 data stored
-! -----------
-!
-IF ( GSTORE ) THEN
- CALL DISTRIBUTE_FLYER(TPFLYER%XX (IN))
- CALL DISTRIBUTE_FLYER(TPFLYER%XY (IN))
- CALL DISTRIBUTE_FLYER(TPFLYER%XZ (IN))
- CALL DISTRIBUTE_FLYER(TPFLYER%XLON(IN))
- CALL DISTRIBUTE_FLYER(TPFLYER%XLAT(IN))
- CALL DISTRIBUTE_FLYER(TPFLYER%XZON(IN))
- CALL DISTRIBUTE_FLYER(TPFLYER%XMER(IN))
- CALL DISTRIBUTE_FLYER(TPFLYER%XW (IN))
- CALL DISTRIBUTE_FLYER(TPFLYER%XP (IN))
- CALL DISTRIBUTE_FLYER(TPFLYER%XTH (IN))
- DO JLOOP=1,SIZE(PR,4)
- CALL DISTRIBUTE_FLYER(TPFLYER%XR (IN,JLOOP))
- END DO
- DO JLOOP=1,SIZE(PSV,4)
- CALL DISTRIBUTE_FLYER(TPFLYER%XSV (IN,JLOOP))
- END DO
- DO JLOOP=1,IKU
- CALL DISTRIBUTE_FLYER(TPFLYER%XRTZ (IN,JLOOP))
- DO JLOOP2=1,SIZE(PR,4)
- CALL DISTRIBUTE_FLYER(TPFLYER%XRZ (IN,JLOOP,JLOOP2))
- ENDDO
- CALL DISTRIBUTE_FLYER(TPFLYER%XFFZ (IN,JLOOP))
- CALL DISTRIBUTE_FLYER(TPFLYER%XCIZ (IN,JLOOP))
- IF (CCLOUD== 'LIMA' ) THEN
- CALL DISTRIBUTE_FLYER(TPFLYER%XCRZ (IN,JLOOP))
- CALL DISTRIBUTE_FLYER(TPFLYER%XCCZ (IN,JLOOP))
- ENDIF
- CALL DISTRIBUTE_FLYER(TPFLYER%XIWCZ (IN,JLOOP))
- CALL DISTRIBUTE_FLYER(TPFLYER%XLWCZ (IN,JLOOP))
- CALL DISTRIBUTE_FLYER(TPFLYER%XCRARE (IN,JLOOP))
- CALL DISTRIBUTE_FLYER(TPFLYER%XCRARE_ATT (IN,JLOOP))
- CALL DISTRIBUTE_FLYER(TPFLYER%XWZ (IN,JLOOP))
- CALL DISTRIBUTE_FLYER(TPFLYER%XZZ (IN,JLOOP))
- END DO
- IF (SIZE(PTKE)>0) CALL DISTRIBUTE_FLYER(TPFLYER%XTKE (IN))
- IF (SIZE(PTS) >0) CALL DISTRIBUTE_FLYER(TPFLYER%XTSRAD(IN))
- IF (LDIAG_IN_RUN) CALL DISTRIBUTE_FLYER(TPFLYER%XTKE_DISS(IN))
- CALL DISTRIBUTE_FLYER(TPFLYER%XZS (IN))
- CALL DISTRIBUTE_FLYER(TPFLYER%XTHW_FLUX(IN))
- CALL DISTRIBUTE_FLYER(TPFLYER%XRCW_FLUX(IN))
- DO JLOOP=1,SIZE(PSV,4)
- CALL DISTRIBUTE_FLYER(TPFLYER%XSVW_FLUX(IN,JLOOP))
- END DO
-END IF
-!
-!----------------------------------------------------------------------------
-!----------------------------------------------------------------------------
-!
+
+
CONTAINS
!
!----------------------------------------------------------------------------
@@ -1486,8 +2051,8 @@ IF (SIZE(PA,1)==2) THEN
JI=1
JJ=1
ELSE
- JI=II
- JJ=IJ
+ JI=II_M
+ JJ=IJ_M
END IF
!
PB = (1.- ZYCOEF) * (1.-ZXCOEF) * ( (1.-ZZCOEF00) * PA(JI ,JJ ,IK00) + ZZCOEF00 * PA(JI ,JJ ,IK00+1)) &
@@ -1508,8 +2073,8 @@ IF (SIZE(PA,1)==2) THEN
JI=1
JJ=1
ELSE
- JI=II
- JJ=IJ
+ JI=II_M
+ JJ=IJ_M
END IF
!
!
@@ -1519,10 +2084,10 @@ DO JK=1,SIZE(PA,3)
PB(JK) = (1.-ZYCOEF) * (1.-ZXCOEF) * PA(JI,JJ,JK) + &
(1.-ZYCOEF) * (ZXCOEF) * PA(JI+1,JJ,JK) + &
(ZYCOEF) * (1.-ZXCOEF) * PA(JI,JJ+1,JK) + &
- (ZYCOEF) * (ZXCOEF) * PA(JI+1,JJ+1,JK)
+ (ZYCOEF) * (ZXCOEF) * PA(JI+1,JJ+1,JK)
ELSE
- PB(JK) = XUNDEF
- END IF
+ PB(JK) = XUNDEF
+ END IF
END DO
!
END FUNCTION FLYER_INTERPZ
@@ -1538,8 +2103,8 @@ IF (SIZE(PA,1)==2) THEN
JI=1
JJ=1
ELSE
- JI=IU
- JJ=IJ
+ JI=II_U
+ JJ=IJ_M
END IF
!
PB = (1.- ZYCOEF) * (1.-ZUCOEF) * ( (1.-ZUCOEF00) * PA(JI ,JJ ,IU00) + ZUCOEF00 * PA(JI ,JJ ,IU00+1)) &
@@ -1561,8 +2126,8 @@ IF (SIZE(PA,1)==2) THEN
JI=1
JJ=1
ELSE
- JI=II
- JJ=IV
+ JI=II_M
+ JJ=IJ_V
END IF
!
PB = (1.- ZVCOEF) * (1.-ZXCOEF) * ( (1.-ZVCOEF00) * PA(JI ,JJ ,IV00) + ZVCOEF00 * PA(JI ,JJ ,IV00+1)) &
@@ -1584,8 +2149,8 @@ IF (SIZE(PA,1)==2) THEN
JI=1
JJ=1
ELSE
- JI=II
- JJ=IJ
+ JI=II_M
+ JJ=IJ_M
END IF
!
PB = (1.- ZYCOEF) * (1.-ZXCOEF) * PA(JI ,JJ ) &
@@ -1639,88 +2204,9 @@ END IF
!
END SUBROUTINE DISTRIBUTE_FLYER_L
!----------------------------------------------------------------------------
-!----------------------------------------------------------------------------
-SUBROUTINE FLYER_CHANGE_MODEL(IMI)
-!
-INTEGER, INTENT(IN) :: IMI ! model index
-!
-INTEGER :: IMK ! kid model index
-INTEGER :: IMODEL ! TPFLYER model index at the beginning of the subroutine
-INTEGER :: IU ! U flux point balloon position (x index)
-INTEGER :: IV ! V flux point balloon position (y index)
-INTEGER :: IU_ABS ! U flux point balloon position (in the model)
-INTEGER :: IV_ABS ! V flux point balloon position (in the model)
-INTEGER :: IXOR ! Origin's coordinates of the extended 2 way subdomain
-INTEGER :: IYOR ! Origin's coordinates of the extended 2 way subdomain
-INTEGER :: IIB ! current processor domain sizes
-INTEGER :: IJB
-INTEGER :: IIE
-INTEGER :: IJE
-!
-!
-IMODEL=TPFLYER%NMODEL
-CALL GET_INDICE_ll (IIB,IJB,IIE,IJE)
-IU=COUNT( PXHAT (:)<=TPFLYER%XX_CUR )
-IV=COUNT( PYHAT (:)<=TPFLYER%XY_CUR )
-ZTHIS_PROC=0.
-IF (IU>=IIB .AND. IU<=IIE .AND. IV>=IJB .AND. IV<=IJE) ZTHIS_PROC=1.
-IF (ZTHIS_PROC .EQ. 1) THEN
- CALL GET_OR_LL('B',IXOR,IYOR)
- IU_ABS=IU + IXOR - 1
- IV_ABS=IV + IYOR - 1
- !
- IF (TPFLYER%NMODEL == IMI) THEN
- !
- ! go to the kid model if the flyer location is inside
- ! ------------------
- !
- DO IMK=IMI+1,NMODEL
- IF (NDAD(IMK) == IMI .AND. &
- IU_ABS>=NXOR_ALL(IMK) .AND. IU_ABS<=NXEND_ALL(IMK) .AND. &
- IV_ABS>=NYOR_ALL(IMK) .AND. IV_ABS<=NYEND_ALL(IMK) ) THEN
- TPFLYER%NMODEL = IMK
- !
- END IF
- END DO
- !
- ELSE
- !
- ! come from the kid model if the flyer location is outside
- ! ------------------
- !
- IMK = TPFLYER%NMODEL
- IF (IU_ABS<NXOR_ALL(IMK) .OR. IU_ABS>NXEND_ALL(IMK) .OR. &
- IV_ABS<NYOR_ALL(IMK) .OR. IV_ABS>NYEND_ALL(IMK) ) THEN
- TPFLYER%NMODEL = IMI
- !
- END IF
- END IF
-END IF
-!
-! send the information to all the processors
-! ----------------------------------------
-!
-CALL DISTRIBUTE_FLYER_N(TPFLYER%NMODEL)
-ZTHIS_PROC=0.
-!
-! print if the model changes
-!---------------------------------
-IF (TPFLYER%NMODEL /= IMODEL) THEN
- IF (NDAD(IMODEL) == TPFLYER%NMODEL) THEN
- WRITE(ILUOUT,*) '-------------------------------------------------------------------'
- WRITE(ILUOUT,*) TPFLYER%CTITLE,' comes from model ',IMODEL,' in model ', &
- TPFLYER%NMODEL,' at ',NINT(TDTCUR%xtime),' sec.'
- WRITE(ILUOUT,*) '-------------------------------------------------------------------'
- ELSE
- WRITE(ILUOUT,*) '-------------------------------------------------------------------'
- WRITE(ILUOUT,*) TPFLYER%CTITLE,' goes from model ',IMODEL,' to model ', &
- TPFLYER%NMODEL,' at ',NINT(TDTCUR%xtime),' sec.'
- WRITE(ILUOUT,*) '-------------------------------------------------------------------'
- ENDIF
-ENDIF
-!
-!
-END SUBROUTINE FLYER_CHANGE_MODEL
-!----------------------------------------------------------------------------
-!----------------------------------------------------------------------------
+
+
END SUBROUTINE AIRCRAFT_BALLOON_EVOL
+!----------------------------------------------------------------------------
+
+END MODULE MODE_AIRCRAFT_BALLOON_EVOL
diff --git a/src/MNH/diag.f90 b/src/MNH/diag.f90
index 27362f835..59e39266e 100644
--- a/src/MNH/diag.f90
+++ b/src/MNH/diag.f90
@@ -91,7 +91,7 @@
! P. Wautelet 07/02/2019: force TYPE to a known value for IO_File_add2list
! P. Wautelet 11/02/2019: added missing use of MODI_CH_MONITOR_n
! P. Wautelet 28/03/2019: use MNHTIME for time measurement variables
-! P. Wautelet 26/07/2019: bug correction: deallocate of zsea and ztown done too early
+! P. Wautelet 26/07/2019: bug correction: deallocate of zsea done too early
! P. Wautelet 13/09/2019: budget: simplify and modernize date/time management
! P. Wautelet 06/07/2021: use FINALIZE_MNH
!-------------------------------------------------------------------------------
@@ -205,7 +205,7 @@ LOGICAL:: GCLOUD_ONLY ! conditionnal radiation computations for
! the only cloudy columns
!
INTEGER :: IIU, IJU, IKU
-REAL, DIMENSION(:,:),ALLOCATABLE :: ZSEA,ZTOWN
+REAL, DIMENSION(:,:),ALLOCATABLE :: ZSEA
REAL, DIMENSION(:,:,:,:),ALLOCATABLE :: ZWETDEPAER
!
TYPE(TFILEDATA),POINTER :: TZDIACFILE => NULL()
@@ -547,21 +547,17 @@ IF ( LAIRCRAFT_BALLOON ) THEN
TDTCUR = TXDTBAL !TDTCUR is used in AIRCRAFT_BALLOON
!
ALLOCATE (ZSEA(SIZE(XRHODJ,1),SIZE(XRHODJ,2)))
- ALLOCATE (ZTOWN(SIZE(XRHODJ,1),SIZE(XRHODJ,2)))
ZSEA(:,:) = 0.
- ZTOWN(:,:)= 0.
- CALL MNHGET_SURF_PARAM_n (PSEA=ZSEA(:,:),PTOWN=ZTOWN(:,:))
- DO ISTEPBAL=1,NTIME_AIRCRAFT_BALLOON,INT(XSTEP_AIRCRAFT_BALLOON)
- CALL AIRCRAFT_BALLOON(XSTEP_AIRCRAFT_BALLOON, &
- XXHAT, XYHAT, XZZ, XMAP, XLONORI, XLATORI, &
- XUT, XVT, XWT, XPABST, XTHT, XRT, XSVT, &
- XTKET, XTSRAD, XRHODREF,XCIT,ZSEA)
+ CALL MNHGET_SURF_PARAM_n (PSEA=ZSEA(:,:))
+ DO ISTEPBAL = 1, NTIME_AIRCRAFT_BALLOON, INT(XSTEP_AIRCRAFT_BALLOON)
+ CALL AIRCRAFT_BALLOON( XSTEP_AIRCRAFT_BALLOON, XZZ, XMAP, XLONORI, XLATORI, XUT, XVT, XWT, &
+ XPABST, XTHT, XRT, XSVT, XTKET, XTSRAD, XRHODREF, XCIT, ZSEA )
TXDTBAL%xtime = TXDTBAL%xtime + XSTEP_AIRCRAFT_BALLOON
CALL DATETIME_CORRECTDATE(TXDTBAL)
TDTCUR = TXDTBAL !TDTCUR is used in AIRCRAFT_BALLOON
- ENDDO
- DEALLOCATE (ZSEA,ZTOWN)
+ END DO
+ DEALLOCATE (ZSEA)
!
TDTCUR = TPDTCUR_SAVE
!
diff --git a/src/MNH/ini_aircraft.f90 b/src/MNH/ini_aircraft.f90
index f6355579a..d0f456f1a 100644
--- a/src/MNH/ini_aircraft.f90
+++ b/src/MNH/ini_aircraft.f90
@@ -118,6 +118,8 @@ ALLOCATE( TAIRCRAFTS(NAIRCRAFTS) )
!Treat aircraft data read in namelist
DO JI = 1, NAIRCRAFTS
+ TAIRCRAFTS(JI)%NID = JI
+
IF ( CTITLE(JI) == '' ) THEN
WRITE( CTITLE(JI), FMT = '( A, I3.3) ') TRIM( CTYPE(JI) ), JI
@@ -139,6 +141,7 @@ DO JI = 1, NAIRCRAFTS
CALL PRINT_MSG( NVERB_WARNING, 'GEN', 'INI_AIRCRAFT', &
'NMODEL is set to 1 at start for a CMODEL="MOB" aircraft (aircraft ' // TRIM( CTITLE(JI) ) // ')' )
END IF
+ IF ( NMODEL_NEST == 1 ) CMODEL(JI) = 'FIX' ! If only one model, FIX and MOB are the same
NMODEL(JI) = 1
ELSE
CMNHMSG(1) = 'invalid CMODEL (' // TRIM( CMODEL(JI) ) // ') for aircraft ' // TRIM( CTITLE(JI) )
@@ -167,10 +170,10 @@ DO JI = 1, NAIRCRAFTS
END IF
TAIRCRAFTS(JI)%TFLYER_TIME%XTSTEP = XTSTEP(JI)
- IF ( NPOS(JI) < 1 ) THEN
- CALL PRINT_MSG( NVERB_ERROR, 'GEN', 'INI_AIRCRAFT', 'NPOS should be at least 1 for aircraft ' // TRIM( CTITLE(JI) ) )
+ IF ( NPOS(JI) < 2 ) THEN
+ CALL PRINT_MSG( NVERB_ERROR, 'GEN', 'INI_AIRCRAFT', 'NPOS should be at least 2 for aircraft ' // TRIM( CTITLE(JI) ) )
END IF
- TAIRCRAFTS(JI)%NSEG = NPOS(JI)-1
+ TAIRCRAFTS(JI)%NPOS = NPOS(JI)
TAIRCRAFTS(JI)%LALTDEF = LALTDEF(JI)
@@ -179,13 +182,13 @@ DO JI = 1, NAIRCRAFTS
// TRIM( CTITLE(JI) ) )
! Allocate trajectory data
- ALLOCATE( TAIRCRAFTS(JI)%XSEGTIME(TAIRCRAFTS(JI)%NSEG ) ); TAIRCRAFTS(JI)%XSEGTIME(:) = XNEGUNDEF
- ALLOCATE( TAIRCRAFTS(JI)%XSEGLAT (TAIRCRAFTS(JI)%NSEG+1) ); TAIRCRAFTS(JI)%XSEGLAT(:) = XNEGUNDEF
- ALLOCATE( TAIRCRAFTS(JI)%XSEGLON (TAIRCRAFTS(JI)%NSEG+1) ); TAIRCRAFTS(JI)%XSEGLON(:) = XNEGUNDEF
+ ALLOCATE( TAIRCRAFTS(JI)%XPOSTIME(TAIRCRAFTS(JI)%NPOS) ); TAIRCRAFTS(JI)%XPOSTIME(:) = XNEGUNDEF
+ ALLOCATE( TAIRCRAFTS(JI)%XPOSLAT (TAIRCRAFTS(JI)%NPOS) ); TAIRCRAFTS(JI)%XPOSLAT(:) = XNEGUNDEF
+ ALLOCATE( TAIRCRAFTS(JI)%XPOSLON (TAIRCRAFTS(JI)%NPOS) ); TAIRCRAFTS(JI)%XPOSLON(:) = XNEGUNDEF
IF ( TAIRCRAFTS(JI)%LALTDEF ) THEN
- ALLOCATE( TAIRCRAFTS(JI)%XSEGP (TAIRCRAFTS(JI)%NSEG+1) ); TAIRCRAFTS(JI)%XSEGP(:) = XNEGUNDEF
+ ALLOCATE( TAIRCRAFTS(JI)%XPOSP (TAIRCRAFTS(JI)%NPOS) ); TAIRCRAFTS(JI)%XPOSP(:) = XNEGUNDEF
ELSE
- ALLOCATE( TAIRCRAFTS(JI)%XSEGZ (TAIRCRAFTS(JI)%NSEG+1) ); TAIRCRAFTS(JI)%XSEGZ(:) = XNEGUNDEF
+ ALLOCATE( TAIRCRAFTS(JI)%XPOSZ (TAIRCRAFTS(JI)%NPOS) ); TAIRCRAFTS(JI)%XPOSZ(:) = XNEGUNDEF
END IF
! Read CSV data (trajectory)
@@ -205,6 +208,7 @@ SUBROUTINE AIRCRAFT_CSV_READ( TPAIRCRAFT, HFILE )
USE MODD_AIRCRAFT_BALLOON, ONLY: TAIRCRAFTDATA
+USE MODE_DATETIME
USE MODE_MSG
IMPLICIT NONE
@@ -215,41 +219,39 @@ CHARACTER(LEN=*), INTENT(IN) :: HFILE !Name of the CSV file with the aircr
CHARACTER(LEN=NMAXLINELGT) :: YSTRING
INTEGER :: ILU ! logical unit of the file
INTEGER :: JI
-REAL :: ZTIME, ZLAT, ZLON, ZALT
-REAL :: ZTIME_OLD
-
-ZTIME_OLD = 0.
+REAL :: ZLAT, ZLON, ZALT
+REAL :: ZTIME
! Open file
OPEN( NEWUNIT = ILU, FILE = HFILE, FORM = 'formatted' )
READ( ILU, END = 101, FMT = '(A)' ) YSTRING ! Reading of header (skip it)
-DO JI = 1, TPAIRCRAFT%NSEG + 1
+DO JI = 1, TPAIRCRAFT%NPOS
! Read aircraft position
READ( ILU, END = 101, FMT = '(A)' ) YSTRING
READ( YSTRING, * ) ZTIME, ZLAT, ZLON, ZALT
- IF ( JI > 1 ) TPAIRCRAFT%XSEGTIME(JI-1) = ZTIME - ZTIME_OLD
- TPAIRCRAFT%XSEGLAT(JI) = ZLAT
- TPAIRCRAFT%XSEGLON(JI) = ZLON
+ TPAIRCRAFT%XPOSTIME(JI) = ZTIME
+ TPAIRCRAFT%XPOSLAT(JI) = ZLAT
+ TPAIRCRAFT%XPOSLON(JI) = ZLON
IF ( TPAIRCRAFT%LALTDEF ) THEN
- TPAIRCRAFT%XSEGP(JI) = ZALT * 100. ! *100 to convert from hPa to Pa
+ TPAIRCRAFT%XPOSP(JI) = ZALT * 100. ! *100 to convert from hPa to Pa
ELSE
- TPAIRCRAFT%XSEGZ(JI) = ZALT
+ TPAIRCRAFT%XPOSZ(JI) = ZALT
END IF
-
- ZTIME_OLD = ZTIME
END DO
101 CONTINUE
CLOSE( ILU )
-IF ( JI < TPAIRCRAFT%NSEG + 1 ) &
+IF ( JI < TPAIRCRAFT%NPOS ) &
CALL PRINT_MSG( NVERB_ERROR, 'GEN', 'AIRCRAFT_CSV_READ', 'Data not found in file ' // TRIM( HFILE ) )
+TPAIRCRAFT%TLAND = TPAIRCRAFT%TLAUNCH + TPAIRCRAFT%XPOSTIME(TPAIRCRAFT%NPOS)
+
END SUBROUTINE AIRCRAFT_CSV_READ
END MODULE MODE_INI_AIRCRAFT
diff --git a/src/MNH/ini_aircraft_balloon.f90 b/src/MNH/ini_aircraft_balloon.f90
index 7509894df..59da978f8 100644
--- a/src/MNH/ini_aircraft_balloon.f90
+++ b/src/MNH/ini_aircraft_balloon.f90
@@ -109,7 +109,6 @@ INTEGER :: ISTORE ! number of storage instants
INTEGER :: ILUOUT ! logical unit
INTEGER :: IRESP ! return code
INTEGER :: JI
-INTEGER :: JSEG ! loop counter
TYPE(TFIELDMETADATA) :: TZFIELD
!
!----------------------------------------------------------------------------
@@ -197,6 +196,7 @@ ENDIF
!
!
allocate( tpflyer%tflyer_time%tpdates(istore) )
+ALLOCATE(TPFLYER%NMODELHIST(ISTORE))
ALLOCATE(TPFLYER%XX (ISTORE))
ALLOCATE(TPFLYER%XY (ISTORE))
ALLOCATE(TPFLYER%XZ (ISTORE))
@@ -236,6 +236,7 @@ ALLOCATE(TPFLYER%XTHW_FLUX(ISTORE))
ALLOCATE(TPFLYER%XRCW_FLUX(ISTORE))
ALLOCATE(TPFLYER%XSVW_FLUX(ISTORE,KSV))
!
+TPFLYER%NMODELHIST = NNEGUNDEF
TPFLYER%XX = XUNDEF
TPFLYER%XY = XUNDEF
TPFLYER%XZ = XUNDEF
@@ -301,6 +302,7 @@ INTEGER(KIND=CDFINT) :: IGROUPID
INTEGER(KIND=CDFINT) :: ISTATUS
INTEGER(KIND=CDFINT), DIMENSION(2) :: IDATA ! Intermediate array to allow merge of 2 MPI broadcasts
#endif
+INTEGER :: IMODEL
LOGICAL :: GREAD ! True if balloon position was read in synchronous file
REAL :: ZLAT ! latitude of the balloon
REAL :: ZLON ! longitude of the balloon
@@ -310,8 +312,11 @@ TYPE(TFILEDATA) :: TZFILE
IF ( IMI /= TPFLYER%NMODEL ) RETURN
+LFLYER = .TRUE.
+
GREAD = .FALSE.
-LFLYER=.TRUE.
+
+CALL SM_XYHAT( PLATOR, PLONOR, TPFLYER%XLATLAUNCH, TPFLYER%XLONLAUNCH, TPFLYER%XXLAUNCH, TPFLYER%XYLAUNCH )
IF ( CPROGRAM == 'MESONH' .OR. CPROGRAM == 'SPAWN ' .OR. CPROGRAM == 'REAL ' ) THEN
! Read the current location in the synchronous file
@@ -490,13 +495,27 @@ IF ( CPROGRAM == 'MESONH' .OR. CPROGRAM == 'SPAWN ' .OR. CPROGRAM == 'REAL ' )
WRITE( CMNHMSG(2), * ) " Lat=", ZLAT, " Lon=", ZLON, " Rho=", TPFLYER%XRHO
END IF
CALL PRINT_MSG( NVERB_INFO, 'GEN', 'INI_LAUNCH' )
-
- TPFLYER%TFLYER_TIME%XTSTEP = MAX ( PTSTEP, TPFLYER%TFLYER_TIME%XTSTEP )
ELSE
! The position is not found, data is not in the synchronous file
! Use the position given in namelist
CALL PRINT_MSG( NVERB_INFO, 'GEN', 'INI_LAUNCH', 'initial location taken from namelist for ' // TRIM( TPFLYER%CTITLE ) )
END IF
+
+ ! Correct timestep if necessary
+ ! This has to be done at first pass (when IMI=1) to have the correct value as soon as possible
+ ! If 'MOB', set balloon store timestep to be at least the timestep of the coarser model (IMI=1) (with higher timestep)
+ ! as the balloon can fly on any model
+ ! If 'FIX', set balloon store timestep to be at least the timestep of its model
+ ! It should also need to be a multiple of the model timestep
+ IF ( IMI == 1 ) THEN
+ IF ( TPFLYER%CMODEL == 'MOB' ) THEN
+ IMODEL = 1
+ ELSE
+ IMODEL = TPFLYER%NMODEL
+ END IF
+
+ CALL FLYER_TIMESTEP_CORRECT( DYN_MODEL(IMODEL)%XTSTEP, TPFLYER )
+ END IF
!
ELSE IF ( CPROGRAM == 'DIAG ' ) THEN
IF ( LAIRCRAFT_BALLOON ) THEN
@@ -513,19 +532,10 @@ ELSE IF ( CPROGRAM == 'DIAG ' ) THEN
CALL PRINT_MSG( NVERB_INFO, 'GEN', 'INI_LAUNCH' )
END IF
!
- TPFLYER%TFLYER_TIME%XTSTEP = MAX (XSTEP_AIRCRAFT_BALLOON , TPFLYER%TFLYER_TIME%XTSTEP )
+ CALL FLYER_TIMESTEP_CORRECT( XSTEP_AIRCRAFT_BALLOON, TPFLYER )
END IF
END IF
-!
-IF ( TPFLYER%XLATLAUNCH == XUNDEF .OR. TPFLYER%XLONLAUNCH == XUNDEF ) THEN
- CMNHMSG(1) = 'Error in balloon initial position (balloon ' // TRIM( TPFLYER%CTITLE ) // ' )'
- CMNHMSG(2) = 'either LATitude or LONgitude is not given'
- CMNHMSG(3) = 'Check your INI_BALLOON routine'
- CALL PRINT_MSG( NVERB_FATAL, 'GEN', 'INI_AIRCRAFT_BALLOON' )
-END IF
-!
-CALL SM_XYHAT( PLATOR, PLONOR, TPFLYER%XLATLAUNCH, TPFLYER%XLONLAUNCH, TPFLYER%XXLAUNCH, TPFLYER%XYLAUNCH )
-!
+
END SUBROUTINE INI_LAUNCH
!----------------------------------------------------------------------------
!----------------------------------------------------------------------------
@@ -533,50 +543,61 @@ SUBROUTINE INI_FLIGHT(KNBR,TPFLYER)
!
INTEGER, INTENT(IN) :: KNBR
CLASS(TAIRCRAFTDATA), INTENT(INOUT) :: TPFLYER
-!
-IF (TPFLYER%CMODEL == 'MOB' .AND. TPFLYER%NMODEL /= 0) TPFLYER%NMODEL=1
-IF (TPFLYER%NMODEL > NMODEL) TPFLYER%NMODEL=0
+
+INTEGER :: IMODEL
+INTEGER :: JSEG ! loop counter
+
IF ( IMI /= TPFLYER%NMODEL ) RETURN
-!
+
LFLYER=.TRUE.
-!
-TPFLYER%TFLYER_TIME%XTSTEP = MAX ( PTSTEP, TPFLYER%TFLYER_TIME%XTSTEP )
-IF (TPFLYER%CTITLE==' ') THEN
- WRITE(TPFLYER%CTITLE,FMT='(A6,I2.2)') TPFLYER%CTYPE,KNBR
-END IF
+! Correct timestep if necessary
+! This has to be done at first pass (when IMI=1) to have the correct value as soon as possible
+! If 'MOB', set balloon store timestep to be at least the timestep of the coarser model (IMI=1) (with higher timestep)
+! as the balloon can fly on any model
+! If 'FIX', set balloon store timestep to be at least the timestep of its model
+! It should also need to be a multiple of the model timestep
+IF ( IMI == 1 ) THEN
+ IF ( TPFLYER%CMODEL == 'MOB' ) THEN
+ IMODEL = 1
+ ELSE
+ IMODEL = TPFLYER%NMODEL
+ END IF
-IF ( TPFLYER%NSEG == 0 ) THEN
- CMNHMSG(1) = 'Error in aircraft flight path (aircraft ' // TRIM( TPFLYER%CTITLE ) // ' )'
- CMNHMSG(2) = 'There is ZERO flight segment defined.'
- CMNHMSG(3) = 'Check your INI_AIRCRAFT routine'
- CALL PRINT_MSG( NVERB_FATAL, 'GEN', 'INI_FLIGHT' )
-END IF
-!
-IF ( ANY(TPFLYER%XSEGLAT(:)==XUNDEF) .OR. ANY(TPFLYER%XSEGLON(:)==XUNDEF) ) THEN
- CMNHMSG(1) = 'Error in aircraft flight path (aircraft ' // TRIM( TPFLYER%CTITLE ) // ' )'
- CMNHMSG(2) = 'either LATitude or LONgitude segment'
- CMNHMSG(3) = 'definiton is not complete.'
- CMNHMSG(4) = 'Check your INI_AIRCRAFT routine'
- CALL PRINT_MSG( NVERB_FATAL, 'GEN', 'INI_FLIGHT' )
+ CALL FLYER_TIMESTEP_CORRECT( DYN_MODEL(IMODEL)%XTSTEP, TPFLYER )
END IF
-!
-ALLOCATE(TPFLYER%XSEGX(TPFLYER%NSEG+1))
-ALLOCATE(TPFLYER%XSEGY(TPFLYER%NSEG+1))
-!
-DO JSEG=1,TPFLYER%NSEG+1
- CALL SM_XYHAT( PLATOR, PLONOR, TPFLYER%XSEGLAT(JSEG), TPFLYER%XSEGLON(JSEG), TPFLYER%XSEGX(JSEG), TPFLYER%XSEGY(JSEG) )
+
+ALLOCATE(TPFLYER%XPOSX(TPFLYER%NPOS))
+ALLOCATE(TPFLYER%XPOSY(TPFLYER%NPOS))
+
+DO JSEG = 1, TPFLYER%NPOS
+ CALL SM_XYHAT( PLATOR, PLONOR, TPFLYER%XPOSLAT(JSEG), TPFLYER%XPOSLON(JSEG), TPFLYER%XPOSX(JSEG), TPFLYER%XPOSY(JSEG) )
END DO
-!
-IF ( ANY(TPFLYER%XSEGTIME(:)==XUNDEF) ) THEN
- CMNHMSG(1) = 'Error in aircraft flight path (aircraft ' // TRIM( TPFLYER%CTITLE ) // ' )'
- CMNHMSG(2) = 'definiton of segment duration is not complete.'
- CMNHMSG(3) = 'Check your INI_AIRCRAFT routine'
- CALL PRINT_MSG( NVERB_FATAL, 'GEN', 'INI_AIRCRAFT_BALLOON' )
-END IF
END SUBROUTINE INI_FLIGHT
!----------------------------------------------------------------------------
+!----------------------------------------------------------------------------
+SUBROUTINE FLYER_TIMESTEP_CORRECT( PTSTEP_MODEL, TPFLYER )
+! Timestep is set to a multiple of the PTSTEP_MODEL value
+REAL, INTENT(IN) :: PTSTEP_MODEL
+CLASS(TFLYERDATA), INTENT(INOUT) :: TPFLYER
+
+REAL :: ZTSTEP_OLD
+
+ZTSTEP_OLD = TPFLYER%TFLYER_TIME%XTSTEP
+
+TPFLYER%TFLYER_TIME%XTSTEP = MAX ( PTSTEP_MODEL, TPFLYER%TFLYER_TIME%XTSTEP )
+TPFLYER%TFLYER_TIME%XTSTEP = NINT( TPFLYER%TFLYER_TIME%XTSTEP / PTSTEP_MODEL ) * PTSTEP_MODEL
+
+IF ( ABS( TPFLYER%TFLYER_TIME%XTSTEP - ZTSTEP_OLD ) > 1E-6 ) THEN
+ WRITE( CMNHMSG(1), '( "Timestep for flyer ", A, " is set to ", EN12.3, " (instead of ", EN12.3, ")" )' ) &
+ TPFLYER%CTITLE, TPFLYER%TFLYER_TIME%XTSTEP, ZTSTEP_OLD
+ CALL PRINT_MSG( NVERB_WARNING, 'GEN', 'INI_LAUNCH' )
+END IF
+
+END SUBROUTINE FLYER_TIMESTEP_CORRECT
+!----------------------------------------------------------------------------
+
!----------------------------------------------------------------------------
!
END SUBROUTINE INI_AIRCRAFT_BALLOON
diff --git a/src/MNH/ini_balloon.f90 b/src/MNH/ini_balloon.f90
index cda0db7c2..3c3c4d6ae 100644
--- a/src/MNH/ini_balloon.f90
+++ b/src/MNH/ini_balloon.f90
@@ -118,10 +118,12 @@ IMPLICIT NONE
INTEGER :: JI
-ALLOCATE( TBALLOONS (NBALLOONS) )
+ALLOCATE( TBALLOONS(NBALLOONS) )
!Treat balloon data read in namelist
DO JI = 1, NBALLOONS
+ TBALLOONS(JI)%NID = JI
+
IF ( CTITLE(JI) == '' ) THEN
WRITE( CTITLE(JI), FMT = '( A, I3.3) ') TRIM( CTYPE(JI) ), JI
@@ -143,6 +145,8 @@ DO JI = 1, NBALLOONS
CALL PRINT_MSG( NVERB_WARNING, 'GEN', 'INI_BALLOON', &
'NMODEL is set to 1 at start for a CMODEL="MOB" balloon (balloon ' // TRIM( CTITLE(JI) ) // ')' )
END IF
+ IF ( NMODEL_NEST == 1 ) CMODEL(JI) = 'FIX' ! If only one model, FIX and MOB are the same
+ ! NMODEL set temporarily to 1. Will be set to the launch model in INI_LAUNCH
NMODEL(JI) = 1
ELSE
CMNHMSG(1) = 'invalid CMODEL (' // TRIM( CMODEL(JI) ) // ') for balloon ' // TRIM( CTITLE(JI) )
diff --git a/src/MNH/modd_aircraft_balloon.f90 b/src/MNH/modd_aircraft_balloon.f90
index 4e71d254d..0ed8002ec 100644
--- a/src/MNH/modd_aircraft_balloon.f90
+++ b/src/MNH/modd_aircraft_balloon.f90
@@ -41,18 +41,23 @@
! ------------
!
!
-use modd_parameters, only: XNEGUNDEF, XUNDEF
+use modd_parameters, only: NNEGUNDEF, XNEGUNDEF, XUNDEF
USE MODD_TYPE_STATPROF, ONLY: TSTATPROFTIME
use modd_type_date, only: date_time
+USE MODE_DATETIME, ONLY: TPREFERENCE_DATE
+
implicit none
save
-!-------------------------------------------------------------------------------------------
-!
+INTEGER, PARAMETER :: NCRASH_NO = 0 ! Not crashed
+INTEGER, PARAMETER :: NCRASH_OUT_HORIZ = 1 ! Flyer is outside of horizontal domain
+INTEGER, PARAMETER :: NCRASH_OUT_LOW = 2 ! Flyer crashed on ground (or sea!)
+INTEGER, PARAMETER :: NCRASH_OUT_HIGH = 3 ! Flyer is too high (outside of domain)
+
LOGICAL :: LFLYER = .FALSE. ! flag to use aircraft/balloons
-!
+
TYPE :: TFLYERDATA
!
!* general information
@@ -62,18 +67,21 @@ TYPE :: TFLYERDATA
! 'MOB' : change od model depends of the
! balloon/aircraft location
INTEGER :: NMODEL = 0 ! model number for each balloon/aircraft
+ INTEGER :: NID = 0 ! Identification number
CHARACTER(LEN=6) :: CTYPE = '' ! flyer type:
! 'RADIOS' : radiosounding balloon
! 'ISODEN' : iso-density balloon
! 'AIRCRA' : aircraft
! 'CVBALL' : Constant Volume balloon
CHARACTER(LEN=10) :: CTITLE = '' ! title or name for the balloon/aircraft
- TYPE(DATE_TIME) :: TLAUNCH ! launch/takeoff date and time
+ TYPE(DATE_TIME) :: TLAUNCH = TPREFERENCE_DATE ! launch/takeoff date and time
LOGICAL :: LCRASH = .FALSE. ! occurence of crash
+ INTEGER :: NCRASH = NCRASH_NO
LOGICAL :: LFLY = .FALSE. ! occurence of flying
!
!* storage monitoring
!
+ LOGICAL :: LSTORE = .FALSE. ! Do we have to store data now
TYPE(TSTATPROFTIME) :: TFLYER_TIME ! Time management for flyer
!
!* current position of the balloon/aircraft
@@ -82,10 +90,11 @@ TYPE :: TFLYERDATA
REAL :: XY_CUR = XNEGUNDEF ! current y
REAL :: XZ_CUR = XNEGUNDEF ! current z (if 'RADIOS' or 'AIRCRA' and 'ALTDEF' = T)
REAL :: XP_CUR = XNEGUNDEF ! current p (if 'AIRCRA' and 'ALTDEF' = F)
+ INTEGER :: NRANK_CUR = NNEGUNDEF ! Rank of the process where the flyer is
!
!* data records
!
-
+ INTEGER, DIMENSION(:), ALLOCATABLE :: NMODELHIST ! List of models where data has been computed
REAL, DIMENSION(:), ALLOCATABLE :: XX ! X(n)
REAL, DIMENSION(:), ALLOCATABLE :: XY ! Y(n)
REAL, DIMENSION(:), ALLOCATABLE :: XZ ! Z(n)
@@ -126,16 +135,16 @@ TYPE, EXTENDS( TFLYERDATA ) :: TAIRCRAFTDATA
!
!* aircraft flight definition
!
- INTEGER :: NSEG = 0 ! number of aircraft flight segments
- INTEGER :: NSEGCURN = 1 ! current flight segment number
- REAL :: XSEGCURT = 0. ! current flight segment time spent
- REAL, DIMENSION(:), ALLOCATABLE :: XSEGLAT ! latitude of flight segment extremities (LEG+1)
- REAL, DIMENSION(:), ALLOCATABLE :: XSEGLON ! longitude of flight segment extremities (LEG+1)
- REAL, DIMENSION(:), ALLOCATABLE :: XSEGX ! X of flight segment extremities (LEG+1)
- REAL, DIMENSION(:), ALLOCATABLE :: XSEGY ! Y of flight segment extremities (LEG+1)
- REAL, DIMENSION(:), ALLOCATABLE :: XSEGP ! pressure of flight segment extremities (LEG+1)
- REAL, DIMENSION(:), ALLOCATABLE :: XSEGZ ! altitude of flight segment extremities (LEG+1)
- REAL, DIMENSION(:), ALLOCATABLE :: XSEGTIME ! duration of flight segments (LEG )
+ INTEGER :: NPOS = 0 ! number of aircraft positions (segment extremities)
+ INTEGER :: NPOSCUR = 1 ! current flight segment number
+ REAL, DIMENSION(:), ALLOCATABLE :: XPOSLAT ! latitude of flight segment extremities (LEG+1)
+ REAL, DIMENSION(:), ALLOCATABLE :: XPOSLON ! longitude of flight segment extremities (LEG+1)
+ REAL, DIMENSION(:), ALLOCATABLE :: XPOSX ! X of flight segment extremities (LEG+1)
+ REAL, DIMENSION(:), ALLOCATABLE :: XPOSY ! Y of flight segment extremities (LEG+1)
+ REAL, DIMENSION(:), ALLOCATABLE :: XPOSP ! pressure of flight segment extremities (LEG+1)
+ REAL, DIMENSION(:), ALLOCATABLE :: XPOSZ ! altitude of flight segment extremities (LEG+1)
+ REAL, DIMENSION(:), ALLOCATABLE :: XPOSTIME ! time since launch (corresponding to flight segments extremities (LEG+1)
+ TYPE(DATE_TIME) :: TLAND = TPREFERENCE_DATE ! landing / end of flight date and time
!
!* aircraft altitude type definition
!
@@ -159,6 +168,9 @@ TYPE, EXTENDS( TFLYERDATA ) :: TBALLOONDATA
REAL :: XINDDRAG = XNEGUNDEF ! induced drag coefficient (i.e. air shifted by the balloon) (if 'CVBALL')
REAL :: XVOLUME = XNEGUNDEF ! volume of the balloon (m3) (if 'CVBALL')
REAL :: XMASS = XNEGUNDEF ! mass of the balloon (kg) (if 'CVBALL')
+
+ TYPE(DATE_TIME) :: TPOS_CUR = TPREFERENCE_DATE ! Time corresponding to the current position (XX_CUR, XY_CUR...)
+
END TYPE TBALLOONDATA
INTEGER :: NAIRCRAFTS = 0 ! Total number of aircrafts
diff --git a/src/MNH/modeln.f90 b/src/MNH/modeln.f90
index 8183346f1..6de3470a4 100644
--- a/src/MNH/modeln.f90
+++ b/src/MNH/modeln.f90
@@ -2112,16 +2112,14 @@ IF (LFLYER) THEN
ALLOCATE(ZSEA(IIU,IJU))
ZSEA(:,:) = 0.
CALL MNHGET_SURF_PARAM_n (PSEA=ZSEA(:,:))
- CALL AIRCRAFT_BALLOON(XTSTEP, &
- XXHAT, XYHAT, XZZ, XMAP, XLONORI, XLATORI, &
- XUT, XVT, XWT, XPABST, XTHT, XRT, XSVT, XTKET, XTSRAD, &
- XRHODREF,XCIT,PSEA=ZSEA(:,:))
+ CALL AIRCRAFT_BALLOON( XTSTEP, XZZ, XMAP, XLONORI, XLATORI, &
+ XUT, XVT, XWT, XPABST, XTHT, XRT, XSVT, XTKET, XTSRAD, &
+ XRHODREF, XCIT, PSEA = ZSEA(:,:) )
DEALLOCATE(ZSEA)
ELSE
- CALL AIRCRAFT_BALLOON(XTSTEP, &
- XXHAT, XYHAT, XZZ, XMAP, XLONORI, XLATORI, &
- XUT, XVT, XWT, XPABST, XTHT, XRT, XSVT, XTKET, XTSRAD, &
- XRHODREF,XCIT)
+ CALL AIRCRAFT_BALLOON( XTSTEP, XZZ, XMAP, XLONORI, XLATORI, &
+ XUT, XVT, XWT, XPABST, XTHT, XRT, XSVT, XTKET, XTSRAD, &
+ XRHODREF, XCIT )
END IF
END IF
--
GitLab