From 569611689af7754c79f9526dfdc600122c97a7bd Mon Sep 17 00:00:00 2001
From: lkugler <lukas.kugler@gmail.com>
Date: Fri, 20 May 2022 17:04:29 +0200
Subject: [PATCH] horizontal shift of initial field
---
dartwrf/create_wbubble_wrfinput.py | 57 ++++++++++++++++++++++++------
1 file changed, 46 insertions(+), 11 deletions(-)
diff --git a/dartwrf/create_wbubble_wrfinput.py b/dartwrf/create_wbubble_wrfinput.py
index f452b34..846cbe4 100644
--- a/dartwrf/create_wbubble_wrfinput.py
+++ b/dartwrf/create_wbubble_wrfinput.py
@@ -1,4 +1,5 @@
import os, sys, shutil
+from re import U
import datetime as dt
import numpy as np
@@ -9,15 +10,19 @@ dx_km = 2
cr = 15 # km horizontal relaxation distance
cz = 2000 # meters vertical relaxation distance
-perturbations = False
-if len(sys.argv) > 1:
- if 'perturb' == sys.argv[1]:
- perturbations = True
+args = sys.argv[1:]
+
+if 'plot' in args:
+ import matplotlib.pyplot as plt
+ fig, ax = plt.subplots(figsize=(5,5))
+
+perturbations = True if 'perturb' in args else False
+print('update state in wrfinput wrfout file to DART background file')
print('perturb wbubble = ', perturbations)
-for iens in range(1, exp.n_ens+1):
- print('update state in wrfinput wrfout file to DART background file')
+for iens in range(1, exp.n_ens+1):
+ print('iens', iens)
wrfin = cluster.wrf_rundir(iens)+'/wrfinput_d01'
# insert wbubble
@@ -31,8 +36,8 @@ for iens in range(1, exp.n_ens+1):
z = z[:, np.newaxis, np.newaxis]
if perturbations:
- cx = (85 + 30*np.random.rand())*dx_km
- cy = (85 + 30*np.random.rand())*dx_km
+ cx = (70 + 60*np.random.rand())*dx_km
+ cy = (70 + 60*np.random.rand())*dx_km
else:
cx = 100*dx_km
cy = 100*dx_km
@@ -44,9 +49,39 @@ for iens in range(1, exp.n_ens+1):
xx, yy = np.meshgrid(dx, dy)
dr = np.sqrt(xx**2 + yy**2)[np.newaxis, :, :]
- pert = 4*np.exp(-(dr/cr)**2)*np.exp(-(z/cz)**2)
+ amplitude = 3
+ if perturbations:
+ amplitude += np.random.rand()*2 - 1 # Uniformly random +/- 1 K
+ increment = amplitude*np.exp(-(dr/cr)**2)*np.exp(-(z/cz)**2)
- ds.variables['T'][0,...] += pert
- ds.variables['THM'][0,...] += pert
+ # now perturbations are centered in domain
+ # to shift it in direction NW
+ shift = 50
+ incr2 = increment.copy()*0.
+ incr2[:, :-shift, :-shift] = increment[:, shift:, shift:] # main area
+ incr2[:, -shift:, :-shift] = increment[:, :shift, shift:] # lower part
+ incr2[:, :-shift, -shift:] = increment[:, shift:, :shift] # right part
+ incr2[:, -shift:, -shift:] = increment[:, :shift, :shift] # bottom right corner
+ if 'plot' in args:
+ pdata = incr2[0,:,:] # select lowest level
+ pdata = pdata - pdata.min()
+ c = next(ax._get_lines.prop_cycler)['color']
+ ax.contour(pdata, levels=[2.,], linewidths=1, colors=[c], zorder=3)
+ ax.contour(pdata, levels=[1.,], linewidths=0.5, colors=[c], zorder=3)
+
+ ds.variables['T'][0,...] += incr2
+ ds.variables['THM'][0,...] += incr2
print(wrfin, 'wbubble inserted.')
+
+if 'plot' in args:
+ ax.set_aspect(1)
+ ax.plot([], [], color='k', lw=1, label='+2 K perturbation')
+ ax.plot([], [], color='k', lw=0.5, label='+1 K perturbation')
+ ax.legend()
+
+ fout = '/home/fs71386/lkugler/data/analysis/'+exp.expname+'/shiftedbubbles.png'
+ os.makedirs(os.path.dirname(fout), exist_ok=True)
+ fig.savefig(fout, dpi=200)
+ print(fout, 'saved.')
+
--
GitLab