Merge pull request #47 from chowland/moist-interp

Interpolation of initial humidity field

.github/workflows/CI.yml

@@ -16,20 +16,22 @@ jobs:
 
     steps:
       - name: Checkout code
-        uses: actions/checkout@v2
+        uses: actions/checkout@v4
 
       # - name: Install GFortran # NB gfortran and build-essential already on VM
       #   run: |
-      #     sudo apt install build-essential gfortran
+      #     sudo apt-get install build-essential gfortran
 
       - name: Install BLAS and LAPACK
-        run: sudo apt install libblas-dev liblapack-dev
+        run: |
+          sudo apt-get update
+          sudo apt-get install libblas-dev liblapack-dev
 
       - name: Install Parallel HDF5
-        run: sudo apt install libhdf5-mpich-dev
+        run: sudo apt-get install libhdf5-openmpi-dev
 
       - name: Install FFTW library
-        run: sudo apt install libfftw3-dev
+        run: sudo apt-get install libfftw3-dev
 
       - name: Build AFiD
         run: make

Makefile

@@ -43,7 +43,7 @@ ifeq ($(MACHINE),PC)
 endif
 ifeq ($(MACHINE),DISCOVERER)
 	ifeq ($(FLAVOUR),GNU)
-		FC = h5pfc -cpp -fdefault-real-8 -fdefault-double-8
+		FC = h5pfc -O2 -cpp -fdefault-real-8 -fdefault-double-8
 		LDFLAGS += -lfftw3 -llapack -ldl
 	else
 		LDFLAGS += -lfftw3 -qmkl=sequential

src/ReadFlowInterp.F90

@@ -6,7 +6,7 @@ subroutine ReadFlowInterp(prow,pcol)
     use input_grids
     use afid_salinity
     use afid_phasefield
-    use afid_moisture, only: humid
+    use afid_moisture, only: humid, InterpInputHum
     use AuxiliaryRoutines
 
     implicit none
@@ -121,6 +121,8 @@ subroutine ReadFlowInterp(prow,pcol)
 
         call InterpInputVel
 
+        if (moist) call InterpInputHum
+
         if (ismaster) write(*,*) "Velocity and T interpolated"
 
         if (multires) then

src/h5tools/MakeMovieXCut.F90

@@ -30,6 +30,7 @@ subroutine Mkmov_xcut
     ic = 1
     if (RayS < 0 .and. RayT < 0) ic = nxm/2
     if (IBM) ic = nxm/2
+    if (moist) ic = nxm/2
 
     ! Record filename as string
     filename = trim("outputdir/flowmov/movie_xcut.h5")

src/moisture.F90

@@ -534,4 +534,51 @@ subroutine CreateMoistH5Groups(filename)
 
 end subroutine CreateMoistH5Groups
 
+!> Interpolate the humidity field provided by `continua_humid.h5`
+!! from a previous grid onto the current grid specified by `bou.in`
+subroutine InterpInputHum
+    use input_grids
+    use HermiteInterpolations, only: interpolate_xyz_old_to_new
+    real, allocatable, dimension(:,:,:) :: qold !< Humidity on the old grid
+    real :: qup     !< Upper boundary condition for humidity
+    real :: qlo     !< Lower boundary condition for humidity
+    integer :: ic, jc
+
+    ! Read field from h5 file and store in qold
+    call AllocateReal3DArray(qold, -1, nxo+1, &
+            xstarto(2)-lvlhalo, xendo(2)+lvlhalo, &
+            xstarto(3)-lvlhalo, xendo(3)+lvlhalo)
+    qold(:,:,:) = 0.d0
+    call HdfReadContinua(nzo, nyo, nxo, &
+        xstarto(2), xendo(2), xstarto(3), xendo(3), 8, &
+        qold(1:nxo,xstarto(2)-lvlhalo:xendo(2)+lvlhalo,xstarto(3)-lvlhalo:xendo(3)+lvlhalo))
+
+    ! Boundary conditions (called after SetHumidityBCs)
+    ! ASSUMES HUMTP, HUMBP ARE CONSTANT)
+    qup = humtp(1,xstart(2),xstart(3))
+    qlo = humbp(1,xstart(2),xstart(3))
+
+    ! Apply boundary conditions in x-halo cells
+    do ic=xstarto(3),xendo(3)
+        do jc=xstarto(2),xendo(2)
+            if (qfixS==1) then
+                qold(0,jc,ic) = 2.0*qlo - qold(1,jc,ic)
+            else
+                qold(0,jc,ic) = qold(1,jc,ic)
+            end if
+            if (qfixN==1) then
+                qold(nxo,jc,ic) = 2.0*qup - qold(nxmo,jc,ic)
+            else
+                qold(nxo,jc,ic) = qold(nxmo,jc,ic)
+            end if
+        end do
+    end do
+    call update_halo(qold, lvlhalo)
+
+    call interpolate_xyz_old_to_new(qold, humid(1:nxm,:,:))
+
+    call DestroyReal3DArray(qold)
+
+end subroutine InterpInputHum
+
 end module afid_moisture

tools/afidtools/__init__.py

@@ -1,3 +1,4 @@
 from .afidtools import *
 from .xdmf import *
-from .boundary_layers import *
+from .boundary_layers import *
+from .moisture import *

tools/drizzle.py

@@ -1,22 +1,24 @@
 import h5py
 from scipy.special import lambertw
-import os
 import afidtools as afid
 import numpy as np
 
-fld = os.environ['WORK']+'/RainyBenard/drizzle_IC'
+fld = '../examples/RainyBenard'
 inputs = afid.InputParams(fld)
 nxm = inputs.nxm
 print(nxm)
 
-saturated_top = False
+mpar = afid.MoistParams(fld)
+alpha = mpar.alpha
+beta = mpar.beta
+gamma = mpar.gamma
 
-alpha = 3.0
-beta = 2.0
-if saturated_top:
-    gamma = beta/(1 - np.exp(-alpha))
-else:
-    gamma = beta
+saturated_top = (mpar.qfixN==1)
+if gamma < 0:
+    if saturated_top:
+        gamma = beta/(1 - np.exp(-alpha))
+    else:
+        gamma = beta
 print(alpha, beta, gamma)
 
 def W(z):
@@ -46,17 +48,18 @@ def W(z):
 
 if not saturated_top:
     # Calculate height above which humidity not at saturation
-    D = alpha*(beta + 1)*np.exp(1-alpha)/(1 + alpha*beta*np.exp(1-alpha))
-    B = beta*D - 1
-    zc = 1 + 1/alpha/(B-beta)
+    A = beta - (1 + gamma)/(1 + gamma*alpha*np.exp(1-alpha))
+    B = alpha*np.exp(1-alpha)*(1 + gamma)/(1 + gamma*alpha*np.exp(1-alpha))
+    zc = 1 - (1 + gamma*alpha*np.exp(1-alpha))/alpha/(1+gamma)
     # Prescribe linear profiles above this level
-    b[z > zc] = beta - 1 + B*(z[z > zc] - 1)
-    q[z > zc] = D*(1 - z[z > zc])
+    b[z > zc] = beta - 1 + A*(z[z > zc] - 1)
+    q[z > zc] = B*(1 - z[z > zc])
 
 import matplotlib.pyplot as plt
 fig, ax = plt.subplots(layout='constrained')
 ax.plot(b, z)
 ax.plot(q, z)
+ax.grid(True)
 fig.savefig('drizzle.png')
 
 with h5py.File(fld+'/drizzle.h5','w') as f: