Line data Source code
1 : module rayleigh_friction
2 :
3 : !---------------------------------------------------------------------------------
4 : ! Module to apply rayleigh friction in region of model top.
5 : ! We specify a decay rate profile that is largest at the model top and
6 : ! drops off vertically using a hyperbolic tangent profile.
7 : ! We compute the tendencies in u and v using an Euler backward scheme.
8 : ! We then apply the negative of the kinetic energy tendency to "s", the dry
9 : ! static energy.
10 : !---------------------------Code history--------------------------------
11 : ! This is a new routine written by Art Mirin in collaboration with Phil Rasch.
12 : ! Initial coding for this version: Art Mirin, May 2007.
13 : ! First CCPP conversion by Andrew Gettelman, Sep 2021
14 : ! Final CCPP completion by Kate T-C, Oct 2024
15 : !---------------------------------------------------------------------------------
16 :
17 : use ccpp_kinds, only: kind_phys
18 :
19 : implicit none
20 :
21 : save
22 : private ! Make default type private to the module
23 : !
24 : ! PUBLIC: interfaces
25 : !
26 : public rayleigh_friction_init
27 : public rayleigh_friction_run
28 :
29 : ! PRIVATE: module variables
30 : integer :: rayk0 = 2 ! vertical level at which rayleigh friction term is centered
31 : real(kind_phys) :: raykrange = 0._kind_phys ! range of rayleigh friction profile
32 : ! if 0, range is set to satisfy x=2 (see below)
33 : real(kind_phys) :: raytau0 = 0._kind_phys ! approximate value of decay time at model top (days)
34 : ! if 0., no rayleigh friction is applied
35 : ! Local
36 : real (kind_phys) :: krange ! range of rayleigh friction profile
37 : real (kind_phys) :: tau0 ! approximate value of decay time at model top
38 : real (kind_phys) :: otau0 ! inverse of tau0
39 : real (kind_phys), allocatable :: otau(:) ! inverse decay time versus vertical level
40 :
41 : ! We apply a profile of the form otau0 * [1 + tanh (x)] / 2 , where
42 : ! x = (k0 - k) / krange. The default is for x to equal 2 at k=1, meaning
43 : ! krange = (k0 - 1) / 2. The default is applied when raykrange is set to 0.
44 : ! If otau0 = 0, no term is applied.
45 :
46 : contains
47 :
48 : !> \section arg_table_rayleigh_friction_init Argument Table
49 : !! \htmlinclude rayleigh_friction_init.html
50 1024 : subroutine rayleigh_friction_init(pver, raytau0_nl, raykrange_nl, rayk0_nl, masterproc, iulog, errmsg, errflg)
51 :
52 : integer, intent(in) :: pver
53 : real (kind_phys), intent(in) :: raytau0_nl
54 : real (kind_phys), intent(in) :: raykrange_nl
55 : integer, intent(in) :: rayk0_nl
56 : logical, intent(in) :: masterproc
57 : integer, intent(in) :: iulog
58 : character(len=512), intent(out) :: errmsg
59 : integer, intent(out) :: errflg
60 :
61 : ! local variables
62 : character(len=*), parameter :: subname = 'rayleigh_friction_init'
63 : real (kind_phys) x
64 : integer k, ierr
65 :
66 1024 : errmsg = ''
67 1024 : errflg = 0
68 1024 : raytau0 = raytau0_nl
69 1024 : raykrange = raykrange_nl
70 1024 : rayk0 = rayk0_nl
71 :
72 1024 : if ( raytau0 > 0._kind_phys) then ! Rayleigh Friction enabled
73 : ! Allocate module data
74 0 : allocate(otau(pver), stat=ierr)
75 0 : if (ierr /= 0) then
76 0 : errflg = ierr
77 0 : errmsg = trim(subname)//': Allocate of otau failed'
78 0 : return
79 : end if
80 :
81 : !-----------------------------------------------------------------------
82 : ! Compute tau array
83 : !-----------------------------------------------------------------------
84 :
85 0 : krange = raykrange
86 0 : if (raykrange .eq. 0._kind_phys) krange = (rayk0 - 1) / 2._kind_phys
87 :
88 0 : tau0 = (86400._kind_phys) * raytau0 ! convert to seconds
89 0 : otau0 = 0._kind_phys
90 0 : if (tau0 .ne. 0._kind_phys) otau0 = 1._kind_phys/tau0
91 :
92 0 : do k = 1, pver
93 0 : x = (rayk0 - k) / krange
94 0 : otau(k) = otau0 * (1 + tanh(x)) / (2._kind_phys)
95 : enddo
96 : end if ! Rayleigh Friction enabled
97 :
98 1024 : if ( masterproc ) then
99 2 : if (raytau0 > 0._kind_phys) then
100 0 : write(iulog,*) 'tuning parameters rayleigh_friction_init: rayk0',rayk0
101 0 : write(iulog,*) 'tuning parameters rayleigh_friction_init: raykrange',raykrange
102 0 : write(iulog,*) 'tuning parameters rayleigh_friction_init: raytau0',raytau0
103 : else
104 2 : write(iulog,*) 'Rayleigh Friction not enabled.'
105 : endif
106 : endif
107 :
108 : end subroutine rayleigh_friction_init
109 :
110 : !=========================================================================================
111 : !> \section arg_table_rayleigh_friction_run Argument Table
112 : !! \htmlinclude rayleigh_friction_run.html
113 63688 : subroutine rayleigh_friction_run(pver, ztodt, u, v, dudt, dvdt, dsdt, errmsg, errflg)
114 :
115 : !------------------------------Arguments--------------------------------
116 : integer, intent(in) :: pver
117 : real(kind_phys), intent(in) :: ztodt !physics timestep
118 : real(kind_phys), intent(in) :: u(:,:)
119 : real(kind_phys), intent(in) :: v(:,:)
120 : real(kind_phys), intent(out) :: dudt(:,:) !tendency_of_eastward_wind
121 : real(kind_phys), intent(out) :: dvdt(:,:) !tendency_of_northward_wind
122 : real(kind_phys), intent(out) :: dsdt(:,:) !heating_rate
123 :
124 : character(len=512), intent(out) :: errmsg
125 : integer, intent(out) :: errflg
126 :
127 : !---------------------------Local storage-------------------------------
128 : character(len=*), parameter :: subname = 'rayleigh_friction_run'
129 : integer :: k ! level
130 : real(kind_phys) :: rztodt ! 1./ztodt
131 : real(kind_phys) :: c1, c2, c3 ! temporary variables
132 : !-----------------------------------------------------------------------
133 :
134 : ! initialize values
135 63688 : errmsg = ''
136 63688 : errflg = 0
137 :
138 63688 : if (raytau0 .eq. 0._kind_phys) return
139 :
140 0 : rztodt = 1._kind_phys/ztodt
141 :
142 : ! u, v and s are modified by rayleigh friction
143 :
144 0 : do k = 1, pver
145 0 : c2 = 1._kind_phys / (1._kind_phys + otau(k)*ztodt)
146 0 : c1 = -otau(k) * c2
147 0 : c3 = 0.5_kind_phys * (1._kind_phys - c2*c2) * rztodt
148 0 : dudt(:,k) = c1 * u(:,k)
149 0 : dvdt(:,k) = c1 * v(:,k)
150 0 : dsdt(:,k) = c3 * (u(:,k)**2 + v(:,k)**2)
151 : enddo
152 :
153 : end subroutine rayleigh_friction_run
154 :
155 : end module rayleigh_friction
156 :
|
