diff --git a/plasim/src/fluxmod.f90 b/plasim/src/fluxmod.f90
index 27e41da5cef6d9f3eeb98de0b5300c6b40fde33c..552ec478d98d467841b9bbf0f58cbd2a784cc473 100644
--- a/plasim/src/fluxmod.f90
+++ b/plasim/src/fluxmod.f90
@@ -389,10 +389,10 @@
 !     (for energy conservation)
 !
 
-      if(ndheat > 0) then
-       zdtdt(:)=-(zun(:)*zun(:)-du(:,NLEV)*du(:,NLEV)                   &
-                 +zvn(:)*zvn(:)-dv(:,NLEV)*dv(:,NLEV))/deltsec2         &
+      zdtdt(:)=-(zun(:)*zun(:)-du(:,NLEV)*du(:,NLEV)                    &
+                +zvn(:)*zvn(:)-dv(:,NLEV)*dv(:,NLEV))/deltsec2          &
      &          *0.5/acpd/(1.+adv*dq(:,NLEV))
+      if(ndheat > 0) then
        dtdt(:,NLEV)=dtdt(:,NLEV)+zdtdt(:)
       endif
 
@@ -450,20 +450,19 @@
 !     entropy/energy diagnostics
 !
       if(nentropy > 0) then
-       if(nentropy > 2) then
-        dentropy(:,31)=-((zun(:)*zun(:)-du(:,NLEV)*du(:,NLEV))/deltsec2 &
-     &                 +(zvn(:)*zvn(:)-dv(:,NLEV)*dv(:,NLEV))/deltsec2) &
-     &                *0.5*dp(:)/ga*dsigma(NLEV)/dt(:,NLEV) 
-       else
-        dentropy(:,31)=-((zun(:)*zun(:)-du(:,NLEV)*du(:,NLEV))/deltsec2 &
-     &                 +(zvn(:)*zvn(:)-dv(:,NLEV)*dv(:,NLEV))/deltsec2) &
-     &                *0.5*dentrop(:)/ga*dsigma(NLEV)/dentrot(:,NLEV) 
-       endif 
+       dentropy(:,31)=zdtdt(:)/dentrot(:,NLEV)*dentrop(:)*dsigma(NLEV)  &
+     &               /ga
+       if(nentro3d > 0) then
+        dentro3d(:,1:NLEM,21)=0.
+        dentro3d(:,NLEV,21)=dentropy(:,31)
+       endif
       endif
       if(nenergy > 0) then
-       denergy(:,21)=-((zun(:)*zun(:)-du(:,NLEV)*du(:,NLEV))/deltsec2   &
-     &               +(zvn(:)*zvn(:)-dv(:,NLEV)*dv(:,NLEV))/deltsec2)   &
-     &              *0.5*dp(:)/ga*dsigma(NLEV) 
+       denergy(:,21)=zdtdt(:)*dp(:)/ga*dsigma(NLEV) 
+       if(nener3d > 0) then
+        dener3d(:,1:NLEM,21)=0.
+        dener3d(:,NLEV,21)=denergy(:,21)
+       endif
       endif
 !
       return
@@ -584,19 +583,21 @@
 !     entropy/energy diagnostics
 !
       if(nentropy > 0) then
-       if(nentropy > 2) then
-        dentropy(:,7)=(ztn(:)-dt(:,NLEV))/deltsec2/dt(:,NLEV)           &
-     &               *acpd*(1.+adv*dq(:,NLEV))*dp(:)/ga*dsigma(NLEV)
-        dentropy(:,14)=dshfl(:)/dt(:,NLEP)
-       else
-        dentropy(:,7)=(ztn(:)-dt(:,NLEV))/deltsec2/dentrot(:,NLEV)      &
+       dentropy(:,7)=(ztn(:)-dt(:,NLEV))/deltsec2/dentrot(:,NLEV)       &
      &        *acpd*(1.+adv*dentroq(:,NLEV))*dentrop(:)/ga*dsigma(NLEV)
-        dentropy(:,14)=dshfl(:)/dt(:,NLEP)
-       endif 
+       dentropy(:,34)=dshfl(:)/dt(:,NLEP)
+       if(nentro3d > 0) then
+        dentro3d(:,1:NLEM,7)=0.
+        dentro3d(:,NLEV,7)=dentropy(:,7)
+       endif
       endif
       if(nenergy > 0) then
        denergy(:,7)=(ztn(:)-dt(:,NLEV))/deltsec2                        &
      &             *acpd*(1.+adv*dq(:,NLEV))*dp(:)/ga*dsigma(NLEV)
+       if(nener3d > 0) then
+        dener3d(:,1:NLEM,7)=0.
+        dener3d(:,NLEV,7)=denergy(:,7)
+       endif  
       endif
 !
       return
@@ -1010,29 +1011,18 @@
        dentropy(:,8)=0.
        dentropy(:,32)=0.
        do jlev=1,NLEV
-        if(nentropy > 2) then
-         dentropy(:,8)=dentropy(:,8)                                    &
-     &                +zdtdt(:,jlev)/dt(:,jlev)                         &
-     &                *acpd*(1.+adv*dq(:,jlev))*dp(:)/ga*dsigma(jlev)
-         dentropy(:,32)=dentropy(:,32)                                  &
-     &                 -((zun(:,jlev)*zun(:,jlev)                       &
+        dentro(:)=zdtdt(:,jlev)/dentrot(:,jlev)                         &
+     &         *acpd*(1.+adv*dentroq(:,jlev))*dentrop(:)/ga*dsigma(jlev) 
+        dentropy(:,8)=dentropy(:,8)+dentro(:)
+        if(nentro3d > 0) dentro3d(:,jlev,8)=dentro(:)
+        dentro(:)=-((zun(:,jlev)*zun(:,jlev)                            &
      &                  -zu(:,jlev)*zu(:,jlev)                          &
      &                  +zvn(:,jlev)*zvn(:,jlev)                        &
      &                  -zv(:,jlev)*zv(:,jlev))/deltsec2                &
      &                  -(zken(:,jlev)-zke(:,jlev))/deltsec2)           &
-     &                 *0.5*dp(:)/ga*dsigma(jlev)/dt(:,jlev)
-        else
-         dentropy(:,8)=dentropy(:,8)                                    &
-     &                +zdtdt(:,jlev)/dentrot(:,jlev)                    &
-     &         *acpd*(1.+adv*dentroq(:,jlev))*dentrop(:)/ga*dsigma(jlev)
-         dentropy(:,32)=dentropy(:,32)                                  &
-     &                 -((zun(:,jlev)*zun(:,jlev)                       &
-     &                  -zu(:,jlev)*zu(:,jlev)                          &
-     &                  +zvn(:,jlev)*zvn(:,jlev)                        &
-     &                  -zv(:,jlev)*zv(:,jlev))/deltsec2                &
-     &                  -(zken(:,jlev)-zke(:,jlev))/deltsec2)           &
-     &                 *0.5*dentrop(:)/ga*dsigma(jlev)/dentrot(:,jlev)  
-        endif
+     &                 *0.5*dentrop(:)/ga*dsigma(jlev)/dentrot(:,jlev)
+        dentropy(:,32)=dentropy(:,32)+dentro(:)
+        if(nentro3d > 0) dentro3d(:,jlev,22)=dentro(:) 
        enddo
       endif
       if(nenergy > 0) then
@@ -1049,6 +1039,16 @@
      &                 -zv(:,jlev)*zv(:,jlev))/deltsec2                 &
      &                 -(zken(:,jlev)-zke(:,jlev))/deltsec2)            &
      &               *0.5*dp(:)/ga*dsigma(jlev)
