1 2! KGEN-generated Fortran source file 3! 4! Filename : mo_psrad_interface.f90 5! Generated at: 2015-02-19 15:30:28 6! KGEN version: 0.4.4 7 8 9 10 MODULE mo_psrad_interface 11 USE mo_spec_sampling, only : read_var_mod5 => kgen_read_var 12 USE mo_kind, ONLY: wp 13 USE mo_rrtm_params, ONLY: nbndlw 14 USE mo_rrtm_params, ONLY: maxinpx 15 USE mo_rrtm_params, ONLY: maxxsec 16 USE mo_lrtm_driver, ONLY: lrtm 17 USE mo_spec_sampling, ONLY: spec_sampling_strategy 18 IMPLICIT NONE 19 PUBLIC lw_strat 20 PUBLIC read_externs_mo_psrad_interface 21 INTEGER, PARAMETER :: kgen_dp = selected_real_kind(15, 307) 22 PUBLIC psrad_interface 23 type, public :: check_t 24 logical :: Passed 25 integer :: numFatal 26 integer :: numTotal 27 integer :: numIdentical 28 integer :: numWarning 29 integer :: VerboseLevel 30 real(kind=kgen_dp) :: tolerance 31 end type check_t 32 TYPE(spec_sampling_strategy), save :: lw_strat 33 !< Spectral sampling strategies for longwave, shortwave 34 INTEGER, parameter :: rng_seed_size = 4 35 CONTAINS 36 37 ! module extern variables 38 39 SUBROUTINE read_externs_mo_psrad_interface(kgen_unit) 40 integer, intent(in) :: kgen_unit 41 call read_var_mod5(lw_strat, kgen_unit) 42 END SUBROUTINE read_externs_mo_psrad_interface 43 44 subroutine kgen_init_check(check,tolerance) 45 type(check_t), intent(inout) :: check 46 real(kind=kgen_dp), intent(in), optional :: tolerance 47 check%Passed = .TRUE. 48 check%numFatal = 0 49 check%numWarning = 0 50 check%numTotal = 0 51 check%numIdentical = 0 52 check%VerboseLevel = 1 53 if(present(tolerance)) then 54 check%tolerance = tolerance 55 else 56 check%tolerance = 1.E-14 57 endif 58 end subroutine kgen_init_check 59 subroutine kgen_print_check(kname, check) 60 character(len=*) :: kname 61 type(check_t), intent(in) :: check 62 write (*,*) 63 write (*,*) TRIM(kname),' KGENPrtCheck: Tolerance for normalized RMS: ',check%tolerance 64 write (*,*) TRIM(kname),' KGENPrtCheck: Number of variables checked: ',check%numTotal 65 write (*,*) TRIM(kname),' KGENPrtCheck: Number of Identical results: ',check%numIdentical 66 write (*,*) TRIM(kname),' KGENPrtCheck: Number of warnings detected: ',check%numWarning 67 write (*,*) TRIM(kname),' KGENPrtCheck: Number of fatal errors detected: ', check%numFatal 68 if (check%numFatal> 0) then 69 write(*,*) TRIM(kname),' KGENPrtCheck: verification FAILED' 70 else 71 write(*,*) TRIM(kname),' KGENPrtCheck: verification PASSED' 72 endif 73 end subroutine kgen_print_check 74 !--------------------------------------------------------------------------- 75 !> 76 !! @brief Sets up (initializes) radation routines 77 !! 78 !! @remarks 79 !! Modify preset variables of module MO_RADIATION which control the 80 !! configuration of the radiation scheme. 81 ! 82 83 !----------------------------------------------------------------------------- 84 !> 85 !! @brief arranges input and calls rrtm sw and lw routines 86 !! 87 !! @par Revision History 88 !! Original Source Rewritten and renamed by B. Stevens (2009-08) 89 !! 90 !! @remarks 91 !! Because the RRTM indexes vertical levels differently than ECHAM a chief 92 !! function of thise routine is to reorder the input in the vertical. In 93 !! addition some cloud physical properties are prescribed, which are 94 !! required to derive cloud optical properties 95 !! 96 !! @par The gases are passed into RRTM via two multi-constituent arrays: 97 !! zwkl and wx_r. zwkl has maxinpx species and wx_r has maxxsec species 98 !! The species are identifed as follows. 99 !! ZWKL [#/cm2] WX_R [#/cm2] 100 !! index = 1 => H20 index = 1 => n/a 101 !! index = 2 => CO2 index = 2 => CFC11 102 !! index = 3 => O3 index = 3 => CFC12 103 !! index = 4 => N2O index = 4 => n/a 104 !! index = 5 => n/a 105 !! index = 6 => CH4 106 !! index = 7 => O2 107 ! 108 109 SUBROUTINE psrad_interface(kbdim, klev, nb_sw, kproma, ktrac, tk_sfc, kgen_unit) 110 integer, intent(in) :: kgen_unit 111 112 ! read interface 113 !interface kgen_read_var 114 ! procedure read_var_real_wp_dim2 115 ! procedure read_var_real_wp_dim1 116 ! procedure read_var_real_wp_dim3 117 ! procedure read_var_integer_4_dim2 118 !end interface kgen_read_var 119 120 121 122 ! verification interface 123 !interface kgen_verify_var 124 ! procedure verify_var_logical 125 ! procedure verify_var_integer 126 ! procedure verify_var_real 127 ! procedure verify_var_character 128 ! procedure verify_var_real_wp_dim2 129 ! procedure verify_var_real_wp_dim1 130 ! procedure verify_var_real_wp_dim3 131 ! procedure verify_var_integer_4_dim2 132 !end interface kgen_verify_var 133 134 INTEGER*8 :: kgen_intvar, start_clock, stop_clock, rate_clock 135 TYPE(check_t):: check_status 136 REAL(KIND=kgen_dp) :: tolerance 137 INTEGER, intent(in) :: kbdim 138 INTEGER, intent(in) :: klev 139 INTEGER, intent(in) :: nb_sw 140 INTEGER, intent(in) :: kproma 141 INTEGER, intent(in) :: ktrac 142 !< aerosol control 143 !< number of longitudes 144 !< first dimension of 2-d arrays 145 !< first dimension of 2-d arrays 146 !< number of levels 147 !< number of tracers 148 !< type of convection 149 !< number of shortwave bands 150 !< land sea mask, land=.true. 151 !< glacier mask, glacier=.true. 152 REAL(KIND=wp), intent(in) :: tk_sfc(kbdim) 153 !< surface emissivity 154 !< mu0 for solar zenith angle 155 !< geopotential above ground 156 !< surface albedo for vis range and dir light 157 !< surface albedo for NIR range and dir light 158 !< surface albedo for vis range and dif light 159 !< surface albedo for NIR range and dif light 160 !< full level pressure in Pa 161 !< half level pressure in Pa 162 !< surface pressure in Pa 163 !< full level temperature in K 164 !< half level temperature in K 165 !< surface temperature in K 166 !< specific humidity in g/g 167 !< specific liquid water content 168 !< specific ice content in g/g 169 !< cloud nuclei concentration 170 !< fractional cloud cover 171 !< total cloud cover in m2/m2 172 !< o3 mass mixing ratio 173 !< co2 mass mixing ratio 174 !< ch4 mass mixing ratio 175 !< n2o mass mixing ratio 176 !< cfc volume mixing ratio 177 !< o2 mass mixing ratio 178 !< tracer mass mixing ratios 179 !< upward LW flux profile, all sky 180 !< upward LW flux profile, clear sky 181 !< downward LW flux profile, all sky 182 !< downward LW flux profile, clear sky 183 !< upward SW flux profile, all sky 184 !< upward SW flux profile, clear sky 185 !< downward SW flux profile, all sky 186 !< downward SW flux profile, clear sky 187 !< Visible (250-680) fraction of net surface radiation 188 !< Downward Photosynthetically Active Radiation (PAR) at surface 189 !< Diffuse fraction of downward surface near-infrared radiation 190 !< Diffuse fraction of downward surface visible radiation 191 !< Diffuse fraction of downward surface PAR 192 ! ------------------------------------------------------------------------------------- 193 !< loop indicies 194 !< index for clear or cloudy 195 REAL(KIND=wp) :: zsemiss (kbdim,nbndlw) 196 REAL(KIND=wp) :: pm_sfc (kbdim) 197 !< LW surface emissivity by band 198 !< pressure thickness in Pa 199 !< surface pressure in mb 200 !< pressure thickness 201 !< scratch array 202 ! 203 ! --- vertically reversed _vr variables 204 ! 205 REAL(KIND=wp) :: cld_frc_vr(kbdim,klev) 206 REAL(KIND=wp) :: aer_tau_lw_vr(kbdim,klev,nbndlw) 207 REAL(KIND=wp) :: pm_fl_vr (kbdim,klev) 208 REAL(KIND=wp) :: tk_fl_vr (kbdim,klev) 209 REAL(KIND=wp) :: tk_hl_vr (kbdim,klev+1) 210 REAL(KIND=wp) :: cld_tau_lw_vr(kbdim,klev,nbndlw) 211 REAL(KIND=wp) :: wkl_vr (kbdim,maxinpx,klev) 212 REAL(KIND=wp) :: wx_vr (kbdim,maxxsec,klev) 213 REAL(KIND=wp) :: col_dry_vr(kbdim,klev) 214 !< number of molecules/cm2 of 215 !< full level pressure [mb] 216 !< half level pressure [mb] 217 !< full level temperature [K] 218 !< half level temperature [K] 219 !< cloud nuclei concentration 220 !< secure cloud fraction 221 !< specific ice water content 222 !< ice water content per volume 223 !< ice water path in g/m2 224 !< specific liquid water content 225 !< liquid water path in g/m2 226 !< liquid water content per 227 !< effective radius of liquid 228 !< effective radius of ice 229 !< number of molecules/cm2 of 230 !< number of molecules/cm2 of 231 !< LW optical thickness of clouds 232 !< extincion 233 !< asymmetry factor 234 !< single scattering albedo 235 !< LW optical thickness of aerosols 236 !< aerosol optical thickness 237 !< aerosol asymmetry factor 238 !< aerosol single scattering albedo 239 REAL(KIND=wp) :: flx_uplw_vr(kbdim,klev+1) 240 REAL(KIND=wp), allocatable :: ref_flx_uplw_vr(:,:) 241 REAL(KIND=wp) :: flx_uplw_clr_vr(kbdim,klev+1) 242 REAL(KIND=wp), allocatable :: ref_flx_uplw_clr_vr(:,:) 243 REAL(KIND=wp) :: flx_dnlw_clr_vr(kbdim,klev+1) 244 REAL(KIND=wp), allocatable :: ref_flx_dnlw_clr_vr(:,:) 245 REAL(KIND=wp) :: flx_dnlw_vr(kbdim,klev+1) 246 REAL(KIND=wp), allocatable :: ref_flx_dnlw_vr(:,:) 247 !< upward flux, total sky 248 !< upward flux, clear sky 249 !< downward flux, total sky 250 !< downward flux, clear sky 251 ! 252 ! Random seeds for sampling. Needs to get somewhere upstream 253 ! 254 INTEGER :: rnseeds(kbdim,rng_seed_size) 255 INTEGER, allocatable :: ref_rnseeds(:,:) 256 ! 257 ! Number of g-points per time step. Determine here to allow automatic array allocation in 258 ! lrtm, srtm subroutines. 259 ! 260 INTEGER :: n_gpts_ts 261 ! 1.0 Constituent properties 262 !-------------------------------- 263 !IBM* ASSERT(NODEPS) 264 ! 265 ! --- control for zero, infintesimal or negative cloud fractions 266 ! 267 ! 268 ! --- main constituent reordering 269 ! 270 !IBM* ASSERT(NODEPS) 271 !IBM* ASSERT(NODEPS) 272 !IBM* ASSERT(NODEPS) 273 ! 274 ! --- CFCs are in volume mixing ratio 275 ! 276 !IBM* ASSERT(NODEPS) 277 ! 278 ! -- Convert to molecules/cm^2 279 ! 280 ! 281 ! 2.0 Surface Properties 282 ! -------------------------------- 283 ! 284 ! 3.0 Particulate Optical Properties 285 ! -------------------------------- 286 ! 287 ! 3.5 Interface for submodels that provide aerosol and/or cloud radiative properties: 288 ! ----------------------------------------------------------------------------------- 289 ! 290 ! 4.0 Radiative Transfer Routines 291 ! -------------------------------- 292 ! 293 ! Seeds for random numbers come from least significant digits of pressure field 294 ! 295 tolerance = 1.E-12 296 CALL kgen_init_check(check_status, tolerance) 297 READ(UNIT=kgen_unit) zsemiss 298 READ(UNIT=kgen_unit) pm_sfc 299 READ(UNIT=kgen_unit) cld_frc_vr 300 READ(UNIT=kgen_unit) aer_tau_lw_vr 301 READ(UNIT=kgen_unit) pm_fl_vr 302 READ(UNIT=kgen_unit) tk_fl_vr 303 READ(UNIT=kgen_unit) tk_hl_vr 304 READ(UNIT=kgen_unit) cld_tau_lw_vr 305 READ(UNIT=kgen_unit) wkl_vr 306 READ(UNIT=kgen_unit) wx_vr 307 READ(UNIT=kgen_unit) col_dry_vr 308 READ(UNIT=kgen_unit) flx_uplw_vr 309 READ(UNIT=kgen_unit) flx_uplw_clr_vr 310 READ(UNIT=kgen_unit) flx_dnlw_clr_vr 311 READ(UNIT=kgen_unit) flx_dnlw_vr 312 READ(UNIT=kgen_unit) rnseeds 313 READ(UNIT=kgen_unit) n_gpts_ts 314 315 !call kgen_read_var(ref_flx_uplw_vr, kgen_unit) 316 !call kgen_read_var(ref_flx_uplw_clr_vr, kgen_unit) 317 !call kgen_read_var(ref_flx_dnlw_clr_vr, kgen_unit) 318 !call kgen_read_var(ref_flx_dnlw_vr, kgen_unit) 319 !call kgen_read_var(ref_rnseeds, kgen_unit) 320 call read_var_real_wp_dim2(ref_flx_uplw_vr, kgen_unit) 321 call read_var_real_wp_dim2(ref_flx_uplw_clr_vr, kgen_unit) 322 call read_var_real_wp_dim2(ref_flx_dnlw_clr_vr, kgen_unit) 323 call read_var_real_wp_dim2(ref_flx_dnlw_vr, kgen_unit) 324 call read_var_integer_4_dim2(ref_rnseeds, kgen_unit) 325 326 ! call to kernel 327 CALL lrtm(kproma, kbdim, klev, pm_fl_vr, pm_sfc, tk_fl_vr, tk_hl_vr, tk_sfc, wkl_vr, wx_vr, col_dry_vr, zsemiss, cld_frc_vr, cld_tau_lw_vr, aer_tau_lw_vr, rnseeds, lw_strat, n_gpts_ts, flx_uplw_vr, flx_dnlw_vr, flx_uplw_clr_vr, flx_dnlw_clr_vr) 328 ! kernel verification for output variables 329 call verify_var_real_wp_dim2("flx_uplw_vr", check_status, flx_uplw_vr, ref_flx_uplw_vr) 330 call verify_var_real_wp_dim2("flx_uplw_clr_vr", check_status, flx_uplw_clr_vr, ref_flx_uplw_clr_vr) 331 call verify_var_real_wp_dim2("flx_dnlw_clr_vr", check_status, flx_dnlw_clr_vr, ref_flx_dnlw_clr_vr) 332 call verify_var_real_wp_dim2("flx_dnlw_vr", check_status, flx_dnlw_vr, ref_flx_dnlw_vr) 333 call verify_var_integer_4_dim2("rnseeds", check_status, rnseeds, ref_rnseeds) 334 CALL kgen_print_check("lrtm", check_status) 335 CALL system_clock(start_clock, rate_clock) 336 DO kgen_intvar=1,100 337 CALL lrtm(kproma, kbdim, klev, pm_fl_vr, pm_sfc, tk_fl_vr, tk_hl_vr, tk_sfc, wkl_vr, wx_vr, col_dry_vr, zsemiss, cld_frc_vr, cld_tau_lw_vr, aer_tau_lw_vr, rnseeds, lw_strat, n_gpts_ts, flx_uplw_vr, flx_dnlw_vr, flx_uplw_clr_vr, flx_dnlw_clr_vr) 338 END DO 339 CALL system_clock(stop_clock, rate_clock) 340 WRITE(*,*) 341 PRINT *, "Elapsed time (sec): ", (stop_clock - start_clock)/REAL(rate_clock*100) 342 ! 343 ! Reset random seeds so SW doesn't depend on what's happened in LW but is also independent 344 ! 345 ! 346 ! Potential pitfall - we're passing every argument but some may not be present 347 ! 348 ! 349 ! 5.0 Post Processing 350 ! -------------------------------- 351 ! 352 ! Lw fluxes are vertically revered but SW fluxes are not 353 ! 354 ! 355 ! 6.0 Interface for submodel diagnosics after radiation calculation: 356 ! ------------------------------------------------------------------ 357 CONTAINS 358 359 ! read subroutines 360 subroutine read_var_real_wp_dim2(var, kgen_unit) 361 integer, intent(in) :: kgen_unit 362 real(kind=wp), intent(out), dimension(:,:), allocatable :: var 363 integer, dimension(2,2) :: kgen_bound 364 logical is_save 365 366 READ(UNIT = kgen_unit) is_save 367 if ( is_save ) then 368 READ(UNIT = kgen_unit) kgen_bound(1, 1) 369 READ(UNIT = kgen_unit) kgen_bound(2, 1) 370 READ(UNIT = kgen_unit) kgen_bound(1, 2) 371 READ(UNIT = kgen_unit) kgen_bound(2, 2) 372 ALLOCATE(var(kgen_bound(2, 1) - kgen_bound(1, 1) + 1, kgen_bound(2, 2) - kgen_bound(1, 2) + 1)) 373 READ(UNIT = kgen_unit) var 374 end if 375 end subroutine 376 subroutine read_var_real_wp_dim1(var, kgen_unit) 377 integer, intent(in) :: kgen_unit 378 real(kind=wp), intent(out), dimension(:), allocatable :: var 379 integer, dimension(2,1) :: kgen_bound 380 logical is_save 381 382 READ(UNIT = kgen_unit) is_save 383 if ( is_save ) then 384 READ(UNIT = kgen_unit) kgen_bound(1, 1) 385 READ(UNIT = kgen_unit) kgen_bound(2, 1) 386 ALLOCATE(var(kgen_bound(2, 1) - kgen_bound(1, 1) + 1)) 387 READ(UNIT = kgen_unit) var 388 end if 389 end subroutine 390 subroutine read_var_real_wp_dim3(var, kgen_unit) 391 integer, intent(in) :: kgen_unit 392 real(kind=wp), intent(out), dimension(:,:,:), allocatable :: var 393 integer, dimension(2,3) :: kgen_bound 394 logical is_save 395 396 READ(UNIT = kgen_unit) is_save 397 if ( is_save ) then 398 READ(UNIT = kgen_unit) kgen_bound(1, 1) 399 READ(UNIT = kgen_unit) kgen_bound(2, 1) 400 READ(UNIT = kgen_unit) kgen_bound(1, 2) 401 READ(UNIT = kgen_unit) kgen_bound(2, 2) 402 READ(UNIT = kgen_unit) kgen_bound(1, 3) 403 READ(UNIT = kgen_unit) kgen_bound(2, 3) 404 ALLOCATE(var(kgen_bound(2, 1) - kgen_bound(1, 1) + 1, kgen_bound(2, 2) - kgen_bound(1, 2) + 1, kgen_bound(2, 3) - kgen_bound(1, 3) + 1)) 405 READ(UNIT = kgen_unit) var 406 end if 407 end subroutine 408 subroutine read_var_integer_4_dim2(var, kgen_unit) 409 integer, intent(in) :: kgen_unit 410 integer(kind=4), intent(out), dimension(:,:), allocatable :: var 411 integer, dimension(2,2) :: kgen_bound 412 logical is_save 413 414 READ(UNIT = kgen_unit) is_save 415 if ( is_save ) then 416 READ(UNIT = kgen_unit) kgen_bound(1, 1) 417 READ(UNIT = kgen_unit) kgen_bound(2, 1) 418 READ(UNIT = kgen_unit) kgen_bound(1, 2) 419 READ(UNIT = kgen_unit) kgen_bound(2, 2) 420 ALLOCATE(var(kgen_bound(2, 1) - kgen_bound(1, 1) + 1, kgen_bound(2, 2) - kgen_bound(1, 2) + 1)) 421 READ(UNIT = kgen_unit) var 422 end if 423 end subroutine 424 425 subroutine verify_var_logical(varname, check_status, var, ref_var) 426 character(*), intent(in) :: varname 427 type(check_t), intent(inout) :: check_status 428 logical, intent(in) :: var, ref_var 429 430 check_status%numTotal = check_status%numTotal + 1 431 IF ( var .eqv. ref_var ) THEN 432 check_status%numIdentical = check_status%numIdentical + 1 433 if(check_status%verboseLevel > 1) then 434 WRITE(*,*) 435 WRITE(*,*) trim(adjustl(varname)), " is IDENTICAL( ", var, " )." 436 endif 437 ELSE 438 if(check_status%verboseLevel > 1) then 439 WRITE(*,*) 440 WRITE(*,*) trim(adjustl(varname)), " is NOT IDENTICAL." 441 if(check_status%verboseLevel > 2) then 442 WRITE(*,*) "KERNEL: ", var 443 WRITE(*,*) "REF. : ", ref_var 444 endif 445 endif 446 check_status%numFatal = check_status%numFatal + 1 447 END IF 448 end subroutine 449 450 subroutine verify_var_integer(varname, check_status, var, ref_var) 451 character(*), intent(in) :: varname 452 type(check_t), intent(inout) :: check_status 453 integer, intent(in) :: var, ref_var 454 455 check_status%numTotal = check_status%numTotal + 1 456 IF ( var == ref_var ) THEN 457 check_status%numIdentical = check_status%numIdentical + 1 458 if(check_status%verboseLevel > 1) then 459 WRITE(*,*) 460 WRITE(*,*) trim(adjustl(varname)), " is IDENTICAL( ", var, " )." 461 endif 462 ELSE 463 if(check_status%verboseLevel > 0) then 464 WRITE(*,*) 465 WRITE(*,*) trim(adjustl(varname)), " is NOT IDENTICAL." 466 if(check_status%verboseLevel > 2) then 467 WRITE(*,*) "KERNEL: ", var 468 WRITE(*,*) "REF. : ", ref_var 469 endif 470 endif 471 check_status%numFatal = check_status%numFatal + 1 472 END IF 473 end subroutine 474 475 subroutine verify_var_real(varname, check_status, var, ref_var) 476 character(*), intent(in) :: varname 477 type(check_t), intent(inout) :: check_status 478 real, intent(in) :: var, ref_var 479 480 check_status%numTotal = check_status%numTotal + 1 481 IF ( var == ref_var ) THEN 482 check_status%numIdentical = check_status%numIdentical + 1 483 if(check_status%verboseLevel > 1) then 484 WRITE(*,*) 485 WRITE(*,*) trim(adjustl(varname)), " is IDENTICAL( ", var, " )." 486 endif 487 ELSE 488 if(check_status%verboseLevel > 0) then 489 WRITE(*,*) 490 WRITE(*,*) trim(adjustl(varname)), " is NOT IDENTICAL." 491 if(check_status%verboseLevel > 2) then 492 WRITE(*,*) "KERNEL: ", var 493 WRITE(*,*) "REF. : ", ref_var 494 endif 495 endif 496 check_status%numFatal = check_status%numFatal + 1 497 END IF 498 end subroutine 499 500 subroutine verify_var_character(varname, check_status, var, ref_var) 501 character(*), intent(in) :: varname 502 type(check_t), intent(inout) :: check_status 503 character(*), intent(in) :: var, ref_var 504 505 check_status%numTotal = check_status%numTotal + 1 506 IF ( var == ref_var ) THEN 507 check_status%numIdentical = check_status%numIdentical + 1 508 if(check_status%verboseLevel > 1) then 509 WRITE(*,*) 510 WRITE(*,*) trim(adjustl(varname)), " is IDENTICAL( ", var, " )." 511 endif 512 ELSE 513 if(check_status%verboseLevel > 0) then 514 WRITE(*,*) 515 WRITE(*,*) trim(adjustl(varname)), " is NOT IDENTICAL." 516 if(check_status%verboseLevel > 2) then 517 WRITE(*,*) "KERNEL: ", var 518 WRITE(*,*) "REF. : ", ref_var 519 end if 520 end if 521 check_status%numFatal = check_status%numFatal + 1 522 END IF 523 end subroutine 524 525 subroutine verify_var_real_wp_dim2(varname, check_status, var, ref_var) 526 character(*), intent(in) :: varname 527 type(check_t), intent(inout) :: check_status 528 real(kind=wp), intent(in), dimension(:,:) :: var 529 real(kind=wp), intent(in), allocatable, dimension(:,:) :: ref_var 530 real(kind=wp) :: nrmsdiff, rmsdiff 531 real(kind=wp), allocatable :: temp(:,:), temp2(:,:) 532 integer :: n 533 534 535 IF ( ALLOCATED(ref_var) ) THEN 536 check_status%numTotal = check_status%numTotal + 1 537 allocate(temp(SIZE(var,dim=1),SIZE(var,dim=2))) 538 allocate(temp2(SIZE(var,dim=1),SIZE(var,dim=2))) 539 IF ( ALL( var == ref_var ) ) THEN 540 check_status%numIdentical = check_status%numIdentical + 1 541 if(check_status%verboseLevel > 1) then 542 WRITE(*,*) 543 WRITE(*,*) "All elements of ", trim(adjustl(varname)), " are IDENTICAL." 544 !WRITE(*,*) "KERNEL: ", var 545 !WRITE(*,*) "REF. : ", ref_var 546 IF ( ALL( var == 0 ) ) THEN 547 if(check_status%verboseLevel > 2) then 548 WRITE(*,*) "All values are zero." 549 end if 550 END IF 551 end if 552 ELSE 553 n = count(var/=ref_var) 554 where(ref_var .NE. 0) 555 temp = ((var-ref_var)/ref_var)**2 556 temp2 = (var-ref_var)**2 557 elsewhere 558 temp = (var-ref_var)**2 559 temp2 = temp 560 endwhere 561 nrmsdiff = sqrt(sum(temp)/real(n)) 562 rmsdiff = sqrt(sum(temp2)/real(n)) 563 564 if(check_status%verboseLevel > 0) then 565 WRITE(*,*) 566 WRITE(*,*) trim(adjustl(varname)), " is NOT IDENTICAL." 567 WRITE(*,*) count( var /= ref_var), " of ", size( var ), " elements are different." 568 if(check_status%verboseLevel > 1) then 569 WRITE(*,*) "Average - kernel ", sum(var)/real(size(var)) 570 WRITE(*,*) "Average - reference ", sum(ref_var)/real(size(ref_var)) 571 endif 572 WRITE(*,*) "RMS of difference is ",rmsdiff 573 WRITE(*,*) "Normalized RMS of difference is ",nrmsdiff 574 end if 575 576 if (nrmsdiff > check_status%tolerance) then 577 check_status%numFatal = check_status%numFatal+1 578 else 579 check_status%numWarning = check_status%numWarning+1 580 endif 581 END IF 582 deallocate(temp,temp2) 583 END IF 584 end subroutine 585 586 subroutine verify_var_real_wp_dim1(varname, check_status, var, ref_var) 587 character(*), intent(in) :: varname 588 type(check_t), intent(inout) :: check_status 589 real(kind=wp), intent(in), dimension(:) :: var 590 real(kind=wp), intent(in), allocatable, dimension(:) :: ref_var 591 real(kind=wp) :: nrmsdiff, rmsdiff 592 real(kind=wp), allocatable :: temp(:), temp2(:) 593 integer :: n 594 595 596 IF ( ALLOCATED(ref_var) ) THEN 597 check_status%numTotal = check_status%numTotal + 1 598 allocate(temp(SIZE(var,dim=1))) 599 allocate(temp2(SIZE(var,dim=1))) 600 IF ( ALL( var == ref_var ) ) THEN 601 check_status%numIdentical = check_status%numIdentical + 1 602 if(check_status%verboseLevel > 1) then 603 WRITE(*,*) 604 WRITE(*,*) "All elements of ", trim(adjustl(varname)), " are IDENTICAL." 605 !WRITE(*,*) "KERNEL: ", var 606 !WRITE(*,*) "REF. : ", ref_var 607 IF ( ALL( var == 0 ) ) THEN 608 if(check_status%verboseLevel > 2) then 609 WRITE(*,*) "All values are zero." 610 end if 611 END IF 612 end if 613 ELSE 614 n = count(var/=ref_var) 615 where(ref_var .NE. 0) 616 temp = ((var-ref_var)/ref_var)**2 617 temp2 = (var-ref_var)**2 618 elsewhere 619 temp = (var-ref_var)**2 620 temp2 = temp 621 endwhere 622 nrmsdiff = sqrt(sum(temp)/real(n)) 623 rmsdiff = sqrt(sum(temp2)/real(n)) 624 625 if(check_status%verboseLevel > 0) then 626 WRITE(*,*) 627 WRITE(*,*) trim(adjustl(varname)), " is NOT IDENTICAL." 628 WRITE(*,*) count( var /= ref_var), " of ", size( var ), " elements are different." 629 if(check_status%verboseLevel > 1) then 630 WRITE(*,*) "Average - kernel ", sum(var)/real(size(var)) 631 WRITE(*,*) "Average - reference ", sum(ref_var)/real(size(ref_var)) 632 endif 633 WRITE(*,*) "RMS of difference is ",rmsdiff 634 WRITE(*,*) "Normalized RMS of difference is ",nrmsdiff 635 end if 636 637 if (nrmsdiff > check_status%tolerance) then 638 check_status%numFatal = check_status%numFatal+1 639 else 640 check_status%numWarning = check_status%numWarning+1 641 endif 642 END IF 643 deallocate(temp,temp2) 644 END IF 645 end subroutine 646 647 subroutine verify_var_real_wp_dim3(varname, check_status, var, ref_var) 648 character(*), intent(in) :: varname 649 type(check_t), intent(inout) :: check_status 650 real(kind=wp), intent(in), dimension(:,:,:) :: var 651 real(kind=wp), intent(in), allocatable, dimension(:,:,:) :: ref_var 652 real(kind=wp) :: nrmsdiff, rmsdiff 653 real(kind=wp), allocatable :: temp(:,:,:), temp2(:,:,:) 654 integer :: n 655 656 657 IF ( ALLOCATED(ref_var) ) THEN 658 check_status%numTotal = check_status%numTotal + 1 659 allocate(temp(SIZE(var,dim=1),SIZE(var,dim=2),SIZE(var,dim=3))) 660 allocate(temp2(SIZE(var,dim=1),SIZE(var,dim=2),SIZE(var,dim=3))) 661 IF ( ALL( var == ref_var ) ) THEN 662 check_status%numIdentical = check_status%numIdentical + 1 663 if(check_status%verboseLevel > 1) then 664 WRITE(*,*) 665 WRITE(*,*) "All elements of ", trim(adjustl(varname)), " are IDENTICAL." 666 !WRITE(*,*) "KERNEL: ", var 667 !WRITE(*,*) "REF. : ", ref_var 668 IF ( ALL( var == 0 ) ) THEN 669 if(check_status%verboseLevel > 2) then 670 WRITE(*,*) "All values are zero." 671 end if 672 END IF 673 end if 674 ELSE 675 n = count(var/=ref_var) 676 where(ref_var .NE. 0) 677 temp = ((var-ref_var)/ref_var)**2 678 temp2 = (var-ref_var)**2 679 elsewhere 680 temp = (var-ref_var)**2 681 temp2 = temp 682 endwhere 683 nrmsdiff = sqrt(sum(temp)/real(n)) 684 rmsdiff = sqrt(sum(temp2)/real(n)) 685 686 if(check_status%verboseLevel > 0) then 687 WRITE(*,*) 688 WRITE(*,*) trim(adjustl(varname)), " is NOT IDENTICAL." 689 WRITE(*,*) count( var /= ref_var), " of ", size( var ), " elements are different." 690 if(check_status%verboseLevel > 1) then 691 WRITE(*,*) "Average - kernel ", sum(var)/real(size(var)) 692 WRITE(*,*) "Average - reference ", sum(ref_var)/real(size(ref_var)) 693 endif 694 WRITE(*,*) "RMS of difference is ",rmsdiff 695 WRITE(*,*) "Normalized RMS of difference is ",nrmsdiff 696 end if 697 698 if (nrmsdiff > check_status%tolerance) then 699 check_status%numFatal = check_status%numFatal+1 700 else 701 check_status%numWarning = check_status%numWarning+1 702 endif 703 END IF 704 deallocate(temp,temp2) 705 END IF 706 end subroutine 707 708 subroutine verify_var_integer_4_dim2(varname, check_status, var, ref_var) 709 character(*), intent(in) :: varname 710 type(check_t), intent(inout) :: check_status 711 integer(kind=4), intent(in), dimension(:,:) :: var 712 integer(kind=4), intent(in), allocatable, dimension(:,:) :: ref_var 713 integer(kind=4) :: nrmsdiff, rmsdiff 714 integer(kind=4), allocatable :: temp(:,:), temp2(:,:) 715 integer :: n 716 717 718 IF ( ALLOCATED(ref_var) ) THEN 719 check_status%numTotal = check_status%numTotal + 1 720 allocate(temp(SIZE(var,dim=1),SIZE(var,dim=2))) 721 allocate(temp2(SIZE(var,dim=1),SIZE(var,dim=2))) 722 IF ( ALL( var == ref_var ) ) THEN 723 check_status%numIdentical = check_status%numIdentical + 1 724 if(check_status%verboseLevel > 1) then 725 WRITE(*,*) 726 WRITE(*,*) "All elements of ", trim(adjustl(varname)), " are IDENTICAL." 727 !WRITE(*,*) "KERNEL: ", var 728 !WRITE(*,*) "REF. : ", ref_var 729 IF ( ALL( var == 0 ) ) THEN 730 if(check_status%verboseLevel > 2) then 731 WRITE(*,*) "All values are zero." 732 end if 733 END IF 734 end if 735 ELSE 736 n = count(var/=ref_var) 737 where(ref_var .NE. 0) 738 temp = ((var-ref_var)/ref_var)**2 739 temp2 = (var-ref_var)**2 740 elsewhere 741 temp = (var-ref_var)**2 742 temp2 = temp 743 endwhere 744 nrmsdiff = sqrt(sum(temp)/real(n)) 745 rmsdiff = sqrt(sum(temp2)/real(n)) 746 747 if(check_status%verboseLevel > 0) then 748 WRITE(*,*) 749 WRITE(*,*) trim(adjustl(varname)), " is NOT IDENTICAL." 750 WRITE(*,*) count( var /= ref_var), " of ", size( var ), " elements are different." 751 if(check_status%verboseLevel > 1) then 752 WRITE(*,*) "Average - kernel ", sum(var)/real(size(var)) 753 WRITE(*,*) "Average - reference ", sum(ref_var)/real(size(ref_var)) 754 endif 755 WRITE(*,*) "RMS of difference is ",rmsdiff 756 WRITE(*,*) "Normalized RMS of difference is ",nrmsdiff 757 end if 758 759 if (nrmsdiff > check_status%tolerance) then 760 check_status%numFatal = check_status%numFatal+1 761 else 762 check_status%numWarning = check_status%numWarning+1 763 endif 764 END IF 765 deallocate(temp,temp2) 766 END IF 767 end subroutine 768 769 END SUBROUTINE psrad_interface 770 END MODULE mo_psrad_interface 771