Skip to content

Commit 42d085f

Browse files
adamrheradamrher
committed
changes to clubb_mf to operate on descending arrays
1 parent f20fbf1 commit 42d085f

2 files changed

Lines changed: 69 additions & 130 deletions

File tree

‎src/physics/cam/clubb_intr.F90‎

Lines changed: 9 additions & 71 deletions
Original file line numberDiff line numberDiff line change
@@ -3473,30 +3473,10 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, &
34733473
rtm_zm_in = zt2zm_api( nzm_clubb, nzt_clubb, ncol, gr, rtm(:ncol,:) )
34743474
thlm_zm_in = zt2zm_api( nzm_clubb, nzt_clubb, ncol, gr, thlm(:ncol,:) )
34753475

3476-
!--------------------------------------- integrate_mf call and flip ---------------------------------------
3477-
! integrate_mf assumes an ascending grid, which is the opposide of the cam grid that
3478-
! clubb_intr now mainly uses, so we need to flip the fields before calling integrate_mf
3479-
!
3480-
! Ideally, integrate_mf would operate in descending mode, then we could remove the flipping.
3476+
!--------------------------------------- integrate_mf call ---------------------------------------
3477+
! integrate_mf expects arguments of individual columns.
34813478
! If the column loop gets pushed into it, we can also avoid the array slicing.
34823479

