fast_thermo0 and statistics: calculate qsat only where used

In fast_thermo0 (faster, vector-friendly thermodynamics), the qsat
value at the end is not accurate in the case the grid cell is not
saturated. Rather than doing one more function call to calculate qsat,
don't store qsat at all in the thermodynamics but calculate it as needed in
the statistics routines. The qsat field is not used anywhere
in the calculations, only for output.

After this commit the qsat 3D field can be removed.

src/modfielddump.f90

@@ -267,14 +267,15 @@ end subroutine initfielddump
 !> if lbinary, collect data to truncated (2 byte) integers, and write them to file
 !> if lnetcdf, write to netCDF (as float32).
   subroutine fielddump
-    use modfields, only : u0,v0,w0,thl0,qt0,ql0,sv0,thv0h,thvh,tmp0,rhof,exnf,qsat
+    use modfields, only : u0,v0,w0,thl0,qt0,ql0,sv0,thv0h,thvh,tmp0,rhof,exnf,presf
     use modsurfdata,only : thls,qts,thvs
     use modglobal, only : imax,i1,ih,jmax,j1,jh,k1,rk3step,dzf, &
                           timee,dt_lim,cexpnr,ifoutput,rtimee,cp,tdn,tup
     use modmpi,    only : myid,cmyidx, cmyidy
     use modstat_nc, only : lnetcdf, writestat_nc
     use modmicrodata, only : iqr, imicro, imicro_none, tuprsg, tdnrsg
     use modraddata, only   :lwu,lwd,swu,swd
+    use modthermodynamics, only: qsat_tab
     implicit none
 
     integer(KIND=selected_int_kind(4)), allocatable :: field(:,:,:)
@@ -451,8 +452,16 @@ subroutine fielddump
 
     if (lnetcdf .and. lhus) vars(:,:,:,ind_hus) = qt0(2:i1:ncoarse,2:j1:ncoarse,klow:khigh) - ql0(2:i1:ncoarse,2:j1:ncoarse,klow:khigh)
 
-    if (lnetcdf .and. lhur) vars(:,:,:,ind_hur) = 100 * (qt0(2:i1:ncoarse,2:j1:ncoarse,klow:khigh) - ql0(2:i1:ncoarse,2:j1:ncoarse,klow:khigh)) / &
-         qsat(2:i1:ncoarse,2:j1:ncoarse,klow:khigh)
+    if (lnetcdf .and. lhur) then
+       do k = klow, khigh
+          do j = 2, j1, ncoarse
+             do i = 2, i1, ncoarse
+                vars(i,j,k,ind_hur) = 100 * (qt0(i,j,k) - ql0(i,j,k)) / &
+                     qsat_tab(tmp0(i,j,k), presf(k))
+             end do
+          end do
+       end do
+    end if
 
     if (lnetcdf .and. ltntr) then
        do k = klow,khigh

src/modgenstat.f90

@@ -481,14 +481,15 @@ end subroutine genstat
   subroutine do_genstat
 
     use modfields, only : u0,v0,w0,thl0,qt0,qt0h,e120, &
-                          ql0,ql0h,thl0h,thv0h,sv0,exnf,exnh,qsat,tmp0
+                          ql0,ql0h,thl0h,thv0h,sv0,exnf,exnh,tmp0,presf
     use modsurfdata,only: thls,qts,svs,ustar,thlflux,qtflux,svflux
     use modsubgriddata,only : ekm, ekh, csz
     use modglobal, only : i1,ih,j1,jh,k1,kmax,nsv,dzf,dzh,rlv,rv,rd,cp, &
                           ijtot,cu,cv,iadv_sv,iadv_kappa,eps1,dxi,dyi,tup,tdn
     use modmpi,    only : comm3d,mpi_sum,mpierr,slabsum,D_MPI_ALLREDUCE
     use advec_kappa, only : halflev_kappa
     use modmicrodata, only: tuprsg, tdnrsg, imicro, imicro_sice, imicro_sice2
+    use modthermodynamics, only: qsat_tab
     implicit none
 
 
@@ -691,7 +692,7 @@ subroutine do_genstat
         do  i=2,i1
           thv0(i,j,k) = (thl0(i,j,k)+rlv*ql0(i,j,k)/(cp*exnf(k))) &
                         *(1+(rv/rd-1)*qt0(i,j,k)-rv/rd*ql0(i,j,k))
-          huravl(k) = huravl(k) + 100 * (qt0(i,j,k) - ql0(i,j,k)) / qsat(i,j,k)
+          huravl(k) = huravl(k) + 100 * (qt0(i,j,k) - ql0(i,j,k)) / qsat_tab(tmp0(i,j,k), presf(k))
         enddo
       enddo
     enddo

src/modradfield.f90

@@ -161,12 +161,13 @@ end subroutine radfield
 
 
   subroutine sample_radfield
-    use modfields, only : rhof,qt0,ql0,tmp0,qsat,u0,v0
+    use modfields, only : rhof,qt0,ql0,tmp0,u0,v0,presf
     use modsurfdata, only: qtflux,thlflux
     use modglobal, only: dzf,tup,tdn,i1,j1,kmax,rlv,cp
     use modraddata, only : lwd,lwu,swd,swu,lwdca,lwuca,swdca,swuca,swdir,swdif,&
                             SW_up_TOA,SW_dn_TOA,LW_up_TOA,&
                             SW_up_ca_TOA,LW_up_ca_TOA
+    use modthermodynamics, only : qsat_tab
 
     implicit none
     integer :: i,j,k
@@ -201,7 +202,7 @@ subroutine sample_radfield
              field_2D_mn (i,j,20) = field_2D_mn (i,j,20) + rhof(k) *  ql0(i,j,k) * dzf(k)
              ilratio=max(0.,min(1.,(tmp0(i,j,k)-tdn)/(tup-tdn)))! cloud water vs cloud ice partitioning
              field_2D_mn (i,j,21) = field_2D_mn (i,j,21) + rhof(k) *  ql0(i,j,k) * dzf(k) *(1-ilratio)
-             field_2D_mn (i,j,22) = field_2D_mn (i,j,22) + rhof(k) *  qsat(i,j,k) * dzf(k)
+             field_2D_mn (i,j,22) = field_2D_mn (i,j,22) + rhof(k) *  qsat_tab(tmp0(i,j,k), presf(k)) * dzf(k)
           end do
        end do
     end do

src/modthermodynamics.f90

@@ -637,7 +637,9 @@ subroutine icethermo0_fast
                 ! The following could
                 ! be done on the fly to save
                 ! precious memory
-                qsat(i,j,k) = qsat_
+
+                !qsat(i,j,k) = qsat_ ! qsat_ is not a good approximation when not saturated
+                                     ! but ql is still good in that case.
                 T = exnf(k)*thl0(i,j,k) + (rlv/cp) * ql
                 tmp0(i,j,k) = T