source: codes/icosagcm/devel/src/dissip/nudging_mod.f90 @ 846

Last change on this file since 846 was 820, checked in by jisesh, 5 years ago

devel : towards nudging (status : relax wind + theta towards initial condition)
File size: 7.2 KB
Line 
1MODULE nudging_mod
2  USE icosa, ONLY : rstd
3  USE grid_param, ONLY : llm, nqdyn
4  USE omp_para, ONLY : ll_begin, ll_end
5  USE domain_mod, ONLY : ndomain, assigned_domain
6  USE dimensions, ONLY : swap_dimensions, u_right, u_lup, u_ldown
7  USE dimensions, ONLY : iim, jjm, ij_begin_ext, ij_end_ext
8  USE geometry, ONLY : swap_geometry
9  USE field_mod
10  IMPLICIT NONE
11  SAVE
12  PRIVATE
13  ! nudging will be active outside a disc of radius 'radius' centered at 'center_lon', 'center lat'.
14  REAL(rstd) :: center_lon, center_lat, nudging_radius, time
15  !$OMP THREADPRIVATE(center_lon, center_lat, nudging_radius, time)
16
17  TYPE(t_field),POINTER :: f_relax_coef_e(:), f_target_ue(:), &
18       f_relax_coef_i(:), f_target_theta_rhodz(:)
19
20  PUBLIC :: init_guided, guided
21 
22CONTAINS
23
24  SUBROUTINE init_guided(f_u,f_theta_rhodz)
25    USE getin_mod, ONLY : getin
26    USE math_const, ONLY : pi
27    USE earth_const, ONLY : scale_factor
28    TYPE(t_field),POINTER :: f_u(:)! initial condition
29    TYPE(t_field),POINTER :: f_theta_rhodz(:)! initial condition
30    REAL(rstd), POINTER :: ue(:,:), target_ue(:,:), coef_e(:) 
31    REAL(rstd), POINTER :: theta_rhodz(:,:,:), target_theta_rhodz(:,:,:), coef_i(:) 
32    INTEGER :: ind
33    ! read DEF keys describing how to relax
34    center_lon=0.
35    CALL getin('nudging_center_lon', center_lon)
36    center_lat=0.
37    CALL getin('nudging_center_lat', center_lat)
38    nudging_radius=0.
39    CALL getin('nudging_radius', nudging_radius)
40    nudging_radius = nudging_radius / scale_factor
41!    time=0.
42!    CALL getin('nudging_time', time)
43
44    ! we should check that radius>0
45
46    !SELECT CASE(TRIM(nudg_name))
47    !    CASE ('f_u')
48        CALL allocate_field(f_relax_coef_e, field_u, type_real, name='nudging_coef_e')
49        CALL allocate_field(f_target_ue, field_u, type_real, llm, name='nudging_target_e')
50    !    CASE ('f_theta_rhodz')
51        CALL allocate_field(f_relax_coef_i, field_t, type_real, name='nudging_coef_i')
52        CALL allocate_field(f_target_theta_rhodz, field_t, type_real, llm,nqdyn, name='nudging_target_theta')
53    !    CASE DEFAULT
54    !END SELECT
55
56    ! compute relax_coef and target_ue
57    center_lon = center_lon * pi/180.
58    center_lat = center_lat * pi/180.
59    DO ind = 1 , ndomain
60        IF (.NOT. assigned_domain(ind)) CYCLE
61            CALL swap_dimensions(ind)
62            CALL swap_geometry(ind) 
63            coef_e = f_relax_coef_e(ind)
64            coef_i = f_relax_coef_i(ind)
65            target_ue = f_target_ue(ind)
66            ue = f_u(ind)
67            target_theta_rhodz = f_target_theta_rhodz(ind)
68            theta_rhodz = f_theta_rhodz(ind)
69            CALL compute_relax_coef(coef_e, coef_i)
70            CALL compute_target_u(ue, target_ue)
71            CALL compute_target_center(theta_rhodz, target_theta_rhodz)
72    END DO
73  END SUBROUTINE init_guided
74
75!----------------------------- Compute relaxation coefficients ------------------------------
76 
77  SUBROUTINE compute_relax_coef(coef_e, coef_i)
78    USE geometry, ONLY : lon_e, lat_e, lon_i, lat_i
79    REAL(rstd), INTENT(OUT) :: coef_e(iim*3*jjm), coef_i(iim*jjm)
80    INTEGER :: l, ij
81    DO ij=ij_begin_ext, ij_end_ext
82       coef_e(ij+u_right) = relax_coef(lon_e(ij+u_right), lat_e(ij+u_right) )
83       coef_e(ij+u_lup)   = relax_coef(lon_e(ij+u_lup),   lat_e(ij+u_lup)   )
84       coef_e(ij+u_ldown) = relax_coef(lon_e(ij+u_ldown), lat_e(ij+u_ldown) )
85       coef_i(ij) = relax_coef(lon_i(ij), lat_i(ij) )
86    END DO
87  END SUBROUTINE compute_relax_coef
88
89  FUNCTION relax_coef(lon,lat)
90    USE spherical_geom_mod, ONLY : dist_lonlat
91    REAL(rstd), INTENT(IN) :: lon,lat
92    REAL(rstd) :: relax_coef, dist, c
93    ! NB : dist is computed on unit sphere
94    CALL dist_lonlat(lon, lat, center_lon, center_lat, dist) 
95    c = tanh((1.-radius*dist/nudging_radius)*20.) ! 1 inside circle, -1 outside
96    relax_coef = .5*(1.+c) ! rescale to [0,1] range ; ! 1 inside circle, 0 outside
97  END FUNCTION relax_coef
98
99!----------------------------- Copy initial condition as target ------------------------------
100
101  SUBROUTINE compute_target_u(ue, target_ue)
102    REAL(rstd), INTENT(OUT)  :: target_ue(iim*3*jjm,llm)
103    REAL(rstd), INTENT(IN) :: ue(iim*3*jjm,llm)
104    INTEGER :: l, ij
105    DO l = ll_begin, ll_end
106       DO ij=ij_begin_ext, ij_end_ext
107          target_ue(ij+u_right,l)=ue(ij+u_right,l)
108          target_ue(ij+u_lup,l)=ue(ij+u_lup,l)
109          target_ue(ij+u_ldown,l)=ue(ij+u_ldown,l)
110       END DO
111    END DO
112  END SUBROUTINE compute_target_u
113
114  SUBROUTINE compute_target_center(theta_rhodz, target_theta_rhodz)
115    REAL(rstd), INTENT(OUT)  :: target_theta_rhodz(iim*jjm,llm,nqdyn)
116    REAL(rstd), INTENT(IN) :: theta_rhodz(iim*jjm,llm,nqdyn)
117    INTEGER :: l, ij, iq
118    DO iq=1, nqdyn
119       DO l = ll_begin, ll_end
120          DO ij=ij_begin_ext, ij_end_ext
121             target_theta_rhodz(ij,l,iq)=theta_rhodz(ij,l,iq)
122          END DO
123       END DO
124    END DO
125  END SUBROUTINE compute_target_center
126
127!----------------------------- Relax towards target ------------------------------
128 
129  SUBROUTINE guided(tt, f_ps, f_theta_rhodz, f_u, f_q)
130    REAL(rstd), INTENT(IN):: tt
131    TYPE(t_field),POINTER :: f_ps(:)
132    TYPE(t_field),POINTER :: f_phis(:)
133    TYPE(t_field),POINTER :: f_theta_rhodz(:)
134    TYPE(t_field),POINTER :: f_u(:) 
135    TYPE(t_field),POINTER :: f_q(:) 
136    REAL(rstd), POINTER :: target_ue(:,:), ue(:,:), coef_e(:)
137    REAL(rstd), POINTER :: target_theta_rhodz(:,:,:), theta_rhodz(:,:,:), coef_i(:) 
138    INTEGER :: ind
139
140    DO ind = 1 , ndomain
141       IF (.NOT. assigned_domain(ind)) CYCLE
142       CALL swap_dimensions(ind)
143       CALL swap_geometry(ind) 
144       coef_e = f_relax_coef_e(ind)
145       target_ue = f_target_ue(ind)
146       ue = f_u(ind)
147       CALL compute_guided_u(coef_e, target_ue, ue)
148       coef_i = f_relax_coef_i(ind)
149       target_theta_rhodz = f_target_theta_rhodz(ind)
150       theta_rhodz = f_theta_rhodz(ind)
151       CALL compute_guided_center(coef_i, target_theta_rhodz, theta_rhodz)
152    END DO
153
154  END SUBROUTINE guided
155
156  SUBROUTINE compute_guided_u(coef_e, target_ue, ue)
157    REAL(rstd), INTENT(IN) :: coef_e(iim*3*jjm)
158    REAL(rstd), INTENT(IN) :: target_ue(iim*3*jjm,llm)
159    REAL(rstd), INTENT(INOUT) :: ue(iim*3*jjm,llm)
160    INTEGER :: l, ij
161    DO l = ll_begin, ll_end
162       DO ij=ij_begin_ext, ij_end_ext
163          ue(ij+u_right,l) = ue(ij+u_right,l)*coef_e(ij+u_right) + &
164               target_ue(ij+u_right,l)*(1.-coef_e(ij+u_right))
165          ue(ij+u_lup,l) = ue(ij+u_lup,l)*coef_e(ij+u_lup) + &
166               target_ue(ij+u_lup,l)*(1.-coef_e(ij+u_lup))
167          ue(ij+u_ldown,l) = ue(ij+u_ldown,l)*coef_e(ij+u_ldown) + &
168               target_ue(ij+u_ldown,l)*(1.-coef_e(ij+u_ldown))
169       END DO
170    END DO
171  END SUBROUTINE compute_guided_u
172
173  SUBROUTINE compute_guided_center(coef_i, target_theta_rhodz, theta_rhodz)
174    REAL(rstd), INTENT(IN) :: coef_i(iim*jjm)
175    REAL(rstd), INTENT(IN) :: target_theta_rhodz(iim*jjm,llm,nqdyn)
176    REAL(rstd), INTENT(INOUT) :: theta_rhodz(iim*jjm,llm,nqdyn)
177    INTEGER :: l, ij, iq
178    DO iq=1, nqdyn
179       DO l = ll_begin, ll_end
180          DO ij=ij_begin_ext, ij_end_ext
181             theta_rhodz(ij,l,iq) = theta_rhodz(ij,l,iq)*coef_i(ij) + &
182                  target_theta_rhodz(ij,l,iq)*(1.-coef_i(ij))
183          END DO
184       END DO
185    END DO
186  END SUBROUTINE compute_guided_center
187 
188END MODULE nudging_mod
Note: See TracBrowser for help on using the repository browser.