summaryrefslogtreecommitdiff
path: root/benchmarks/CUDA/WP
diff options
context:
space:
mode:
authorTor Aamodt <[email protected]>2010-07-15 18:09:46 -0800
committerTor Aamodt <[email protected]>2010-07-15 18:09:46 -0800
commit69f2911e04ffb1b19eef1fafb8c040af271f656e (patch)
tree231d3b6bdc3a202f7c255bfcf7bf2c36e32cee9e /benchmarks/CUDA/WP
creating branch for adding support for CUDA 3.x and Fermi
[git-p4: depot-paths = "//depot/gpgpu_sim_research/fermi/distribution/": change = 6829]
Diffstat (limited to 'benchmarks/CUDA/WP')
-rw-r--r--benchmarks/CUDA/WP/README62
-rw-r--r--benchmarks/CUDA/WP/README.GPGPU-Sim4
-rw-r--r--benchmarks/CUDA/WP/SAMPLE_SESSION192
-rw-r--r--benchmarks/CUDA/WP/body_inline.h741
-rw-r--r--benchmarks/CUDA/WP/compare_snaps.F83
-rw-r--r--benchmarks/CUDA/WP/cublas.h2
-rw-r--r--benchmarks/CUDA/WP/d2s.h36
-rw-r--r--benchmarks/CUDA/WP/data/wsm5_in_010bin0 -> 4530080 bytes
-rw-r--r--benchmarks/CUDA/WP/debug.m435
-rw-r--r--benchmarks/CUDA/WP/libmassv.F385
-rw-r--r--benchmarks/CUDA/WP/makefile138
-rw-r--r--benchmarks/CUDA/WP/microclock.c33
-rw-r--r--benchmarks/CUDA/WP/module_mp_wsm5.F1731
-rw-r--r--benchmarks/CUDA/WP/s2d.h36
-rw-r--r--benchmarks/CUDA/WP/spt.h150
-rwxr-xr-xbenchmarks/CUDA/WP/spt.pl219
-rw-r--r--benchmarks/CUDA/WP/timings_2007120524
-rw-r--r--benchmarks/CUDA/WP/util.h265
-rw-r--r--benchmarks/CUDA/WP/util4.cu46
-rw-r--r--benchmarks/CUDA/WP/wsm5.cu500
-rw-r--r--benchmarks/CUDA/WP/wsm5_constants.h92
-rw-r--r--benchmarks/CUDA/WP/wsm5_gpu.cu783
22 files changed, 5557 insertions, 0 deletions
diff --git a/benchmarks/CUDA/WP/README b/benchmarks/CUDA/WP/README
new file mode 100644
index 0000000..04213b1
--- /dev/null
+++ b/benchmarks/CUDA/WP/README
@@ -0,0 +1,62 @@
+WSM5 (WRF microphysics) test driver for host CPU and GPU, complete with test data sets.
+
+To compile:
+ Edit makefile for your system.
+ (Variables XXX and YYY control number of threads per block. See settings
+ in makefile for settings for different GPUs)
+ make
+ (make clean also works)
+
+ Makes 2 execs: wsm5_driver_vanilla (host cpu) and wsm5_driver_chocolate (gpu)
+
+To run:
+ CD to run and type runit
+
+Contents:
+
+This directory:
+ Source & Build:
+ makefile Builds CUDA and Original versions (also CUDA Emulated)
+ module_mp_wsm5.F Contains standalone driver and orig WSM5 source
+ wsm5.cu WSM5 CUDA code, Host side
+ wsm5_gpu.cu WSM5 CUDA code, GPU side
+ libmassv.F MASS compatibility lib
+ SAMPLE_SESSION Commands and their output for build and run for GTX 5600 GPU
+ using intel compilers on qp.ncsa.uiuc.edu
+
+ Includes:
+ wsm5_constants.h Constants (give good output agreement Intel/Gfortran)
+ s2d.h Copy file prec. to code prec. (might be same)
+ d2s.h Copy from code prec. to driver prec. (ditto)
+ debug.m4 M4 macros for debug out from the GPU code if enabled
+ spt.h Header (supports my directives)
+
+ Utilities:
+ compare_snaps.F Comp CUDA v HOST output and gen difference plot data
+ microclock.c Microsecond timer
+ spt.pl Perl script
+
+Subdirectory: run
+ match.c Compare 2 floats and report # of places of agreement.
+ runit Script to run codes and compare
+ wsm5_in_001 Ten gracious sets of input and output
+ wsm5_in_002
+ wsm5_in_003
+ wsm5_in_004
+ wsm5_in_005
+ wsm5_in_006
+ wsm5_in_007
+ wsm5_in_008
+ wsm5_in_009
+ wsm5_in_010
+ wsm5_out_001
+ wsm5_out_002
+ wsm5_out_003
+ wsm5_out_004
+ wsm5_out_005
+ wsm5_out_006
+ wsm5_out_007
+ wsm5_out_008
+ wsm5_out_009
+ wsm5_out_010
+
diff --git a/benchmarks/CUDA/WP/README.GPGPU-Sim b/benchmarks/CUDA/WP/README.GPGPU-Sim
new file mode 100644
index 0000000..d1d223c
--- /dev/null
+++ b/benchmarks/CUDA/WP/README.GPGPU-Sim
@@ -0,0 +1,4 @@
+ensure GPGPUSIM_ROOT/lib is in LD_LIBRARY_PATH
+
+make
+echo "10 ./data/" | ./gpgpu_ptx_sim__wsm5_standalone
diff --git a/benchmarks/CUDA/WP/SAMPLE_SESSION b/benchmarks/CUDA/WP/SAMPLE_SESSION
new file mode 100644
index 0000000..bf3f8d5
--- /dev/null
+++ b/benchmarks/CUDA/WP/SAMPLE_SESSION
@@ -0,0 +1,192 @@
+Script started on Mon 01 Sep 2008 02:10:27 PM CDT
+$ tcsh
+$ setenv XXX 8
+$ setenv YYY 8
+$ make clean
+\rm -f *.o *.cu.c x.cu y.cu wsm5_driver_* *.mod *.f90 x.ptx
+rm: No match.
+make: *** [clean] Error 1
+$ make
+/lib/cpp -C -P libmassv.F > libmassv.f90
+ifort -c -O3 -align all -FR -convert big_endian -fno-alias -fno-fnalias libmassv.f90
+icc -c -w -O3 -ip -c -DMKX=28 microclock.c
+/lib/cpp -C -P -DDEBUG_I=59 -DDEBUG_J=45 -DDEBUG_K=1 -DRWORDSIZE=4 module_mp_wsm5.F > module_mp_wsm5.f90
+ifort -c -O3 -align all -FR -convert big_endian -fno-alias -fno-fnalias module_mp_wsm5.f90
+fortcom: Info: module_mp_wsm5.f90, line 382: This statement function has not been used. [CONDEN]
+ viscos, xka, venfac, conden, diffac, &
+---------------------------------^
+ifort -o wsm5_driver_vanilla module_mp_wsm5.o libmassv.o microclock.o /usr/local/cuda/lib/libcublasemu.so
+m4 wsm5.cu | sed "s/float/float/g" > y.cu
+nvcc -DCUDA -ccbin /usr/bin -cuda -DDEBUG_I=59 -DDEBUG_J=45 -DDEBUG_K=1 -DXXX=8 -DYYY=8 -DMKX=28 --host-compilation 'C++' --use_fast_math y.cu
+"y.cu", line 146: warning: variable "m" was declared but never referenced
+ int i, dc, m ;
+ ^
+
+"y.cu", line 324: warning: variable "k" was declared but never referenced
+ { int k ;
+ ^
+
+"y.cu", line 216: warning: variable "i" was declared but never referenced
+ int i, j, k ;
+ ^
+
+"y.cu", line 216: warning: variable "j" was declared but never referenced
+ int i, j, k ;
+ ^
+
+"y.cu", line 216: warning: variable "k" was declared but never referenced
+ int i, j, k ;
+ ^
+
+"y.cu", line 217: warning: variable "bigbuf" was set but never used
+ float *bigbuf ;
+ ^
+
+"y.cu", line 218: warning: variable "s1" was declared but never referenced
+ int s, e, s1, e1, s2, e2 ;
+ ^
+
+"y.cu", line 218: warning: variable "e1" was declared but never referenced
+ int s, e, s1, e1, s2, e2 ;
+ ^
+
+"y.cu", line 473: warning: missing return statement at end of non-void
+ function "get_wsm5_gpu_levels_"
+ }
+ ^
+
+"y.cu", line 146: warning: variable "m" was declared but never referenced
+ int i, dc, m ;
+ ^
+
+"y.cu", line 324: warning: variable "k" was declared but never referenced
+ { int k ;
+ ^
+
+"y.cu", line 216: warning: variable "i" was declared but never referenced
+ int i, j, k ;
+ ^
+
+"y.cu", line 216: warning: variable "j" was declared but never referenced
+ int i, j, k ;
+ ^
+
+"y.cu", line 216: warning: variable "k" was declared but never referenced
+ int i, j, k ;
+ ^
+
+"y.cu", line 217: warning: variable "bigbuf" was set but never used
+ float *bigbuf ;
+ ^
+
+"y.cu", line 218: warning: variable "s1" was declared but never referenced
+ int s, e, s1, e1, s2, e2 ;
+ ^
+
+"y.cu", line 218: warning: variable "e1" was declared but never referenced
+ int s, e, s1, e1, s2, e2 ;
+ ^
+
+"y.cu", line 473: warning: missing return statement at end of non-void
+ function "get_wsm5_gpu_levels_"
+ }
+ ^
+
+/bin/mv y.cu.cpp wsm5.cu.cpp
+icc -w -O3 -ip -c -c wsm5.cu.cpp
+m4 wsm5_gpu.cu | spt.pl | sed "s/float/float/g" > x.cu
+nvcc -DCUDA -ccbin /usr/bin -cuda -DDEBUG_I=59 -DDEBUG_J=45 -DDEBUG_K=1 -DXXX=8 -DYYY=8 -DMKX=28 --host-compilation 'C++' --use_fast_math x.cu
+"x.cu", line 364: warning: variable "fall1_kp1" was set but never used
+ float den_k, falk1_k, falk1_kp1, fall1_k, fall1_kp1, delz_k, delz_kp1 ;
+ ^
+
+"x.cu", line 365: warning: variable "fall2_kp1" was set but never used
+ float falk2_k, falk2_kp1, fall2_k, fall2_kp1 ;
+ ^
+
+"x.cu", line 57: warning: variable "xni_k" was declared but never referenced
+ float t_k, q_k, qr_k, qc_k, qs_k, qi_k, qs1_k, qs2_k, cpm_k, xl_k, xni_k, w1_k, w2_k, w3_k ;
+ ^
+
+"util.h", line 111: warning: type qualifier on return type is meaningless
+ extern __attribute__((__device__)) const struct Float4 _Zplf6Float4(const float, const struct Float4);
+ ^
+
+"util.h", line 120: warning: type qualifier on return type is meaningless
+ extern __attribute__((__device__)) const struct Float4 _Zmlf6Float4(const float, const struct Float4);
+ ^
+
+"util.h", line 129: warning: type qualifier on return type is meaningless
+ extern __attribute__((__device__)) const struct Float4 _Zdvf6Float4(const float, const struct Float4);
+ ^
+
+"util.h", line 138: warning: type qualifier on return type is meaningless
+ extern __attribute__((__device__)) const struct Float4 _Zmif6Float4(const float, const struct Float4);
+ ^
+
+"util.h", line 111: warning: type qualifier on return type is meaningless
+ __attribute__((__device__)) const struct Float4 _Zplf6Float4(
+ ^
+
+"util.h", line 120: warning: type qualifier on return type is meaningless
+ __attribute__((__device__)) const struct Float4 _Zmlf6Float4(
+ ^
+
+"util.h", line 129: warning: type qualifier on return type is meaningless
+ __attribute__((__device__)) const struct Float4 _Zdvf6Float4(
+ ^
+
+"util.h", line 138: warning: type qualifier on return type is meaningless
+ __attribute__((__device__)) const struct Float4 _Zmif6Float4(
+ ^
+
+"x.cu", line 364: warning: variable "fall1_kp1" was set but never used
+ float den_k, falk1_k, falk1_kp1, fall1_k, fall1_kp1, delz_k, delz_kp1 ;
+ ^
+
+"x.cu", line 365: warning: variable "fall2_kp1" was set but never used
+ float falk2_k, falk2_kp1, fall2_k, fall2_kp1 ;
+ ^
+
+"x.cu", line 57: warning: variable "xni_k" was declared but never referenced
+ float t_k, q_k, qr_k, qc_k, qs_k, qi_k, qs1_k, qs2_k, cpm_k, xl_k, xni_k, w1_k, w2_k, w3_k ;
+ ^
+
+/bin/mv x.cu.cpp wsm5_gpu.cu.cpp
+icc -w -O3 -ip -c -c wsm5_gpu.cu.cpp
+/lib/cpp -C -P -DDEBUG_I=59 -DDEBUG_J=45 -DDEBUG_K=1 -DRUN_ON_GPU -DRWORDSIZE=4 module_mp_wsm5.F > module_mp_wsm5.f90
+ifort -c -O3 -align all -FR -convert big_endian -fno-alias -fno-fnalias module_mp_wsm5.f90
+fortcom: Info: module_mp_wsm5.f90, line 364: This statement function has not been used. [CONDEN]
+ viscos, xka, venfac, conden, diffac, &
+---------------------------------^
+ifort -o wsm5_driver_chocolate module_mp_wsm5.o wsm5.cu.o wsm5_gpu.cu.o libmassv.o microclock.o /usr/local/cuda/lib/libcublas.so
+$ cd run
+$ runit
+vanilla
+chocolate
+vanilla: Call to WSM5 on host for step 10 is 195998 microseconds
+chocolate:Call to wsm5_gpu (not including data xfer): 12083 microseconds
+chocolate:Call to wsm5_gpu (including data xfer): 22516 microseconds
+chocolate: Call to WSM5 on host for step 10 is 35985 microseconds
+running compare_snaps
+ 1 74 1 61 1 28
+ -3 78 -3 65 1 28
+ 2 72 2 59 1 27
+ 1 74 1 61 1 28
+ -3 78 -3 65 1 28
+ 2 72 2 59 1 27
+fort.88 fort.98 differ: byte 376, line 23
+$ paste fort.88 fort.98 | head
+ 71 58 th 71 58 th
+ 286.2785 286.2785
+ 286.2001 286.2001
+ 286.0901 286.0901
+ 285.9821 285.9821
+ 285.9263 285.9263
+ 285.9068 285.9068
+ 285.8943 285.8943
+ 285.9103 285.9103
+ 285.9970 285.9970
+$ exit
+
+Script done on Mon 01 Sep 2008 02:16:58 PM CDT
diff --git a/benchmarks/CUDA/WP/body_inline.h b/benchmarks/CUDA/WP/body_inline.h
new file mode 100644
index 0000000..5da0941
--- /dev/null
+++ b/benchmarks/CUDA/WP/body_inline.h
@@ -0,0 +1,741 @@
+#if defined(DEVICEEMU) && defined(DEBUGOUTPUT)
+for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
+kDIAGOUTPUT1(t)
+}
+for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
+kDIAGOUTPUT1(q)
+}
+for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
+kDIAGOUTPUT1(qc)
+}
+for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
+kDIAGOUTPUT1(qi)
+}
+for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
+kDIAGOUTPUT1(qr)
+}
+for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
+kDIAGOUTPUT1(qs)
+}
+for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
+kDIAGOUTPUT1(den)
+}
+for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
+kDIAGOUTPUT1(p)
+}
+for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
+kDIAGOUTPUT1(delz)
+}
+for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
+kDIAGOUTPUT1(cpm)
+}
+for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
+kDIAGOUTPUT1(xl)
+}
+#endif
+
+// 585 !----------------------------------------------------------------
+// 586 ! initialize the large scale variables
+ mstep = 1 ;
+
+ ttp=t0c+0.01 ;
+ dldt=cvap-cliq ;
+ xa=-dldt/rv ;
+ xb=xa+hvap/(rv*ttp) ;
+ dldti=cvap-cice ;
+ xai=-dldti/rv ;
+ xbi=xai+hsub/(rv*ttp) ;
+
+ for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
+
+ pp = p[k] ;
+ tt = t[k] ;
+ tr = ttp/tt ;
+ ltr = log(tr) ;
+
+ qq=psat*exp(ltr*(xa)+xb*(1.-tr)) ;
+ qq=ep2*qq/(pp-qq) ;
+ qs1[k] = MAX(qq,qmin) ;
+ rh1[k] = MAX( q[k]/qs1[k],qmin) ;
+
+ if( tt < ttp ) {
+ qq=psat*exp(ltr*(xai)+xbi*(1.-tr)) ;
+ } else {
+ qq=psat*exp(ltr*(xa)+xb*(1.-tr)) ;
+ }
+ qq = ep2 * qq / (pp - qq) ;
+ qs2[k] = MAX(qq,qmin) ;
+ rh2[k] = MAX(q[k]/qs2[k],qmin) ;
+
+ }
+
+ for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
+ prevp[k] = 0. ;
+ psdep[k] = 0. ;
+ praut[k] = 0. ;
+ psaut[k] = 0. ;
+ pracw[k] = 0. ;
+ psaci[k] = 0. ;
+ psacw[k] = 0. ;
+ pigen[k] = 0. ;
+ pidep[k] = 0. ;
+ pcond[k] = 0. ;
+ psmlt[k] = 0. ;
+ psevp[k] = 0. ;
+ falk1[k] = 0. ;
+ falk2[k] = 0. ;
+ fall1[k] = 0. ;
+ fall2[k] = 0. ;
+ fallc[k] = 0. ;
+ falkc[k] = 0. ;
+ xni[k] = 1.e3 ;
+ }
+
+#define LAMDAR(x,y) sqrt(sqrt(pidn0r/((x)*(y))))
+#define LAMDAS(x,y,z) sqrt(sqrt(pidn0s*(z)/((x)*(y))))
+// calculate mstep for this colum
+
+
+ for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
+ float supcol = t0c - t[k] ;
+#ifdef DEVICEEMU
+if ( ig == IDEBUG && jg == JDEBUG && k+1 == KDEBUG ) fprintf(stderr,"ZAP t0c %25.17e\n",t0c) ;
+if ( ig == IDEBUG && jg == JDEBUG && k+1 == KDEBUG ) fprintf(stderr,"ZAP supcol %25.17e\n",supcol) ;
+#endif
+DIAGOUTPUT1(t)
+ n0sfac[k] = MAX(MIN(exp(alpha*supcol),n0smax/n0s),1.) ;
+ if ( qr[k] <= qcrmin ) {
+ rsloper[k] = rslopermax ;
+ rslopebr[k] = rsloperbmax ;
+ rslope2r[k] = rsloper2max ;
+ rslope3r[k] = rsloper3max ;
+ } else {
+DIAGOUTPUT1(qr)
+DIAGOUTPUT1(den)
+ rsloper[k] = 1./LAMDAR(qr[k],den[k]) ;
+DIAGOUTPUT1(rsloper)
+ rslopebr[k] = exp(log(rsloper[k])*bvtr) ;
+ rslope2r[k] = rsloper[k] * rsloper[k] ;
+ rslope3r[k] = rslope2r[k] * rsloper[k] ;
+ }
+ if ( qs[k] <= qcrmin ) {
+ rslopes[k] = rslopesmax ;
+DIAGOUTPUT1(rslopes) ;
+ rslopebs[k] = rslopesbmax ;
+DIAGOUTPUT1(rslopebs) ;
+ rslope2s[k] = rslopes2max ;
+ rslope3s[k] = rslopes3max ;
+ } else {
+DIAGOUTPUT1(qs) ;
+DIAGOUTPUT1(den) ;
+DIAGOUTPUT1(n0sfac) ;
+ rslopes[k] = 1./LAMDAS(qs[k],den[k],n0sfac[k]) ;
+DIAGOUTPUT1(rslopes) ;
+ rslopebs[k] = exp(log(rslopes[k])*bvts) ;
+DIAGOUTPUT1(rslopebs) ;
+ rslope2s[k] = rslopes[k] * rslopes[k] ;
+ rslope3s[k] = rslope2s[k] * rslopes[k] ;
+ }
+ denfac[k] = sqrt(den0/den[k]) ;
+ w1[k] = pvtr*rslopebr[k]*denfac[k]/delz[k] ;
+ w2[k] = pvts*rslopebs[k]*denfac[k]/delz[k] ;
+
+DIAGOUTPUT1(w1)
+DIAGOUTPUT1(rslopebr)
+DIAGOUTPUT1(w2)
+DIAGOUTPUT1(rslopebs)
+DIAGOUTPUT1(denfac)
+DIAGOUTPUT1(delz)
+
+ w = MAX(w1[k],w2[k]) ;
+ numdt = MAX(trunc(w*dtcld+.5+.5),1) ;
+ if ( numdt >= mstep ) mstep = numdt ;
+//-------------------------------------------------------------
+// Ni: ice crystal number concentration [HDC 5c]
+//-------------------------------------------------------------
+ float temp = (den[k]*MAX(qi[k],qmin)) ;
+ temp = sqrt(sqrt(temp*temp*temp)) ;
+#ifdef DEBUGDEBUG
+ xni[k] = 1.e3 ;
+#else
+ xni[k] = MIN(MAX(5.38e7*temp,1.e3),1.e6) ;
+#endif
+ }
+ rmstep = 1./mstep ;
+
+ for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
+ fall1[k] = 0. ;
+ fall2[k] = 0. ;
+ }
+
+// diffus(x,y) = 8.794e-5 * exp(log(x)*(1.81)) / y ! 8.794e-5*x**1.81/y
+// viscos(x,y) = 1.496e-6 * (x*sqrt(x)) /(x+120.)/y ! 1.496e-6*x**1.5/(x+120.)/y
+// xka(x,y) = 1.414e3*viscos(x,y)*y
+// diffac(a,b,c,d,e) = d*a*a/(xka(c,d)*rv*c*c)+1./(e*diffus(c,b))
+// venfac(a,b,c) = exp(log((viscos(b,c)/diffus(b,a)))*((.3333333))) &
+// /sqrt(viscos(b,c))*sqrt(sqrt(den0/c))
+
+#define DIFFUS(x,y) (8.794e-5 * exp(log(x)*(1.81)) / (y))
+#define VISCOS(x,y) (1.496e-6 * ((x)*sqrt(x)) /((x)+120.)/(y))
+#define XKA(x,y) (1.414e3*VISCOS((x),(y))*(y))
+#define DIFFAC(a,b,c,d,e) ((d)*(a)*(a)/(XKA((c),(d))*rv*(c)*(c))+1./((e)*DIFFUS((c),(b))))
+#define VENFAC(a,b,c) (exp(log((VISCOS((b),(c))/DIFFUS((b),(a))))*((.3333333)))/sqrt(VISCOS((b),(c)))*sqrt(sqrt(den0/(c))))
+#define CONDEN(a,b,c,d,e) ((MAX((b),qmin)-(c))/(1.+(d)*(d)/(rv*(e))*(c)/((a)*(a))))
+
+ for ( n = 1 ; n <= mstep ; n++ ) {
+ k = kpe - 1 ;
+ falk1[k] = den[k]*qr[k]*w1[k]*rmstep ;
+ fall1[k] += falk1[k] ;
+ falk2[k] = den[k]*qs[k]*w2[k]*rmstep ;
+ fall2[k] += falk2[k] ;
+ qr[k] = MAX(qr[k]-falk1[k]*dtcldden,0.) ;
+ qs[k] = MAX(qs[k]-falk2[k]*dtcldden,0.) ;
+ for ( k = kpe-2 ; k >= kps-1 ; k-- ) {
+ falk1[k] = den[k]*qr[k]*w1[k]*rmstep ;
+ fall1[k] += falk1[k] ;
+ falk2[k] = den[k]*qs[k]*w2[k]*rmstep ;
+ fall2[k] += falk2[k] ;
+ dtcldden = dtcld/den[k] ;
+ rdelz = 1./delz[k] ;
+DIAGOUTPUT1i(loop) ;
+DIAGOUTPUT1i(mstep) ;
+DIAGOUTPUT1i(n) ;
+DIAGOUTPUT1(qr) ;
+DIAGOUTPUT1(falk1) ;
+DIAGOUTPUT11(falk1) ;
+DIAGOUTPUT1(delz) ;
+DIAGOUTPUT11(delz) ;
+ qr[k] = MAX(qr[k]-
+ (falk1[k]-falk1[k+1]*delz[k+1]*rdelz)*
+ dtcldden,0.) ;
+DIAGOUTPUT1(qr) ;
+DIAGOUTPUT1(qs) ;
+DIAGOUTPUT1(falk2) ;
+DIAGOUTPUT1(w2) ;
+DIAGOUTPUT11(falk2) ;
+ qs[k] = MAX(qs[k]-
+ (falk2[k]-falk2[k+1]*delz[k+1]*rdelz)*
+ dtcldden,0.) ;
+DIAGOUTPUT1(qs) ;
+ }
+
+ for ( k = kpe-1 ; k >= kps-1 ; k-- ) {
+DIAGOUTPUT1(t) ;
+DIAGOUTPUT1(qs) ;
+ if ( t[k] > t0c && qs[k] > 0.) {
+ xlf = xlf0 ;
+ w3[k] = VENFAC(p[k],t[k],den[k]) ;
+ coeres = rslope2s[k]*sqrt(rslopes[k]*rslopebs[2]) ;
+ psmlt[k] = XKA(t[k],den[k])/xlf*(t0c-t[k])*pi/2.
+ *n0sfac[k]*(precs1*rslope2s[k]+precs2
+ *w3[k]*coeres) ;
+ psmlt[k] = MIN(MAX(psmlt[k]*dtcld*rmstep,-qs[k]*rmstep),0.) ;
+ qs[k] += psmlt[k] ;
+DIAGOUTPUT1i(mstep) ;
+DIAGOUTPUT1i(n) ;
+DIAGOUTPUT1(qr) ;
+DIAGOUTPUT1(psmlt) ;
+ qr[k] -= psmlt[k] ;
+DIAGOUTPUT1(qr) ;
+
+ t[k] += xlf/CPMCAL(q[k])*psmlt[k] ;
+ }
+ }
+ }
+
+//---------------------------------------------------------------
+// Vice [ms-1] : fallout of ice crystal [HDC 5a]
+//---------------------------------------------------------------
+ mstep = 1 ;
+ numdt = 1 ;
+ for ( k = kpe-1 ; k >= kps-1 ; k-- ) {
+ if (qi[k] <= 0.) {
+ w2[k] = 0. ;
+ } else {
+ xmi = den[k]*qi[k]/xni[k] ;
+ diameter = MAX(MIN(dicon * sqrt(xmi),dimax), 1.e-25) ;
+ w1[k] = 1.49e4*exp(log(diameter)*(1.31)) ;
+ w2[k] = w1[k]/delz[k] ;
+ }
+ numdt = MAX( trunc(w2[k]*dtcld+.5+.5),1) ;
+ if(numdt > mstep) mstep = numdt ;
+ }
+ rmstep = 1./mstep ;
+
+ for ( n = 1 ; n <= mstep ; n++ ) {
+ k = kpe - 1 ;
+ falkc[k] = den[k]*qi[k]*w2[k]*rmstep ;
+ fallc[k] = fallc[k]+falkc[k] ;
+ qi[k] = MAX(qi[k]-falkc[k]*dtcld/den[k],0.) ;
+ for ( k = kpe-2 ; k >= kps-1 ; k-- ) {
+ falkc[k] = den[k]*qi[k]*w2[k]*rmstep ;
+ fallc[k] = fallc[k]+falkc[k] ;
+ qi[k] = MAX(qi[k]-(falkc[k]-falkc[k+1]
+ *delz[k+1]/delz[k])*dtcld/den[k],0.) ;
+ }
+ }
+ fallsum = fall1[1]+fall2[1]+fallc[1] ;
+ fallsum_qsi = fall2[1]+fallc[1] ;
+ rainncv = 0. ;
+ if(fallsum > 0.) {
+ rainncv = fallsum*delz[1]/denr*dtcld*1000. ;
+ rain = fallsum*delz[1]/denr*dtcld*1000. + rain ;
+ }
+ snowncv = 0. ;
+ if(fallsum_qsi > 0.) {
+ snowncv = fallsum_qsi*delz[0]/denr*dtcld*1000. ;
+ snow = fallsum_qsi*delz[0]/denr*dtcld*1000. + snow ;
+ }
+ sr = 0. ;
+ if ( fallsum > 0. ) sr = fallsum_qsi*delz[0]/denr*dtcld*1000./(rainncv+1.e-12) ;
+
+//---------------------------------------------------------------
+// pimlt: instantaneous melting of cloud ice [HL A47] [RH83 A28]
+// (T>T0: I->C)
+//---------------------------------------------------------------
+ for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
+ float supcol = t0c-t[k] ;
+ xlf = xls-xl[k] ;
+ if( supcol < 0. ) xlf = xlf0 ;
+ if( supcol < 0 && qi[k] > 0. ) {
+ qc[k] = qc[k] + qi[k] ;
+ t[k] = t[k] - xlf/cpm[k]*qi[k] ;
+ qi[k] = 0. ;
+ }
+//---------------------------------------------------------------
+// pihmf: homogeneous freezing of cloud water below -40c [HL A45]
+// (T<-40C: C->I)
+//---------------------------------------------------------------
+ if( supcol > 40. && qc[k] > 0. ) {
+ qi[k] = qi[k] + qc[k] ;
+ t[k] = t[k] + xlf/cpm[k]*qc[k] ;
+ qc[k] = 0. ;
+ }
+//---------------------------------------------------------------
+// pihtf: heterogeneous freezing of cloud water [HL A44]
+// (T0>T>-40C: C->I)
+//---------------------------------------------------------------
+ if ( supcol > 0. && qc[k] > 0.) {
+ float pfrzdtc = MIN(pfrz1*(exp(pfrz2*supcol)-1.)
+ *den[k]/denr/xncr*qc[k]*qc[k]*dtcld,qc[k]) ;
+ qi[k] = qi[k] + pfrzdtc ;
+ t[k] = t[k] + xlf/cpm[k]*pfrzdtc ;
+ qc[k] = qc[k]-pfrzdtc ;
+ }
+//---------------------------------------------------------------
+// psfrz: freezing of rain water [HL A20] [LFO 45]
+// (T<T0, R->S)
+//---------------------------------------------------------------
+ if( supcol > 0. && qr[k] > 0. ) {
+ float temp = rsloper[k] ;
+ temp = temp*temp*temp*temp*temp*temp*temp ;
+ float pfrzdtr = MIN(20.*(pi*pi)*pfrz1*n0r*denr/den[k]
+ *(exp(pfrz2*supcol)-1.)*temp*dtcld,
+ qr[k]) ;
+ qs[k] = qs[k] + pfrzdtr ;
+ t[k] = t[k] + xlf/cpm[k]*pfrzdtr ;
+ qr[k] = qr[k]-pfrzdtr ;
+DIAGOUTPUT1(qr) ;
+ }
+ }
+
+//----------------------------------------------------------------
+// rsloper: reverse of the slope parameter of the rain(m)
+// xka: thermal conductivity of air(jm-1s-1k-1)
+// work1: the thermodynamic term in the denominator associated with
+// heat conduction and vapor diffusion
+// (ry88, y93, h85)
+// work2: parameter associated with the ventilation effects(y93)
+
+ for ( k = kps-1 ; k <= kpe - 1 ; k++ ) {
+ float supcol = t0c - t[k] ;
+ n0sfac[k] = MAX(MIN(exp(alpha*supcol),n0smax/n0s),1.) ;
+ if ( qr[k] <= qcrmin ) {
+ rsloper[k] = rslopermax ;
+ rslopebr[k] = rsloperbmax ;
+ rslope2r[k] = rsloper2max ;
+ rslope3r[k] = rsloper3max ;
+ } else {
+ rsloper[k] = 1./(sqrt(sqrt(pidn0r/((qr[k])*(den[k]))))) ;
+DIAGOUTPUT1(rsloper) ;
+DIAGOUTPUT1(qr) ;
+DIAGOUTPUT1(den) ;
+ rslopebr[k] = exp(log(rsloper[k])*bvtr) ;
+ rslope2r[k] = rsloper[k] * rsloper[k] ;
+ rslope3r[k] = rslope2r[k] * rsloper[k] ;
+ }
+ if ( qs[k] <= qcrmin ) {
+ rslopes[k] = rslopesmax ;
+ rslopebs[k] = rslopesbmax ;
+ rslope2s[k] = rslopes2max ;
+ rslope3s[k] = rslopes3max ;
+ } else {
+ rslopes[k] = 1./(sqrt(sqrt(pidn0s*(n0sfac[k])/((qs[k])*(den[k]))))) ;
+ rslopebs[k] = exp(log(rslopes[k])*bvts) ;
+ rslope2s[k] = rslopes[k] * rslopes[k] ;
+ rslope3s[k] = rslope2s[k] * rslopes[k] ;
+ }
+ }
+
+ for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
+ w1[k] = DIFFAC(xl[k],p[k],t[k],den[k],qs1[k]) ;
+ w2[k] = DIFFAC(xls,p[k],t[k],den[k],qs2[k]) ;
+ w3[k] = VENFAC(p[k],t[k],den[k]) ;
+ }
+
+//
+//===============================================================
+//
+// warm rain processes
+//
+// - follows the processes in RH83 and LFO except for autoconcersion
+//
+//===============================================================
+//
+ for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
+ float supsat = MAX(q[k],qmin)-qs1[k] ;
+ float satdt = supsat/dtcld ;
+//---------------------------------------------------------------
+// praut: auto conversion rate from cloud to rain [HDC 16]
+// (C->R)
+//---------------------------------------------------------------
+ if(qc[k] > qc0) {
+ praut[k] = qck1*exp(log(qc[k])*((7./3.))) ;
+ praut[k] = MIN(praut[k],qc[k]/dtcld) ;
+ }
+//---------------------------------------------------------------
+// pracw: accretion of cloud water by rain [HL A40] [LFO 51]
+// (C->R)
+//---------------------------------------------------------------
+ if(qr[k] > qcrmin && qc[k] > qmin) {
+ pracw[k] = MIN(pacrr*rslope3r[k]*rslopebr[k]
+ *qc[k]*denfac[k],qc[k]/dtcld) ;
+ }
+//---------------------------------------------------------------
+// prevp: evaporation/condensation rate of rain [HDC 14]
+// (V->R or R->V)
+//---------------------------------------------------------------
+ if(qr[k] > 0.) {
+ coeres = rslope2r[k]*sqrt(rsloper[k]*rslopebr[k]) ;
+ prevp[k] = (rh1[k]-1.)*(precr1*rslope2r[k]
+ +precr2*w3[k]*coeres)/w1[k] ;
+DIAGOUTPUT1(prevp) ;
+DIAGOUTPUT1(qr) ;
+DIAGOUTPUT1(rsloper) ;
+DIAGOUTPUT1(rslope2r) ;
+DIAGOUTPUT1(rslopebr) ;
+DIAGOUTPUT1(w1) ;
+DIAGOUTPUT1(rh1) ;
+ if(prevp[k] < 0.) {
+ prevp[k] = MAX(prevp[k],-qr[k]/dtcld) ;
+ prevp[k] = MAX(prevp[k],satdt/2) ;
+ } else {
+ prevp[k] = MIN(prevp[k],satdt/2) ;
+ }
+ }
+ }
+
+//
+//===============================================================
+//
+// cold rain processes
+//
+// - follows the revised ice microphysics processes in HDC
+// - the processes same as in RH83 and RH84 and LFO behave
+// following ice crystal hapits defined in HDC, inclduing
+// intercept parameter for snow (n0s), ice crystal number
+// concentration (ni), ice nuclei number concentration
+// (n0i), ice diameter (d)
+//
+//===============================================================
+//
+ rdtcld = 1./dtcld ;
+ for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
+ float supcol = t0c-t[k] ;
+ float supsat = MAX(q[k],qmin)-qs2[k] ;
+ float satdt = supsat/dtcld ;
+ int ifsat = 0 ;
+//-------------------------------------------------------------
+// Ni: ice crystal number concentraiton [HDC 5c]
+//-------------------------------------------------------------
+ float temp = (den[k]*MAX(qi[k],qmin)) ;
+ temp = sqrt(sqrt(temp*temp*temp)) ;
+ xni[k] = MIN(MAX(5.38e7*temp,1.e3),1.e6) ;
+ float eacrs = exp(0.07*(-supcol)) ;
+//
+ if(supcol > 0) {
+ if(qs[k] > qcrmin && qi[k] > qmin) {
+ xmi = den[k]*qi[k]/xni[k] ;
+ diameter = MIN(dicon * sqrt(xmi),dimax) ;
+ vt2i = 1.49e4*pow(diameter,(float)1.31) ;
+ vt2s = pvts*rslopebs[k]*denfac[k] ;
+//-------------------------------------------------------------
+// psaci: Accretion of cloud ice by rain [HDC 10]
+// (T<T0: I->S)
+//-------------------------------------------------------------
+ acrfac = 2.*rslope3s[k]+2.*diameter*rslope2s[k]
+ +diameter*diameter*rslopes[k] ;
+ psaci[k] = pi*qi[k]*eacrs*n0s*n0sfac[k]
+ *abs(vt2s-vt2i)*acrfac*.25 ;
+ }
+//-------------------------------------------------------------
+// psacw: Accretion of cloud water by snow [HL A7] [LFO 24]
+// (T<T0: C->S, and T>=T0: C->R)
+//-------------------------------------------------------------
+ if(qs[k] > qcrmin && qc[k] > qmin) {
+ psacw[k] = MIN(pacrc*n0sfac[k]*rslope3s[k]
+ *rslopebs[k]*qc[k]*denfac[k]
+ ,qc[k]*rdtcld) ;
+ }
+//-------------------------------------------------------------
+// pidep: Deposition/Sublimation rate of ice [HDC 9]
+// (T<T0: V->I or I->V)
+//-------------------------------------------------------------
+ if(qi[k] > 0 && ifsat != 1) {
+ xmi = den[k]*qi[k]/xni[k] ;
+ diameter = dicon * sqrt(xmi) ;
+ pidep[k] = 4.*diameter*xni[k]*(rh2[k]-1.)/w2[k] ;
+ supice = satdt-prevp[k] ;
+ if(pidep[k] < 0.) {
+ pidep[k] = MAX(MAX(pidep[k],satdt*.5),supice) ;
+ pidep[k] = MAX(pidep[k],-qi[k]*rdtcld) ;
+ } else {
+ pidep[k] = MIN(MIN(pidep[k],satdt*.5),supice) ;
+ }
+ if(abs(prevp[k]+pidep[k]) >= abs(satdt)) ifsat = 1 ;
+ }
+ }
+//-------------------------------------------------------------
+// psdep: deposition/sublimation rate of snow [HDC 14]
+// (V->S or S->V)
+//-------------------------------------------------------------
+ if( qs[k] > 0. && ifsat != 1) {
+ coeres = rslope2s[k]*sqrt(rslopes[k]*rslopebs[k]) ;
+ psdep[k] = (rh2[k]-1.)*n0sfac[k]
+ *(precs1*rslope2s[k]+precs2
+ *w3[k]*coeres)/w2[k] ;
+ supice = satdt-prevp[k]-pidep[k] ;
+ if(psdep[k] < 0.) {
+ psdep[k] = MAX(psdep[k],-qs[k]*rdtcld) ;
+ psdep[k] = MAX(MAX(psdep[k],satdt*.5),supice) ;
+ } else {
+ psdep[k] = MIN(MIN(psdep[k],satdt*.5),supice) ;
+ }
+ if(abs(prevp[k]+pidep[k]+psdep[k]) >= abs(satdt))
+ ifsat = 1 ;
+ }
+//-------------------------------------------------------------
+// pigen: generation(nucleation) of ice from vapor [HL A50] [HDC 7-8]
+// (T<T0: V->I)
+//-------------------------------------------------------------
+ if(supcol > 0) {
+ if(supsat > 0 && ifsat != 1) {
+ supice = satdt-prevp[k]-pidep[k]-psdep[k] ;
+ xni0 = 1.e3*exp(0.1*supcol) ;
+ roqi0 = 4.92e-11*exp(log(xni0)*(1.33));
+ pigen[k] = MAX(0.,(roqi0/den[k]-MAX(qi[k],0.))
+ *rdtcld) ;
+ pigen[k] = MIN(MIN(pigen[k],satdt),supice) ;
+ }
+//
+//-------------------------------------------------------------
+// psaut: conversion(aggregation) of ice to snow [HDC 12]
+// (T<T0: I->S)
+//-------------------------------------------------------------
+ if(qi[k] > 0.) {
+ qimax = roqimax/den[k] ;
+ psaut[k] = MAX(0.,(qi[k]-qimax)*rdtcld) ;
+ }
+ }
+//-------------------------------------------------------------
+// psevp: Evaporation of melting snow [HL A35] [RH83 A27]
+// (T>T0: S->V)
+//-------------------------------------------------------------
+ if(supcol < 0.) {
+ if(qs[k] > 0. && rh1[k] < 1.) {
+ psevp[k] = psdep[k]*w2[k]/w1[k] ;
+ } // asked Jimy about this, 11.6.07, JM
+ psevp[k] = MIN(MAX(psevp[k],-qs[k]*rdtcld),0.) ;
+ }
+ }
+
+
+//
+//
+//----------------------------------------------------------------
+// check mass conservation of generation terms and feedback to the
+// large scale
+//
+ for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
+ if(t[k]<=t0c) {
+//
+// cloud water
+//
+ value = MAX(qmin,qc[k]) ;
+ source = (praut[k]+pracw[k]+psacw[k])*dtcld ;
+ if (source > value) {
+ factor = value/source ;
+ praut[k] = praut[k]*factor ;
+ pracw[k] = pracw[k]*factor ;
+ psacw[k] = psacw[k]*factor ;
+ }
+//
+// cloud ice
+//
+ value = MAX(qmin,qi[k]) ;
+ source = (psaut[k]+psaci[k]-pigen[k]-pidep[k])*dtcld ;
+ if (source > value) {
+ factor = value/source ;
+ psaut[k] = psaut[k]*factor ;
+ psaci[k] = psaci[k]*factor ;
+ pigen[k] = pigen[k]*factor ;
+ pidep[k] = pidep[k]*factor ;
+ }
+//
+ w3[k]=-(prevp[k]+psdep[k]+pigen[k]+pidep[k]) ;
+// update
+DIAGOUTPUT1(q) ;
+DIAGOUTPUT1(prevp) ;
+DIAGOUTPUT1(psdep) ;
+DIAGOUTPUT1(pigen) ;
+DIAGOUTPUT1(pidep) ;
+ q[k] = q[k]+w3[k]*dtcld ;
+DIAGOUTPUT1(q) ;
+ qc[k] = MAX(qc[k]-(praut[k]+pracw[k]+psacw[k])*dtcld,0.) ;
+ qr[k] = MAX(qr[k]+(praut[k]+pracw[k]+prevp[k])*dtcld,0.) ;
+ qi[k] = MAX(qi[k]-(psaut[k]+psaci[k]-pigen[k]-pidep[k])*dtcld,0.) ;
+DIAGOUTPUT1(qs)
+ qs[k] = MAX(qs[k]+(psdep[k]+psaut[k]+psaci[k]+psacw[k])*dtcld,0.) ;
+DIAGOUTPUT1(qs)
+ xlf = xls-xl[k] ;
+ xlwork2 = -xls*(psdep[k]+pidep[k]+pigen[k])-xl[k]*prevp[k]-xlf*psacw[k] ;
+ t[k] = t[k]-xlwork2/cpm[k]*dtcld ;
+ } else {
+//
+// cloud water
+//
+ value = MAX(qmin,qc[k]) ;
+ source=(praut[k]+pracw[k]+psacw[k])*dtcld ;
+ if (source > value) {
+ factor = value/source ;
+ praut[k] = praut[k]*factor ;
+ pracw[k] = pracw[k]*factor ;
+ psacw[k] = psacw[k]*factor ;
+ }
+//
+// snow
+//
+ value = MAX(qcrmin,qs[k]) ;
+ source=(-psevp[k])*dtcld ;
+ if (source > value) {
+ factor = value/source ;
+ psevp[k] = psevp[k]*factor ;
+ }
+ w3[k]=-(prevp[k]+psevp[k]) ;
+// update
+DIAGOUTPUT1(q) ;
+DIAGOUTPUT1(prevp) ;
+DIAGOUTPUT1(psdep) ;
+DIAGOUTPUT1(pigen) ;
+DIAGOUTPUT1(pidep) ;
+ q[k] = q[k]+w3[k]*dtcld ;
+DIAGOUTPUT1(q) ;
+ qc[k] = MAX(qc[k]-(praut[k]+pracw[k]+psacw[k])*dtcld,0.) ;
+ qr[k] = MAX(qr[k]+(praut[k]+pracw[k]+prevp[k] +psacw[k])*dtcld,0.) ;
+DIAGOUTPUT1(qs)
+DIAGOUTPUT1(psevp)
+
+#ifdef DEVICEEMU
+if (ig == IDEBUG && jg == JDEBUG && k+1 == KDEBUG ) fprintf(stderr,"%8s %25.17e\n","ZAP p*dt",psevp[k]*dtcld) ;
+if (ig == IDEBUG && jg == JDEBUG && k+1 == KDEBUG ) fprintf(stderr,"%8s %25.17e\n","ZAP q+p*dt",qs[k]+psevp[k]*dtcld) ;
+#endif
+ qs[k] = MAX(qs[k]+psevp[k]*dtcld,0.) ;
+DIAGOUTPUT1(qs)
+ xlf = xls-xl[k] ;
+ xlwork2 = -xl[k]*(prevp[k]+psevp[k]) ;
+ t[k] = t[k]-xlwork2/cpm[k]*dtcld ;
+ }
+ }
+DIAGOUTPUT2(qs)
+//
+// Inline expansion for fpvs
+ cvap = cpv ;
+ ttp=t0c+0.01 ;
+ dldt=cvap-cliq ;
+ xa=-dldt/rv ;
+ xb=xa+hvap/(rv*ttp) ;
+ dldti=cvap-cice ;
+ xai=-dldti/rv ;
+ xbi=xai+hsub/(rv*ttp) ;
+ for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
+ tr=ttp/t[k] ;
+ qs1[k]=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr)) ;
+ qs1[k] = ep2 * qs1[k] / (p[k] - qs1[k]) ;
+ qs1[k] = MAX(qs1[k],qmin) ;
+ }
+//
+//----------------------------------------------------------------
+// pcond: condensational/evaporational rate of cloud water [HL A46] [RH83 A6]
+// if there exists additional water vapor condensated/if
+// evaporation of cloud water is not enough to remove subsaturation
+//
+ for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
+ w1[k] = ((MAX(q[k],qmin)-(qs1[k])) /
+ (1.+(xl[k])*(xl[k])/(rv*(cpm[k]))*(qs1[k])/((t[k])*(t[k])))) ;
+ // w3[k] = qc[k]+w1[k] ; NOT USED
+ pcond[k] = MIN(MAX(w1[k]/dtcld,0.),MAX(q[k],0.)/dtcld) ;
+ if(qc[k] > 0. && w1[k] < 0.) {
+ pcond[k] = MAX(w1[k],-qc[k])/dtcld ;
+ }
+DIAGOUTPUT1(q) ;
+DIAGOUTPUT1(pcond) ;
+DIAGOUTPUT1(qs1) ;
+ q[k] = q[k]-pcond[k]*dtcld ;
+DIAGOUTPUT1(q) ;
+ qc[k] = MAX(qc[k]+pcond[k]*dtcld,0.) ;
+ t[k] = t[k]+pcond[k]*xl[k]/cpm[k]*dtcld ;
+ }
+//
+//
+//----------------------------------------------------------------
+// padding for small values
+//
+ for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
+ if(qc[k] <= qmin) qc[k] = 0.0 ;
+ if(qi[k] <= qmin) qi[k] = 0.0 ;
+ }
+
+//////////// end of loop ////////////////
+for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
+kDIAGOUTPUT1(t)
+}
+for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
+kDIAGOUTPUT1(q)
+}
+for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
+kDIAGOUTPUT1(qc)
+}
+for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
+kDIAGOUTPUT1(qi)
+}
+for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
+kDIAGOUTPUT1(qr)
+}
+for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
+kDIAGOUTPUT1(qs)
+}
+for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
+kDIAGOUTPUT1(den)
+}
+for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
+kDIAGOUTPUT1(p)
+}
+for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
+kDIAGOUTPUT1(delz)
+}
+for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
+kDIAGOUTPUT1(cpm)
+}
+for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
+kDIAGOUTPUT1(xl)
+}
diff --git a/benchmarks/CUDA/WP/compare_snaps.F b/benchmarks/CUDA/WP/compare_snaps.F
new file mode 100644
index 0000000..cedbd39
--- /dev/null
+++ b/benchmarks/CUDA/WP/compare_snaps.F
@@ -0,0 +1,83 @@
+! file format
+!write(46) ids,ide,jds,jde,kds,kde,ims,ime,jms,jme,kms,kme,its,ite,jts,jte,kts,kte
+!write(46) th(its:ite,kts:kte,jts:jte)
+!write(46) q(its:ite,kts:kte,jts:jte)
+!write(46) qc(its:ite,kts:kte,jts:jte)
+!write(46) qr(its:ite,kts:kte,jts:jte)
+!write(46) qi(its:ite,kts:kte,jts:jte)
+!write(46) qs(its:ite,kts:kte,jts:jte)
+!!INOUT 2D
+!write(46) rain(its:ite,jts:jte)
+!write(46) rainncv(its:ite,jts:jte)
+!write(46) sr(its:ite,jts:jte)
+!write(46) snow(its:ite,jts:jte)
+!write(46) snowncv(its:ite,jts:jte)
+
+ program compare_snaps
+ implicit none
+ integer, parameter :: unit1 = 31
+ integer, parameter :: unit2 = 32
+ integer ids,ide,jds,jde,kds,kde,ims,ime,jms,jme,kms,kme,its,ite,jts,jte,kts,kte
+
+ open(unit=unit1,file="snap_out_010",status="old",form="unformatted",err=22)
+ open(unit=unit2,file="snap_gpu_010",status="old",form="unformatted",err=23)
+
+ read(unit=unit1) ids,ide,jds,jde,kds,kde,ims,ime,jms,jme,kms,kme,its,ite,jts,jte,kts,kte
+ write(*,*)ids,ide,jds,jde,kds,kde
+ write(*,*)ims,ime,jms,jme,kms,kme
+ write(*,*)its,ite,jts,jte,kts,kte
+ read(unit=unit2) ids,ide,jds,jde,kds,kde,ims,ime,jms,jme,kms,kme,its,ite,jts,jte,kts,kte
+ write(*,*)ids,ide,jds,jde,kds,kde
+ write(*,*)ims,ime,jms,jme,kms,kme
+ write(*,*)its,ite,jts,jte,kts,kte
+
+ call read_and_output3( 'th', unit1, 88 )
+ call read_and_output3( 'th', unit2, 98 )
+ call read_and_output3( 'q', unit1, 88 )
+ call read_and_output3( 'q', unit2, 98 )
+ call read_and_output3( 'qc', unit1, 88 )
+ call read_and_output3( 'qc', unit2, 98 )
+ call read_and_output3( 'qr', unit1, 88 )
+ call read_and_output3( 'qr', unit2, 98 )
+ call read_and_output3( 'qi', unit1, 88 )
+ call read_and_output3( 'qi', unit2, 98 )
+ call read_and_output3( 'qs', unit1, 88 )
+ call read_and_output3( 'qs', unit2, 98 )
+
+ close(unit1)
+ close(unit2)
+ stop
+22 write(*,*)'error opening ',unit1
+ stop
+23 write(*,*)'error opening ',unit2
+ stop
+ end
+
+ subroutine read_and_output3 ( varname,iu, ou )
+ implicit none
+ integer iu,ou
+ character*(*) varname
+! 2 72 2 59 1 27
+ integer, parameter :: ix = 71
+ integer, parameter :: jx = 58
+ integer, parameter :: kx = 27
+ integer i,j,k
+ real s3(ix,kx,jx)
+
+ read(iu)s3
+ write(ou,*)ix,jx,' ',trim(varname)
+ do j = 1,jx
+ do i = 1,ix
+ write(ou,*)s3(i,1,j)
+ enddo
+ enddo
+ return
+ end
+
+
+
+
+
+
+
+
diff --git a/benchmarks/CUDA/WP/cublas.h b/benchmarks/CUDA/WP/cublas.h
new file mode 100644
index 0000000..4453b30
--- /dev/null
+++ b/benchmarks/CUDA/WP/cublas.h
@@ -0,0 +1,2 @@
+//intentially empty file for running cpp in this directory
+//otherwise, nvcc gets the cublas.h file from the cuda dir
diff --git a/benchmarks/CUDA/WP/d2s.h b/benchmarks/CUDA/WP/d2s.h
new file mode 100644
index 0000000..973f550
--- /dev/null
+++ b/benchmarks/CUDA/WP/d2s.h
@@ -0,0 +1,36 @@
+
+delt_s = delt
+g_s = g
+rd_s = rd
+rv_s = rv
+t0c_s = t0c
+den0_s = den0
+cpd_s = cpd
+cpv_s = cpv
+ep1_s = ep1
+ep2_s = ep2
+qmin_s = qmin
+XLS_s = XLS
+XLV0_s = XLV0
+XLF0_s = XLF0
+cliq_s = cliq
+cice_s = cice
+psat_s = psat
+denr_s = denr
+
+th_s(its:ite,kts:kte,jts:jte) = th(its:ite,kts:kte,jts:jte)
+pii_s(its:ite,kts:kte,jts:jte) = pii(its:ite,kts:kte,jts:jte)
+q_s(its:ite,kts:kte,jts:jte) = q(its:ite,kts:kte,jts:jte)
+qc_s(its:ite,kts:kte,jts:jte) = qc(its:ite,kts:kte,jts:jte)
+qi_s(its:ite,kts:kte,jts:jte) = qi(its:ite,kts:kte,jts:jte)
+qr_s(its:ite,kts:kte,jts:jte) = qr(its:ite,kts:kte,jts:jte)
+qs_s(its:ite,kts:kte,jts:jte) = qs(its:ite,kts:kte,jts:jte)
+den_s(its:ite,kts:kte,jts:jte) = den(its:ite,kts:kte,jts:jte)
+p_s(its:ite,kts:kte,jts:jte) = p(its:ite,kts:kte,jts:jte)
+delz_s(its:ite,kts:kte,jts:jte) = delz(its:ite,kts:kte,jts:jte)
+rain_s(its:ite,jts:jte) = rain(its:ite,jts:jte)
+rainncv_s(its:ite,jts:jte) = rainncv(its:ite,jts:jte)
+sr_s(its:ite,jts:jte) = sr(its:ite,jts:jte)
+snow_s(its:ite,jts:jte) = snow(its:ite,jts:jte)
+snowncv_s(its:ite,jts:jte) = snowncv(its:ite,jts:jte)
+
diff --git a/benchmarks/CUDA/WP/data/wsm5_in_010 b/benchmarks/CUDA/WP/data/wsm5_in_010
new file mode 100644
index 0000000..d1d502b
--- /dev/null
+++ b/benchmarks/CUDA/WP/data/wsm5_in_010
Binary files differ
diff --git a/benchmarks/CUDA/WP/debug.m4 b/benchmarks/CUDA/WP/debug.m4
new file mode 100644
index 0000000..0d6ba17
--- /dev/null
+++ b/benchmarks/CUDA/WP/debug.m4
@@ -0,0 +1,35 @@
+define( DIAGOUTPUT1, `
+#if defined(DEVICEEMU) && defined(DEBUGOUTPUT)
+ if (ig==IDEBUG&&jg==JDEBUG&&k+1==KDEBUG) fprintf(stderr,"ZAP %8s %25.17e \n", "$1", $1[k] );
+#endif
+')
+define( DIAGOUTPUT1i, `
+#if defined(DEVICEEMU) && defined(DEBUGOUTPUT)
+ if (ig==IDEBUG&&jg==JDEBUG&&k+1==KDEBUG) fprintf(stderr,"ZAP %8s %20d \n", "$1", $1 );
+#endif
+')
+define( DIAGOUTPUT11, `
+#if defined(DEVICEEMU) && defined(DEBUGOUTPUT)
+ if (ig==IDEBUG&&jg==JDEBUG&&k+1==KDEBUG) fprintf(stderr,"ZAP %8s %25.17e \n", "$1", $1[k+1] );
+#endif
+')
+define( DIAGOUTPUT2, `
+#if defined(DEVICEEMU) && defined(DEBUGOUTPUT)
+ if (ig==IDEBUG&&jg==JDEBUG) fprintf(stderr,"ZAP %8s %25.17e \n", "$1", $1[KDEBUG-1] );
+#endif
+')
+define( kDIAGOUTPUT1, `
+#if defined(DEVICEEMU) && defined(DEBUGOUTPUT)
+ if (ig==IDEBUG&&jg==JDEBUG) fprintf(stderr,"ZAP %8s %25.17e \n", "$1", $1[k] );
+#endif
+')
+define( kDIAGOUTPUT1i, `
+#if defined(DEVICEEMU) && defined(DEBUGOUTPUT)
+ if (ig==IDEBUG&&jg==JDEBUG) fprintf(stderr,"ZAP %8s %20d \n", "$1", $1 );
+#endif
+')
+define( kDIAGOUTPUT11, `
+#if defined(DEVICEEMU) && defined(DEBUGOUTPUT)
+ if (ig==IDEBUG&&jg==JDEBUG) fprintf(stderr,"ZAP %8s %25.17e \n", "$1", $1[k+1] );
+#endif
+')
diff --git a/benchmarks/CUDA/WP/libmassv.F b/benchmarks/CUDA/WP/libmassv.F
new file mode 100644
index 0000000..344946b
--- /dev/null
+++ b/benchmarks/CUDA/WP/libmassv.F
@@ -0,0 +1,385 @@
+! IBM libmassv compatibility library
+!
+
+#ifndef NATIVE_MASSV
+ subroutine vdiv(z,x,y,n)
+ double precision x(*),y(*),z(*)
+ do 10 j=1,n
+ z(j)=x(j)/y(j)
+ 10 continue
+ return
+ end
+
+ subroutine vsdiv(z,x,y,n)
+ real x(*),y(*),z(*)
+ do 10 j=1,n
+ z(j)=x(j)/y(j)
+ 10 continue
+ return
+ end
+
+ subroutine vexp(y,x,n)
+ double precision x(*),y(*)
+ do 10 j=1,n
+ y(j)=exp(x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vsexp(y,x,n)
+ real x(*),y(*)
+ do 10 j=1,n
+ y(j)=exp(x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vlog(y,x,n)
+ double precision x(*),y(*)
+ do 10 j=1,n
+ y(j)=log(x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vslog(y,x,n)
+ real x(*),y(*)
+ do 10 j=1,n
+ y(j)=log(x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vrec(y,x,n)
+ double precision x(*),y(*)
+ do 10 j=1,n
+ y(j)=1.d0/x(j)
+ 10 continue
+ return
+ end
+
+ subroutine vsrec(y,x,n)
+ real x(*),y(*)
+ do 10 j=1,n
+ y(j)=1.d0/x(j)
+ 10 continue
+ return
+ end
+
+ subroutine vrsqrt(y,x,n)
+ double precision x(*),y(*)
+ do 10 j=1,n
+ y(j)=1.d0/sqrt(x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vsrsqrt(y,x,n)
+ real x(*),y(*)
+ do 10 j=1,n
+ y(j)=1.d0/sqrt(x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vsincos(x,y,z,n)
+ double precision x(*),y(*),z(*)
+ do 10 j=1,n
+ x(j)=sin(z(j))
+ y(j)=cos(z(j))
+ 10 continue
+ return
+ end
+
+ subroutine vssincos(x,y,z,n)
+ real x(*),y(*),z(*)
+ do 10 j=1,n
+ x(j)=sin(z(j))
+ y(j)=cos(z(j))
+ 10 continue
+ return
+ end
+
+ subroutine vsqrt(y,x,n)
+ double precision x(*),y(*)
+ do 10 j=1,n
+ y(j)=sqrt(x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vssqrt(y,x,n)
+ real x(*),y(*)
+ do 10 j=1,n
+ y(j)=sqrt(x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vtan(y,x,n)
+ double precision x(*),y(*)
+ do 10 j=1,n
+ y(j)=tan(x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vstan(y,x,n)
+ real x(*),y(*)
+ do 10 j=1,n
+ y(j)=tan(x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vatan2(z,y,x,n)
+ double precision x(*),y(*),z(*)
+ do 10 j=1,n
+ z(j)=atan2(y(j),x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vsatan2(z,y,x,n)
+ real x(*),y(*),z(*)
+ do 10 j=1,n
+ z(j)=atan2(y(j),x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vasin(y,x,n)
+ double precision x(*),y(*)
+ do 10 j=1,n
+ y(j)=asin(x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vsin(y,x,n)
+ double precision x(*),y(*)
+ do 10 j=1,n
+ y(j)=sin(x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vssin(y,x,n)
+ real x(*),y(*)
+ do 10 j=1,n
+ y(j)=sin(x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vacos(y,x,n)
+ double precision x(*),y(*)
+ do 10 j=1,n
+ y(j)=acos(x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vcos(y,x,n)
+ double precision x(*),y(*)
+ do 10 j=1,n
+ y(j)=cos(x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vscos(y,x,n)
+ real x(*),y(*)
+ do 10 j=1,n
+ y(j)=cos(x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vcosisin(y,x,n)
+ complex*16 y(*)
+ double precision x(*)
+ do 10 j=1,n
+ y(j)=dcmplx(cos(x(j)),sin(x(j)))
+ 10 continue
+ return
+ end
+
+ subroutine vscosisin(y,x,n)
+ complex*8 y(*)
+ real x(*)
+ do 10 j=1,n
+ y(j)= cmplx(cos(x(j)),sin(x(j)))
+ 10 continue
+ return
+ end
+
+ subroutine vdint(y,x,n)
+ double precision x(*),y(*)
+ do 10 j=1,n
+! y(j)=dint(x(j))
+ y(j)=int(x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vdnint(y,x,n)
+ double precision x(*),y(*)
+ do 10 j=1,n
+! y(j)=dnint(x(j))
+ y(j)=nint(x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vlog10(y,x,n)
+ double precision x(*),y(*)
+ do 10 j=1,n
+ y(j)=log10(x(j))
+ 10 continue
+ return
+ end
+
+! subroutine vlog1p(y,x,n)
+! double precision x(*),y(*)
+! interface
+! double precision function log1p(%val(x))
+! double precision x
+! end function log1p
+! end interface
+! do 10 j=1,n
+! y(j)=log1p(x(j))
+! 10 continue
+! return
+! end
+
+ subroutine vcosh(y,x,n)
+ double precision x(*),y(*)
+ do 10 j=1,n
+ y(j)=cosh(x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vsinh(y,x,n)
+ double precision x(*),y(*)
+ do 10 j=1,n
+ y(j)=sinh(x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vtanh(y,x,n)
+ double precision x(*),y(*)
+ do 10 j=1,n
+ y(j)=tanh(x(j))
+ 10 continue
+ return
+ end
+
+! subroutine vexpm1(y,x,n)
+! double precision x(*),y(*)
+! interface
+! double precision function expm1(%val(x))
+! double precision x
+! end function expm1
+! end interface
+! do 10 j=1,n
+! y(j)=expm1(x(j))
+! 10 continue
+! return
+! end
+
+
+ subroutine vsasin(y,x,n)
+ real x(*),y(*)
+ do 10 j=1,n
+ y(j)=asin(x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vsacos(y,x,n)
+ real x(*),y(*)
+ do 10 j=1,n
+ y(j)=acos(x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vscosh(y,x,n)
+ real x(*),y(*)
+ do 10 j=1,n
+ y(j)=cosh(x(j))
+ 10 continue
+ return
+ end
+
+! subroutine vsexpm1(y,x,n)
+! real x(*),y(*)
+! interface
+! double precision function expm1(%val(x))
+! double precision x
+! end function expm1
+! end interface
+! do 10 j=1,n
+! y(j)=expm1(real(x(j),8))
+! 10 continue
+! return
+! end
+
+ subroutine vslog10(y,x,n)
+ real x(*),y(*)
+ do 10 j=1,n
+ y(j)=log10(x(j))
+ 10 continue
+ return
+ end
+
+! subroutine vslog1p(y,x,n)
+! real x(*),y(*)
+! interface
+! double precision function log1p(%val(x))
+! double precision x
+! end function log1p
+! end interface
+! do 10 j=1,n
+! y(j)=log1p(real(x(j),8))
+! 10 continue
+! return
+! end
+
+
+ subroutine vssinh(y,x,n)
+ real x(*),y(*)
+ do 10 j=1,n
+ y(j)=sinh(x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vstanh(y,x,n)
+ real x(*),y(*)
+ do 10 j=1,n
+ y(j)=tanh(x(j))
+ 10 continue
+ return
+ end
+#endif
+
+ subroutine vspow(z,y,x,n)
+ real x(*),y(*),z(*)
+ do 10 j=1,n
+ z(j)=y(j)**x(j)
+ 10 continue
+ return
+ end
+
+ subroutine vpow(z,y,x,n)
+ double precision x(*),y(*),z(*)
+ do 10 j=1,n
+ z(j)=y(j)**x(j)
+ 10 continue
+ return
+ end
+
diff --git a/benchmarks/CUDA/WP/makefile b/benchmarks/CUDA/WP/makefile
new file mode 100644
index 0000000..835d0a1
--- /dev/null
+++ b/benchmarks/CUDA/WP/makefile
@@ -0,0 +1,138 @@
+#SHELL = csh
+##################### intel compiler #############
+#CC = icc
+#FC = ifort
+#CFLAGS = -w -O3 -ip -c
+#FCFLAGS = -align all -FR -convert big_endian -fno-alias -fno-fnalias
+#OPT = -O3
+#LD = $(FC)
+#LDOPT =
+# settings for GTX 280
+#XXX = 32
+#YYY = 8
+# settings for 5600 Quadro and GTX 8800
+XXX = 8
+YYY = 8
+##################### gcc/gfortran ###############
+CC = gcc
+FC = gfortran
+FCFLAGS = -fconvert=big-endian -frecord-marker=4 -ffree-form
+OPT = -O3
+LD = $(FC)
+LDOPT =
+##################################################
+# do not change this definition, change the one further down
+FLOAT = float
+RWORDSIZE = 4
+##################### cuda location ##############
+# eces-shell
+#CUDALIBPATH = ~/emu/cuda/lib
+# ncsa
+#CUDALIBPATH = /usr/local/cuda/lib
+#LIBCUBLAS = $(CUDALIBPATH)/libcublas.so
+#LIBCUBLASEMU = $(CUDALIBPATH)/libcublasemu.so
+
+
+######################## THIS SECTION YOU CAN CHANGE ##################
+#
+# Hard coded number of levels (35 for conus, 28 for jan00)
+MKX = 28
+
+# uncomment this to use FLOAT4 data type (optimization)
+#FLOAT4 = -DFLOAT_4=4
+
+# this must always be defined but they do not do anything unless
+# other settings are set
+DEBUG_I = 59
+DEBUG_J = 45
+DEBUG_K = 1
+DEBUGOUTPUT = -DDEBUG_I=$(DEBUG_I) -DDEBUG_J=$(DEBUG_J) -DDEBUG_K=$(DEBUG_K) $(FLOAT4)
+
+# uncomment to run on emulator instead of the device
+#DEVICEEMU = -DDEVICEEMU
+#DEVICEEMU_NVCC = -deviceemu $(DEVICEEMU)
+#LIBCUBLAS = $(LIBCUBLASEMU)
+
+# uncomment to output detailed debug data output
+# must have DEVICEEMU settings uncommented above
+#DEBUGOUTPUT = -DDEBUGOUTPUT -DDEBUG_I=$(DEBUG_I) -DDEBUG_J=$(DEBUG_J) -DDEBUG_K=$(DEBUG_K) $(FLOAT4)
+
+# uncomment to allow settings to force closer agreement
+#DEBUGDEBUG = -DDEBUGDEBUG
+
+# uncomment for to promote to 8 byte floats
+# note, if you do this without DEVICE EMU above, compiler will complain it does not have enough shared mem
+#PROMOTE = -DPROMOTE
+#FLOAT = double
+#FCFLAGS = -fconvert=big-endian -frecord-marker=4 -fdefault-real-8
+#RWORDSIZE = 8
+######################## END OF SECTION YOU CAN CHANGE ################
+
+#NVCC = nvcc -DCUDA
+NVCC = nvcc -DCUDA #--ptxas-options=-v
+PHASE = -cuda -arch sm_11 --compiler-options -fno-strict-aliasing
+NVOPT = $(DEVICEEMU_NVCC) $(PROMOTE) $(DEBUGDEBUG) $(DEBUGOUTPUT) \
+ -DXXX=$(XXX) -DYYY=$(YYY) -DMKX=$(MKX) --host-compilation 'C++' --use_fast_math
+
+ROOTBINDIR ?= bin
+BINDIR ?= $(ROOTBINDIR)
+ROOTOBJDIR ?= obj
+LIBDIR := $(NVIDIA_CUDA_SDK_LOCATION)/lib
+COMMONDIR := $(NVIDIA_CUDA_SDK_LOCATION)/common
+GPGPUSIM_ROOT ?= ../../..
+INTERMED_FILES := *.cpp*.i *.cpp*.ii *.cu.c *.cudafe*.* *.fatbin.c *.cubin *.hash
+
+GPGPUSIM_CL:=$(shell p4 changes -m 1 //...\#have | awk '{print $$2}')
+NVOPENCC_VER:=$(shell $$CUDAHOME/open64/bin/nvopencc --version 2>&1 | awk '/GPGPU-Sim/ {printf("_nvopencc_CL%d", $$3);}')
+GPGPULINK = -L$(GPGPUSIM_ROOT)/lib/ -lcudart -L$(NVIDIA_CUDA_SDK_LOCATION)/C/lib/ -lcutil -lm -lz -ldl -lGL $(NEWLIBDIR) $(LIB)
+# /usr/lib64/gcc/x86_64-pc-linux-gnu/4.1.2/libstdc++.so.6
+.SUFFIXES :
+
+all : vanilla chocolate compare_snaps
+
+wsm5.f.cu : wsm5.cu
+ m4 wsm5.cu | sed "s/float/$(FLOAT)/g" > wsm5.f.cu
+
+wsm5.cu.o : wsm5.f.cu
+ $(NVCC) $(PHASE) $(NVOPT) wsm5.f.cu
+ $(CC) $(CFLAGS) -o wsm5.cu.o -c wsm5.f.cu.cpp
+
+
+wsm5_gpu.f.cu : wsm5_gpu.cu spt.h
+ m4 wsm5_gpu.cu | ./spt.pl | sed "s/float/$(FLOAT)/g" > wsm5_gpu.f.cu
+
+wsm5_gpu.cu.o : wsm5_gpu.f.cu
+ $(NVCC) $(PHASE) $(NVOPT) wsm5_gpu.f.cu
+ $(CC) $(CFLAGS) -o wsm5_gpu.cu.o -c wsm5_gpu.f.cu.cpp
+
+
+libmassv.o : libmassv.F
+ $(FC) -E -C -P libmassv.F > libmassv.f90
+ $(FC) -c $(OPT) $(FCFLAGS) libmassv.f90
+
+microclock.o : microclock.c
+ $(CC) -c $(CFLAGS) -DMKX=$(MKX) microclock.c
+
+vanilla : module_mp_wsm5.F libmassv.o microclock.o
+ $(FC) -E -C -P $(DEBUGDEBUG) $(DEVICEEMU) $(DEBUGOUTPUT) -DRWORDSIZE=$(RWORDSIZE) module_mp_wsm5.F > module_mp_wsm5.f90
+ $(FC) -c $(OPT) $(FCFLAGS) module_mp_wsm5.f90
+ $(LD) -o wsm5_driver_vanilla $(LDOPT) module_mp_wsm5.o libmassv.o microclock.o $(LIBCUBLASEMU) $(GPGPULINK)
+
+chocolate : module_mp_wsm5.F wsm5.cu.o wsm5_gpu.cu.o libmassv.o microclock.o
+ $(FC) -E -C -P $(DEBUGDEBUG) $(DEVICEEMU) $(DEBUGOUTPUT) -DRUN_ON_GPU -DRWORDSIZE=$(RWORDSIZE) module_mp_wsm5.F > module_mp_wsm5.f90
+ $(FC) -c $(OPT) $(FCFLAGS) module_mp_wsm5.f90
+ $(LD) -o gpgpu_ptx_sim__wsm5_standalone $(LDOPT) module_mp_wsm5.o wsm5.cu.o wsm5_gpu.cu.o libmassv.o microclock.o $(LIBCUBLAS) $(GPGPULINK)
+
+compare_snaps : compare_snaps.F
+ cp compare_snaps.F compare_snaps.f90
+ $(FC) -o compare_snaps $(FCFLAGS) compare_snaps.f90
+ rm -f compare_snaps.f90
+
+clean :
+ rm -f *.o *.cu.cpp *.f.cu wsm5_driver_* *.mod *.f90 gpgpu_ptx_sim__wsm5_standalone* compare_snaps gpgpu_inst_stats.txt snap_gpu_010 gpgpu_inst_stats.txt *.log *.ptx
+ make clean_imm
+clean_imm :
+ rm -f $(INTERMED_FILES)
+
+tar:
+ tar cvf wsm5gpu_`date +"%Y%m%d"`.tar *.cu *.F *.h *.m4 *.pl makefile
diff --git a/benchmarks/CUDA/WP/microclock.c b/benchmarks/CUDA/WP/microclock.c
new file mode 100644
index 0000000..3cfc9d7
--- /dev/null
+++ b/benchmarks/CUDA/WP/microclock.c
@@ -0,0 +1,33 @@
+#ifndef CRAY
+# ifdef NOUNDERSCORE
+# define RSL_INTERNAL_MICROCLOCK rsl_internal_microclock
+# else
+# ifdef F2CSTYLE
+# define RSL_INTERNAL_MICROCLOCK rsl_internal_microclock__
+# else
+# define RSL_INTERNAL_MICROCLOCK rsl_internal_microclock_
+# endif
+# endif
+#endif
+#include <sys/time.h>
+
+RSL_INTERNAL_MICROCLOCK ()
+{
+ struct timeval tb ;
+ struct timezone tzp ;
+ int isec ; /* seconds */
+ int usec ; /* microseconds */
+ int msecs ;
+ gettimeofday( &tb, &tzp ) ;
+ isec = tb.tv_sec ;
+ usec = tb.tv_usec ;
+ msecs = 1000000 * isec + usec ;
+ return(msecs) ;
+}
+
+c_pow_ ( float * a, float * b )
+{
+ *a = pow( *a , *b ) ;
+}
+
+
diff --git a/benchmarks/CUDA/WP/module_mp_wsm5.F b/benchmarks/CUDA/WP/module_mp_wsm5.F
new file mode 100644
index 0000000..f7132ba
--- /dev/null
+++ b/benchmarks/CUDA/WP/module_mp_wsm5.F
@@ -0,0 +1,1731 @@
+#define IDEBUG (DEBUG_I)
+#define JDEBUG (DEBUG_J)
+#define KDEBUG (DEBUG_K)
+
+#if defined(DEVICEEMU) && defined(DEBUGOUTPUT)
+# define DIAGOUTPUT1(b,c) if ( lat.EQ.JDEBUG.and.k.eq.KDEBUG.and.i.EQ.IDEBUG ) write(0,'("ZAP ",a8,e25.17)') b,c
+# define DIAGOUTPUT1i(b,c) if ( lat.EQ.JDEBUG.and.k.eq.KDEBUG.and.i.EQ.IDEBUG ) write(0,'("ZAP ",a8,I20)') b,c
+# define DIAGOUTPUT2(b,c) if ( lat.EQ.JDEBUG ) write(0,'("ZAP ",a8,e25.17)') b,c
+# define kDIAGOUTPUT1(b,c) if ( lat.EQ.JDEBUG.and.i.EQ.IDEBUG ) write(0,'("ZAP ",a8,e25.17)') b,c
+# define kDIAGOUTPUT1i(b,c) if ( lat.EQ.JDEBUG.and.i.EQ.IDEBUG ) write(0,'("ZAP ",a8,I20)') b,c
+#else
+# define DIAGOUTPUT1(b,c)
+# define DIAGOUTPUT1i(b,c)
+# define DIAGOUTPUT2(b,c)
+# define kDIAGOUTPUT1(b,c)
+# define kDIAGOUTPUT1i(b,c)
+#endif
+
+#define STANDALONE
+!#define OUTPUT_SNAPSHOTS
+!#define INPUT_SNAPSHOTS
+!#define DIAGS
+!#define RUN_ON_GPU
+
+#ifdef STANDALONE
+# define TENTIN
+# define TENTINOUT
+# define OPTY
+#else
+# define TENTIN ,INTENT(IN)
+# define TENTINOUT ,INTENT(INOUT)
+# define OPTY ,OPTIONAL
+#endif
+
+#if ( RWORDSIZE == 4 )
+# define VREC vsrec
+# define VSQRT vssqrt
+#else
+# define VREC vrec
+# define VSQRT vsqrt
+#endif
+
+!Including inline expansion statistical function
+MODULE module_mp_wsm5
+!
+!
+ REAL, PARAMETER, PRIVATE :: dtcldcr = 120.
+ REAL, PARAMETER, PRIVATE :: n0r = 8.e6
+ REAL, PARAMETER, PRIVATE :: avtr = 841.9
+ REAL, PARAMETER, PRIVATE :: bvtr = 0.8
+ REAL, PARAMETER, PRIVATE :: r0 = .8e-5 ! 8 microm in contrast to 10 micro m
+ REAL, PARAMETER, PRIVATE :: peaut = .55 ! collection efficiency
+ REAL, PARAMETER, PRIVATE :: xncr = 3.e8 ! maritime cloud in contrast to 3.e8 in tc80
+ REAL, PARAMETER, PRIVATE :: xmyu = 1.718e-5 ! the dynamic viscosity kgm-1s-1
+ REAL, PARAMETER, PRIVATE :: avts = 11.72
+ REAL, PARAMETER, PRIVATE :: bvts = .41
+ REAL, PARAMETER, PRIVATE :: n0smax = 1.e11 ! t=-90C unlimited
+ REAL, PARAMETER, PRIVATE :: lamdarmax = 8.e4
+ REAL, PARAMETER, PRIVATE :: lamdasmax = 1.e5
+ REAL, PARAMETER, PRIVATE :: lamdagmax = 6.e4
+ REAL, PARAMETER, PRIVATE :: betai = .6
+ REAL, PARAMETER, PRIVATE :: xn0 = 1.e-2
+ REAL, PARAMETER, PRIVATE :: dicon = 11.9
+ REAL, PARAMETER, PRIVATE :: di0 = 12.9e-6
+ REAL, PARAMETER, PRIVATE :: dimax = 500.e-6
+ REAL, PARAMETER, PRIVATE :: n0s = 2.e6 ! temperature dependent n0s
+ REAL, PARAMETER, PRIVATE :: alpha = .12 ! .122 exponen factor for n0s
+ REAL, PARAMETER, PRIVATE :: pfrz1 = 100.
+ REAL, PARAMETER, PRIVATE :: pfrz2 = 0.66
+ REAL, PARAMETER, PRIVATE :: qcrmin = 1.e-9
+ REAL, PARAMETER, PRIVATE :: t40c = 233.16
+ REAL, PARAMETER, PRIVATE :: eacrc = 1.0
+ REAL, SAVE :: &
+ qc0, qck1,bvtr1,bvtr2,bvtr3,bvtr4,g1pbr,&
+ g3pbr,g4pbr,g5pbro2,pvtr,eacrr,pacrr, &
+ precr1,precr2,xm0,xmmax,roqimax,bvts1, &
+ bvts2,bvts3,bvts4,g1pbs,g3pbs,g4pbs, &
+ g5pbso2,pvts,pacrs,precs1,precs2,pidn0r,&
+ pidn0s,xlv1,pacrc, &
+ rslopermax,rslopesmax,rslopegmax, &
+ rsloperbmax,rslopesbmax,rslopegbmax, &
+ rsloper2max,rslopes2max,rslopeg2max, &
+ rsloper3max,rslopes3max,rslopeg3max
+
+CONTAINS
+!===================================================================
+!
+#ifndef STANDALONE
+ SUBROUTINE wsm5(th, q, qc, qr, qi, qs &
+ ,den, pii, p, delz &
+ ,delt,g, cpd, cpv, rd, rv, t0c &
+ ,ep1, ep2, qmin &
+ ,XLS, XLV0, XLF0, den0, denr &
+ ,cliq,cice,psat &
+ ,rain, rainncv &
+ ,snow, snowncv &
+ ,sr &
+ ,ids,ide, jds,jde, kds,kde &
+ ,ims,ime, jms,jme, kms,kme &
+ ,its,ite, jts,jte, kts,kte &
+ )
+#else
+ SUBROUTINE wsm5 ( &
+ ids,ide, jds,jde, kds,kde &
+ ,ims,ime, jms,jme, kms,kme &
+ ,its,ite, jts,jte, kts,kte &
+ )
+#endif
+
+!-------------------------------------------------------------------
+ IMPLICIT NONE
+!-------------------------------------------------------------------
+!
+! This code is a 5-class mixed ice microphyiscs scheme (WSM5) of the WRF
+! Single-Moment MicroPhyiscs (WSMMP). The WSMMP assumes that ice nuclei
+! number concentration is a function of temperature, and seperate assumption
+! is developed, in which ice crystal number concentration is a function
+! of ice amount. A theoretical background of the ice-microphysics and related
+! processes in the WSMMPs are described in Hong et al. (2004).
+! Production terms in the WSM6 scheme are described in Hong and Lim (2006).
+! All units are in m.k.s. and source/sink terms in kgkg-1s-1.
+!
+! WSM5 cloud scheme
+!
+! Coded by Song-You Hong (Yonsei Univ.)
+! Jimy Dudhia (NCAR) and Shu-Hua Chen (UC Davis)
+! Summer 2002
+!
+! Implemented by Song-You Hong (Yonsei Univ.) and Jimy Dudhia (NCAR)
+! Summer 2003
+!
+! Reference) Hong, Dudhia, Chen (HDC, 2004) Mon. Wea. Rev.
+! Rutledge, Hobbs (RH83, 1983) J. Atmos. Sci.
+! Hong and Lim (HL, 2006) J. Korean Meteor. Soc.
+!
+ INTEGER TENTIN :: ids,ide, jds,jde, kds,kde , &
+ ims,ime, jms,jme, kms,kme , &
+ its,ite, jts,jte, kts,kte
+ REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) TENTINOUT :: &
+ th, &
+ q, &
+ qc, &
+ qi, &
+ qr, &
+ qs
+ REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) TENTIN :: &
+ den, &
+ pii, &
+ p, &
+ delz
+ REAL TENTIN :: &
+ delt, &
+ g, &
+ rd, &
+ rv, &
+ t0c, &
+ den0, &
+ cpd, &
+ cpv, &
+ ep1, &
+ ep2, &
+ qmin, &
+ XLS, &
+ XLV0, &
+ XLF0, &
+ cliq, &
+ cice, &
+ psat, &
+ denr
+ REAL, DIMENSION( ims:ime , jms:jme ) TENTINOUT :: &
+ rain, &
+ rainncv, &
+ sr
+
+ REAL, DIMENSION( ims:ime , jms:jme ) OPTY TENTINOUT :: &
+ snow, &
+ snowncv
+
+! LOCAL VAR
+ INTEGER :: ids2,ide2, jds2,jde2, kds2,kde2
+ REAL, DIMENSION( its:ite , kts:kte ) :: t
+ REAL, DIMENSION( its:ite , kts:kte, 2 ) :: qci, qrs
+ INTEGER :: i,j,k
+
+ REAL*4 delt_s,g_s,rd_s,rv_s,t0c_s,den0_s,cpd_s,cpv_s,ep1_s,ep2_s
+ REAL*4 qmin_s,XLS_s,XLV0_s,XLF0_s,cliq_s,cice_s,psat_s,denr_s
+ REAL*4 th_s(its:ite,kts:kte,jts:jte)
+ REAL*4 pii_s(its:ite,kts:kte,jts:jte)
+ REAL*4 q_s(its:ite,kts:kte,jts:jte)
+ REAL*4 qc_s(its:ite,kts:kte,jts:jte)
+ REAL*4 qi_s(its:ite,kts:kte,jts:jte)
+ REAL*4 qr_s(its:ite,kts:kte,jts:jte)
+ REAL*4 qs_s(its:ite,kts:kte,jts:jte)
+ REAL*4 den_s(its:ite,kts:kte,jts:jte)
+ REAL*4 p_s(its:ite,kts:kte,jts:jte)
+ REAL*4 delz_s(its:ite,kts:kte,jts:jte)
+ REAL*4 rain_s(its:ite,jts:jte)
+ REAL*4 rainncv_s(its:ite,jts:jte)
+ REAL*4 sr_s(its:ite,jts:jte)
+ REAL*4 snow_s(its:ite,jts:jte)
+ REAL*4 snowncv_s(its:ite,jts:jte)
+
+!-------------------------------------------------------------------
+CHARACTER*80 fname
+#ifdef STANDALONE
+INTEGER :: thisstep
+CHARACTER*80 :: pathtofile
+COMMON /wsm5_driver_block/ thisstep, pathtofile
+INTEGER s, e
+INTEGER, EXTERNAL :: RSL_INTERNAL_MICROCLOCK
+#else
+INTEGER, SAVE :: thisstep = 0
+#endif
+
+!INOUT 3D
+
+thisstep = thisstep + 1
+
+
+#ifdef OUTPUT_SNAPSHOTS
+write(fname,'(A,"wsm5_in_",i3.3)')trim(pathtofile),thisstep
+open(45,file=fname,form='UNFORMATTED')
+write(45) ids,ide,jds,jde,kds,kde,ims,ime,jms,jme,kms,kme,its,ite,jts,jte,kts,kte
+!IN 0D
+write(45) delt_s,g_s,rd_s,rv_s,t0c_s,den0_s,cpd_s,cpv_s,ep1_s,ep2_s,qmin_s,XLS_s,XLV0_s,XLF0_s,cliq_s,cice_s,psat_s,denr_s
+!INOUT 3D
+write(45) th_s(its:ite,kts:kte,jts:jte)
+write(45) q_s(its:ite,kts:kte,jts:jte)
+write(45) qc_s(its:ite,kts:kte,jts:jte)
+write(45) qr_s(its:ite,kts:kte,jts:jte)
+write(45) qi_s(its:ite,kts:kte,jts:jte)
+write(45) qs_s(its:ite,kts:kte,jts:jte)
+!INOUT 2D
+write(45) rain_s(its:ite,jts:jte)
+write(45) rainncv_s(its:ite,jts:jte)
+write(45) sr_s(its:ite,jts:jte)
+write(45) snow_s(its:ite,jts:jte)
+write(45) snowncv_s(its:ite,jts:jte)
+!IN 3D
+write(45) den_s(its:ite,kts:kte,jts:jte)
+write(45) pii_s(its:ite,kts:kte,jts:jte)
+write(45) p_s(its:ite,kts:kte,jts:jte)
+write(45) delz_s(its:ite,kts:kte,jts:jte)
+close(45)
+#endif
+
+#if defined (STANDALONE)
+write(fname,'(A,"wsm5_in_",i3.3)')trim(pathtofile),thisstep
+open(45,file=fname,form='UNFORMATTED')
+read(45) ids,ide,jds,jde,kds,kde,ims,ime,jms,jme,kms,kme,its,ite,jts,jte,kts,kte
+!IN 0D
+read(45) delt_s,g_s,rd_s,rv_s,t0c_s,den0_s,cpd_s,cpv_s,ep1_s,ep2_s,qmin_s,XLS_s,XLV0_s,XLF0_s,cliq_s,cice_s,psat_s,denr_s
+!INOUT 3D
+read(45) th_s(its:ite,kts:kte,jts:jte)
+read(45) q_s(its:ite,kts:kte,jts:jte)
+read(45) qc_s(its:ite,kts:kte,jts:jte)
+read(45) qr_s(its:ite,kts:kte,jts:jte)
+read(45) qi_s(its:ite,kts:kte,jts:jte)
+read(45) qs_s(its:ite,kts:kte,jts:jte)
+!INOUT 2D
+read(45) rain_s(its:ite,jts:jte)
+read(45) rainncv_s(its:ite,jts:jte)
+read(45) sr_s(its:ite,jts:jte)
+read(45) snow_s(its:ite,jts:jte)
+read(45) snowncv_s(its:ite,jts:jte)
+!IN 3D
+read(45) den_s(its:ite,kts:kte,jts:jte)
+read(45) pii_s(its:ite,kts:kte,jts:jte)
+read(45) p_s(its:ite,kts:kte,jts:jte)
+read(45) delz_s(its:ite,kts:kte,jts:jte)
+close(45)
+#endif
+
+#ifdef DIAGS
+open(80,file='th0',form='FORMATTED')
+open(81,file='q0',form='FORMATTED')
+open(82,file='qc0',form='FORMATTED')
+open(83,file='qr0',form='FORMATTED')
+open(84,file='qi0',form='FORMATTED')
+open(85,file='qs0',form='FORMATTED')
+do k = kts,kte
+write(80,*)ite-its+1,jte-jts+1,' k '
+write(81,*)ite-its+1,jte-jts+1,' k '
+write(82,*)ite-its+1,jte-jts+1,' k '
+write(83,*)ite-its+1,jte-jts+1,' k '
+write(84,*)ite-its+1,jte-jts+1,' k '
+write(85,*)ite-its+1,jte-jts+1,' k '
+do j = jts,jte
+do i = its,ite
+write(80,*)th(i,k,j)
+write(81,*)q(i,k,j)
+write(82,*)qc(i,k,j)
+write(83,*)qr(i,k,j)
+write(84,*)qi(i,k,j)
+write(85,*)qs(i,k,j)
+enddo
+enddo
+enddo
+close(80)
+close(81)
+close(82)
+close(83)
+close(84)
+close(85)
+#endif
+
+#if !defined(INPUT_SNAPSHOTS) || defined (STANDALONE)
+#ifdef STANDALONE
+#endif
+
+#include "s2d.h"
+
+#ifdef RUN_ON_GPU
+
+ CALL wsm5_gpu_init(0,1,0)
+
+s = RSL_INTERNAL_MICROCLOCK()
+
+write(0,*)'calling wsm5_host',its,ite,jts,jte,kts,kte
+ CALL wsm5_host ( &
+ th(its:ite,kts:kte,jts:jte), pii(its:ite,kts:kte,jts:jte) &
+ ,q(its:ite,kts:kte,jts:jte), qc(its:ite,kts:kte,jts:jte) &
+ ,qi(its:ite,kts:kte,jts:jte), qr(its:ite,kts:kte,jts:jte) &
+ ,qs(its:ite,kts:kte,jts:jte), den(its:ite,kts:kte,jts:jte) &
+ ,p(its:ite,kts:kte,jts:jte), delz(its:ite,kts:kte,jts:jte) &
+ ,delt &
+ ,rain(its:ite,jts:jte),rainncv(its:ite,jts:jte) &
+ ,sr(its:ite,jts:jte) &
+ ,snow(its:ite,jts:jte),snowncv(its:ite,jts:jte) &
+ ,its, ite, jts, jte, kts, kte &
+ ,its, ite, jts, jte, kts, kte &
+ ,its, ite, jts, jte, kts, kte &
+ )
+write(0,*)'back from wsm5_host'
+
+
+e = RSL_INTERNAL_MICROCLOCK()
+#else
+
+s = RSL_INTERNAL_MICROCLOCK()
+ DO j=jts,jte
+ DO k=kts,kte
+ DO i=its,ite
+ t(i,k)=th(i,k,j)*pii(i,k,j)
+ qci(i,k,1) = qc(i,k,j)
+ qci(i,k,2) = qi(i,k,j)
+ qrs(i,k,1) = qr(i,k,j)
+ qrs(i,k,2) = qs(i,k,j)
+ ENDDO
+ ENDDO
+ CALL wsm52D(t, q(ims,kms,j), qci, qrs &
+ ,den(ims,kms,j) &
+ ,p(ims,kms,j), delz(ims,kms,j) &
+ ,delt,g, cpd, cpv, rd, rv, t0c &
+ ,ep1, ep2, qmin &
+ ,XLS, XLV0, XLF0, den0, denr &
+ ,cliq,cice,psat &
+ ,j &
+ ,rain(ims,j),rainncv(ims,j) &
+ ,sr(ims,j) &
+ ,ids,ide, jds,jde, kds,kde &
+ ,ims,ime, jms,jme, kms,kme &
+ ,its,ite, jts,jte, kts,kte &
+ ,snow(ims,j),snowncv(ims,j) &
+ )
+ DO K=kts,kte
+ DO I=its,ite
+ th(i,k,j)=t(i,k)/pii(i,k,j)
+ qc(i,k,j) = qci(i,k,1)
+ qi(i,k,j) = qci(i,k,2)
+ qr(i,k,j) = qrs(i,k,1)
+ qs(i,k,j) = qrs(i,k,2)
+ ENDDO
+ ENDDO
+ ENDDO
+e = RSL_INTERNAL_MICROCLOCK()
+
+#endif
+
+#include "d2s.h"
+
+#ifdef STANDALONE
+write(0,*)'Call to WSM5 on host for step ',thisstep,' is ',e-s,'microseconds'
+#endif
+#else
+
+write(fname,'(A,"wsm5_out_",i3.3)')trim(pathtofile),thisstep
+open(46,file=fname,form='UNFORMATTED')
+!INOUT 3D
+read(46) ids,ide,jds,jde,kds,kde,ims,ime,jms,jme,kms,kme,its,ite,jts,jte,kts,kte
+read(46) th_s(its:ite,kts:kte,jts:jte)
+read(46) q_s(its:ite,kts:kte,jts:jte)
+read(46) qc_s(its:ite,kts:kte,jts:jte)
+read(46) qr_s(its:ite,kts:kte,jts:jte)
+read(46) qi_s(its:ite,kts:kte,jts:jte)
+read(46) qs_s(its:ite,kts:kte,jts:jte)
+!INOUT 2D
+read(46) rain_s(its:ite,jts:jte)
+read(46) rainncv_s(its:ite,jts:jte)
+read(46) sr_s(its:ite,jts:jte)
+read(46) snow_s(its:ite,jts:jte)
+read(46) snowncv_s(its:ite,jts:jte)
+close(46)
+#endif
+
+#if defined(OUTPUT_SNAPSHOTS) || defined(STANDALONE)
+write(fname,'(A,"wsm5_out_",i3.3)')trim(pathtofile),thisstep
+# if defined(STANDALONE)
+# if ( defined(RUN_ON_GPU) )
+write(fname,'("snap_gpu_",i3.3)')thisstep
+# else
+write(fname,'("snap_out_",i3.3)')thisstep
+# endif
+# endif
+open(46,file=fname,form='UNFORMATTED')
+!INOUT 3D
+write(46) ids,ide,jds,jde,kds,kde,ims,ime,jms,jme,kms,kme,its,ite,jts,jte,kts,kte
+write(46) th_s(its:ite,kts:kte,jts:jte)
+write(46) q_s(its:ite,kts:kte,jts:jte)
+write(46) qc_s(its:ite,kts:kte,jts:jte)
+write(46) qr_s(its:ite,kts:kte,jts:jte)
+write(46) qi_s(its:ite,kts:kte,jts:jte)
+write(46) qs_s(its:ite,kts:kte,jts:jte)
+!INOUT 2D
+write(46) rain_s(its:ite,jts:jte)
+write(46) rainncv_s(its:ite,jts:jte)
+write(46) sr_s(its:ite,jts:jte)
+write(46) snow_s(its:ite,jts:jte)
+write(46) snowncv_s(its:ite,jts:jte)
+close(46)
+#endif
+
+#ifdef DIAGS
+open(80,file='th',form='FORMATTED')
+open(81,file='q',form='FORMATTED')
+open(82,file='qc',form='FORMATTED')
+open(83,file='qr',form='FORMATTED')
+open(84,file='qi',form='FORMATTED')
+open(85,file='qs',form='FORMATTED')
+do k = kts,kte
+write(80,*)ite-its+1,jte-jts+1,' k '
+write(81,*)ite-its+1,jte-jts+1,' k '
+write(82,*)ite-its+1,jte-jts+1,' k '
+write(83,*)ite-its+1,jte-jts+1,' k '
+write(84,*)ite-its+1,jte-jts+1,' k '
+write(85,*)ite-its+1,jte-jts+1,' k '
+do j = jts,jte
+do i = its,ite
+write(80,*)th(i,k,j)
+write(81,*)q(i,k,j)
+write(82,*)qc(i,k,j)
+write(83,*)qr(i,k,j)
+write(84,*)qi(i,k,j)
+write(85,*)qs(i,k,j)
+enddo
+enddo
+enddo
+close(80)
+close(81)
+close(82)
+close(83)
+close(84)
+close(85)
+#endif
+
+ END SUBROUTINE wsm5
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!! ORIGINAL SUBROUTINE IN WRF (WITH A FEW CHANGES FOR STANDALONE !!!!
+!!!! AND SOME HOOKS FOR VERIFICATION WITH RESPECT TO BASELINE) !!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+!===================================================================
+!
+ SUBROUTINE wsm52D(t, q, qci, qrs, den, p, delz &
+ ,delt,g, cpd, cpv, rd, rv, t0c &
+ ,ep1, ep2, qmin &
+ ,XLS, XLV0, XLF0, den0, denr &
+ ,cliq,cice,psat &
+ ,lat &
+ ,rain,rainncv &
+ ,sr &
+ ,ids,ide, jds,jde, kds,kde &
+ ,ims,ime, jms,jme, kms,kme &
+ ,its,ite, jts,jte, kts,kte &
+ ,snow,snowncv &
+ )
+!-------------------------------------------------------------------
+ IMPLICIT NONE
+!-------------------------------------------------------------------
+ INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
+ ims,ime, jms,jme, kms,kme , &
+ its,ite, jts,jte, kts,kte, &
+ lat
+ REAL, DIMENSION( its:ite , kts:kte ), &
+ INTENT(INOUT) :: &
+ t
+ REAL, DIMENSION( its:ite , kts:kte, 2 ), &
+ INTENT(INOUT) :: &
+ qci, &
+ qrs
+
+ REAL, DIMENSION( ims:ime , kms:kme ), &
+ INTENT(INOUT) :: &
+ q
+ REAL, DIMENSION( ims:ime , kms:kme ), &
+ INTENT(IN ) :: &
+ den, &
+ p, &
+ delz
+ REAL, INTENT(IN ) :: delt, &
+ g, &
+ cpd, &
+ cpv, &
+ t0c, &
+ den0, &
+ rd, &
+ rv, &
+ ep1, &
+ ep2, &
+ qmin, &
+ XLS, &
+ XLV0, &
+ XLF0, &
+ cliq, &
+ cice, &
+ psat, &
+ denr
+ REAL, DIMENSION( ims:ime ), &
+ INTENT(INOUT) :: rain, &
+ rainncv, &
+ sr
+
+ REAL, DIMENSION( ims:ime ), OPTIONAL, &
+ INTENT(INOUT) :: snow, &
+ snowncv
+
+! LOCAL VAR
+ REAL, DIMENSION( its:ite , kts:kte , 2) :: &
+ rh, qs, rslope, rslope2, rslope3, rslopeb, &
+ falk, fall, work1
+ REAL, DIMENSION( its:ite , kts:kte ) :: &
+ falkc, work1c, work2c, fallc
+ REAL, DIMENSION( its:ite , kts:kte ) :: &
+ praut, psaut, prevp, psdep, pracw, psaci, psacw, &
+ pigen, pidep, pcond, xl, cpm, work2, psmlt, psevp, denfac, xni,&
+ n0sfac
+! variables for optimization
+ REAL, DIMENSION( its:ite ) :: tvec1
+ INTEGER, DIMENSION( its:ite ) :: mstep, numdt
+ REAL, DIMENSION(its:ite) :: rmstep
+ REAL dtcldden, rdelz, rdtcld
+ LOGICAL, DIMENSION( its:ite ) :: flgcld
+ REAL :: pi, &
+ cpmcal, xlcal, lamdar, lamdas, diffus, &
+ viscos, xka, venfac, conden, diffac, &
+ x, y, z, a, b, c, d, e, &
+ qdt, holdrr, holdrs, supcol, pvt, &
+ coeres, supsat, dtcld, xmi, eacrs, satdt, &
+ vt2i,vt2s,acrfac, &
+ qimax, diameter, xni0, roqi0, &
+ fallsum, fallsum_qsi, xlwork2, factor, source, &
+ value, xlf, pfrzdtc, pfrzdtr, supice
+ REAL :: temp
+ REAL :: holdc, holdci
+ INTEGER :: i, j, k, mstepmax, &
+ iprt, latd, lond, loop, loops, ifsat, n
+! Temporaries used for inlining fpvs function
+ REAL :: dldti, xb, xai, tr, xbi, xa, hvap, cvap, hsub, dldt, ttp
+!
+!=================================================================
+! compute internal functions
+!
+ cpmcal(x) = cpd*(1.-max(x,qmin))+max(x,qmin)*cpv
+ xlcal(x) = xlv0-xlv1*(x-t0c)
+!----------------------------------------------------------------
+! size distributions: (x=mixing ratio, y=air density):
+! valid for mixing ratio > 1.e-9 kg/kg.
+!
+! Optimizatin : A**B => exp(log(A)*(B))
+ lamdar(x,y)= sqrt(sqrt(pidn0r/(x*y))) ! (pidn0r/(x*y))**.25
+ lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y))) ! (pidn0s*z/(x*y))**.25
+!
+!----------------------------------------------------------------
+! diffus: diffusion coefficient of the water vapor
+! viscos: kinematic viscosity(m2s-1)
+ diffus(x,y) = 8.794d-5 * exp(log(x)*(1.81d0)) / y ! 8.794e-5*x**1.81/y
+ viscos(x,y) = 1.496d-6 * (x*sqrt(x)) /(x+120.d0)/y ! 1.496e-6*x**1.5/(x+120.)/y
+ xka(x,y) = 1.414d3*viscos(x,y)*y
+ diffac(a,b,c,d,e) = d*a*a/(xka(c,d)*rv*c*c)+1./(e*diffus(c,b))
+ venfac(a,b,c) = exp(log((viscos(b,c)/diffus(b,a)))*((.3333333d0))) &
+ /sqrt(viscos(b,c))*sqrt(sqrt(den0/c))
+ conden(a,b,c,d,e) = (max(b,qmin)-c)/(1.+d*d/(rv*e)*c/(a*a))
+!
+!
+ pi = 4. * atan(1.)
+!
+!----------------------------------------------------------------
+! paddint 0 for negative values generated by dynamics
+!
+ do k = kts, kte
+ do i = its, ite
+ qci(i,k,1) = max(qci(i,k,1),0.0)
+ qrs(i,k,1) = max(qrs(i,k,1),0.0)
+ qci(i,k,2) = max(qci(i,k,2),0.0)
+ qrs(i,k,2) = max(qrs(i,k,2),0.0)
+ enddo
+ enddo
+!
+!----------------------------------------------------------------
+! latent heat for phase changes and heat capacity. neglect the
+! changes during microphysical process calculation
+! emanuel(1994)
+!
+
+
+!#include "con.h"
+
+ do k = kts, kte
+ do i = its, ite
+DIAGOUTPUT1("q_start",q(i,k))
+ cpm(i,k) = cpmcal(q(i,k))
+ xl(i,k) = xlcal(t(i,k))
+ enddo
+ enddo
+!
+!----------------------------------------------------------------
+! compute the minor time steps.
+!
+ loops = max(nint(delt/dtcldcr),1)
+ dtcld = delt/loops
+ if(delt.le.dtcldcr) dtcld = delt
+!
+ do loop = 1,loops
+
+! THESE LOOPS ARE ADDED TO OUTPUT THE MAIN INPUT FIELDS
+! AT THE START OF THE MAIN LOOP IN WSM5. THE MACROS ARE
+! ONLY ENABLED
+
+do i = its, ite
+do k = kts, kte
+kDIAGOUTPUT1("t_loop",t(i,k))
+enddo
+do k = kts, kte
+kDIAGOUTPUT1("q_loop",q(i,k))
+enddo
+do k = kts, kte
+kDIAGOUTPUT1("qc_loop",qci(i,k,1))
+enddo
+do k = kts, kte
+kDIAGOUTPUT1("qi_loop",qci(i,k,2))
+enddo
+do k = kts, kte
+kDIAGOUTPUT1("qr_loop",qrs(i,k,1))
+enddo
+do k = kts, kte
+kDIAGOUTPUT1("qs_loop",qrs(i,k,2))
+enddo
+do k = kts, kte
+kDIAGOUTPUT1("den_loop",den(i,k))
+enddo
+do k = kts, kte
+kDIAGOUTPUT1("p_loop",p(i,k))
+enddo
+do k = kts, kte
+kDIAGOUTPUT1("delz_loop",delz(i,k))
+enddo
+do k = kts, kte
+kDIAGOUTPUT1("cpm_loop",cpm(i,k))
+enddo
+do k = kts, kte
+kDIAGOUTPUT1("xl_loop",xl(i,k))
+enddo
+enddo
+
+!
+!----------------------------------------------------------------
+! initialize the large scale variables
+!
+ do i = its, ite
+ mstep(i) = 1
+ flgcld(i) = .true.
+ enddo
+!
+! do k = kts, kte
+! do i = its, ite
+! denfac(i,k) = sqrt(den0/den(i,k))
+! enddo
+! enddo
+ do k = kts, kte
+ CALL VREC( tvec1(its), den(its,k), ite-its+1)
+ do i = its, ite
+ tvec1(i) = tvec1(i)*den0
+ enddo
+ CALL VSQRT( denfac(its,k), tvec1(its), ite-its+1)
+ enddo
+!
+! Inline expansion for fpvs
+! qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
+! qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
+ hsub = xls
+ hvap = xlv0
+ cvap = cpv
+ ttp=t0c+0.01
+ dldt=cvap-cliq
+ xa=-dldt/rv
+ xb=xa+hvap/(rv*ttp)
+ dldti=cvap-cice
+ xai=-dldti/rv
+ xbi=xai+hsub/(rv*ttp)
+ do k = kts, kte
+ do i = its, ite
+ tr=ttp/t(i,k)
+ qs(i,k,1)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
+ qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
+ qs(i,k,1) = max(qs(i,k,1),qmin)
+ rh(i,k,1) = max(q(i,k) / qs(i,k,1),qmin)
+ if(t(i,k).lt.ttp) then
+ qs(i,k,2)=psat*exp(log(tr)*(xai))*exp(xbi*(1.-tr))
+ else
+ qs(i,k,2)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
+ endif
+ qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
+ qs(i,k,2) = max(qs(i,k,2),qmin)
+ rh(i,k,2) = max(q(i,k) / qs(i,k,2),qmin)
+ enddo
+ enddo
+!
+!----------------------------------------------------------------
+! initialize the variables for microphysical physics
+!
+!
+ do k = kts, kte
+ do i = its, ite
+ prevp(i,k) = 0.
+ psdep(i,k) = 0.
+ praut(i,k) = 0.
+ psaut(i,k) = 0.
+ pracw(i,k) = 0.
+ psaci(i,k) = 0.
+ psacw(i,k) = 0.
+ pigen(i,k) = 0.
+ pidep(i,k) = 0.
+ pcond(i,k) = 0.
+ psmlt(i,k) = 0.
+ psevp(i,k) = 0.
+ falk(i,k,1) = 0.
+ falk(i,k,2) = 0.
+ fall(i,k,1) = 0.
+ fall(i,k,2) = 0.
+ fallc(i,k) = 0.
+ falkc(i,k) = 0.
+ xni(i,k) = 1.e3
+ enddo
+ enddo
+!
+!----------------------------------------------------------------
+! compute the fallout term:
+! first, vertical terminal velosity for minor loops
+!
+ do k = kts, kte
+ do i = its, ite
+ supcol = t0c-t(i,k)
+DIAGOUTPUT1("t0c",t0c)
+DIAGOUTPUT1("supcol",supcol)
+DIAGOUTPUT1("t(i,k)",t(i,k))
+!---------------------------------------------------------------
+! n0s: Intercept parameter for snow [m-4] [HDC 6]
+!---------------------------------------------------------------
+ n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
+ if(qrs(i,k,1).le.qcrmin)then
+ rslope(i,k,1) = rslopermax
+ rslopeb(i,k,1) = rsloperbmax
+ rslope2(i,k,1) = rsloper2max
+ rslope3(i,k,1) = rsloper3max
+ else
+DIAGOUTPUT1("qrs(i,k,2)",qrs(i,k,1))
+DIAGOUTPUT1("den(i,k)",den(i,k))
+ rslope(i,k,1) = 1./lamdar(qrs(i,k,1),den(i,k))
+DIAGOUTPUT1("rslope",rslope(i,k,1))
+ rslopeb(i,k,1) = exp(log(rslope(i,k,1))*(bvtr))
+ rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
+ rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
+ endif
+ if(qrs(i,k,2).le.qcrmin)then
+ rslope(i,k,2) = rslopesmax
+DIAGOUTPUT1("rslope",rslope(i,k,2))
+ rslopeb(i,k,2) = rslopesbmax
+DIAGOUTPUT1("rslopeb",rslopeb(i,k,2))
+ rslope2(i,k,2) = rslopes2max
+ rslope3(i,k,2) = rslopes3max
+ else
+DIAGOUTPUT1("qrs(i,k,2)",qrs(i,k,2))
+DIAGOUTPUT1("den(i,k)",den(i,k))
+DIAGOUTPUT1("n0sfac(i,k)",n0sfac(i,k))
+ rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k),n0sfac(i,k))
+DIAGOUTPUT1("rslope",rslope(i,k,2))
+ rslopeb(i,k,2) = exp(log(rslope(i,k,2))*(bvts))
+DIAGOUTPUT1("rslopeb",rslopeb(i,k,2))
+ rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
+ rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
+ endif
+!-------------------------------------------------------------
+! Ni: ice crystal number concentraiton [HDC 5c]
+!-------------------------------------------------------------
+! xni(i,k) = min(max(5.38e7*(den(i,k) &
+! *max(qci(i,k,2),qmin))**0.75,1.e3),1.e6)
+ temp = (den(i,k)*max(qci(i,k,2),qmin))
+ temp = sqrt(sqrt(temp*temp*temp))
+#ifdef DEBUGDEBUG
+ xni(i,k) = 1000.
+#else
+ xni(i,k) = min(max(5.38d7*temp,1.d3),1.d6)
+#endif
+ enddo
+ enddo
+!
+ mstepmax = 1
+ numdt = 1
+ do k = kte, kts, -1
+ do i = its, ite
+ work1(i,k,1) = pvtr*rslopeb(i,k,1)*denfac(i,k)/delz(i,k)
+ work1(i,k,2) = pvts*rslopeb(i,k,2)*denfac(i,k)/delz(i,k)
+DIAGOUTPUT1("work1_1",work1(i,k,1))
+DIAGOUTPUT1("rslopeb_1",rslopeb(i,k,1))
+DIAGOUTPUT1("work1_2",work1(i,k,2))
+DIAGOUTPUT1("rslopeb_2",rslopeb(i,k,2))
+DIAGOUTPUT1("denfac",denfac(i,k))
+DIAGOUTPUT1("delz",delz(i,k))
+ numdt(i) = max(nint(max(work1(i,k,1),work1(i,k,2))*dtcld+.5),1)
+ if(numdt(i).ge.mstep(i)) mstep(i) = numdt(i)
+ enddo
+ enddo
+#if FLOAT_4 == 4
+ mstep = 4 ! hard code this to match gpu version (must fix this)
+#endif
+ do i = its, ite
+ if(mstepmax.le.mstep(i)) mstepmax = mstep(i)
+ rmstep(i) = 1./mstep(i)
+ enddo
+!
+ do n = 1, mstepmax
+ k = kte
+ do i = its, ite
+ if(n.le.mstep(i)) then
+ falk(i,k,1) = den(i,k)*qrs(i,k,1)*work1(i,k,1)*rmstep(i)
+ falk(i,k,2) = den(i,k)*qrs(i,k,2)*work1(i,k,2)*rmstep(i)
+ fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
+ fall(i,k,2) = fall(i,k,2)+falk(i,k,2)
+ dtcldden = dtcld/den(i,k)
+ qrs(i,k,1) = max(qrs(i,k,1)-falk(i,k,1)*dtcldden,0.)
+ qrs(i,k,2) = max(qrs(i,k,2)-falk(i,k,2)*dtcldden,0.)
+ endif
+ enddo
+ do k = kte-1, kts, -1
+ do i = its, ite
+ if(n.le.mstep(i)) then
+ falk(i,k,1) = den(i,k)*qrs(i,k,1)*work1(i,k,1)*rmstep(i)
+ falk(i,k,2) = den(i,k)*qrs(i,k,2)*work1(i,k,2)*rmstep(i)
+ fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
+ fall(i,k,2) = fall(i,k,2)+falk(i,k,2)
+ dtcldden = dtcld/den(i,k)
+ rdelz = 1./delz(i,k)
+DIAGOUTPUT1i("loop",loop)
+DIAGOUTPUT1i("mstepi",mstep(i))
+DIAGOUTPUT1i("n",n)
+DIAGOUTPUT1("qr_090",qrs(i,k,1))
+DIAGOUTPUT1("falk",falk(i,k,1))
+DIAGOUTPUT1("falk1",falk(i,k+1,1))
+DIAGOUTPUT1("delz",delz(i,k))
+DIAGOUTPUT1("delz1",delz(i,k+1))
+ qrs(i,k,1) = max(qrs(i,k,1)-(falk(i,k,1)-falk(i,k+1,1) &
+ *delz(i,k+1)*rdelz)*dtcldden,0.)
+DIAGOUTPUT1("qr_091",qrs(i,k,1))
+DIAGOUTPUT1("qs_092",qrs(i,k,2))
+DIAGOUTPUT1("falk2",falk(i,k,2))
+DIAGOUTPUT1("work2",work1(i,k,2))
+DIAGOUTPUT1("falk21",falk(i,k+1,2))
+ qrs(i,k,2) = max(qrs(i,k,2)-(falk(i,k,2)-falk(i,k+1,2) &
+ *delz(i,k+1)*rdelz)*dtcldden,0.)
+DIAGOUTPUT1("qs_093",qrs(i,k,2))
+ endif
+ enddo
+ enddo
+ do k = kte, kts, -1
+ do i = its, ite
+ if(n.le.mstep(i)) then
+DIAGOUTPUT1("t",t(i,k))
+DIAGOUTPUT1("qs",qrs(i,k,2))
+ if(t(i,k).gt.t0c.and.qrs(i,k,2).gt.0.) then
+!----------------------------------------------------------------
+! psmlt: melting of snow [HL A33] [RH83 A25]
+! (T>T0: S->R)
+!----------------------------------------------------------------
+ xlf = xlf0
+! work2(i,k)= venfac(p(i,k),t(i,k),den(i,k))
+ work2(i,k)= (exp(log(((1.496e-6*((t(i,k))*sqrt(t(i,k))) &
+ /((t(i,k))+120.)/(den(i,k)))/(8.794e-5 &
+ *exp(log(t(i,k))*(1.81))/p(i,k)))) &
+ *((.3333333)))/sqrt((1.496e-6*((t(i,k)) &
+ *sqrt(t(i,k)))/((t(i,k))+120.)/(den(i,k)))) &
+ *sqrt(sqrt(den0/(den(i,k)))))
+ coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
+! psmlt(i,k) = xka(t(i,k),den(i,k))/xlf*(t0c-t(i,k))*pi/2. &
+! *n0sfac(i,k)*(precs1*rslope2(i,k,2)+precs2 &
+! *work2(i,k)*coeres)
+ psmlt(i,k) = &
+(1.414e3*(1.496e-6 * ((t(i,k))*sqrt(t(i,k))) /((t(i,k))+120.)/(den(i,k)) )*(den(i,k)))&
+ /xlf*(t0c-t(i,k))*pi/2. &
+ *n0sfac(i,k)*(precs1*rslope2(i,k,2)+precs2 &
+ *work2(i,k)*coeres)
+ psmlt(i,k) = min(max(psmlt(i,k)*dtcld/mstep(i), &
+ -qrs(i,k,2)/mstep(i)),0.)
+ qrs(i,k,2) = qrs(i,k,2) + psmlt(i,k)
+DIAGOUTPUT1i("mstepi",mstep(i))
+DIAGOUTPUT1i("n",n)
+DIAGOUTPUT1("qr_100",qrs(i,k,1))
+DIAGOUTPUT1("psmlt",psmlt(i,k))
+ qrs(i,k,1) = qrs(i,k,1) - psmlt(i,k)
+DIAGOUTPUT1("qr_101",qrs(i,k,1))
+ t(i,k) = t(i,k) + xlf/cpm(i,k)*psmlt(i,k)
+ endif
+ endif
+ enddo
+ enddo
+ enddo
+
+#if FLOAT_4 != 4
+
+!---------------------------------------------------------------
+! Vice [ms-1] : fallout of ice crystal [HDC 5a]
+!---------------------------------------------------------------
+ mstepmax = 1
+ mstep = 1
+ numdt = 1
+ do k = kte, kts, -1
+ do i = its, ite
+ if(qci(i,k,2).le.0.) then
+ work2c(i,k) = 0.
+ else
+ xmi = den(i,k)*qci(i,k,2)/xni(i,k)
+! diameter = min(dicon * sqrt(xmi),dimax)
+ diameter = max(min(dicon * sqrt(xmi),dimax), 1.e-25)
+ work1c(i,k) = 1.49e4*exp(log(diameter)*(1.31))
+ work2c(i,k) = work1c(i,k)/delz(i,k)
+ endif
+ numdt(i) = max(nint(work2c(i,k)*dtcld+.5),1)
+ if(numdt(i).ge.mstep(i)) mstep(i) = numdt(i)
+ enddo
+ enddo
+ do i = its, ite
+ if(mstepmax.le.mstep(i)) mstepmax = mstep(i)
+ enddo
+!
+ do n = 1, mstepmax
+ k = kte
+ do i = its, ite
+ if(n.le.mstep(i)) then
+ falkc(i,k) = den(i,k)*qci(i,k,2)*work2c(i,k)/mstep(i)
+ fallc(i,k) = fallc(i,k)+falkc(i,k)
+ qci(i,k,2) = max(qci(i,k,2)-falkc(i,k)*dtcld/den(i,k),0.)
+ endif
+ enddo
+ do k = kte-1, kts, -1
+ do i = its, ite
+ if(n.le.mstep(i)) then
+ falkc(i,k) = den(i,k)*qci(i,k,2)*work2c(i,k)/mstep(i)
+ fallc(i,k) = fallc(i,k)+falkc(i,k)
+ qci(i,k,2) = max(qci(i,k,2)-(falkc(i,k)-falkc(i,k+1) &
+ *delz(i,k+1)/delz(i,k))*dtcld/den(i,k),0.)
+ endif
+ enddo
+ enddo
+ enddo
+!
+!
+!----------------------------------------------------------------
+! rain (unit is mm/sec;kgm-2s-1: /1000*delt ===> m)==> mm for wrf
+!
+ do i = its, ite
+ fallsum = fall(i,1,1)+fall(i,1,2)+fallc(i,1)
+ fallsum_qsi = fall(i,1,2)+fallc(i,1)
+ rainncv(i) = 0.
+ if(fallsum.gt.0.) then
+ rainncv(i) = fallsum*delz(i,1)/denr*dtcld*1000.
+ rain(i) = fallsum*delz(i,1)/denr*dtcld*1000. + rain(i)
+ endif
+ IF ( PRESENT (snowncv) .AND. PRESENT (snow)) THEN
+ snowncv(i) = 0.
+ if(fallsum_qsi.gt.0.) then
+ snowncv(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000.
+ snow(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. + snow(i)
+ endif
+ ENDIF
+ sr(i) = 0.
+ if(fallsum.gt.0.)sr(i)=fallsum_qsi*delz(i,kts)/denr*dtcld*1000./(rainncv(i)+1.e-12)
+ enddo
+!
+!---------------------------------------------------------------
+! pimlt: instantaneous melting of cloud ice [HL A47] [RH83 A28]
+! (T>T0: I->C)
+!---------------------------------------------------------------
+ do k = kts, kte
+ do i = its, ite
+ supcol = t0c-t(i,k)
+ xlf = xls-xl(i,k)
+ if(supcol.lt.0.) xlf = xlf0
+ if(supcol.lt.0.and.qci(i,k,2).gt.0.) then
+ qci(i,k,1) = qci(i,k,1) + qci(i,k,2)
+ t(i,k) = t(i,k) - xlf/cpm(i,k)*qci(i,k,2)
+ qci(i,k,2) = 0.
+ endif
+!---------------------------------------------------------------
+! pihmf: homogeneous freezing of cloud water below -40c [HL A45]
+! (T<-40C: C->I)
+!---------------------------------------------------------------
+ if(supcol.gt.40..and.qci(i,k,1).gt.0.) then
+ qci(i,k,2) = qci(i,k,2) + qci(i,k,1)
+ t(i,k) = t(i,k) + xlf/cpm(i,k)*qci(i,k,1)
+ qci(i,k,1) = 0.
+ endif
+!---------------------------------------------------------------
+! pihtf: heterogeneous freezing of cloud water [HL A44]
+! (T0>T>-40C: C->I)
+!---------------------------------------------------------------
+ if(supcol.gt.0..and.qci(i,k,1).gt.0.) then
+! pfrzdtc = min(pfrz1*(exp(pfrz2*supcol)-1.) &
+! *den(i,k)/denr/xncr*qci(i,k,1)**2*dtcld,qci(i,k,1))
+ pfrzdtc = min(pfrz1*(exp(pfrz2*supcol)-1.) &
+ *den(i,k)/denr/xncr*qci(i,k,1)*qci(i,k,1)*dtcld,qci(i,k,1))
+ qci(i,k,2) = qci(i,k,2) + pfrzdtc
+ t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtc
+ qci(i,k,1) = qci(i,k,1)-pfrzdtc
+ endif
+!---------------------------------------------------------------
+! psfrz: freezing of rain water [HL A20] [LFO 45]
+! (T<T0, R->S)
+!---------------------------------------------------------------
+ if(supcol.gt.0..and.qrs(i,k,1).gt.0.) then
+! pfrzdtr = min(20.*pi**2*pfrz1*n0r*denr/den(i,k) &
+! *(exp(pfrz2*supcol)-1.)*rslope(i,k,1)**7*dtcld, &
+! qrs(i,k,1))
+ temp = rslope(i,k,1)
+ temp = temp*temp*temp*temp*temp*temp*temp
+ pfrzdtr = min(20.*(pi*pi)*pfrz1*n0r*denr/den(i,k) &
+ *(exp(pfrz2*supcol)-1.)*temp*dtcld, &
+ qrs(i,k,1))
+ qrs(i,k,2) = qrs(i,k,2) + pfrzdtr
+ t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtr
+ qrs(i,k,1) = qrs(i,k,1)-pfrzdtr
+DIAGOUTPUT1("qr_BBB",qrs(i,k,1))
+ endif
+ enddo
+ enddo
+!
+!----------------------------------------------------------------
+! rsloper: reverse of the slope parameter of the rain(m)
+! xka: thermal conductivity of air(jm-1s-1k-1)
+! work1: the thermodynamic term in the denominator associated with
+! heat conduction and vapor diffusion
+! (ry88, y93, h85)
+! work2: parameter associated with the ventilation effects(y93)
+!
+ do k = kts, kte
+ do i = its, ite
+ if(qrs(i,k,1).le.qcrmin)then
+ rslope(i,k,1) = rslopermax
+ rslopeb(i,k,1) = rsloperbmax
+ rslope2(i,k,1) = rsloper2max
+ rslope3(i,k,1) = rsloper3max
+ else
+! rslope(i,k,1) = 1./lamdar(qrs(i,k,1),den(i,k))
+ rslope(i,k,1) = 1./(sqrt(sqrt(pidn0r/((qrs(i,k,1))*(den(i,k))))))
+DIAGOUTPUT1("rslope_set",rslope(i,k,1))
+DIAGOUTPUT1("rslope_qr",qrs(i,k,1))
+DIAGOUTPUT1("den",den(i,k))
+ rslopeb(i,k,1) = exp(log(rslope(i,k,1))*(bvtr))
+ rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
+ rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
+ endif
+ if(qrs(i,k,2).le.qcrmin)then
+ rslope(i,k,2) = rslopesmax
+ rslopeb(i,k,2) = rslopesbmax
+ rslope2(i,k,2) = rslopes2max
+ rslope3(i,k,2) = rslopes3max
+ else
+! rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k),n0sfac(i,k))
+ rslope(i,k,2) = 1./(sqrt(sqrt(pidn0s*(n0sfac(i,k))/((qrs(i,k,2))*(den(i,k))))))
+ rslopeb(i,k,2) = exp(log(rslope(i,k,2))*(bvts))
+ rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
+ rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
+ endif
+ enddo
+ enddo
+!
+ do k = kts, kte
+ do i = its, ite
+ work1(i,k,1) = diffac(xl(i,k),p(i,k),t(i,k),den(i,k),qs(i,k,1))
+! work1(i,k,1) = &
+! ((((den(i,k))*(xl(i,k))*(xl(i,k))) * ((t(i,k))+120.) * (den(i,k))) &
+! / &
+! ( 1.414e3 * (1.496e-6 * ((t(i,k))*sqrt(t(i,k)))) * (den(i,k)) * &
+! (rv*(t(i,k))*(t(i,k))))) &
+! + &
+! p(i,k) / ( (qs(i,k,1)) * ( 8.794e-5 * exp(log(t(i,k))*(1.81)) ) )
+ work1(i,k,2) = diffac(xls,p(i,k),t(i,k),den(i,k),qs(i,k,2))
+! work1(i,k,2) = &
+! ( &
+! (((den(i,k))*(xls)*(xls))*((t(i,k))+120.)*(den(i,k))) &
+! / &
+! ( &
+! 1.414e3 * (1.496e-6 * ((t(i,k))*sqrt(t(i,k)))) * (den(i,k)) * &
+! (rv*(t(i,k))*(t(i,k))) &
+! ) &
+! + &
+! p(i,k) &
+! / &
+! ( qs(i,k,2) * (8.794e-5 * exp(log(t(i,k))*(1.81)))) &
+! )
+ work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
+! work2(i,k) = &
+! ( &
+! exp(.3333333*log( &
+! ((1.496e-6 * ((t(i,k))*sqrt(t(i,k))))*p(i,k)) &
+! / &
+! (((t(i,k))+120.)*den(i,k)*(8.794e-5 * exp(log(t(i,k))*(1.81)))) &
+! )) &
+! * &
+! sqrt(sqrt(den0/(den(i,k)))) &
+! ) &
+! / &
+! sqrt( &
+! (1.496e-6 * ((t(i,k))*sqrt(t(i,k)))) &
+! / &
+! ( &
+! ((t(i,k))+120.) * den(i,k) &
+! ) &
+! )
+ ENDDO
+ ENDDO
+
+!
+!===============================================================
+!
+! warm rain processes
+!
+! - follows the processes in RH83 and LFO except for autoconcersion
+!
+!===============================================================
+!
+ do k = kts, kte
+ do i = its, ite
+ supsat = max(q(i,k),qmin)-qs(i,k,1)
+ satdt = supsat/dtcld
+!---------------------------------------------------------------
+! praut: auto conversion rate from cloud to rain [HDC 16]
+! (C->R)
+!---------------------------------------------------------------
+ if(qci(i,k,1).gt.qc0) then
+ praut(i,k) = qck1*exp(log(qci(i,k,1))*((7./3.)))
+ praut(i,k) = min(praut(i,k),qci(i,k,1)/dtcld)
+ endif
+!---------------------------------------------------------------
+! pracw: accretion of cloud water by rain [HL A40] [LFO 51]
+! (C->R)
+!---------------------------------------------------------------
+ if(qrs(i,k,1).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
+ pracw(i,k) = min(pacrr*rslope3(i,k,1)*rslopeb(i,k,1) &
+ *qci(i,k,1)*denfac(i,k),qci(i,k,1)/dtcld)
+ endif
+!---------------------------------------------------------------
+! prevp: evaporation/condensation rate of rain [HDC 14]
+! (V->R or R->V)
+!---------------------------------------------------------------
+ if(qrs(i,k,1).gt.0.) then
+ coeres = rslope2(i,k,1)*sqrt(rslope(i,k,1)*rslopeb(i,k,1))
+ prevp(i,k) = (rh(i,k,1)-1.)*(precr1*rslope2(i,k,1) &
+ +precr2*work2(i,k)*coeres)/work1(i,k,1)
+DIAGOUTPUT1("prevpset",prevp(i,k))
+DIAGOUTPUT1("prevpset_qrs",qrs(i,k,1))
+DIAGOUTPUT1("prevpset_rslope",rslope(i,k,1))
+DIAGOUTPUT1("prevpset_rslope2",rslope2(i,k,1))
+DIAGOUTPUT1("prevpset_rslopeb",rslopeb(i,k,1))
+DIAGOUTPUT1("prevpset_work1",work1(i,k,1))
+DIAGOUTPUT1("prevpset_rh",rh(i,k,1))
+ if(prevp(i,k).lt.0.) then
+ prevp(i,k) = max(prevp(i,k),-qrs(i,k,1)/dtcld)
+ prevp(i,k) = max(prevp(i,k),satdt/2)
+ else
+ prevp(i,k) = min(prevp(i,k),satdt/2)
+ endif
+ endif
+ enddo
+ enddo
+!
+!===============================================================
+!
+! cold rain processes
+!
+! - follows the revised ice microphysics processes in HDC
+! - the processes same as in RH83 and RH84 and LFO behave
+! following ice crystal hapits defined in HDC, inclduing
+! intercept parameter for snow (n0s), ice crystal number
+! concentration (ni), ice nuclei number concentration
+! (n0i), ice diameter (d)
+!
+!===============================================================
+!
+ rdtcld = 1./dtcld
+ do k = kts, kte
+ do i = its, ite
+ supcol = t0c-t(i,k)
+ supsat = max(q(i,k),qmin)-qs(i,k,2)
+ satdt = supsat/dtcld
+ ifsat = 0
+!-------------------------------------------------------------
+! Ni: ice crystal number concentraiton [HDC 5c]
+!-------------------------------------------------------------
+! xni(i,k) = min(max(5.38e7*(den(i,k) &
+! *max(qci(i,k,2),qmin))**0.75,1.e3),1.e6)
+ temp = (den(i,k)*max(qci(i,k,2),qmin))
+ temp = sqrt(sqrt(temp*temp*temp))
+ xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
+ eacrs = exp(0.07*(-supcol))
+!
+ if(supcol.gt.0) then
+ if(qrs(i,k,2).gt.qcrmin.and.qci(i,k,2).gt.qmin) then
+ xmi = den(i,k)*qci(i,k,2)/xni(i,k)
+ diameter = min(dicon * sqrt(xmi),dimax)
+ vt2i = 1.49e4*diameter**1.31
+ vt2s = pvts*rslopeb(i,k,2)*denfac(i,k)
+!-------------------------------------------------------------
+! psaci: Accretion of cloud ice by rain [HDC 10]
+! (T<T0: I->S)
+!-------------------------------------------------------------
+ acrfac = 2.*rslope3(i,k,2)+2.*diameter*rslope2(i,k,2) &
+ +diameter**2*rslope(i,k,2)
+ psaci(i,k) = pi*qci(i,k,2)*eacrs*n0s*n0sfac(i,k) &
+ *abs(vt2s-vt2i)*acrfac/4.
+ endif
+ endif
+!-------------------------------------------------------------
+! psacw: Accretion of cloud water by snow [HL A7] [LFO 24]
+! (T<T0: C->S, and T>=T0: C->R)
+!-------------------------------------------------------------
+ if(qrs(i,k,2).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
+ psacw(i,k) = min(pacrc*n0sfac(i,k)*rslope3(i,k,2) &
+ *rslopeb(i,k,2)*qci(i,k,1)*denfac(i,k) &
+! ,qci(i,k,1)/dtcld)
+ ,qci(i,k,1)*rdtcld)
+ endif
+
+ if(supcol .gt. 0) then
+!-------------------------------------------------------------
+! pidep: Deposition/Sublimation rate of ice [HDC 9]
+! (T<T0: V->I or I->V)
+!-------------------------------------------------------------
+ if(qci(i,k,2).gt.0.and.ifsat.ne.1) then
+ xmi = den(i,k)*qci(i,k,2)/xni(i,k)
+ diameter = dicon * sqrt(xmi)
+ pidep(i,k) = 4.*diameter*xni(i,k)*(rh(i,k,2)-1.)/work1(i,k,2)
+ supice = satdt-prevp(i,k)
+ if(pidep(i,k).lt.0.) then
+! pidep(i,k) = max(max(pidep(i,k),satdt/2),supice)
+! pidep(i,k) = max(pidep(i,k),-qci(i,k,2)/dtcld)
+ pidep(i,k) = max(max(pidep(i,k),satdt*.5),supice)
+ pidep(i,k) = max(pidep(i,k),-qci(i,k,2)*rdtcld)
+ else
+! pidep(i,k) = min(min(pidep(i,k),satdt/2),supice)
+ pidep(i,k) = min(min(pidep(i,k),satdt*.5),supice)
+ endif
+ if(abs(prevp(i,k)+pidep(i,k)).ge.abs(satdt)) ifsat = 1
+ endif
+!-------------------------------------------------------------
+! psdep: deposition/sublimation rate of snow [HDC 14]
+! (V->S or S->V)
+!-------------------------------------------------------------
+ if(qrs(i,k,2).gt.0..and.ifsat.ne.1) then
+ coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
+ psdep(i,k) = (rh(i,k,2)-1.)*n0sfac(i,k) &
+ *(precs1*rslope2(i,k,2)+precs2 &
+ *work2(i,k)*coeres)/work1(i,k,2)
+ supice = satdt-prevp(i,k)-pidep(i,k)
+ if(psdep(i,k).lt.0.) then
+! psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)/dtcld)
+! psdep(i,k) = max(max(psdep(i,k),satdt/2),supice)
+ psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)*rdtcld)
+ psdep(i,k) = max(max(psdep(i,k),satdt*.5),supice)
+ else
+! psdep(i,k) = min(min(psdep(i,k),satdt/2),supice)
+ psdep(i,k) = min(min(psdep(i,k),satdt*.5),supice)
+ endif
+ if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)).ge.abs(satdt)) &
+ ifsat = 1
+ endif
+!-------------------------------------------------------------
+! pigen: generation(nucleation) of ice from vapor [HL A50] [HDC 7-8]
+! (T<T0: V->I)
+!-------------------------------------------------------------
+ if(supsat.gt.0.and.ifsat.ne.1) then
+ supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)
+ xni0 = 1.e3*exp(0.1*supcol)
+ roqi0 = 4.92e-11*exp(log(xni0)*(1.33))
+ pigen(i,k) = max(0.,(roqi0/den(i,k)-max(qci(i,k,2),0.)) &
+! /dtcld)
+ *rdtcld)
+ pigen(i,k) = min(min(pigen(i,k),satdt),supice)
+ endif
+!
+!-------------------------------------------------------------
+! psaut: conversion(aggregation) of ice to snow [HDC 12]
+! (T<T0: I->S)
+!-------------------------------------------------------------
+ if(qci(i,k,2).gt.0.) then
+ qimax = roqimax/den(i,k)
+! psaut(i,k) = max(0.,(qci(i,k,2)-qimax)/dtcld)
+ psaut(i,k) = max(0.,(qci(i,k,2)-qimax)*rdtcld)
+ endif
+ endif
+!-------------------------------------------------------------
+! psevp: Evaporation of melting snow [HL A35] [RH83 A27]
+! (T>T0: S->V)
+!-------------------------------------------------------------
+ if(supcol.lt.0.) then
+ if(qrs(i,k,2).gt.0..and.rh(i,k,1).lt.1.) &
+ psevp(i,k) = psdep(i,k)*work1(i,k,2)/work1(i,k,1)
+! psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)/dtcld),0.)
+ psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)*rdtcld),0.)
+ endif
+ enddo
+ enddo
+!
+!
+!----------------------------------------------------------------
+! check mass conservation of generation terms and feedback to the
+! large scale
+!
+ do k = kts, kte
+ do i = its, ite
+ if(t(i,k).le.t0c) then
+!
+! cloud water
+!
+ value = max(qmin,qci(i,k,1))
+ source = (praut(i,k)+pracw(i,k)+psacw(i,k))*dtcld
+ if (source.gt.value) then
+ factor = value/source
+ praut(i,k) = praut(i,k)*factor
+ pracw(i,k) = pracw(i,k)*factor
+ psacw(i,k) = psacw(i,k)*factor
+ endif
+!
+! cloud ice
+!
+ value = max(qmin,qci(i,k,2))
+ source = (psaut(i,k)+psaci(i,k)-pigen(i,k)-pidep(i,k))*dtcld
+ if (source.gt.value) then
+ factor = value/source
+ psaut(i,k) = psaut(i,k)*factor
+ psaci(i,k) = psaci(i,k)*factor
+ pigen(i,k) = pigen(i,k)*factor
+ pidep(i,k) = pidep(i,k)*factor
+ endif
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! added from WRF 3.1.0
+! rain
+!
+!
+ value = max(qmin,qrs(i,k,1))
+ source = (-praut(i,k)-pracw(i,k)-prevp(i,k))*dtcld
+ if (source.gt.value) then
+ factor = value/source
+ praut(i,k) = praut(i,k)*factor
+ pracw(i,k) = pracw(i,k)*factor
+ prevp(i,k) = prevp(i,k)*factor
+ endif
+!
+! snow
+!
+ value = max(qmin,qrs(i,k,2))
+ source = (-psdep(i,k)-psaut(i,k)-psaci(i,k)-psacw(i,k))*dtcld
+ if (source.gt.value) then
+ factor = value/source
+ psdep(i,k) = psdep(i,k)*factor
+ psaut(i,k) = psaut(i,k)*factor
+ psaci(i,k) = psaci(i,k)*factor
+ psacw(i,k) = psacw(i,k)*factor
+ endif
+!
+! end addition from WRF 3.1.0
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+ work2(i,k)=-(prevp(i,k)+psdep(i,k)+pigen(i,k)+pidep(i,k))
+! update
+DIAGOUTPUT1("q_pre-update1",q(i,k))
+DIAGOUTPUT1("prevp1",prevp(i,k))
+DIAGOUTPUT1("psdep1",psdep(i,k))
+DIAGOUTPUT1("pigen1",pigen(i,k))
+DIAGOUTPUT1("pidep1",pidep(i,k))
+ q(i,k) = q(i,k)+work2(i,k)*dtcld
+DIAGOUTPUT1("q_update1",q(i,k))
+ qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k) &
+ +psacw(i,k))*dtcld,0.)
+ qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k) &
+ +prevp(i,k))*dtcld,0.)
+ qci(i,k,2) = max(qci(i,k,2)-(psaut(i,k)+psaci(i,k) &
+ -pigen(i,k)-pidep(i,k))*dtcld,0.)
+DIAGOUTPUT1("qrs_690",qrs(i,k,2))
+ qrs(i,k,2) = max(qrs(i,k,2)+(psdep(i,k)+psaut(i,k) &
+ +psaci(i,k)+psacw(i,k))*dtcld,0.)
+DIAGOUTPUT1("qrs_700",qrs(i,k,2))
+ xlf = xls-xl(i,k)
+ xlwork2 = -xls*(psdep(i,k)+pidep(i,k)+pigen(i,k)) &
+ -xl(i,k)*prevp(i,k)-xlf*psacw(i,k)
+ t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
+ else
+!
+! cloud water
+!
+ value = max(qmin,qci(i,k,1))
+ source=(praut(i,k)+pracw(i,k)+psacw(i,k))*dtcld
+ if (source.gt.value) then
+ factor = value/source
+ praut(i,k) = praut(i,k)*factor
+ pracw(i,k) = pracw(i,k)*factor
+ psacw(i,k) = psacw(i,k)*factor
+ endif
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! added from WRF 3.1.0
+! rain
+!
+ value = max(qmin,qrs(i,k,1))
+ source = (-praut(i,k)-pracw(i,k)-prevp(i,k)-psacw(i,k))*dtcld
+ if (source.gt.value) then
+ factor = value/source
+ praut(i,k) = praut(i,k)*factor
+ pracw(i,k) = pracw(i,k)*factor
+ prevp(i,k) = prevp(i,k)*factor
+ psacw(i,k) = psacw(i,k)*factor
+ endif
+!
+! end addition from WRF 3.1.0
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! snow
+!
+ value = max(qcrmin,qrs(i,k,2))
+ source=(-psevp(i,k))*dtcld
+ if (source.gt.value) then
+ factor = value/source
+ psevp(i,k) = psevp(i,k)*factor
+ endif
+ work2(i,k)=-(prevp(i,k)+psevp(i,k))
+! update
+DIAGOUTPUT1("q_pre-update2",q(i,k))
+DIAGOUTPUT1("prevp2",prevp(i,k))
+DIAGOUTPUT1("psdep2",psdep(i,k))
+DIAGOUTPUT1("pigen2",pigen(i,k))
+DIAGOUTPUT1("pidep2",pidep(i,k))
+ q(i,k) = q(i,k)+work2(i,k)*dtcld
+DIAGOUTPUT1("q_update2",q(i,k))
+ qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k) &
+ +psacw(i,k))*dtcld,0.)
+ qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k) &
+ +prevp(i,k) +psacw(i,k))*dtcld,0.)
+DIAGOUTPUT1("qrs_691",qrs(i,k,2))
+DIAGOUTPUT1("psevp",psevp(i,k))
+DIAGOUTPUT1("p*dt",psevp(i,k)*dtcld)
+DIAGOUTPUT1("q+p*dt",qrs(i,k,2)+psevp(i,k)*dtcld)
+ qrs(i,k,2) = max(qrs(i,k,2)+psevp(i,k)*dtcld,0.)
+DIAGOUTPUT1("qrs_692",qrs(i,k,2))
+ xlf = xls-xl(i,k)
+ xlwork2 = -xl(i,k)*(prevp(i,k)+psevp(i,k))
+ t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
+ endif
+ enddo
+ enddo
+DIAGOUTPUT2("qs_800",qrs(IDEBUG,KDEBUG,2))
+!
+! Inline expansion for fpvs
+! qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
+! qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
+ hsub = xls
+ hvap = xlv0
+ cvap = cpv
+ ttp=t0c+0.01
+ dldt=cvap-cliq
+ xa=-dldt/rv
+ xb=xa+hvap/(rv*ttp)
+ dldti=cvap-cice
+ xai=-dldti/rv
+ xbi=xai+hsub/(rv*ttp)
+ do k = kts, kte
+ do i = its, ite
+ tr=ttp/t(i,k)
+ qs(i,k,1)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
+ qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
+ qs(i,k,1) = max(qs(i,k,1),qmin)
+ tr=ttp/t(i,k)
+ if(t(i,k).lt.ttp) then
+ qs(i,k,2)=psat*exp(log(tr)*(xai))*exp(xbi*(1.-tr))
+ else
+ qs(i,k,2)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
+ endif
+ qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
+ qs(i,k,2) = max(qs(i,k,2),qmin)
+ enddo
+ enddo
+!
+!----------------------------------------------------------------
+! pcond: condensational/evaporational rate of cloud water [HL A46] [RH83 A6]
+! if there exists additional water vapor condensated/if
+! evaporation of cloud water is not enough to remove subsaturation
+!
+ do k = kts, kte
+ do i = its, ite
+! work1(i,k,1) = conden(t(i,k),q(i,k),qs(i,k,1),xl(i,k),cpm(i,k))
+ work1(i,k,1) = ((max(q(i,k),qmin)-(qs(i,k,1)))/ &
+ (1.+(xl(i,k))*(xl(i,k))/(rv*(cpm(i,k)))*(qs(i,k,1))/((t(i,k))*(t(i,k)))))
+ work2(i,k) = qci(i,k,1)+work1(i,k,1)
+ pcond(i,k) = min(max(work1(i,k,1)/dtcld,0.),max(q(i,k),0.)/dtcld)
+ if(qci(i,k,1).gt.0..and.work1(i,k,1).lt.0.) &
+ pcond(i,k) = max(work1(i,k,1),-qci(i,k,1))/dtcld
+DIAGOUTPUT1("q_a",q(i,k))
+DIAGOUTPUT1("pcond_a",pcond(i,k))
+DIAGOUTPUT1("qs_1_a",qs(i,k,1))
+ q(i,k) = q(i,k)-pcond(i,k)*dtcld
+DIAGOUTPUT1("q_b",q(i,k))
+ qci(i,k,1) = max(qci(i,k,1)+pcond(i,k)*dtcld,0.)
+ t(i,k) = t(i,k)+pcond(i,k)*xl(i,k)/cpm(i,k)*dtcld
+ enddo
+ enddo
+!
+!
+!----------------------------------------------------------------
+! padding for small values
+!
+ do k = kts, kte
+ do i = its, ite
+ if(qci(i,k,1).le.qmin) qci(i,k,1) = 0.0
+ if(qci(i,k,2).le.qmin) qci(i,k,2) = 0.0
+ enddo
+ enddo
+#endif
+
+do i = its, ite
+do k = kts, kte
+kDIAGOUTPUT1("t_bot",t(i,k))
+enddo
+do k = kts, kte
+kDIAGOUTPUT1("q_bot",q(i,k))
+enddo
+do k = kts, kte
+kDIAGOUTPUT1("qc_bot",qci(i,k,1))
+enddo
+do k = kts, kte
+kDIAGOUTPUT1("qi_bot",qci(i,k,2))
+enddo
+do k = kts, kte
+kDIAGOUTPUT1("qr_bot",qrs(i,k,1))
+enddo
+do k = kts, kte
+kDIAGOUTPUT1("qs_bot",qrs(i,k,2))
+enddo
+do k = kts, kte
+kDIAGOUTPUT1("den_bot",den(i,k))
+enddo
+do k = kts, kte
+kDIAGOUTPUT1("p_bot",p(i,k))
+enddo
+do k = kts, kte
+kDIAGOUTPUT1("delz_bot",delz(i,k))
+enddo
+do k = kts, kte
+kDIAGOUTPUT1("cpm_bot",cpm(i,k))
+enddo
+do k = kts, kte
+kDIAGOUTPUT1("xl_bot",xl(i,k))
+enddo
+enddo
+ enddo ! big loops
+ END SUBROUTINE wsm52d
+! ...................................................................
+ REAL FUNCTION rgmma(x)
+!-------------------------------------------------------------------
+ IMPLICIT NONE
+!-------------------------------------------------------------------
+! rgmma function: use infinite product form
+ REAL :: euler
+ PARAMETER (euler=0.577215664901532)
+ REAL :: x, y
+ INTEGER :: i
+ if(x.eq.1.)then
+ rgmma=0.
+ else
+ rgmma=x*exp(euler*x)
+ do i=1,10000
+ y=float(i)
+ rgmma=rgmma*(1.000+x/y)*exp(-x/y)
+ enddo
+ rgmma=1./rgmma
+ endif
+ END FUNCTION rgmma
+!
+!--------------------------------------------------------------------------
+ REAL FUNCTION fpvs(t,ice,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c)
+!--------------------------------------------------------------------------
+ IMPLICIT NONE
+!--------------------------------------------------------------------------
+ REAL t,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c,dldt,xa,xb,dldti, &
+ xai,xbi,ttp,tr
+ INTEGER ice
+! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+ ttp=t0c+0.01
+ dldt=cvap-cliq
+ xa=-dldt/rv
+ xb=xa+hvap/(rv*ttp)
+ dldti=cvap-cice
+ xai=-dldti/rv
+ xbi=xai+hsub/(rv*ttp)
+ tr=ttp/t
+ if(t.lt.ttp.and.ice.eq.1) then
+ fpvs=psat*exp(log(tr)*(xai))*exp(xbi*(1.-tr))
+ else
+ fpvs=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
+ endif
+! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+ END FUNCTION fpvs
+!-------------------------------------------------------------------
+ SUBROUTINE wsm5init(den0,denr,dens,cl,cpv,allowed_to_read)
+!-------------------------------------------------------------------
+ IMPLICIT NONE
+!-------------------------------------------------------------------
+!.... constants which may not be tunable
+ REAL TENTIN :: den0,denr,dens,cl,cpv
+ LOGICAL TENTIN :: allowed_to_read
+ REAL :: pi
+!
+ pi = 4.*atan(1.)
+ xlv1 = cl-cpv
+!
+ qc0 = 4./3.*pi*denr*r0**3*xncr/den0 ! 0.419e-3 -- .61e-3
+ qck1 = .104*9.8*peaut/(xncr*denr)**(1./3.)/xmyu*den0**(4./3.) ! 7.03
+!
+ bvtr1 = 1.+bvtr
+ bvtr2 = 2.5+.5*bvtr
+ bvtr3 = 3.+bvtr
+ bvtr4 = 4.+bvtr
+ g1pbr = rgmma(bvtr1) ! 0.9312320
+ g3pbr = rgmma(bvtr3) ! 4.690781
+ g4pbr = rgmma(bvtr4) ! 17.81741
+ g5pbro2 = rgmma(bvtr2) ! 1.8273
+ pvtr = avtr*g4pbr/6.
+ eacrr = 1.0
+ pacrr = pi*n0r*avtr*g3pbr*.25*eacrr
+ precr1 = 2.*pi*n0r*.78
+ precr2 = 2.*pi*n0r*.31*avtr**.5*g5pbro2
+ xm0 = (di0/dicon)**2
+ xmmax = (dimax/dicon)**2
+ roqimax = 2.08e22*dimax**8
+!
+ bvts1 = 1.+bvts
+ bvts2 = 2.5+.5*bvts
+ bvts3 = 3.+bvts
+ bvts4 = 4.+bvts
+ g1pbs = rgmma(bvts1) !.8875
+ g3pbs = rgmma(bvts3) ! 3.011540
+ g4pbs = rgmma(bvts4) ! 10.26537
+ g5pbso2 = rgmma(bvts2) ! 1.550308
+ pvts = avts*g4pbs/6.
+ pacrs = pi*n0s*avts*g3pbs*.25
+ precs1 = 4.*n0s*.65
+ precs2 = 4.*n0s*.44*avts**.5*g5pbso2
+ pidn0r = pi*denr*n0r
+ pidn0s = pi*dens*n0s
+ pacrc = pi*n0s*avts*g3pbs*.25*eacrc
+!
+ rslopermax = 1./lamdarmax
+ rslopesmax = 1./lamdasmax
+ rsloperbmax = rslopermax ** bvtr
+ rslopesbmax = rslopesmax ** bvts
+ rsloper2max = rslopermax * rslopermax
+ rslopes2max = rslopesmax * rslopesmax
+ rsloper3max = rsloper2max * rslopermax
+ rslopes3max = rslopes2max * rslopesmax
+!
+ END SUBROUTINE wsm5init
+END MODULE module_mp_wsm5
+
+#ifdef STANDALONE
+
+ PROGRAM wsm_driver
+ USE module_mp_wsm5
+ IMPLICIT NONE
+ INTEGER thisstep
+ CHARACTER*80 pathtofile
+ COMMON /wsm5_driver_block/ thisstep, pathtofile
+ CHARACTER*80 fname
+ INTEGER ids,ide,jds,jde,kds,kde,ims,ime,jms,jme,kms,kme,its,ite,jts,jte,kts,kte
+ REAL , PARAMETER :: r_d = 287.
+ REAL , PARAMETER :: cp = 7.*r_d/2.
+ REAL , PARAMETER :: cv = cp-r_d
+ REAL , PARAMETER :: cliq = 4190.
+ REAL , PARAMETER :: rhowater = 1000.
+ REAL , PARAMETER :: rhosnow = 100.
+ REAL , PARAMETER :: rhoair0 = 1.28
+
+ print *,'Step number? Path to File? '
+ read(*,'(I3.2,A)') thisstep, pathtofile
+ write(fname,'(A,"wsm5_in_",i3.3)')trim(pathtofile),thisstep
+ print*,'opening ',trim(fname)
+write(0,*)__LINE__
+ open(45,file=fname,form='UNFORMATTED')
+write(0,*)__LINE__
+ read(45) ids,ide,jds,jde,kds,kde,ims,ime,jms,jme,kms,kme,its,ite,jts,jte,kts,kte
+write(0,*)__LINE__
+ close(45)
+write(0,*)__LINE__
+ thisstep = thisstep - 1
+write(0,*)'thisstep ',thisstep
+write(0,*)ids,ide,jds,jde,kds,kde,ims,ime,jms,jme,kms,kme,its,ite,jts,jte,kts,kte
+ CALL wsm5init(rhoair0,rhowater,rhosnow,cliq,cv,.FALSE.)
+ CALL wsm5( ids,ide,jds,jde,kds,kde,ims,ime,jms,jme,kms,kme,its,ite,jts,jte,kts,kte )
+
+ stop
+ END PROGRAM wsm_driver
+
+#endif
+
diff --git a/benchmarks/CUDA/WP/s2d.h b/benchmarks/CUDA/WP/s2d.h
new file mode 100644
index 0000000..f270d05
--- /dev/null
+++ b/benchmarks/CUDA/WP/s2d.h
@@ -0,0 +1,36 @@
+
+delt = delt_s
+g = g_s
+rd = rd_s
+rv = rv_s
+t0c = t0c_s
+den0 = den0_s
+cpd = cpd_s
+cpv = cpv_s
+ep1 = ep1_s
+ep2 = ep2_s
+qmin = qmin_s
+XLS = XLS_s
+XLV0 = XLV0_s
+XLF0 = XLF0_s
+cliq = cliq_s
+cice = cice_s
+psat = psat_s
+denr = denr_s
+
+th(its:ite,kts:kte,jts:jte) = th_s(its:ite,kts:kte,jts:jte)
+pii(its:ite,kts:kte,jts:jte) = pii_s(its:ite,kts:kte,jts:jte)
+q(its:ite,kts:kte,jts:jte) = q_s(its:ite,kts:kte,jts:jte)
+qc(its:ite,kts:kte,jts:jte) = qc_s(its:ite,kts:kte,jts:jte)
+qi(its:ite,kts:kte,jts:jte) = qi_s(its:ite,kts:kte,jts:jte)
+qr(its:ite,kts:kte,jts:jte) = qr_s(its:ite,kts:kte,jts:jte)
+qs(its:ite,kts:kte,jts:jte) = qs_s(its:ite,kts:kte,jts:jte)
+den(its:ite,kts:kte,jts:jte) = den_s(its:ite,kts:kte,jts:jte)
+p(its:ite,kts:kte,jts:jte) = p_s(its:ite,kts:kte,jts:jte)
+delz(its:ite,kts:kte,jts:jte) = delz_s(its:ite,kts:kte,jts:jte)
+rain(its:ite,jts:jte) = rain_s(its:ite,jts:jte)
+rainncv(its:ite,jts:jte) = rainncv_s(its:ite,jts:jte)
+sr(its:ite,jts:jte) = sr_s(its:ite,jts:jte)
+snow(its:ite,jts:jte) = snow_s(its:ite,jts:jte)
+snowncv(its:ite,jts:jte) = snowncv_s(its:ite,jts:jte)
+
diff --git a/benchmarks/CUDA/WP/spt.h b/benchmarks/CUDA/WP/spt.h
new file mode 100644
index 0000000..27a1467
--- /dev/null
+++ b/benchmarks/CUDA/WP/spt.h
@@ -0,0 +1,150 @@
+// macros and whatnot for translator
+
+#ifdef CUDA
+
+#if 0
+Types
+
+1. Array stored in device memory
+
+
+ 1.a fully dimensional
+ Example 3D
+ name: qc
+ how allocated: argument
+ dimensionality: ims:ime,kms:kme,jms:jme
+ index as: P3(ti,k,tj)
+
+ Example 2D
+ name: qc
+ how allocated: argument
+ dimensionality: ims:ime,jms:jme
+ index as: P2(ti,k,tj)
+
+
+ 1.b 1 dimensional (vertical only), local storage
+ Example
+ name: w3
+ how allocated: local
+ dimensionality: constant (MKX)
+ index as: w3[k]
+
+2. Array stored in shared memory
+
+ 2.a fully dimensional
+
+ 3D
+ Example
+ name: qc_s
+ how allocated: as offsets into sm[SM_SIZE]
+ dimensionality: bx*by*kx
+ index as: S3(ti,k,tj)
+
+ 2D
+ Example
+ name: <none>_s
+ how allocated: as offsets into sm[SM_SIZE]
+ dimensionality: bx*by
+ index as: S2(ti,tj)
+
+#endif
+
+#ifndef MKX
+ -- intentional syntax error -- need a defined constant MKX that is static number of levels --
+#endif
+
+#ifdef PROMOTE
+# define float double
+#endif
+
+#define SM_SIZE (0x1000-0xd4)
+#define MAX_THREADS_PER_BLOCK 512
+
+#define bi blockIdx.x
+#define bj blockIdx.y
+#define bx blockDim.x
+#define by blockDim.y
+#define ti threadIdx.x
+#define tj threadIdx.y
+
+# define ix (ime-ims+1)
+# define jx (jme-jms+1)
+# define kx (kme-kms+1)
+
+
+// basic indexing macros. indices are always given as global indices
+// in undecompsed Domain(ids:ide,jds:jde)
+//
+// That is, given IJ (global index), the global Index mapped to
+// a local index on a Patch(0:nx-1,0:ny-1) in Device Memory as:
+//
+// I - (ips-ims) + nx * ( J - (jps-jms) )
+//
+// where ips is the global index of the start of the patch (the -1 is
+// for translating from WRF fortran indices).
+//
+// The global index I is mapped to a local index on a GPU Block's
+// shared memory (0:bx-1, 0:by-1) as:
+//
+// I - (ips-ims) - bi * bx + by * ( J - (jps-jms) - bj * by )
+//
+// Where bi is the index into the GPU Block, and bx is the
+// GPU Block Width.
+
+// global to patch index converter
+#define GtoP(i,p,P) ((i)-(p)+(P))
+#define GtoB(i,n,N,p,P) ((i)-(p)+(P)-(n)*(N))
+
+// thread index to local memory index = i + bi * bx + ips - ims
+#define TtoP(i,a,b,c,d) ((i)+(a)*(b)+(c)-(d))
+
+//#define MAX(x,y) ((x)>(y)?(x):(y))
+//#define MIN(x,y) ((x)<(y)?(x):(y))
+#define MAX(x,y) max(x,y)
+#define MIN(x,y) min(x,y)
+
+// basic indexing macros
+//#define I2(i,j,m) ((i)+((j)*(m)))
+//#define I3(i,j,m,k,n) (I2(i,j,m)+((k)*(m)*(n)))
+#define I2(i,j,m) ((i)+(__mul24((j),(m))))
+#define I3(i,j,m,k,n) (I2(i,j,m)+(__mul24((k),__mul24((m),(n)))))
+
+#if (FLOAT_4 != 4)
+
+// index into a patch stored on device memory - 1
+# define P2(i,j) I2(TtoP(i,bi,bx,ips,ims),TtoP(j,bj,by,jps,jms),ime-ims+1)
+# define P3(i,k,j) I3(TtoP(i,bi,bx,ips,ims),k,ime-ims+1,TtoP(j,bj,by,jps,jms),kme-kms+1)
+// index into a block stored on shared memory
+# define S2(i,j) I2(i,j,bx)
+//# define S3(i,k,j) I3(i,k,bx,j,kme-kms+1)
+//# define S3(i,k,j) I3(i,j,bx,k,by)
+# define S3(i,k,j) I3(k,i,kx,j,bx)
+ // Local arrays in device mem
+# define LOCDM(a,s) float * a ; cudaMalloc( (void**) & a , (s)*sizeof(float)) ;
+ // Local scratch arrays in shared memory
+# define LOCSM(a,s) __shared__ float * a ; a = &(sm[isize]) ; isize += (s) ;
+#define ig (TtoP(ti,bi,bx,ips,ims))
+#define jg (TtoP(tj,bj,by,jps,jms))
+
+#else
+
+// index into a patch stored on device memory - 4
+# define P2(i,j) (I2(TtoP((i)*4,bi,(bx)*4,ips,ims),TtoP((j),bj,by,jps,jms),ime-ims+1)/4)
+# define P3(i,k,j) (I3(TtoP((i)*4,bi,(bx)*4,ips,ims),k,ime-ims+1,TtoP((j),bj,by,jps,jms),kme-kms+1)/4)
+// index into a block stored on shared memory
+# define S2(i,j) (I2(i*4,j,(bx)*4)/4)
+# define S3(i,k,j) (I3(i*4,k,(bx)*4,j,kme-kms+1)/4)
+ // Local arrays in device mem
+# define LOCDM(a,s) float4 * a ; cudaMalloc( (void**) & a , (s)*sizeof(float4)) ;
+ // Local scratch arrays in shared memory
+# define LOCSM(a,s) __shared__ Float4 * a ; a = &(sm[isize]) ; isize += (s)*4 ;
+#define ig (TtoP((ti)*4,bi,(bx)*4,ips,ims))
+#define jg (TtoP(tj,bj,by,jps,jms))
+
+#endif
+
+
+#endif
+
+
+
diff --git a/benchmarks/CUDA/WP/spt.pl b/benchmarks/CUDA/WP/spt.pl
new file mode 100755
index 0000000..7dc272a
--- /dev/null
+++ b/benchmarks/CUDA/WP/spt.pl
@@ -0,0 +1,219 @@
+#!/usr/bin/perl
+
+# 0 line number
+# 1 keyword
+# 2 action
+
+open DIRECTIVES, "| cat -n | grep '//_def_' | sed 's/ .*\\\/\\\/_def_/ / ' > directives__" or die "writing directives" ;
+open TEMP1, "> tmp1__" or die "making copy of code" ;
+
+while (<STDIN>) {
+ print DIRECTIVES ;
+ print TEMP1 ;
+}
+close DIRECTIVES ;
+close TEMP1 ;
+
+open DEBUG, "> debug__" or die ;
+## first pass, preprocess the directives
+
+open DIRECTIVES, "< directives__" or die ;
+while (<DIRECTIVES>) {
+ print TEMP1 ;
+ $line = $_ ;
+ @t = split( ' ',$line ) ;
+ $keyword = $t[1] ;
+ $action = $t[2] ;
+ @actionlist = split( ';', $action ) ;
+ foreach $act ( @actionlist ) {
+ @dim_vlist = split ( ':', $act ) ;
+ $dim = $dim_vlist[0] ; $vlist = $dim_vlist[1] ;
+ foreach $v ( split( ',', $vlist ) ) {
+ $vars{$v} = $v ;
+ $dimensionality{$v} = $dim ;
+ if ( $keyword eq "arg" ) {
+ $key{$v} = $keyword ;
+ if ( $dim eq "ikj" ) { $ikj_args{$v} = $v ; }
+ if ( $dim eq "ij" ) { $ij_args{$v} = $v ; }
+ }
+ if ( $keyword eq "local" ) {
+ $key{$v} = $keyword ;
+ if ( $dim eq "k" ) { $k_local{$v} = $v ; }
+ }
+ if ( $keyword eq "register" ) {
+ $key{$v} = $keyword ;
+ if ( $dim eq "0" ) { $register{$v} = $v ; }
+ }
+ if ( $keyword eq "copy_up_mem" ) {
+ if ( $key{$v} ne "arg" ) {
+ print "//warning: copy_up_mem of $v when $v is not arg.\n" ;
+ } else {
+ $copy_up_mem{$v} = $dim ;
+ if ( $dim eq "ikj" ) { $ikj_shared{$v} = $v ; }
+ if ( $dim eq "ij" ) { $ij_shared{$v} = $v ; }
+ }
+ }
+ if ( $keyword eq "shared_mem_local" ) {
+ $key{$v} = $keyword ;
+ $shared_mem_local{$v} = $dim ;
+ if ( $dim eq "ikj" ) { $ikj_shared{$v} = $v ; }
+ if ( $dim eq "ij" ) { $ij_shared{$v} = $v ; }
+ }
+ if ( $keyword eq "copy_down_mem" ) {
+ if ( $key{$v} ne "arg" ) {
+ print "//warning: copy_down_mem of $v when $v is not arg.\n" ;
+ } else {
+ $copy_down_mem{$v} = $dim ;
+ if ( $dim eq "ikj" ) { $ikj_shared{$v} = $v ; }
+ if ( $dim eq "ij" ) { $ij_shared{$v} = $v ; }
+ }
+ }
+ }
+ }
+}
+close DIRECTIVES ;
+
+## seond pass, modify the code
+## and preprocess deferred directives
+
+$spton=0 ;
+
+open TEMP1, "< tmp1__" or die ;
+while (<TEMP1>) {
+ $line = $_ ;
+ # toggle on and off between SPTSTART and SPTSTOP
+ if ( $line =~ "SPTSTART" ) { $spton = 1 ; }
+ elsif ( $line =~ "SPTSTOP" ) { $spton = 0 ; }
+ if ( $spton == 1 ) {
+
+ # handle copy_up_mem and copy_down_mem directives in line
+ if ( $line =~ m/\/\/\s*_def_\s+copy_up_mem\s/ ) {
+ @t = split( ' ',$line ) ;
+ $action = $t[2] ;
+ @dim_vlist = split ( ':', $action ) ;
+ $vlist = $dim_vlist[1] ;
+ foreach $v ( split( ',', $vlist ) ) {
+ print "LOCSM(${v}_s,bx*by*kx) ;\n" ;
+ }
+ print "{ int k ; \n" ;
+ foreach $v ( split( ',', $vlist ) ) {
+ print "for(k=kps-1;k<kpe;k++){${v}_s[S3(ti,k,tj)]=${v}[P3(ti,k,tj)];}\n" ;
+ }
+ print "}\n" ;
+ }
+ elsif ( $line =~ m/\/\/\s*_def_\s+register\s/ ) {
+ @t = split( ' ',$line ) ;
+ $action = $t[2] ;
+ @dim_vlist = split ( ':', $action ) ;
+ $vlist = $dim_vlist[1] ;
+ foreach $v ( split( ',', $vlist ) ) {
+ print "float ${v}_reg ;\n" ;
+ }
+ }
+ elsif ( $line =~ m/\/\/\s*_def_\s+shared_mem_local\s/ ) {
+ @t = split( ' ',$line ) ;
+ $action = $t[2] ;
+ @dim_vlist = split ( ':', $action ) ;
+ $vlist = $dim_vlist[1] ;
+ foreach $v ( split( ',', $vlist ) ) {
+ print "LOCSM(${v}_s,bx*by*kx) ;\n" ;
+ }
+ }
+ elsif ( $line =~ m/\/\/\s*_def_\s+copy_down_mem\s/ ) {
+ @t = split( ' ',$line ) ;
+ $action = $t[2] ;
+ @dim_vlist = split ( ':', $action ) ;
+ $vlist = $dim_vlist[1] ;
+ print "{ int k ; \n" ;
+ foreach $v ( split( ',', $vlist ) ) {
+ print "for(k=kps-1;k<kpe;k++){${v}[P3(ti,k,tj)]=${v}_s[S3(ti,k,tj)];}\n" ;
+ }
+ print "}\n" ;
+ }
+ elsif ( $line =~ m/\/\/\s*_def_\s+local\s/ ) {
+ $line = $_ ;
+ @t = split( ' ',$line ) ;
+ $keyword = $t[1] ;
+ $action = $t[2] ;
+ @actionlist = split( ';', $action ) ;
+ foreach $act ( @actionlist ) {
+ @dim_vlist = split ( ':', $act ) ;
+ $dim = $dim_vlist[0] ; $vlist = $dim_vlist[1] ;
+ foreach $v ( split( ',', $vlist ) ) {
+ if ( $dim eq "k" ) {
+ print "#if (FLOAT_4 == 4)\n" ;
+ print " Float4 ${v}[MKX] ; \n" ;
+ print "#else\n" ;
+ print " float ${v}[MKX] ; \n" ;
+ print "#endif\n" ;
+ }
+ }
+ }
+ }
+ # otherwise do not touch lines with // in them
+ elsif ( ! ($line =~ m/\/\//) ) {
+ @t = split( /\W+/,$line ) ;
+ %seen = "" ;
+ foreach $token ( @t ) {
+ if ( ! $seen{$token} ) {
+ $seen{$token} = $token ;
+ foreach $v ( keys %vars ) {
+ if ( "$v" eq "$token" ) {
+ $dim = $dimensionality{$v} ;
+ $keyw = $key{$v} ;
+ $nodex=0 ;
+ if ( $keyw eq "arg" || $keyw eq "shared_mem_local" ) {
+ if ( $copy_up_mem{$v} || $shared_mem_local{$v} ) {
+ if ( $dim eq "ikj" ) { $orig = $v."\\[\(.*\?\)\\]" ; $repl1 = $v."_s"."A|S3(ti," ; $repl2 = ",tj)B|" ; }
+ elsif ( $dim eq "ij" ) { $orig = $v ; $repl1 = $v."_s"."A|S2(ti,tj)B|" ; $repl2 = "" ; }
+ } else {
+ if ( $dim eq "ikj" ) { $orig = $v."\\[\(.*\?\)\\]" ; $repl1 = $v."A|P3(ti," ; $repl2 = ",tj)B|" ; }
+ elsif ( $dim eq "ij" ) { $orig = $v ; $repl1 = $v."A|P2(ti,tj)B|" ; $repl2 = "" ; }
+ }
+ } elsif ( $keyw eq "register" ) {
+ if ( $dim eq "0" ) { $orig = $v."\\[\(.*\?\)\\]" ; $repl1 = $v."_reg" ; $repl2 = "" ; $nodex=1}
+ } elsif ( $keyw eq "local" ) {
+ if ( $dim eq "k" ) { $orig = $v."\\[\(.*\?\)\\]" ; $repl1 = $v."A|" ; $repl2 = "B|" ; }
+ }
+# these repetitions are to handle multiple instances of the
+# variable being indexed differently on the same line.
+ if ( $nodex == 0 ) {
+ $line =~ s/(\W)$orig(\W)/$1$repl1$2$repl2$3/g ;
+ $line =~ s/(\W)$orig(\W)/$1$repl1$2$repl2$3/g ;
+ $line =~ s/(\W)$orig(\W)/$1$repl1$2$repl2$3/g ;
+ $line =~ s/(\W)$orig(\W)/$1$repl1$2$repl2$3/g ;
+ $line =~ s/(\W)$orig(\W)/$1$repl1$2$repl2$3/g ;
+ $line =~ s/(\W)$orig(\W)/$1$repl1$2$repl2$3/g ;
+ $line =~ s/(\W)$orig(\W)/$1$repl1$2$repl2$3/g ;
+ $line =~ s/(\W)$orig(\W)/$1$repl1$2$repl2$3/g ;
+ $line =~ s/(\W)$orig(\W)/$1$repl1$2$repl2$3/g ;
+ } else {
+ $line =~ s/(\W)$orig(\W)/$1$repl1$repl2$3/g ;
+ $line =~ s/(\W)$orig(\W)/$1$repl1$repl2$3/g ;
+ $line =~ s/(\W)$orig(\W)/$1$repl1$repl2$3/g ;
+ $line =~ s/(\W)$orig(\W)/$1$repl1$repl2$3/g ;
+ $line =~ s/(\W)$orig(\W)/$1$repl1$repl2$3/g ;
+ $line =~ s/(\W)$orig(\W)/$1$repl1$repl2$3/g ;
+ $line =~ s/(\W)$orig(\W)/$1$repl1$repl2$3/g ;
+ $line =~ s/(\W)$orig(\W)/$1$repl1$repl2$3/g ;
+ }
+
+ $line =~ s/A\|/[/g ;
+ $line =~ s/B\|/]/g ;
+ }
+ }
+ }
+ }
+ }
+ }
+ print $line ;
+}
+
+close TEMP1 ;
+close DEBUG ;
+
+unlink "directives__" ;
+unlink "tmp1__" ;
+unlink "debug__" ;
+
+
diff --git a/benchmarks/CUDA/WP/timings_20071205 b/benchmarks/CUDA/WP/timings_20071205
new file mode 100644
index 0000000..d0f8a5c
--- /dev/null
+++ b/benchmarks/CUDA/WP/timings_20071205
@@ -0,0 +1,24 @@
+vanilla: Call to WSM5 on host for step 10 is 1266066 microseconds
+chocolate:Call to wsm5_gpu (not including data xfer): 117187 microseconds
+chocolate:Call to wsm5_gpu (including data xfer): 128778 microseconds
+vanilla: Call to WSM5 on host for step 10 is 1265687 microseconds
+chocolate:Call to wsm5_gpu (not including data xfer): 117192 microseconds
+chocolate:Call to wsm5_gpu (including data xfer): 128576 microseconds
+chocolate: Call to WSM5 on host for step 10 is 146600 microseconds
+vanilla: Call to WSM5 on host for step 10 is 1265479 microseconds
+chocolate:Call to wsm5_gpu (not including data xfer): 117198 microseconds
+chocolate:Call to wsm5_gpu (including data xfer): 128639 microseconds
+chocolate: Call to WSM5 on host for step 10 is 146507 microseconds
+vanilla: Call to WSM5 on host for step 10 is 1261755 microseconds
+chocolate:Call to wsm5_gpu (not including data xfer): 117124 microseconds
+chocolate:Call to wsm5_gpu (including data xfer): 128550 microseconds
+chocolate: Call to WSM5 on host for step 10 is 146404 microseconds
+vanilla: Call to WSM5 on host for step 10 is 1266827 microseconds
+chocolate:Call to wsm5_gpu (not including data xfer): 117189 microseconds
+chocolate:Call to wsm5_gpu (including data xfer): 128725 microseconds
+chocolate: Call to WSM5 on host for step 10 is 146660 microseconds
+vanilla: Call to WSM5 on host for step 10 is 1260294 microseconds
+chocolate:Call to wsm5_gpu (not including data xfer): 124672 microseconds
+chocolate:Call to wsm5_gpu (including data xfer): 146497 microseconds
+chocolate: Call to WSM5 on host for step 10 is 164375 microseconds
+
diff --git a/benchmarks/CUDA/WP/util.h b/benchmarks/CUDA/WP/util.h
new file mode 100644
index 0000000..56cb493
--- /dev/null
+++ b/benchmarks/CUDA/WP/util.h
@@ -0,0 +1,265 @@
+ class Float4 {
+ public:
+ float x, y, z, w;
+
+ __device__ const Float4
+ operator+(const Float4& iv) const {
+ Float4 rv ;
+ rv.x = x + iv.x ;
+ rv.y = y + iv.y ;
+ rv.z = z + iv.z ;
+ rv.w = w + iv.w ;
+ return Float4( rv ) ;
+ }
+ __device__ const Float4
+ operator*(const Float4& iv) const {
+ Float4 rv ;
+ rv.x = x * iv.x ;
+ rv.y = y * iv.y ;
+ rv.z = z * iv.z ;
+ rv.w = w * iv.w ;
+ return Float4( rv ) ;
+ }
+ __device__ const Float4
+ operator/(const Float4& iv) const {
+ Float4 rv ;
+ rv.x = x / iv.x ;
+ rv.y = y / iv.y ;
+ rv.z = z / iv.z ;
+ rv.w = w / iv.w ;
+ return Float4( rv ) ;
+ }
+ __device__ const Float4
+ operator-(const Float4& iv) const {
+ Float4 rv ;
+ rv.x = x - iv.x ;
+ rv.y = y - iv.y ;
+ rv.z = z - iv.z ;
+ rv.w = w - iv.w ;
+ return Float4( rv ) ;
+ }
+
+ __device__ const Float4
+ operator+(const float iv) const {
+ Float4 rv ;
+ rv.x = x + iv ;
+ rv.y = y + iv ;
+ rv.z = z + iv ;
+ rv.w = w + iv ;
+ return Float4( rv ) ;
+ }
+ __device__ const Float4
+ operator*(const float iv) const {
+ Float4 rv ;
+ rv.x = x * iv ;
+ rv.y = y * iv ;
+ rv.z = z * iv ;
+ rv.w = w * iv ;
+ return Float4( rv ) ;
+ }
+ __device__ const Float4
+ operator/(const float iv) const {
+ Float4 rv ;
+ rv.x = x / iv ;
+ rv.y = y / iv ;
+ rv.z = z / iv ;
+ rv.w = w / iv ;
+ return Float4( rv ) ;
+ }
+ __device__ const Float4
+ operator-(const float iv) const {
+ Float4 rv ;
+ rv.x = x - iv ;
+ rv.y = y - iv ;
+ rv.z = z - iv ;
+ rv.w = w - iv ;
+ return Float4( rv ) ;
+ }
+ __device__ const Float4
+ operator-() const {
+ Float4 rv ;
+ rv.x = -x ;
+ rv.y = -y ;
+ rv.z = -z ;
+ rv.w = -w ;
+ return Float4( rv ) ;
+ }
+
+ __device__ void operator=(const float iv) {
+ x = iv ;
+ y = iv ;
+ z = iv ;
+ w = iv ;
+ }
+
+ __device__ void operator+=(const Float4 iv) {
+ x += iv.x ;
+ y += iv.y ;
+ z += iv.z ;
+ w += iv.w ;
+ }
+ __device__ void operator-=(const Float4 iv) {
+ x -= iv.x ;
+ y -= iv.y ;
+ z -= iv.z ;
+ w -= iv.w ;
+ }
+
+ };
+
+ __device__ const Float4
+ operator+( const float iv1, const Float4 iv2 ) {
+ Float4 rv ;
+ rv.x = iv1 + iv2.x ;
+ rv.y = iv1 + iv2.y ;
+ rv.z = iv1 + iv2.z ;
+ rv.w = iv1 + iv2.w ;
+ return Float4( rv ) ;
+ }
+ __device__ const Float4
+ operator*( const float iv1, const Float4 iv2 ) {
+ Float4 rv ;
+ rv.x = iv1 * iv2.x ;
+ rv.y = iv1 * iv2.y ;
+ rv.z = iv1 * iv2.z ;
+ rv.w = iv1 * iv2.w ;
+ return Float4( rv ) ;
+ }
+ __device__ const Float4
+ operator/( const float iv1, const Float4 iv2 ) {
+ Float4 rv ;
+ rv.x = iv1 / iv2.x ;
+ rv.y = iv1 / iv2.y ;
+ rv.z = iv1 / iv2.z ;
+ rv.w = iv1 / iv2.w ;
+ return Float4( rv ) ;
+ }
+ __device__ const Float4
+ operator-( const float iv1, const Float4 iv2 ) {
+ Float4 rv ;
+ rv.x = iv1 - iv2.x ;
+ rv.y = iv1 - iv2.y ;
+ rv.z = iv1 - iv2.z ;
+ rv.w = iv1 - iv2.w ;
+ return Float4( rv ) ;
+ }
+
+__device__ Float4 max ( const Float4 a , const Float4 b )
+{
+ Float4 c ;
+ c.x = (a.x>b.x)?a.x:b.x;
+ c.y = (a.y>b.y)?a.y:b.y;
+ c.z = (a.z>b.z)?a.z:b.z;
+ c.w = (a.w>b.w)?a.w:b.w;
+ return(c) ;
+}
+__device__ Float4 max ( const float a , const Float4 b )
+{
+ Float4 c ;
+ c.x = (a>b.x)?a:b.x;
+ c.y = (a>b.y)?a:b.y;
+ c.z = (a>b.z)?a:b.z;
+ c.w = (a>b.w)?a:b.w;
+ return(c) ;
+}
+__device__ Float4 max ( const Float4 a , const float b )
+{
+ Float4 c ;
+ c.x = (a.x>b)?a.x:b;
+ c.y = (a.y>b)?a.y:b;
+ c.z = (a.z>b)?a.z:b;
+ c.w = (a.w>b)?a.w:b;
+ return(c) ;
+}
+//__device__ float max ( const float a , const float b )
+//{
+// return(a>b)?a:b) ;
+//}
+
+__device__ Float4 min ( const Float4 a , const Float4 b )
+{
+ Float4 c ;
+ c.x = (a.x<b.x)?a.x:b.x;
+ c.y = (a.y<b.y)?a.y:b.y;
+ c.z = (a.z<b.z)?a.z:b.z;
+ c.w = (a.w<b.w)?a.w:b.w;
+ return(c) ;
+}
+__device__ Float4 min ( const float a , const Float4 b )
+{
+ Float4 c ;
+ c.x = (a<b.x)?a:b.x;
+ c.y = (a<b.y)?a:b.y;
+ c.z = (a<b.z)?a:b.z;
+ c.w = (a<b.w)?a:b.w;
+ return(c) ;
+}
+__device__ Float4 min ( const Float4 a , const float b )
+{
+ Float4 c ;
+ c.x = (a.x<b)?a.x:b;
+ c.y = (a.y<b)?a.y:b;
+ c.z = (a.z<b)?a.z:b;
+ c.w = (a.w<b)?a.w:b;
+ return(c) ;
+}
+
+__device__ Float4 trunc ( const Float4 a )
+{
+ Float4 c ;
+ c.x = trunc(a.x) ;
+ c.y = trunc(a.y) ;
+ c.z = trunc(a.z) ;
+ c.w = trunc(a.w) ;
+ return(c) ;
+}
+
+__device__ Float4 log ( const Float4 a )
+{
+ Float4 c ;
+ c.x = log(a.x) ; c.y = log(a.y) ; c.z = log(a.z) ; c.w = log(a.w) ;
+ return(c) ;
+}
+
+__device__ Float4 exp ( const Float4 a )
+{
+ Float4 c ;
+ c.x = exp(a.x) ; c.y = exp(a.y) ; c.z = exp(a.z) ; c.w = exp(a.w) ;
+ return(c) ;
+}
+
+__device__ Float4 sqrt ( const Float4 a )
+{
+ Float4 c ;
+ c.x = sqrt(a.x) ; c.y = sqrt(a.y) ; c.z = sqrt(a.z) ; c.w = sqrt(a.w) ;
+ return(c) ;
+}
+
+#if 0
+ int main() {
+
+ Float4 a, b, c ;
+
+ a.x = 0. ; a.y = 1. ; a.z = 2. ; a.w = 3. ;
+ b.x = 0. ; b.y = 1. ; b.z = 2. ; b.w = 3. ;
+
+ c = 2. + a ;
+ fprintf(stderr,"%f %f %f %f\n",a.x,a.y,a.z,a.w) ;
+ fprintf(stderr,"%f %f %f %f\n",b.x,b.y,b.z,b.w) ;
+ fprintf(stderr,"%f %f %f %f\n",c.x,c.y,c.z,c.w) ;
+ c = 2. * b ;
+ fprintf(stderr,"%f %f %f %f\n",a.x,a.y,a.z,a.w) ;
+ fprintf(stderr,"%f %f %f %f\n",b.x,b.y,b.z,b.w) ;
+ fprintf(stderr,"%f %f %f %f\n",c.x,c.y,c.z,c.w) ;
+ c = 2. - b ;
+ fprintf(stderr,"%f %f %f %f\n",a.x,a.y,a.z,a.w) ;
+ fprintf(stderr,"%f %f %f %f\n",b.x,b.y,b.z,b.w) ;
+ fprintf(stderr,"%f %f %f %f\n",c.x,c.y,c.z,c.w) ;
+ c = 2. / b ;
+ fprintf(stderr,"%f %f %f %f\n",a.x,a.y,a.z,a.w) ;
+ fprintf(stderr,"%f %f %f %f\n",b.x,b.y,b.z,b.w) ;
+ fprintf(stderr,"%f %f %f %f\n",c.x,c.y,c.z,c.w) ;
+
+ }
+#endif
+
diff --git a/benchmarks/CUDA/WP/util4.cu b/benchmarks/CUDA/WP/util4.cu
new file mode 100644
index 0000000..4b629d4
--- /dev/null
+++ b/benchmarks/CUDA/WP/util4.cu
@@ -0,0 +1,46 @@
+#ifndef PREPASS
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include "cublas.h"
+#endif
+
+__device__ float4 max4 ( const float4 a , const float4 b )
+{
+ float4 c ;
+ c.x = (a.x>b.x)?a.x:b.x;
+ c.y = (a.y>b.y)?a.y:b.y;
+ c.z = (a.z>b.z)?a.z:b.z;
+ c.w = (a.w>b.w)?a.w:b.w;
+ return(c) ;
+}
+__device__ float4 min4 ( const float4 a , const float4 b )
+{
+ float4 c ;
+ c.x = (a.x<b.x)?a.x:b.x;
+ c.y = (a.y<b.y)?a.y:b.y;
+ c.z = (a.z<b.z)?a.z:b.z;
+ c.w = (a.w<b.w)?a.w:b.w;
+ return(c) ;
+}
+
+__device__ float4 log4 ( const float4 a )
+{
+ float4 c ;
+ c.x = log(a.x) ; c.y = log(a.y) ; c.z = log(a.z) ; c.w = log(a.w) ;
+ return(c) ;
+}
+
+__device__ float4 exp4 ( const float4 a )
+{
+ float4 c ;
+ c.x = exp(a.x) ; c.y = exp(a.y) ; c.z = exp(a.z) ; c.w = exp(a.w) ;
+ return(c) ;
+}
+
+__device__ float4 sqrt4 ( const float4 a )
+{
+ float4 c ;
+ c.x = sqrt(a.x) ; c.y = sqrt(a.y) ; c.z = sqrt(a.z) ; c.w = sqrt(a.w) ;
+ return(c) ;
+}
diff --git a/benchmarks/CUDA/WP/wsm5.cu b/benchmarks/CUDA/WP/wsm5.cu
new file mode 100644
index 0000000..894e937
--- /dev/null
+++ b/benchmarks/CUDA/WP/wsm5.cu
@@ -0,0 +1,500 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include "cublas.h"
+
+#define IDEBUG 12
+#define JDEBUG 0
+
+#ifndef CRAY
+# ifdef NOUNDERSCORE
+# define WSM5_HOST wsm5_host
+# define WSM5_HOST_2 wsm5_host_2
+# define WSM5_GPU_INIT wsm5_gpu_init
+# define GET_WSM5_GPU_LEVELS get_wsm5_gpu_levels
+# else
+# ifdef F2CSTYLE
+# define WSM5_HOST wsm5_host__
+# define WSM5_HOST_2 wsm5_host_2__
+# define WSM5_GPU_INIT wsm5_gpu_init__
+# define GET_WSM5_GPU_LEVELS get_wsm5_gpu_levels__
+# else
+# define WSM5_HOST wsm5_host_
+# define WSM5_HOST_2 wsm5_host_2_
+# define WSM5_GPU_INIT wsm5_gpu_init_
+# define GET_WSM5_GPU_LEVELS get_wsm5_gpu_levels_
+# endif
+# endif
+#endif
+
+#define I2(i,j,m) ((i)+((j)*(m)))
+#define I3(i,j,m,k,n) (I2(i,j,m)+((k)*(m)*(n)))
+
+#if 1
+# define TODEV(A,s) float *A##_d;cudaMalloc((void**)&A##_d,((s))*sizeof(float));cudaMemcpy(A##_d,A,(s)*sizeof(float),cudaMemcpyHostToDevice);
+# define FROMDEV(A,s) cudaMemcpy(A,A##_d,(s)*sizeof(float),cudaMemcpyDeviceToHost);
+# define CLNUP(A) cudaFree(A##_d)
+#else
+# define TODEV(A,s) s1=rsl_internal_microclock_() ; float *A##_d;cudaMalloc((void**)&A##_d,((s))*sizeof(float));cudaMemcpy(A##_d,A,(s)*sizeof(float),cudaMemcpyHostToDevice); e1=rsl_internal_microclock_() ; fprintf(stderr,"TODEV %d\n",e1-s1)
+# define FROMDEV(A,s) s1=rsl_internal_microclock_() ; cudaMemcpy(A,A##_d,(s)*sizeof(float),cudaMemcpyDeviceToHost); e1=rsl_internal_microclock_() ; fprintf(stderr,"FROMDEV %d\n",e1-s1)
+# define CLNUP(A) s1=rsl_internal_microclock_() ; cudaFree(A##_d) ; e1=rsl_internal_microclock_() ; fprintf(stderr,"Free %d\n",e1-s1)
+#endif
+
+#if FLOAT_4==4
+#define TODEV2(A)s1=rsl_internal_microclock_();\
+float*A##_d;\
+cudaMalloc((void**)&A##_d,(dipe*djpe*sizeof(float)));\
+for(j=*jps-1;j<=*jpe-1;j++){\
+ for(i=*ips-1;i<=*ipe-1;i++){\
+ bigbuf[I2(i-*ips+1,j-*jps+1,dipe)]=\
+ A[I2(i-*ims+1,j-*jms+1,(*ime-*ims+1))];\
+}}\
+cudaMemcpy(A##_d,bigbuf,(dipe*djpe)*sizeof(float),cudaMemcpyHostToDevice);\
+e1=rsl_internal_microclock_();fprintf(stderr,"TODEV2 %d\n",e1-s1);
+
+#define TODEV3(A)s1=rsl_internal_microclock_();\
+float*A##_d;\
+cudaMalloc((void**)&A##_d,(dipe*djpe*dkpe*sizeof(float)));\
+for(j=*jps-1;j<=*jpe-1;j++){\
+ for(k=*kps-1;k<=*kpe-1;k++){\
+ for(i=*ips-1;i<=*ipe-1;i++){\
+ bigbuf[I3(i-*ips+1,k-*kps+1,dipe,j-*jps+1,dkpe)]=\
+ A[I3(i-*ims+1,k-*kms+1,*ime-*ims+1,j-*jms+1,*kme-*kms+1)];\
+}}}\
+cudaMemcpy(A##_d,bigbuf,(dipe*djpe*dkpe)*sizeof(float),cudaMemcpyHostToDevice);\
+e1=rsl_internal_microclock_();fprintf(stderr,"TODEV3 %d\n",e1-s1);
+
+// for debugging only
+#define TODEV3a(A)s1=rsl_internal_microclock_();\
+float*A##_d;\
+cudaMalloc((void**)&A##_d,(dipe*djpe*dkpe*sizeof(float)));\
+for(j=*jps-1;j<=*jpe-1;j++){\
+ for(k=*kps-1;k<=*kpe-1;k++){\
+ for(i=*ips-1;i<=*ipe-1;i++){\
+ bigbuf[I3(i-*ips+1,k-*kps+1,dipe,j-*jps+1,dkpe)]=\
+ A[I3(i-*ims+1,k-*kms+1,*ime-*ims+1,j-*jms+1,*kme-*kms+1)];\
+if (i==*ips-1){\
+ fprintf(stderr,"There %d %d %d (%d)| %d %d %d (%d)| %f\n",\
+ i-*ips+1,k-*kps+1,j-*jps+1,I3(i-*ips+1,k-*kps+1,dipe,j-*jps+1,dkpe),\
+ i-*ims+1,k-*kms+1,j-*jms+1,I3(i-*ims+1,k-*kms+1,*ime-*ims+1,j-*jms+1,*kme-*kms+1),\
+ bigbuf[I3(i-*ips+1,k-*kps+1,dipe,j-*jps+1,dkpe)]);\
+ A[I3(i-*ims+1,k-*kms+1,*ime-*ims+1,j-*jms+1,*kme-*kms+1)]=199.;\
+}\
+}}}\
+cudaMemcpy(A##_d,bigbuf,(dipe*djpe*dkpe)*sizeof(float),cudaMemcpyHostToDevice);\
+e1=rsl_internal_microclock_();fprintf(stderr,"TODEV3 %d\n",e1-s1);
+
+
+#define FROMDEV2(A) s1=rsl_internal_microclock_();\
+cudaMemcpy(bigbuf,A##_d,dipe*djpe*sizeof(float),cudaMemcpyDeviceToHost);\
+for(j=*jps-1;j<=*jpe-1;j++){\
+ for(i=*ips-1;i<=*ipe-1;i++){\
+ A[I2(i-*ims+1,j-*jms+1,(*ime-*ims+1))]=\
+ bigbuf[I2(i-*ips+1,j-*jps+1,dipe)];\
+}}\
+e1=rsl_internal_microclock_() ; fprintf(stderr,"FROMDEV2 %d\n",e1-s1);
+
+#define FROMDEV3(A) s1=rsl_internal_microclock_();\
+cudaMemcpy(bigbuf,A##_d,dipe*djpe*dkpe*sizeof(float),cudaMemcpyDeviceToHost);\
+for(j=*jps-1;j<=*jpe-1;j++){\
+ for(k=*kps-1;k<=*kpe-1;k++){\
+ for(i=*ips-1;i<=*ipe-1;i++){\
+ A[I3(i-*ims+1,k-*kms+1,*ime-*ims+1,j-*jms+1,*kme-*kms+1)]=\
+ bigbuf[I3(i-*ips+1,k-*kps+1,dipe,j-*jps+1,dkpe)];\
+}}}\
+e1=rsl_internal_microclock_();fprintf(stderr,"FROMDEV3 %d\n",e1-s1);
+#else
+# define TODEV3(A) TODEV(A,d3)
+# define TODEV2(A) TODEV(A,d2)
+# define FROMDEV3(A) FROMDEV(A,d3)
+# define FROMDEV2(A) FROMDEV(A,d2)
+#endif
+
+extern "C" int rsl_internal_microclock_() ;
+
+extern __global__ void wsm5_gpu (
+ float *th, float *pii //_def_ arg ikj:th,pii
+ ,float *q //_def_ arg ikj:q
+ ,float *qc,float *qi,float *qr,float *qs //_def_ arg ikj:qc,qi,qr,qs
+ ,float *den, float *p, float *delz //_def_ arg ikj:den,p,delz
+#ifdef DEBUGGAL_ARRAY
+,float *debuggal //_def_ arg ikj:debuggal
+#endif
+ ,float *rain,float *rainncv //_def_ arg ij:rain,rainncv
+ ,float *sr //_def_ arg ij:sr
+ ,float *snow,float *snowncv //_def_ arg ij:snow,snowncv
+ ,float delt
+,float* retvals
+ ,int ids, int ide, int jds, int jde, int kds, int kde
+ ,int ims, int ime, int jms, int jme, int kms, int kme
+ ,int ips, int ipe, int jps, int jpe, int kps, int kpe
+ ) ;
+
+extern "C" {
+
+int gethostname(char *name, size_t len);
+void bzero(void *s, size_t n);
+char *strcpy(char *dest, const char *src);
+
+#define MAXDEVICES 4
+#define MAXNODES 16
+int
+WSM5_GPU_INIT ( int * myproc , int * nproc, int * mydevice )
+{
+ float x, *x_d ;
+ int s, e ;
+ int i, dc, m ;
+ cudaError_t cerr ;
+ char hostname[64] ;
+ struct cudaDeviceProp dp ;
+// manage devices if multiheaded
+ cudaGetDeviceCount( &dc ) ;
+ if ( dc > MAXDEVICES )
+ { fprintf(stderr, "warning: more than %d devices on node (%d)\n", MAXDEVICES, dc ) ; dc = MAXDEVICES ; }
+ fprintf(stderr,"Number of devices on this node: %d\n", dc) ;
+
+ // i = *myproc % dc ;
+
+ i = *mydevice ;
+ if ( dc > 0 )
+ {
+ if ( cerr = cudaSetDevice( i ) ) {
+ fprintf(stderr," non-zero cerr %d\n",cerr) ;
+ }
+ }
+ gethostname( hostname, 64 ) ;
+ fprintf(stderr,"Setting device %02d for task %03d on host %s\n",i,*myproc,hostname) ;
+
+ if ( cerr = cudaGetDeviceProperties( &dp, i ) ) {
+ fprintf(stderr,"Device %02d: cerr = %d\n", cerr) ;
+ } else {
+ fprintf(stderr,"Device %02d: name %s\n",i,dp.name) ;
+ fprintf(stderr,"Device %02d: mem %d\n",i,dp.totalGlobalMem) ;
+ fprintf(stderr,"Device %02d: smem %d\n",i,dp.sharedMemPerBlock) ;
+ fprintf(stderr,"Device %02d: nreg %d\n",i,dp.regsPerBlock) ;
+ fprintf(stderr,"Device %02d: warp %d\n",i,dp.warpSize) ;
+ fprintf(stderr,"Device %02d: pitch %d\n",i,dp.memPitch) ;
+ fprintf(stderr,"Device %02d: maxthrds %d\n",i,dp.maxThreadsPerBlock) ;
+ fprintf(stderr,"Device %02d: maxtdim %d %d %d\n",i,dp.maxThreadsDim[0]
+ ,dp.maxThreadsDim[1]
+ ,dp.maxThreadsDim[2]) ;
+ fprintf(stderr,"Device %02d: maxgdim %d %d %d\n",i,dp.maxGridSize[0]
+ ,dp.maxGridSize[1]
+ ,dp.maxGridSize[2]) ;
+ fprintf(stderr,"Device %02d: clock %d\n",i,dp.clockRate) ;
+ fprintf(stderr,"Device %02d: talign %d\n",i,dp.textureAlignment) ;
+ }
+
+// do a dummy init to get things going
+ s=rsl_internal_microclock_() ;
+ cudaMalloc((void **)&x_d,sizeof(float)) ;
+ cudaMemcpy(x_d,&x,sizeof(float),cudaMemcpyHostToDevice) ;
+ cudaFree(x_d) ;
+ e=rsl_internal_microclock_() ;
+ fprintf(stderr,"wsm5_init: %d\n",e-s) ;
+ return(0) ;
+}
+
+int
+WSM5_HOST (
+ float *th, float *pii
+ ,float *q
+ ,float *qc, float *qi, float *qr, float *qs
+ ,float *den, float *p, float *delz
+#ifdef DEBUGGAL_ARRAY
+,float *debuggal
+#endif
+ ,float *delt
+ ,float *rain,float *rainncv
+ ,float *sr
+ ,float *snow,float *snowncv
+ ,int *ids, int *ide, int *jds, int *jde, int *kds, int *kde
+ ,int *ims, int *ime, int *jms, int *jme, int *kms, int *kme
+ ,int *ips, int *ipe, int *jps, int *jpe, int *kps, int *kpe
+ )
+{
+ int i, j, k ;
+ float *bigbuf ;
+ int s, e, s1, e1, s2, e2 ;
+ int d3 = (*ime-*ims+1) * (*jme-*jms+1) * (*kme-*kms+1) ;
+ int d2 = (*ime-*ims+1) * (*jme-*jms+1) ;
+
+
+//fprintf(stderr,"d3 = %d\n",d3) ;
+//fprintf(stderr,"d2 = %d\n",d2) ;
+
+#if FLOAT_4 == 4
+ int dips = 0 ; int dipe = (((*ipe-*ips+1+3)/4)*4) ; // round up four
+#else
+ int dips = 0 ; int dipe = (*ipe-*ips+1) ;
+#endif
+ int djps = 0 ; int djpe = (*jpe-*jps+1) ;
+ int dkps = 0 ; int dkpe = (*kpe-*kps+1) ;
+
+ bigbuf = (float *)malloc( dipe * djpe * dkpe * sizeof(float) ) ;
+
+//fprintf(stderr,"ids %d ide %d jds %d jde %d kds %d kde %d\n",*ids,*ide,*jds,*jde,*kds,*kde) ;
+//fprintf(stderr,"ims %d ime %d jms %d jme %d kms %d kme %d\n",*ims,*ime,*jms,*jme,*kms,*kme) ;
+//fprintf(stderr,"ips %d ipe %d jps %d jpe %d kps %d kpe %d\n",*ips,*ipe,*jps,*jpe,*kps,*kpe) ;
+//fprintf(stderr,"dipe %d djpe %d dkpe %d\n",dipe,djpe,dkpe) ;
+
+ s = rsl_internal_microclock_() ;
+ TODEV3(th) ;
+ TODEV3(pii) ;
+ TODEV3(q) ;
+ TODEV3(qc) ;
+ TODEV3(qi) ;
+ TODEV3(qr) ;
+ TODEV3(qs) ;
+ TODEV3(den) ;
+ TODEV3(p) ;
+ TODEV3(delz) ;
+#ifdef DEBUGGAL_ARRAY
+//TODEV3(debuggal) ;
+#endif
+ TODEV2(rain) ;
+ TODEV2(rainncv) ;
+ TODEV2(sr) ;
+ TODEV2(snow) ;
+ TODEV2(snowncv) ;
+float retvals[100] ;
+{ int k ;
+for (k=0 ;k<*kme-*kms+1;k++) {retvals[k] = 0.; }
+}
+TODEV(retvals,(*kme-*kms+1)) ;
+
+ int remx, remy ; // remainder?
+
+ remx = (*ipe-*ips+1) % XXX != 0 ? 1 : 0 ;
+ remy = (*jpe-*jps+1) % YYY != 0 ? 1 : 0 ;
+
+ dim3 dimBlock( XXX , YYY ) ;
+// fprintf(stderr,"ipe ips remx jpe jps remy %d %d %d %d %d %d\n",*ipe,*ips,remx,*jpe,*jps,remy) ;
+ dim3 dimGrid ( (*ipe-*ips+1) / XXX + remx , (*jpe-*jps+1) / YYY + remy ) ;
+
+ fprintf(stderr,"Call to wsm5_gpu: block dims %d %d\n",dimBlock.x,dimBlock.y) ;
+ fprintf(stderr,"Call to wsm5_gpu: grid dims %d %d\n",dimGrid.x,dimGrid.y) ;
+
+#if 1
+//fprintf(stderr,"calling wsm5_gpu \n") ;
+//fprintf(stderr,"d %d %d %d %d %d %d\n",dips+1 , (*ipe-*ips+1) , djps+1 , (*jpe-*jps+1) , dkps+1 , (*kpe-*kps+1)) ;
+//fprintf(stderr,"m %d %d %d %d %d %d\n",dips+1 , dipe , djps+1 , djpe , dkps+1 , dkpe ) ;
+//fprintf(stderr,"p %d %d %d %d %d %d\n",dips+1 , dipe , djps+1 , djpe , dkps+1 , dkpe ) ;
+
+ s2 = rsl_internal_microclock_() ;
+ wsm5_gpu <<< dimGrid, dimBlock >>> (
+ th_d, pii_d, q_d, qc_d, qi_d, qr_d, qs_d, den_d, p_d, delz_d
+#ifdef DEBUGGAL_ARRAY
+,debuggal_d
+#endif
+ ,rain_d,rainncv_d
+ ,sr_d
+ ,snow_d,snowncv_d
+ ,*delt
+,retvals_d
+ ,dips+1 , (*ipe-*ips+1) , djps+1 , (*jpe-*jps+1) , dkps+1 , (*kpe-*kps+1)
+ ,dips+1 , dipe , djps+1 , djpe , dkps+1 , dkpe
+ ,dips+1 , dipe , djps+1 , djpe , dkps+1 , dkpe
+ ) ;
+ cudaThreadSynchronize() ;
+ e2 = rsl_internal_microclock_() ;
+ fprintf(stderr,"Call to wsm5_gpu (not including data xfer): %d microseconds\n",e2-s2) ;
+#endif
+
+ FROMDEV3(th) ;
+ FROMDEV3(pii) ;
+ FROMDEV3(q) ;
+ FROMDEV3(qc) ;
+ FROMDEV3(qi) ;
+ FROMDEV3(qr) ;
+ FROMDEV3(qs) ;
+#ifdef DEBUGGAL_ARRAY
+FROMDEV3(debuggal) ;
+#endif
+ FROMDEV2(rain) ;
+ FROMDEV2(rainncv) ;
+ FROMDEV2(sr) ;
+ FROMDEV2(snow) ;
+ FROMDEV2(snowncv) ;
+ e = rsl_internal_microclock_() ;
+//fprintf(stderr,"retrieving retvals %d\n",*kme-*kms+1) ;
+FROMDEV(retvals,(*kme-*kms+1)) ;
+ fprintf(stderr,"Call to wsm5_gpu (including data xfer): %d microseconds\n",e-s) ;
+
+{ int k ;
+//for (k=0 ;k<*kme-*kms+1;k++) {fprintf(stderr,"retvals %d %f\n",k,retvals[k]) ;}
+//for (k=0 ;k<5;k++) {fprintf(stderr,"retvals %d %f\n",k,retvals[k]) ;}
+}
+
+ CLNUP(th) ;
+ CLNUP(pii) ;
+ CLNUP(q) ;
+ CLNUP(qc) ;
+ CLNUP(qi) ;
+ CLNUP(qr) ;
+ CLNUP(qs) ;
+ CLNUP(den) ;
+ CLNUP(p) ;
+ CLNUP(delz) ;
+#ifdef DEBUGGAL_ARRAY
+CLNUP(debuggal) ;
+#endif
+ CLNUP(rain) ;
+ CLNUP(rainncv) ;
+ CLNUP(sr) ;
+ CLNUP(snow) ;
+ CLNUP(snowncv) ;
+CLNUP(retvals) ;
+
+ return(0) ;
+}
+
+#if 0
+static int first_wsm5_host_2=1 ;
+// 3d
+static float * th_h ;
+static float * pii_h ;
+static float * q_h ;
+static float * qc_h ;
+static float * qi_h ;
+static float * qr_h ;
+static float * qs_h ;
+static float * den_h ;
+static float * p_h ;
+static float * delz_h ;
+// 2d
+static float * rain_h ;
+static float * rainncv_h ;
+static float * sr_h ;
+static float * snow_h ;
+static float * snowncv_h ;
+
+// idea here is to copy the data into pinned (paged-locked) mem for faster xfer
+int
+WSM5_HOST_2 (
+ float *th, float *pii
+ ,float *q
+ ,float *qc, float *qi, float *qr, float *qs
+ ,float *den, float *p, float *delz
+ ,float *delt
+ ,float *rain,float *rainncv
+ ,float *sr
+ ,float *snow,float *snowncv
+ ,int *ids, int *ide, int *jds, int *jde, int *kds, int *kde
+ ,int *ims, int *ime, int *jms, int *jme, int *kms, int *kme
+ ,int *ips, int *ipe, int *jps, int *jpe, int *kps, int *kpe
+ )
+{
+ int i,j,k ;
+ float *ptr ;
+ int d3 = (*ipe-*ips+1) * (*jpe-*jps+1) * (*kpe-*kps+1) ;
+ int d2 = (*ipe-*ips+1) * (*jpe-*jps+1) ;
+
+ if ( first_wsm5_host_2 == 1 ) {
+ cudaMallocHost( (void **)&th_h , d3*sizeof(float) ) ; //3d
+ cudaMallocHost( (void **)&pii_h , d3*sizeof(float) ) ;
+ cudaMallocHost( (void **)&q_h , d3*sizeof(float) ) ;
+ cudaMallocHost( (void **)&qc_h , d3*sizeof(float) ) ;
+ cudaMallocHost( (void **)&qi_h , d3*sizeof(float) ) ;
+ cudaMallocHost( (void **)&qr_h , d3*sizeof(float) ) ;
+ cudaMallocHost( (void **)&qs_h , d3*sizeof(float) ) ;
+ cudaMallocHost( (void **)&den_h , d3*sizeof(float) ) ;
+ cudaMallocHost( (void **)&p_h , d3*sizeof(float) ) ;
+ cudaMallocHost( (void **)&delz_h , d3*sizeof(float) ) ;
+ cudaMallocHost( (void **)&rain_h , d2*sizeof(float) ) ; //2d
+ cudaMallocHost( (void **)&rainncv_h , d2*sizeof(float) ) ;
+ cudaMallocHost( (void **)&sr_h , d2*sizeof(float) ) ;
+ cudaMallocHost( (void **)&snow_h , d2*sizeof(float) ) ;
+ cudaMallocHost( (void **)&snowncv_h , d2*sizeof(float) ) ;
+ first_wsm5_host_2 = 0 ;
+ }
+
+#define PIN3(A) ptr=A##_h;for(j=*jps;j<=*jpe;j++){for(k=*kps;k<=*kpe;k++){for(i=*ips;i<=*ipe;i++){*ptr++=A [I3(i-*ims,k-*kms,*ime-*ims+1,j-*jms,*kme-*kms+1)];}}};
+#define PIN2(A) ptr=A##_h;for(j=*jps;j<=*jpe;j++);for(i=*ips;i<=*ipe;i++){*ptr++=A [I2(i-*ims,j-*jms,*ime-*ims+1)];};
+#define UNPIN3(A) ptr=A##_h;for(j=*jps;j<=*jpe;j++);for(k=*kps;k<=*kpe;k++);for(i=*ips;i<=*ipe;i++){A [I3(i-*ims,k-*kms,*ime-*ims+1,j-*jms,*kme-*kms+1)]=*ptr++;};
+#define UNPIN2(A) ptr=A##_h;for(j=*jps;j<=*jpe;j++);for(i=*ips;i<=*ipe;i++){A [I2(i-*ims,j-*jms,*ime-*ims+1)]=*ptr++;};
+
+ PIN3(th) ;
+ PIN3(th) ;
+ PIN3(pii) ;
+ PIN3(q) ;
+ PIN3(qc) ;
+ PIN3(qi) ;
+ PIN3(qr) ;
+ PIN3(qs) ;
+ PIN3(den) ;
+ PIN3(p) ;
+ PIN3(delz) ;
+ PIN2(rain) ;
+ PIN2(rainncv) ;
+ PIN2(sr) ;
+ PIN2(snow) ;
+ PIN2(snowncv) ;
+
+ WSM5_HOST (
+ th_h, pii_h
+ ,q_h
+ ,qc_h, qi_h, qr_h, qs_h
+ ,den_h, p_h, delz_h
+ ,delt
+ ,rain_h,rainncv_h
+ ,sr_h
+ ,snow_h,snowncv_h
+ ,ids, ide, jds, jde, kds, kde
+ ,ips, ipe, jps, jpe, kps, kpe
+ ,ips, ipe, jps, jpe, kps, kpe
+ ) ;
+
+
+ UNPIN3(th) ;
+ UNPIN3(th) ;
+ UNPIN3(pii) ;
+ UNPIN3(q) ;
+ UNPIN3(qc) ;
+ UNPIN3(qi) ;
+ UNPIN3(qr) ;
+ UNPIN3(qs) ;
+ UNPIN3(den) ;
+ UNPIN3(p) ;
+ UNPIN3(delz) ;
+ UNPIN2(rain) ;
+ UNPIN2(rainncv) ;
+ UNPIN2(sr) ;
+ UNPIN2(snow) ;
+ UNPIN2(snowncv) ;
+
+}
+#endif
+
+int
+GET_WSM5_GPU_LEVELS ( int * retval )
+{
+ *retval = MKX ; /* MKX is hard coded value set in the makefile */
+}
+}
+
+#if 0
+main( int argc, char **argv )
+{
+ float *th ; float *pii ; float *q ;
+ float *qc; float *qi; float *qr; float *qs ;
+ float *den; float *p; float *delz ;
+ float *delt ;
+ float *rain;float *rainncv ;
+ float *sr ;
+ float *snow;float *snowncv ;
+ int *ids; int *ide; int *jds; int *jde; int *kds; int *kde ;
+ int *ims; int *ime; int *jms; int *jme; int *kms; int *kme ;
+ int *ips; int *ipe; int *jps; int *jpe; int *kps; int *kpe ;
+ WSM5_HOST (
+ th, pii, q, qc, qi, qr, qs, den, p, delz
+ ,rain,rainncv
+ ,sr
+ ,snow,snowncv
+ ,delt
+ ,ids, ide, jds, jde, kds, kde
+ ,ims, ime, jms, jme, kms, kme
+ ,ips, ipe, jps, jpe, kps, kpe
+ ) ;
+}
+#endif
diff --git a/benchmarks/CUDA/WP/wsm5_constants.h b/benchmarks/CUDA/WP/wsm5_constants.h
new file mode 100644
index 0000000..5c435dc
--- /dev/null
+++ b/benchmarks/CUDA/WP/wsm5_constants.h
@@ -0,0 +1,92 @@
+// WSM5 Constants
+
+#if 1
+# define epsilon 1.e-15
+# define r_d 287.
+# define rhoair0 1.28
+# define rhosnow 100.
+# define dens rhosnow
+# define rhowater 1000.
+# define svpt0 .27314999389648438e+03
+# define xlv 2.5e6
+#endif
+
+#define g 0.981000041961670E+01
+#define r_v 0.461600006103516E+03
+#define rv r_v
+#define cice 0.210600000000000E+04
+#define cliq 0.419000000000000E+04
+#define denr 0.100000000000000E+04
+#define den0 0.127999997138977E+01
+#define xlf0 0.350000000000000E+06
+#define xlv0 0.250000000000000E+07
+#define xls 0.285000000000000E+07
+#define t0c 0.273149993896484E+03
+#define qmin 0.100000000362749E-14
+#define ep1 0.608362436294556E+00
+#define ep2 0.621750414371490E+00
+#define psat 0.610780029296875E+03
+#define alpha 0.120000000000000E+00
+#define n0smax 0.100000000000000E+12
+#define n0s 0.200000000000000E+07
+#define n0r 0.800000000000000E+07
+#define qcrmin 0.100000000000000E-08
+#define avtr 0.841900000000000E+03
+#define bvtr 0.800000000000000E+00
+#define g1pbr 0.931232915622909E+00
+#define g3pbr 0.469078683336385E+01
+#define g4pbr 0.178173289058329E+02
+#define g5pbro2 0.182658695197891E+01
+#define avts 0.117200000000000E+02
+#define bvts 0.410000000000000E+00
+#define g1pbs 0.886676521690526E+00
+#define g3pbs 0.301156382231086E+01
+#define g4pbs 0.102654190601850E+02
+#define g5pbso2 1.550308
+#define r0 0.800000000000000E-05
+#define peaut 0.550000000000000E+00
+#define xncr 0.300000000000000E+09
+#define xmyu 0.171800000000000E-04
+#define lamdarmax 0.800000000000000E+05
+#define lamdasmax 0.100000000000000E+06
+#define lamdagmax 0.600000000000000E+05
+#define pi 0.314159265358979E+01
+#define dicon 0.119000000000000E+02
+#define dimax 0.500000000000000E-03
+#define pfrz1 0.100000000000000E+03
+#define pfrz2 0.660000000000000E+00
+#define eacrr 0.100000000000000E+01
+#define eacrc 0.100000000000000E+01
+
+ double cpv = 4.*r_v ;
+ double cp = 7.*r_d/2. ;
+ double cv = cp-r_d ;
+ double cpd = cp ;
+
+ //double ep_1 = r_v/r_d-1. ;
+ //double ep_2 = r_d/r_v ;
+ double pvtr = avtr*g4pbr/6. ;
+ double pvts = avts*g4pbs/6. ;
+ double xlv1 = cliq - cv ;
+
+ double rslopermax = 1./lamdarmax ;
+ double rslopesmax = .10000000000000001e-04 ; // 1./lamdasmax ;
+ double rsloperbmax = 0.11954406247375457E-03 ; // exp(log(rslopermax) * bvtr) ;
+ double rslopesbmax = .89125093813374589e-02 ; // exp(log(rslopesmax) * bvts) ;
+ double rsloper2max = rslopermax * rslopermax ;
+ double rslopes2max = rslopesmax * rslopesmax ;
+ double rsloper3max = rsloper2max * rslopermax ;
+ double rslopes3max = rslopes2max * rslopesmax ;
+
+ double pidn0r = pi*denr*n0r ;
+ double pidn0s = pi*dens*n0s ;
+
+ double precs1 = 4.*n0s*.65 ;
+ double precs2 = 4.*n0s*.44*sqrt(avts)*g5pbso2 ;
+ double qc0 = 4./3.*pi*denr*(r0*r0*r0)*xncr/den0 ;
+ double qck1 = .104*9.8*peaut/pow((xncr*denr),(1./3.))/xmyu*pow(den0,(4./3.)) ;
+ double precr1 = 2.*pi*n0r*.78 ;
+ double precr2 = 2.*pi*n0r*.31*sqrt(avtr)*g5pbro2 ;
+ double pacrr = pi*n0r*avtr*g3pbr*.25*eacrr ;
+ double pacrc = pi*n0s*avts*g3pbs*.25*eacrc ;
+ double roqimax = 2.08e22*pow(dimax,8) ;
diff --git a/benchmarks/CUDA/WP/wsm5_gpu.cu b/benchmarks/CUDA/WP/wsm5_gpu.cu
new file mode 100644
index 0000000..7000cb3
--- /dev/null
+++ b/benchmarks/CUDA/WP/wsm5_gpu.cu
@@ -0,0 +1,783 @@
+#define REWORK_FALL
+#define REWORK_PART2
+// wsm5_gpu.cu gets preprocessed by spt.pl, which handles the _def_ directives before it is compiled
+
+#ifndef PREPASS
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include "cublas.h"
+#endif
+
+#define IDEBUG ((DEBUG_I)-2)
+#define JDEBUG ((DEBUG_J)-2)
+#define KDEBUG (DEBUG_K)
+
+// this is an M4 include
+include(debug.m4)
+
+//SPTSTART
+
+#include "spt.h"
+
+#include "util.h"
+
+# define float float
+
+
+
+__global__ void wsm5_gpu (
+ float *th, float *pii //_def_ arg ikj:th,pii
+ ,float *q //_def_ arg ikj:q
+ ,float *qc,float *qi,float *qr,float *qs //_def_ arg ikj:qc,qi,qr,qs
+ ,float *den, float *p, float *delz //_def_ arg ikj:den,p,delz
+#ifdef DEBUGAL_ARRAY
+,float *debuggal //_def_ arg ikj:debuggal
+#endif
+ ,float *rain,float *rainncv //_def_ arg ij:rain,rainncv
+ ,float *sr //_def_ arg ij:sr
+ ,float *snow,float *snowncv //_def_ arg ij:snow,snowncv
+ ,float delt
+,float* retvals
+ ,int ids, int ide, int jds, int jde, int kds, int kde
+ ,int ims, int ime, int jms, int jme, int kms, int kme
+ ,int ips, int ipe, int jps, int jpe, int kps, int kpe
+ )
+{
+
+ float xlf, xmi, acrfac, vt2i, vt2s, supice, diameter ;
+ float roqi0, xni0, qimax, value, source, factor, xlwork2 ;
+ float t_k, q_k, qr_k, qc_k, qs_k, qi_k, qs1_k, qs2_k, cpm_k, xl_k, xni_k, w1_k, w2_k, w3_k ;
+
+#define hsub xls
+#define hvap xlv0
+#define cvap cpv
+ float ttp ;
+ float dldt ;
+ float xa ;
+ float xb ;
+ float dldti ;
+ float xai ;
+ float xbi ;
+
+ //_def_ local k:qs1,qs2,rh1,rh2
+
+#ifdef DEBUGAL_ARRAY
+ debuggal[0] = 999.00 ;
+#endif
+
+if ( ig < ide-ids+1 && jg < jde-jds+1 ) {
+
+
+ int k ;
+
+#include "wsm5_constants.h"
+
+ //_def_ local k:t
+ //_def_ local k:cpm,xl
+
+ for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
+ t[k] = th[k] * pii[k] ;
+ }
+
+ for( k=kps-1 ;k<=kpe-1;k++) {
+ if ( qc[k] < 0. ) { qc[k] = 0. ; }
+ if ( qi[k] < 0. ) { qi[k] = 0. ; }
+ if ( qr[k] < 0. ) { qr[k] = 0. ; }
+ if ( qs[k] < 0. ) { qs[k] = 0. ; }
+ }
+
+// 564 !----------------------------------------------------------------
+// 565 ! latent heat for phase changes and heat capacity. neglect the
+// 566 ! changes during microphysical process calculation
+// 567 ! emanuel(1994)
+
+#define CPMCAL(x) (cpd*(1.-max(x,qmin))+max(x,qmin)*cpv)
+#define XLCAL(x) (xlv0-xlv1*((x)-t0c))
+
+ for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
+ cpm[k] = CPMCAL(q[k]) ;
+ xl[k] = XLCAL(t[k]) ;
+ }
+
+// 576 !----------------------------------------------------------------
+// 577 ! compute the minor time steps.
+
+ float dtcldcr = 120. ;
+ int loops = delt/dtcldcr+.5 ;
+
+ loops = MAX(loops,1) ;
+ float dtcld = delt/loops ;
+ if ( delt <= dtcldcr) dtcld = delt ;
+
+ int loop ;
+
+
+ for ( loop = 1 ; loop <= loops ; loop++ ) {
+// 585 !----------------------------------------------------------------
+// 586 ! initialize the large scale variables
+ int mstep = 1 ;
+
+ ttp=t0c+0.01 ;
+ dldt=cvap-cliq ;
+ xa=-dldt/rv ;
+ xb=xa+hvap/(rv*ttp) ;
+ dldti=cvap-cice ;
+ xai=-dldti/rv ;
+ xbi=xai+hsub/(rv*ttp) ;
+
+
+ float tr, ltr, tt, pp, qq ;
+
+ for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
+
+ pp = p[k] ;
+ tt = t[k] ;
+ tr = ttp/tt ;
+ ltr = log(tr) ;
+
+ qq=psat*exp(ltr*(xa)+xb*(1.-tr)) ;
+ qq=ep2*qq/(pp-qq) ;
+ qs1[k] = MAX(qq,qmin) ;
+ rh1[k] = MAX( q[k]/qs1[k],qmin) ;
+
+ if( tt < ttp ) {
+ qq=psat*exp(ltr*(xai)+xbi*(1.-tr)) ;
+ } else {
+ qq=psat*exp(ltr*(xa)+xb*(1.-tr)) ;
+ }
+ qq = ep2 * qq / (pp - qq) ;
+ qs2[k] = MAX(qq,qmin) ;
+ rh2[k] = MAX(q[k]/qs2[k],qmin) ;
+
+ }
+
+ //_def_ register 0:prevp,psdep,praut,psaut,pracw,psaci,psacw,pigen,pidep,pcond,psmlt,psevp
+ //_def_ local k:xni
+
+ for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
+ xni[k] = 1.e3 ;
+ }
+
+// diffus(x,y) = 8.794e-5 * exp(log(x)*(1.81)) / y ! 8.794e-5*x**1.81/y
+// viscos(x,y) = 1.496e-6 * (x*sqrt(x)) /(x+120.)/y ! 1.496e-6*x**1.5/(x+120.)/y
+// xka(x,y) = 1.414e3*viscos(x,y)*y
+// diffac(a,b,c,d,e) = d*a*a/(xka(c,d)*rv*c*c)+1./(e*diffus(c,b))
+// venfac(a,b,c) = exp(log((viscos(b,c)/diffus(b,a)))*((.3333333))) &
+// /sqrt(viscos(b,c))*sqrt(sqrt(den0/c))
+
+#define DIFFUS(x,y) (8.794e-5 * exp(log(x)*(1.81)) / (y))
+#define VISCOS(x,y) (1.496e-6 * ((x)*sqrt(x)) /((x)+120.)/(y))
+#define XKA(x,y) (1.414e3*VISCOS((x),(y))*(y))
+#define DIFFAC(a,b,c,d,e) ((d)*(a)*(a)/(XKA((c),(d))*rv*(c)*(c))+1./((e)*DIFFUS((c),(b))))
+#define VENFAC(a,b,c) (exp(log((VISCOS((b),(c))/DIFFUS((b),(a))))*((.3333333)))*rsqrt(VISCOS((b),(c)))*sqrt(sqrt(den0/(c))))
+#define CONDEN(a,b,c,d,e) ((MAX((b),qmin)-(c))/(1.+(d)*(d)/(rv*(e))*(c)/((a)*(a))))
+
+#define LAMDAR(x,y) sqrt(sqrt(pidn0r/((x)*(y))))
+#define LAMDAS(x,y,z) sqrt(sqrt(pidn0s*(z)/((x)*(y))))
+
+// calculate mstep for this colum
+
+ //_def_ local k:rsloper,rslopebr,rslope2r,rslope3r
+ //_def_ local k:rslopes,rslopebs,rslope2s,rslope3s
+ //_def_ local k:denfac
+ //_def_ local k:n0sfac
+ //_def_ local k:w1,w2,w3
+
+
+ float w ;
+ float rmstep ;
+ int numdt ;
+ for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
+ float supcol = t0c - t[k] ;
+ n0sfac[k] = MAX(MIN(exp(alpha*supcol),n0smax/n0s),1.) ;
+ if ( qr[k] <= qcrmin ) {
+ rsloper[k] = rslopermax ;
+ rslopebr[k] = rsloperbmax ;
+ rslope2r[k] = rsloper2max ;
+ rslope3r[k] = rsloper3max ;
+ } else {
+ rsloper[k] = 1./LAMDAR(qr[k],den[k]) ;
+ rslopebr[k] = exp(log(rsloper[k])*bvtr) ;
+ rslope2r[k] = rsloper[k] * rsloper[k] ;
+ rslope3r[k] = rslope2r[k] * rsloper[k] ;
+ }
+ if ( qs[k] <= qcrmin ) {
+ rslopes[k] = rslopesmax ;
+ rslopebs[k] = rslopesbmax ;
+ rslope2s[k] = rslopes2max ;
+ rslope3s[k] = rslopes3max ;
+ } else {
+ rslopes[k] = 1./LAMDAS(qs[k],den[k],n0sfac[k]) ;
+ rslopebs[k] = exp(log(rslopes[k])*bvts) ;
+ rslope2s[k] = rslopes[k] * rslopes[k] ;
+ rslope3s[k] = rslope2s[k] * rslopes[k] ;
+ }
+ denfac[k] = sqrt(den0/den[k]) ;
+ w1[k] = pvtr*rslopebr[k]*denfac[k]/delz[k] ;
+ w2[k] = pvts*rslopebs[k]*denfac[k]/delz[k] ;
+
+ w = MAX(w1[k],w2[k]) ;
+ numdt = MAX((int)trunc(w*dtcld+.5+.5),1) ;
+ if ( numdt >= mstep ) mstep = numdt ;
+//-------------------------------------------------------------
+// Ni: ice crystal number concentration [HDC 5c]
+//-------------------------------------------------------------
+ float temp = (den[k]*MAX(qi[k],qmin)) ;
+ temp = sqrt(sqrt(temp*temp*temp)) ;
+#ifdef DEBUGDEBUG
+ xni[k] = 1.e3 ;
+#else
+ xni[k] = MIN(MAX(5.38e7*temp,1.e3),1.e6) ;
+#endif
+ }
+ rmstep = 1./mstep ;
+
+ int n ;
+ float dtcldden, coeres, rdelz ;
+
+
+ float den_k, falk1_k, falk1_kp1, fall1_k, fall1_kp1, delz_k, delz_kp1 ;
+ float falk2_k, falk2_kp1, fall2_k, fall2_kp1 ;
+
+ for ( n = 1 ; n <= mstep ; n++ ) {
+ k = kpe - 1 ;
+ den_k = den[k] ;
+ falk1_kp1 = den_k*qr[k]*w1[k]*rmstep ;
+ fall1_kp1 = falk1_kp1 ;
+ falk2_kp1 = den_k*qs[k]*w2[k]*rmstep ;
+ fall2_kp1 = falk2_kp1 ;
+ dtcldden = dtcld/den_k ;
+ qr[k] = MAX(qr[k]-falk1_kp1*dtcldden,0.0) ;
+ qs[k] = MAX(qs[k]-falk2_kp1*dtcldden,0.0) ;
+ delz_kp1 = delz[k] ;
+ for ( k = kpe-2 ; k >= kps-1 ; k-- ) {
+ den_k = den[k] ;
+ falk1_k = den_k*qr[k]*w1[k]*rmstep ;
+ fall1_k = falk1_k ;
+ falk2_k = den_k*qs[k]*w2[k]*rmstep ;
+ fall2_k = falk2_k ;
+ dtcldden = dtcld/den_k ;
+ delz_k = delz[k] ;
+ rdelz = 1./delz_k ;
+ qr[k] = MAX(qr[k]- (falk1_k-falk1_kp1*delz_kp1*rdelz)* dtcldden,0.) ;
+ qs[k] = MAX(qs[k]- (falk2_k-falk2_kp1*delz_kp1*rdelz)* dtcldden,0.) ;
+ delz_kp1 = delz_k ;
+ falk1_kp1 = falk1_k ;
+ fall1_kp1 = fall1_k ;
+ falk2_kp1 = falk2_k ;
+ fall2_kp1 = fall2_k ;
+ }
+
+ for ( k = kpe-1 ; k >= kps-1 ; k-- ) {
+ if ( t[k] > t0c && qs[k] > 0.) {
+ xlf = xlf0 ;
+ w3[k] = VENFAC(p[k],t[k],den[k]) ;
+ coeres = rslope2s[k]*sqrt(rslopes[k]*rslopebs[2]) ;
+ psmlt[k] = XKA(t[k],den[k])/xlf*(t0c-t[k])*pi/2.
+ *n0sfac[k]*(precs1*rslope2s[k]+precs2
+ *w3[k]*coeres) ;
+ psmlt[k] = MIN(MAX(psmlt[k]*dtcld*rmstep,-qs[k]*rmstep),0.) ;
+ qs[k] += psmlt[k] ;
+ qr[k] -= psmlt[k] ;
+ t[k] += xlf/CPMCAL(q[k])*psmlt[k] ;
+ }
+ }
+ }
+
+//---------------------------------------------------------------
+// Vice [ms-1] : fallout of ice crystal [HDC 5a]
+//---------------------------------------------------------------
+ mstep = 1 ;
+ numdt = 1 ;
+ for ( k = kpe-1 ; k >= kps-1 ; k-- ) {
+ if (qi[k] <= 0.) {
+ w2[k] = 0. ;
+ } else {
+ xmi = den[k]*qi[k]/xni[k] ;
+ diameter = MAX(MIN(dicon * sqrt(xmi),dimax), 1.e-25) ;
+ w1[k] = 1.49e4*exp(log(diameter)*(1.31)) ;
+ w2[k] = w1[k]/delz[k] ;
+ }
+ numdt = MAX( (int) trunc(w2[k]*dtcld+.5+.5),1) ;
+ if(numdt > mstep) mstep = numdt ;
+ }
+ rmstep = 1./mstep ;
+
+ float falkc_k, falkc_kp1, fallc_k, fallc_kp1 ;
+ for ( n = 1 ; n <= mstep ; n++ ) {
+ k = kpe - 1 ;
+ den_k = den[k] ;
+ falkc_kp1 = den_k*qi[k]*w2[k]*rmstep ;
+ fallc_kp1 = fallc_kp1+falkc_kp1 ;
+ qi[k] = MAX(qi[k]-falkc_kp1*dtcld/den_k,0.) ;
+ delz_kp1 = delz[k] ;
+ for ( k = kpe-2 ; k >= kps-1 ; k-- ) {
+ den_k = den[k] ;
+ falkc_k = den_k*qi[k]*w2[k]*rmstep ;
+ fallc_k = fallc_k+falkc_k ;
+ delz_k = delz[k] ;
+ qi[k] = MAX(qi[k]-(falkc_k-falkc_kp1
+ *delz_kp1/delz_k)*dtcld/den_k,0.) ;
+ delz_kp1 = delz_k ;
+ falkc_kp1 = falkc_k ;
+ fallc_kp1 = fallc_k ;
+ }
+ }
+ float fallsum = fall1_k+fall2_k+fallc_k ;
+ float fallsum_qsi = fall2_k+fallc_k ;
+
+ rainncv = 0. ;
+ if(fallsum > 0.) {
+ rainncv = fallsum*delz[1]/denr*dtcld*1000. ;
+ rain = fallsum*delz[1]/denr*dtcld*1000. + rain ;
+ }
+ snowncv = 0. ;
+ if(fallsum_qsi > 0.) {
+ snowncv = fallsum_qsi*delz[0]/denr*dtcld*1000. ;
+ snow = fallsum_qsi*delz[0]/denr*dtcld*1000. + snow ;
+ }
+ sr = 0. ;
+ if ( fallsum > 0. ) sr = fallsum_qsi*delz[0]/denr*dtcld*1000./(rainncv+1.e-12) ;
+
+//---------------------------------------------------------------
+// pimlt: instantaneous melting of cloud ice [HL A47] [RH83 A28]
+// (T>T0: I->C)
+//---------------------------------------------------------------
+
+
+ for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
+
+ // note -- many of these are turned into scalars of form name_reg by _def_ above
+ // so that they will be stored in registers
+ prevp[k] = 0. ;
+ psdep[k] = 0. ;
+ praut[k] = 0. ;
+ psaut[k] = 0. ;
+ pracw[k] = 0. ;
+ psaci[k] = 0. ;
+ psacw[k] = 0. ;
+ pigen[k] = 0. ;
+ pidep[k] = 0. ;
+ pcond[k] = 0. ;
+ psevp[k] = 0. ;
+
+ q_k = q[k] ;
+ t_k = t[k] ;
+ qr_k = qr[k] ;
+ qc_k = qc[k] ;
+ qs_k = qs[k] ;
+ qi_k = qi[k] ;
+ qs1_k = qs1[k] ;
+ qs2_k = qs2[k] ;
+ cpm_k = cpm[k] ;
+ xl_k = xl[k] ;
+
+ float supcol = t0c-t_k ;
+ xlf = xls-xl_k ;
+ if( supcol < 0. ) xlf = xlf0 ;
+ if( supcol < 0 && qi_k > 0. ) {
+ qc_k = qc_k + qi_k ;
+ t_k = t_k - xlf/cpm_k*qi_k ;
+ qi_k = 0. ;
+ }
+//---------------------------------------------------------------
+// pihmf: homogeneous freezing of cloud water below -40c [HL A45]
+// (T<-40C: C->I)
+//---------------------------------------------------------------
+ if( supcol > 40. && qc_k > 0. ) {
+ qi_k = qi_k + qc_k ;
+ t_k = t_k + xlf/cpm_k*qc_k ;
+ qc_k = 0. ;
+ }
+//---------------------------------------------------------------
+// pihtf: heterogeneous freezing of cloud water [HL A44]
+// (T0>T>-40C: C->I)
+//---------------------------------------------------------------
+ if ( supcol > 0. && qc_k > 0.) {
+ float pfrzdtc = MIN(pfrz1*(exp(pfrz2*supcol)-1.)
+ *den[k]/denr/xncr*qc_k*qc_k*dtcld,qc_k) ;
+ qi_k = qi_k + pfrzdtc ;
+ t_k = t_k + xlf/cpm_k*pfrzdtc ;
+ qc_k = qc_k-pfrzdtc ;
+ }
+//---------------------------------------------------------------
+// psfrz: freezing of rain water [HL A20] [LFO 45]
+// (T<T0, R->S)
+//---------------------------------------------------------------
+ if( supcol > 0. && qr_k > 0. ) {
+ float temp = rsloper[k] ;
+ temp = temp*temp*temp*temp*temp*temp*temp ;
+ float pfrzdtr = MIN(20.*(pi*pi)*pfrz1*n0r*denr/den[k]
+ *(exp(pfrz2*supcol)-1.)*temp*dtcld,
+ qr_k) ;
+ qs_k = qs_k + pfrzdtr ;
+ t_k = t_k + xlf/cpm_k*pfrzdtr ;
+ qr_k = qr_k-pfrzdtr ;
+ }
+
+//----------------------------------------------------------------
+// rsloper: reverse of the slope parameter of the rain(m)
+// xka: thermal conductivity of air(jm-1s-1k-1)
+// work1: the thermodynamic term in the denominator associated with
+// heat conduction and vapor diffusion
+// (ry88, y93, h85)
+// work2: parameter associated with the ventilation effects(y93)
+
+ n0sfac[k] = MAX(MIN(exp(alpha*supcol),n0smax/n0s),1.) ;
+ if ( qr_k <= qcrmin ) {
+ rsloper[k] = rslopermax ;
+ rslopebr[k] = rsloperbmax ;
+ rslope2r[k] = rsloper2max ;
+ rslope3r[k] = rsloper3max ;
+ } else {
+ rsloper[k] = 1./(sqrt(sqrt(pidn0r/((qr_k)*(den[k]))))) ;
+ rslopebr[k] = exp(log(rsloper[k])*bvtr) ;
+ rslope2r[k] = rsloper[k] * rsloper[k] ;
+ rslope3r[k] = rslope2r[k] * rsloper[k] ;
+ }
+ if ( qs_k <= qcrmin ) {
+ rslopes[k] = rslopesmax ;
+ rslopebs[k] = rslopesbmax ;
+ rslope2s[k] = rslopes2max ;
+ rslope3s[k] = rslopes3max ;
+ } else {
+ rslopes[k] = 1./(sqrt(sqrt(pidn0s*(n0sfac[k])/((qs_k)*(den[k]))))) ;
+ rslopebs[k] = exp(log(rslopes[k])*bvts) ;
+ rslope2s[k] = rslopes[k] * rslopes[k] ;
+ rslope3s[k] = rslope2s[k] * rslopes[k] ;
+ }
+
+ w1_k = DIFFAC(xl_k,p[k],t_k,den[k],qs1_k) ;
+ w2_k = DIFFAC(xls,p[k],t_k,den[k],qs2_k) ;
+ w3_k = VENFAC(p[k],t_k,den[k]) ;
+
+//
+//===============================================================
+//
+// warm rain processes
+//
+// - follows the processes in RH83 and LFO except for autoconcersion
+//
+//===============================================================
+//
+ float supsat = MAX(q_k,qmin)-qs1_k ;
+ float satdt = supsat/dtcld ;
+//---------------------------------------------------------------
+// praut: auto conversion rate from cloud to rain [HDC 16]
+// (C->R)
+//---------------------------------------------------------------
+ if(qc_k > qc0) {
+ praut[k] = qck1*exp(log(qc_k)*((7./3.))) ;
+ praut[k] = MIN(praut[k],qc_k/dtcld) ;
+ }
+//---------------------------------------------------------------
+// pracw: accretion of cloud water by rain [HL A40] [LFO 51]
+// (C->R)
+//---------------------------------------------------------------
+ if(qr_k > qcrmin && qc_k > qmin) {
+ pracw[k] = MIN(pacrr*rslope3r[k]*rslopebr[k]
+ *qc_k*denfac[k],qc_k/dtcld) ;
+ }
+//---------------------------------------------------------------
+// prevp: evaporation/condensation rate of rain [HDC 14]
+// (V->R or R->V)
+//---------------------------------------------------------------
+ if(qr_k > 0.) {
+ coeres = rslope2r[k]*sqrt(rsloper[k]*rslopebr[k]) ;
+ prevp[k] = (rh1[k]-1.)*(precr1*rslope2r[k]
+ +precr2*w3_k*coeres)/w1_k ;
+ if(prevp[k] < 0.) {
+ prevp[k] = MAX(prevp[k],-qr_k/dtcld) ;
+ prevp[k] = MAX(prevp[k],satdt/2) ;
+ } else {
+ prevp[k] = MIN(prevp[k],satdt/2) ;
+ }
+ }
+
+//
+//===============================================================
+//
+// cold rain processes
+//
+// - follows the revised ice microphysics processes in HDC
+// - the processes same as in RH83 and RH84 and LFO behave
+// following ice crystal hapits defined in HDC, inclduing
+// intercept parameter for snow (n0s), ice crystal number
+// concentration (ni), ice nuclei number concentration
+// (n0i), ice diameter (d)
+//
+//===============================================================
+//
+ float rdtcld = 1./dtcld ;
+ supsat = MAX(q_k,qmin)-qs2_k ;
+ satdt = supsat/dtcld ;
+ int ifsat = 0 ;
+//-------------------------------------------------------------
+// Ni: ice crystal number concentraiton [HDC 5c]
+//-------------------------------------------------------------
+ float temp = (den[k]*MAX(qi_k,qmin)) ;
+ temp = sqrt(sqrt(temp*temp*temp)) ;
+ xni[k] = MIN(MAX(5.38e7*temp,1.e3),1.e6) ;
+ float eacrs = exp(0.07*(-supcol)) ;
+//-------------------------------------------------------------
+// psacw: Accretion of cloud water by snow [HL A7] [LFO 24]
+// (T<T0: C->S, and T>=T0: C->R)
+//-------------------------------------------------------------
+ if(qs_k > qcrmin && qc_k > qmin) {
+ psacw[k] = MIN(pacrc*n0sfac[k]*rslope3s[k]
+ *rslopebs[k]*qc_k*denfac[k]
+ ,qc_k*rdtcld) ;
+ }
+//
+ if(supcol > 0) {
+ if(qs_k > qcrmin && qi_k > qmin) {
+ xmi = den[k]*qi_k/xni[k] ;
+ diameter = MIN(dicon * sqrt(xmi),dimax) ;
+ vt2i = 1.49e4*pow(diameter,(float)1.31) ;
+ vt2s = pvts*rslopebs[k]*denfac[k] ;
+//-------------------------------------------------------------
+// psaci: Accretion of cloud ice by rain [HDC 10]
+// (T<T0: I->S)
+//-------------------------------------------------------------
+ acrfac = 2.*rslope3s[k]+2.*diameter*rslope2s[k]
+ +diameter*diameter*rslopes[k] ;
+ psaci[k] = pi*qi_k*eacrs*n0s*n0sfac[k]
+ *abs(vt2s-vt2i)*acrfac*.25 ;
+ }
+//-------------------------------------------------------------
+// pidep: Deposition/Sublimation rate of ice [HDC 9]
+// (T<T0: V->I or I->V)
+//-------------------------------------------------------------
+ if(qi_k > 0 && ifsat != 1) {
+ xmi = den[k]*qi_k/xni[k] ;
+ diameter = dicon * sqrt(xmi) ;
+ pidep[k] = 4.*diameter*xni[k]*(rh2[k]-1.)/w2_k ;
+ supice = satdt-prevp[k] ;
+ if(pidep[k] < 0.) {
+ pidep[k] = MAX(MAX(pidep[k],satdt*.5),supice) ;
+ pidep[k] = MAX(pidep[k],-qi_k*rdtcld) ;
+ } else {
+ pidep[k] = MIN(MIN(pidep[k],satdt*.5),supice) ;
+ }
+ if(abs(prevp[k]+pidep[k]) >= abs(satdt)) ifsat = 1 ;
+ }
+//-------------------------------------------------------------
+// psdep: deposition/sublimation rate of snow [HDC 14]
+// (V->S or S->V)
+//-------------------------------------------------------------
+ if( qs_k > 0. && ifsat != 1) {
+ coeres = rslope2s[k]*sqrt(rslopes[k]*rslopebs[k]) ;
+ psdep[k] = (rh2[k]-1.)*n0sfac[k]
+ *(precs1*rslope2s[k]+precs2
+ *w3_k*coeres)/w2_k ;
+ supice = satdt-prevp[k]-pidep[k] ;
+ if(psdep[k] < 0.) {
+ psdep[k] = MAX(psdep[k],-qs_k*rdtcld) ;
+ psdep[k] = MAX(MAX(psdep[k],satdt*.5),supice) ;
+ } else {
+ psdep[k] = MIN(MIN(psdep[k],satdt*.5),supice) ;
+ }
+ if(abs(prevp[k]+pidep[k]+psdep[k]) >= abs(satdt))
+ ifsat = 1 ;
+ }
+//-------------------------------------------------------------
+// pigen: generation(nucleation) of ice from vapor [HL A50] [HDC 7-8]
+// (T<T0: V->I)
+//-------------------------------------------------------------
+ if(supsat > 0 && ifsat != 1) {
+ supice = satdt-prevp[k]-pidep[k]-psdep[k] ;
+ xni0 = 1.e3*exp(0.1*supcol) ;
+ roqi0 = 4.92e-11*exp(log(xni0)*(1.33));
+ pigen[k] = MAX(0.,(roqi0/den[k]-MAX(qi_k,0.))
+ *rdtcld) ;
+ pigen[k] = MIN(MIN(pigen[k],satdt),supice) ;
+ }
+//
+//-------------------------------------------------------------
+// psaut: conversion(aggregation) of ice to snow [HDC 12]
+// (T<T0: I->S)
+//-------------------------------------------------------------
+ if(qi_k > 0.) {
+ qimax = roqimax/den[k] ;
+ psaut[k] = MAX(0.,(qi_k-qimax)*rdtcld) ;
+ }
+ }
+//-------------------------------------------------------------
+// psevp: Evaporation of melting snow [HL A35] [RH83 A27]
+// (T>T0: S->V)
+//-------------------------------------------------------------
+ if(supcol < 0.) {
+ if(qs_k > 0. && rh1[k] < 1.) {
+ psevp[k] = psdep[k]*w2_k/w1_k ;
+ } // asked Jimy about this, 11.6.07, JM
+ psevp[k] = MIN(MAX(psevp[k],-qs_k*rdtcld),0.) ;
+ }
+
+
+//
+//
+//----------------------------------------------------------------
+// check mass conservation of generation terms and feedback to the
+// large scale
+//
+ if(t_k<=t0c) {
+//
+// cloud water
+//
+ value = MAX(qmin,qc_k) ;
+ source = (praut[k]+pracw[k]+psacw[k])*dtcld ;
+ if (source > value) {
+ factor = value/source ;
+ praut[k] = praut[k]*factor ;
+ pracw[k] = pracw[k]*factor ;
+ psacw[k] = psacw[k]*factor ;
+ }
+//
+// cloud ice
+//
+ value = MAX(qmin,qi_k) ;
+ source = (psaut[k]+psaci[k]-pigen[k]-pidep[k])*dtcld ;
+ if (source > value) {
+ factor = value/source ;
+ psaut[k] = psaut[k]*factor ;
+ psaci[k] = psaci[k]*factor ;
+ pigen[k] = pigen[k]*factor ;
+ pidep[k] = pidep[k]*factor ;
+ }
+
+//
+// rain (added for WRFV3.0.1)
+//
+ value = MAX(qmin,qr_k) ;
+ source = (-praut[k]+pracw[k]-prevp[k])*dtcld ;
+ if (source > value) {
+ factor = value/source ;
+ praut[k] = praut[k]*factor ;
+ pracw[k] = pracw[k]*factor ;
+ prevp[k] = prevp[k]*factor ;
+ }
+//
+// snow (added for WRFV3.0.1)
+//
+ value = MAX(qmin,qs_k) ;
+ source = (-psdep[k]+psaut[k]-psaci[k]-psacw[k])*dtcld ;
+ if (source > value) {
+ factor = value/source ;
+ psdep[k] = psdep[k]*factor ;
+ psaut[k] = psaut[k]*factor ;
+ psaci[k] = psaci[k]*factor ;
+ psacw[k] = psacw[k]*factor ;
+ }
+// (end added for WRFV3.0.1)
+
+//
+ w3_k=-(prevp[k]+psdep[k]+pigen[k]+pidep[k]) ;
+// update
+ q_k = q_k+w3_k*dtcld ;
+ qc_k = MAX(qc_k-(praut[k]+pracw[k]+psacw[k])*dtcld,0.) ;
+ qr_k = MAX(qr_k+(praut[k]+pracw[k]+prevp[k])*dtcld,0.) ;
+ qi_k = MAX(qi_k-(psaut[k]+psaci[k]-pigen[k]-pidep[k])*dtcld,0.) ;
+ qs_k = MAX(qs_k+(psdep[k]+psaut[k]+psaci[k]+psacw[k])*dtcld,0.) ;
+ xlf = xls-xl_k ;
+ xlwork2 = -xls*(psdep[k]+pidep[k]+pigen[k])-xl_k*prevp[k]-xlf*psacw[k] ;
+ t_k = t_k-xlwork2/cpm_k*dtcld ;
+ } else {
+//
+// cloud water
+//
+ value = MAX(qmin,qc_k) ;
+ source=(praut[k]+pracw[k]+psacw[k])*dtcld ;
+ if (source > value) {
+ factor = value/source ;
+ praut[k] = praut[k]*factor ;
+ pracw[k] = pracw[k]*factor ;
+ psacw[k] = psacw[k]*factor ;
+ }
+//
+// rain (added for WRFV3.0.1)
+//
+ value = MAX(qmin,qr_k) ;
+ source = (-praut[k]-pracw[k]-prevp[k]-psacw[k])*dtcld ;
+ if (source > value) {
+ factor = value/source ;
+ praut[k] = praut[k]*factor ;
+ pracw[k] = pracw[k]*factor ;
+ prevp[k] = prevp[k]*factor ;
+ psacw[k] = psacw[k]*factor ;
+ }
+// (end added for WRFV3.0.1)
+//
+// snow
+//
+ value = MAX(qcrmin,qs_k) ;
+ source=(-psevp[k])*dtcld ;
+ if (source > value) {
+ factor = value/source ;
+ psevp[k] = psevp[k]*factor ;
+ }
+ w3_k=-(prevp[k]+psevp[k]) ;
+// update
+ q_k = q_k+w3_k*dtcld ;
+ qc_k = MAX(qc_k-(praut[k]+pracw[k]+psacw[k])*dtcld,0.) ;
+ qr_k = MAX(qr_k+(praut[k]+pracw[k]+prevp[k] +psacw[k])*dtcld,0.) ;
+ qs_k = MAX(qs_k+psevp[k]*dtcld,0.) ;
+ xlf = xls-xl_k ;
+ xlwork2 = -xl_k*(prevp[k]+psevp[k]) ;
+ t_k = t_k-xlwork2/cpm_k*dtcld ;
+ }
+//
+// Inline expansion for fpvs
+ cvap = cpv ;
+ ttp=t0c+0.01 ;
+ dldt=cvap-cliq ;
+ xa=-dldt/rv ;
+ xb=xa+hvap/(rv*ttp) ;
+ dldti=cvap-cice ;
+ xai=-dldti/rv ;
+ xbi=xai+hsub/(rv*ttp) ;
+ tr=ttp/t_k ;
+ qs1_k=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr)) ;
+ qs1_k = ep2 * qs1_k / (p[k] - qs1_k) ;
+ qs1_k = MAX(qs1_k,qmin) ;
+//
+//----------------------------------------------------------------
+// pcond: condensational/evaporational rate of cloud water [HL A46] [RH83 A6]
+// if there exists additional water vapor condensated/if
+// evaporation of cloud water is not enough to remove subsaturation
+//
+ w1_k = ((MAX(q_k,qmin)-(qs1_k)) /
+ (1.+(xl_k)*(xl_k)/(rv*(cpm_k))*(qs1_k)/((t_k)*(t_k)))) ;
+ // w3_k = qc_k+w1_k ; NOT USED
+ pcond[k] = MIN(MAX(w1_k/dtcld,0.),MAX(q_k,0.)/dtcld) ;
+ if(qc_k > 0. && w1_k < 0.) {
+ pcond[k] = MAX(w1_k,-qc_k)/dtcld ;
+ }
+ q_k = q_k-pcond[k]*dtcld ;
+ qc_k = MAX(qc_k+pcond[k]*dtcld,0.) ;
+ t_k = t_k+pcond[k]*xl_k/cpm_k*dtcld ;
+//
+//
+//----------------------------------------------------------------
+// padding for small values
+//
+ if(qc_k <= qmin) qc_k = 0.0 ;
+ if(qi_k <= qmin) qi_k = 0.0 ;
+
+ q[k] = q_k ;
+ t[k] = t_k ;
+ qr[k] = qr_k ;
+ qc[k] = qc_k ;
+ qs[k] = qs_k ;
+ qi[k] = qi_k ;
+ qs1[k] = qs1_k ;
+
+ }
+ }
+ for ( k = kps-1 ; k <= kpe-1 ; k++ ) {
+ th[k] = t[k] / pii[k] ;
+ }
+ } // guard
+}
+
+
generated by cgit v1.3 (git 2.53.0) at 2026-07-11 04:12:08 +0000