coupla.f 13.1 KB
Newer Older
1
2
* METRo : Model of the Environment and Temperature of Roads
* METRo is Free and is proudly provided by the Government of Canada
Miguel Tremblay's avatar
Miguel Tremblay committed
3
* Copyright (C) Her Majesty The Queen in Right of Canada, Environment Canada, 2006
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36

*  Questions or bugs report: metro@ec.gc.ca
*  METRo repository: https://gna.org/projects/metro/
*  Documentation: http://documentation.wikia.com/wiki/METRo
*
*
* Code contributed by:
*  Louis-Philippe Crevier - Canadian meteorological center
*  Miguel Tremblay - Canadian meteorological center
*
*  $LastChangedDate$
*  $LastChangedRevision$
*
************************************************************************
*  This program is free software; you can redistribute it and/or modify
*  it under the terms of the GNU General Public License as published by
*  the Free Software Foundation; either version 2 of the License, or
*  (at your option) any later version.
*
*  This program is distributed in the hope that it will be useful,
*  but WITHOUT ANY WARRANTY; without even the implied warranty of
*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
*  GNU General Public License for more details.
*
*  You should have received a copy of the GNU General Public License
*  along with this program; if not, write to the Free Software
*  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*
*
*


***
Miguel Tremblay's avatar
Miguel Tremblay committed
37
38
*     Subroutine COUPLA: Do the coupling between the atmospheric forecast
*                        and the observations of the RWIS station.
39
*     Sous-routine COUPLA: effectue le couplage des previsions  avec 
40
41
42
43
44
*                          les observations des stations meteo-routieres
*
*     Auteur /  Author: Louis-Philippe Crevier
*     Date: Aout 1999 / August 1999
***
45
      SUBROUTINE COUPLA ( FS, FI, P0, TA , QA , VA , TYP, PR,
Miguel Tremblay's avatar
   
Miguel Tremblay committed
46
     *                    PVC , iref, NTP, NTP2, ITP, TSO, FLAT,
47
     *                    FCOR, ALN, ALR, FP,
48
     *                    FSCORR   , FICORR  , ER1, ER2, 
49
     *                    FAIL, EPSILON, Z0, Z0T, ZU, ZT, ECHEC,
50
     *                    dpRA, dpSN, npRC, dpRT, dpIR, dpSF, dpFV,
Miguel Tremblay's avatar
   
Miguel Tremblay committed
51
52
     *                    dpFC, dpFA, dpG, dpBB, dpFP, 
     *                    dpCapacity, dpConductivity)
53
54
55
56
57

      IMPLICIT NONE
      INTEGER i, j
      INTEGER Nl, n
      REAL DT
58
      INTEGER DTMAX
François Fortin's avatar
François Fortin committed
59
60
      REAL PADDING
      COMMON /BUFFER_SIZE/ DTMAX, Nl, DT, PADDING, n
61
62
63
64
65

***                 ***
*     DEFINITIONS     *
***                 ***
***
Miguel Tremblay's avatar
Miguel Tremblay committed
66
*     Input
67
*     -------
Miguel Tremblay's avatar
Miguel Tremblay committed
68
69
70
71
72
73
74
75
76
77
78
79
*     FS: Incident solar flux (W)
*     FI: Incident infra-red flux (W)
*     P0: Surface pressure (Pa)
*     TA: Air temperature (C) at level ZT
*     QA: Specific humidity (g/kg) at level ZT
*     VA: Wind (m/s) at level ZU
*     TYP: Probable type of precipitation ( 1 -> rain, 2 -> snow )
*     PR: Precipitation rate (m/s)
*     PVC: Road condition (0 [dry] ou 1 [wet])
*     iref: Number of levels in the grid
*     NTP: Index for the start of coupling
*     NTP2: Index for the end of coupling
80
*     TSO: Road surface temperature
Miguel Tremblay's avatar
Miguel Tremblay committed
81
82
83
84
85
86
87
*     Z0: Roughness length (m)
*     Z0T: Roughness length (m)
*     ZU: Height of level of the wind forecast (m)
*     ZT:  Height of level of the atmospheric forecast (m)
*     EPSILON: Road emissivity
*     TS0: Target temperature for the end of coupling (C)
*     FCOR: Coriolis factor
88
89
90
91
92
93
94
95
*     dpRA: Liquid accumulation on the road
*     dpSN: Solid (snow/ice) accumulation on the road
*     dpRC: Road condition
*     dpRT: Road temperature
*     dpIR: Infra-red flux
*     dpSF: Solar flux
*     dpFV: Vapor flux
*     dpFC: Sensible heat
96
*     dpFA: Anthropogenic flux
97
98
99
*     dpG:  Ground flux 
*     dpBB: Black body radiation
*     dpFP: Phase change energy
Miguel Tremblay's avatar
Miguel Tremblay committed
100
101
102
103
*     ALN: Snow Albedo 
*     ALR: Road Albedo 
*     FP: Frozing point (C)
*     FLAT: Road or bridge
Miguel Tremblay's avatar
   
