The TELEMAC-MASCARET system  trunk
thermal_khione.f
Go to the documentation of this file.
1 ! *********************
2  MODULE thermal_khione
3 ! *********************
4 !
5 !***********************************************************************
6 ! KHIONE V7P3
7 !***********************************************************************
8 !
9 !brief Module containing all subroutines to deal with the physics
10 !+ of thermal exchanges, at a node level.
11 !+ To be joined up with THERMAL_WAQTEL.
12 !
13 !history F. SOUILLE (EDF)
14 !+ 30/09/2019
15 !+ V8P0
16 !+ Added coefficients for full heat budget
17 !+ Added freezing temperature
18 !
19 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
20 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
21 !
23  IMPLICIT NONE
24 !
25  PRIVATE
27  & daynum
28 !
29 !=======================================================================
30 !
31 ! 1) THERMAL FLUXES
32 !
33  CONTAINS
34 !
35 !=======================================================================
36 !
37 ! *************************
38  SUBROUTINE thermal_fluxes
39 ! *************************
40 !
41  &(tair,twat,tfrz,tdew,cc,visb,wind,pluie,sumph,phcl,phri,phps,
42  & phib,phie,phih,phip,anfem,dt,at,mardat,martim,lambd0,cice)
43 !
44 !***********************************************************************
45 ! RICE-2D V7P2 11/11/2016
46 !***********************************************************************
47 !
48 !brief
49 !
50 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
51 !| ANFEM |-->| CONCENTRATION OF SURFACE ICE PARTICLES
52 !| AT |-->| CURRENT TIME
53 !| CC |-->| CLOUD COVER
54 !| DN |-->| CURRENT TIME
55 !| DT |-->| TIME STEP
56 !| LAMBD0 |-->| LATITUDE OF ORIGIN POINT (KEYWORD, IN DEGREES)
57 !| MARDAT |-->| DATE (YEAR, MONTH,DAY)
58 !| MARTIM |-->| TIME (HOUR, MINUTE,SECOND)
59 !| PHCL |<->| SOLAR RAD (FLUX) REACHING SURFACE, UNDER CLEAR SKY
60 !| PHIB |<->| EFFECTIVE BACK RADIATION (OPEN WATER OR ICE)
61 !| PHIE |<->| EVAPORATIVE HEAT TRANSFER
62 !| PHIH |<->| CONVECTIVE HEAT TRANSFER
63 !| PHIP |<->| HEAT TRANSFER DUE TO PRECIPITATION
64 !| PHPS |<->| NET SOLAR RADIATION (FLUX) AFTER REFLEXION
65 !| PHRI |<->| SOLAR RAD (FLUX) REACHING SURFACE, UNDER CLOUDY SKY
66 !| PLUIE |-->| RAIN
67 !| SUMPH |<->| NET SUM OF ALL THERMAL FLUXES
68 !| T1 |-->| STARTING HOUR FOR SOLAR RADIATION CALCULATION (HRS)
69 !| T2 |-->| ENDING HOUR FOR SOLAR RADIATION CALCULATION (HRS)
70 !| TAIR |-->| AIR TEMPERATURE
71 !| TDEW |-->| DEWPOINT TEMPERATURE
72 !| TFRZ |-->| FREEZING TEMPERATURE
73 !| TWAT |-->| WATER TEMPERATURE
74 !| VISB |-->| VISIBILITY
75 !| WIND |-->| WIND SPEED EFFECT ON ICE
76 !| CICE |-->| EXCHANGES WITH OPEN WATER (0) OR ICE COVER (1)
77 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
78 !
80  & lin_watair,cst_watair,lin_iceair,cst_iceair,
81  & coef_phib,coef_phie,coef_phih,coef_phip,sgma,
82  & cp_eau, atmoexch
83  USE meteo_telemac, ONLY: windz
84 !
85  IMPLICIT NONE
86 !
87 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
88 !
89  INTEGER, INTENT(IN) :: MARDAT(3),MARTIM(3),CICE
90 !
91  DOUBLE PRECISION, INTENT(IN) :: TWAT,ANFEM,TFRZ
92  DOUBLE PRECISION, INTENT(IN) :: TAIR,TDEW,CC,VISB,WIND,PLUIE
93  DOUBLE PRECISION, INTENT(IN) :: DT,AT
94  DOUBLE PRECISION, INTENT(INOUT) :: SUMPH,PHCL,PHRI
95  DOUBLE PRECISION, INTENT(INOUT) :: PHPS,PHIB,PHIE,PHIH,PHIP
96  DOUBLE PRECISION, INTENT(IN) :: LAMBD0
97 !
98 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
99 !
100 ! DAY NUMBER, ORBITAL CORRECTION
101  INTEGER IYEAR,IMONTH,IDAY,IHOUR,IMIN,ISEC
102  INTEGER NDLT,I
103 !
104  DOUBLE PRECISION DN, DAY,DAYREEL,NDAYS,DTHR
105 !
106  DOUBLE PRECISION PHBA,PHBC,PHBR,PHBW
107  DOUBLE PRECISION ES1,EA1,EPINA,AKN,VA,ASV
108  DOUBLE PRECISION PCL,PRI,PPS,TAK,TSK,TDK
109  DOUBLE PRECISION T1,T2,T10,T20
110 !
111 !-----------------------------------------------------------------------
112 !
113  iyear = mardat(1)
114  imonth = mardat(2)
115  iday = mardat(3)
116  ihour = martim(1)
117  imin = martim(2)
118  isec = martim(3)
119 !
120 !-----------------------------------------------------------------------
121 !
122 ! PHCL: SOLAR RAD (FLUX) REACHING SURFACE, UNDER CLEAR SKY
123 ! PHRI: SOLAR RAD (FLUX) REACHING SURFACE, UNDER CLOUDY SKY
124 ! PHPS: NET SOLAR RAD (FLUX) AFTER REFLECTION
125 ! PHIB: BACK RADIATION HEAT TRANSFER
126 ! PHIE: EVAPORATIVE HEAT TRANSFER
127 ! PHIH: CONVECTIVE HEAT TRANSFER
128 ! PHIP: HEAT TRANSFER DUE TO PRECIPITATION
129  sumph = 0.d0
130 !
131 !-----------------------------------------------------------------------
132 !
133 ! DAY NUMBER, ORBITAL CORRECTION
134  day = daynum(iyear,imonth,iday,ihour,imin,isec)
135  & + floor(at/86400.d0) ! DAY FROM JAN.1
136  ndays = 365.d0 + REAL(leap(iyear))
137  dayreel = modulo(day, ndays) ! DAY NUMBER IN DATE
138 !
139  dn = day ! INT # OF DAYS FROM JAN 1 = DAYREEL
140  dthr = dt / 3600.d0
141 
142  phcl = 0.d0
143  phri = 0.d0
144  phps = 0.d0
145 
146  phib = 0.d0
147  phie = 0.d0
148  phih = 0.d0
149  phip = 0.d0
150 
151 ! MSEC = NDAYS * 86400 + IHOUR * 3600 + IMIN * 60 + ISEC
152 !
153 ! T10 = FRACTION # OF HOURS FROM 0:00 @ CURRENT TIME (TSUM1)
154 ! T20 = FRACTION # OF HOURS FROM 0:00 AFTER DTHR
155 !
156  t10 = ihour + modulo(at,86400.d0)/3600.d0
157  t20 = t10 + dthr
158  t1 = t10
159 !
160  IF (dthr.GE.1.d0) THEN ! IF TIME INTERVAL > 1.0 HR (RARE)
161  ndlt = int(dthr + 0.001d0)
162  DO i = 1,ndlt
163  t2 = t1 + 1.d0
164  IF (t2.GT.t20) THEN
165  CALL solar(dn,t1,t20,cc,pcl,pri,pps,cice,lambd0)
166  phcl = phcl + pcl ! ADD FRACTION OF HOUR LEFT
167  phri = phri + pri
168  phps = phps + pps
169  EXIT
170  ELSE
171  CALL solar(dn,t1,t2,cc,pcl,pri,pps,cice,lambd0)
172  phcl = phcl + pcl ! ADD HOUR INCREMENTS
173  phri = phri + pri
174  phps = phps + pps
175  ENDIF
176  t1 = t2
177  ENDDO
178  ELSE
179  CALL solar(dn,t10,t20,cc,pcl,pri,pps,cice,lambd0)
180  phcl = phcl + pcl ! IF TIME INTERVAL A FRACTION OF AN HOUR
181  phri = phri + pri
182  phps = phps + pps
183  ENDIF
184 !
185  phps = phps / dthr ! NET SOLAR RAD AT SURFACE, W/M^2
186 !
187 !-----------------------------------------------------------------------
188 !
189 ! LINEAR HEAT TRANSFER FOR AIR-WATER INTERFACE ONLY
190 !
191  IF( atmoexch.EQ.0 ) THEN
192  IF(cice.EQ.1) THEN ! ICE
193  phih = cst_iceair + ( tair-tfrz )*lin_iceair
194  ELSE ! OPEN WATER
195  phih = cst_watair + ( tair-twat )*lin_watair
196  ENDIF
197  sumph = phps + phih
198 !
199 !-----------------------------------------------------------------------
200 !
201 ! FULL BUDGET FOR AIR-WATER INTERFACE
202 !
203  ELSEIF( atmoexch.EQ.1 ) THEN
204 !
205  tak = tair + 273.16d0 ! TAK = AIR TEMPERATURE
206  tdk = tdew + 273.16d0 ! TDK = DEW POINT
207 !
208  IF (cice.EQ.0) THEN
209  tsk = twat + 273.16d0 ! TSK = WATER TEMP
210  ELSE
211  tsk = 273.16d0
212  ENDIF
213 !
214 ! ~~~~~~~~~~~~~~
215 ! BACK RADIATION
216 ! ~~~~~~~~~~~~~~
217  IF (cice.EQ.0) THEN
218 !
219 ! FOR WATER SURFACE, ES1 FOUND USING WATER TEMP (TSK), (4.32)
220  es1 = 7.95357242d10*exp((-18.1972839d0*373.16d0/tsk)+
221  & 5.02808d0*log(373.16d0/tsk)-20242.1852d0*
222  & exp(-26.1205253d0*tsk/373.16d0)+
223  & 58.0691913d0*exp(-8.039282d0*373.16d0/tsk))
224  ea1 = 7.95357242d10*exp((-18.1972839d0*373.16d0/tdk)+
225  & 5.02808d0*log(373.16d0/tdk)-20242.1852d0*
226  & exp(-26.1205253d0*tdk/373.16d0)+
227  & 58.0691913d0*exp(-8.039282d0*373.16d0/tdk))
228  ELSE
229 !
230 ! FOR ICE SURFACE, ES1 FOUND USING AIR TEMP (TAK), (4.33)
231 ! USUALLY CONSIDERED WHEN ANFEM(I) > 0.5,
232 ! MORE ICE THAN WATER AT SURFACE
233  es1 = 5.75185606d10*exp((-20.947031d0*273.16d0/tak)-
234  & 3.56654d0*log(273.16d0/tak)-2.01889049d0/273.16d0*tak)
235  ea1 = 5.75185606d10*exp((-20.947031d0*273.16d0/tdk)-
236  & 3.56654d0*log(273.16d0/tdk)-2.01889049d0/273.16d0*tdk)
237  ENDIF
238 !
239 ! EMISSIVITY OF ATMOSPHERE = PHBA (INCL. EFFECTS OF CLOUDS)
240  IF (tair.LT.0.d0) THEN
241  ! (4.31) SATTERLUND (1979)
242  epina = 1.08d0 * (1.d0 - exp(-ea1 ** (tak / 2016.d0)))
243  ELSE
244  ! TK 03-2010 IDSO-JACKSON(1969)
245  epina = 1.d0-0.261d0*exp(-0.000777d0*(273.16d0-tak)**2)
246  ENDIF
247  ! (4.30)
248  phbc = epina * sgma * tak ** 4
249  ! = ATMOSPHERIC RADIATION UNDER CLOUDY SKY (TELEMAC 2D)
250  phba = phbc * (1.d0 + 0.0017d0 * cc **2)
251 !
252 ! EMISSIVITY OF RIVER SURFACE = PHBW
253  ! (4.29), TSK = WATER TEMP (ALWAYS)
254  phbw = 0.97d0 * sgma * tsk ** 4
255 
256 ! REFLECTED LONG WAVE RADIATION = PHBR
257  ! (4.36)
258  phbr = 0.03d0 * phba
259 !
260 ! EFFECTIVE BACK RADIATION
261  phib = phba - phbr - phbw
262  phib = coef_phib*phib
263 !
264 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
265 ! EVAPORATIVE HEAT FLUX AND CONVECTIVE HEAT FLUX
266 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
267  IF (cice.EQ.0) THEN
268  akn = 8.d0 + 0.35d0 * (twat - tair)
269  va = (2.d0/windz) ** 0.15d0 * wind
270 ! WINDZ = HEIGHT WIND VELOCITY IS MEASURED (*.PAR)
271 !
272 ! EVAPORATION
273  phie = -(1.56d0*akn + 6.08d0 * va) * (es1-ea1)*4.1855d0/8.64d0
274  phie = coef_phie*phie
275  ELSE
276  akn = 8.d0 + 0.35d0 * (0.d0 - tair)
277  va = (2.d0/windz) ** 0.15d0 * wind
278  ENDIF
279 !
280 ! CONDUCTIVE HEAT TRANSFER
281  phih = (akn + 3.9d0*va) * (tak-tsk)*4.1855d0/8.64d0
282  phih = coef_phih*phih
283 !
284 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
285 ! HEAT TRANSFER DUE TO PRECIPITATION
286 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
287  IF (visb.GT.0.d0.AND.visb.LT.1.d0) THEN
288 ! MIN VALUE = 1.0 KM, PHIP = 380.28 W/M^2
289  asv = 78.5d0 / 86400.d0
290  ELSE
291 ! AT VISB = 10 KM, PHIP = 1.604 W/M^2
292  asv = 78.5d0 / 86400.d0 * visb ** (-2.375d0)
293  ENDIF
294 !
295 ! NO SNOW
296  IF(visb.LT.1.d-5) THEN ! NO SNOW
297  asv = 0.d0
298  ENDIF
299 !
300 ! SNOW FALL
301  IF(cice.EQ.0) THEN
302  phip = asv * (cp_ice * (tair - twat) - lh_ice)
303  ELSE
304  phip = 0.d0
305  ENDIF
306 !
307 ! RAIN FALL
308  IF(pluie.GT.0.d0) THEN
309  asv = pluie/3600.d0
310  IF(cice.EQ.0) THEN ! WATER
311  IF(tair.LT.0.d0) THEN
312  phip = - asv*cp_eau*(0.d0 - twat) ! HEAT LOSS
313  ELSE
314  phip = - asv*cp_eau*(tair - twat) ! HEAT GAIN
315  ENDIF
316  ELSE ! ICE
317  IF(tair.GT.0) THEN
318  phip = - asv*cp_eau*(tair - 0.d0) ! HEAT GAIN
319  ELSE
320  phip = 0.d0
321  ENDIF
322  ENDIF
323  ENDIF
324  phip = coef_phip*phip
325 !
326 ! ~~~~~~~~~~~~~~~~~~~~
327 ! SUMMATION AND OUTPUT
328 ! ~~~~~~~~~~~~~~~~~~~~
329  sumph = phps + phib + phie + phih + phip
330 !
331  ENDIF
332 !
333 !-----------------------------------------------------------------------
334 !
335  RETURN
336  END SUBROUTINE thermal_fluxes
337 !
338 !=======================================================================
339 !
340 ! ************************
341  SUBROUTINE icover_growth
342 ! ************************
343 !
344  &(twat,tmelt,sumph,anfem,thifems,hwi,dt)
345 !
346 !***********************************************************************
347 ! RICE-2D V7P3 11/11/2016
348 !***********************************************************************
349 !
350 !brief
351 !
352 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
353 !| DT |-->| TIME STEP
354 !| HWI |-->| WATERICE HEAT COEFFICIENT
355 !| SUMPH |-->| NET HEAT EXCHANGE FLUX WITH THE ATMOSPHERE
356 !| ANFEM |<->| CONCENTRATION OF SURFACE ICE PARTICLES
357 !| THIFEMS |<->| SOLID ICE THICKNESS
358 !| TMELT |-->| FREEZING POINT OF WATER
359 !| TWAT |-->| WATER TEMPERATURE
360 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
361 !
363 !
364  IMPLICIT NONE
365 !
366 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
367 !
368  DOUBLE PRECISION, INTENT(IN) :: HWI
369  DOUBLE PRECISION, INTENT(IN) :: TWAT,TMELT,SUMPH
370  DOUBLE PRECISION, INTENT(IN) :: DT
371  DOUBLE PRECISION, INTENT(INOUT) :: ANFEM,THIFEMS
372 !
373 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
374 !
375 !
376  DOUBLE PRECISION AUX,DH
377 !
378 !-----------------------------------------------------------------------
379 !
380  aux = dt/(rho_ice*lh_ice)
381  dh = - aux*sumph/(thifems*lin_iceair/tc_bi + 1.d0)
382  & - aux*hwi*(twat-tmelt)
383 !
384  IF( dh.GT.(-thifems) ) THEN
385  thifems = thifems + dh
386  ELSE
387  thifems = 0.d0
388  ENDIF
389 !
390 !-----------------------------------------------------------------------
391 !
392 ! NOT ALLOWING NEGATIVE ICE COVER THICKNESS
393 !
394  IF( thifems.LE.0.d0 ) thifems = 0.d0
395  IF( thifems.LE.0.d0 ) anfem = 0.d0
396 !
397 !-----------------------------------------------------------------------
398 !
399  RETURN
400  END SUBROUTINE icover_growth
401 !
402 !=======================================================================
403 !
404 ! ****************
405  SUBROUTINE solar
406 ! ****************
407  &(dn,t1,t2,cc,phcl,phri,phps,cice,lambd0)
408 !
409 !***********************************************************************
410 ! KHIONE V7P3
411 !***********************************************************************
412 !
413 !brief Same as SOLRAD within EXCHANGE_WITH_ATMOSPHERE (WAQTEL)
414 !+ TODO: Combined with SOLAR
415 !
416 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
417 !| CC |-->| CLOUD COVER, IN THENTHS, 0-10 (/!\ NOT IN OCTAS)
418 !| CICE |-->| ICE CONDITION, OPEN WATER : 0 ; ICE : 1
419 !| DN |-->| CURRENT TIME (THE NUMBER OF DAYS FROM JAN. 1)
420 !| T1 |-->| STARTING HOUR FOR SOLAR RADIATION CALCULATION (HRS)
421 !| T2 |-->| ENDING HOUR FOR SOLAR RADIATION CALCULATION (HRS)
422 !| PHCL |<->| SOLAR RAD (FLUX) REACHING SURFACE, UNDER CLEAR SKY
423 !| PHPS |<->| NET SOLAR RADIATION (FLUX) AFTER REFLEXION
424 !| PHRI |<->| SOLAR RAD (FLUX) REACHING SURFACE, UNDER CLOUDY SKY
425 !| LAMBD0 |-->| LATITUDE OF ORIGIN POINT (NORTH +, SOUTH -, DEGREES)
426 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
427 !
429  USE meteo_telemac, ONLY: modelz,alphsd,alphrd
430 !
431  IMPLICIT NONE
432 !
433 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
434 !
435  INTEGER , INTENT(IN) :: CICE
436  DOUBLE PRECISION, INTENT(IN) :: CC
437  DOUBLE PRECISION, INTENT(IN) :: DN,T1,T2
438  DOUBLE PRECISION, INTENT(INOUT) :: PHCL,PHRI,PHPS
439  DOUBLE PRECISION, INTENT(IN) :: LAMBD0
440 !
441 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
442 !
443  DOUBLE PRECISION :: PI
444  DOUBLE PRECISION :: AA, ALPHA,ALPHR,ALPHS,ALHA
445  DOUBLE PRECISION R,BB,AMO,AM,CR,DELTA,EO,ET,HRIS,HSS,
446  & papo,phiso,phrr,phs,rhs1,rhs2,sia,alr,hs1,hs2,hset
447 !
448 !-----------------------------------------------------------------------
449 !
450 ! DELTA = SOLAR DECLINATION OF THE SUN;
451 ! HSS = LOCAL HOUR ANGLE OF THE SUNSET, RADIANS;
452 ! HSR = LOCAL HOUR ANGLE OF THE SUNRISE, RADIANS;
453 ! ALMST = LOCAL MEAN SOLAR TIME IN HRS;
454 ! ALST = STANDARD TIME OF THE TIME ZONE, IN HRS, COUNTED FROM MIDNIGHT, 0-24.0;
455 ! ALSM = LONGITUDE OF THE STANDARD MERIDEAN, IN DEGREE; (*.PAR FILE)
456 ! ALLM = LONGITUDE OF THE LOCAL MERIDEAN, IN DEGREES; (*.PAR FILE)
457 ! ETA = -1, FOR WEST; + FOR EAST; (*.PAR FILE)
458 ! MODELZ = ELEVATION ABOVE SEA LEVEL, M. (*.PAR FILE)
459 ! ALPHSD = SUN EXIT ANGLE, 180 FOR HORIZONTAL; (*.PAR FILE)
460 ! ALPHRD = SUN EMISION ANGLE, 0 DEGREE FOR HORIZONTAL; (*.PAR FILE)
461 ! NM = NUMBER OF MONTH
462 ! ND = DAY NUMBER IN THE DATE
463 !
464 !-----------------------------------------------------------------------
465 !
466  phcl = 0.d0
467  phri = 0.d0
468  phps = 0.d0
469 !
470  pi = 4.d0*atan(1.d0)
471 
472 ! LOCAL GEOGRAPHIC LATITUDE, CONVERT TO RADIANS
473  phs = lambd0 * pi / 180.d0 ! LATITUDE
474  alphs = alphsd * pi / 180.d0 ! EXIT ANGLE, RADIANS
475  alphr = alphrd * pi / 180.d0 ! EMISSION ANGLE, RADIANS
476 
477 ! COOPER 1969 SOLAR DECLINATION, IN RADIANS (4.6)
478  delta = 23.45d0*pi/180.d0*sin(360.d0*(284.d0+dn)/365.d0*pi/180.d0)
479 
480 ! DIFFERENCE BETWEEN TRUE SOLAR TIME AND MEAN SOLAR TIME
481  ! (4.3)
482  r = 2.d0 * pi * (dn - 1.d0) / 365.d0
483 ! DUFFIE AND BECKMAN 1959, ECCENTRICITY CORRECTION FACTION OF THE EARTH ORBIT
484  ! (4.5)
485  eo = 1.d0 + 0.033d0 * cos(2.d0 * pi / 365.d0 * dn)
486 ! EQUATION OF TIME, IN HRS (4.2)
487  et = 3.8197d0 * (0.000075d0 + 0.001868d0 * cos(r)
488  & -0.032077d0*sin(r)-0.014615d0*cos(2.d0*r)-0.04089d0*sin(2.d0*r))
489 
490 ! HOUR ANGLE AT SUNRISE, RADIANS (4.9)
491  hss = acos( -tan(phs) * tan(delta) )
492  IF (hss.LT.0.d0) THEN
493  hss=-hss
494  ENDIF
495 ! SUN SET, HRS (4.11)
496  hset = 12.d0 + hss * 12.d0 / pi - (pi - alphs) / pi * 12.d0
497 ! SUN RISE, HRS (4.10)
498  hris = 12.d0 - hss * 12.d0 / pi + alphr / pi * 12.d0
499 ! CALCULATE HOUR ANGLE, TIME CORRECTION (4.4)
500  hs1 = t1 - etadir / 15.d0 * (alsm - allm) + et ! HOURS
501  rhs1 = (12.d0 - hs1) * pi / 12.d0 ! RADIANS
502  hs2 = t2 - etadir / 15.d0 * (alsm - allm) + et ! HOURS
503  rhs2 = (12.d0 - hs2) * pi / 12.d0 ! RADIANS
504 ! RHS1/RHS2 IS # OF RADIANS OFF OF NOON, + BEFORE, - AFTER
505 
506 ! NET SOLAR RADIATION, PHPS
507  phcl = 0.d0 ! CLEAR SKY SOLAR RADIATION
508  phri = 0.d0 ! INCL. CLOUD EFFECTS, IF ANY
509  phps = 0.d0 ! NET SOLAR, AFTER REFLECTION AT EARTH SURFACE
510 
511 ! NO SUN, BEFORE SUNRISE OR AFTER SUNSET
512  IF (hs1.GT.hset.AND.hs2.GT.hset) THEN
513  RETURN
514  ENDIF
515  IF (hs1.LT.hris.AND.hs2.LT.hris) THEN
516  RETURN
517  ENDIF
518 
519  IF (hs1.LT.hris .AND. hs2.GT.hris) THEN
520  rhs1=(12.d0-hris)*pi/12.d0
521  hs1=hris
522  ENDIF
523  IF (hs2.GT.hset.AND.hs1.LT.hset) THEN
524  rhs2=(12.0-hset)*pi/12.d0
525  hs2=hset
526  ENDIF
527 
528 ! AVERAGE SOLAR TIME ANGLE
529  alha = (rhs1 + rhs2) / 2.d0 ! RADIANS
530 ! TOTAL SOLAR RADIATION T1 TO T2, PER UNIT AREA, (4.13)-RADIANS, (4.14)-HOURS
531  phiso = 12.d0/pi*sio*eo*((rhs1-rhs2)*sin(delta)*sin(phs)+
532  & (sin(rhs1)-sin(rhs2))*cos(delta)*cos(phs))
533 
534  sia = sin(delta) * sin(phs) + cos(delta) * cos(phs) * cos(alha)
535 ! LINE ABOVE -- (4.7)
536  alpha = asin( sia ) * 180.d0 / pi ! CONVERT TO DEGREES
537  ! (4.17)
538  amo = 1.d0 / (sia + 0.15d0 * (alpha + 3.885d0) ** (-1.253d0))
539  ! (4.18)
540  papo = exp(-0.0001184d0 * modelz)
541  ! (4.16)
542  am = amo * papo
543  ! (4.15)
544  am = 0.99d0 - 0.17d0 * am
545  ! ENERGY FLUX, REACHING GROUND UNDER CLEAR SKY
546  phcl = am * phiso
547 
548  IF (phcl.LT.0.d0) phcl=0.d0
549  ! (4.19)
550  phri = phcl * ( 1.d0 - 0.0065d0 * cc **2 )
551 ! SOLAR RADIATION REACHING THE EARTH UNDER CLOUDY SKIES (ABOVE)
552 
553  IF (cice.EQ.0) THEN ! IF ALBE = 0, "OPEN WATER" CONDITIONS
554  IF (phcl.EQ.0.d0) THEN
555  phrr=0.d0
556  ELSE ! ESTIMATE REFLECTIVITY OF OPEN WATER SURFACE
557  ! (4.25) - CLOUDINESS RATIO
558  cr = max((1.d0 - phri / phcl),0.d0)
559  aa = 2.2d0 +cr**0.7d0/4.d0 - (cr**0.7d0 - 0.4d0)**2 / 0.16d0
560  bb = -1.02d0+cr**0.7d0/16.d0 + (cr**0.7d0 - 0.4d0)**2 / 0.64d0
561  ! ESTIMATE OPEN WATER ALBEDO (4.22-4.24)
562  alr = aa * alpha ** bb
563  ! AMOUNT OF SOLAR RADIATION REFLECTED
564  phrr = alr * phri
565  ENDIF
566 
567 ! USE THE ALBEDO SPECIFIED TO INCLUDE REFLECTIVE EFFECTS OF ICE PRESENT
568  ELSE
569  ! (4.20)
570  phrr = albe * phri
571  IF (phcl.EQ.0.) phrr=0.d0
572  ENDIF
573  ! NET SOLAR RADIATION BETWEEN T1-T2, (4.20)
574  phps = phri - phrr
575 
576  RETURN
577  END SUBROUTINE solar
578 !
579 !=======================================================================
580 !
581 ! ********************************************
582  DOUBLE PRECISION FUNCTION waterice_heat_coef
583 ! ********************************************
584  &(h,u,twat,tfrz)
585 !
586 !***********************************************************************
587 ! KHIONE V7P3
588 !***********************************************************************
589 !
590 !brief Computes the heat exchange coefficient with the ice cover
591 !
592 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
593 !| H |-->| PROPAGATION DEPTH
594 !| U |-->| VELOCITY MAGNITUDE
595 !| TWAT |-->| WATER TEMPERATURE
596 !| TFRZ |-->| FREEZING POINT FOR WATER
597 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
598 !
600 !
601  IMPLICIT NONE
602 !
603 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
604 !
605  DOUBLE PRECISION, INTENT(IN) :: H,U,TWAT,TFRZ
606 !
607 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
608 !
609  DOUBLE PRECISION RE, DH
610 !
611 !-----------------------------------------------------------------------
612 !
613 ! REYNOLDS NUMBER
614  dh = 4.d0*h
615  re = u * dh / xnu
616 !
617 ! DRY BANKS
618  IF(dh.LE.0.001d0.AND.dh.GT.0.d0) THEN
619  waterice_heat_coef = 0.d0 !1394.D0
620  RETURN
621  ENDIF
622 !
623 ! LAMINAR FLOW
624  IF(re.LT.2200.d0) THEN
626 !
627 ! TURBULENT FLOW
628  ELSEIF( twat.GT.tfrz ) THEN ! WATER TEMP.>TFRZ
629  waterice_heat_coef = cwi1*u**0.8d0/dh**0.2d0 ! CWI1 = 1448
630  ELSE ! WATER TEMP.<TFRZ
631  waterice_heat_coef = ciw1*u**0.8d0/dh**0.2d0 ! CIW1 = 1118
632  ENDIF
633 !
634  RETURN
635 !
636 !-----------------------------------------------------------------------
637 !
638  END FUNCTION waterice_heat_coef
639 !
640 ! *********************
641  INTEGER FUNCTION leap
642 ! *********************
643 !
644  &(iyear)
645 !
646 !***********************************************************************
647 ! TELEMAC-3D V6P2 27/06/2012
648 !***********************************************************************
649 !
650 !brief DETERMINES WHETHER IYEAR IS A LEAP YEAR
651 !+ DESCRIPTION - RETURNS 1 IF IYEAR IS A LEAP YEAR, 0 OTHERWISE
652 !+
653 !
654 !history C. GUILBAUD (SOGREAH)
655 !+ JUNE 2001
656 !+ V6P0?
657 !+
658 !
659 !history C.-T. PHAM (EDF-LNHE)
660 !+ 27/06/2012
661 !+ V6P2
662 !+ Introduction into EXCHANGE_WITH_ATMOSPHERE module + INTENT
663 !+
664 !
665 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
666 !| IYEAR |-->| INDEX OF YEAR
667 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
668 !
669  IMPLICIT NONE
670 !
671 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
672 !
673  INTEGER, INTENT(IN) :: IYEAR
674 !
675 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
676 !
677  IF( mod(iyear,4).EQ.0.AND.
678  & (mod(iyear,100).NE.0.OR.mod(iyear,400).EQ.0)) THEN
679  leap = 1
680  ELSE
681  leap = 0
682  ENDIF
683 !
684 !-----------------------------------------------------------------------
685 !
686  RETURN
687  END FUNCTION leap
688 
689 ! ********************************
690  DOUBLE PRECISION FUNCTION daynum
691 ! ********************************
692 !
693  &(iyear,imonth,iday,ihour,imin,isec)
694 !
695 !***********************************************************************
696 ! TELEMAC-3D V6P2 27/06/2012
697 !***********************************************************************
698 !
699 !brief RETURNS DAY NUMBER OF THE YEAR (FRACTIONAL)
700 !+
701 !
702 !history C. GUILBAUD (SOGREAH)
703 !+ JUNE 2001
704 !+ V6P0?
705 !+
706 !
707 !history C.-T. PHAM (EDF-LNHE)
708 !+ 27/06/2012
709 !+ V6P2
710 !+ Introduction into EXCHANGE_WITH_ATMOSPHERE module
711 !+ Change for type of result (double precision, not integer)
712 !+ + REAL conversion + addition of seconds ISEC
713 !+
714 !
715 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
716 !| IDAY |-->| INDEX OF DAY
717 !| IHOUR |-->| INDEX OF HOUR
718 !| IMIN |-->| INDEX OF MINUTE
719 !| IMONTH |-->| INDEX OF MONTH
720 !| ISEC |-->| INDEX OF SECOND
721 !| IYEAR |-->| INDEX OF YEAR
722 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
723 !
724  IMPLICIT NONE
725 !
726 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
727 !
728  INTEGER, INTENT(IN) :: IYEAR,IMONTH,IDAY,IHOUR,IMIN,ISEC
729 !
730 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
731 !
732 !
733 ! INTEGER LEAP
734 ! EXTERNAL LEAP
735 !
736  INTEGER MONTH(12)
737  parameter( month=(/0,31,59,90,120,151,181,212,243,273,304,334/) )
738 !
739 !-----------------------------------------------------------------------
740 !
741  daynum = REAL(month(imonth)+iday)
742  & + REAL(ihour)/24.D0+REAL(imin)/1440.D0+REAL(isec)/86400.D0
743  IF(imonth.GT.2) daynum = daynum + REAL(leap(iyear))
744 !
745 !-----------------------------------------------------------------------
746 !
747  RETURN
748  END FUNCTION daynum
749 
750 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
751 !
752  END MODULE thermal_khione
subroutine solar(DN, T1, T2, CC, PHCL, PHRI, PHPS, CICE, LAMBD0)
double precision function, public daynum(IYEAR, IMONTH, IDAY, IHOUR, IMIN, ISEC)
double precision etadir
integer function, public leap(IYEAR)
double precision cp_ice
double precision lh_ice
subroutine, public thermal_fluxes(TAIR, TWAT, TFRZ, TDEW, CC, VISB, WIND, PLUIE, SUMPH, PHCL, PHRI, PHPS, PHIB, PHIE, PHIH, PHIP, ANFEM, DT, AT, MARDAT, MARTIM, LAMBD0, CICE)
double precision tc_wt
double precision lin_iceair
double precision, public alphsd
double precision function, public waterice_heat_coef(H, U, TWAT, TFRZ)
double precision, public alphrd
double precision, public modelz
double precision rho_ice
subroutine, public icover_growth(TWAT, TMELT, SUMPH, ANFEM, THIFEMS, HWI, DT)
double precision, public windz
double precision tc_bi