Commit 4c6fae37 authored by Miguel Tremblay's avatar Miguel Tremblay
Browse files

Conversion of road condition from double to long.

Improve comparison handling in fortran code, using < 1.0D-5 instead of ".eq. 0" while comparing different data type.
parent 638c96ca
......@@ -228,8 +228,9 @@ class Metro_model(Metro_module):
lIR = wf_interpolated_data.get_matrix_col('IR').tolist()
npFA = wf_interpolated_data.get_matrix_col('FA')
lFA = npFA.tolist()
lPI = wf_interpolated_data.get_matrix_col('PI').tolist()
lSC = wf_interpolated_data.get_matrix_col('SC').tolist()
lPI = wf_interpolated_data.get_matrix_col('PI').astype(numpy.int32).tolist()
lSC = wf_interpolated_data.get_matrix_col('SC').astype(numpy.int32).tolist()
# Number of 30 seconds step.
npFT = wf_interpolated_data.get_matrix_col('Time')
......@@ -306,14 +307,13 @@ class Metro_model(Metro_module):
bEchec = []
macadam.Do_Metro(bFlat,\
fLat, fLon, \
lLayerThick, nNbrOfLayer, lLayerType, \
lAT, lQP, lWS, lAP, lSF, lIR, lFA, lPI, lSC,\
lAT_obs,lST_obs, lSST_obs, \
lAH, lTime_obs, lSWO, bNoObs,\
fDeltaTMetroObservation, nLenObservation, \
nNbrTimeSteps, bSilent, \
macadam.Do_Metro(bFlat, fLat, fLon, lLayerThick, \
nNbrOfLayer, lLayerType, lAT, lQP, \
lWS, lAP, lSF, lIR, \
lFA, lPI, lSC, lAT_obs, \
lST_obs, lSST_obs, lAH, lTime_obs, \
lSWO, bNoObs, fDeltaTMetroObservation, \
nLenObservation, nNbrTimeSteps, bSilent, \
dSstDepth, bDeepTemp, dDeepTemp)
bEchec = (macadam.get_echec())[0]
# Check if the execution of the model was a succes:
......
......@@ -111,10 +111,10 @@
* dpConductivity: Thermic conductivy of the road at every level
***
LOGICAL FLAT
INTEGER iref, ir40, NTP, NTFM
INTEGER iref, ir40, NTP, NTFM, TYP(DTMAX)
DOUBLE PRECISION FS(DTMAX), FI(DTMAX), P0(DTMAX)
DOUBLE PRECISION TA(DTMAX), QA(DTMAX), VA(DTMAX)
DOUBLE PRECISION TYP(DTMAX), PR(DTMAX)
DOUBLE PRECISION PR(DTMAX)
DOUBLE PRECISION ITP(n)
DOUBLE PRECISION FCOR
DOUBLE PRECISION ALN, ALR, FP
......@@ -162,7 +162,7 @@
EXTERNAL CHKDIV
INTEGER next, now
REAL T(n,2), G(0:n)
DOUBLE PRECISION FGD, AL, RA, M
DOUBLE PRECISION FGD, AL, RA, M, DIV
REAL QG, RHO, TSK
REAL CL, PR1, PR2, DX, PRG, FZ
DOUBLE PRECISION COFS, COFI
......@@ -310,7 +310,8 @@
* humidite at the suface
call SRFHUM ( QG, CL, ER1, ER2, TSK, P0(i), QA(i), FP )
* Air density at the surface
ier = CHKDIV ( RGASD * FOTVT ( TSK , QG ), "balanc.ftn", 196 )
DIV = RGASD * FOTVT ( TSK , QG )
ier = CHKDIV (DIV, "balanc.f", 196 )
if ( ier .eq. 1 ) then
WRITE(*,*) "Echec dans balanc.f"
ECHEC = .true.
......
......@@ -105,11 +105,10 @@
* dpConductivity: Thermic conductivity of the road at every level
***
LOGICAL FLAT
INTEGER iref, NTP, NTP2
INTEGER iref, NTP, NTP2, TYP(DTMAX), PVC(DTMAX)
DOUBLE PRECISION FS(DTMAX), FI(DTMAX), TA(DTMAX)
DOUBLE PRECISION QA(DTMAX), PR(DTMAX)
DOUBLE PRECISION VA(DTMAX), P0(DTMAX)
DOUBLE PRECISION TYP(DTMAX), PVC(DTMAX)
DOUBLE PRECISION FCOR
DOUBLE PRECISION ALN, ALR, TSO, FP
DOUBLE PRECISION EPSILON, ZU, ZT, Z0, Z0T
......@@ -152,7 +151,7 @@
LOGICAL DOWN
REAL T(n,2), G(0:n)
REAL QG, RHO, TSK, COUDUR
DOUBLE PRECISION AL, M, RA
DOUBLE PRECISION AL, M, RA, DIV
REAL CL, PR1, PR2, DX, PRG, FZ
REAL coeff1, coeff2, cotemp, COFS, COFI
***
......@@ -326,9 +325,11 @@
* Humidity at the surface
call SRFHUM ( QG, CL, ER1, ER2, TSK, P0(i), QA(i), FP )
* Air density at ground
ier = CHKDIV ( RGASD * FOTVT ( TSK , QG ), "coupla.ftn", 209 )
DIV = RGASD * FOTVT ( TSK , QG )
ier = CHKDIV (DIV, "coupla.f", 209 )
if ( ier .eq. 1 ) then
ECHEC = .true.
WRITE(*,*) "Echec in coupla.f"
return
end if
RHO = REAL(P0(i) / ( RGASD * FOTVT ( TSK , QG ) ))
......@@ -395,7 +396,7 @@
FSCORR = 0.0
FAIL = .true.
WRITE(*,*) "More than 25 iterations in the coupling phase "
else if ( abs(min(99.99,max(-99.99,T(1,now)))) .eq. 99.99 .or.
else if ( T(1,now) < -99.99 .or. T(1,now) > 99.99 .or.
* FAIL ) then
FICORR = 0.0
FSCORR = 0.0
......@@ -405,8 +406,8 @@
FAIL = .true.
WRITE(*,*) " Temperature over 99 degrees, aborting in coupling"
else if ( mod(iter,15) .eq. 0 .and.
* min(REAL(max(T(1,now),REAL(TSO-0.1))),REAL(TSO+0.1))
* .ne. T(1,now) ) then
* abs(min(REAL(max(T(1,now),REAL(TSO-0.1))),REAL(TSO+0.1))
* - T(1,now)) > 1.0D-5 ) then
coeff1 = 0.0
coeff2 = 1.5*coeff2
go to 11
......
......@@ -177,13 +177,13 @@
* Auteur/Author: Louis-Philippe Crevier
* Date: Octobre 2000
***
FUNCTION CHKDIV ( DIV, FICH, LIGNE )
FUNCTION CHKDIV (DIV, FICH, LIGNE )
IMPLICIT NONE
INTEGER CHKDIV
INTEGER a
CHARACTER*(*) FICH, LIGNE
DOUBLE PRECISION DIV
if ( DIV .eq. 0.0 ) then
if ( abs(DIV) < 1.0D-5 ) then
a=index(FICH,' ')-1
FICH=FICH(1:a)//', '//LIGNE
CHKDIV = 1
......
......@@ -49,8 +49,9 @@
***
* Entrees
* -------
INTEGER TYP
REAL TS, FZ
DOUBLE PRECISION TYP, FP, PR
DOUBLE PRECISION FP, PR
***
* Sorties
* -------
......@@ -66,6 +67,8 @@
* --------
***
DOUBLE PRECISION RATIO
INTEGER RAIN, SNOW
***
* Declarations des constantes physique
* et des fonctions thermodynamiques
......@@ -87,36 +90,39 @@
* Procedure
* =========
RATIO = 0.0
if ( TYP.eq.2. .and. TS.le.FP ) then
* To be used to identify TYP
RAIN =1
SNOW =2
if ( TYP .eq. SNOW .and. TS.le.FP ) then
* neige sur route a T < FP
PR1 = 0.0
PR2 = REAL(1e3 * PR)
PRG = 0.0
else if ( TYP.eq.1. .and. TS.gt.FP ) then
else if ( TYP .eq. RAIN .and. TS.gt.FP ) then
* pluie sur route a T > FP
PR1 = REAL(1e3 * PR)
PR2 = 0.0
PRG = 0.0
else if ( TYP.eq.2..and. TS.gt.FP+FZ ) then
else if ( TYP .eq. SNOW .and. TS.gt.FP+FZ ) then
* neige sur route a T > FZ; neige fond au contact de la route
PR1 = REAL(1e3 * PR)
PR2 = 0.0
PRG = REAL(- (PR * CHLF * RAUW))
else if ( TYP.eq.1. .and. TS.le.FP-FZ ) then
else if ( TYP .eq. RAIN .and. TS.le.FP-FZ ) then
* pluie sur route a T < FZ; pluie gele au contact de la route
PR1 = 0.0
PR2 = REAL(1e3 * PR)
PRG = REAL(PR * CHLF * RAUW)
else if ( TYP.eq.2. .and.
* max(min(TS,REAL(FP+FZ)),REAL(FP)).eq.TS ) then
else if ( TYP .eq. SNOW .and.
* abs(max(min(TS,REAL(FP+FZ)),REAL(FP)) - TS) < 1.0D-5 ) then
* neige sur route a FZ > T > 0; une partie de la neige fond au
* contact de la route
RATIO = 0.25*sin(PI*(TS-FP)/FZ)+0.75
PR1 = REAL(RATIO * 1e3 * PR)
PR2 = REAL((1.0-RATIO) * 1e3 * PR)
PRG = REAL(- (RATIO * PR * CHLF * RAUW))
else if ( TYP.eq.1. .and.
* min(max(TS,REAL(FP-FZ)),REAL(FP)).eq.TS ) then
else if ( TYP .eq. RAIN .and.
* abs(min(max(TS,REAL(FP-FZ)),REAL(FP)) - TS) < 1.0D-5 ) then
* pluie sur route a -FZ < T < 0; une partie de la pluie gele
* au contact de la route
RATIO = 0.25*sin(PI*(TS-FP)/FZ)+0.75
......@@ -333,13 +339,13 @@
* -----------------------------------------------------------
ETAT = 6
else if ( ER1.lt.CUTOFF .and. ER2.lt.CUTOFF ) then
if ( FV1.gt.0.0 .and. PR1 .eq. 0.0 ) then
if ( FV1 > 0.0 .and. abs(PR1) < 1.0D-5 ) then
* -----------------------------------------------------------
* 5. -> Dew
* -----------------------------------------------------------
ETAT = 5
else if ( ( FV2 .gt. 0.0 .and. PR2 .eq. 0.0 ) .or.
else if ( ( FV2 .gt. 0.0 .and. abs(PR2) < 1.0D-5 ) .or.
* DX .lt. 0.0 ) then
* -----------------------------------------------------------
* 7. -> *** Black ice warning ***
......
......@@ -99,7 +99,7 @@ static struct doubleStruct stLT; /* Level temperature */
[I double* dpFI : interpolated visible flux]
[I double* dpFA : interpolated anthropogenic flux]
[I double* TYP : Type of precipitation: 0 = nada, 1= liquid, 2=solid]
[I double* dpRC : interpolated road condition. 0 = dry, 1=wet]
[I long* npRC : interpolated road condition. 0 = dry, 1=wet]
[I double dpTAO : interpolated observed air temperature]
[I double* dpRTO : interpolated observed road temperature]
[I double* dpDTO : interpolated observed deep road temperature]
......@@ -149,7 +149,7 @@ Miguel Tremblay May 2004
void Do_Metro( BOOL bFlat, double dMLat, double dMLon, double* dpZones,\
long nNbrOfZone, long* npMateriau, double* dpTA, double* dpQP,\
double* dpFF, double* dpPS, double* dpFS, double* dpFI, \
double* dpFA, double* dpTYP, double* dpRC, double* dpTAO, \
double* dpFA, long* npTYP, long* npRC, double* dpTAO, \
double* dpRTO, double* dpDTO, double* dpAH, double* dpTimeO,\
long* npSwo, BOOL* bpNoObs, double dDeltaT,\
long nLenObservation, long nNbrTimeSteps, BOOL bSilent,\
......@@ -284,7 +284,7 @@ void Do_Metro( BOOL bFlat, double dMLat, double dMLon, double* dpZones,\
&bDeepTemp, &dDeepTemp);
nNtp = - nDeltaTIndice + nNtdcl;
nNtp2 = nLenObservation - nDeltaTIndice;
f77name(coupla)(dpFS, dpFI, dpPS, dpTA, dpAH, dpFF, dpTYP, dpQP, dpRC, \
f77name(coupla)(dpFS, dpFI, dpPS, dpTA, dpAH, dpFF, npTYP, dpQP, npRC, \
&stTemperatureDepth.nSize, &nNtp, &nNtp2, dpItp, \
&(dpRTO[nLenObservation]), &bFlat, &dFCorr, \
&dAln, &dAlr, &dFp, &dFsCorr, &dFiCorr, &dEr1, &dEr2, \
......@@ -323,8 +323,8 @@ void Do_Metro( BOOL bFlat, double dMLat, double dMLon, double* dpZones,\
&nIR40, &bFlat, npSwo, dpCapacity, dpConductivity);
nNtp = 0 + nNtdcl;
nNtp2 = nLenObservation - nDeltaTIndice;
f77name(coupla)(dpFS, dpFI, dpPS, dpTA, dpAH, dpFF, dpTYP, dpQP, \
dpRC, &stTemperatureDepth.nSize, &nNtp, &nNtp2, dpItp,\
f77name(coupla)(dpFS, dpFI, dpPS, dpTA, dpAH, dpFF, npTYP, dpQP, \
npRC, &stTemperatureDepth.nSize, &nNtp, &nNtp2, dpItp,\
&(dpRTO[nLenObservation]), &bFlat, &dFCorr, \
&dAln, &dAlr, &dFp, &dFsCorr, &dFiCorr, &dEr1, &dEr2,\
&bFail, &dEpsilon, &dZ0, &dZ0t, &dZu, &dZt, stEc.plArray,\
......@@ -348,7 +348,7 @@ void Do_Metro( BOOL bFlat, double dMLat, double dMLon, double* dpZones,\
}/* End else observation complete */
/************ roadcast **************************************************/
f77name(balanc)(dpFS, dpFI, dpPS, dpTA, dpAH, dpFF, dpTYP, dpQP,\
f77name(balanc)(dpFS, dpFI, dpPS, dpTA, dpAH, dpFF, npTYP, dpQP,\
&stTemperatureDepth.nSize, \
&nIR40, &nNtp2, &nNbrTimeSteps, dpItp, &bFlat, &dFCorr,\
&dAln, &dAlr, &dFp, &dFsCorr, &dFiCorr, &dEr1,\
......
......@@ -39,7 +39,7 @@
void Do_Metro(BOOL bFlat, double dMLat, double dMLon, double* dpZones, \
long nNbrOfZone, long* npMateriau, double* dpTA, double* dpQP,\
double* dpFF, double* dpPS, double* dpFsPy, double* dpFI, \
double* dpFA, double* dpTYP, double* dpRc, double* dpTAO, \
double* dpFA, long* npTYP, long* npRc, double* dpTAO, \
double* dpRTO, double* dpDTO, double* dpAH, double* dpTimeO,\
long* npSWO, BOOL* bpNoObs, double dDeltaT, \
long nLenObservation, long nNbrTimeSteps, BOOL bSilent,\
......@@ -61,7 +61,7 @@ extern void initial_(double* dpItp , double* dpRTO, double* dpDTO, double* dpTAO
long* nLenObservation, long* nSize, long* nIR40, \
BOOL* bFlat, long* npSwo, double* dpCapacity, double* dpConductivity);
extern void coupla_(double* dpFS, double* dpFI, double* dpPS, double* dpTA, double* dpAH, \
double* dpFF, double* dpTYP, double* dpQP, double* dpRC, \
double* dpFF, long* npTYP, double* dpQP, long* npRC, \
long* nSize, long* nNtp, long* nNtp2, double* dpItp, double* dpRTO, \
BOOL* bFlat, double* dFCorr, double* dAln, double* dAlr, \
double* dFp, double* dFsCorr, double* dFiCorr, double* dEr1, double* dEr2, \
......@@ -73,7 +73,7 @@ extern void coupla_(double* dpFS, double* dpFI, double* dpPS, double* dpTA, doub
extern void balanc_(double* dpFS, double* dpFI, double* dpPS, double* dpTA, double* dpAH, \
double* dpFF, double* dpTYP, double * dpQP, \
double* dpFF, long* npTYP, double * dpQP, \
long* nSize, long* nIR40, long* nNtp2, long* nNbrTimeSteps, double* dpItp,\
BOOL* bFlat, double* dFCorr, double* dAln, double* dAlr, \
double* dFp, double* dFsCorr, double* dFiCorr, double* dEr1, \
......
......@@ -393,7 +393,7 @@ int print_args(char **argv) {
extern void Do_Metro(long, double, double, double*,\
long, long*, double*, double*,\
double*, double*, double*, double*,\
double*, double*, double*, double*,\
double*, long*, long*, double*,\
double*, double*, double*, double*,\
long*, long*, double,\
long, long, long, \
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment