Remove obsolescent and non-standard features in physics modules

phys/module_cu_mskf.F

@@ -1762,8 +1762,8 @@ subroutine mskf_mphy(su,    qu,   mu,   du,   cmel, cmei, zf,  pm,  te,   qe, ep
        mtime=deltat/900._r8
        mtimec=deltat/900._r8
 
-       mtime = AMAX1(1.0,mtime)   !TWG remove time scale limitation from CAM5
-       mtimec = AMAX1(1.0,mtimec)    
+       mtime = max(1.0_r8,mtime)   !TWG remove time scale limitation from CAM5
+       mtimec = max(1.0_r8,mtimec)    
 
 ! conservation of qc
 
@@ -5160,19 +5160,19 @@ SUBROUTINE MSKF_eta_PARA (I, J,                           &
 
                JK = KX-KQ+1
 !               print *,'kf qliq=', QLIQ(KQ)
-                 QLIQ(KQ) = amax1(0.0,zmqliq(1,JK))
-                 QICE(KQ) = amax1(0.0,zmqice(1,JK))
+                 QLIQ(KQ) = max(0._r8,zmqliq(1,JK))
+                 QICE(KQ) = max(0._r8,zmqice(1,JK))
 !TWG 06/14/16
-                 QRAIN(KQ) = amax1(0.0,zmqrain(1,JK))
-                 QSNOW(KQ) = amax1(0.0,zmqsnow(1,JK))
-                 NLIQ(KQ) = amax1(0.0,ncmp(1,JK))
-                 NICE(KQ) = amax1(0.0,nimp(1,JK))
-                 NRAIN(KQ) = amax1(0.0,nrmp(1,JK))
-                 NSNOW(KQ) = amax1(0.0,nsmp(1,JK))
-                 CCN(KQ) = amax1(0.0,zmccn(1,JK))
-                 EFFCH(KQ) = MAX(2.49, MIN(effc(1,JK), 50.))
-                 EFFIH(KQ) = MAX(4.99, MIN(effi(1,JK), 125.))
-                 EFFSH(KQ) = MAX(9.99, MIN(effs(1,JK), 999.))
+                 QRAIN(KQ) = max(0._r8,zmqrain(1,JK))
+                 QSNOW(KQ) = max(0._r8,zmqsnow(1,JK))
+                 NLIQ(KQ) = max(0._r8,ncmp(1,JK))
+                 NICE(KQ) = max(0._r8,nimp(1,JK))
+                 NRAIN(KQ) = max(0._r8,nrmp(1,JK))
+                 NSNOW(KQ) = max(0._r8,nsmp(1,JK))
+                 CCN(KQ) = max(0._r8,zmccn(1,JK))
+                 EFFCH(KQ) = max(2.49_r8, min(effc(1,JK), 50._r8))
+                 EFFIH(KQ) = max(4.99_r8, min(effi(1,JK), 125._r8))
+                 EFFSH(KQ) = max(9.99_r8, min(effs(1,JK), 999._r8))
 ! END TWG        
                   DETLQ(KQ)= QLIQ(KQ)*UDR(KQ)
                   DETIC(KQ)= QICE(KQ)*UDR(KQ)

phys/module_cu_scalesas.F

@@ -4222,7 +4222,7 @@ subroutine mfshalcnv(im,ix,km,delt,delp,prslp,psp,phil,ql,        &
         if(cnvflg(i)) then
           if (gdx(i) < dxcrt) then
             scaldfunc(i) = (1.-sigmagfm(i)) * (1.-sigmagfm(i))
-            scaldfunc(i) = max(min(scaldfunc(i), 1.0), 0._kind_phys)
+            scaldfunc(i) = max(min(scaldfunc(i), 1._kind_phys), 0._kind_phys)
           else
             scaldfunc(i) = 1.0
           endif

phys/module_dust_emis.F

@@ -142,7 +142,7 @@ subroutine bulk_dust_emis (ktau,dt,num_soil_layers,u_phy,v_phy,          &
 
 ! increase erodability where the surface albedo is high to account better for real deserts
       if (erodin .gt. 1.E-8 .AND. albbck(i,j).gt.0.175 .and. vegfra(i,j).lt.12.5) then
-         erodin = MIN(0.5, erodin + 0.1*albbck(i,j))
+         erodin = min(0.5d0, dble(erodin + 0.1*albbck(i,j)))
       endif
 
 !  volumetric soil moisture over porosity
@@ -158,14 +158,14 @@ subroutine bulk_dust_emis (ktau,dt,num_soil_layers,u_phy,v_phy,          &
 ! Case of surface dry enough to erode
       IF (gwet < 0.5) THEN  !  Pete's modified value
 !     IF (gwet < 0.2) THEN
-         u_ts = MAX(0.0D+0,u_ts0*(1.2D+0+2.0D-1*LOG10(MAX(1.0D-3, gwet))))
+         u_ts = max(0.0d0,dble(u_ts0*(1.2d0+2.0d-1*log10(max(1.0d-3, dble(gwet))))))
       ELSE
 ! Case of wet surface, no erosion
          u_ts = 100.0
       END IF
       srce = frac_s(n)*erodin*dxy  ! (m2)
 !     srce = 1.1*erodin*dxy  ! (m2)
-      dsrc = MAX(0.0, ch_dust(n,month)*srce*w10m*w10m *(w10m - u_ts)*dt) ! (kg)
+      dsrc = max(0.0d0, dble(ch_dust(n,month)*srce*w10m*w10m *(w10m - u_ts)*dt)) ! (kg)
 
 ! unit change from kg/timestep/cell to ug/m2/s
 !   totalemis=((dsrc)/dt)*converi/dxy

phys/module_mp_fast_sbm.F

@@ -1978,7 +1978,7 @@ SUBROUTINE JERSUPSAT_KS (DEL1,DEL2,DEL1N,DEL2N, &
   IRI = 1
   IPI = 1
 
-	IF(AMAX1(RW,PW,RI,PI)<=RW_PW_RI_PI_MIN) THEN
+	IF(max(RW,PW,RI,PI)<=RW_PW_RI_PI_MIN) THEN
 
     RW = 0.0
     IRW = 0

phys/module_mp_milbrandt2mom.F

@@ -190,7 +190,7 @@ real FUNCTION gamma(xx)
   gammadp= exp(gammadp)
 
 !!GEM:
-  gamma  = sngl(gammadp)
+  gamma = gammadp
 
  END FUNCTION gamma
 !======================================================================!
@@ -260,7 +260,7 @@ real FUNCTION gammln(xx)
   enddo
 
 !! GEM:
-  gammln= sngl( tmp+log(stp*ser/x)  )
+  gammln=tmp+log(stp*ser/x)
 
  END FUNCTION gammln
 !======================================================================!
@@ -3051,7 +3051,7 @@ SUBROUTINE mp_milbrandt2mom_main(WZ,T,Q,QC,QR,QI,QN,QG,QH,NC,NR,NY,NN,NG,NH,PS,
            LAMr     = 1./iLAMr
           !note: The following coding of 'No_r=...' prevents overflow:
           !No_r     = NR(i,k)*LAMr**(1.+alpha_r))*iGR31
-           No_r     = sngl(dble(NR(i,k))*dble(LAMr)**dble(1.+alpha_r))*iGR31
+           No_r     = dble(NR(i,k))*dble(LAMr)**dble(1.+alpha_r)*iGR31
           !note: There is an error in MY05a_eq(8) for VENTx (corrected in code)
            VENTr    = Avx*GR32*iLAMr**cexr5 + Bvx*ScTHRD*sqrt(gam*afr*iMUkin)*GR17*iLAMr**cexr6
            ABw      = CHLC**2/(Ka*RGASV*T(i,k)**2)+1./(DE(i,k)*(QSW(i,k))*Cdiff)

phys/module_mp_p3.F

@@ -354,7 +354,7 @@ subroutine p3_init(lookup_file_dir,nCat,model,stat,abort_on_err)
  if (global_status == STATUS_ERROR) then
     if (err_abort) then
        print*,'Stopping in P3 init'
-       call flush(6)
+       flush(6)
        stop
     endif
     return
@@ -430,7 +430,7 @@ subroutine p3_init(lookup_file_dir,nCat,model,stat,abort_on_err)
    if (global_status == STATUS_ERROR) then
       if (err_abort) then
          print*,'Stopping in P3 init'
-         call flush(6)
+         flush(6)
          stop
       endif
       return

phys/module_mp_thompson.F

@@ -3008,7 +3008,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
            !           -tpc_wev(idx_d, idx_c, idx_n)*orho*odt)
             prw_vcd(k) = MAX(DBLE(-rc(k)*0.99*orho*odt), prw_vcd(k))
             pnc_wcd(k) = MAX(DBLE(-nc(k)*0.99*orho*odt),                &
-                       -tnc_wev(idx_d, idx_c, idx_n)*orho*odt)
+                         DBLE(-tnc_wev(idx_d, idx_c, idx_n)*orho*odt))
 
            endif
           else

phys/module_ra_goddard.F

@@ -5254,7 +5254,7 @@ subroutine swrad ( np,cosz, pl,ta,wa,oa, fcld,ict,icb, &
 !dir$ vector aligned
  DO ic=1,irestrict
  if(lmask(ic) .eqv. .true.) then
-       df_sub(ic,k) = max(df(ic,k) - df(ic,k-1), 0.)  !df for each layer (remove negative df_sub) 
+       df_sub(ic,k) = max(df(ic,k) - df(ic,k-1), 0._fp_kind)  !df for each layer (remove negative df_sub) 
  endif
  ENDDO
  enddo
@@ -5275,7 +5275,7 @@ subroutine swrad ( np,cosz, pl,ta,wa,oa, fcld,ict,icb, &
 !dir$ vector aligned
  DO ic=1,irestrict
  if(lmask(ic) .eqv. .true.) then
-    flc(ic,k)=max(flc(ic,k)-df_clr(ic,k),0.)  !this filter is for small cosine zenith angle.
+    flc(ic,k)=max(flc(ic,k)-df_clr(ic,k),0._fp_kind)  !this filter is for small cosine zenith angle.
  endif
  ENDDO
  enddo
@@ -5316,7 +5316,7 @@ subroutine swrad ( np,cosz, pl,ta,wa,oa, fcld,ict,icb, &
            i_tau = min(max(int(cld_alb*10.)+1,1),10) !1~10
            ratio = ratio_lut(i_tau,i_cos)
        else !use computed clear and cloudy flux ratio (not fast_overcast)
-           ratio = max(0.01, min(1.,(flx(ic,k)/flc(ic,k))))
+           ratio = max(0.01_fp_kind, min(1._fp_kind,(flx(ic,k)/flc(ic,k))))
        endif
        df_sub(ic,k)  = df_sub(ic,k)*ratio  !compute cloudy-sky df_sub
     enddo !k
@@ -5340,8 +5340,8 @@ subroutine swrad ( np,cosz, pl,ta,wa,oa, fcld,ict,icb, &
 !dir$ vector aligned
  DO ic=1,irestrict
  if(lmask(ic) .eqv. .true.) then
-    flx(ic,k)  = max(flx(ic,k)-df_cld(ic,k) , 0.)  !this max is for small cosz
-    flxd(ic,k) = max(flxd(ic,k)-df_cld(ic,k), 0.)  !this max is for small cosz
+    flx(ic,k)  = max(flx(ic,k)-df_cld(ic,k) , 0._fp_kind)  !this max is for small cosz
+    flxd(ic,k) = max(flxd(ic,k)-df_cld(ic,k), 0._fp_kind)  !this max is for small cosz
     flxu(ic,k) = flx(ic,k)-flxd(ic,k)
  endif
  ENDDO
@@ -5361,7 +5361,7 @@ subroutine swrad ( np,cosz, pl,ta,wa,oa, fcld,ict,icb, &
    if (  (fdirir(ic)-df_cld(ic,np+1)) >= 0. ) then ! normal
     fdirir(ic)=fdirir(ic)-df_cld(ic,np+1)  ! updated 
    else ! if negative, it also reduces diffuse component.
-    fdifir(ic) = MAX(0., fdifir(ic) + (fdirir(ic)-df_cld(ic,np+1)) )
+    fdifir(ic) = max(0._fp_kind, fdifir(ic) + (fdirir(ic)-df_cld(ic,np+1)) )
     fdirir(ic)=0.
    endif
 !
@@ -5717,7 +5717,7 @@ subroutine sw_uvpar (np,wh,oh,dp, &
 
 !-----for direct incident radiation
 !     the effective layer optical properties. eqs. (6.2)-(6.4)  !modify max & min
-           tautob(ic)=tausto(ic) + max(tauclb(ic,k,ib),0.e0)
+           tautob(ic)=tausto(ic) + max(tauclb(ic,k,ib),0._fp_kind)
            ssatob(ic)=min(max((ssatau(ic)+ssacl(ic,k,ib)*tauclb(ic,k,ib))/tautob(ic) , ssa_min), ssa_max)  !SSA of cloud is unity
            asytob(ic)=min(max((asysto(ic)+asycl(ic,k,ib)*ssacl(ic,k,ib)*tauclb(ic,k,ib)) &
                          /(   max(ssatob(ic)*tautob(ic),const_tiny)   ),asy_min), asy_max)
@@ -6176,7 +6176,7 @@ subroutine sw_ir (np,wh,dp, &
 !-----compute reflectance and transmittance of the cloudy portion of a layer
 
 !-----for direct incident radiation. eqs.(6.2)-(6.4)
-           tautob(ic,k)=tausto(ic,k)+max(tauclb(ic,k),0.e0)
+           tautob(ic,k)=tausto(ic,k)+max(tauclb(ic,k),0._fp_kind)
            ssatob(ic,k)=min(max((ssatau(ic,k)+ssacl(ic,k,iv)*tauclb(ic,k))/tautob(ic,k),ssa_min),ssa_max)
            asytob(ic,k)=min(max((asysto(ic,k)+asycl(ic,k,iv)*ssacl(ic,k,iv)*tauclb(ic,k)) &
                       /(  max(ssatob(ic,k)*tautob(ic,k),const_tiny)  ),ssa_min),ssa_max)
@@ -6187,7 +6187,7 @@ subroutine sw_ir (np,wh,dp, &
              ssatof(ic,k)=ssatob(ic,k)
              asytof(ic,k)=asytob(ic,k)
            else
-             tautof(ic,k)=tausto(ic,k)+max(tauclf(ic,k),0.e0)
+             tautof(ic,k)=tausto(ic,k)+max(tauclf(ic,k),0._fp_kind)
              ssatof(ic,k)=min(max((ssatau(ic,k)+ssacl(ic,k,iv)*tauclf(ic,k))/tautof(ic,k),ssa_min),ssa_max) 
              asytof(ic,k)=min(max((asysto(ic,k)+asycl(ic,k,iv)*ssacl(ic,k,iv)*tauclf(ic,k)) &
                         /(  max(ssatof(ic,k)*tautof(ic,k),const_tiny)  ),asy_min),asy_max)
@@ -6529,7 +6529,7 @@ subroutine cloud_scale (np,cosz,fcld,taucld,ict,icb, &
 !-----normalize cloud cover following eq. (7.8)
            fa=fcld(ic,k)/cc(ic,kk)
 !-----table look-up
-           taux=min(taux,32.)
+           taux=min(taux,32._fp_kind)
            fm=cosz(ic)/dm
            ft=(log10(taux)-t1)/dt
            fa=fa/da
@@ -6555,8 +6555,8 @@ subroutine cloud_scale (np,cosz,fcld,taucld,ict,icb, &
            xai=xai+(-caib(im,it,ia-1)*(1.-fa)+ &
            caib(im,it,ia+1)*(1.+fa))*fa*.5+caib(im,it,ia)*(1.-fa*fa)
            xai= xai-2.*caib(im,it,ia)
-           xai=max(xai,0.0)
-           xai=min(xai,1.0)
+           xai=max(xai,0._fp_kind)
+           xai=min(xai,1._fp_kind)
            tauclb(ic,k) = taux*xai
 !-----scale cloud optical thickness for diffuse radiation following eq. (7.4)
 !     the scaling factor, xai, is a function of the cloud optical
@@ -6566,8 +6566,8 @@ subroutine cloud_scale (np,cosz,fcld,taucld,ict,icb, &
            xai=xai+(-caif(it,ia-1)*(1.-fa)+ &
             caif(it,ia+1)*(1.+fa))*fa*.5+caif(it,ia)*(1.-fa*fa)
            xai= xai-caif(it,ia)
-           xai=max(xai,0.0)
-           xai=min(xai,1.0)
+           xai=max(xai,0._fp_kind)
+           xai=min(xai,1._fp_kind)
            tauclf(ic,k) = taux*xai
          endif
        endif
@@ -10427,11 +10427,11 @@ subroutine lwrad ( np, emiss, tb, ts, ict, icb,&
 !new
      dp   (ic,k) = pl(ic,k+1)-pl(ic,k)
      dh2o (ic,k) = 1.02*wa(ic,k)*dp(ic,k)
-     dh2o (ic,k) = max(dh2o (ic,k),1.e-30)
+     dh2o (ic,k) = max(dh2o (ic,k),1.e-30_fp_kind)
      do3  (ic,k) = 476.*oa(ic,k)*dp(ic,k)
-     do3  (ic,k) = max(do3 (ic,k),1.e-30)
+     do3  (ic,k) = max(do3 (ic,k),1.e-30_fp_kind)
      dco2 (ic,k) = 789.*co2*dp(ic,k)
-     dco2 (ic,k) = max(dco2 (ic,k),1.e-30)
+     dco2 (ic,k) = max(dco2 (ic,k),1.e-30_fp_kind)
      dch4 (ic,k) = 789.*ch4*dp(ic,k)
      dn2o (ic,k) = 789.*n2o*dp(ic,k)
      df11 (ic,k) = 789.*cfc11*dp(ic,k)
@@ -10576,7 +10576,7 @@ subroutine lwrad ( np, emiss, tb, ts, ict, icb,&
               taux=taux*(1.-ww*ff)
 !-----compute cloud diffuse transmittance. it is approximated by using
 !     a diffusivity factor of 1.66.
-              tauxa=max(0.,1.66*taux)
+              tauxa=max(0._fp_kind,1.66*taux)
               tcldlyr(ic,k)=0.
               if(tauxa.lt.80.) tcldlyr(ic,k)=exp(-tauxa)
            endif
@@ -10885,8 +10885,8 @@ subroutine lwrad ( np, emiss, tb, ts, ict, icb,&
        if (k2 .eq. k1+1) then
 !dir$ vector aligned
        DO ic=1,irestrict
-         yy=min(0.999,trant(ic))
-         yy=max(0.001,yy)
+         yy=min(0.999_fp_kind,trant(ic))
+         yy=max(0.001_fp_kind,yy)
 !-hmhj use log instead of alog for default intrinsic function
          xx=(blevel(ic,k1)-blevel(ic,k2))/ log(yy)
          bu=(blevel(ic,k1)-blevel(ic,k2)*yy)/(1.0-yy)+xx
@@ -11652,9 +11652,9 @@ subroutine tablup(k1,k2,np,nx,nh,sabs,spre,stem,w1,p1, &
 !-----normalize we and pe
         pe=(log10(x2)-p1)/dpe
 !-----restrict the magnitudes of the normalized we and pe.
-        we=min(we,REAL(nh-1))
-        pe=max(pe,0.0)
-        pe=min(pe,REAL(nx-1))
+        we=min(we,real(nh-1,kind=fp_kind))
+        pe=max(pe,0._fp_kind)
+        pe=min(pe,real(nx-1,kind=fp_kind))
 !-----assign iw and ip and compute the distance of we and pe
 !     from iw and ip.
         iw=int(we+1.0)
@@ -11684,8 +11684,8 @@ subroutine tablup(k1,k2,np,nx,nh,sabs,spre,stem,w1,p1, &
         t2 = ca*(1.-fw) + cb*fw
 !-----update the total transmittance between levels k1 and k2
         tran(ic)= (ax + (t1+t2*x3) * x3)*tran(ic)
-        tran(ic)=min(tran(ic),0.9999999)
-        tran(ic)=max(tran(ic),0.0000001)
+        tran(ic)=min(tran(ic),0.9999999_fp_kind)
+        tran(ic)=max(tran(ic),0.0000001_fp_kind)
        else
         tran(ic)=0.9999999
        endif

phys/module_shcu_deng.F

@@ -4478,7 +4478,7 @@ SUBROUTINE TPMIX2(p,thes,tu,qu,qliq,qice,qnewlq,qnewic,XLV1,XLV0,ktau,i,j,nk,tra
        IF(IPTB.GE.220 .OR. IPTB.LE.1 .OR. ITHTB.GE.250 .OR. ITHTB.LE.1)THEN
          write(98,'(a,6i5,2f9.2,i3)')'*** OUT OF BOUNDS ***', ktau,i,j,nk,IPTB, ITHTB, &
                     p/100.0, thes, tracker
-         call flush(98)
+         flush(98)
        ENDIF
 !
       t00=ttab(ithtb  ,iptb  )
@@ -4796,7 +4796,7 @@ SUBROUTINE TPMIX2DD(p,thes,ts,qs,ktau,i,j,nk)
       ithtb=int(tth)+1
        IF(IPTB.GE.220 .OR. IPTB.LE.1 .OR. ITHTB.GE.250 .OR. ITHTB.LE.1)THEN
          write(97,'(a,6i5,2f9.2)')'*** OUT OF BOUNDS ***', ktau,i,j,nk,IPTB, ITHTB, p/100., thes
-         call flush(97)
+         flush(97)
        ENDIF
 !
       t00=ttab(ithtb  ,iptb  )