+        if(nener3d > 0) then
+         dener3d(:,jlev,8)=zdtdt(:,jlev)                                &
+     &                   *acpd*(1.+adv*dq(:,jlev))*dp(:)/ga*dsigma(jlev)
+         dener3d(:,jlev,22)=-((zun(:,jlev)*zun(:,jlev)                  &
+     &                        -zu(:,jlev)*zu(:,jlev)                    &
+     &                        +zvn(:,jlev)*zvn(:,jlev)                  &
+     &                        -zv(:,jlev)*zv(:,jlev))/deltsec2          &
+     &                       -(zken(:,jlev)-zke(:,jlev))/deltsec2)      &
+     &                     *0.5*dp(:)/ga*dsigma(jlev)
+        endif   
        enddo
       endif
 !
diff --git a/plasim/src/miscmod.f90 b/plasim/src/miscmod.f90
index 3e815c0558db9e32b54c025e1665ffe4b8c8c34d..61b2db33e20c95d7a59d6ecf7ab8d7746ffc2ec1 100644
--- a/plasim/src/miscmod.f90
+++ b/plasim/src/miscmod.f90
@@ -277,17 +277,20 @@
 !     entropy/energy diagnostics
 !
       if(nentropy > 0) then
-       if(nentropy > 2) then
-        dentropy(:,6)=zdtdt(:)/dt(:,1)                                  &
-     &               *acpd*(1.+adv*dq(:,1))*dp(:)/ga*dsigma(1)
-       else
-        dentropy(:,6)=zdtdt(:)/dentrot(:,1)                             &
+       dentropy(:,6)=zdtdt(:)/dentrot(:,1)                              &
      &              *acpd*(1.+adv*dentroq(:,1))*dentrop(:)/ga*dsigma(1)
+       if(nentro3d > 0) then
+        dentro3d(:,2:NLEV,6)=0.
+        dentro3d(:,1,6)=dentropy(:,6) 
        endif
       endif
       if(nenergy > 0) then
        denergy(:,6)=zdtdt(:)                                            &
      &              *acpd*(1.+adv*dq(:,1))*dp(:)/ga*dsigma(1)
+       if(nener3d > 0) then
+        dener3d(:,2:NLEV,6)=0.
+        dener3d(:,1,6)=denergy(:,6)
+       endif
       endif
 !
       return
diff --git a/plasim/src/outmod.f90 b/plasim/src/outmod.f90
index 8c5ab08e92adc7eaa514766e6647a6daf72a1215..6b95e11ad910b04bb17716ce89a694ef612dda48 100644
--- a/plasim/src/outmod.f90
+++ b/plasim/src/outmod.f90
@@ -541,11 +541,20 @@
 !     **************************************
 
       if(nentropy > 0) then
-       do jdiag=1,33
+       do jdiag=1,36
         jcode=319+jdiag
+        if(jcode == 333) cycle                      !333 is reserved
         call writegp(40,dentropy(1,jdiag),jcode,0)
        enddo
       end if
+      if(nentro3d > 0) then
+       do jdiag=1,23
+        jcode=419+jdiag
+        do jlev=1,NLEV
+         call writegp(40,dentro3d(1,jlev,jdiag),jcode,jlev)
+        enddo
+       enddo
+      end if
 
 !     *************************************
 !     * energy diagnostics if switched on *
@@ -557,6 +566,14 @@
         call writegp(40,denergy(1,jdiag),jcode,0)
        enddo
       end if
+      if(nener3d > 0) then
+       do jdiag=1,28
+        jcode=459+jdiag
+        do jlev=1,NLEV
+         call writegp(40,dener3d(1,jlev,jdiag),jcode,jlev)
+        enddo
+       enddo
+      end if
 !
       return
       end
diff --git a/plasim/src/plasim.f90 b/plasim/src/plasim.f90
index e05ac193300841f9320a11ae5d478048ec529622..69c988395257cc244c6e27f66302ad2ed1621e56 100644
--- a/plasim/src/plasim.f90
+++ b/plasim/src/plasim.f90
@@ -249,7 +249,9 @@ plasimversion = "https://github.com/Edilbert/PLASIM/ : 15-Dec-2015"
       call mpbci(ndiagsp   ) ! switch for franks spectral diagnostics
       call mpbci(ndiagcf   ) ! switch for cloud forcing diagnostics
       call mpbci(nentropy  ) ! switch for entropy diagnostics
+      call mpbci(nentro3d  ) ! switch for 3d entropy diagnostics
       call mpbci(nenergy  )  ! switch for energy diagnostics
+      call mpbci(nener3d  )  ! switch for 3d energy diagnostics
       call mpbci(ndiaggp3d ) ! no of 3d gp diagnostic arrays
       call mpbci(ndiaggp2d ) ! no of 2d gp diagnostic arrays
       call mpbci(ndiagsp3d ) ! no of 3d sp diagnostic arrays
@@ -350,16 +352,25 @@ plasimversion = "https://github.com/Edilbert/PLASIM/ : 15-Dec-2015"
        dclforc(:,:)=0.
       end if
       if(nentropy > 0) then
-       allocate(dentropy(NHOR,33))
+       allocate(dentropy(NHOR,36))
        allocate(dentrot(NHOR,NLEV))
        allocate(dentroq(NHOR,NLEV))
        allocate(dentrop(NHOR))
+       allocate(dentro(NHOR))
        dentropy(:,:)=0.
       end if
+      if(nentro3d > 0) then
+       allocate(dentro3d(NHOR,NLEV,23))
+       dentro3d(:,:,:)=0.
+      end if
       if(nenergy > 0) then
        allocate(denergy(NHOR,28))
        denergy(:,:)=0.
       end if
+      if(nener3d > 0) then
+       allocate(dener3d(NHOR,NLEV,28))
+       dener3d(:,:,:)=0.
+      end if
 
       call legini
 
@@ -602,8 +613,12 @@ plasimversion = "https://github.com/Edilbert/PLASIM/ : 15-Dec-2015"
       if(ndiaggp3d > 0) deallocate(dgp3d)
       if(ndiagsp3d > 0) deallocate(dsp3d)
       if(ndiagcf   > 0) deallocate(dclforc)
-      if(nentropy  > 0) deallocate(dentropy,dentrop,dentrot,dentroq)
+      if(nentropy  > 0) then 
+       deallocate(dentropy,dentrop,dentrot,dentroq,dentro)
+      endif
       if(nenergy   > 0) deallocate(denergy)
+      if(nener3d   > 0) deallocate(dener3d)
+      if(nentro3d  > 0) deallocate(dentro3d)
 !
 !     close output file
 !
@@ -947,7 +962,7 @@ plasimversion = "https://github.com/Edilbert/PLASIM/ : 15-Dec-2015"
                    , ndel    , ndheat  , ndiag   , ndiagcf , ndiaggp    &
                    , ndiaggp2d , ndiaggp3d                              &
                    , ndiagsp   , ndiagsp2d , ndiagsp3d                  &
-                   , ndl     , nentropy, neqsig  , nflux                &
+                   , ndl     , nentropy, nentro3d, neqsig  , nflux      &
                    , ngui    , nguidbg , nhdiff  , nhordif , nkits      &
                    , noutput    &
                    , npackgp , npacksp , nperpetual        , nprhor     &
@@ -962,7 +977,7 @@ plasimversion = "https://github.com/Edilbert/PLASIM/ : 15-Dec-2015"
                    , tdissd  , tdissz  , tdisst  , tdissq  , tgr        &
                    , psurf                                              &
                    , restim  , t0      , tfrc                           &
-                   , sigh     , nenergy  , nsponge  , dampsp
+                   , sigh    , nenergy , nener3d , nsponge , dampsp
 !
 !     preset namelist parameter according to model set up
 !
@@ -2064,11 +2079,11 @@ plasimversion = "https://github.com/Edilbert/PLASIM/ : 15-Dec-2015"
        do jlev=1,NLEV
         dentrot(:,jlev)=ct*(gt(:,jlev)+t0(jlev))
         dentroq(:,jlev)=gq(:,jlev)*psurf/dentrop(:)