Miguel Tremblay committed
104
105
*     dpCapacity: Thermic capacity of the road at every level
*     dpConductivity: Thermic conductivity of the road at every level
106
107
***
      LOGICAL FLAT
108
      INTEGER iref, NTP, NTP2, TYP(DTMAX), PVC(DTMAX)
109
110
      DOUBLE PRECISION  FS(DTMAX), FI(DTMAX), TA(DTMAX)
      DOUBLE PRECISION QA(DTMAX), PR(DTMAX)
111
112
      DOUBLE PRECISION VA(DTMAX), P0(DTMAX)
      DOUBLE PRECISION FCOR
113
      DOUBLE PRECISION ALN, ALR, TSO, FP
114
115
116
117
118
119
120
      DOUBLE PRECISION EPSILON, ZU, ZT, Z0, Z0T
      DOUBLE PRECISION dpSN(DTMAX), dpRA(DTMAX)
      DOUBLE PRECISION dpIR(DTMAX), dpSF(DTMAX)
      DOUBLE PRECISION dpFC(DTMAX), dpFA(DTMAX)
      DOUBLE PRECISION dpG(DTMAX), dpBB(DTMAX)
      DOUBLE PRECISION dpRT(DTMAX), dpFV(DTMAX)
      DOUBLE PRECISION dpFP(DTMAX)
Miguel Tremblay's avatar
   
Miguel Tremblay committed
121
      DOUBLE PRECISION dpCapacity(DTMAX), dpConductivity(DTMAX)
122
      INTEGER npRC(DTMAX)
123
      INTEGER nCheckBefore, nCheckAfter, n30SecondsStepsIn3Hours
124
***
Miguel Tremblay's avatar
Miguel Tremblay committed
125
*     I/O
126
*     ---------------
Miguel Tremblay's avatar
Miguel Tremblay committed
127
*     ITP: Road temperature profile  (C, iref levels != 0)
128
129
130
131
***
      LOGICAL ECHEC
      DOUBLE PRECISION ITP(n)
***
Miguel Tremblay's avatar
Miguel Tremblay committed
132
*     Output
133
*     -------
Miguel Tremblay's avatar
Miguel Tremblay committed
134
135
136
137
138
139
*     FAIL: Tells if the coupling have failed or not
*     TS: Time series for the surface temperature
*     FSCORR: Coupling Coefficient for the solar flux
*     FICORR: Coupling Coefficient for the infra-red flux
*     ER1: Water quantity of the road (mm)
*     ER2: Water Equivalent for ice/snow on the road (mm)
140
141
142
143
144
145
146
147
148
149
150
***
      LOGICAL FAIL
      DOUBLE PRECISION  FSCORR, FICORR, ER1, ER2
***
*     Internes
*     --------
***
      CHARACTER*250 outfmt
      INTEGER next, now, iter
      LOGICAL DOWN
      REAL T(n,2), G(0:n)
151
      REAL QG, RHO, TSK, COUDUR
152
      DOUBLE PRECISION AL, M, RA, DIV
153
154
155
      REAL CL, PR1, PR2, DX, PRG, FZ
      REAL coeff1, coeff2, cotemp, COFS, COFI
