@@ -417,246 +417,3 @@ function external_forcing_tendency!(Yₜ, Y, p, t, ::ISDACForcing)
417
417
# total specific humidity
418
418
@. Yₜ. c. ρq_tot += Y. c. ρ * ᶜdqtdt_nudging
419
419
end
420
-
421
- # # For ERA5 reanalysis forcing data
422
- # function external_forcing_cache(Y, external_forcing::ERA5Forcing, params)
423
- # FT = Spaces.undertype(axes(Y.c))
424
- # ᶜdTdt_fluc = similar(Y.c, FT)
425
- # ᶜdqtdt_fluc = similar(Y.c, FT)
426
- # ᶜdTdt_hadv = similar(Y.c, FT)
427
- # ᶜdqtdt_hadv = similar(Y.c, FT)
428
- # ᶜdTdt_rad = similar(Y.c, FT)
429
- # ᶜT_nudge = similar(Y.c, FT)
430
- # ᶜqt_nudge = similar(Y.c, FT)
431
- # ᶜu_nudge = similar(Y.c, FT)
432
- # ᶜv_nudge = similar(Y.c, FT)
433
- # ᶜinv_τ_wind = similar(Y.c, FT)
434
- # ᶜinv_τ_scalar = similar(Y.c, FT)
435
- # ᶜls_subsidence = similar(Y.c, FT)
436
- # insolation = similar(Fields.level(Y.c.ρ, 1), FT)
437
- # cos_zenith = similar(Fields.level(Y.c.ρ, 1), FT)
438
-
439
- # (; external_forcing_file, cfsite_number) = external_forcing
440
-
441
- # NC.Dataset(external_forcing_file, "r") do ds
442
-
443
- # function setvar!(cc_field, varname, colidx, zc_gcm, zc_forcing)
444
- # parent(cc_field[colidx]) .= interp_vertical_prof(
445
- # zc_gcm,
446
- # zc_forcing,
447
- # era5_driven_profile(ds.group[cfsite_number], varname),
448
- # )
449
- # end
450
-
451
- # function setvar_subsidence!(
452
- # cc_field,
453
- # varname,
454
- # colidx,
455
- # zc_gcm,
456
- # zc_forcing,
457
- # params,
458
- # )
459
- # parent(cc_field[colidx]) .= interp_vertical_prof(
460
- # zc_gcm,
461
- # zc_forcing,
462
- # era5_driven_profile(ds.group[cfsite_number], varname) ./
463
- # .-(era5_driven_profile(ds.group[cfsite_number], "rho")) ./
464
- # CAP.grav(params),
465
- # )
466
- # end
467
-
468
- # function set_insolation!(cc_field)
469
- # parent(cc_field) .= mean(
470
- # ds.group[cfsite_number]["rsdt"][:] ./
471
- # ds.group[cfsite_number]["coszen"][:],
472
- # )
473
- # end
474
-
475
- # function set_cos_zenith!(cc_field)
476
- # parent(cc_field) .= ds.group[cfsite_number]["coszen"][1]
477
- # end
478
-
479
- # zc_forcing = era5_height(ds.group[cfsite_number])
480
- # Fields.bycolumn(axes(Y.c)) do colidx
481
-
482
- # zc_gcm = Fields.coordinate_field(Y.c).z[colidx]
483
-
484
- # setvar!(ᶜdTdt_hadv, "tntha", colidx, zc_gcm, zc_forcing)
485
- # setvar!(ᶜdqtdt_hadv, "tnhusha", colidx, zc_gcm, zc_forcing)
486
- # #setvar!(ᶜdTdt_rad, "tntr", colidx, zc_gcm, zc_forcing)
487
- # setvar_subsidence!(
488
- # ᶜls_subsidence,
489
- # "wap",
490
- # colidx,
491
- # zc_gcm,
492
- # zc_forcing,
493
- # params,
494
- # )
495
- # # GCM states, used for nudging + vertical eddy advection
496
- # setvar!(ᶜT_nudge, "ta", colidx, zc_gcm, zc_forcing)
497
- # setvar!(ᶜqt_nudge, "hus", colidx, zc_gcm, zc_forcing)
498
- # setvar!(ᶜu_nudge, "ua", colidx, zc_gcm, zc_forcing)
499
- # setvar!(ᶜv_nudge, "va", colidx, zc_gcm, zc_forcing)
500
-
501
- # # vertical eddy advection (Shen et al., 2022; eqn. 9,10)
502
- # # sum of two terms to give total tendency. First term:
503
- # setvar!(ᶜdTdt_fluc, "tntva", colidx, zc_gcm, zc_forcing)
504
- # setvar!(ᶜdqtdt_fluc, "tnhusva", colidx, zc_gcm, zc_forcing)
505
- # # second term:
506
- # eddy_vert_fluctuation!(ᶜdTdt_fluc, ᶜT_nudge, ᶜls_subsidence)
507
- # eddy_vert_fluctuation!(ᶜdqtdt_fluc, ᶜqt_nudge, ᶜls_subsidence)
508
-
509
- # set_insolation!(insolation)
510
- # set_cos_zenith!(cos_zenith)
511
-
512
- # @. ᶜinv_τ_wind[colidx] = 1 / (6 * 3600)
513
- # @. ᶜinv_τ_scalar[colidx] = compute_gcm_driven_scalar_inv_τ(zc_gcm)
514
- # end
515
- # end
516
-
517
- # return (;
518
- # ᶜdTdt_fluc,
519
- # ᶜdqtdt_fluc,
520
- # ᶜdTdt_hadv,
521
- # ᶜdqtdt_hadv,
522
- # #ᶜdTdt_rad,
523
- # ᶜT_nudge,
524
- # ᶜqt_nudge,
525
- # ᶜu_nudge,
526
- # ᶜv_nudge,
527
- # ᶜinv_τ_wind,
528
- # ᶜinv_τ_scalar,
529
- # ᶜls_subsidence,
530
- # insolation,
531
- # cos_zenith,
532
- # )
533
- # end
534
-
535
- function external_forcing_tendency! (
536
- Yₜ,
537
- Y,
538
- p,
539
- t,
540
- :: Union{ERA5Forcing, GCMForcing} ,
541
- )
542
- # horizontal advection, vertical fluctuation, nudging, subsidence (need to add),
543
- (; params) = p
544
- thermo_params = CAP. thermodynamics_params (params)
545
- (; ᶜspecific, ᶜts, ᶜh_tot) = p. precomputed
546
- (;
547
- ᶜdTdt_fluc,
548
- ᶜdqtdt_fluc,
549
- ᶜdTdt_hadv,
550
- ᶜdqtdt_hadv,
551
- # ᶜdTdt_rad,
552
- ᶜT_nudge,
553
- ᶜqt_nudge,
554
- ᶜu_nudge,
555
- ᶜv_nudge,
556
- ᶜls_subsidence,
557
- ᶜinv_τ_wind,
558
- ᶜinv_τ_scalar,
559
- ) = p. external_forcing
560
-
561
- ᶜlg = Fields. local_geometry_field (Y. c)
562
- ᶜuₕ_nudge = p. scratch. ᶜtemp_C12
563
- @. ᶜuₕ_nudge = C12 (Geometry. UVVector (ᶜu_nudge, ᶜv_nudge), ᶜlg)
564
- @. Yₜ. c. uₕ -= (Y. c. uₕ - ᶜuₕ_nudge) * ᶜinv_τ_wind
565
-
566
- # nudging tendency
567
- ᶜdTdt_nudging = p. scratch. ᶜtemp_scalar
568
- ᶜdqtdt_nudging = p. scratch. ᶜtemp_scalar_2
569
- @. ᶜdTdt_nudging =
570
- - (TD. air_temperature (thermo_params, ᶜts) - ᶜT_nudge) * ᶜinv_τ_scalar
571
- @. ᶜdqtdt_nudging = - (ᶜspecific. q_tot - ᶜqt_nudge) * ᶜinv_τ_scalar
572
-
573
- ᶜdTdt_sum = p. scratch. ᶜtemp_scalar
574
- ᶜdqtdt_sum = p. scratch. ᶜtemp_scalar_2
575
- @. ᶜdTdt_sum = ᶜdTdt_hadv + ᶜdTdt_nudging + ᶜdTdt_fluc
576
- @. ᶜdqtdt_sum = ᶜdqtdt_hadv + ᶜdqtdt_nudging + ᶜdqtdt_fluc
577
-
578
- T_0 = TD. Parameters. T_0 (thermo_params)
579
- Lv_0 = TD. Parameters. LH_v0 (thermo_params)
580
- cv_v = TD. Parameters. cv_v (thermo_params)
581
- R_v = TD. Parameters. R_v (thermo_params)
582
- # total energy
583
- @. Yₜ. c. ρe_tot +=
584
- Y. c. ρ * (
585
- TD. cv_m (thermo_params, ᶜts) * ᶜdTdt_sum +
586
- (
587
- cv_v * (TD. air_temperature (thermo_params, ᶜts) - T_0) + Lv_0 -
588
- R_v * T_0
589
- ) * ᶜdqtdt_sum
590
- )
591
- # total specific humidity
592
- @. Yₜ. c. ρq_tot += Y. c. ρ * ᶜdqtdt_sum
593
-
594
- # # subsidence -->
595
- ᶠls_subsidence³ = p. scratch. ᶠtemp_CT3
596
- @. ᶠls_subsidence³ =
597
- ᶠinterp (ᶜls_subsidence * CT3 (unit_basis_vector_data (CT3, ᶜlg)))
598
- subsidence! (
599
- Yₜ. c. ρe_tot,
600
- Y. c. ρ,
601
- ᶠls_subsidence³,
602
- ᶜh_tot,
603
- Val {:first_order} (),
604
- )
605
- subsidence! (
606
- Yₜ. c. ρq_tot,
607
- Y. c. ρ,
608
- ᶠls_subsidence³,
609
- ᶜspecific. q_tot,
610
- Val {:first_order} (),
611
- )
612
-
613
- # needed to address top boundary condition for forcings. Otherwise upper portion of domain is anomalously cold
614
- ρe_tot_top = Fields. level (Yₜ. c. ρe_tot, Spaces. nlevels (axes (Y. c)))
615
- @. ρe_tot_top = 0.0
616
-
617
- ρq_tot_top = Fields. level (Yₜ. c. ρq_tot, Spaces. nlevels (axes (Y. c)))
618
- @. ρq_tot_top = 0.0
619
- # <-- subsidence
620
-
621
- return nothing
622
- end
623
-
624
-
625
- function external_height (ds, :: ERA5Forcing )
626
- vec (mean (ds[" z" ][:], dims = 2 ))
627
- end
628
-
629
- function external_height (ds, :: GCMForcing )
630
- vec (mean (ds[" zg" ][:, :], dims = 2 ))
631
- end
632
-
633
- function external_driven_profile_tmean (ds, varname, :: GCMForcing )
634
- vec (mean (ds[varname][:, :], dims = 2 ))
635
- end
636
-
637
- function external_driven_profile_tmean (ds, varname, :: ERA5Forcing )
638
- ds[varname][:]
639
- end
640
-
641
-
642
- external_forcing_cache (Y, atmos:: AtmosModel , params) =
643
- external_forcing_cache (Y, atmos. external_forcing, params)
644
-
645
- external_forcing_cache (Y, external_forcing:: Nothing , params) = (;)
646
- function external_forcing_cache (
647
- Y,
648
- external_forcing:: Union{GCMForcing, ERA5Forcing} ,
649
- params,
650
- )
651
- # varname = "z" if external_forcing isa ERA5Forcing else "zg"
652
- height_tv = TimeVaryingInput (
653
- axes (Y. c),
654
- " myfile.nc" ,
655
- alpha,
656
- preprocess_func = (x) -> 1 / x,
657
- )
658
- height = similar (Y. c, FT)
659
-
660
- FT = Spaces. undertype (axes (Y. c))
661
- ᶜdTdt_fluc = similar (Y. c, FT)
662
- end
0 commit comments