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

Last change on this file since 1010 was 989, checked in by jisesh, 5 years ago

devel: Added nudging outside the boundary of the domain in addidion to inside nudging
File size: 14.7 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, nudging_relaxation_time, &
15       nudging_relaxation_time_in,nudging_relaxation_time_out
16  !$OMP THREADPRIVATE(center_lon, center_lat, nudging_radius, nudging_relaxation_time,nudging_relaxation_time_in,nudging_relaxation_time_out)
17  TYPE(t_field),POINTER :: f_relax_coef_e(:), f_target_ue(:), &
18       f_relax_coef_i(:), f_target_theta_rhodz(:),f_target_ps(:)
19  CHARACTER(LEN=255),SAVE :: guided_nudging_field
20  INTEGER,SAVE :: nudging_time
21
22  PUBLIC :: init_guided, guided
23 
24CONTAINS
25
26  SUBROUTINE init_guided(f_u,f_theta_rhodz,f_ps)
27    USE getin_mod, ONLY : getin
28    USE math_const, ONLY : pi
29    USE earth_const, ONLY : scale_factor
30    TYPE(t_field),POINTER :: f_u(:)! initial condition
31    TYPE(t_field),POINTER :: f_theta_rhodz(:)! initial condition
32    TYPE(t_field),POINTER :: f_ps(:)! initial condition
33    REAL(rstd), POINTER :: ue(:,:), target_ue(:,:), coef_e(:) 
34    REAL(rstd), POINTER :: theta_rhodz(:,:,:), target_theta_rhodz(:,:,:), coef_i(:) 
35    REAL(rstd), POINTER :: ps2(:), target_ps(:)
36    INTEGER :: ind
37    ! read DEF keys describing how to relax
38    center_lon=0.
39    CALL getin('nudging_center_lon', center_lon)
40    center_lat=0.
41    CALL getin('nudging_center_lat', center_lat)
42    nudging_radius=0.
43    CALL getin('nudging_radius', nudging_radius)
44    nudging_radius = nudging_radius / scale_factor
45    ! we should check that radius>0
46
47    nudging_time=0
48    CALL getin('nudging_time', nudging_time)
49    nudging_time = nudging_time/scale_factor
50    !nudging_time=0
51
52    nudging_relaxation_time_in = 0.
53    nudging_relaxation_time_out = 0.
54    CALL getin('nudging_relaxation_time_in', nudging_relaxation_time_in)
55    CALL getin('nudging_relaxation_time_out', nudging_relaxation_time_out)
56    nudging_relaxation_time_in = nudging_relaxation_time_in/scale_factor
57    nudging_relaxation_time_out = nudging_relaxation_time_out/scale_factor
58    IF(nudging_relaxation_time_in .EQ. 0. .AND. nudging_relaxation_time_out .EQ. 0.) THEN
59        PRINT*," nudging_relaxation_time_in and nudging_relaxation_time_out values can not be 0. at the same time"
60    ENDIF
61
62    CALL getin("guided_nudging_field",guided_nudging_field)
63
64    SELECT CASE(TRIM(guided_nudging_field))
65        CASE('wind','temperature','ps','wind_temperature','wind_temperature_ps')
66            CALL allocate_field(f_relax_coef_e, field_u, type_real, name='nudging_coef_e')
67            CALL allocate_field(f_target_ue, field_u, type_real, llm, name='nudging_target_e')
68            CALL allocate_field(f_relax_coef_i, field_t, type_real, name='nudging_coef_i')
69            CALL allocate_field(f_target_theta_rhodz, field_t, type_real, llm,nqdyn, name='nudging_target_theta')
70            CALL allocate_field(f_target_ps, field_t, type_real, name='nudging_target_ps')
71        CASE DEFAULT
72                 PRINT*,"Bad selector for varaible init_guided_nudging>",TRIM(guided_nudging_field)
73                 PRINT*,"options are <wind>, <temperature>,<ps>,<wind_temperature>,<wind_temperature_ps>"
74                 STOP
75    END SELECT
76    ! compute relax_coef and target_ue
77    center_lon = center_lon * pi/180.
78    center_lat = center_lat * pi/180.
79    DO ind = 1 , ndomain
80        IF (.NOT. assigned_domain(ind)) CYCLE
81            CALL swap_dimensions(ind)
82            CALL swap_geometry(ind) 
83            coef_e = f_relax_coef_e(ind)
84            coef_i = f_relax_coef_i(ind)
85            target_ue = f_target_ue(ind)
86            ue = f_u(ind)
87            target_theta_rhodz = f_target_theta_rhodz(ind)
88            theta_rhodz = f_theta_rhodz(ind)
89            target_ps = f_target_ps(ind)
90            ps2 = f_ps(ind)
91            CALL compute_relax_coef(coef_e, coef_i)
92            CALL compute_target_u(ue, target_ue)
93            CALL compute_target_center(theta_rhodz, target_theta_rhodz)
94            CALL compute_target_center2(ps2, target_ps)
95    END DO
96  END SUBROUTINE init_guided
97
98!----------------------------- Compute relaxation coefficients ------------------------------
99 
100  SUBROUTINE compute_relax_coef(coef_e, coef_i)
101    USE geometry, ONLY : lon_e, lat_e, lon_i, lat_i
102    REAL(rstd), INTENT(OUT) :: coef_e(iim*3*jjm), coef_i(iim*jjm)
103    INTEGER :: l, ij
104    DO ij=ij_begin_ext, ij_end_ext
105       coef_e(ij+u_right) = relax_coef(lon_e(ij+u_right), lat_e(ij+u_right) )
106       coef_e(ij+u_lup)   = relax_coef(lon_e(ij+u_lup),   lat_e(ij+u_lup)   )
107       coef_e(ij+u_ldown) = relax_coef(lon_e(ij+u_ldown), lat_e(ij+u_ldown) )
108       coef_i(ij) = relax_coef(lon_i(ij), lat_i(ij) )
109    END DO
110  END SUBROUTINE compute_relax_coef
111
112  FUNCTION relax_coef(lon,lat)
113    USE spherical_geom_mod, ONLY : dist_lonlat
114    USE time_mod!, ONLY : dt ! time step for a full RK step
115    REAL(rstd), INTENT(IN) :: lon,lat
116    REAL(rstd) :: relax_coef, dist, c
117    ! NB : dist is computed on unit sphere
118    IF(nudging_radius>0.) THEN
119       CALL dist_lonlat(lon, lat, center_lon, center_lat, dist) 
120       c = tanh((1.-radius*dist/nudging_radius)*20.) ! 1 inside circle, -1 outside
121       IF (c>0.99) c=1
122       IF (c<-0.99) c=-1
123       c = .5*(1.+c) ! rescale to [0,1] range ; ! 1 inside circle, 0 outside
124    ELSE
125       c = 0.
126    END IF
127
128    nudging_relaxation_time = c*nudging_relaxation_time_in + (1.-c)*nudging_relaxation_time_out 
129    relax_coef = 1. - dt/(dt + nudging_relaxation_time)
130  END FUNCTION relax_coef
131
132!----------------------------- Copy initial condition as target ------------------------------
133
134  SUBROUTINE compute_target_u(ue, target_ue)
135    REAL(rstd), INTENT(OUT)  :: target_ue(iim*3*jjm,llm)
136    REAL(rstd), INTENT(IN) :: ue(iim*3*jjm,llm)
137    INTEGER :: l, ij
138    DO l = ll_begin, ll_end
139       DO ij=ij_begin_ext, ij_end_ext
140          target_ue(ij+u_right,l)=ue(ij+u_right,l)
141          target_ue(ij+u_lup,l)=ue(ij+u_lup,l)
142          target_ue(ij+u_ldown,l)=ue(ij+u_ldown,l)
143       END DO
144    END DO
145  END SUBROUTINE compute_target_u
146
147  SUBROUTINE compute_target_center(theta_rhodz, target_theta_rhodz)
148    REAL(rstd), INTENT(OUT)  :: target_theta_rhodz(iim*jjm,llm,nqdyn)
149    REAL(rstd), INTENT(IN) :: theta_rhodz(iim*jjm,llm,nqdyn)
150    INTEGER :: l, ij, iq
151    DO iq=1, nqdyn
152       DO l = ll_begin, ll_end
153          DO ij=ij_begin_ext, ij_end_ext
154             target_theta_rhodz(ij,l,iq)=theta_rhodz(ij,l,iq)
155          END DO
156       END DO
157    END DO
158  END SUBROUTINE compute_target_center
159
160  SUBROUTINE compute_target_center2(theta_rhodz, target_theta_rhodz)
161    REAL(rstd), INTENT(OUT)  :: target_theta_rhodz(iim*jjm,nqdyn)
162    REAL(rstd), INTENT(IN) :: theta_rhodz(iim*jjm,nqdyn)
163    INTEGER :: ij, iq
164    DO iq=1, nqdyn
165          DO ij=ij_begin_ext, ij_end_ext
166             target_theta_rhodz(ij,iq)=theta_rhodz(ij,iq)
167          END DO
168    END DO
169  END SUBROUTINE compute_target_center2
170
171!----------------------------- Relax towards target ------------------------------
172 
173  SUBROUTINE guided(tt, f_ps, f_theta_rhodz, f_u, f_q)
174    USE xios_mod
175    USE theta2theta_rhodz_mod
176    USE wind_mod
177    USE transfert_mod
178    USE time_mod, ONLY : dt ! time step for a full RK step
179    USE vertical_remap_mod
180    USE compute_pression_mod, ONLY : pression_mid
181    REAL(rstd), INTENT(IN):: tt
182    TYPE(t_field),POINTER :: f_ps(:)
183    TYPE(t_field),POINTER,SAVE :: f_pmid_target(:)
184    TYPE(t_field),POINTER :: f_phis(:)
185    TYPE(t_field),POINTER :: f_theta_rhodz(:)
186    TYPE(t_field),POINTER :: f_u(:) 
187    TYPE(t_field),POINTER :: f_q(:)
188    TYPE(t_field),POINTER, SAVE :: f_T_guided(:),  f_ulon_guided(:), f_ulat_guided(:),f_T_guided_interp(:), &
189                                   f_ulon_guided_interp(:),f_ulat_guided_interp(:),f_ps_guided(:)
190    REAL(rstd), POINTER :: target_ue(:,:), ue(:,:), coef_e(:)
191    REAL(rstd), POINTER :: target_theta_rhodz(:,:,:), theta_rhodz(:,:,:), coef_i(:) 
192    REAL(rstd), POINTER :: target_ps(:), ps2(:)
193    INTEGER :: ind
194   
195    IF (abs(MOD(tt,REAL(nudging_time))-dt) < 1.0D-2) THEN
196      CALL allocate_field(f_ps_guided, field_t, type_real, name='nudging_ps')
197      CALL xios_read_field("ps_guided_read",f_ps_guided)
198      CALL transfert_request(f_ps_guided,req_i0)
199      CALL transfert_request(f_ps_guided,req_i1)
200      CALL allocate_field(f_T_guided, field_t, type_real, llm, name='nudging_T')
201      CALL allocate_field(f_T_guided_interp, field_t, type_real, llm, name='nudging_T')
202      CALL allocate_field(f_ulon_guided_interp, field_t, type_real, llm, name='nudging_T')
203      CALL allocate_field(f_ulat_guided_interp, field_t, type_real, llm, name='nudging_T')
204      CALL allocate_field(f_ulon_guided, field_t, type_real, llm, name='nudging_ulon')
205      CALL allocate_field(f_ulat_guided, field_t, type_real, llm, name='nudging_ulat')
206      CALL allocate_field(f_pmid_target,field_t,type_real,llm,name='target_pressure') 
207      CALL xios_read_field("T_guided_read",f_T_guided)
208      CALL pression_mid(f_ps, f_pmid_target) 
209      CALL vertical_remap_extdata(f_T_guided,f_pmid_target,f_T_guided_interp)
210      CALL transfert_request(f_T_guided,req_i0)
211      CALL transfert_request(f_T_guided_interp,req_i0)
212      CALL transfert_request(f_T_guided,req_i1)
213      CALL transfert_request(f_T_guided_interp,req_i1)
214      CALL temperature2theta_rhodz(f_ps,f_T_guided_interp,f_target_theta_rhodz)
215      CALL xios_read_field("ulon_guided_read",f_ulon_guided)
216      CALL xios_read_field("ulat_guided_read",f_ulat_guided)
217      CALL vertical_remap_extdata(f_ulon_guided,f_pmid_target,f_ulon_guided_interp)
218      CALL vertical_remap_extdata(f_ulat_guided,f_pmid_target,f_ulat_guided_interp)
219      CALL transfert_request(f_ulon_guided,req_i0)
220      CALL transfert_request(f_ulon_guided_interp,req_i0)
221      CALL transfert_request(f_ulat_guided,req_i0)
222      CALL transfert_request(f_ulat_guided_interp,req_i0)
223      CALL transfert_request(f_ulon_guided,req_i1)
224      CALL transfert_request(f_ulon_guided_interp,req_i1)
225      CALL transfert_request(f_ulat_guided,req_i1)
226      CALL transfert_request(f_ulat_guided_interp,req_i1)
227      !CALL xios_write_field("ulat_guided_out",f_ulat_guided)
228      CALL xios_write_field("ulat_guided_out",f_T_guided_interp) ! FIX
229      CALL xios_write_field("ulon_guided_out",f_ulon_guided)
230      CALL xios_write_field("T_guided_out",f_T_guided)
231      CALL xios_write_field("PS_guided_out",f_ps_guided) 
232      CALL ulonlat2un(f_ulon_guided_interp, f_ulat_guided,f_target_ue) 
233      CALL deallocate_field(f_T_guided)
234      CALL deallocate_field(f_T_guided_interp)
235      CALL deallocate_field(f_ulon_guided)
236      CALL deallocate_field(f_ulon_guided_interp)
237      CALL deallocate_field(f_ulat_guided)
238      CALL deallocate_field(f_ulat_guided_interp)
239      CALL deallocate_field(f_pmid_target)
240    ENDIF
241
242
243    DO ind = 1 , ndomain
244       IF (.NOT. assigned_domain(ind)) CYCLE
245       CALL swap_dimensions(ind)
246       CALL swap_geometry(ind) 
247       
248       SELECT CASE(TRIM(guided_nudging_field))
249         CASE ('wind')
250           coef_e = f_relax_coef_e(ind)
251           target_ue = f_target_ue(ind)
252           ue = f_u(ind)
253           CALL compute_guided_u(coef_e, target_ue, ue)
254         CASE ('temperature')
255           coef_i = f_relax_coef_i(ind)
256           target_theta_rhodz = f_target_theta_rhodz(ind)
257           theta_rhodz = f_theta_rhodz(ind)
258           CALL compute_guided_center(coef_i, target_theta_rhodz, theta_rhodz)
259         CASE ('ps')
260           coef_i = f_relax_coef_i(ind)
261           target_ps = f_ps_guided(ind)
262           ps2 = f_ps(ind)
263           CALL compute_guided_center2(coef_i, target_ps, ps2)
264         CASE ('wind_temperature')
265           coef_e = f_relax_coef_e(ind)
266           target_ue = f_target_ue(ind)
267           ue = f_u(ind)
268           CALL compute_guided_u(coef_e, target_ue, ue)
269           coef_i = f_relax_coef_i(ind)
270           target_theta_rhodz = f_target_theta_rhodz(ind)
271           theta_rhodz = f_theta_rhodz(ind)
272           CALL compute_guided_center(coef_i, target_theta_rhodz, theta_rhodz)
273         CASE ('wind_temperature_ps')
274           coef_e = f_relax_coef_e(ind)
275           target_ue = f_target_ue(ind)
276           ue = f_u(ind)
277           CALL compute_guided_u(coef_e, target_ue, ue)
278           coef_i = f_relax_coef_i(ind)
279           target_theta_rhodz = f_target_theta_rhodz(ind)
280           theta_rhodz = f_theta_rhodz(ind)
281           CALL compute_guided_center(coef_i, target_theta_rhodz, theta_rhodz)
282           target_ps = f_ps_guided(ind)
283           ps2 = f_ps(ind)
284           CALL compute_guided_center2(coef_i, target_ps, ps2)
285         CASE DEFAULT
286                 PRINT*,"Bad selector for varaible <guided_nudging_field>",TRIM(guided_nudging_field)
287                 PRINT*,"options are <wind>, <temperature>,<ps>,<wind_temperature>,<wind_temperature_ps>"
288                 STOP
289       END SELECT
290    END DO
291
292  END SUBROUTINE guided
293
294  SUBROUTINE compute_guided_u(coef_e, target_ue, ue)
295    REAL(rstd), INTENT(IN) :: coef_e(iim*3*jjm)
296    REAL(rstd), INTENT(IN) :: target_ue(iim*3*jjm,llm)
297    REAL(rstd), INTENT(INOUT) :: ue(iim*3*jjm,llm)
298    INTEGER :: l, ij
299    DO l = ll_begin, ll_end
300       DO ij=ij_begin_ext, ij_end_ext
301          ue(ij+u_right,l) = ue(ij+u_right,l)*coef_e(ij+u_right) + &
302               target_ue(ij+u_right,l)*(1.-coef_e(ij+u_right))
303          ue(ij+u_lup,l) = ue(ij+u_lup,l)*coef_e(ij+u_lup) + &
304               target_ue(ij+u_lup,l)*(1.-coef_e(ij+u_lup))
305          ue(ij+u_ldown,l) = ue(ij+u_ldown,l)*coef_e(ij+u_ldown) + &
306               target_ue(ij+u_ldown,l)*(1.-coef_e(ij+u_ldown))
307       END DO
308    END DO
309  END SUBROUTINE compute_guided_u
310
311  SUBROUTINE compute_guided_center(coef_i, target_theta_rhodz, theta_rhodz)
312    REAL(rstd), INTENT(IN) :: coef_i(iim*jjm)
313    REAL(rstd), INTENT(IN) :: target_theta_rhodz(iim*jjm,llm,nqdyn)
314    REAL(rstd), INTENT(INOUT) :: theta_rhodz(iim*jjm,llm,nqdyn)
315    INTEGER :: l, ij, iq
316    DO iq=1, nqdyn
317       DO l = ll_begin, ll_end
318          DO ij=ij_begin_ext, ij_end_ext
319             theta_rhodz(ij,l,iq) = theta_rhodz(ij,l,iq)*coef_i(ij) + &
320                  target_theta_rhodz(ij,l,iq)*(1.-coef_i(ij))
321          END DO
322       END DO
323    END DO
324  END SUBROUTINE compute_guided_center
325
326  SUBROUTINE compute_guided_center2(coef_i, target_theta_rhodz, theta_rhodz)
327    REAL(rstd), INTENT(IN) :: coef_i(iim*jjm)
328    REAL(rstd), INTENT(IN) :: target_theta_rhodz(iim*jjm,nqdyn)
329    REAL(rstd), INTENT(INOUT) :: theta_rhodz(iim*jjm,nqdyn)
330    INTEGER :: ij, iq
331    DO iq=1, nqdyn
332          DO ij=ij_begin_ext, ij_end_ext
333             theta_rhodz(ij,iq) = theta_rhodz(ij,iq)*coef_i(ij) + &
334                  target_theta_rhodz(ij,iq)*(1.-coef_i(ij))
335          END DO
336    END DO
337  END SUBROUTINE compute_guided_center2
338 
339END MODULE nudging_mod
Note: See TracBrowser for help on using the repository browser.