-        dentropy(:,1)=dentropy(:,1)                                     &
-     &       +acpd*(1.+adv*dentroq(:,jlev))*dentrop(:)/ga*dsigma(jlev)  &
-     &       *log(dentrot(:,jlev)*(zp00/(dentrop(:)*sigma(jlev)))**akap)
+        dentro(:)=acpd*(1.+adv*dentroq(:,jlev))*dentrop(:)*dsigma(jlev) &
+     &    *log(dentrot(:,jlev)*(zp00/(dentrop(:)*sigma(jlev)))**akap)/ga
+        dentropy(:,1)=dentropy(:,1)+dentro(:)
+        if(nentro3d > 0) dentro3d(:,jlev,1)=dentro(:)
        enddo
-       if(nentropy==2) dentrot(:,:)=1.
       endif
 !
 !     save u, v, and ps (at time t) for tracer transport
@@ -2530,54 +2545,29 @@ plasimversion = "https://github.com/Edilbert/PLASIM/ : 15-Dec-2015"
 !
 !     entropy/energy diagnostics
 !
-      if(nentropy > 0) then
+      if(nenergy > 0 .or. nentropy > 0) then
        allocate(ztt(NESP,NLEV))
        allocate(zttgp(NHOR,NLEV))
-       call mpgallsp(ztt,stt,NLEV)
-       ztt(:,:)=ztt(:,:)*ct*ww
-       call sp2fl(ztt,zttgp,NLEV)
-       call fc2gp(zttgp,NLON,NLPP*NLEV)
-       deallocate(ztt)
-       dentropy(:,2)=0.
-       do jlev=1,NLEV
-        dentropy(:,2)=dentropy(:,2)                                     &
-     &         +zttgp(:,jlev)/dentrot(:,jlev)                           &
-     &         *acpd*(1.+adv*dentroq(:,jlev))*dentrop(:)/ga*dsigma(jlev)
-       enddo
-       deallocate(zttgp)
-      endif
-      if(nenergy > 0) then
-       allocate(ztt(NESP,NLEV))
        allocate(zst(NESP,NLEV))
-       allocate(zsd(NESP,NLEV))
-       allocate(zsz(NESP,NLEV))
        if (nqspec == 1) allocate(zsq(NESP,NLEV))
        allocate(zsp(NESP))
-       allocate(zugp(NHOR,NLEV))
-       allocate(zvgp(NHOR,NLEV))
        allocate(ztgp(NHOR,NLEV))
        allocate(zqgp(NHOR,NLEV))
        allocate(zqmgp(NHOR,NLEV))
        allocate(zpgp(NHOR))
        allocate(zpmgp(NHOR))
-       allocate(zekin(NHOR,NLEV))
-       allocate(zepot(NHOR,NLEV))
-       allocate(zttgp(NHOR,NLEV))
-       call mpgallsp(zst,atm,NLEV) 
-       call mpgallsp(zsd,adm,NLEV)
-       call mpgallsp(zsz,azm,NLEV)
-       if (nqspec == 1) call mpgallsp(zsq,aqm,NLEV)
        call mpgallsp(ztt,stt,NLEV)
+       call mpgallsp(zst,atm,NLEV)
+       if (nqspec == 1) call mpgallsp(zsq,aqm,NLEV)
        call mpgallsp(zsp,apm,1)
-       call dv2uv(zsd,zsz,zugp,zvgp)
+       ztt(:,:)=ztt(:,:)*ct*ww
+       call sp2fl(ztt,zttgp,NLEV)
        if (nqspec == 1) call sp2fl(sq,zqgp,NLEV)
        if (nqspec == 1) call sp2fl(zsq,zqmgp,NLEV)
        call sp2fl(zst,ztgp,NLEV)
-       call sp2fl(ztt,zttgp,NLEV)
        call sp2fl(zsp,zpmgp,1)
        call sp2fl(sp,zpgp,1)
-       call fc2gp(zugp,NLON,NLPP*NLEV)
-       call fc2gp(zvgp,NLON,NLPP*NLEV)
+       call fc2gp(zttgp,NLON,NLPP*NLEV)
        if (nqspec == 1) then
           call fc2gp(zqgp,NLON,NLPP*NLEV)
           call fc2gp(zqmgp,NLON,NLPP*NLEV)
@@ -2585,16 +2575,61 @@ plasimversion = "https://github.com/Edilbert/PLASIM/ : 15-Dec-2015"
           zqgp(:,:) = 0.0
           zqmgp(:,:) = 0.0
        endif
-       call fc2gp(zttgp,NLON,NLPP*NLEV)
        call fc2gp(ztgp,NLON,NLPP*NLEV)
        call fc2gp(zpgp,NLON,NLPP)
        call fc2gp(zpmgp,NLON,NLPP)
        zpmgp(:)=psurf*exp(zpmgp(:))
        zpgp(:)=psurf*exp(zpgp(:))
-       zttgp(:,:)=ct*ww*zttgp(:,:)
        do jlev=1,NLEV
         ztgp(:,jlev)=ct*(ztgp(:,jlev)+t0(jlev))
         if (nqspec == 1) zqmgp(:,jlev)=zqmgp(:,jlev)*psurf/zpmgp(:)
+        if (nqspec == 1) zqgp(:,jlev)=zqgp(:,jlev)*psurf/zpgp(:)
+       enddo
+       deallocate(ztt)
+       deallocate(zst)
+       if (nqspec == 1) deallocate(zsq)
+       deallocate(zsp)
+      endif       
+      if(nentropy > 0) then
+       dentropy(:,2)=0.
+       dentropy(:,35)=0.
+       dentropy(:,36)=0.
+       do jlev=1,NLEV
+        dentro(:)=zttgp(:,jlev)/dentrot(:,jlev)                         &
+     &         *acpd*(1.+adv*dentroq(:,jlev))*dentrop(:)/ga*dsigma(jlev)&
+     &         -gascon*(1.+(1./rdbrv-1.)*dentroq(:,jlev))               &
+     &         *(zpgp(:)-zpmgp(:))/deltsec2/dentrop(:)                  &
+     &         *dentrop(:)*dsigma(jlev)/ga
+        dentropy(:,2)=dentropy(:,2)+dentro(:)
+        if(nentro3d > 0) dentro3d(:,jlev,2)=dentro(:)
+        dentro(:)=zttgp(:,jlev)/dentrot(:,jlev)                         &
+     &         *acpd*(1.+adv*dentroq(:,jlev))*dentrop(:)/ga*dsigma(jlev)
+        dentropy(:,35)=dentropy(:,35)+dentro(:)
+        dentropy(:,36)=dentropy(:,36)                                   &
+     &                +((acpd*(1.+adv*zqgp(:,jlev))                     & 
+     &                  *log(ztgp(:,jlev)+zttgp(:,jlev)*deltsec2)       &
+     &                  -gascon*(1.+(1./rdbrv-1.)*zqgp(:,jlev))         &
+     &                  *log(zpgp(:)))*zpgp(:)                          &
+     &                 -(acpd*(1.+adv*zqmgp(:,jlev))                    &
+     &                  *log(ztgp(:,jlev))                              &
+     &                  -gascon*(1.+(1./rdbrv-1.)*zqmgp(:,jlev))        &
+     &                  *log(zpmgp(:)))*zpmgp(:))                       &
+     &                /deltsec2/ga*dsigma(jlev) 
+       enddo
+      endif
+      if(nenergy > 0) then
+       allocate(zsd(NESP,NLEV))
+       allocate(zsz(NESP,NLEV))
+       allocate(zugp(NHOR,NLEV))
+       allocate(zvgp(NHOR,NLEV))
+       allocate(zekin(NHOR,NLEV))
+       allocate(zepot(NHOR,NLEV))
+       call mpgallsp(zsd,adm,NLEV)
+       call mpgallsp(zsz,azm,NLEV)
+       call dv2uv(zsd,zsz,zugp,zvgp)
+       call fc2gp(zugp,NLON,NLPP*NLEV)
+       call fc2gp(zvgp,NLON,NLPP*NLEV)
+       do jlev=1,NLEV
         zekin(:,jlev)=0.5*(zugp(:,jlev)*zugp(:,jlev)                    &
      &                    +zvgp(:,jlev)*zvgp(:,jlev))*cv*cv*rcsq(:)     &
      &               *zpmgp(:)   
