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

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

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