1
2! KGEN-generated Fortran source file
3!
4! Filename    : prim_advance_mod.F90
5! Generated at: 2015-04-12 19:37:49
6! KGEN version: 0.4.9
7
8
9
10    MODULE prim_advance_mod
11        USE kgen_utils_mod, ONLY : kgen_dp, check_t, kgen_init_check, kgen_print_check
12    USE element_mod, ONLY : kgen_read_mod9 => kgen_read
13    USE element_mod, ONLY : kgen_verify_mod9 => kgen_verify
14        ! _EXTERNAL
15        IMPLICIT NONE
16        PRIVATE
17        PUBLIC compute_and_apply_rhs
18        CONTAINS
19
20        ! write subroutines
21        ! No subroutines
22        ! No module extern variables
23
24
25
26
27
28
29
30
31
32        !
33        ! phl notes: output is stored in first argument. Advances from 2nd argument using tendencies evaluated at 3rd rgument:
34        ! phl: for offline winds use time at 3rd argument (same as rhs currently)
35        !
36
37        SUBROUTINE compute_and_apply_rhs(elem, kgen_unit)
38                USE kgen_utils_mod, ONLY : kgen_dp, check_t, kgen_init_check, kgen_print_check
39            ! ===================================
40            ! compute the RHS, accumulate into u(np1) and apply DSS
41            !
42            !           u(np1) = u(nm1) + dt2*DSS[ RHS(u(n0)) ]
43            !
44            ! This subroutine is normally called to compute a leapfrog timestep
45            ! but by adjusting np1,nm1,n0 and dt2, many other timesteps can be
46            ! accomodated.  For example, setting nm1=np1=n0 this routine will
47            ! take a forward euler step, overwriting the input with the output.
48            !
49            !    qn0 = timelevel used to access Qdp() in order to compute virtual Temperature
50            !          qn0=-1 for the dry case
51            !
52            ! if  dt2<0, then the DSS'd RHS is returned in timelevel np1
53            !
54            ! Combining the RHS and DSS pack operation in one routine
55            ! allows us to fuse these two loops for more cache reuse
56            !
57            ! Combining the dt advance and DSS unpack operation in one routine
58            ! allows us to fuse these two loops for more cache reuse
59            !
60            ! note: for prescribed velocity case, velocity will be computed at
61            ! "real_time", which should be the time of timelevel n0.
62            !
63            !
64            ! ===================================
65            USE kinds, ONLY: real_kind
66            USE dimensions_mod, ONLY: np
67            USE dimensions_mod, ONLY: nlev
68            USE element_mod, ONLY: element_t
69            USE prim_si_mod, ONLY: preq_hydrostatic
70            IMPLICIT NONE
71            integer, intent(in) :: kgen_unit
72            INTEGER*8 :: kgen_intvar, start_clock, stop_clock, rate_clock
73            TYPE(check_t):: check_status
74            REAL(KIND=kgen_dp) :: tolerance
75            TYPE(element_t), intent(inout), target :: elem(:)
76            ! weighting for eta_dot_dpdn mean flux
77            ! local
78            ! surface pressure for current tiime level
79            REAL(KIND=real_kind), pointer, dimension(:,:,:) :: phi
80            REAL(KIND=real_kind), pointer :: ref_phi(:,:,:) => NULL()
81            REAL(KIND=real_kind), dimension(np,np,nlev) :: t_v
82            ! half level vertical velocity on p-grid
83            ! temporary field
84            ! generic gradient storage
85            ! v.grad(T)
86            ! kinetic energy + PHI term
87            ! lat-lon coord version
88            ! vorticity
89            REAL(KIND=real_kind), dimension(np,np,nlev) :: p ! pressure
90            REAL(KIND=real_kind), dimension(np,np,nlev) :: dp ! delta pressure
91            ! inverse of delta pressure
92            ! temperature vertical advection
93            ! v.grad(p)
94            ! half level pressures on p-grid
95            ! velocity vertical advection
96            INTEGER :: ie
97            !JMD  call t_barrierf('sync_compute_and_apply_rhs', hybrid%par%comm)
98                tolerance = 1.E-14
99                CALL kgen_init_check(check_status, tolerance)
100                CALL kgen_read_real_real_kind_dim3_ptr(phi, kgen_unit)
101                READ(UNIT=kgen_unit) t_v
102                READ(UNIT=kgen_unit) p
103                READ(UNIT=kgen_unit) dp
104                READ(UNIT=kgen_unit) ie
105
106                CALL kgen_read_real_real_kind_dim3_ptr(ref_phi, kgen_unit)
107
108
109                ! call to kernel
110                CALL preq_hydrostatic(phi, elem(ie)%state%phis, t_v, p, dp)
111                ! kernel verification for output variables
112                CALL kgen_verify_real_real_kind_dim3_ptr( "phi", check_status, phi, ref_phi)
113                CALL kgen_print_check("preq_hydrostatic", check_status)
114                CALL system_clock(start_clock, rate_clock)
115                DO kgen_intvar=1,10
116                    CALL preq_hydrostatic(phi, elem(ie) % state % phis, t_v, p, dp)
117                END DO
118                CALL system_clock(stop_clock, rate_clock)
119                WRITE(*,*)
120                PRINT *, "Elapsed time (sec): ", (stop_clock - start_clock)/REAL(rate_clock*10)
121            ! =============================================================
122            ! Insert communications here: for shared memory, just a single
123            ! sync is required
124            ! =============================================================
125        CONTAINS
126
127        ! write subroutines
128            SUBROUTINE kgen_read_real_real_kind_dim3_ptr(var, kgen_unit, printvar)
129                INTEGER, INTENT(IN) :: kgen_unit
130                CHARACTER(*), INTENT(IN), OPTIONAL :: printvar
131                real(KIND=real_kind), INTENT(OUT), POINTER, DIMENSION(:,:,:) :: var
132                LOGICAL :: is_true
133                INTEGER :: idx1,idx2,idx3
134                INTEGER, DIMENSION(2,3) :: kgen_bound
135
136                READ(UNIT = kgen_unit) is_true
137
138                IF ( is_true ) THEN
139                    READ(UNIT = kgen_unit) kgen_bound(1, 1)
140                    READ(UNIT = kgen_unit) kgen_bound(2, 1)
141                    READ(UNIT = kgen_unit) kgen_bound(1, 2)
142                    READ(UNIT = kgen_unit) kgen_bound(2, 2)
143                    READ(UNIT = kgen_unit) kgen_bound(1, 3)
144                    READ(UNIT = kgen_unit) kgen_bound(2, 3)
145                    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))
146                    READ(UNIT = kgen_unit) var
147                    IF ( PRESENT(printvar) ) THEN
148                        PRINT *, "** " // printvar // " **", var
149                    END IF
150                END IF
151            END SUBROUTINE kgen_read_real_real_kind_dim3_ptr
152
153
154        ! verify subroutines
155            SUBROUTINE kgen_verify_real_real_kind_dim3_ptr( varname, check_status, var, ref_var)
156                character(*), intent(in) :: varname
157                type(check_t), intent(inout) :: check_status
158                real(KIND=real_kind), intent(in), DIMENSION(:,:,:), POINTER :: var, ref_var
159                real(KIND=real_kind) :: nrmsdiff, rmsdiff
160                real(KIND=real_kind), allocatable, DIMENSION(:,:,:) :: temp, temp2
161                integer :: n
162                IF ( ASSOCIATED(var) ) THEN
163                check_status%numTotal = check_status%numTotal + 1
164                IF ( ALL( var == ref_var ) ) THEN
165
166                    check_status%numIdentical = check_status%numIdentical + 1
167                    if(check_status%verboseLevel > 1) then
168                        WRITE(*,*)
169                        WRITE(*,*) "All elements of ", trim(adjustl(varname)), " are IDENTICAL."
170                        !WRITE(*,*) "KERNEL: ", var
171                        !WRITE(*,*) "REF.  : ", ref_var
172                        IF ( ALL( var == 0 ) ) THEN
173                            if(check_status%verboseLevel > 2) then
174                                WRITE(*,*) "All values are zero."
175                            end if
176                        END IF
177                    end if
178                ELSE
179                    allocate(temp(SIZE(var,dim=1),SIZE(var,dim=2),SIZE(var,dim=3)))
180                    allocate(temp2(SIZE(var,dim=1),SIZE(var,dim=2),SIZE(var,dim=3)))
181
182                    n = count(var/=ref_var)
183                    where(abs(ref_var) > check_status%minvalue)
184                        temp  = ((var-ref_var)/ref_var)**2
185                        temp2 = (var-ref_var)**2
186                    elsewhere
187                        temp  = (var-ref_var)**2
188                        temp2 = temp
189                    endwhere
190                    nrmsdiff = sqrt(sum(temp)/real(n))
191                    rmsdiff = sqrt(sum(temp2)/real(n))
192
193                    if(check_status%verboseLevel > 0) then
194                        WRITE(*,*)
195                        WRITE(*,*) trim(adjustl(varname)), " is NOT IDENTICAL."
196                        WRITE(*,*) count( var /= ref_var), " of ", size( var ), " elements are different."
197                        if(check_status%verboseLevel > 1) then
198                            WRITE(*,*) "Average - kernel ", sum(var)/real(size(var))
199                            WRITE(*,*) "Average - reference ", sum(ref_var)/real(size(ref_var))
200                        endif
201                        WRITE(*,*) "RMS of difference is ",rmsdiff
202                        WRITE(*,*) "Normalized RMS of difference is ",nrmsdiff
203                    end if
204
205                    if (nrmsdiff > check_status%tolerance) then
206                        check_status%numFatal = check_status%numFatal+1
207                    else
208                        check_status%numWarning = check_status%numWarning+1
209                    endif
210
211                    deallocate(temp,temp2)
212                END IF
213                END IF
214            END SUBROUTINE kgen_verify_real_real_kind_dim3_ptr
215
216        END SUBROUTINE compute_and_apply_rhs
217        !TRILINOS
218
219
220    END MODULE prim_advance_mod
221