@@ -2606,14 +2641,17 @@ plasimversion = "https://github.com/Edilbert/PLASIM/ : 15-Dec-2015"
        denergy(:,27)=0.
        denergy(:,26)=0.
        denergy(:,2)=0.
+       denergy(:,1)=0.
        do jlev=1,NLEV
-        if (nqspec == 1) zqgp(:,jlev)=zqgp(:,jlev)*psurf/zpgp(:)
         zekin(:,jlev)=0.5*(zugp(:,jlev)*zugp(:,jlev)                    &
      &                    +zvgp(:,jlev)*zvgp(:,jlev))*cv*cv*rcsq(:)     &
      &               *zpgp(:)                                           &
      &               -zekin(:,jlev)
         denergy(:,27)=denergy(:,27)                                     &
      &               -zekin(:,jlev)/deltsec2/ga*dsigma(jlev)
+        denergy(:,1)=denergy(:,1)                                       &
+     &              +(ztgp(:,jlev)+zttgp(:,jlev)*deltsec2)              &
+     &              *acpd*(1.+adv*zqgp(:,jlev))*zpgp(:)*dsigma(jlev)/ga
         denergy(:,2)=denergy(:,2)                                       &
      &              +zttgp(:,jlev)                                      &
      &              *acpd*(1.+adv*zqgp(:,jlev))*zpgp(:)/ga*dsigma(jlev)
@@ -2621,22 +2659,30 @@ plasimversion = "https://github.com/Edilbert/PLASIM/ : 15-Dec-2015"
      &               +((ztgp(:,jlev)+zttgp(:,jlev)*deltsec2)            &
      &                 *acpd*(1.+adv*zqgp(:,jlev))*zpgp(:)              &
      &                -zepot(:,jlev))/deltsec2/ga*dsigma(jlev)
+        if(nener3d > 0) then 
+         dener3d(:,jlev,27)=-zekin(:,jlev)/deltsec2/ga*dsigma(jlev)
+         dener3d(:,jlev,1)=(ztgp(:,jlev)+zttgp(:,jlev)*deltsec2)        &
+     &               *acpd*(1.+adv*zqgp(:,jlev))*zpgp(:)*dsigma(jlev)/ga
+         dener3d(:,jlev,2)=zttgp(:,jlev)                                &
+     &              *acpd*(1.+adv*zqgp(:,jlev))*zpgp(:)/ga*dsigma(jlev)
+         dener3d(:,jlev,26)=((ztgp(:,jlev)+zttgp(:,jlev)*deltsec2)      &
+     &                     *acpd*(1.+adv*zqgp(:,jlev))*zpgp(:)          &
+     &                     -zepot(:,jlev))/deltsec2/ga*dsigma(jlev)
+        endif  
        enddo
-       deallocate(ztt)
-       deallocate(zst)
        deallocate(zsd)
        deallocate(zsz)
-       if (nqspec == 1) deallocate(zsq)
-       deallocate(zsp)
        deallocate(zugp)
        deallocate(zvgp)
+       deallocate(zekin)
+       deallocate(zepot)
+      endif
+      if(nenergy > 0 .or. nentropy > 0) then
+       deallocate(zttgp)
        deallocate(zqgp)
        deallocate(zqmgp)
        deallocate(zpgp)
        deallocate(zpmgp)
-       deallocate(zekin)
-       deallocate(zepot)
-       deallocate(zttgp)
        deallocate(ztgp)
       endif
 !
@@ -2852,25 +2898,16 @@ plasimversion = "https://github.com/Edilbert/PLASIM/ : 15-Dec-2015"
        dentropy(:,4)=0.
        dentropy(:,33)=0.
        do jlev=1,NLEV
-        if(nentropy > 2) then
-        dentropy(:,4)=dentropy(:,4)                                     &
-     &         +dtdt(:,jlev)/dt(:,jlev)                                 &
-     &         *acpd*(1.+adv*dq(:,jlev))*dp(:)/ga*dsigma(jlev)
-        zfac=-1.*ww*tfrc(jlev)/(1.+tfrc(jlev)*delt2)
-        dentropy(:,33)=dentropy(:,33)                                   &
-     &                +((du(:,jlev)+dudt(:,jlev)*deltsec2)**2           &
-     &                 +(dv(:,jlev)+dvdt(:,jlev)*deltsec2)**2)          &
-     &                *zfac*dp(:)/ga*dsigma(jlev)/dt(:,jlev)
-        else
-        dentropy(:,4)=dentropy(:,4)                                     &
-     &         +dtdt(:,jlev)/dentrot(:,jlev)                            &
+        dentro(:)=dtdt(:,jlev)/dentrot(:,jlev)                          &
      &         *acpd*(1.+adv*dentroq(:,jlev))*dentrop(:)/ga*dsigma(jlev)
+        dentropy(:,4)=dentropy(:,4)+dentro(:)
+        if(nentro3d > 0) dentro3d(:,jlev,4)=dentro(:)
         zfac=-1.*ww*tfrc(jlev)/(1.+tfrc(jlev)*delt2)
-        dentropy(:,33)=dentropy(:,33)                                   &
-     &                +((du(:,jlev)+dudt(:,jlev)*deltsec2)**2           &
-     &                 +(dv(:,jlev)+dvdt(:,jlev)*deltsec2)**2)          &
-     &                *zfac*dentrop(:)/ga*dsigma(jlev)/dentrot(:,jlev)
-        endif
+        dentro(:)=((du(:,jlev)+dudt(:,jlev)*deltsec2)**2                &
+     &             +(dv(:,jlev)+dvdt(:,jlev)*deltsec2)**2)              &
+     &             *zfac*dentrop(:)/ga*dsigma(jlev)/dentrot(:,jlev)
+        dentropy(:,33)=dentropy(:,33)+dentro(:)
+        if(nentro3d > 0) dentro3d(:,jlev,23)=dentro(:)
        enddo
       endif
       if(nenergy > 0) then
@@ -2879,6 +2916,10 @@ plasimversion = "https://github.com/Edilbert/PLASIM/ : 15-Dec-2015"
         denergy(:,4)=denergy(:,4)                                       &
      &              +dtdt(:,jlev)                                       &
      &              *acpd*(1.+adv*dq(:,jlev))*dp(:)/ga*dsigma(jlev)
+        if(nener3d > 0) then
+         dener3d(:,jlev,4)=dtdt(:,jlev)                                 &
+     &                   *acpd*(1.+adv*dq(:,jlev))*dp(:)/ga*dsigma(jlev)
+        endif
        enddo
       endif
 !
@@ -3207,15 +3248,10 @@ plasimversion = "https://github.com/Edilbert/PLASIM/ : 15-Dec-2015"
        deallocate(ztt)
        dentropy(:,3)=0.
        do jlev=1,NLEV
-        if(nentropy > 2) then
-        dentropy(:,3)=dentropy(:,3)                                     &
-     &         +zttgp(:,jlev)/dt(:,jlev)                                &
-     &         *acpd*(1.+adv*dq(:,jlev))*dp(:)/ga*dsigma(jlev)
-        else
-        dentropy(:,3)=dentropy(:,3)                                     &
-     &         +zttgp(:,jlev)/dentrot(:,jlev)                           &
+        dentro(:)=zttgp(:,jlev)/dentrot(:,jlev)                         &
      &         *acpd*(1.+adv*dentroq(:,jlev))*dentrop(:)/ga*dsigma(jlev)
-        endif
+        dentropy(:,3)=dentropy(:,3)+dentro(:)
+        if(nentro3d > 0) dentro3d(:,jlev,3)=dentro(:) 
        enddo
        deallocate(zttgp)
       endif
