Browse Source

Merge branch 'feature/fortran_intel' into develop

lon0
Bill Ladwig 8 years ago
parent
commit
99beada350
  1. 46
      fortran/rip_cape.f90
  2. 18
      fortran/wrf_cloud_fracf.f90
  3. 4
      fortran/wrf_rip_phys_routines.f90
  4. 26
      fortran/wrf_user.f90
  5. 4
      fortran/wrf_user_dbz.f90
  6. 18
      fortran/wrf_user_latlon_routines.f90
  7. 2
      fortran/wrf_vinterp.f90

46
fortran/rip_cape.f90

@ -374,17 +374,17 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& @@ -374,17 +374,17 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
G/(RD*tvirtual(tmk(i,j,kpar1), qvp(i,j,kpar1)))
p2 = MIN(prs(i,j,kpar1)+.5d0*pavg, prsf(i,j,mkzh))
p1 = p2 - pavg
totthe = 0.d0
totqvp = 0.d0
totprs = 0.d0
totthe = 0.D0
totqvp = 0.D0
totprs = 0.D0
DO k = mkzh,2,-1
IF (prsf(i,j,k) .LE. p1) EXIT !GOTO 35
IF (prsf(i,j,k-1) .GE. p2) CYCLE !GOTO 34
p = prs(i,j,k)
pup = prsf(i,j,k)
pdn = prsf(i,j,k-1)
q = MAX(qvp(i,j,k),1.d-15)
th = tmk(i,j,k)*(1000.d0/p)**(GAMMA*(1.d0 + GAMMAMD*q))
q = MAX(qvp(i,j,k),1.D-15)
th = tmk(i,j,k)*(1000.D0/p)**(GAMMA*(1.D0 + GAMMAMD*q))
pp1 = MAX(p1,pdn)
pp2 = MIN(p2,pup)
IF (pp2 .GT. pp1) THEN
@ -398,7 +398,7 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& @@ -398,7 +398,7 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
! 35 CONTINUE
qvppari = totqvp/totprs
tmkpari = (totthe/totprs)*&
(prs(i,j,kpar1)/1000.d0)**(GAMMA*(1.d0+GAMMAMD*qvp(i,j,kpar1)))
(prs(i,j,kpar1)/1000.D0)**(GAMMA*(1.D0+GAMMAMD*qvp(i,j,kpar1)))
END IF
DO kpar = kpar1, kpar2
@ -412,13 +412,13 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& @@ -412,13 +412,13 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
END IF
prspari = prs(i,j,kpar)
ghtpari = ght(i,j,kpar)
gammam = GAMMA * (1.d0 + GAMMAMD*qvppari)
cpm = CP * (1.d0 + CPMD*qvppari)
gammam = GAMMA * (1.D0 + GAMMAMD*qvppari)
cpm = CP * (1.D0 + CPMD*qvppari)
e = MAX(1.d-20,qvppari*prspari/(EPS + qvppari))
e = MAX(1.D-20,qvppari*prspari/(EPS + qvppari))
tlcl = TLCLC1/(LOG(tmkpari**TLCLC2/e) - TLCLC3) + TLCLC4
ethpari = tmkpari*(1000.d0/prspari)**(GAMMA*(1.d0 + GAMMAMD*qvppari))*&
EXP((THTECON1/tlcl - THTECON2)*qvppari*(1.d0 + THTECON3*qvppari))
ethpari = tmkpari*(1000.D0/prspari)**(GAMMA*(1.D0 + GAMMAMD*qvppari))*&
EXP((THTECON1/tlcl - THTECON2)*qvppari*(1.D0 + THTECON3*qvppari))
zlcl = ghtpari + (tmkpari - tlcl)/(G/cpm)
! Calculate buoyancy and relative height of lifted parcel at
@ -476,13 +476,13 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& @@ -476,13 +476,13 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
buoy(kk) = G*(tvlift - tvenv)/tvenv
zrel(kk) = ghtlift - ghtpari
IF ((kk .GT. 1) .AND. (buoy(kk)*buoy(kk-1) .LT. 0.0d0)) THEN
IF ((kk .GT. 1) .AND. (buoy(kk)*buoy(kk-1) .LT. 0.0D0)) THEN
! Parcel ascent curve crosses sounding curve, so create a new level
! in the bottom-up array at the crossing.
kk = kk + 1
buoy(kk) = buoy(kk-1)
zrel(kk) = zrel(kk-1)
buoy(kk-1) = 0.d0
buoy(kk-1) = 0.D0
zrel(kk-1) = zrel(kk-2) + buoy(kk-2)/&
(buoy(kk-2) - buoy(kk))*(zrel(kk) - zrel(kk-2))
END IF
@ -511,11 +511,11 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& @@ -511,11 +511,11 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
! Get the accumulated buoyant energy from the parcel's starting
! point, at all levels up to the top level.
benaccum(1) = 0.0d0
benaccum(1) = 0.0D0
benamin = 9d9
DO k = 2,kmax
dz = zrel(k) - zrel(k-1)
benaccum(k) = benaccum(k-1) + .5d0*dz*(buoy(k-1) + buoy(k))
benaccum(k) = benaccum(k-1) + .5D0*dz*(buoy(k-1) + buoy(k))
IF (benaccum(k) .LT. benamin) THEN
benamin = benaccum(k)
END IF
@ -527,7 +527,7 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& @@ -527,7 +527,7 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
elfound = .FALSE.
DO k = kmax,klcl,-1
IF (buoy(k) .GE. 0.d0) THEN
IF (buoy(k) .GE. 0.D0) THEN
! k of equilibrium level
kel = k
elfound = .TRUE.
@ -549,8 +549,8 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& @@ -549,8 +549,8 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
! to a more appropriate missing value, which is passed into this
! routine as cmsg.
! cape(i,j,kpar) = -0.1d0
! cin(i,j,kpar) = -0.1d0
! cape(i,j,kpar) = -0.1D0
! cin(i,j,kpar) = -0.1D0
IF (.NOT. elfound) THEN
cape(i,j,kpar) = cmsg
cin(i,j,kpar) = cmsg
@ -572,7 +572,7 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& @@ -572,7 +572,7 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
! inhibition (cin).
! First get the lfc according to the above definition.
benamin = 9d9
benamin = 9D9
klfc = kmax
DO k = klcl,kel
IF (benaccum(k) .LT. benamin) THEN
@ -583,8 +583,8 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& @@ -583,8 +583,8 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
! Now we can assign values to cape and cin
cape(i,j,kpar) = MAX(benaccum(kel)-benamin, 0.1d0)
cin(i,j,kpar) = MAX(-benamin, 0.1d0)
cape(i,j,kpar) = MAX(benaccum(kel)-benamin, 0.1D0)
cin(i,j,kpar) = MAX(-benamin, 0.1D0)
! cin is uninteresting when cape is small (< 100 j/kg), so set
! cin to -0.1 (see note about missing values in v6.1.0) in
@ -595,8 +595,8 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& @@ -595,8 +595,8 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
! to a more appropriate missing value, which is passed into this
! routine as cmsg.
! IF (cape(i,j,kpar).lt.100.d0) cin(i,j,kpar) = -0.1d0
IF (cape(i,j,kpar) .LT. 100.d0) cin(i,j,kpar) = cmsg
! IF (cape(i,j,kpar).lt.100.D0) cin(i,j,kpar) = -0.1D0
IF (cape(i,j,kpar) .LT. 100.D0) cin(i,j,kpar) = cmsg
! 102 CONTINUE
END DO

