update microphysics / thermodynamics for imicro_bulk3

For imicro_bulk3, the saturation adjustment in thermodynamics should be
done with respect to liquid water only. This is done by initializing
the vapor pressure table using liquid water values (for bulk3)
and ice/liquid mix values for all other cases.  Move the
initialization of those tables from modglobal to modthermodynamics,
since the initialization needs imicro and importing that into
modglobal would create a cyclic dependency.

Report an error for running bulk3 with old thermodynamics,
since that has not been updated with the liquid-only special case.

src/modglobal.f90

@@ -235,7 +235,6 @@ subroutine initglobal
 
     integer :: advarr(4)
     real phi, colat, silat, omega, omega_gs
-    real :: ilratio
     integer :: k, n, m, ierr
     character(80) chmess
 
@@ -324,23 +323,6 @@ subroutine initglobal
 
     ! Global constants
 
-
-
-    ! esatltab(m) gives the saturation vapor pressure over water at T corresponding to m
-    ! esatitab(m) is the same over ice
-    ! esatmtab(m) is interpolated between the ice and liquid values with ilratio
-    ! http://www.radiativetransfer.org/misc/atmlabdoc/atmlab/h2o/thermodynamics/e_eq_water_mk.html
-    ! Murphy and Koop 2005 parameterization formula.
-    do m=1,2000
-    ttab(m)=150.+0.2*m
-    esatltab(m)=exp(54.842763-6763.22/ttab(m)-4.21*log(ttab(m))+0.000367*ttab(m)+&
-         tanh(0.0415*(ttab(m)-218.8))*(53.878-1331.22/ttab(m)-9.44523*log(ttab(m))+ 0.014025*ttab(m)))
-
-    esatitab(m)=exp(9.550426-5723.265/ttab(m)+3.53068*log(ttab(m))-0.00728332*ttab(m))
-    ilratio = max(0.,min(1.,(ttab(m)-tdn)/(tup-tdn)))
-    esatmtab(m) = ilratio*esatltab(m) + (1-ilratio)*esatitab(m)
-    end do
-
     mygamma251(-100)=0.
     mygamma21(-100)=0.
     do m=-99,4000

src/modmicrophysics.f90

@@ -37,7 +37,7 @@ module modmicrophysics
 contains
   subroutine initmicrophysics
     use modmpi,   only :myid,comm3d,mpi_integer,mpi_logical, D_MPI_BCAST
-    use modglobal,only :ifnamopt,fname_options,nsv,checknamelisterror
+    use modglobal,only :ifnamopt,fname_options,nsv,checknamelisterror,lfast_thermo
     use modbulkmicro, only : initbulkmicro
     use modsimpleice, only : initsimpleice
     use modsimpleice2, only : initsimpleice2
@@ -59,6 +59,13 @@ subroutine initmicrophysics
       call checknamelisterror(ierr, ifnamopt, 'NAMMICROPHYSICS')
       write(6 ,NAMMICROPHYSICS)
       close(ifnamopt)
+
+      if(imicro == imicro_bulk3) then
+         if (.not. lfast_thermo) then
+            STOP 'Bulkmicro 3 requires lfast_thermo (special treatment, liquid-only saturation adjustment)'
+            ! Bulkmicro 3 special treatment of saturation adjustment is currently only impmemented in the fast_thermo scheme
+         end if
+      end if
     end if
 
     call D_MPI_BCAST(imicro,   1, 0,comm3d,ierr)

src/modthermodynamics.f90

@@ -45,14 +45,36 @@ module modthermodynamics
 
 !> Allocate and initialize arrays
   subroutine initthermodynamics
-    use modglobal, only : ih,i1,jh,j1,k1
-
+    use modglobal, only : ih,i1,jh,j1,k1,tdn,tup,esatltab,esatitab,esatmtab,ttab
+    use modmicrodata, only: imicro,imicro_bulk3
     implicit none
+    real :: ilratio
+    integer :: m
 
     allocate(th0av(k1))
     allocate(thv0(2-ih:i1+ih,2-jh:j1+jh,k1))
     th0av = 0.
 
+    ! esatltab(m) gives the saturation vapor pressure over water at T corresponding to m
+    ! esatitab(m) is the same over ice
+    ! esatmtab(m) is interpolated between the ice and liquid values with ilratio
+    ! http://www.radiativetransfer.org/misc/atmlabdoc/atmlab/h2o/thermodynamics/e_eq_water_mk.html
+    ! Murphy and Koop 2005 parameterization formula.
+    do m=1,2000
+       ttab(m)=150.+0.2*m
+       esatltab(m)=exp(54.842763-6763.22/ttab(m)-4.21*log(ttab(m))+0.000367*ttab(m)+&
+            tanh(0.0415*(ttab(m)-218.8))*(53.878-1331.22/ttab(m)-9.44523*log(ttab(m))+ 0.014025*ttab(m)))
+
+       esatitab(m)=exp(9.550426-5723.265/ttab(m)+3.53068*log(ttab(m))-0.00728332*ttab(m))
+       ilratio = max(0.,min(1.,(ttab(m)-tdn)/(tup-tdn)))
+       if(imicro.eq.imicro_bulk3) then
+          ! bulkmicro3 thermodynamics is for liquid only, ice is explicitely accounted for separately.
+          esatmtab(m) = esatltab(m)
+       else
+          ! for all other microphysics, saturation is w.r.t. liquid and ice
+          esatmtab(m) = ilratio*esatltab(m) + (1-ilratio)*esatitab(m)
+       end if
+    end do
   end subroutine initthermodynamics
 
 !> Do moist thermodynamics.
@@ -527,6 +549,7 @@ pure function esat_tab(T) result(es)
     real(field_r) :: tlo, thi, es
 
     ! interpolated ice-liquid saturation vapor pressure from table
+    ! note if imicto==imicro_bulk3, the table is for liquid only
     tlonr=int((T-150.)*5.)
     tlo = 150. + 0.2*tlonr
     thi = tlo + 0.2