@@ -3232,6 +3268,10 @@ plasimversion = "https://github.com/Edilbert/PLASIM/ : 15-Dec-2015"
         denergy(:,3)=denergy(:,3)                                       &
      &              +zttgp(:,jlev)                                      &
      &              *acpd*(1.+adv*dq(:,jlev))*dp(:)/ga*dsigma(jlev)
+        if(nener3d > 0) then
+         dener3d(:,jlev,3)=zttgp(:,jlev)                                &
+     &                   *acpd*(1.+adv*dq(:,jlev))*dp(:)/ga*dsigma(jlev)
+        endif 
        enddo
 
        deallocate(zttgp)
@@ -3409,15 +3449,10 @@ plasimversion = "https://github.com/Edilbert/PLASIM/ : 15-Dec-2015"
        deallocate(ztt)
        dentropy(:,5)=0.
        do jlev=1,NLEV
-        if(nentropy > 2) then
-         dentropy(:,5)=dentropy(:,5)                                    &
-     &                +zttgp(:,jlev)/dt(:,jlev)                         &
-     &                *acpd*(1.+adv*dq(:,jlev))*dp(:)/ga*dsigma(jlev)
-        else
-         dentropy(:,5)=dentropy(:,5)                                    &
-     &                +zttgp(:,jlev)/dentrot(:,jlev)                    &
+        dentro(:)=zttgp(:,jlev)/dentrot(:,jlev)                         &
      &         *acpd*(1.+adv*dentroq(:,jlev))*dentrop(:)/ga*dsigma(jlev)
-        endif
+        dentropy(:,5)=dentropy(:,5)+dentro(:)
+        if(nentro3d > 0) dentro3d(:,jlev,5)=dentro(:)
        enddo
        deallocate(zttgp)
       endif
@@ -3433,6 +3468,10 @@ plasimversion = "https://github.com/Edilbert/PLASIM/ : 15-Dec-2015"
         denergy(:,5)=denergy(:,5)                                       &
      &              +zttgp(:,jlev)                                      &
      &              *acpd*(1.+adv*dq(:,jlev))*dp(:)/ga*dsigma(jlev)
+        if(nener3d > 0) then
+         dener3d(:,jlev,5)=zttgp(:,jlev)                                &
+     &                  *acpd*(1.+adv*dq(:,jlev))*dp(:)/ga*dsigma(jlev)
+        endif 
        enddo
        call mpgallsp(ztt,zsttd,NLEV)
        ztt(:,:)=ztt(:,:)*ct*ww
@@ -3443,6 +3482,10 @@ plasimversion = "https://github.com/Edilbert/PLASIM/ : 15-Dec-2015"
         denergy(:,24)=denergy(:,24)                                     &
      &              +zttgp(:,jlev)                                      &
      &              *acpd*(1.+adv*dq(:,jlev))*dp(:)/ga*dsigma(jlev)
+        if(nener3d > 0) then
+         dener3d(:,jlev,24)=zttgp(:,jlev)                               &
+     &                   *acpd*(1.+adv*dq(:,jlev))*dp(:)/ga*dsigma(jlev)
+        endif 
        enddo
        deallocate(ztt)
        deallocate(zttgp)
@@ -3612,6 +3655,10 @@ plasimversion = "https://github.com/Edilbert/PLASIM/ : 15-Dec-2015"
         denergy(:,23)=denergy(:,23)                                     &
      &               +zdtdt(:,jlev)                                     &
      &               *acpd*(1.+adv*dq(:,jlev))*dp(:)/ga*dsigma(jlev)
+        if(nener3d > 0) then
+         dener3d(:,jlev,23)=zdtdt(:,jlev)                               &
+     &                   *acpd*(1.+adv*dq(:,jlev))*dp(:)/ga*dsigma(jlev)
+        endif
        enddo
        call mpgallsp(zstf2,zstt2,NLEV)
        zstf2(:,:)=zstf2(:,:)*ct*ww
@@ -3622,6 +3669,10 @@ plasimversion = "https://github.com/Edilbert/PLASIM/ : 15-Dec-2015"
         denergy(:,25)=denergy(:,25)                                     &
      &               +zdtdt(:,jlev)                                     &
      &               *acpd*(1.+adv*dq(:,jlev))*dp(:)/ga*dsigma(jlev)
+        if(nener3d > 0) then
+         dener3d(:,jlev,25)=zdtdt(:,jlev)                               &
+     &                   *acpd*(1.+adv*dq(:,jlev))*dp(:)/ga*dsigma(jlev)
+        endif
        enddo
 
       endif
diff --git a/plasim/src/plasimmod.f90 b/plasim/src/plasimmod.f90
index d8dffa551ea26536fb8a5a4dba2dd98dfa2f8bae..a2a29ed6fa668127942094dc36697075df84b49d 100644
--- a/plasim/src/plasimmod.f90
+++ b/plasim/src/plasimmod.f90
@@ -182,7 +182,9 @@
       integer :: nperpetual = 0 ! radiation day for perpetual integration
       integer :: n_sea_points=0 ! number of sea points on grid
       integer :: nentropy = 0   ! switch for entropy diagnostics
+      integer :: nentro3d = 0   ! switch for 3d entropy diagnostics
       integer :: nenergy  = 0   ! switch for energy diagnostics
+      integer :: nener3d  = 0   ! switch for 3d energy diagnostics
       integer :: ndheat   = 1   ! switch for heating due to momentum dissipation
       integer :: nseedlen = 0   ! length of random seed (set by lib call)
       integer :: nsela    = 0   ! enable (1) or disable (0) Semi Lagrangian Advection
@@ -383,12 +385,15 @@
 
       real, allocatable :: dgp2d(:,:),dsp2d(:,:)     ! 2-d diagnostics
       real, allocatable :: dgp3d(:,:,:),dsp3d(:,:,:) ! 3-d diagnostics
-      real, allocatable :: dclforc(:,:)  ! cloud forcing diagnostics
-      real, allocatable :: dentropy(:,:) ! entropy diagnostics
-      real, allocatable :: denergy(:,:)  ! energy diagnostics
-      real, allocatable :: dentrop(:)    ! ps for entropy diagnostics
-      real, allocatable :: dentrot(:,:)  ! t for entropy diagnostics
-      real, allocatable :: dentroq(:,:)  ! q for entropy diagnostics
+      real, allocatable :: dclforc(:,:)   ! cloud forcing diagnostics
+      real, allocatable :: dentropy(:,:)  ! entropy diagnostics
+      real, allocatable :: dentro3d(:,:,:)! entropy diagnostics 3d
+      real, allocatable :: denergy(:,:)   ! energy diagnostics
+      real, allocatable :: dener3d(:,:,:) ! energy diagnostics 3d
+      real, allocatable :: dentrop(:)     ! ps for entropy diagnostics
+      real, allocatable :: dentrot(:,:)   ! t for entropy diagnostics
+      real, allocatable :: dentroq(:,:)   ! q for entropy diagnostics
+      real, allocatable :: dentro(:)      ! 2d entropy for diagnostics
 
 !
 !     accumulated output
diff --git a/plasim/src/radmod.f90 b/plasim/src/radmod.f90
index b90a9e292050650f1e3e1449ea9f535587be581b..12dc6484ae61de6d64c153d1990e7922d764c3d1 100644
--- a/plasim/src/radmod.f90
+++ b/plasim/src/radmod.f90
@@ -574,152 +574,91 @@
        dentropy(:,17)=dswfl(:,NLEP)/dt(:,NLEP)
        dentropy(:,27)=dftd(:,NLEP)/dt(:,NLEP)
        dentropy(:,28)=dftu(:,NLEP)/dt(:,NLEP)
-       dentropy(:,29)=0.
-       dentropy(:,30)=0.
-       do jlev=1,NLEV
-        dentropy(:,29)=dentropy(:,29)+dftu0(:,jlev)/dentrot(:,jlev)
-        dentropy(:,30)=dentropy(:,30)+dftd0(:,jlev)/dentrot(:,jlev)
-       enddo
-       allocate(zdtdte(NHOR,NLEV))
        dentropy(:,9)=0.