***
Miguel Tremblay's avatar
Miguel Tremblay committed
156
*     Variables of FLXSURFZ
157
158
159
160
161
*     ---------------------
***
      REAL CMU, CTU, RIB,  ILMO, FTEMP, FVAP
      REAL H, UE, LZZ0, LZZ0T, FM, FH
***
Miguel Tremblay's avatar
Miguel Tremblay committed
162
*     Initialisation of common CLELOG
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
*     -------------------------------
***

      REAL  TCDK
      REAL CPD, CPV, RGASD, RGASV, TRPL, RAUW, EPS1
      REAL DELTA, CAPPA, TGL, CONSOL, GRAV, RAYT, STEFAN, PI
      REAL OMEGA, EPS2
      REAL KNAMS, STLO, KARMAN, RIC, CHLC, CHLF, DG2RAD

      COMMON/CTSPHY/  CPD, CPV, RGASD, RGASV, TRPL, TCDK, RAUW,
     $                EPS1, EPS2, DELTA, CAPPA, TGL, CONSOL,
     $                GRAV, RAYT, STEFAN, PI, OMEGA,
     $                KNAMS, STLO, KARMAN, RIC, CHLC, CHLF, DG2RAD


178
179
      REAL TTT,QQQ
      REAL FOTVT
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
      LOGICAL bSilent
      COMMON /SILENT/ bSilent

      INTRINSIC SIGN

*   DEFINITION DES FONCTIONS THERMODYNAMIQUES DE BASE
*   POUR LES CONSTANTES, UTILISER LE COMMON /CTESDYN/
*     NOTE: TOUTES LES FONCTIONS TRAVAILLENT AVEC LES UNITES S.I.
*           I.E. TTT EN DEG K, PRS EN PA, QQQ EN KG/KG
*          *** N. BRUNET - MAI 90 ***
*          * REVISION 01 - MAI 94 - N. BRUNET
*                          NOUVELLE VERSION POUR FAIBLES PRESSIONS
*
*      FONCTION CALCULANT TEMP VIRT. (TVI) DE TEMP (TTT) ET HUM SP (QQQ)
      FOTVT(TTT,QQQ) = TTT * (1.0 + DELTA*QQQ)
*
*
***
*
*     Procedure
*     =========
      if( .not. bSilent) then
         WRITE(*,*) "DEBUT COUPLA"
      end if
*****
Miguel Tremblay's avatar
Miguel Tremblay committed
205
**     Initialisation of parametres
206
207
208
209
210
211
212
213
**     ++++++++++++++++++++++++++++
      FAIL = .FALSE.
      DOWN = .FALSE.
      FZ = 0.1
      iter = 0
      coeff1 = 0.
      coeff2 = 1.
      ILMO = 1.
Miguel Tremblay's avatar
Miguel Tremblay committed
214
**    Compute the length of coupling
215
      COUDUR = real(NTP2-NTP+1)*DT/3.6e3
Miguel Tremblay's avatar
Typo    
Miguel Tremblay committed
216
*     Initialization of albedo variation due to snow presence
217
218
*     ---------------------------------------------------------------
*     |                                                             |
Miguel Tremblay's avatar
Miguel Tremblay committed
219
*     |  See balanc.F for the description of the function used      |
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
*     |                                                             |
*     ---------------------------------------------------------------
      M = (ALN - ALR) / 5.
 11   next = 1
      now = 2
      H = 300.0
      ER1 = 0.0
      ER2 = 0.0
      iter = iter + 1
      G(0) = 0.0
      do j=1,iref
         T(j,now)=REAL(ITP(j))
         T(j,next)=REAL(ITP(j))
         G(j) = 0.
      end do

************************************************
*    Criterium: Solar flux must be non zero
*     during the 3 hours BEFORE the roadcast and 
*     1 hour AFTER the begining of the roadcast.
240
241
242
243
*       3*60*2 = 360
      n30SecondsStepsIn3Hours = 360
      if (NTP2 - n30SecondsStepsIn3Hours > NTP) then
         nCheckBefore = NTP2-n30SecondsStepsIn3Hours 