18
fortran/wrf_cloud_fracf.f90

@ -29,11 +29,11 @@ SUBROUTINE DCLOUDFRAC(pres, rh, lowc, midc, highc, nz, ns, ew) @@ -29,11 +29,11 @@ SUBROUTINE DCLOUDFRAC(pres, rh, lowc, midc, highc, nz, ns, ew)
DO k = 1,nz-1
IF (k .GE. kclo .AND. k .LT. kcmi) THEN
lowc(i,j) = AMAX1(rh(i,j,k), lowc(i,j))
lowc(i,j) = MAX(rh(i,j,k), lowc(i,j))
ELSE IF (k .GE. kcmi .AND. k .LT. kchi) THEN ! mid cloud
midc(i,j) = AMAX1(rh(i,j,k), midc(i,j))
midc(i,j) = MAX(rh(i,j,k), midc(i,j))
ELSE if (k .GE. kchi) THEN ! high cloud
highc(i,j) = AMAX1(rh(i,j,k), highc(i,j))
highc(i,j) = MAX(rh(i,j,k), highc(i,j))
END IF
END DO
@ -41,12 +41,12 @@ SUBROUTINE DCLOUDFRAC(pres, rh, lowc, midc, highc, nz, ns, ew) @@ -41,12 +41,12 @@ SUBROUTINE DCLOUDFRAC(pres, rh, lowc, midc, highc, nz, ns, ew)
midc(i,j) = 4.0*midc(i,j)/100. - 3.0
highc(i,j) = 2.5*highc(i,j)/100. - 1.5
lowc(i,j) = amin1(lowc(i,j), 1.0)
lowc(i,j) = amax1(lowc(i,j), 0.0)
midc(i,j) = amin1(midc(i,j), 1.0)
midc(i,j) = amax1(midc(i,j), 0.0)
highc(i,j) = amin1(highc(i,j), 1.0)
highc(i,j) = amax1(highc(i,j), 0.0)
lowc(i,j) = MIN(lowc(i,j), 1.0)
lowc(i,j) = MAX(lowc(i,j), 0.0)
midc(i,j) = MIN(midc(i,j), 1.0)
midc(i,j) = MAX(midc(i,j), 0.0)
highc(i,j) = MIN(highc(i,j), 1.0)
highc(i,j) = MAX(highc(i,j), 0.0)
END DO
END DO