-       do jlev=1,NLEV
-        jlep=jlev+1
-        zdtdte(:,jlev)=-ga*(dlwfl(:,jlep)-dlwfl(:,jlev))                &
-     &                /(dsigma(jlev)*dp(:)*acpd*(1.+ADV*dq(:,jlev)))
-        if(nentropy > 2) then
-        dentropy(:,9)=dentropy(:,9)                                     &
-     &         +zdtdte(:,jlev)/dt(:,jlev)                               &
-     &         *acpd*(1.+adv*dq(:,jlev))*dp(:)/ga*dsigma(jlev)
-        else
-        dentropy(:,9)=dentropy(:,9)                                     &
-     &         +zdtdte(:,jlev)/dentrot(:,jlev)                          &
-     &         *acpd*(1.+adv*dentroq(:,jlev))*dentrop(:)/ga*dsigma(jlev)
-        endif
-       enddo
+       dentropy(:,10)=0.
        dentropy(:,21)=0.
-       do jlev=1,NLEV
-        jlep=jlev+1
-        zdtdte(:,jlev)=-ga*(dftd(:,jlep)-dftd(:,jlev))                  &
-     &                /(dsigma(jlev)*dp(:)*acpd*(1.+ADV*dq(:,jlev)))
-        if(nentropy > 2) then
-        dentropy(:,21)=dentropy(:,21)                                   &
-     &         +zdtdte(:,jlev)/dt(:,jlev)                               &
-     &         *acpd*(1.+adv*dq(:,jlev))*dp(:)/ga*dsigma(jlev)
-        else
-        dentropy(:,21)=dentropy(:,21)                                   &
-     &         +zdtdte(:,jlev)/dentrot(:,jlev)                          &
-     &         *acpd*(1.+adv*dentroq(:,jlev))*dentrop(:)/ga*dsigma(jlev)
-        endif
-       enddo
        dentropy(:,22)=0.
-       do jlev=1,NLEV
-        jlep=jlev+1
-        zdtdte(:,jlev)=-ga*(dftu(:,jlep)-dftu(:,jlev))                  &
-     &                /(dsigma(jlev)*dp(:)*acpd*(1.+ADV*dq(:,jlev)))
-        if(nentropy > 2) then
-        dentropy(:,22)=dentropy(:,22)                                   &
-     &         +zdtdte(:,jlev)/dt(:,jlev)                               &
-     &         *acpd*(1.+adv*dq(:,jlev))*dp(:)/ga*dsigma(jlev)
-        else
-        dentropy(:,22)=dentropy(:,22)                                   &
-     &         +zdtdte(:,jlev)/dentrot(:,jlev)                          &
-     &         *acpd*(1.+adv*dentroq(:,jlev))*dentrop(:)/ga*dsigma(jlev)
-        endif
-       enddo
        dentropy(:,23)=0.
-       do jlev=1,NLEV
-        jlep=jlev+1
-        zdtdte(:,jlev)=-ga*(dftue1(:,jlep)-dftue1(:,jlev))              &
-     &                /(dsigma(jlev)*dp(:)*acpd*(1.+ADV*dq(:,jlev)))
-        if(nentropy > 2) then
-        dentropy(:,23)=dentropy(:,23)                                   &
-     &         +zdtdte(:,jlev)/dt(:,jlev)                               &
-     &         *acpd*(1.+adv*dq(:,jlev))*dp(:)/ga*dsigma(jlev)
-        else
-        dentropy(:,23)=dentropy(:,23)                                   &
-     &         +zdtdte(:,jlev)/dentrot(:,jlev)                          &
-     &         *acpd*(1.+adv*dentroq(:,jlev))*dentrop(:)/ga*dsigma(jlev)
-        endif
-       enddo
        dentropy(:,24)=0.
-       do jlev=1,NLEV
-        jlep=jlev+1
-        zdtdte(:,jlev)=-ga*(dftue2(:,jlep)-dftue2(:,jlev))              &
-     &                /(dsigma(jlev)*dp(:)*acpd*(1.+ADV*dq(:,jlev)))
-        if(nentropy > 2) then
-        dentropy(:,24)=dentropy(:,24)                                   &
-     &         +zdtdte(:,jlev)/dt(:,jlev)                               &
-     &         *acpd*(1.+adv*dq(:,jlev))*dp(:)/ga*dsigma(jlev)
-        else
-        dentropy(:,24)=dentropy(:,24)                                   &
-     &         +zdtdte(:,jlev)/dentrot(:,jlev)                          &
-     &         *acpd*(1.+adv*dentroq(:,jlev))*dentrop(:)/ga*dsigma(jlev)
-        endif
-       enddo
        dentropy(:,26)=0.
+       dentropy(:,29)=0.
+       dentropy(:,30)=0.
        do jlev=1,NLEV
-        jlep=jlev+1
-        zdtdte(:,jlev)=-ga*(dftue2(:,jlep)-dftue2(:,jlev))              &
-     &                /(dsigma(jlev)*dp(:)*acpd*(1.+ADV*dq(:,jlev)))
-        dentropy(:,26)=dentropy(:,26)                                   &
-     &              +zdtdte(:,jlev)                                     &
-     &              *acpd*(1.+adv*dq(:,jlev))*dp(:)/ga*dsigma(jlev)
+        jlep=jlev+1  
+        dentro(:)=dftu0(:,jlev)/dentrot(:,jlev)
+        dentropy(:,29)=dentropy(:,29)+dentro(:)
+        if(nentro3d > 0) dentro3d(:,jlev,19)=dentro(:)
+        dentro(:)=dftd0(:,jlev)/dentrot(:,jlev)
+        dentropy(:,30)=dentropy(:,30)+dentro(:)
+        if(nentro3d > 0) dentro3d(:,jlev,20)=dentro(:)
+        dentro(:)=-ga*(dlwfl(:,jlep)-dlwfl(:,jlev))                     &
+     &           /(dsigma(jlev)*dp(:)*acpd*(1.+ADV*dq(:,jlev)))         &
+     &         *acpd*(1.+adv*dentroq(:,jlev))*dentrop(:)/ga*dsigma(jlev)&
+     &         /dentrot(:,jlev)
+        dentropy(:,9)=dentropy(:,9)+dentro(:)
+        if(nentro3d > 0) dentro3d(:,jlev,9)=dentro(:)
+        dentro(:)=-ga*(dswfl(:,jlep)-dswfl(:,jlev))                     &
+     &           /(dsigma(jlev)*dp(:)*acpd*(1.+ADV*dq(:,jlev)))         &
+     &         *acpd*(1.+adv*dentroq(:,jlev))*dentrop(:)/ga*dsigma(jlev)&
+     &         /dentrot(:,jlev)
+        dentropy(:,10)=dentropy(:,10)+dentro(:)
+        if(nentro3d > 0) dentro3d(:,jlev,10)=dentro(:)
+        dentro(:)=-ga*(dftd(:,jlep)-dftd(:,jlev))                       &
+     &           /(dsigma(jlev)*dp(:)*acpd*(1.+ADV*dq(:,jlev)))         &
+     &         *acpd*(1.+adv*dentroq(:,jlev))*dentrop(:)/ga*dsigma(jlev)&
+     &         /dentrot(:,jlev) 
+        dentropy(:,21)=dentropy(:,21)+dentro(:)
+        if(nentro3d > 0) dentro3d(:,jlev,15)=dentro(:)
+        dentro(:)=-ga*(dftu(:,jlep)-dftu(:,jlev))                       &
+     &           /(dsigma(jlev)*dp(:)*acpd*(1.+ADV*dq(:,jlev)))         &
+     &         *acpd*(1.+adv*dentroq(:,jlev))*dentrop(:)/ga*dsigma(jlev)&
+     &         /dentrot(:,jlev)
+        dentropy(:,22)=dentropy(:,22)+dentro(:)
+        if(nentro3d > 0) dentro3d(:,jlev,16)=dentro(:)
+        dentro(:)=-ga*(dftue1(:,jlep)-dftue1(:,jlev))                   &
+     &           /(dsigma(jlev)*dp(:)*acpd*(1.+ADV*dq(:,jlev)))         &
+     &         *acpd*(1.+adv*dentroq(:,jlev))*dentrop(:)/ga*dsigma(jlev)&
+     &         /dentrot(:,jlev)
+        dentropy(:,23)=dentropy(:,23)+dentro(:)
+        if(nentro3d > 0) dentro3d(:,jlev,17)=dentro(:) 
+        dentro(:)=-ga*(dftue2(:,jlep)-dftue2(:,jlev))                   &
+     &           /(dsigma(jlev)*dp(:)*acpd*(1.+ADV*dq(:,jlev)))         &
+     &         *acpd*(1.+adv*dentroq(:,jlev))*dentrop(:)/ga*dsigma(jlev)&
+     &         /dentrot(:,jlev)
+        dentropy(:,24)=dentropy(:,24)+dentro(:)
+        if(nentro3d > 0) dentro3d(:,jlev,18)=dentro(:)
+        dentropy(:,26)=dentropy(:,26)+dentro(:)*dentrot(:,jlev)
        enddo
        dentropy(:,25)=(dftu(:,NLEP)+dentropy(:,26))/dt(:,NLEP)
        dentropy(:,26)=-dentropy(:,26)/dt(:,NLEP)
