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

Last change on this file was 1026, checked in by dubos, 4 years ago

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