4
fortran/wrf_rip_phys_routines.f90

@ -70,11 +70,11 @@ SUBROUTINE WETBULBCALC(prs, tmk, qvp, twb, nx, ny, nz, psafile, errstat, errmsg) @@ -70,11 +70,11 @@ SUBROUTINE WETBULBCALC(prs, tmk, qvp, twb, nx, ny, nz, psafile, errstat, errmsg)
DO k=1,nz
DO j=1,ny
DO i=1,nx
q = DMAX1(qvp(i,j,k), 1.D-15)
q = MAX(qvp(i,j,k), 1.D-15)
t = tmk(i,j,k)
p = prs(i,j,k)/100.
e = q*p/(EPS + q)
tlcl = TLCLC1/(DLOG(t**TLCLC2/e) - TLCLC3) + TLCLC4
tlcl = TLCLC1/(LOG(t**TLCLC2/e) - TLCLC3) + TLCLC4
eth = t*(1000./p)**(GAMMA*(1. + GAMMAMD*q))*&
EXP((THTECON1/tlcl - THTECON2)*q*(1. + THTECON3*q))

26
fortran/wrf_user.f90

@ -137,8 +137,8 @@ SUBROUTINE DZSTAG(znew, nx, ny, nz, z, nxz, nyz ,nzz, terrain) @@ -137,8 +137,8 @@ SUBROUTINE DZSTAG(znew, nx, ny, nz, z, nxz, nyz ,nzz, terrain)
DO k = 1,nz
DO j = 1,ny
DO i = 1,nx
ii = MIN0(i,nxz)
im1 = MAX0(i-1,1)
ii = MIN(i,nxz)
im1 = MAX(i-1,1)
znew(i,j,k) = 0.5D0*(z(ii,j,k) + z(im1,j,k))
END DO
END DO
@ -147,8 +147,8 @@ SUBROUTINE DZSTAG(znew, nx, ny, nz, z, nxz, nyz ,nzz, terrain) @@ -147,8 +147,8 @@ SUBROUTINE DZSTAG(znew, nx, ny, nz, z, nxz, nyz ,nzz, terrain)
ELSE IF (ny .GT. nyz) THEN
DO k = 1,nz
DO j = 1,NY
jj = MIN0(j,nyz)
jm1 = MAX0(j-1,1)
jj = MIN(j,nyz)
jm1 = MAX(j-1,1)
DO i = 1,nx
znew(i,j,k) = 0.5D0*(z(i,jj,k) + z(i,jm1,k))
END DO
@ -196,8 +196,8 @@ SUBROUTINE DINTERP2DXY(v3d, v2d, xy, nx, ny, nz, nxy) @@ -196,8 +196,8 @@ SUBROUTINE DINTERP2DXY(v3d, v2d, xy, nx, ny, nz, nxy)
REAL(KIND=8) :: w11, w12, w21, w22, wx, wy
DO ij = 1,nxy
i = MAX0(1,MIN0(nx-1,INT(xy(1,ij)+1)))
j = MAX0(1,MIN0(ny-1,INT(xy(2,ij)+1)))
i = MAX(1,MIN(nx-1,INT(xy(1,ij)+1)))
j = MAX(1,MIN(ny-1,INT(xy(2,ij)+1)))
wx = DBLE(i+1) - (xy(1,ij)+1)
wy = DBLE(j+1) - (xy(2,ij)+1)
w11 = wx*wy
@ -255,7 +255,7 @@ SUBROUTINE DINTERP1D(v_in, v_out, z_in, z_out, vmsg, nz_in, nz_out) @@ -255,7 +255,7 @@ SUBROUTINE DINTERP1D(v_in, v_out, z_in, z_out, vmsg, nz_in, nz_out)
DO WHILE ((.NOT. interp) .AND. (kp .GE. 2))
IF (((z_in(kp-im) .LE. height) .AND. (z_in(kp-ip) .GT. height))) THEN
w2 = (height - z_in(kp-im))/(z_in(kp-ip) - z_in(kp-im))
w1 = 1.d0 - w2
w1 = 1.D0 - w2
v_out(k) = w1*v_in(kp-im) + w2*v_in(kp-ip)
interp = .TRUE.
END IF
@ -608,9 +608,9 @@ SUBROUTINE DCOMPUTERH(qv, p, t, rh, nx) @@ -608,9 +608,9 @@ SUBROUTINE DCOMPUTERH(qv, p, t, rh, nx)
es = EZERO*EXP(ESLCON1*(temperature - CELKEL)/(temperature - ESLCON2))
! qvs = ep_2*es/(pressure-es)
qvs = EPS*es/(0.01D0*pressure - (1.D0 - EPS)*es)
! rh = 100*amax1(1., qv(i)/qvs)
! rh = 100*MAX(1., qv(i)/qvs)
! rh(i) = 100.*qv(i)/qvs
rh(i) = 100.D0*DMAX1(DMIN1(qv(i)/qvs, 1.0D0), 0.0D0)
rh(i) = 100.D0*MAX(MIN(qv(i)/qvs, 1.0D0), 0.0D0)
END DO
RETURN
@ -645,7 +645,7 @@ SUBROUTINE DGETIJLATLONG(lat_array, long_array, lat, longitude, ii, jj, nx, ny, @@ -645,7 +645,7 @@ SUBROUTINE DGETIJLATLONG(lat_array, long_array, lat, longitude, ii, jj, nx, ny,
DO i = 1,nx
latd = (lat_array(i,j) - lat)**2
longd = (long_array(i,j) - longitude)**2
! longd = dmin1((long_array(i,j)-longitude)**2, &
! longd = MIN((long_array(i,j)-longitude)**2, &
! (long_array(i,j)+longitude)**2)
dist = SQRT(latd + longd)
IF (dist_min .GT. dist) THEN
@ -792,12 +792,12 @@ SUBROUTINE DCOMPUTETD(td, pressure, qv_in, nx) @@ -792,12 +792,12 @@ SUBROUTINE DCOMPUTETD(td, pressure, qv_in, nx)
INTEGER :: i
DO i = 1,nx
qv = DMAX1(qv_in(i), 0.D0)
qv = MAX(qv_in(i), 0.D0)
! vapor pressure
tdc = qv*pressure(i)/(.622D0 + qv)
! avoid problems near zero
tdc = DMAX1(tdc, 0.001D0)
tdc = MAX(tdc, 0.001D0)
td(i) = (243.5D0*LOG(tdc) - 440.8D0)/(19.48D0 - LOG(tdc))
END DO
@ -834,7 +834,7 @@ SUBROUTINE DCOMPUTEICLW(iclw, pressure, qc_in, nx, ny, nz) @@ -834,7 +834,7 @@ SUBROUTINE DCOMPUTEICLW(iclw, pressure, qc_in, nx, ny, nz)
DO j = 3,ny - 2
DO i = 3,nx - 2
DO k = 1,nz
qclw = DMAX1(qc_in(i,j,k), 0.D0)
qclw = MAX(qc_in(i,j,k), 0.D0)
IF (k.EQ.1) THEN
dp = pressure(i,j,k-1) - pressure(i,j,k)
ELSE IF (k.EQ.nz) then

4
fortran/wrf_user_dbz.f90

@ -136,8 +136,8 @@ SUBROUTINE CALCDBZ(prs, tmk, qvp, qra, qsn, qgr, sn0, ivarint, iliqskin, dbz, nx @@ -136,8 +136,8 @@ SUBROUTINE CALCDBZ(prs, tmk, qvp, qra, qsn, qgr, sn0, ivarint, iliqskin, dbz, nx
! Calculate variable intercept parameters
IF (ivarint .EQ. 1) THEN
temp_c = DMIN1(-0.001D0, tmk(i,j,k)-CELKEL)
sonv = DMIN1(2.0D8, 2.0D6*EXP(-0.12D0*temp_c))
temp_c = MIN(-0.001D0, tmk(i,j,k)-CELKEL)
sonv = MIN(2.0D8, 2.0D6*EXP(-0.12D0*temp_c))
gonv = gon
IF (qgr(i,j,k) .GT. R1) THEN

18
fortran/wrf_user_latlon_routines.f90

@ -142,14 +142,14 @@ SUBROUTINE DLLTOIJ(map_proj, truelat1, truelat2, stdlon, lat1, lon1,& @@ -142,14 +142,14 @@ SUBROUTINE DLLTOIJ(map_proj, truelat1, truelat2, stdlon, lat1, lon1,&
! the rsw tag.
rsw = 0.D0
IF (lat1 .NE. 0.D0) THEN
rsw = (DLOG(TAN(0.5D0*((lat1 + 90.D0)*RAD_PER_DEG))))/dlon
rsw = (LOG(TAN(0.5D0*((lat1 + 90.D0)*RAD_PER_DEG))))/dlon
END IF
deltalon = lon - lon1
IF (deltalon .LT. -180.D0) deltalon = deltalon + 360.D0
IF (deltalon .GT. 180.D0) deltalon = deltalon - 360.D0
i = knowni + (deltalon/(dlon*DEG_PER_RAD))
j = knownj + (DLOG(TAN(0.5D0*((lat + 90.D0)*RAD_PER_DEG))))/dlon - rsw
j = knownj + (LOG(TAN(0.5D0*((lat + 90.D0)*RAD_PER_DEG))))/dlon - rsw
! ps
ELSE IF (map_proj .EQ. 2) THEN
@ -183,9 +183,9 @@ SUBROUTINE DLLTOIJ(map_proj, truelat1, truelat2, stdlon, lat1, lon1,& @@ -183,9 +183,9 @@ SUBROUTINE DLLTOIJ(map_proj, truelat1, truelat2, stdlon, lat1, lon1,&
END IF
IF (ABS(truelat1 - truelat2) .GT. 0.1D0) THEN
cone = (DLOG(COS(truelat1*RAD_PER_DEG))-DLOG(COS(truelat2*RAD_PER_DEG)))/&
(DLOG(TAN((90.D0 - ABS(truelat1))*RAD_PER_DEG*0.5D0))-&
DLOG(TAN((90.D0 - ABS(truelat2))*RAD_PER_DEG*0.5D0)))
cone = (LOG(COS(truelat1*RAD_PER_DEG))-LOG(COS(truelat2*RAD_PER_DEG)))/&
(LOG(TAN((90.D0 - ABS(truelat1))*RAD_PER_DEG*0.5D0))-&
LOG(TAN((90.D0 - ABS(truelat2))*RAD_PER_DEG*0.5D0)))
ELSE
cone = SIN(ABS(truelat1)*RAD_PER_DEG)
END IF
@ -358,7 +358,7 @@ SUBROUTINE DIJTOLL(map_proj, truelat1, truelat2, stdlon, lat1, lon1,& @@ -358,7 +358,7 @@ SUBROUTINE DIJTOLL(map_proj, truelat1, truelat2, stdlon, lat1, lon1,&
! the rsw tag.
rsw = 0.D0
IF (lat1 .NE. 0.D0) THEN
rsw = (DLOG(TAN(0.5D0*((lat1 + 90.D0)*RAD_PER_DEG))))/dlon
rsw = (LOG(TAN(0.5D0*((lat1 + 90.D0)*RAD_PER_DEG))))/dlon
END IF
lat = 2.0D0*ATAN(EXP(dlon*(rsw + aj - knownj)))*DEG_PER_RAD - 90.D0
@ -418,9 +418,9 @@ SUBROUTINE DIJTOLL(map_proj, truelat1, truelat2, stdlon, lat1, lon1,& @@ -418,9 +418,9 @@ SUBROUTINE DIJTOLL(map_proj, truelat1, truelat2, stdlon, lat1, lon1,&
END IF
IF (ABS(truelat1 - truelat2) .GT. 0.1D0) THEN
cone = (DLOG(COS(truelat1*RAD_PER_DEG)) - DLOG(COS(truelat2*RAD_PER_DEG)))/&
(DLOG(TAN((90.D0 - ABS(truelat1))*RAD_PER_DEG*0.5D0)) - &
DLOG(TAN((90.D0 - ABS(truelat2))*RAD_PER_DEG*0.5D0)))
cone = (LOG(COS(truelat1*RAD_PER_DEG)) - LOG(COS(truelat2*RAD_PER_DEG)))/&
(LOG(TAN((90.D0 - ABS(truelat1))*RAD_PER_DEG*0.5D0)) - &
LOG(TAN((90.D0 - ABS(truelat2))*RAD_PER_DEG*0.5D0)))
ELSE
cone = SIN(ABS(truelat1)*RAD_PER_DEG)
END IF

2
fortran/wrf_vinterp.f90

@ -369,7 +369,7 @@ SUBROUTINE wrf_vintrp(datain, dataout, pres, tk, qvp, ght, terrain,& @@ -369,7 +369,7 @@ SUBROUTINE wrf_vintrp(datain, dataout, pres, tk, qvp, ght, terrain,&
ELSE !else for checking above ground
ptarget = psurfsm - 150.D0
dpmin = 1.e4
dpmin = 1.E4
DO k=1,nz
ripk = nz-k+1
dp = ABS((pres(i,j,ripk) * 0.01D0) - ptarget)

Loading…
Cancel
Save