-!
-       dentropy(:,10)=0.
-       do jlev=1,NLEV
-        jlep=jlev+1
-        zdtdte(:,jlev)=-ga*(dswfl(:,jlep)-dswfl(:,jlev))                &
-     &                /(dsigma(jlev)*dp(:)*acpd*(1.+ADV*dq(:,jlev)))
-        if(nentropy > 2) then
-        dentropy(:,10)=dentropy(:,10)                                   &
-     &         +zdtdte(:,jlev)/dt(:,jlev)                               &
-     &         *acpd*(1.+adv*dq(:,jlev))*dp(:)/ga*dsigma(jlev)
-        else
-        dentropy(:,10)=dentropy(:,10)                                   &
-     &         +zdtdte(:,jlev)/dentrot(:,jlev)                          &
-     &         *acpd*(1.+adv*dentroq(:,jlev))*dentrop(:)/ga*dsigma(jlev)
-        endif
-       enddo
-       deallocate(zdtdte)
       endif
       if(nenergy > 0) then
        allocate(zdtdte(NHOR,NLEV))
        denergy(:,9)=0.
-       do jlev=1,NLEV
-        jlep=jlev+1
-        denergy(:,9)=denergy(:,9)-(dlwfl(:,jlep)-dlwfl(:,jlev))  
-       enddo
+       denergy(:,10)=0.
        denergy(:,17)=0.
-       do jlev=1,NLEV
-        jlep=jlev+1
-        denergy(:,17)=denergy(:,17)-(dftd(:,jlep)-dftd(:,jlev)) 
-       enddo
        denergy(:,18)=0.
-       do jlev=1,NLEV
-        jlep=jlev+1
-        denergy(:,18)=denergy(:,18)-(dftu(:,jlep)-dftu(:,jlev))
-       enddo
        denergy(:,19)=0.
-       do jlev=1,NLEV
-        jlep=jlev+1
-        denergy(:,19)=denergy(:,19)-(dftue1(:,jlep)-dftue1(:,jlev))
-       enddo
        denergy(:,20)=0.
-       do jlev=1,NLEV
-        jlep=jlev+1
-        denergy(:,20)=denergy(:,20)-(dftue2(:,jlep)-dftue2(:,jlev)) 
-       enddo
        denergy(:,28)=0.
-       do jlev=1,NLEV
-        denergy(:,28)=denergy(:,28)+dt(:,jlev)*dsigma(jlev)
-       enddo
-       denergy(:,10)=0.
        do jlev=1,NLEV
         jlep=jlev+1
+        denergy(:,9)=denergy(:,9)-(dlwfl(:,jlep)-dlwfl(:,jlev))  
         denergy(:,10)=denergy(:,10)-(dswfl(:,jlep)-dswfl(:,jlev))
+        denergy(:,17)=denergy(:,17)-(dftd(:,jlep)-dftd(:,jlev)) 
+        denergy(:,18)=denergy(:,18)-(dftu(:,jlep)-dftu(:,jlev))
+        denergy(:,19)=denergy(:,19)-(dftue1(:,jlep)-dftue1(:,jlev))
+        denergy(:,20)=denergy(:,20)-(dftue2(:,jlep)-dftue2(:,jlev)) 
+        denergy(:,28)=denergy(:,28)+dt(:,jlev)*dsigma(jlev)
+        if(nener3d > 0) then
+         dener3d(:,jlev,9)=-(dlwfl(:,jlep)-dlwfl(:,jlev))
+         dener3d(:,jlev,10)=-(dswfl(:,jlep)-dswfl(:,jlev))
+         dener3d(:,jlev,17)=-(dftd(:,jlep)-dftd(:,jlev))
+         dener3d(:,jlev,18)=-(dftu(:,jlep)-dftu(:,jlev))
+         dener3d(:,jlev,19)=-(dftue1(:,jlep)-dftue1(:,jlev))
+         dener3d(:,jlev,20)=-(dftue2(:,jlep)-dftue2(:,jlev))
+         dener3d(:,jlev,28)=dt(:,jlev)*dsigma(jlev)
+        endif
        enddo
        deallocate(zdtdte)
       endif
diff --git a/plasim/src/rainmod.f90 b/plasim/src/rainmod.f90
index 4d22cf1536ff7e6612778b3212160974d947478f..776c6f9d7ae657119d5b4cad579be9b946c2d6be 100644
--- a/plasim/src/rainmod.f90
+++ b/plasim/src/rainmod.f90
@@ -312,11 +312,17 @@
       if(nentropy > 0) then
        dentropy(:,11)=0.
       endif
+      if(nentro3d > 0) then
+       dentro3d(:,:,11)=0.
+      endif
       if(nenergy > 0) then
        denergy(:,11)=0.
        denergy(:,15)=0.
       endif
-
+      if(nener3d > 0) then
+       dener3d(:,:,11)=0.
+       dener3d(:,:,15)=0.
+      endif
 !
 !*    initialize
 !
@@ -495,20 +501,19 @@
 !     entropy/energy diagnostics
 !
        if(nentropy > 0) then
-        if(nentropy > 2) then
-         dentropy(:,11)=dentropy(:,11)                                  &
-     &                 +zdtdt(:)/dt(:,jlev)                             &
-     &                 *acpd*(1.+adv*dq(:,jlev))*dp(:)/ga*dsigma(jlev)
-        else
-         dentropy(:,11)=dentropy(:,11)                                  &
-     &                 +zdtdt(:)/dentrot(:,jlev)                        &
-     &         *acpd*(1.+adv*dentroq(:,jlev))*dentrop(:)/ga*dsigma(jlev)
-        endif
+        dentro(:)=zdtdt(:)/dentrot(:,jlev)                              &
+     &         *acpd*(1.+adv*dentroq(:,jlev))*dentrop(:)/ga*dsigma(jlev) 
+        dentropy(:,11)=dentropy(:,11)+dentro(:)
+        if(nentro3d > 0) dentro3d(:,jlev,11)=dentro(:)
        endif
        if(nenergy > 0) then
         denergy(:,11)=denergy(:,11)                                     &
      &               +zdtdt(:)                                          &
      &               *acpd*(1.+adv*dq(:,jlev))*dp(:)/ga*dsigma(jlev)
+        if(nener3d > 0) then 
+         dener3d(:,jlev,11)=zdtdt(:)                                    &
+     &                   *acpd*(1.+adv*dq(:,jlev))*dp(:)/ga*dsigma(jlev)
+        endif
         do jhor=1,NHOR
          if(zt(jhor) > TMELT) then
           zzal=als
