summaryrefslogtreecommitdiff
path: root/benchmarks/CUDA/WP
diff options
context:
space:
mode:
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_010bin4530080 -> 0 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, 0 insertions, 5557 deletions
diff --git a/benchmarks/CUDA/WP/README b/benchmarks/CUDA/WP/README
deleted file mode 100644
index 04213b1..0000000
--- a/benchmarks/CUDA/WP/README
+++ /dev/null
@@ -1,62 +0,0 @@
-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
deleted file mode 100644
index d1d223c..0000000
--- a/benchmarks/CUDA/WP/README.GPGPU-Sim
+++ /dev/null
@@ -1,4 +0,0 @@
-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
deleted file mode 100644
index bf3f8d5..0000000
--- a/benchmarks/CUDA/WP/SAMPLE_SESSION
+++ /dev/null
@@ -1,192 +0,0 @@
-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
deleted file mode 100644
index 5da0941..0000000
--- a/benchmarks/CUDA/WP/body_inline.h
+++ /dev/null
@@ -1,741 +0,0 @@
-#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
deleted file mode 100644
index cedbd39..0000000
--- a/benchmarks/CUDA/WP/compare_snaps.F
+++ /dev/null
@@ -1,83 +0,0 @@
-! 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
deleted file mode 100644
index 4453b30..0000000
--- a/benchmarks/CUDA/WP/cublas.h
+++ /dev/null
@@ -1,2 +0,0 @@
-//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
deleted file mode 100644
index 973f550..0000000
--- a/benchmarks/CUDA/WP/d2s.h
+++ /dev/null
@@ -1,36 +0,0 @@
-
-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
deleted file mode 100644
index d1d502b..0000000
--- a/benchmarks/CUDA/WP/data/wsm5_in_010
+++ /dev/null
Binary files differ
diff --git a/benchmarks/CUDA/WP/debug.m4 b/benchmarks/CUDA/WP/debug.m4
deleted file mode 100644
index 0d6ba17..0000000
--- a/benchmarks/CUDA/WP/debug.m4
+++ /dev/null
@@ -1,35 +0,0 @@
-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
deleted file mode 100644
index 344946b..0000000
--- a/benchmarks/CUDA/WP/libmassv.F
+++ /dev/null
@@ -1,385 +0,0 @@
-! 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
deleted file mode 100644
index 835d0a1..0000000
--- a/benchmarks/CUDA/WP/makefile
+++ /dev/null
@@ -1,138 +0,0 @@
-#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
deleted file mode 100644
index 3cfc9d7..0000000
--- a/benchmarks/CUDA/WP/microclock.c
+++ /dev/null
@@ -1,33 +0,0 @@
-#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
deleted file mode 100644
index f7132ba..0000000
--- a/benchmarks/CUDA/WP/module_mp_wsm5.F
+++ /dev/null
@@ -1,1731 +0,0 @@
-#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
deleted file mode 100644
index f270d05..0000000
--- a/benchmarks/CUDA/WP/s2d.h
+++ /dev/null
@@ -1,36 +0,0 @@
-
-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
deleted file mode 100644
index 27a1467..0000000
--- a/benchmarks/CUDA/WP/spt.h
+++ /dev/null
@@ -1,150 +0,0 @@
-// 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
deleted file mode 100755
index 7dc272a..0000000
--- a/benchmarks/CUDA/WP/spt.pl
+++ /dev/null
@@ -1,219 +0,0 @@
-#!/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
deleted file mode 100644
index d0f8a5c..0000000
--- a/benchmarks/CUDA/WP/timings_20071205
+++ /dev/null
@@ -1,24 +0,0 @@
-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
deleted file mode 100644
index 56cb493..0000000
--- a/benchmarks/CUDA/WP/util.h
+++ /dev/null
@@ -1,265 +0,0 @@
- 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
deleted file mode 100644
index 4b629d4..0000000
--- a/benchmarks/CUDA/WP/util4.cu
+++ /dev/null
@@ -1,46 +0,0 @@
-#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
deleted file mode 100644
index 894e937..0000000
--- a/benchmarks/CUDA/WP/wsm5.cu
+++ /dev/null
@@ -1,500 +0,0 @@
-#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
deleted file mode 100644
index 5c435dc..0000000
--- a/benchmarks/CUDA/WP/wsm5_constants.h
+++ /dev/null
@@ -1,92 +0,0 @@
-// 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
deleted file mode 100644
index 7000cb3..0000000
--- a/benchmarks/CUDA/WP/wsm5_gpu.cu
+++ /dev/null
@@ -1,783 +0,0 @@
-#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-10 23:34:02 +0000