244
245
246
      else
         nCheckBefore = NTP
      end if
247
248
      if (NTP2 + n30SecondsStepsIn3Hours/3 < DTMAX) then
         nCheckAfter = NTP2+n30SecondsStepsIn3Hours/3
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
      else
         nCheckAfter = DTMAX
      end if

      if(FS(nCheckBefore) > 0.0 .and.
     *   FS(NTP2) > 100.0 .and. 
     *   FS(nCheckAfter) > 0.0) then
         COFI = 1.0
         COFS = coeff2
      else
*        Nuit
         COFI = coeff2
         COFS = 1.0
      end if

*     +++++++++++++++++++++++++++++++++
Miguel Tremblay's avatar
Miguel Tremblay committed
265
*     Entry of the coupling loop
266
267
*     +++++++++++++++++++++++++++++++++
      do i = NTP+1,NTP2
Miguel Tremblay's avatar
Miguel Tremblay committed
268
*        Adjustement of the road conditions
269
270
271
*        +++++++++++++++++++++++++++++++++++++
         ER1 = PVC(i)*ER1
         ER2 = PVC(i)*ER2
Miguel Tremblay's avatar
Miguel Tremblay committed
272
*        Preparation of differents terms of balance
273
274
275
*        ++++++++++++++++++++++++++++++++++++++++++
*        road temp in Kelvins         
         TSK = T(1,now) + TCDK
Miguel Tremblay's avatar
Miguel Tremblay committed
276
*        Black body radiation 
277
         RA = EPSILON*STEFAN*TSK**4
Miguel Tremblay's avatar
Miguel Tremblay committed
278
*        Humidity at the surface
279
         call SRFHUM  ( QG, CL, ER1, ER2, TSK, P0(i), QA(i), FP )
Miguel Tremblay's avatar
Miguel Tremblay committed
280
*        Air density at ground
281
         DIV = RGASD * FOTVT ( TSK , QG )
282
         if (  abs(DIV) < 1.0D-5 ) then
283
            ECHEC = .true.
284
            WRITE(*,*) "Echec in coupla.f"
285
286
            return
         end if
287
         RHO = REAL(P0(i) / ( RGASD * FOTVT ( TSK , QG ) ))
Miguel Tremblay's avatar
Miguel Tremblay committed
288
*        Energy used/free by the melting snow/freezing rain
289
290
         call VERGLAS ( TYP(i), T(1,now), FP, FZ, PR(i),
     *                  PR1, PR2, PRG, FAIL )
Miguel Tremblay's avatar
Miguel Tremblay committed
291
*        Heat flux coefficients
292
293
294
295
296
         call FLXSURFZ( CMU  , CTU , RIB  , FTEMP, FVAP ,
     *                  ILMO, UE   , FCOR , TA(i)+TCDK,
     *                  QA(i), ZU, ZT  , VA(i), TSK  , QG   ,
     *                  H    , Z0  , Z0T  , LZZ0 , LZZ0T,
     *                  FM   , FH  ,  1,  1     )
Miguel Tremblay's avatar
Miguel Tremblay committed
297
*        Effective albedo 
298
         AL = 1.0 - max( ALR,min( ALN,M*ER2+ALR-M ) )
Miguel Tremblay's avatar
Miguel Tremblay committed
299
*        Energetic balance computation
300
*        +++++++++++++++++++++++++++
301
         G(0) = REAL(max(0.0,REAL(COFS*AL*FS(i))) - RA 
302
303
     *        + COFI*EPSILON*FI(i)
     *        + RHO*CTU*( CPD*(TA(i)-T(1,now))
304
     *        + CL*(QA(i)-QG)) + PRG + dpFA(i))
305
306
307
308
309
310
311
312
313
*        Output information         
         dpIR(i) = COFI*EPSILON*FI(i)
         dpSF(i) = COFS*AL*FS(i)
         dpFV(i) = RHO*CTU*(CL*(QA(i)-QG))
         dpBB(i) = - RA
         dpFC(i) =  RHO*CTU*( CPD*(TA(i)-T(1,now)))
         dpFP(i) = PRG
         dpG(i) = G(0)