@@ -517,9 +522,12 @@
          endif
          denergy(jhor,15)=denergy(jhor,15)+zzal*zdqdt(jhor)             &
      &                                    *dsigma(jlev)*dp(jhor)/ga
+         if(nener3d > 0) then
+          dener3d(jhor,jlev,15)=zzal*zdqdt(jhor)                        &
+     &                         *dsigma(jlev)*dp(jhor)/ga
+         endif
         enddo
        endif
-
       enddo
 
       return
@@ -1371,32 +1379,26 @@
 !
       if(nentropy > 0) then
        dentropy(:,12)=0.
+       if(nentro3d > 0) dentro3d(:,:,12)=0.
        do jlev=1,NLEV
-        if(nentropy > 2) then
-        where(ilift(:) >= jlev .and. itop(:) <= jlev .and. zi(:) > 0.   &
-     &       .and. ishallow(:) == 0)
-         dentropy(:,12)=dentropy(:,12)                                  &
-     &         +zat(:)*zdtdt(:,jlev)/deltsec2/dt(:,jlev)                &
-     &         *acpd*(1.+adv*dq(:,jlev))*dp(:)/ga*dsigma(jlev)
-        endwhere
-        else 
         where(ilift(:) >= jlev .and. itop(:) <= jlev .and. zi(:) > 0.   &
      &       .and.ishallow(:) == 0)
-         dentropy(:,12)=dentropy(:,12)                                  &
-     &         +zat(:)*zdtdt(:,jlev)/deltsec2/dentrot(:,jlev)           &
+         dentro(:)=+zat(:)*zdtdt(:,jlev)/deltsec2/dentrot(:,jlev)       &
      &         *acpd*(1.+adv*dentroq(:,jlev))*dentrop(:)/ga*dsigma(jlev)
+         dentropy(:,12)=dentropy(:,12)+dentro(:)
         endwhere
+        if(nentro3d > 0) then
+         where(ilift(:) >= jlev .and. itop(:) <= jlev .and. zi(:) > 0.  &
+     &        .and.ishallow(:) == 0)
+          dentro3d(:,jlev,12)=dentro(:)
+         endwhere
         endif
         if(nshallow == 1) then
-         if(nentropy > 2) then
-          dentropy(:,12)=dentropy(:,12)                                 &
-     &                  +zdtdts(:,jlev)/dt(:,jlev)                      &
-     &                  *acpd*(1.+adv*dq(:,jlev))*dp(:)/ga*dsigma(jlev)
-         else
-          dentropy(:,12)=dentropy(:,12)                                 &
-     &                  +zdtdts(:,jlev)/dentrot(:,jlev)                 &
+         dentro(:)=zdtdts(:,jlev)/dentrot(:,jlev)                       &
      &        *acpd*(1.+adv*dentroq(:,jlev))*dentrop(:)/ga*dsigma(jlev)
-         endif
+         dentropy(:,12)=dentropy(:,12)+dentro(:)
+         if(nentro3d > 0) dentro3d(:,jlev,12)=dentro3d(:,jlev,12)       &
+     &                                       +dentro(:)
         endif
        enddo
       endif
@@ -1422,6 +1424,27 @@
      &                *dsigma(jlev)*dp(:)/ga
         endif
        enddo
+       if(nener3d > 0) then
+        dener3d(:,:,12)=0.
+        dener3d(:,:,16)=0.
+        do jlev=1,NLEV
+         where(ilift(:) >= jlev .and. itop(:) <= jlev .and. zi(:) > 0.  &
+     &        .and.ishallow(:) == 0)
+          dener3d(:,jlev,12)=zat(:)*zdtdt(:,jlev)/deltsec2              &
+     &                   *acpd*(1.+adv*dq(:,jlev))*dp(:)/ga*dsigma(jlev)
+         endwhere
+         dener3d(:,jlev,16)=alv*(dqdt(:,jlev)-zdqdtold(:,jlev))         &
+     &                     *dsigma(jlev)*dp(:)/ga
+         if(nshallow == 1) then
+          dener3d(:,jlev,12)=dener3d(:,jlev,12)                         &
+     &                      +zdtdts(:,jlev)/deltsec2                    &
+     &                   *acpd*(1.+adv*dq(:,jlev))*dp(:)/ga*dsigma(jlev)
+          dener3d(:,jlev,16)=dener3d(:,jlev,16)                         &
+     &                      +alv*zdqdts(:,jlev)                         &
+     &                      *dsigma(jlev)*dp(:)/ga
+         endif
+        enddo
+       endif
       endif
 !
       return
@@ -2046,15 +2069,10 @@
       if(nentropy > 0) then
        dentropy(:,13)=0.
        do jlev=1,NLEV
-        if(nentropy > 2) then
-         dentropy(:,13)=dentropy(:,13)                                  &
-     &                 +zdtdt(:,jlev)/dt(:,jlev)                        &
-     &                 *acpd*(1.+adv*dq(:,jlev))*dp(:)/ga*dsigma(jlev)
-        else
-         dentropy(:,13)=dentropy(:,13)                                  &
-     &                 +zdtdt(:,jlev)/dentrot(:,jlev)                   &
-     &         *acpd*(1.+adv*dentroq(:,jlev))*dentrop(:)/ga*dsigma(jlev)
-        endif
+        dentro(:)=zdtdt(:,jlev)/dentrot(:,jlev)                         &
+     &         *acpd*(1.+adv*dentroq(:,jlev))*dentrop(:)/ga*dsigma(jlev) 
+        dentropy(:,13)=dentropy(:,13)+dentro(:)
+        if(nentro3d > 0) dentro3d(:,jlev,13)=dentro(:)
        enddo
       endif
       if(nenergy > 0) then
@@ -2063,6 +2081,10 @@
         denergy(:,13)=denergy(:,13)                                     &
      &               +zdtdt(:,jlev)                                     &
      &               *acpd*(1.+adv*dq(:,jlev))*dp(:)/ga*dsigma(jlev)
+        if(nener3d > 0) then
+         dener3d(:,jlev,13)=zdtdt(:,jlev)                               &
+     &                   *acpd*(1.+adv*dq(:,jlev))*dp(:)/ga*dsigma(jlev)
+        endif
        enddo
       endif
 !
@@ -2240,15 +2262,10 @@
       if(nentropy > 0) then
        dentropy(:,20)=0.
        do jlev=1,NLEV
-        if(nentropy > 2) then
-         dentropy(:,20)=dentropy(:,20)                                  &
-     &                 +zdtdt(:,jlev)/dt(:,jlev)                        &
-     &                 *acpd*(1.+adv*dq(:,jlev))*dp(:)/ga*dsigma(jlev)
-        else
-         dentropy(:,20)=dentropy(:,20)                                  &
-     &                 +zdtdt(:,jlev)/dentrot(:,jlev)                   &
+        dentro(:)=zdtdt(:,jlev)/dentrot(:,jlev)                         &
      &         *acpd*(1.+adv*dentroq(:,jlev))*dentrop(:)/ga*dsigma(jlev)
-        endif
+        dentropy(:,20)=dentropy(:,20)+dentro(:)
+        if(nentro3d > 0) dentro3d(:,jlev,14)=dentro(:)
        enddo
       endif
       if(nenergy > 0) then
@@ -2257,6 +2274,10 @@
          denergy(:,14)=denergy(:,14)                                    &
      &                +zdtdt(:,jlev)                                    &
      &                *acpd*(1.+adv*dq(:,jlev))*dp(:)/ga*dsigma(jlev)
+         if(nener3d > 0) then
+          dener3d(:,jlev,14)=zdtdt(:,jlev)                              &
+     &                   *acpd*(1.+adv*dq(:,jlev))*dp(:)/ga*dsigma(jlev)
+         endif 
        enddo
       endif