source: palm/trunk/SOURCE/biometeorology_mod.f90 @ 3361

Last change on this file since 3361 was 3361, checked in by knoop, 3 years ago

Introduced global constant rd_d_rv=0.622

File size: 27.7 KB
Line 
1!> @file biometeorology.f90
2!------------------------------------------------------------------------------!
3! This file is part of PALM-4U.
4!
5! PALM-4U is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM-4U is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 2018, Institute of Computer Science,Academy of Sciences, Prague
18!
19! Calculation of PET:
20! Copyright 2016, Deutscher Wetterdienst (DWD) /
21! German Meteorological Service (DWD)
22!------------------------------------------------------------------------------!
23!
24! Current revisions:
25! -----------------
26!
27!
28! Former revisions:
29! -----------------
30! $Id$
31! Initial revision
32!
33!
34!
35! Authors:
36! --------
37! @author Jaroslav Resler <resler@cs.cas.cz>
38! @author Dominik Froehlich <dominik.froehlich@dwd.de>
39!
40!------------------------------------------------------------------------------!
41
42MODULE biometeorology_mod
43!
44!-- Load required variables from existing modules
45    USE arrays_3d,                                                             &
46        ONLY:  hyp, p, pt, q, u, v, w
47
48    USE basic_constants_and_equations_mod,                                     &
49        ONLY:  rd_d_rv
50
51    USE kinds  !< to set precision of INTEGER and REAL arrays according to PALM
52
53    USE radiation_model_mod,  &
54        ONLY: ix, iy, iz, id, mrt_nlevels, mrt_include_sw,                     &
55              mrtinsw, mrtinlw, mrtbl, nmrtbl, radiation
56
57
58    IMPLICIT NONE
59
60    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  bio_mrt              !< biometeorology mean radiant temperature for each MRT box
61    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  bio_mrt_av           !< time average mean radiant temperature for each MRT box
62    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  bio_pet              !< PhysiologiCALLy Equivalent Temperature (PET) for each MRT box
63    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  bio_pet_av           !< time average PhysiologiCALLy Equivalent Temperature (PET) for each MRT box
64
65
66    REAL(wp), PARAMETER :: sigma_sb       = 5.67037321E-8_wp,       & !< Stefan-Boltzmann constant
67                           t_zero = -273.15_wp,                     & !< temperature 0K in Celsius
68                           human_absorb = 0.7_wp,                   & !< SW absorbtivity of human body
69                           human_emiss = 0.97_wp                      !< emissivity of human body (Lindberg 2008)
70
71
72    INTERFACE biometeorology_3d_data_averaging
73       MODULE PROCEDURE biometeorology_3d_data_averaging
74    END INTERFACE biometeorology_3d_data_averaging
75
76    INTERFACE biometeorology_check_data_output
77       MODULE PROCEDURE biometeorology_check_data_output
78    END INTERFACE biometeorology_check_data_output
79
80    INTERFACE biometeorology_data_output_3d
81       MODULE PROCEDURE biometeorology_data_output_3d
82    END INTERFACE biometeorology_data_output_3d
83
84    INTERFACE biometeorology_define_netcdf_grid
85       MODULE PROCEDURE biometeorology_define_netcdf_grid
86    END INTERFACE biometeorology_define_netcdf_grid
87
88
89    SAVE
90
91    PRIVATE
92
93!
94!-- Public functions
95    PUBLIC biometeorology_3d_data_averaging, biometeorology_check_data_output, &
96           biometeorology_data_output_3d, biometeorology_define_netcdf_grid
97
98!
99!-- Public variables and constants / NEEDS SORTING
100!   PUBLIC
101
102
103 CONTAINS
104
105
106
107
108!------------------------------------------------------------------------------!
109! Description:
110! ------------
111!> Check data output for biometeorology model
112!------------------------------------------------------------------------------!
113    SUBROUTINE biometeorology_check_data_output( var, unit, i, ilen, k )
114 
115 
116       USE control_parameters,                                                 &
117           ONLY: data_output, message_string
118
119       IMPLICIT NONE
120
121       CHARACTER (LEN=*) ::  unit     !<
122       CHARACTER (LEN=*) ::  var      !<
123
124       INTEGER(iwp) :: i
125       INTEGER(iwp) :: ilen
126       INTEGER(iwp) :: k
127
128       SELECT CASE ( TRIM( var ) )
129
130          CASE ( 'bio_mrt', 'bio_pet'  )
131             IF ( .NOT.  radiation ) THEN
132                message_string = 'output of "' // TRIM( var ) // '" require'&
133                                 // 's radiation = .TRUE.'
134                CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
135             ENDIF
136             IF ( mrt_nlevels == 0 ) THEN
137                message_string = 'output of "' // TRIM( var ) // '" require'&
138                                 // 's mrt_nlevels > 0'
139                CALL message( 'check_parameters', 'PA0510', 1, 2, 0, 6, 0 )
140             ENDIF
141             IF ( TRIM( var ) == 'bio_mrt' ) THEN
142                unit = 'K'
143             ELSE
144                unit = 'W m-2'
145             ENDIF
146
147          CASE DEFAULT
148             unit = 'illegal'
149
150       END SELECT
151
152    END SUBROUTINE biometeorology_check_data_output
153
154
155
156!------------------------------------------------------------------------------!
157!
158! Description:
159! ------------
160!> Subroutine for averaging 3D data
161!------------------------------------------------------------------------------!
162SUBROUTINE biometeorology_3d_data_averaging( mode, variable )
163
164    USE control_parameters
165
166    USE indices
167
168    USE kinds
169
170    IMPLICIT NONE
171
172    CHARACTER (LEN=*) ::  mode    !<
173    CHARACTER (LEN=*) :: variable !<
174
175    INTEGER(iwp) ::  i !<
176    INTEGER(iwp) ::  j !<
177    INTEGER(iwp) ::  k !<
178    INTEGER(iwp) ::  l, m !< index of current surface element
179
180    REAL(wp)     ::  mrt, pet
181
182    IF ( mode == 'allocate' )  THEN
183
184       SELECT CASE ( TRIM( variable ) )
185             CASE ( 'bio_mrt' )
186                IF ( .NOT. ALLOCATED( bio_mrt_av ) )  THEN
187                   ALLOCATE( bio_mrt_av(nmrtbl) )
188                ENDIF
189                bio_mrt_av = 0.0_wp
190
191             CASE ( 'bio_pet' )
192                IF ( .NOT. ALLOCATED( bio_pet_av ) )  THEN
193                   ALLOCATE( bio_pet_av(nmrtbl) )
194                ENDIF
195                bio_pet_av = 0.0_wp
196
197          CASE DEFAULT
198             CONTINUE
199
200       END SELECT
201
202    ELSEIF ( mode == 'sum' )  THEN
203
204       SELECT CASE ( TRIM( variable ) )
205
206          CASE ( 'bio_mrt' )
207             IF ( ALLOCATED( bio_mrt_av ) )  THEN
208
209                 IF ( nmrtbl > 0 )  THEN
210                    IF ( mrt_include_sw )  THEN
211                       bio_mrt_av(:) = bio_mrt_av(:) + ((human_absorb*mrtinsw(:) &
212                           + human_emiss*mrtinlw(:)) / (human_emiss*sigma_sb)) ** .25_wp
213                    ELSE
214                       bio_mrt_av(:) = bio_mrt_av(:) + &
215                                       (human_emiss * mrtinlw(:) / sigma_sb) ** .25_wp
216                    ENDIF
217                 ENDIF
218             ENDIF
219
220          CASE ( 'bio_pet' )
221             IF ( ALLOCATED( bio_pet_av ) )  THEN
222                DO  l = 1, nmrtbl
223                   i = mrtbl(ix,l)
224                   j = mrtbl(iy,l)
225                   k = mrtbl(iz,l)
226
227                   IF ( mrt_include_sw )  THEN
228                      mrt = ((human_absorb*mrtinsw(l) + human_emiss*mrtinlw(l)) &
229                             / (human_emiss*sigma_sb)) ** .25_wp
230                   ELSE
231                      mrt = mrt + (human_emiss * mrtinlw(l) / sigma_sb) ** .25_wp
232                   ENDIF
233
234                   CALL calculate_pet_static(                                      &
235                           pt(k,j,i) * (hyp(k) / 100000.0_wp )**0.286_wp + t_zero, &  !< Air temperature (°C)
236                           q(k,j,i) * hyp(k) / ( q(k,j,i) + rd_d_rv ) / 100._wp,   &  !< Vapor pressure (hPa)
237                           SQRT( MAX( ( ( u(k,j,i) + u(k,j,i+1) ) * 0.5_wp )**2 +  &
238                                      ( ( v(k,j,i) + v(k,j+1,i) ) * 0.5_wp )**2 +  &
239                                      ( ( w(k,j,i) + w(k-1,j,i) ) * 0.5_wp )**2,   &
240                                      0.01_wp ) ),                                 &  !< Wind speed (at scalar gridpoint) (m/s)
241                           mrt + t_zero,                                           &  !< Mean radiant temperature (°C)
242                           (hyp(k) + p(k,j,i)) * 0.01_wp,                          &  !< Air pressure (hPa)
243                           pet )                                                      !< output variable of PET
244
245                    bio_pet_av(l) = bio_pet_av(l) + pet
246                ENDDO
247             ENDIF
248
249          CASE DEFAULT
250             CONTINUE
251
252       END SELECT
253
254    ELSEIF ( mode == 'average' )  THEN
255
256       SELECT CASE ( TRIM( variable ) )
257
258          CASE ( 'bio_mrt' )
259             IF ( ALLOCATED( bio_mrt_av ) )  THEN
260                bio_mrt_av(:) = bio_mrt_av(:) / REAL( average_count_3d, KIND=wp )
261             ENDIF
262
263          CASE ( 'bio_pet' )
264             IF ( ALLOCATED( bio_pet_av ) )  THEN
265                bio_pet_av(:) = bio_pet_av(:) / REAL( average_count_3d, KIND=wp )
266             ENDIF
267
268       END SELECT
269
270    ENDIF
271
272END SUBROUTINE biometeorology_3d_data_averaging
273
274
275!------------------------------------------------------------------------------!
276!
277! Description:
278! ------------
279!> Subroutine defining appropriate grid for netcdf variables.
280!> It is called out from subroutine netcdf.
281!------------------------------------------------------------------------------!
282 SUBROUTINE biometeorology_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
283
284    IMPLICIT NONE
285
286    CHARACTER (LEN=*), INTENT(IN)  ::  var         !<
287    LOGICAL, INTENT(OUT)           ::  found       !<
288    CHARACTER (LEN=*), INTENT(OUT) ::  grid_x      !<
289    CHARACTER (LEN=*), INTENT(OUT) ::  grid_y      !<
290    CHARACTER (LEN=*), INTENT(OUT) ::  grid_z      !<
291
292    found  = .TRUE.
293
294
295!
296!-- Check for the grid
297    SELECT CASE ( TRIM( var ) )
298
299       CASE ( 'bio_mrt', 'bio_pet' )
300          grid_x = 'x'
301          grid_y = 'y'
302          grid_z = 'zu'
303
304       CASE DEFAULT
305          found  = .FALSE.
306          grid_x = 'none'
307          grid_y = 'none'
308          grid_z = 'none'
309
310        END SELECT
311
312 END SUBROUTINE biometeorology_define_netcdf_grid
313
314
315!------------------------------------------------------------------------------!
316!
317! Description:
318! ------------
319!> Subroutine defining 3D output variables
320!------------------------------------------------------------------------------!
321 SUBROUTINE biometeorology_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
322 
323
324    USE indices
325
326    USE kinds
327
328
329    IMPLICIT NONE
330
331    CHARACTER (LEN=*) ::  variable !<
332
333    INTEGER(iwp) ::  av          !<
334    INTEGER(iwp) ::  i, j, k, l  !<
335    INTEGER(iwp) ::  nzb_do      !<
336    INTEGER(iwp) ::  nzt_do      !<
337
338    LOGICAL      ::  found       !<
339
340    REAL(wp)     ::  fill_value = -999.0_wp    !< value for the _FillValue attribute
341
342    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !<
343   
344    REAL(wp)     ::  mrt, pet
345
346    found = .TRUE.
347
348
349    SELECT CASE ( TRIM( variable ) )
350
351
352        CASE ( 'bio_mrt' )
353
354            local_pf = REAL( fill_value, KIND = wp )
355            DO  l = 1, nmrtbl
356                i = mrtbl(ix,l)
357                j = mrtbl(iy,l)
358                k = mrtbl(iz,l)
359                IF ( mrt_include_sw )  THEN
360                    mrt = ((human_absorb*mrtinsw(l) + human_emiss*mrtinlw(l)) &
361                           / (human_emiss*sigma_sb)) ** .25_wp
362                ELSE
363                    mrt = (human_emiss*mrtinlw(l) / sigma_sb) ** .25_wp
364                ENDIF
365                local_pf(i,j,k) = mrt
366            ENDDO
367
368        CASE ( 'bio_pet' )
369            local_pf = REAL( fill_value, KIND = wp )
370            IF ( av == 0 )  THEN
371                DO  l = 1, nmrtbl
372                    i = mrtbl(ix,l)
373                    j = mrtbl(iy,l)
374                    k = mrtbl(iz,l)
375
376                    IF ( mrt_include_sw )  THEN
377                        mrt = ((human_absorb*mrtinsw(l) + human_emiss*mrtinlw(l)) &
378                               / (human_emiss*sigma_sb)) ** .25_wp
379                    ELSE
380                        mrt = mrt + (human_emiss*mrtinlw(l) / sigma_sb) ** .25_wp
381                    ENDIF
382
383                    CALL calculate_pet_static(                                      &
384                       pt(k,j,i) * (hyp(k) / 100000.0_wp )**0.286_wp + t_zero, &  !< Air temperature (°C)
385                       q(k,j,i) * hyp(k) / ( q(k,j,i) + rd_d_rv ) / 100.0_wp,  &  !< Vapor pressure (hPa)
386                       SQRT( MAX( ( ( u(k,j,i) + u(k,j,i+1) ) * 0.5_wp )**2 +  &
387                                  ( ( v(k,j,i) + v(k,j+1,i) ) * 0.5_wp )**2 +  &
388                                  ( ( w(k,j,i) + w(k-1,j,i) ) * 0.5_wp )**2,   &
389                                  0.01_wp ) ),                                 &  !< Wind speed (at scalar gridpoint) (m/s)
390                       mrt + t_zero,                                           &  !< Mean radiant temperature (°C)
391                       (hyp(k) + p(k,j,i)) * 0.01_wp,                          &  !< Air pressure (hPa)
392                       pet )                                                      !< output variable of PET
393
394               local_pf(i,j,k) = pet
395
396            ENDDO
397         ELSE
398            IF ( ALLOCATED( bio_pet_av ) ) THEN
399               DO  l = 1, nmrtbl
400                  i = mrtbl(ix,l)
401                  j = mrtbl(iy,l)
402                  k = mrtbl(iz,l)
403                  local_pf(i,j,k) = bio_pet_av(l)
404               ENDDO
405            ENDIF
406         ENDIF
407
408       CASE DEFAULT
409          found = .FALSE.
410
411    END SELECT
412
413
414 END SUBROUTINE biometeorology_data_output_3d
415
416
417
418!------------------------------------------------------------------------------!
419!
420! Description:
421! ------------
422!> PhysiologiCALLy Equivalent Temperature (PET),
423!  stationary (calculated based on MEMI),
424!  Subroutine based on PETBER vers. 1.5.1996 by P. Hoeppe
425!
426!  Input arguments:
427! ----------------
428!  - ta         : Air temperature               (°C)   REAL(wp)
429!  - tmrt       : Mean radiant temperature      (°C)   REAL(wp)
430!  - v          : Wind speed                    (m/s)   REAL(wp)
431!  - vpa        : Vapor pressure                (hPa)   REAL(wp)
432!  - p          : Air pressure                  (hPa)   REAL(wp)
433!
434!  Output arguments:
435! ----------------
436!  - tx         : PET                           (°C)   REAL(wp)
437!------------------------------------------------------------------------------!
438
439 SUBROUTINE calculate_pet_static(                                              &
440!-- Meteorological input
441  ta, vpa, v, tmrt, p,                                                         &
442!-- Output variables
443  tx )  !,                                                                     &
444!-- Configure sample person (optional)
445!  age, mbody, ht, work, eta, icl, fcl, pos, sex )
446
447   IMPLICIT NONE
448
449   REAL(wp), INTENT( IN ) :: ta, tmrt, v, vpa, p
450
451   REAL(wp), INTENT ( OUT ) :: tx
452
453!  REAL(wp), INTENT ( in ), optional :: age, mbody, ht, work, eta, icl, fcl
454   REAL(wp) :: age, mbody, ht, work, eta, icl, fcl
455!  INTEGER(iwp), INTENT ( in ), optional :: pos, sex
456   INTEGER(iwp) :: pos, sex
457
458   REAL(wp) :: acl, adu, aeff, cair, cb, emcl, emsk, ere, erel, esw,           &
459          evap, facl, feff, food, h, po, rdcl, rdsk, rob, rtv,              &
460          sigm, vpts,                                                          &
461   !  former intent (out)
462   !  - tsk        : Skin temperature              (°C)    real
463   !  - tcl        : Clothing temperature          (°C)    real
464   !  - ws         :                                       real
465   !  - wetsk      : Fraction of wet skin                  real
466          tsk, tcl, wetsk
467
468
469!--      Person data
470!   IF ( .NOT. present( age ) )   age   = 35.
471!   IF ( .NOT. present( mbody ) ) mbody = 75.
472!   IF ( .NOT. present( ht ) )    ht    =  1.75
473!   IF ( .NOT. present( work ) )  work  = 80.
474!   IF ( .NOT. present( eta ) )   eta   =  0.
475!   IF ( .NOT. present( icl ) )   icl   =  0.9
476!   IF ( .NOT. present( fcl ) )   fcl   =  1.15
477!   IF ( .NOT. present( pos ) )   pos   =  1
478!   IF ( .NOT. present( sex ) )   sex   =  1
479   
480   age   = 35.
481   mbody = 75.
482   ht    =  1.75
483   work  = 80.
484   eta   =  0.
485   icl   =  0.9
486   fcl   =  1.15
487   pos   =  1
488   sex   =  1
489
490!-- constants
491   po   = 1013.25  !< preassure at sea level
492   ! p    = 1013.25   !< local preassure (hPa), now defined as input variable
493   rob  =    1.06
494   cb   = 3.64 * 1000.
495   food =    0.
496   emsk =    0.99
497   emcl =    0.95
498   evap = 2.42 * 10. ** 6.
499   sigm = 5.67 * 10. **(-8.)
500
501
502!   write(9,*) 'Call calculate_pet_static(ta=', ta, ', vpa=', vpa, &
503!              ', v=', v, ', tmrt=', tmrt, ', p=', p
504!   flush(9)
505!-- call subfunctions
506   CALL  INKOERP ( age, cair, eta, ere, erel, evap, h, ht, mbody,              &
507         p, rtv, sex, ta, vpa, work )
508
509
510   CALL BERECH ( acl, adu, aeff, cair, cb, emcl, emsk,                 &
511        ere, erel, esw, evap, facl, fcl, feff, food, h, ht, icl,               &
512        mbody, p, po, rdcl, rdsk, rob, sex, sigm,                              &
513        ta, tcl, tmrt, tsk, v, vpa, vpts, wetsk )
514
515
516   CALL PET ( acl, adu, aeff, cair, emcl, emsk, esw, evap,                     &
517        facl, feff, h, p, po, rdcl, rdsk,                                      &
518        rtv, sigm, ta, tcl, tsk, tx, vpts, wetsk )
519
520
521 END SUBROUTINE calculate_pet_static
522
523
524!------------------------------------------------------------------------------!
525! Description:
526! ------------
527!> Calculate internal energy ballance
528!------------------------------------------------------------------------------!
529 SUBROUTINE inkoerp( age, cair, eta, ere, erel, evap, h, ht, mbody,            &
530   &     p, rtv, sex, ta, vpa, work )
531
532
533   REAL(wp) :: age, cair, eta, ere, erel, eres, evap, h, ht, mbody,            &
534   &           met, p, rtv, ta, tex, vpa, vpex, work
535
536   INTEGER(iwp) :: sex
537
538   IF ( sex .EQ. 1 ) THEN
539     met = 3.45 * mbody ** ( 3. / 4. ) * (1. + 0.004 *          &
540           ( 30. - age) + 0.010 * ( ( ht * 100. /                     &
541           ( mbody ** ( 1. / 3. ) ) ) - 43.4 ) )
542   ELSE IF ( sex .EQ. 2 ) THEN
543     met = 3.19 * mbody ** ( 3. / 4. ) * ( 1. + 0.004 *         &
544           ( 30. - age ) + 0.018 * ( ( ht * 100. / ( mbody **         &
545           ( 1. / 3. ) ) ) - 42.1 ) )
546   END IF
547   met = work + met
548
549   h = met * (1. - eta)
550
551!-- SENSIBLE RESPIRATION ENERGY
552
553   cair = 1.01 * 1000.
554   tex  = 0.47 * ta + 21.0
555   rtv  = 1.44 * 10. ** (-6.) * met
556   eres = cair * (ta - tex) * rtv
557
558!-- LATENT RESPIRATION ENERGY
559
560   vpex = 6.11 * 10. ** ( 7.45 * tex / ( 235. + tex ) )
561   erel = 0.623 * evap / p * ( vpa - vpex ) * rtv
562
563!-- SUM OF RESULTS
564
565   ere = eres + erel
566   RETURN
567 END SUBROUTINE inkoerp
568
569
570!------------------------------------------------------------------------------!
571! Description:
572! ------------
573!> Calculate heat gain or loss
574!------------------------------------------------------------------------------!
575 SUBROUTINE BERECH( acl, adu, aeff, cair, cb, emcl, emsk,              &
576        ere, erel, esw, evap, facl, fcl, feff, food, h, ht, icl,               &
577        mbody, p, po, rdcl, rdsk, rob, sex, sigm,                              &
578        ta, tcl, tmrt, tsk, v, vpa, vpts, wetsk )
579
580
581   REAL(wp) :: acl, adu, aeff, c(0:10), cair, cb, cbare,                       &
582        cclo, csum, di, ed, emcl, emsk, enbal,                                 &
583        enbal2, ere, erel, esw, eswdif, eswphy, eswpot,                        &
584        evap, facl, fcl, fec, feff, food, h, hc, he, ht, htcl, icl,            &
585        mbody, p, po, r1, r2, rbare, rcl,                                      &
586        rclo, rclo2, rdcl, rdsk, rob, rsum, sigm, sw, swf, swm,                &
587        ta, tbody, tcl, tcore(1:7), tmrt, tsk, v, vb, vb1, vb2,                &
588        vpa, vpts, wetsk, wd, wr, ws, wsum, xx, y
589
590   INTEGER(iwp) :: count1, count3, j, sex
591   logical :: skipIncreaseCount
592
593   wetsk = 0.
594   adu = 0.203 * mbody ** 0.425 * ht ** 0.725
595
596   hc = 2.67 + ( 6.5 * v ** 0.67)
597   hc   = hc * (p /po) ** 0.55
598   feff = 0.725  !< Posture: 0.725 for stading, 0.696 for sitting
599   facl = (- 2.36 + 173.51 * icl - 100.76 * icl * icl + 19.28                  &
600         * (icl ** 3.)) / 100.
601
602   IF ( facl .GT. 1. )   facl = 1.
603   rcl = ( icl / 6.45) / facl
604   IF ( icl .GE. 2. )  y  = 1.
605
606   IF ( ( icl .GT. 0.6 ) .AND. ( icl .LT. 2. ) )  y = ( ht - 0.2 ) / ht
607   IF ( ( icl .LE. 0.6 ) .AND. ( icl .GT. 0.3 ) ) y = 0.5
608   IF ( ( icl .LE. 0.3 ) .AND. ( icl .GT. 0. ) )  y = 0.1
609
610   r2   = adu * (fcl - 1. + facl) / (2. * 3.14 * ht * y)
611   r1   = facl * adu / (2. * 3.14 * ht * y)
612
613   di   = r2 - r1
614
615!-- SKIN TEMPERATURE
616
617   DO j = 1,7
618
619     tsk    = 34.
620     count1 = 0
621     tcl    = ( ta + tmrt + tsk ) / 3.
622     count3 = 1
623     enbal2 = 0.
624
625     DO
626       acl   = adu * facl + adu * ( fcl - 1. )
627       rclo2 = emcl * sigm * ( ( tcl + 273.2 )**4. -                           &
628         ( tmrt + 273.2 )** 4. ) * feff
629       htcl  = 6.28 * ht * y * di / ( rcl * LOG( r2 / r1 ) * acl )
630       tsk   = 1. / htcl * ( hc * ( tcl - ta ) + rclo2 ) + tcl
631
632!--    RADIATION SALDO
633
634       aeff  = adu * feff
635       rbare = aeff * ( 1. - facl ) * emsk * sigm *                            &
636         ( ( tmrt + 273.2 )** 4. - ( tsk + 273.2 )** 4. )
637       rclo  = feff * acl * emcl * sigm *                                      &
638         ( ( tmrt + 273.2 )** 4. - ( tcl + 273.2 )** 4. )
639       rsum  = rbare + rclo
640
641!--    CONVECTION
642
643       cbare = hc * ( ta - tsk ) * adu * ( 1. - facl )
644       cclo  = hc * ( ta - tcl ) * acl
645       csum  = cbare + cclo
646
647 !--   CORE TEMPERATUR
648
649       c(0)  = h + ere
650       c(1)  = adu * rob * cb
651       c(2)  = 18. - 0.5 * tsk
652       c(3)  = 5.28 * adu * c(2)
653       c(4)  = 0.0208 * c(1)
654       c(5)  = 0.76075 * c(1)
655       c(6)  = c(3) - c(5) - tsk * c(4)
656       c(7)  = - c(0) * c(2) - tsk * c(3) + tsk * c(5)
657       c(8)  = c(6) * c(6) - 4. * c(4) * c(7)
658       c(9)  = 5.28 * adu - c(5) - c(4) * tsk
659       c(10) = c(9) * c(9) - 4. * c(4) *                                       &
660         ( c(5) * tsk - c(0) - 5.28 * adu * tsk )
661
662       IF ( tsk .EQ. 36. ) tsk = 36.01
663       tcore(7) = c(0) / ( 5.28 * adu + c(1) * 6.3 / 3600. ) + tsk
664       tcore(3) = c(0) / ( 5.28 * adu + ( c(1) * 6.3 / 3600. ) /               &
665         ( 1. + 0.5 * ( 34. - tsk ) ) ) + tsk
666       IF ( c(10) .GE. 0.) THEN
667         tcore(6) = ( - c(9) - c(10)**0.5 ) / ( 2. * c(4) )
668         tcore(1) = ( - c(9) + c(10)**0.5 ) / ( 2. * c(4) )
669       END IF
670       ! 22
671       IF ( c(8) .GE. 0. ) THEN
672         tcore(2) = ( - c(6) + ABS( c(8) ) ** 0.5 ) / ( 2. * c(4) )
673         tcore(5) = ( - c(6) - ABS( c(8) ) ** 0.5 ) / ( 2. * c(4) )
674         tcore(4) = c(0) / ( 5.28 * adu + c(1) * 1. / 40. ) + tsk
675       END IF
676       ! 24
677
678!--    TRANSPIRATION
679
680       tbody = 0.1 * tsk + 0.9 * tcore(j)
681       swm   = 304.94 * ( tbody - 36.6 ) * adu / 3600000.
682       vpts  = 6.11 * 10.**( 7.45 * tsk / ( 235. + tsk ) )
683
684       IF ( tbody .LE. 36.6 ) swm = 0.
685       ! swf = 0.7 * swm
686
687       IF ( sex .EQ. 1 ) sw = swm
688       IF ( sex .EQ. 2 ) sw = 0.7 * swm  ! swf
689       eswphy = - sw * evap
690       he     = 0.633 * hc / ( p * cair )
691       fec    = 1. / ( 1. + 0.92 * hc * rcl )
692       eswpot = he * ( vpa - vpts ) * adu * evap * fec
693       wetsk  = eswphy / eswpot
694
695       IF ( wetsk .GT. 1. ) wetsk = 1.
696
697       eswdif = eswphy - eswpot
698
699       IF ( eswdif .LE. 0. ) esw = eswpot
700       IF ( eswdif .GT. 0. ) esw = eswphy
701       IF ( esw  .GT. 0. )   esw = 0.
702
703!--    DIFFUSION
704
705       rdsk = 0.79 * 10. ** 7.
706       rdcl = 0.
707       ed   = evap / ( rdsk + rdcl ) * adu * ( 1. - wetsk ) * ( vpa - vpts )
708
709!--    MAX VB
710
711       vb1 = 34. - tsk
712       vb2 = tcore(j) - 36.6
713
714       IF ( vb2 .LT. 0. ) vb2 = 0.
715       IF ( vb1 .LT. 0. ) vb1 = 0.
716       vb = ( 6.3 + 75. * vb2 ) / ( 1. + 0.5 * vb1 )
717
718!--    ENERGY BALLANCE
719
720       enbal = h + ed + ere + esw + csum + rsum + food
721
722
723!--    CLOTHING TEMPERATURE
724
725       xx = 0.001
726       IF ( count1 .EQ. 0 ) xx = 1.
727       IF ( count1 .EQ. 1 ) xx = 0.1
728       IF ( count1 .EQ. 2 ) xx = 0.01
729       IF ( count1 .EQ. 3 ) xx = 0.001
730
731       IF ( enbal .GT. 0. ) tcl = tcl + xx
732       IF ( enbal .LT. 0. ) tcl = tcl - xx
733
734       skipIncreaseCount = .FALSE.
735       IF ( ( (enbal .LE. 0.) .AND. (enbal2 .GT. 0.) ) .OR.                    &
736          ( ( enbal .GE. 0. ) .AND. ( enbal2 .LT. 0. ) ) ) THEN
737         skipIncreaseCount = .TRUE.
738       ELSE
739         enbal2 = enbal
740         count3 = count3 + 1
741       END IF
742       
743       IF ( ( count3 .GT. 200 ) .OR. skipIncreaseCount ) THEN
744         IF ( count1 .LT. 3 ) THEN
745           count1 = count1 + 1
746           enbal2 = 0.
747         ELSE
748           EXIT
749         END IF
750       END IF
751     END DO
752
753     IF ( count1 .EQ. 3 ) THEN
754       SELECT CASE ( j )
755         CASE ( 2, 5) 
756           IF ( .NOT. ( ( tcore(j) .GE. 36.6 ) .AND.                        &
757             ( tsk .LE. 34.050 ) ) ) CYCLE
758         CASE ( 6, 1 )
759           IF ( c(10) .LT. 0. ) CYCLE
760           IF ( .NOT. ( ( tcore(j) .GE. 36.6 ) .AND.                        &
761             ( tsk .GT. 33.850 ) ) ) CYCLE
762         CASE ( 3 )
763           IF ( .NOT. ( ( tcore(j) .LT. 36.6 ) .AND.                        &
764             ( tsk .LE. 34.000 ) ) ) CYCLE
765         CASE ( 7 )
766           IF ( .NOT. ( ( tcore(j) .LT. 36.6 ) .AND.                        &
767             ( tsk .GT. 34.000 ) ) ) CYCLE
768         CASE default  ! = CASE ( 4 ), does actually nothing
769       END SELECT
770     END IF
771
772     IF ( ( j .NE. 4 ) .AND. ( vb .GE. 91. ) ) CYCLE
773     IF ( ( j .EQ. 4 ) .AND. ( vb .LT. 89. ) ) CYCLE
774     IF ( vb .GT. 90.) vb = 90.
775
776!--  LOSSES BY WATER
777
778     ws = sw * 3600. * 1000.
779     IF ( ws .GT. 2000. ) ws = 2000.
780     wd = ed / evap * 3600. * ( -1000. )
781     wr = erel / evap * 3600. * ( -1000. )
782
783     wsum = ws + wr + wd
784
785     RETURN
786   END DO
787   RETURN
788 END SUBROUTINE Berech
789
790
791 
792!------------------------------------------------------------------------------!
793! Description:
794! ------------
795!> Calculate PET
796!------------------------------------------------------------------------------!
797 SUBROUTINE PET ( acl, adu, aeff, cair, emcl, emsk, esw, evap,                 &
798     facl, feff, h, p, po, rdcl, rdsk, rtv, sigm, ta, tcl, tsk, tx, vpts, wetsk)
799
800   REAL ( wp )  :: acl, adu, aeff, cair, cbare, cclo, csum, ed,                &
801      emcl, emsk, enbal, enbal2, ere, erel, eres, esw, evap,                   &
802      facl, feff, h, hc, p, po, rbare, rclo, rdcl, rdsk, rsum,                 &
803      rtv, sigm, ta, tcl, tex, tsk, tx, vpex, vpts, wetsk, xx
804
805   INTEGER ( iwp ) :: count1
806
807   tx = ta
808   enbal2 = 0.
809
810    DO count1 = 0, 3
811     DO
812       hc = 2.67 + 6.5 * 0.1 ** 0.67
813       hc = hc * ( p / po ) ** 0.55
814
815!--    Radiation
816
817       aeff  = adu * feff
818       rbare = aeff * ( 1. - facl ) * emsk * sigm *                            &
819           ( ( tx + 273.2 ) ** 4. - ( tsk + 273.2 ) ** 4. )
820       rclo  = feff * acl * emcl * sigm *                                      &
821           ( ( tx + 273.2 ) ** 4. - ( tcl + 273.2 ) ** 4. )
822       rsum  = rbare + rclo
823
824!--    Covection
825
826       cbare = hc * ( tx - tsk ) * adu * ( 1. - facl )
827       cclo  = hc * ( tx - tcl ) * acl
828       csum  = cbare + cclo
829
830!--    Diffusion
831
832       ed = evap / ( rdsk + rdcl ) * adu * ( 1. - wetsk ) * ( 12. - vpts )
833
834!--    Respiration
835
836       tex  = 0.47 * tx + 21.
837       eres = cair * ( tx - tex ) * rtv
838       vpex = 6.11 * 10. ** ( 7.45 * tex / ( 235. + tex ) )
839       erel = 0.623 * evap / p * ( 12. - vpex ) * rtv
840       ere  = eres + erel
841
842!--    Energy ballance
843
844       enbal = h + ed + ere + esw + csum + rsum
845
846!--    Iteration concerning ta
847
848       IF ( count1 .EQ. 0 )  xx = 1.
849       IF ( count1 .EQ. 1 )  xx = 0.1
850       IF ( count1 .EQ. 2 )  xx = 0.01
851       IF ( count1 .EQ. 3 )  xx = 0.001
852       IF ( enbal .GT. 0. )  tx = tx - xx
853       IF ( enbal .LT. 0. )  tx = tx + xx
854       IF ( ( enbal .LE. 0. ) .AND. ( enbal2 .GT. 0. ) ) EXIT
855       IF ( ( enbal .GE. 0. ) .AND. ( enbal2 .LT. 0. ) ) EXIT
856
857       enbal2 = enbal
858     END DO
859   END DO
860   RETURN
861 END SUBROUTINE PET
862
863END MODULE biometeorology_mod
Note: See TracBrowser for help on using the repository browser.