3483-
dz_g = dz_g(:,nzt_clubb:1:-1)
3484-
p_in_Pa = p_in_Pa(:,nzt_clubb:1:-1)
3485-
invrs_exner_zt = invrs_exner_zt(:,nzt_clubb:1:-1)
3486-
um = um(:,nzt_clubb:1:-1)
3487-
vm = vm(:,nzt_clubb:1:-1)
3488-
thlm = thlm(:,nzt_clubb:1:-1)
3489-
rtm = rtm(:,nzt_clubb:1:-1)
3490-
3491-
thv_ds_zt = thv_ds_zt(:,nzt_clubb:1:-1)
3492-
3493-
! Flip zm inputs
3494-
zi_g = zi_g(:,nzm_clubb:1:-1)
3495-
p_in_Pa_zm = p_in_Pa_zm(:,nzm_clubb:1:-1)
3496-
invrs_exner_zm = invrs_exner_zm(:,nzm_clubb:1:-1)
3497-
thlm_zm_in = thlm_zm_in(:,nzm_clubb:1:-1)
3498-
rtm_zm_in = rtm_zm_in(:,nzm_clubb:1:-1)
3499-
35003480
do i = 1, ncol
35013481
call integrate_mf( nzm_clubb, nzt_clubb, dz_g(i,:), zi_g(i,:), p_in_Pa_zm(i,:), invrs_exner_zm(i,:), & ! input
35023482
p_in_Pa(i,:), invrs_exner_zt(i,:), & ! input
@@ -3510,56 +3490,14 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, &
35103490
mf_dry_u(i,:), mf_moist_u(i,:), & ! output - plume diagnostics
35113491
mf_dry_v(i,:), mf_moist_v(i,:), & ! output - plume diagnostics
35123492
mf_moist_qc(i,:), & ! output - plume diagnostics
3513-
s_ae(i,:), s_aw(i,:), & ! output - plume diagnostics
3514-
s_awthl(i,:), s_awqt(i,:), & ! output - plume diagnostics
3515-
s_awql(i,:), s_awqi(i,:), & ! output - plume diagnostics
3516-
s_awu(i,:), s_awv(i,:), & ! output - plume diagnostics
3517-
mf_thlflx(i,:), mf_qtflx(i,:) ) ! output - variables needed for solver
3493+
s_ae(i,:), s_aw(i,:), & ! output - plume diagnostics
3494+
s_awthl(i,:), s_awqt(i,:), & ! output - plume diagnostics
3495+
s_awql(i,:), s_awqi(i,:), & ! output - plume diagnostics
3496+
s_awu(i,:), s_awv(i,:), & ! output - plume diagnostics
3497+
mf_thlflx(i,:), mf_qtflx(i,:) ) ! output - variables needed for solver
35183498
end do
3519-
3520-
! Flip zt inputs back
3521-
dz_g = dz_g(:,nzt_clubb:1:-1)
3522-
p_in_Pa = p_in_Pa(:,nzt_clubb:1:-1)
3523-
invrs_exner_zt = invrs_exner_zt(:,nzt_clubb:1:-1)
3524-
um = um(:,nzt_clubb:1:-1)
3525-
vm = vm(:,nzt_clubb:1:-1)
3526-
thlm = thlm(:,nzt_clubb:1:-1)
3527-
rtm = rtm(:,nzt_clubb:1:-1)
3528-
3529-
thv_ds_zt = thv_ds_zt(:,nzt_clubb:1:-1)
3530-
3531-
! Flip zm inputs back
3532-
zi_g = zi_g(:,nzm_clubb:1:-1)
3533-
p_in_Pa_zm = p_in_Pa_zm(:,nzm_clubb:1:-1)
3534-
invrs_exner_zm = invrs_exner_zm(:,nzm_clubb:1:-1)
3535-
thlm_zm_in = thlm_zm_in(:,nzm_clubb:1:-1)
3536-
rtm_zm_in = rtm_zm_in(:,nzm_clubb:1:-1)
3537-
3538-
! Flip clubb_mf output, since it assumes an ascending grid currently
3539-
mf_dry_a = mf_dry_a(:,nzm_clubb:1:-1)
3540-
mf_moist_a = mf_moist_a(:,nzm_clubb:1:-1)
3541-
mf_dry_w = mf_dry_w(:,nzm_clubb:1:-1)
3542-
mf_moist_w = mf_moist_w(:,nzm_clubb:1:-1)
3543-
mf_dry_qt = mf_dry_qt(:,nzm_clubb:1:-1)
3544-
mf_moist_qt = mf_moist_qt(:,nzm_clubb:1:-1)
3545-
mf_dry_thl = mf_dry_thl(:,nzm_clubb:1:-1)
3546-
mf_moist_thl = mf_moist_thl(:,nzm_clubb:1:-1)
3547-
mf_dry_u = mf_dry_u(:,nzm_clubb:1:-1)
3548-
mf_moist_u = mf_moist_u(:,nzm_clubb:1:-1)
3549-
mf_dry_v = mf_dry_v(:,nzm_clubb:1:-1)
3550-
mf_moist_v = mf_moist_v(:,nzm_clubb:1:-1)
3551-
mf_moist_qc = mf_moist_qc(:,nzm_clubb:1:-1)
3552-
s_ae = s_ae(:,nzm_clubb:1:-1)
3553-
s_aw = s_aw(:,nzm_clubb:1:-1)
3554-
s_awthl = s_awthl(:,nzm_clubb:1:-1)
3555-
s_awqt = s_awqt(:,nzm_clubb:1:-1)
3556-
s_awql = s_awql(:,nzm_clubb:1:-1)
3557-
s_awqi = s_awqi(:,nzm_clubb:1:-1)
3558-
s_awu = s_awu(:,nzm_clubb:1:-1)
3559-
s_awv = s_awv(:,nzm_clubb:1:-1)
3560-
mf_thlflx = mf_thlflx(:,nzm_clubb:1:-1)
3561-
mf_qtflx = mf_qtflx(:,nzm_clubb:1:-1)
3562-
!--------------------------------------- END integrate_mf call and flip ---------------------------------------
3499+
3500+
!--------------------------------------- END integrate_mf call ---------------------------------------
35633501

35643502
! pass MF turbulent advection term as CLUBB explicit forcing term
35653503
do k = 1, nzt_clubb

‎src/physics/cam/clubb_mf.F90‎