Miguel Tremblay's avatar
   
Miguel Tremblay committed
314
315
         G(1) = REAL(dpConductivity(1) * ( T(2,now) - T(1,now) ))
         T(1,next) = REAL(T(1,now)+DT*(dpCapacity(1)*( G(1) - G(0)) ))
316
317

         DX = 0.0
Miguel Tremblay's avatar
Miguel Tremblay committed
318
*        Compute the evolution of temperature in the ground
319
*        +++++++++++++++++++++++++++++++++++++++++++++++++
Miguel Tremblay's avatar
   
Miguel Tremblay committed
320
321
         call TSEVOL ( T, iref, now, G, FLAT, TA(i), 
     *        dpCapacity, dpConductivity)
Miguel Tremblay's avatar
Miguel Tremblay committed
322
*        Balance of accumulation at the ground
323
324
*        ++++++++++++++++++++++++++++++
         call RODCON ( ER1, ER2, RHO  , CTU, CL , FP , FZ , 
325
     *                 T(1,now), QA(i), QG , PR1, PR2, PRG, DX , 
326
327
328
329
330
331
     *                 npRC(i))

         dpRT(i) = T(1,next)
         dpRA(i) = ER1
         dpSN(i) = ER2

Miguel Tremblay's avatar
Miguel Tremblay committed
332
*        Inversion of indices
333
334
335
336
337
338
339
340
341
342
343
*        +++++++++++++++++++++
         next = 3 - next
         now = 3 - now
      end do
*      if ( VERBOSE ) then
         FSCORR = COFS - 1.0
         FICORR = COFI - 1.0
         write(outfmt,223) iter, COUDUR,
     *        min(99.99,max(-99.99,T(1,now))),
     *        min(99.99,max(-99.99,REAL(TSO))),COFI,COFS
*      end if
Miguel Tremblay's avatar
Miguel Tremblay committed
344
*     Check if other iterations are necessary
345
346
347
348
349
*     +++++++++++++++++++++++++++++++++++++++++++
      if ( iter .eq. 25 ) then
         FICORR = 0.0
         FSCORR = 0.0
         FAIL = .true.
350
         WRITE(*,*) "More than 25 iterations in the coupling phase "
351
      else if ( T(1,now) < -99.99 .or. T(1,now) > 99.99 .or.
352
353
354
355
356
357
358
     *        FAIL ) then
         FICORR = 0.0
         FSCORR = 0.0
         write(outfmt,223) iter, COUDUR,
     *        min(99.99,max(-99.99,T(1,now))),
     *        min(99.99,max(-99.99,REAL(TSO))),COFI,COFS
         FAIL = .true.
359
         WRITE(*,*) " Temperature over 99 degrees, aborting in coupling"
360
      else if ( mod(iter,15) .eq. 0 .and.
361
362
     *        abs(min(REAL(max(T(1,now),REAL(TSO-0.1))),REAL(TSO+0.1)) 
     *    - T(1,now)) > 1.0D-5 ) then
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
         coeff1 = 0.0
         coeff2 = 1.5*coeff2
         go to 11
      else if ( T(1,now) - TSO .gt. 0.1 ) then
         coeff2 = coeff2 - 0.5*(coeff2-coeff1)
         DOWN = .true.
         go to 11
      else if ( TSO - T(1,now) .gt. 0.1 ) then
        if ( DOWN ) then
           cotemp = coeff1
           coeff1 = coeff2
           coeff2 = coeff2 + 0.5*(coeff2-cotemp)
         else
            cotemp = coeff1
            coeff1 = coeff2
            coeff2 = 2.0*coeff2 - cotemp
         end if
         go to 11
      else
         FSCORR = COFS - 1.0
         FICORR = COFI - 1.0
         FAIL = .false.
         do j=1,iref
            ITP(j) = T(j,now)
         end do
      end if
 223  format(1x,i3,3x,f5.2,3x,f6.2,3x,f6.2,f9.5,f9.5)

      if( .not. bSilent) then
         WRITE(*,*) "FIN COUPLA"
      end if

      RETURN
      END