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

Last change on this file since 3337 was 3337, checked in by kanani, 3 years ago

reintegrate branch resler to trunk

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