Lines changed: 60 additions & 59 deletions
Original file line numberDiff line numberDiff line change
@@ -281,7 +281,7 @@ subroutine integrate_mf( nzm, nzt, dzt, zm, p_zm, iexner_zm,
281281
zcb = zcb_unset
282282

283283
pblh = max(pblh,pblhmin)
284-
wthv = wthl+zvir*thv(1)*wqt
284+
wthv = wthl+zvir*thv(nzt)*wqt
285285

286286
! if surface buoyancy is positive then do mass-flux
287287
if ( wthv > 0._r8 ) then
@@ -290,28 +290,29 @@ subroutine integrate_mf( nzm, nzt, dzt, zm, p_zm, iexner_zm,
290290
! overide stochastic entrainment with fixent
291291
ent(:,:) = fixent
292292
else
293-
293+
294294
! get entrainment coefficient, dz/L0
295295
do i=1,clubb_mf_nup
296-
do k=1,nzt
296+
!do k=1,nzt
297+
do k=nzt,1,-1
297298
entf(k,i) = dzt(k) / clubb_mf_L0
298299
enddo
299300
enddo
300-
301+
301302
! get poisson, P(dz/L0)
302-
call poisson( nzt, clubb_mf_nup, entf, enti, u(1:4))
303-
303+
call poisson( nzt, clubb_mf_nup, entf, enti, u(nzt:nzt-3:-1))
304+
304305
! get entrainment, ent=ent0/dz*P(dz/L0)
305306
do i=1,clubb_mf_nup
306-
do k=1,nzt
307+
do k=nzt,1,-1
307308
ent(k,i) = real( enti(k,i))*clubb_mf_ent0/dzt(k)
308309
enddo
309310
enddo
310-
311+
311312
end if
312313

313314
! get surface conditions
314-
wstar = max( wstarmin, (gravit/thv(1)*wthv*pblh)**(1._r8/3._r8) )
315+
wstar = max( wstarmin, (gravit/thv(nzt)*wthv*pblh)**(1._r8/3._r8) )
315316
qstar = wqt / wstar
316317
thvstar = wthv / wstar
317318

@@ -327,105 +328,105 @@ subroutine integrate_mf( nzm, nzt, dzt, zm, p_zm, iexner_zm,
327328
wlv = wmin + (wmax-wmin) / (real(clubb_mf_nup,r8)) * (real(i-1, r8))
328329
wtv = wmin + (wmax-wmin) / (real(clubb_mf_nup,r8)) * real(i,r8)
329330

330-
upw(1,i) = 0.5_r8 * (wlv+wtv)
331-
upa(1,i) = 0.5_r8 * erf( wtv/(sqrt(2._r8)*sigmaw) ) &
331+
upw(nzm,i) = 0.5_r8 * (wlv+wtv)
332+
upa(nzm,i) = 0.5_r8 * erf( wtv/(sqrt(2._r8)*sigmaw) ) &
332333
- 0.5_r8 * erf( wlv/(sqrt(2._r8)*sigmaw) )
333334

334-
upu(1,i) = u(1)
335-
upv(1,i) = v(1)
335+
upu(nzm,i) = u(nzt)
336+
upv(nzm,i) = v(nzt)
336337

337-
upqt(1,i) = qt(1) + cwqt * upw(1,i) * sigmaqt/sigmaw
338-
upthv(1,i) = thv(1) + cwthv * upw(1,i) * sigmathv/sigmaw
339-
upthl(1,i) = upthv(1,i) / (1._r8+zvir*upqt(1,i))
338+
upqt(nzm,i) = qt(nzt) + cwqt * upw(nzm,i) * sigmaqt/sigmaw
339+
upthv(nzm,i) = thv(nzt) + cwthv * upw(nzm,i) * sigmathv/sigmaw
340+
upthl(nzm,i) = upthv(nzm,i) / (1._r8+zvir*upqt(nzm,i))
340341

341342
! get cloud, lowest momentum level
342343
if (do_condensation) then
343-
call condensation_mf(upqt(1,i), upthl(1,i), p_zm(1), iexner_zm(1), &
344+
call condensation_mf(upqt(nzm,i), upthl(nzm,i), p_zm(nzm), iexner_zm(nzm), &
344345
thvn, qcn, thn, qln, qin, qsn, lmixn)
345-
upthv(1,i) = thvn
346-
upqc(1,i) = qcn
347-
upql(1,i) = qln
348-
upqi(1,i) = qin
349-
upqs(1,i) = qsn
350-
if (qcn > 0._r8) zcb(i) = zm(1)
346+
upthv(nzm,i) = thvn
347+
upqc(nzm,i) = qcn
348+
upql(nzm,i) = qln
349+
upqi(nzm,i) = qin
350+
upqs(nzm,i) = qsn
351+
if (qcn > 0._r8) zcb(i) = zm(nzm)
351352
else
352353
! assume no cldliq
353-
upqc(1,i) = 0._r8
354+
upqc(nzm,i) = 0._r8
354355
end if
355356

356357
enddo
357358

358359
! get updraft properties
359360
do i=1,clubb_mf_nup
360-
do k=1,nzm-1
361-
361+
do k=nzm,2,-1
362+
362363
! get microphysics, autoconversion
363364
if (do_precip .and. upqc(k,i) > 0._r8) then
364-
call precip_mf(upqs(k,i),upqt(k,i),upw(k,i),dzt(k),zm(k+1)-zcb(i),supqt)
365+
call precip_mf(upqs(k,i),upqt(k,i),upw(k,i),dzt(k-1),zm(k-1)-zcb(i),supqt)
365366

366-
supthl = -1._r8*lmixn*supqt*iexner_zt(k)/cpair
367+
supthl = -1._r8*lmixn*supqt*iexner_zt(k-1)/cpair
367368
else
368369
supqt = 0._r8
369370
supthl = 0._r8
370371
end if
371372

372373
! integrate updraft
373-
entexp = exp(-ent(k,i)*dzt(k))
374-
entexpu = exp(-ent(k,i)*dzt(k)/3._r8)
374+
entexp = exp(-ent(k-1,i)*dzt(k-1))
375+
entexpu = exp(-ent(k-1,i)*dzt(k-1)/3._r8)
375376

376-
qtn = qt(k) *(1._r8-entexp ) + upqt (k,i)*entexp + supqt
377-
thln = thl(k)*(1._r8-entexp ) + upthl(k,i)*entexp + supthl
378-
un = u(k) *(1._r8-entexpu) + upu (k,i)*entexpu
379-
vn = v(k) *(1._r8-entexpu) + upv (k,i)*entexpu
377+
qtn = qt(k-1) *(1._r8-entexp ) + upqt (k,i)*entexp + supqt
378+
thln = thl(k-1)*(1._r8-entexp ) + upthl(k,i)*entexp + supthl
379+
un = u(k-1) *(1._r8-entexpu) + upu (k,i)*entexpu
380+
vn = v(k-1) *(1._r8-entexpu) + upv (k,i)*entexpu
380381

381382
! get cloud, momentum levels
382383
if (do_condensation) then
383-
call condensation_mf(qtn, thln, p_zm(k+1), iexner_zm(k+1), &
384+
call condensation_mf(qtn, thln, p_zm(k-1), iexner_zm(k-1), &
384385
thvn, qcn, thn, qln, qin, qsn, lmixn)
385-
if (zcb(i).eq.zcb_unset .and. qcn > 0._r8) zcb(i) = zm(k+1)
386+
if (zcb(i).eq.zcb_unset .and. qcn > 0._r8) zcb(i) = zm(k-1)
386387
else
387388
thvn = thln*(1._r8+zvir*qtn)
388389
end if
389390

390391
! get buoyancy
391-
B=gravit*(0.5_r8*(thvn + upthv(k,i))/thv(k)-1._r8)
392+
B=gravit*(0.5_r8*(thvn + upthv(k,i))/thv(k-1)-1._r8)
392393
if (debug) then
393394
if ( masterproc ) then
394395
write(iulog,*) "B(k,i), k, i ", B, k, i
395396
end if
396397
end if
397398

398399
! get wn^2
399-
wp = wb*ent(k,i)
400+
wp = wb*ent(k-1,i)
400401
if (wp==0._r8) then
401-
wn2 = upw(k,i)**2._r8+2._r8*wa*B*dzt(k)
402+
wn2 = upw(k,i)**2._r8+2._r8*wa*B*dzt(k-1)
402403
else
403-
entw = exp(-2._r8*wp*dzt(k))
404-
wn2 = entw*upw(k,i)**2._r8+wa*B/(wb*ent(k,i))*(1._r8-entw)
404+
entw = exp(-2._r8*wp*dzt(k-1))
405+
wn2 = entw*upw(k,i)**2._r8+wa*B/(wb*ent(k-1,i))*(1._r8-entw)
405406
end if
406407

407408
if (wn2>0._r8) then
408-
upw(k+1,i) = sqrt(wn2)
409-
upthv(k+1,i) = thvn
410-
upthl(k+1,i) = thln
411-
upqt(k+1,i) = qtn
412-
upqc(k+1,i) = qcn
413-
upqs(k+1,i) = qsn
414-
upu(k+1,i) = un
415-
upv(k+1,i) = vn
416-
upa(k+1,i) = upa(k,i)
417-
upql(k+1,i) = qln
418-
upqi(k+1,i) = qin
419-
upqv(k+1,i) = qtn - qcn
409+
upw(k-1,i) = sqrt(wn2)
410+
upthv(k-1,i) = thvn
411+
upthl(k-1,i) = thln
412+
upqt(k-1,i) = qtn
413+
upqc(k-1,i) = qcn
414+
upqs(k-1,i) = qsn
415+
upu(k-1,i) = un
416+
upv(k-1,i) = vn
417+
upa(k-1,i) = upa(k,i)
418+
upql(k-1,i) = qln
419+
upqi(k-1,i) = qin
420+
upqv(k-1,i) = qtn - qcn
420421
else
421422
exit
422423
end if
423424
enddo
424425
enddo
425426

426427
! writing updraft properties for output
427-
do k=1,nzm
428-
428+
do k=nzm,1,-1
429+
429430
! first sum over all i-updrafts
430431
do i=1,clubb_mf_nup
431432
if (upqc(k,i)>0._r8) then
@@ -478,7 +479,7 @@ subroutine integrate_mf( nzm, nzt, dzt, zm, p_zm, iexner_zm,
478479

479480
enddo
480481

481-
do k=1,nzm
482+
do k=nzm,1,-1
482483
do i=1,clubb_mf_nup
483484
ae (k) = ae (k) - upa(k,i)
484485
aw (k) = aw (k) + upa(k,i)*upw(k,i)
@@ -492,7 +493,7 @@ subroutine integrate_mf( nzm, nzt, dzt, zm, p_zm, iexner_zm,
492493
enddo
493494
enddo
494495

495-
do k=1,nzm
496+
do k=nzm,1,-1
496497
thlflx(k)= awthl(k) - aw(k)*thl_zm(k)
497498
qtflx(k)= awqt(k) - aw(k)*qt_zm(k)
498499
enddo
@@ -643,7 +644,7 @@ subroutine poisson(nz,nup,lambda,poi,state)
643644
! Set seed
644645
kiss_gen = ShrKissRandGen(tmpseed)
645646

646-
do i=1,nz
647+
do i=nz,1,-1
647648
do j=1,nup
648649
call knuth(kiss_gen,lambda(i,j),poi(i,j))
649650
enddo
@@ -681,4 +682,4 @@ subroutine knuth(kiss_gen,lambda,kout)
681682
end subroutine knuth
682683

683684
end module clubb_mf
684-
685+

0 commit comments

Comments
 (0)