Ignore:
Timestamp:
Sep 24, 2018 3:42:55 PM (3 years ago)
Author:
knoop
Message:

Modularization of all bulk cloud physics code components

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/radiation_model_mod.f90

    r3272 r3274  
    2828! -----------------
    2929! $Id$
     30! Modularization of all bulk cloud physics code components
     31!
     32! 3272 2018-09-24 10:16:32Z suehring
    3033! - split direct and diffusion shortwave radiation using RRTMG rather than using
    3134!   calc_diffusion_radiation, in case of RRTMG
     
    103106!
    104107! 3117 2018-07-11 09:59:11Z maronga
    105 ! Bugfix: water vapor was not transfered to RRTMG when cloud_physics = .F.
     108! Bugfix: water vapor was not transfered to RRTMG when bulk_cloud_model = .F.
    106109! Bugfix: changed the calculation of RRTMG boundary conditions (Mohamed Salim)
    107110! Bugfix: dry residual atmosphere is replaced by standard RRTMG atmosphere
     
    437440 
    438441    USE arrays_3d,                                                             &
    439         ONLY:  dzw, hyp, nc, pt, q, ql, zu, zw
     442        ONLY:  dzw, hyp, nc, pt, q, ql, zu, zw, exner, d_exner
     443
     444    USE basic_constants_and_equations_mod,                                     &
     445        ONLY:  c_p, g, lv_d_cp, l_v, pi, r_d, rho_l, solar_constant,           &
     446               barometric_formula
    440447
    441448    USE calc_mean_profile_mod,                                                 &
    442449        ONLY:  calc_mean_profile
    443450
    444     USE cloud_parameters,                                                      &
    445         ONLY:  cp, l_d_cp, l_v, r_d, rho_l
    446 
    447     USE constants,                                                             &
    448         ONLY:  pi
    449 
    450451    USE control_parameters,                                                    &
    451         ONLY:  cloud_droplets, cloud_physics, coupling_char, dz, g,            &
     452        ONLY:  cloud_droplets, coupling_char, dz,                              &
    452453               humidity,                                                       &
    453454               initializing_actions, io_blocks, io_group,                      &
    454455               latitude, longitude, large_scale_forcing, lsf_surf,             &
    455                message_string, microphysics_morrison, plant_canopy, pt_surface,&
     456               message_string, plant_canopy, pt_surface,&
    456457               rho_surface, surface_pressure, time_since_reference_point,      &
    457458               urban_surface, land_surface, end_time, spinup_time, dt_spinup
     
    475476    USE kinds
    476477
    477     USE microphysics_mod,                                                      &
    478         ONLY:  na_init, nc_const, sigma_gc
     478    USE bulk_cloud_model_mod,                                                  &
     479        ONLY:  bulk_cloud_model, microphysics_morrison, na_init, nc_const, sigma_gc
    479480
    480481#if defined ( __netcdf )
     
    561562                                                         /)
    562563
     564    REAL(wp), PARAMETER :: sigma_sb       = 5.67037321E-8_wp !< Stefan-Boltzmann constant
     565
    563566    INTEGER(iwp) :: albedo_type  = 9999999, & !< Albedo surface type
    564567                    dots_rad     = 0          !< starting index for timeseries output
     
    578581                                                        !< will be considered. However fewer SVFs are expected.
    579582                radiation_interactions_on = .TRUE.      !< namelist flag to force RTM activiation regardless to vertical urban/land surface and trees
    580 
    581 
    582     REAL(wp), PARAMETER :: sigma_sb       = 5.67037321E-8_wp,       & !< Stefan-Boltzmann constant
    583                            solar_constant = 1368.0_wp                 !< solar constant at top of atmosphere
    584583
    585584    REAL(wp) :: albedo = 9999999.9_wp,           & !< NAMELIST alpha
     
    22742273!
    22752274!--       Initialize RRTMG
    2276           IF ( lw_radiation )  CALL rrtmg_lw_ini ( cp )
    2277           IF ( sw_radiation )  CALL rrtmg_sw_ini ( cp )
     2275          IF ( lw_radiation )  CALL rrtmg_lw_ini ( c_p )
     2276          IF ( sw_radiation )  CALL rrtmg_sw_ini ( c_p )
    22782277
    22792278!
     
    23352334
    23362335       INTEGER(iwp) ::  l         !< running index for surface orientation
    2337 
    2338        REAL(wp)     ::  exn       !< Exner functions at surface
    2339        REAL(wp)     ::  exn1      !< Exner functions at first grid level or at urban layer top
    23402336       REAL(wp)     ::  pt1       !< potential temperature at first grid level or mean value at urban layer top
    23412337       REAL(wp)     ::  pt1_l     !< potential temperature at first grid level or mean value at urban layer top at local subdomain
     
    23552351!
    23562352!--    Calculate value of the Exner function at model surface
    2357        exn = (surface_pressure / 1000.0_wp )**0.286_wp
    23582353!
    23592354!--    In case averaged radiation is used, calculate mean temperature and
     
    23612356       IF ( average_radiation ) THEN   
    23622357          pt1   = 0.0_wp
    2363           IF ( cloud_physics  .OR.  cloud_droplets )  ql1   = 0.0_wp
     2358          IF ( bulk_cloud_model  .OR.  cloud_droplets )  ql1   = 0.0_wp
    23642359
    23652360          pt1_l = SUM( pt(nzut,nys:nyn,nxl:nxr) )
    2366           IF ( cloud_physics  .OR.  cloud_droplets  )  ql1_l = SUM( ql(nzut,nys:nyn,nxl:nxr) )
     2361          IF ( bulk_cloud_model  .OR.  cloud_droplets  )  ql1_l = SUM( ql(nzut,nys:nyn,nxl:nxr) )
    23672362
    23682363#if defined( __parallel )     
    23692364          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    23702365          CALL MPI_ALLREDUCE( pt1_l, pt1, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
    2371           IF ( cloud_physics  .OR.  cloud_droplets )                            &
     2366          IF ( bulk_cloud_model  .OR.  cloud_droplets )                            &
    23722367             CALL MPI_ALLREDUCE( ql1_l, ql1, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
    23732368#else
    23742369          pt1 = pt1_l
    2375           IF ( cloud_physics  .OR.  cloud_droplets )  ql1 = ql1_l
     2370          IF ( bulk_cloud_model  .OR.  cloud_droplets )  ql1 = ql1_l
    23762371#endif
    2377  
    2378           exn1 = ( hyp(nzut) / 100000.0_wp )**0.286_wp
    2379           IF ( cloud_physics  .OR.  cloud_droplets  )  pt1 = pt1 + l_d_cp / exn1 * ql1
     2372
     2373          IF ( bulk_cloud_model  .OR.  cloud_droplets  )  pt1 = pt1 + lv_d_cp / exner(nzut) * ql1
    23802374!
    23812375!--       Finally, divide by number of grid points
     
    24182412                k = nzut
    24192413
    2420                 exn1 = ( hyp(k+1) / 100000.0_wp )**0.286_wp
    2421 
    24222414                surf%rad_sw_in  = solar_constant * sky_trans * zenith(0)
    24232415                surf%rad_sw_out = albedo_urb * surf%rad_sw_in
    24242416               
    2425                 surf%rad_lw_in  = 0.8_wp * sigma_sb * (pt1 * exn1)**4
     2417                surf%rad_lw_in  = 0.8_wp * sigma_sb * (pt1 * exner(k+1))**4
    24262418
    24272419                surf%rad_lw_out = emissivity_urb * sigma_sb * (t_rad_urb)**4   &
     
    24432435                   j = surf%j(m)
    24442436                   k = surf%k(m)
    2445 
    2446                    exn1 = (hyp(k) / 100000.0_wp )**0.286_wp
    24472437
    24482438                   surf%rad_sw_in(m) = solar_constant * sky_trans * zenith(0)
     
    24692459                                        )                                      &
    24702460                                        * sigma_sb                             &
    2471                                         * ( surf%pt_surface(m) * exn )**4
     2461                                        * ( surf%pt_surface(m) * exner(nzb) )**4
    24722462
    24732463                   surf%rad_lw_out_change_0(m) =                               &
     
    24792469                                        surf%emissivity(ind_wat_win,m)         &
    24802470                                      ) * 3.0_wp * sigma_sb                    &
    2481                                       * ( surf%pt_surface(m) * exn )** 3
    2482 
    2483 
    2484                    IF ( cloud_physics  .OR.  cloud_droplets  )  THEN
    2485                       pt1 = pt(k,j,i) + l_d_cp / exn1 * ql(k,j,i)
    2486                       surf%rad_lw_in(m)  = 0.8_wp * sigma_sb * (pt1 * exn1)**4
     2471                                      * ( surf%pt_surface(m) * exner(nzb) )** 3
     2472
     2473
     2474                   IF ( bulk_cloud_model  .OR.  cloud_droplets  )  THEN
     2475                      pt1 = pt(k,j,i) + lv_d_cp / exner(k) * ql(k,j,i)
     2476                      surf%rad_lw_in(m)  = 0.8_wp * sigma_sb * (pt1 * exner(k))**4
    24872477                   ELSE
    2488                       surf%rad_lw_in(m)  = 0.8_wp * sigma_sb * (pt(k,j,i) * exn1)**4
     2478                      surf%rad_lw_in(m)  = 0.8_wp * sigma_sb * (pt(k,j,i) * exner(k))**4
    24892479                   ENDIF
    24902480
     
    25242514       INTEGER(iwp) ::  l         !< running index for surface orientation
    25252515
    2526        REAL(wp)     ::  exn       !< Exner functions at surface
    2527        REAL(wp)     ::  exn1      !< Exner functions at first grid level
    25282516       REAL(wp)     ::  pt1       !< potential temperature at first grid level or mean value at urban layer top
    25292517       REAL(wp)     ::  pt1_l     !< potential temperature at first grid level or mean value at urban layer top at local subdomain
     
    25312519       REAL(wp)     ::  ql1_l     !< liquid water mixing ratio at first grid level or mean value at urban layer top at local subdomain
    25322520
    2533        TYPE(surf_type), POINTER ::  surf !< pointer on respective surface type, used to generalize routine   
    2534 
    2535 !
    2536 !--    Calculate value of the Exner function
    2537        exn = (surface_pressure / 1000.0_wp )**0.286_wp
     2521       TYPE(surf_type), POINTER ::  surf !< pointer on respective surface type, used to generalize routine
     2522
    25382523!
    25392524!--    In case averaged radiation is used, calculate mean temperature and
     
    25412526       IF ( average_radiation ) THEN   
    25422527          pt1   = 0.0_wp
    2543           IF ( cloud_physics  .OR.  cloud_droplets )  ql1   = 0.0_wp
     2528          IF ( bulk_cloud_model  .OR.  cloud_droplets )  ql1   = 0.0_wp
    25442529
    25452530          pt1_l = SUM( pt(nzut,nys:nyn,nxl:nxr) )
    2546           IF ( cloud_physics  .OR.  cloud_droplets )  ql1_l = SUM( ql(nzut,nys:nyn,nxl:nxr) )
     2531          IF ( bulk_cloud_model  .OR.  cloud_droplets )  ql1_l = SUM( ql(nzut,nys:nyn,nxl:nxr) )
    25472532
    25482533#if defined( __parallel )     
    25492534          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    25502535          CALL MPI_ALLREDUCE( pt1_l, pt1, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
    2551           IF ( cloud_physics  .OR.  cloud_droplets )                             &
     2536          IF ( bulk_cloud_model  .OR.  cloud_droplets )                             &
    25522537             CALL MPI_ALLREDUCE( ql1_l, ql1, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
    25532538#else
    25542539          pt1 = pt1_l
    2555           IF ( cloud_physics  .OR.  cloud_droplets )  ql1 = ql1_l
     2540          IF ( bulk_cloud_model  .OR.  cloud_droplets )  ql1 = ql1_l
    25562541#endif
    2557           IF ( cloud_physics  .OR.  cloud_droplets )  pt1 = pt1 + l_d_cp / exn1 * ql1
     2542          IF ( bulk_cloud_model  .OR.  cloud_droplets )  pt1 = pt1 + lv_d_cp / exner(nzut+1) * ql1
    25582543!
    25592544!--       Finally, divide by number of grid points
     
    25952580             IF ( average_radiation ) THEN
    25962581
    2597                 ! set height above canopy
    2598                 k = nzut
    2599 
    26002582                surf%rad_net = net_radiation
    2601 ! MS: Wyh k + 1 ?
    2602                 exn1 = (hyp(k+1) / 100000.0_wp )**0.286_wp
    2603 
    2604                 surf%rad_lw_in  = 0.8_wp * sigma_sb * (pt1 * exn1)**4
     2583
     2584                surf%rad_lw_in  = 0.8_wp * sigma_sb * (pt1 * exner(nzut+1))**4
    26052585
    26062586                surf%rad_lw_out = emissivity_urb * sigma_sb * (t_rad_urb)**4   &
     
    26362616                   surf%rad_net(m) = net_radiation
    26372617
    2638                    exn1 = (hyp(k) / 100000.0_wp )**0.286_wp
    2639 
    2640                    IF ( cloud_physics  .OR.  cloud_droplets )  THEN
    2641                       pt1 = pt(k,j,i) + l_d_cp / exn1 * ql(k,j,i)
    2642                       surf%rad_lw_in(m)  = 0.8_wp * sigma_sb * (pt1 * exn1)**4
     2618                   IF ( bulk_cloud_model  .OR.  cloud_droplets )  THEN
     2619                      pt1 = pt(k,j,i) + lv_d_cp / exner(k) * ql(k,j,i)
     2620                      surf%rad_lw_in(m)  = 0.8_wp * sigma_sb * (pt1 * exner(k))**4
    26432621                   ELSE
    26442622                      surf%rad_lw_in(m)  = 0.8_wp * sigma_sb *                 &
    2645                                              ( pt(k,j,i) * exn1 )**4
     2623                                             ( pt(k,j,i) * exner(k) )**4
    26462624                   ENDIF
    26472625
     
    26562634                                        )                                      &
    26572635                                      * sigma_sb                               &
    2658                                       * ( surf%pt_surface(m) * exn )**4
     2636                                      * ( surf%pt_surface(m) * exner(nzb) )**4
    26592637
    26602638                   surf%rad_sw_in(m) = ( surf%rad_net(m) - surf%rad_lw_in(m)   &
     
    29542932          rrtm_tlev(0,nzb+1) = t_rad_urb
    29552933
    2956           IF ( cloud_physics )  THEN
     2934          IF ( bulk_cloud_model )  THEN
    29572935
    29582936             CALL calc_mean_profile( ql, 54 )
     
    29622940             DO k = nzb+1, nzt+1
    29632941                rrtm_tlay(0,k) = pt_av(k) * ( (hyp(k) ) / 100000._wp       &
    2964                                  )**.286_wp + l_d_cp * ql_av(k)
     2942                                 )**.286_wp + lv_d_cp * ql_av(k)
    29652943                rrtm_h2ovmr(0,k) = mol_mass_air_d_wv * (q_av(k) - ql_av(k))
    29662944             ENDDO
     
    30132991          rrtm_icld   = 0
    30142992
    3015           IF ( cloud_physics )  THEN
     2993          IF ( bulk_cloud_model )  THEN
    30162994             DO k = nzb+1, nzt+1
    30172995                rrtm_cliqwp(0,k) =  ql_av(k) * 1000._wp *                   &
     
    31333111!--             Prepare profiles of temperature and H2O volume mixing ratio
    31343112                DO  m = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
    3135                    rrtm_tlev(0,nzb+1) = surf_lsm_h%pt_surface(m) * &
    3136                                     ( surface_pressure / 1000.0_wp )**0.286_wp
     3113                   rrtm_tlev(0,nzb+1) = surf_lsm_h%pt_surface(m) * exner(nzb)
    31373114                ENDDO
    31383115                DO  m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
    3139                    rrtm_tlev(0,nzb+1) = surf_usm_h%pt_surface(m) * &
    3140                                     ( surface_pressure / 1000.0_wp )**0.286_wp
     3116                   rrtm_tlev(0,nzb+1) = surf_usm_h%pt_surface(m) * exner(nzb)
    31413117                ENDDO
    31423118
    31433119
    3144                 IF ( cloud_physics )  THEN
     3120                IF ( bulk_cloud_model )  THEN
    31453121                   DO k = nzb+1, nzt+1
    3146                       rrtm_tlay(0,k) = pt(k,j,i) * ( (hyp(k) ) / 100000.0_wp   &
    3147                                        )**0.286_wp + l_d_cp * ql(k,j,i)
     3122                      rrtm_tlay(0,k) = pt(k,j,i) * exner(k)                    &
     3123                                        + lv_d_cp * ql(k,j,i)
    31483124                      rrtm_h2ovmr(0,k) = mol_mass_air_d_wv * (q(k,j,i) - ql(k,j,i))
    31493125                   ENDDO
    31503126                ELSEIF ( cloud_droplets )  THEN
    31513127                   DO k = nzb+1, nzt+1
    3152                       rrtm_tlay(0,k) = pt(k,j,i) * ( (hyp(k) ) / 100000.0_wp   &
    3153                                        )**0.286_wp + l_d_cp * ql(k,j,i)
     3128                      rrtm_tlay(0,k) = pt(k,j,i) * exner(k)                     &
     3129                                        + lv_d_cp * ql(k,j,i)
    31543130                      rrtm_h2ovmr(0,k) = mol_mass_air_d_wv * q(k,j,i)
    31553131                   ENDDO
    31563132                ELSE
    31573133                   DO k = nzb+1, nzt+1
    3158                       rrtm_tlay(0,k) = pt(k,j,i) * ( (hyp(k) ) / 100000.0_wp   &
    3159                                        )**0.286_wp
     3134                      rrtm_tlay(0,k) = pt(k,j,i) * exner(k)
    31603135                   ENDDO
    31613136
     
    32023177                rrtm_icld   = 0
    32033178
    3204                 IF ( cloud_physics  .OR.  cloud_droplets )  THEN
     3179                IF ( bulk_cloud_model  .OR.  cloud_droplets )  THEN
    32053180                   DO k = nzb+1, nzt+1
    32063181                      rrtm_cliqwp(0,k) =  ql(k,j,i) * 1000.0_wp *              &
     
    32143189!
    32153190!--                      Calculate cloud droplet effective radius
    3216                          IF ( cloud_physics )  THEN
     3191                         IF ( bulk_cloud_model )  THEN
    32173192!
    32183193!--                         Calculete effective droplet radius. In case of using
     
    32853260                               surf_lsm_h%frac(ind_wat_win,m)   *              &
    32863261                               surf_lsm_h%emissivity(ind_wat_win,m)
    3287                    rrtm_tsfc = surf_lsm_h%pt_surface(m) *                      &
    3288                                        (surface_pressure / 1000.0_wp )**0.286_wp
     3262                   rrtm_tsfc = surf_lsm_h%pt_surface(m) * exner(nzb)
    32893263                ENDDO             
    32903264                DO  m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
     
    32953269                               surf_usm_h%frac(ind_wat_win,m)   *              &
    32963270                               surf_usm_h%emissivity(ind_wat_win,m)
    3297                    rrtm_tsfc = surf_usm_h%pt_surface(m) *                      &
    3298                                        (surface_pressure / 1000.0_wp )**0.286_wp
     3271                   rrtm_tsfc = surf_usm_h%pt_surface(m) * exner(nzb)
    32993272                ENDDO
    33003273!
     
    39853958       ALLOCATE ( rrtm_plev(0:0,nzb+1:nzt_rad+2)   )
    39863959
    3987        t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp
     3960       t_surface = pt_surface * exner(nzb)
    39883961       DO k = nzb+1, nzt+1
    39893962          rrtm_play(0,k) = hyp(k) * 0.01_wp
    3990           rrtm_plev(0,k) = surface_pressure * ( (t_surface - g/cp * zw(k-1)) / &
    3991                          t_surface )**(1.0_wp/0.286_wp)
     3963          rrtm_plev(0,k) = barometric_formula(zw(k-1),                         &
     3964                              pt_surface * exner(nzb), &
     3965                              surface_pressure )
    39923966       ENDDO
    39933967
     
    44064380 SUBROUTINE radiation_tendency_ij ( i, j, tend )
    44074381
    4408     USE cloud_parameters,                                                      &
    4409         ONLY:  pt_d_t
    4410 
    44114382    IMPLICIT NONE
    44124383
     
    44214392       DO k = nzb+1, nzt+1
    44224393          tend(k,j,i) = tend(k,j,i) + (rad_lw_hr(k,j,i) + rad_sw_hr(k,j,i))    &
    4423                                          * pt_d_t(k) * d_seconds_hour
     4394                                         * d_exner(k) * d_seconds_hour
    44244395       ENDDO
    44254396#endif
     
    44374408 SUBROUTINE radiation_tendency ( tend )
    44384409
    4439     USE cloud_parameters,                                                      &
    4440         ONLY:  pt_d_t
    4441 
    44424410    USE indices,                                                               &
    44434411        ONLY:  nxl, nxr, nyn, nys
     
    44574425             DO k = nzb+1, nzt+1
    44584426                tend(k,j,i) = tend(k,j,i) + ( rad_lw_hr(k,j,i)                 &
    4459                                           +  rad_sw_hr(k,j,i) ) * pt_d_t(k)    &
     4427                                          +  rad_sw_hr(k,j,i) ) * d_exner(k)   &
    44604428                                          * d_seconds_hour
    44614429             ENDDO
     
    45254493#if ! defined( __nopointer )
    45264494     IF ( plant_canopy )  THEN
    4527          pchf_prep(:) = r_d * (hyp(nzub:nzut) / 100000.0_wp)**0.286_wp &
    4528                      / (cp * hyp(nzub:nzut) * dx*dy*dz(1)) !< equals to 1 / (rho * c_p * Vbox * T)
    4529          pctf_prep(:) = r_d * (hyp(nzub:nzut) / 100000.0_wp)**0.286_wp &
     4495         pchf_prep(:) = r_d * exner(nzub:nzut)                                &
     4496                     / (c_p * hyp(nzub:nzut) * dx*dy*dz(1)) !< equals to 1 / (rho * c_p * Vbox * T)
     4497         pctf_prep(:) = r_d * exner(nzub:nzut)                                &
    45304498                     / (l_v * hyp(nzub:nzut) * dx*dy*dz(1))
    45314499     ENDIF
     
    50325000!         k                = surf_lsm_h%k(m)
    50335001!         pt_surf_urb_l    = pt_surf_urb_l + surf_lsm_h%pt_surface(m)            &
    5034 !                            * ( hyp(k) / 100000.0_wp )**0.286_wp
     5002!                            * exner(k)
    50355003!         count_surfaces_l = count_surfaces_l + 1.0_wp
    50365004!      ENDDO
     
    50385006!         k                = surf_usm_h%k(m)
    50395007!         pt_surf_urb_l    = pt_surf_urb_l + surf_usm_h%pt_surface(m)            &
    5040 !                            * ( hyp(k) / 100000.0_wp )**0.286_wp
     5008!                            * exner(k)
    50415009!         count_surfaces_l = count_surfaces_l + 1.0_wp
    50425010!      ENDDO
     
    50455013!            k                = surf_lsm_v(l)%k(m)
    50465014!            pt_surf_urb_l    = pt_surf_urb_l + surf_lsm_v(l)%pt_surface(m)      &
    5047 !                            * ( hyp(k) / 100000.0_wp )**0.286_wp
     5015!                            * exner(k)
    50485016!            count_surfaces_l = count_surfaces_l + 1.0_wp
    50495017!         ENDDO
    50505018!         DO  m = 1, surf_usm_v(l)%ns
    50515019!            k                = surf_usm_v(l)%k(m)
    5052 !            pt_surf_urb_l    = pt_surf_urb_l + surf_usm_v(l)%pt_surface(m)      &
    5053 !                            * ( hyp(k) / 100000.0_wp )**0.286_wp
     5020!            pt_surf_urb_l    = pt_surf_urb_l + surf_usm_v(l)%pt_surface(m) * exner(k)
    50545021!            count_surfaces_l = count_surfaces_l + 1.0_wp
    50555022!         ENDDO
Note: See TracChangeset for help on using the changeset viewer.