1 2! KGEN-generated Fortran source file 3! 4! Filename : prim_advance_mod.F90 5! Generated at: 2015-03-16 09:25:31 6! KGEN version: 0.4.5 7 8 9 10 MODULE prim_advance_mod 11 USE hybvcoord_mod, ONLY : kgen_read_mod5 => kgen_read 12 ! _EXTERNAL 13 IMPLICIT NONE 14 PRIVATE 15 INTEGER, PARAMETER :: kgen_dp = selected_real_kind(15, 307) 16 PUBLIC compute_and_apply_rhs 17 type, public :: check_t 18 logical :: Passed 19 integer :: numFatal 20 integer :: numTotal 21 integer :: numIdentical 22 integer :: numWarning 23 integer :: VerboseLevel 24 real(kind=kgen_dp) :: tolerance 25 end type check_t 26 CONTAINS 27 28 ! write subroutines 29 ! No subroutines 30 ! No module extern variables 31 subroutine kgen_init_check(check,tolerance) 32 type(check_t), intent(inout) :: check 33 real(kind=kgen_dp), intent(in), optional :: tolerance 34 check%Passed = .TRUE. 35 check%numFatal = 0 36 check%numWarning = 0 37 check%numTotal = 0 38 check%numIdentical = 0 39 check%VerboseLevel = 1 40 if(present(tolerance)) then 41 check%tolerance = tolerance 42 else 43 check%tolerance = 1.E-14 44 endif 45 end subroutine kgen_init_check 46 subroutine kgen_print_check(kname, check) 47 character(len=*) :: kname 48 type(check_t), intent(in) :: check 49 write (*,*) 50 write (*,*) TRIM(kname),' KGENPrtCheck: Tolerance for normalized RMS: ',check%tolerance 51 write (*,*) TRIM(kname),' KGENPrtCheck: Number of variables checked: ',check%numTotal 52 write (*,*) TRIM(kname),' KGENPrtCheck: Number of Identical results: ',check%numIdentical 53 write (*,*) TRIM(kname),' KGENPrtCheck: Number of warnings detected: ',check%numWarning 54 write (*,*) TRIM(kname),' KGENPrtCheck: Number of fatal errors detected: ', check%numFatal 55 if (check%numFatal> 0) then 56 write(*,*) TRIM(kname),' KGENPrtCheck: verification FAILED' 57 else 58 write(*,*) TRIM(kname),' KGENPrtCheck: verification PASSED' 59 endif 60 end subroutine kgen_print_check 61 62 63 64 65 66 67 68 69 70 ! 71 ! phl notes: output is stored in first argument. Advances from 2nd argument using tendencies evaluated at 3rd rgument: 72 ! phl: for offline winds use time at 3rd argument (same as rhs currently) 73 ! 74 75 SUBROUTINE compute_and_apply_rhs(hvcoord, kgen_unit) 76 ! =================================== 77 ! compute the RHS, accumulate into u(np1) and apply DSS 78 ! 79 ! u(np1) = u(nm1) + dt2*DSS[ RHS(u(n0)) ] 80 ! 81 ! This subroutine is normally called to compute a leapfrog timestep 82 ! but by adjusting np1,nm1,n0 and dt2, many other timesteps can be 83 ! accomodated. For example, setting nm1=np1=n0 this routine will 84 ! take a forward euler step, overwriting the input with the output. 85 ! 86 ! qn0 = timelevel used to access Qdp() in order to compute virtual Temperature 87 ! qn0=-1 for the dry case 88 ! 89 ! if dt2<0, then the DSS'd RHS is returned in timelevel np1 90 ! 91 ! Combining the RHS and DSS pack operation in one routine 92 ! allows us to fuse these two loops for more cache reuse 93 ! 94 ! Combining the dt advance and DSS unpack operation in one routine 95 ! allows us to fuse these two loops for more cache reuse 96 ! 97 ! note: for prescribed velocity case, velocity will be computed at 98 ! "real_time", which should be the time of timelevel n0. 99 ! 100 ! 101 ! =================================== 102 USE kinds, ONLY: real_kind 103 USE dimensions_mod, ONLY: np 104 USE dimensions_mod, ONLY: nlev 105 USE hybvcoord_mod, ONLY: hvcoord_t 106 USE prim_si_mod, ONLY: preq_omega_ps 107 IMPLICIT NONE 108 integer, intent(in) :: kgen_unit 109 INTEGER*8 :: kgen_intvar, start_clock, stop_clock, rate_clock 110 TYPE(check_t):: check_status 111 REAL(KIND=kgen_dp) :: tolerance 112 TYPE(hvcoord_t), intent(in) :: hvcoord 113 ! weighting for eta_dot_dpdn mean flux 114 ! local 115 ! surface pressure for current tiime level 116 REAL(KIND=real_kind), dimension(np,np,nlev) :: omega_p 117 REAL(KIND=real_kind) :: ref_omega_p(np,np,nlev) 118 REAL(KIND=real_kind), dimension(np,np,nlev) :: divdp 119 ! half level vertical velocity on p-grid 120 ! temporary field 121 ! generic gradient storage 122 ! 123 ! 124 ! v.grad(T) 125 ! kinetic energy + PHI term 126 ! lat-lon coord version 127 ! gradient(p - p_met) 128 ! vorticity 129 REAL(KIND=real_kind), dimension(np,np,nlev) :: p ! pressure 130 ! delta pressure 131 ! inverse of delta pressure 132 ! temperature vertical advection 133 REAL(KIND=real_kind), dimension(np,np,nlev) :: vgrad_p ! v.grad(p) 134 ! half level pressures on p-grid 135 ! velocity vertical advection 136 !JMD call t_barrierf('sync_compute_and_apply_rhs', hybrid%par%comm) 137 tolerance = 1.E-14 138 CALL kgen_init_check(check_status, tolerance) 139 READ(UNIT=kgen_unit) omega_p 140 READ(UNIT=kgen_unit) divdp 141 READ(UNIT=kgen_unit) p 142 READ(UNIT=kgen_unit) vgrad_p 143 144 READ(UNIT=kgen_unit) ref_omega_p 145 146 ! call to kernel 147 CALL preq_omega_ps(omega_p, hvcoord, p, vgrad_p, divdp) 148 ! kernel verification for output variables 149 CALL kgen_verify_real_real_kind_dim3( "omega_p", check_status, omega_p, ref_omega_p) 150 CALL kgen_print_check("preq_omega_ps", check_status) 151 CALL system_clock(start_clock, rate_clock) 152 DO kgen_intvar=1,10 153 CALL preq_omega_ps(omega_p, hvcoord, p, vgrad_p, divdp) 154 END DO 155 CALL system_clock(stop_clock, rate_clock) 156 WRITE(*,*) 157 PRINT *, "Elapsed time (sec): ", (stop_clock - start_clock)/REAL(rate_clock*10) 158 ! ============================================================= 159 ! Insert communications here: for shared memory, just a single 160 ! sync is required 161 ! ============================================================= 162 CONTAINS 163 164 ! write subroutines 165 SUBROUTINE kgen_read_real_real_kind_dim3(var, kgen_unit) 166 INTEGER, INTENT(IN) :: kgen_unit 167 real(KIND=real_kind), INTENT(OUT), ALLOCATABLE, DIMENSION(:,:,:) :: var 168 LOGICAL :: is_true 169 INTEGER :: idx1,idx2,idx3 170 INTEGER, DIMENSION(2,3) :: kgen_bound 171 172 READ(UNIT = kgen_unit) is_true 173 174 IF ( is_true ) THEN 175 READ(UNIT = kgen_unit) kgen_bound(1, 1) 176 READ(UNIT = kgen_unit) kgen_bound(2, 1) 177 READ(UNIT = kgen_unit) kgen_bound(1, 2) 178 READ(UNIT = kgen_unit) kgen_bound(2, 2) 179 READ(UNIT = kgen_unit) kgen_bound(1, 3) 180 READ(UNIT = kgen_unit) kgen_bound(2, 3) 181 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)) 182 READ(UNIT = kgen_unit) var 183 END IF 184 END SUBROUTINE kgen_read_real_real_kind_dim3 185 186 187 subroutine kgen_verify_logical(varname, check_status, var, ref_var) 188 character(*), intent(in) :: varname 189 type(check_t), intent(inout) :: check_status 190 logical, intent(in) :: var, ref_var 191 192 check_status%numTotal = check_status%numTotal + 1 193 IF ( var .eqv. ref_var ) THEN 194 check_status%numIdentical = check_status%numIdentical + 1 195 if(check_status%verboseLevel > 1) then 196 WRITE(*,*) 197 WRITE(*,*) trim(adjustl(varname)), " is IDENTICAL( ", var, " )." 198 endif 199 ELSE 200 if(check_status%verboseLevel > 1) then 201 WRITE(*,*) 202 WRITE(*,*) trim(adjustl(varname)), " is NOT IDENTICAL." 203 if(check_status%verboseLevel > 2) then 204 WRITE(*,*) "KERNEL: ", var 205 WRITE(*,*) "REF. : ", ref_var 206 endif 207 endif 208 check_status%numFatal = check_status%numFatal + 1 209 END IF 210 end subroutine 211 212 subroutine kgen_verify_integer(varname, check_status, var, ref_var) 213 character(*), intent(in) :: varname 214 type(check_t), intent(inout) :: check_status 215 integer, intent(in) :: var, ref_var 216 217 check_status%numTotal = check_status%numTotal + 1 218 IF ( var == ref_var ) THEN 219 check_status%numIdentical = check_status%numIdentical + 1 220 if(check_status%verboseLevel > 1) then 221 WRITE(*,*) 222 WRITE(*,*) trim(adjustl(varname)), " is IDENTICAL( ", var, " )." 223 endif 224 ELSE 225 if(check_status%verboseLevel > 0) then 226 WRITE(*,*) 227 WRITE(*,*) trim(adjustl(varname)), " is NOT IDENTICAL." 228 if(check_status%verboseLevel > 2) then 229 WRITE(*,*) "KERNEL: ", var 230 WRITE(*,*) "REF. : ", ref_var 231 endif 232 endif 233 check_status%numFatal = check_status%numFatal + 1 234 END IF 235 end subroutine 236 237 subroutine kgen_verify_real(varname, check_status, var, ref_var) 238 character(*), intent(in) :: varname 239 type(check_t), intent(inout) :: check_status 240 real, intent(in) :: var, ref_var 241 242 check_status%numTotal = check_status%numTotal + 1 243 IF ( var == ref_var ) THEN 244 check_status%numIdentical = check_status%numIdentical + 1 245 if(check_status%verboseLevel > 1) then 246 WRITE(*,*) 247 WRITE(*,*) trim(adjustl(varname)), " is IDENTICAL( ", var, " )." 248 endif 249 ELSE 250 if(check_status%verboseLevel > 0) then 251 WRITE(*,*) 252 WRITE(*,*) trim(adjustl(varname)), " is NOT IDENTICAL." 253 if(check_status%verboseLevel > 2) then 254 WRITE(*,*) "KERNEL: ", var 255 WRITE(*,*) "REF. : ", ref_var 256 endif 257 endif 258 check_status%numFatal = check_status%numFatal + 1 259 END IF 260 end subroutine 261 262 subroutine kgen_verify_character(varname, check_status, var, ref_var) 263 character(*), intent(in) :: varname 264 type(check_t), intent(inout) :: check_status 265 character(*), intent(in) :: var, ref_var 266 267 check_status%numTotal = check_status%numTotal + 1 268 IF ( var == ref_var ) THEN 269 check_status%numIdentical = check_status%numIdentical + 1 270 if(check_status%verboseLevel > 1) then 271 WRITE(*,*) 272 WRITE(*,*) trim(adjustl(varname)), " is IDENTICAL( ", var, " )." 273 endif 274 ELSE 275 if(check_status%verboseLevel > 0) then 276 WRITE(*,*) 277 WRITE(*,*) trim(adjustl(varname)), " is NOT IDENTICAL." 278 if(check_status%verboseLevel > 2) then 279 WRITE(*,*) "KERNEL: ", var 280 WRITE(*,*) "REF. : ", ref_var 281 end if 282 end if 283 check_status%numFatal = check_status%numFatal + 1 284 END IF 285 end subroutine 286 287 subroutine kgen_verify_real_real_kind_dim3(varname, check_status, var, ref_var) 288 character(*), intent(in) :: varname 289 type(check_t), intent(inout) :: check_status 290 real(kind=real_kind), intent(in), dimension(:,:,:) :: var, ref_var 291 !real(kind=real_kind), intent(in), dimension(:,:,:) :: ref_var 292 real(kind=real_kind) :: nrmsdiff, rmsdiff 293 real(kind=real_kind), allocatable :: temp(:,:,:), temp2(:,:,:) 294 integer :: n 295 296 297 IF ( .TRUE. ) THEN 298 check_status%numTotal = check_status%numTotal + 1 299 allocate(temp(SIZE(var,dim=1),SIZE(var,dim=2),SIZE(var,dim=3))) 300 allocate(temp2(SIZE(var,dim=1),SIZE(var,dim=2),SIZE(var,dim=3))) 301 IF ( ALL( var == ref_var ) ) THEN 302 check_status%numIdentical = check_status%numIdentical + 1 303 if(check_status%verboseLevel > 1) then 304 WRITE(*,*) 305 WRITE(*,*) "All elements of ", trim(adjustl(varname)), " are IDENTICAL." 306 !WRITE(*,*) "KERNEL: ", var 307 !WRITE(*,*) "REF. : ", ref_var 308 IF ( ALL( var == 0 ) ) THEN 309 if(check_status%verboseLevel > 2) then 310 WRITE(*,*) "All values are zero." 311 end if 312 END IF 313 end if 314 ELSE 315 n = count(var/=ref_var) 316 where(ref_var .NE. 0) 317 temp = ((var-ref_var)/ref_var)**2 318 temp2 = (var-ref_var)**2 319 elsewhere 320 temp = (var-ref_var)**2 321 temp2 = temp 322 endwhere 323 nrmsdiff = sqrt(sum(temp)/real(n)) 324 rmsdiff = sqrt(sum(temp2)/real(n)) 325 326 if(check_status%verboseLevel > 0) then 327 WRITE(*,*) 328 WRITE(*,*) trim(adjustl(varname)), " is NOT IDENTICAL." 329 WRITE(*,*) count( var /= ref_var), " of ", size( var ), " elements are different." 330 if(check_status%verboseLevel > 1) then 331 WRITE(*,*) "Average - kernel ", sum(var)/real(size(var)) 332 WRITE(*,*) "Average - reference ", sum(ref_var)/real(size(ref_var)) 333 endif 334 WRITE(*,*) "RMS of difference is ",rmsdiff 335 WRITE(*,*) "Normalized RMS of difference is ",nrmsdiff 336 end if 337 338 if (nrmsdiff > check_status%tolerance) then 339 check_status%numFatal = check_status%numFatal+1 340 else 341 check_status%numWarning = check_status%numWarning+1 342 endif 343 END IF 344 deallocate(temp,temp2) 345 END IF 346 end subroutine 347 348 END SUBROUTINE compute_and_apply_rhs 349 350 !TRILINOS 351 352 353 END MODULE prim_advance_mod 354