-
Notifications
You must be signed in to change notification settings - Fork 2
/
initialisation.f90
513 lines (451 loc) · 17.4 KB
/
initialisation.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
!>@author
!>Paul Connolly, The University of Manchester
!>@brief
!>initialisation for the thermal cloud model and solves the
!>hydrostatic equation for pressure:
!>\f$ \frac{\partial P}{\partial z} = - \frac{P}{R_aT}g \f$
module initialisation
use numerics_type
! use variables
private
public :: calc_profile_2d, allocate_arrays
contains
!>@author
!>Paul J. Connolly, The University of Manchester
!>@brief
!>allocate arrays for q_type and q_init
!>@param[in] nq number of q fields
!>@param[in] n_levels number of levels for reading in sounding
!>@param[inout] q_type: integer array
!>@param[inout] q_init: logical array
!>@param[inout] q_read: real array
subroutine allocate_arrays(nq,n_levels,q_type,q_init,q_read)
use numerics_type
implicit none
integer(i4b), intent(in) :: nq, n_levels
integer(i4b), dimension(:), allocatable, intent(inout) :: q_type
logical, dimension(:), allocatable, intent(inout) :: q_init
real(wp), dimension(:,:), allocatable, intent(inout) :: q_read
! local variables:
integer(i4b) :: AllocateStatus
! allocate arrays
allocate( q_type(1:nq), STAT = AllocateStatus)
if (AllocateStatus /= 0) STOP "*** Not enough memory ***"
allocate( q_init(1:nq), STAT = AllocateStatus)
if (AllocateStatus /= 0) STOP "*** Not enough memory ***"
allocate( q_read(1:nq,1:n_levels), STAT = AllocateStatus)
if (AllocateStatus /= 0) STOP "*** Not enough memory ***"
end subroutine allocate_arrays
!>@author
!>Paul J. Connolly, The University of Manchester
!>@brief
!>interpolates the sounding to the grid
!>@param[in] nq number of q fields
!>@param[in] nprec number of precipitation arrays
!>@param[in] n_levels number levels for sounding
!>@param[in] psurf surface pressure
!>@param[in] tsurf surface temperature
!>@param[in] t_cbase cloud base temperature
!>@param[in] t_ctop cloud top temperature
!>@param[in] adiabatic_prof: flag if we want an adiabatic profile
!>@param[in] adiabatic_frac: fraction of adiabatic liquid water content in cloud
!>@param[in] q_type flag for type of q-field
!>@param[in] q_init flag for whether q-field initialised
!>@param[in] z_read vertical levels for sounding
!>@param[in] theta_read theta on vertical levels for sounding
!>@param[in] q_read q fields on vertical levels for sounding
!>@param[in] ip number of vertical levels of grid
!>@param[in] kp number of vertical levels of grid
!>@param[in] o_halo number of extra grid levels required for advection
!>@param[in] dx horizontal resolution of grid
!>@param[in] dz vertical resolution of grid
!>@param[in] dx2 horizontal resolution of grid
!>@param[in] dz2 vertical resolution of grid
!>@param[inout] q, qold,
!>@param[in] iqv,iqc,inc,iqi,ini
!>@param[inout] precip, theta, th_old,
!> pressure, x,xn,z,zn, temperature, rho,u,w, delsq, vis
!>@param[in] drop_num_init: flag to initialise number of drops where liquid water>0
!>@param[in] number conc of drops #/kg
!>@param[in] ice_init: flag to initialise ice crystals in model
!>@param[in] number conc of ice crystals #/kg
!>@param[in] mass of a single ice crystal kg.
!>@param[inout] zbase, ztop
!>@param[in] microphysics_flag: flag for kind of microphysics used
!>@param[in] above_cloud: flag for assumption above cloud
subroutine calc_profile_2d(nq,nprec,n_levels,psurf,tsurf,t_cbase, &
t_ctop, adiabatic_prof, adiabatic_frac,q_type,q_init, &
z_read,theta_read,q_read, &
ip,kp,o_halo,dx,dz,dx2,dz2,q,qold, &
iqv,iqc,inc,iqi,ini,&
precip,theta,th_old, p,x,xn,z,zn,t,rho,u,w,&
delsq, vis, &
drop_num_init, num_drop, &
ice_init, num_ice, mass_ice, &
zbase,ztop, &
microphysics_flag, above_cloud)
use numerics_type
use numerics, only : find_pos, poly_int, vode_integrate, zeroin
use constants
use variables, only : theta_surf, theta_q_sat, w_cb, t1old, p111
!use micro_module, only : svp_liq
implicit none
! inputs
integer(i4b), intent(in) :: n_levels, nq,o_halo, nprec
real(wp), dimension(n_levels), intent(in) :: z_read, theta_read
real(wp), dimension(nq,n_levels), intent(in) :: q_read
integer(i4b), dimension(nq), intent(in) :: q_type
logical, dimension(nq), intent(in) :: q_init
integer(i4b), intent(in) :: ip, kp, microphysics_flag, above_cloud
real(wp), intent(in) :: dx, dz, psurf, tsurf, t_cbase, t_ctop
logical, intent(in) :: adiabatic_prof, ice_init, drop_num_init
real(wp), intent(in) :: adiabatic_frac
real(wp), intent(in) :: num_drop, num_ice, mass_ice
! inouts
real(wp), intent(inout) :: zbase, ztop
real(wp), dimension(:,:), allocatable, intent(inout) :: theta, th_old, &
p, t, rho,u, w,delsq, &
vis
real(wp), dimension(:), allocatable, intent(inout) :: x, z,xn,zn, dx2, dz2
real(wp), dimension(:,:,:), allocatable, intent(inout) :: q, qold, precip
integer(i4b), intent(in) :: iqv, iqc, inc, iqi, ini
! local variables:
integer(i4b) :: i,j, iloc, AllocateStatus, istore,istore2
real(wp) :: var, dummy
! variables for odesolver:
real(wp), dimension(1) :: z1, p11,p22
real(wp) :: htry,hmin,eps2,p1,p2,p_ctop,z11,z22,theta1
! allocate arrays
allocate( precip(1:kp,1:ip,1:nprec), STAT = AllocateStatus)
if (AllocateStatus /= 0) STOP "*** Not enough memory ***"
allocate( q(-o_halo+1:kp+o_halo,-o_halo+1:ip+o_halo,1:nq), STAT = AllocateStatus)
if (AllocateStatus /= 0) STOP "*** Not enough memory ***"
allocate( qold(-o_halo+1:kp+o_halo,-o_halo+1:ip+o_halo,1:nq), STAT = AllocateStatus)
if (AllocateStatus /= 0) STOP "*** Not enough memory ***"
allocate( theta(-o_halo+1:kp+o_halo,-o_halo+1:ip+o_halo), STAT = AllocateStatus)
if (AllocateStatus /= 0) STOP "*** Not enough memory ***"
allocate( th_old(-o_halo+1:kp+o_halo,-o_halo+1:ip+o_halo), STAT = AllocateStatus)
if (AllocateStatus /= 0) STOP "*** Not enough memory ***"
allocate( p(-o_halo+1:kp+o_halo,-o_halo+1:ip+o_halo), STAT = AllocateStatus)
if (AllocateStatus /= 0) STOP "*** Not enough memory ***"
allocate( x(-o_halo+1:ip+o_halo), STAT = AllocateStatus)
if (AllocateStatus /= 0) STOP "*** Not enough memory ***"
allocate( z(-o_halo+1:kp+o_halo), STAT = AllocateStatus)
if (AllocateStatus /= 0) STOP "*** Not enough memory ***"
allocate( dx2(-o_halo+1:ip+o_halo), STAT = AllocateStatus)
if (AllocateStatus /= 0) STOP "*** Not enough memory ***"
allocate( dz2(-o_halo+1:kp+o_halo), STAT = AllocateStatus)
if (AllocateStatus /= 0) STOP "*** Not enough memory ***"
allocate( xn(-o_halo+1:ip+o_halo), STAT = AllocateStatus)
if (AllocateStatus /= 0) STOP "*** Not enough memory ***"
allocate( zn(-o_halo+1:kp+o_halo), STAT = AllocateStatus)
if (AllocateStatus /= 0) STOP "*** Not enough memory ***"
allocate( t(-o_halo+1:kp+o_halo,-o_halo+1:ip+o_halo), STAT = AllocateStatus)
if (AllocateStatus /= 0) STOP "*** Not enough memory ***"
allocate( rho(-o_halo+1:kp+o_halo,-o_halo+1:ip+o_halo), STAT = AllocateStatus)
if (AllocateStatus /= 0) STOP "*** Not enough memory ***"
allocate( u(-o_halo+1:kp+o_halo,-o_halo+1:ip+o_halo), STAT = AllocateStatus)
if (AllocateStatus /= 0) STOP "*** Not enough memory ***"
allocate( w(-o_halo+1:kp+o_halo,-o_halo+1:ip+o_halo), STAT = AllocateStatus)
if (AllocateStatus /= 0) STOP "*** Not enough memory ***"
allocate( delsq(1:kp,1:ip), STAT = AllocateStatus)
if (AllocateStatus /= 0) STOP "*** Not enough memory ***"
allocate( vis(-o_halo+1:kp+o_halo,-o_halo+1:ip+o_halo), STAT = AllocateStatus)
if (AllocateStatus /= 0) STOP "*** Not enough memory ***"
precip=0._wp
! set up vertical level array
z=dz*(/(i,i=-o_halo,kp+o_halo-1)/) +dz
zn=z-0.5_wp*dz
! set up horizontal level array
x=dx*(/(i,i=-o_halo,ip+o_halo-1)/)!-0.5_wp*dx
x=x-0.5_wp*dx
xn=x-0.5_wp*dx
dx2=dx
dz2=dz
q=0._wp
if (adiabatic_prof) then
! ++++ Dry adiabat
! calculate the dry adiabat:
theta_surf=tsurf*(1.e5_wp/psurf)**(ra/cp)
p1=psurf
z1=0._wp
p2=1.e5_wp*(t_cbase/theta_surf)**(cp/ra)
htry=p2-psurf
eps2=1.e-5_wp
call vode_integrate(z1,p1,p2,eps2,htry,hmin,hydrostatic1)
p1=p2 ! p1 is the cloud-base pressure, z1 is the cb height
zbase=z1(1)-dz*0.5_wp
! integrate going downwards - dry adiabatic layer
p(1,:)=psurf
do i=1,-o_halo+2,-1
z11=z(i)
z22=z(i-1)
p11=p(i,1)
htry=-dz
hmin=1.e-2_wp
! finds the pressure on different levels:
call vode_integrate(p11,z11,z22,eps2,htry,hmin,hydrostatic1a)
p(i-1,:)=p11(1)
enddo
! integrate going upwards - dry adiabatic layer
p(1,:)=psurf
do i=1,kp+o_halo-1
z11=z(i)
z22=z(i+1)
if(z22.gt.z1(1)) exit ! gt cloud-base height
p11=p(i,1)
htry=dz
hmin=1.e-2_wp
! finds the pressure on different levels:
call vode_integrate(p11,z11,z22,eps2,htry,hmin,hydrostatic1a)
p(i+1,:)=p11(1)
enddo
istore=i-1
! adiabatic temperature
t(-o_halo+1:istore,:)=theta_surf*(p(-o_halo+1:istore,:)/1.e5_wp)**(ra/cp)
! adiabatic vapour mixing ratio - at all levels below CB:
q(-o_halo+1:istore,:,iqv)=eps1*svp_liq(t_cbase)/(p1-svp_liq(t_cbase))
! ---- Dry adiabat done
! ++++ Moist adiabat
! now calculate the moist adiabat
w_cb=eps1*svp_liq(t_cbase)/(p1-svp_liq(t_cbase))
theta_q_sat=t_cbase*(1.e5_wp/p1)**(ra/cp)*exp(lv*w_cb/cp/t_cbase)
! theta_q_sat is conserved. Use it to calculate the new temperature
t1old=t_ctop
! calculates the pressure at t1old, conserving theta_q_sat:
p_ctop=zeroin(p1,3000._wp,calc_theta_q2,1.e-5_wp)
! stop
! integrate going upwards - moist adiabatic layer
t1old=t_ctop
istore=istore+1
! now find the temperature
p111=p(istore,1)
t(istore,:)=theta_surf*( &
(p(istore,:)-dz*p(istore,:)/ra/t_ctop)/1.e5_wp)**(ra/cp) ! a temperature colder than
! next level
! calculates the temperature at p111, conserving theta_q_sat
! this is at cloud-base (i.e. istore)
t(istore,:)=zeroin(1.01_wp*t(istore,1),t_ctop,calc_theta_q,1.e-5_wp)
t1old=t(istore,1)
do i=istore,kp+o_halo-1
z22=z(i+1)
if (i.eq.istore) then
z11=z1(1)
p11=p1
else
z11=z(i)
p11=p(i,1)
endif
htry=dz
hmin=1.e-2_wp
!print *,t1old,p11,z11,z22
call vode_integrate(p11,z11,z22,eps2,htry,hmin,hydrostatic2a)
p(i+1,:)=p11(1)
t(i+1,:)=theta_surf*(p(i+1,:)/1.e5_wp)**(ra/cp)
t1old=t(i,1)
p111=p(i+1,1)
t(i+1,:)=zeroin(t(i+1,1),t1old*1.01_wp,calc_theta_q,1.e-5_wp)
if(t(i+1,1).lt.t_ctop) exit
enddo
istore2=i-1
ztop=zn(istore2) !(z(istore2-1)+z(istore2))*0.5_wp
do i=istore,istore2
q(i,:,iqv)=eps1*svp_liq(t(i,1))/ &
(p(i,1)-svp_liq(t(i,1)))
q(i,:,iqc)=adiabatic_frac* &
max(q(1,1,iqv)-eps1*svp_liq(t(i,1))/(p(i,1)-svp_liq(t(i,1))) ,0._wp)
if(drop_num_init .and. &
((microphysics_flag .eq. 2) .or. (microphysics_flag .eq. 3))) then
q(i,:,inc) = num_drop
endif
enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! integrate going upwards - above cloud-top
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
select case (above_cloud)
case (0) ! dry above, conserving theta
theta1=t_ctop*(1.e5_wp/p_ctop)**(ra/cp)
istore2=istore2+1
do i=istore2,kp+o_halo-1
z11=z(i)
z22=z(i+1)
p11=p(i,1)
htry=dz
hmin=1.e-2_wp
call vode_integrate(p11,z11,z22,eps2,htry,hmin,hydrostatic1a)
p(i+1,:)=p11(1)
enddo
t(istore2:kp+o_halo,:)=theta1*(p(istore2:kp+o_halo,:)/1.e5_wp)**(ra/cp)
case (1) ! rh=0.95 above, conserving theta_q
do i=istore2,kp+o_halo-1
z11=z(i)
z22=z(i+1)
p11=p(i,1)
htry=dz
hmin=1.e-2_wp
call vode_integrate(p11,z11,z22,eps2,htry,hmin,hydrostatic2a)
p(i+1,:)=p11(1)
t(i+1,:)=theta_surf*(p(i+1,:)/1.e5_wp)**(ra/cp)
t1old=t(i,1)
p111=p(i+1,1)
t(i+1,:)=zeroin(t(i+1,1),t1old*1.01_wp,calc_theta_q,1.e-5_wp)
q(i,:,iqv)=0.95_wp*eps1*svp_liq(t(i,1))/ &
(p(i,1)-svp_liq(t(i,1)))
enddo
case default
print *,'error above_cloud: ',above_cloud
stop
end select
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! initialise ice crystals
if(ice_init .and. (microphysics_flag .eq. 1)) then
where(t(istore:istore2,:).lt.ttr)
q(istore:istore2,:,iqi)=num_ice*mass_ice
q(istore:istore2,:,ini)=num_ice
end where
endif
else
! use linear interpolation to put sounding on grid:
do i=-o_halo+1,kp+o_halo
iloc=find_pos(z_read(1:n_levels),z(i))
iloc=min(n_levels-1,iloc)
iloc=max(1,iloc)
! linear interp theta
call poly_int(z_read(iloc:iloc+1), theta_read(iloc:iloc+1), &
min(z(i),z_read(n_levels)), var,dummy)
theta(i,:)=var
! q-fields
do j=1,nq
if(q_init(j)) then
! linear interp q fields
call poly_int(z_read(iloc:iloc+1), q_read(j,iloc:iloc+1), &
min(z(i),z_read(n_levels)), var, dummy)
q(i,:,j)=var
else
if(q_type(j).eq.2) then
q(i,:,j) = .0_wp
else
q(i,:,j) = 0.0_wp
endif
endif
enddo
enddo
p(1,:)=psurf
do i=1,kp+o_halo
! solve the hydrostatic equation to get pressure and temperature
t(i,:)=theta(i,:)*(p(i,:)/psurf)**(ra/cp)
if(i.lt.(kp+o_halo)) then
p(i+1,:)=p(i,:)-dz*p(i,:)/(ra*t(i,:))*grav
endif
enddo
do i=0,-o_halo+1,-1
! solve the hydrostatic equation to get pressure and temperature
p(i,:)=p(i+1,:)+dz*p(i+1,:)/(ra*t(i+1,:))*grav
t(i,:)=theta(i,:)*(p(i,:)/psurf)**(ra/cp)
enddo
endif
u(:,:)=0._wp
w(:,:)=0._wp
theta=t*(1.e5_wp/p)**(ra/cp)
rho=p/(ra*t)
end subroutine calc_profile_2d
subroutine hydrostatic1(p,z,dzdp)
use numerics_type
use constants
use variables, only : theta_surf
implicit none
real(wp), intent(in) :: p
real(wp), dimension(:), intent(in) :: z
real(wp), dimension(:), intent(out) :: dzdp
real(wp) :: t
t=theta_surf*(p/1.e5_wp)**(ra/cp)
dzdp(1)=-(ra*t) / (grav*p)
end subroutine hydrostatic1
subroutine hydrostatic1a(z,p,dpdz)
use numerics_type
use constants
use variables, only : theta_surf
implicit none
real(wp), intent(in) :: z
real(wp), dimension(:), intent(in) :: p
real(wp), dimension(:), intent(out) :: dpdz
real(wp) :: t
t=theta_surf*(p(1)/1.e5_wp)**(ra/cp)
dpdz(1)=-(grav*p(1)) / (ra*t)
end subroutine hydrostatic1a
subroutine hydrostatic2(p,z,dzdp)
use numerics_type
use numerics, only : zeroin
use constants
use variables, only : theta_surf,theta_q_sat, w_cb, t1old, p111
implicit none
real(wp), intent(in) :: p
real(wp), dimension(:), intent(in) :: z
real(wp), dimension(:), intent(out) :: dzdp
real(wp) :: t
p111=p
t=theta_surf*(p111/1.e5_wp)**(ra/cp)
t=zeroin(t,t1old*1.01_wp,calc_theta_q,1.e-5_wp)
! print *,'hi',t,calc_theta_q(t)
! find the temperature by iteration
dzdp(1)=-(ra*t) / (grav*p)
end subroutine hydrostatic2
subroutine hydrostatic2a(z,p,dpdz)
use numerics_type
use numerics, only : zeroin
use constants
use variables, only : theta_surf,theta_q_sat, w_cb, t1old, p111
implicit none
real(wp), intent(in) :: z
real(wp), dimension(:), intent(in) :: p
real(wp), dimension(:), intent(out) :: dpdz
real(wp) :: t
p111=p(1)
t=theta_surf*(p111/1.e5_wp)**(ra/cp)
t=zeroin(t,t1old*1.01_wp,calc_theta_q,1.e-5_wp)
! print *,'hi',t,calc_theta_q(t)
! find the temperature by iteration
dpdz(1)=-(grav*p(1))/(ra*t)
end subroutine hydrostatic2a
function calc_theta_q(t111)
use numerics_type
use constants
use variables, only : theta_q_sat, p111
implicit none
real(wp), intent(in) :: t111
real(wp) :: calc_theta_q
real(wp) :: ws
ws=eps1*svp_liq(t111)/(p111-svp_liq(t111))
calc_theta_q=t111*(1.e5_wp/p111)**(ra/cp)*exp(lv*ws/cp/t111)-theta_q_sat
end function calc_theta_q
function calc_theta_q2(p)
use numerics_type
use constants
use variables, only : theta_q_sat, t1old
implicit none
real(wp), intent(in) :: p
real(wp) :: calc_theta_q2
real(wp) :: ws
ws=eps1*svp_liq(t1old)/(p-svp_liq(t1old))
calc_theta_q2=t1old*(1e5_wp/p)**(ra/cp)*exp(lv*ws/cp/t1old)-theta_q_sat
end function calc_theta_q2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! saturation vapour pressure over liquid !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!>@author
!>Paul J. Connolly, The University of Manchester
!>@brief
!>calculates the saturation vapour pressure over liquid water according to buck fit
!>@param[in] t: temperature
!>@return svp_liq: saturation vapour pressure over liquid water
function svp_liq(t)
use numerics_type
use constants, only : ttr
implicit none
real(wp), intent(in) :: t
real(wp) :: svp_liq
svp_liq = 100._wp*6.1121_wp* &
exp((18.678_wp - (t-ttr)/ 234.5_wp)* &
(t-ttr)/(257.14_wp + (t-ttr)))
end function svp_liq
end module initialisation