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

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

devel: Nudging external data and vertical interpolatiion
File size: 10.8 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  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)
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    REAL(rstd), POINTER :: ue(:,:), target_ue(:,:), coef_e(:) 
33    REAL(rstd), POINTER :: theta_rhodz(:,:,:), target_theta_rhodz(:,:,:), coef_i(:) 
34    INTEGER :: ind
35    ! read DEF keys describing how to relax
36    center_lon=0.
37    CALL getin('nudging_center_lon', center_lon)
38    center_lat=0.
39    CALL getin('nudging_center_lat', center_lat)
40    nudging_radius=0.
41    CALL getin('nudging_radius', nudging_radius)
42    nudging_radius = nudging_radius / scale_factor
43
44    nudging_time=0.
45    CALL getin('nudging_time', nudging_time)
46    nudging_time = nudging_time/scale_factor
47    ! we should check that radius>0
48    CALL getin("guided_nudging_field",guided_nudging_field)
49
50    SELECT CASE(TRIM(guided_nudging_field))
51        CASE ('wind')
52            CALL allocate_field(f_relax_coef_e, field_u, type_real, name='nudging_coef_e')
53            CALL allocate_field(f_target_ue, field_u, type_real, llm, name='nudging_target_e')
54        CASE ('temperature')
55            CALL allocate_field(f_relax_coef_i, field_t, type_real, name='nudging_coef_i')
56            CALL allocate_field(f_target_theta_rhodz, field_t, type_real, llm,nqdyn, name='nudging_target_theta')
57        CASE ('wind_temperature')
58            CALL allocate_field(f_relax_coef_e, field_u, type_real, name='nudging_coef_e')
59            CALL allocate_field(f_target_ue, field_u, type_real, llm, name='nudging_target_e')
60            CALL allocate_field(f_relax_coef_i, field_t, type_real, name='nudging_coef_i')
61            CALL allocate_field(f_target_theta_rhodz, field_t, type_real, llm,nqdyn, name='nudging_target_theta')
62        CASE DEFAULT
63                 PRINT*,"Bad selector for varaible init_guided_nudging>",TRIM(guided_nudging_field)
64                 PRINT*," option are <none>, <wind>, <temperature>,<wind_temperature>"
65                 STOP
66    END SELECT
67
68    ! compute relax_coef and target_ue
69    center_lon = center_lon * pi/180.
70    center_lat = center_lat * pi/180.
71    DO ind = 1 , ndomain
72        IF (.NOT. assigned_domain(ind)) CYCLE
73            CALL swap_dimensions(ind)
74            CALL swap_geometry(ind) 
75            coef_e = f_relax_coef_e(ind)
76            coef_i = f_relax_coef_i(ind)
77            target_ue = f_target_ue(ind)
78            ue = f_u(ind)
79            target_theta_rhodz = f_target_theta_rhodz(ind)
80            theta_rhodz = f_theta_rhodz(ind)
81            CALL compute_relax_coef(coef_e, coef_i)
82            CALL compute_target_u(ue, target_ue)
83            CALL compute_target_center(theta_rhodz, target_theta_rhodz)
84    END DO
85  END SUBROUTINE init_guided
86
87!----------------------------- Compute relaxation coefficients ------------------------------
88 
89  SUBROUTINE compute_relax_coef(coef_e, coef_i)
90    USE geometry, ONLY : lon_e, lat_e, lon_i, lat_i
91    REAL(rstd), INTENT(OUT) :: coef_e(iim*3*jjm), coef_i(iim*jjm)
92    INTEGER :: l, ij
93    DO ij=ij_begin_ext, ij_end_ext
94       coef_e(ij+u_right) = relax_coef(lon_e(ij+u_right), lat_e(ij+u_right) )
95       coef_e(ij+u_lup)   = relax_coef(lon_e(ij+u_lup),   lat_e(ij+u_lup)   )
96       coef_e(ij+u_ldown) = relax_coef(lon_e(ij+u_ldown), lat_e(ij+u_ldown) )
97       coef_i(ij) = relax_coef(lon_i(ij), lat_i(ij) )
98    END DO
99  END SUBROUTINE compute_relax_coef
100
101  FUNCTION relax_coef(lon,lat)
102    USE spherical_geom_mod, ONLY : dist_lonlat
103    REAL(rstd), INTENT(IN) :: lon,lat
104    REAL(rstd) :: relax_coef, dist, c
105    ! NB : dist is computed on unit sphere
106    CALL dist_lonlat(lon, lat, center_lon, center_lat, dist) 
107    c = tanh((1.-radius*dist/nudging_radius)*20.) ! 1 inside circle, -1 outside
108    IF (c>0.99) c=1
109    IF (c<-0.99) c=-1
110    relax_coef = .5*(1.+c) ! rescale to [0,1] range ; ! 1 inside circle, 0 outside
111  END FUNCTION relax_coef
112
113!----------------------------- Copy initial condition as target ------------------------------
114
115  SUBROUTINE compute_target_u(ue, target_ue)
116    REAL(rstd), INTENT(OUT)  :: target_ue(iim*3*jjm,llm)
117    REAL(rstd), INTENT(IN) :: ue(iim*3*jjm,llm)
118    INTEGER :: l, ij
119    DO l = ll_begin, ll_end
120       DO ij=ij_begin_ext, ij_end_ext
121          target_ue(ij+u_right,l)=ue(ij+u_right,l)
122          target_ue(ij+u_lup,l)=ue(ij+u_lup,l)
123          target_ue(ij+u_ldown,l)=ue(ij+u_ldown,l)
124       END DO
125    END DO
126  END SUBROUTINE compute_target_u
127
128  SUBROUTINE compute_target_center(theta_rhodz, target_theta_rhodz)
129    REAL(rstd), INTENT(OUT)  :: target_theta_rhodz(iim*jjm,llm,nqdyn)
130    REAL(rstd), INTENT(IN) :: theta_rhodz(iim*jjm,llm,nqdyn)
131    INTEGER :: l, ij, iq
132    DO iq=1, nqdyn
133       DO l = ll_begin, ll_end
134          DO ij=ij_begin_ext, ij_end_ext
135             target_theta_rhodz(ij,l,iq)=theta_rhodz(ij,l,iq)
136          END DO
137       END DO
138    END DO
139  END SUBROUTINE compute_target_center
140
141!----------------------------- Relax towards target ------------------------------
142 
143  SUBROUTINE guided(tt, f_ps, f_theta_rhodz, f_u, f_q)
144    USE xios_mod
145    USE theta2theta_rhodz_mod
146    USE wind_mod
147    USE transfert_mod
148    USE time_mod
149    USE pression_mod
150    USE vertical_remap_mod
151    REAL(rstd), INTENT(IN):: tt
152    TYPE(t_field),POINTER :: f_ps(:)
153    TYPE(t_field),POINTER,SAVE :: f_pmid_target(:)
154    TYPE(t_field),POINTER :: f_phis(:)
155    TYPE(t_field),POINTER :: f_theta_rhodz(:)
156    TYPE(t_field),POINTER :: f_u(:) 
157    TYPE(t_field),POINTER :: f_q(:)
158    TYPE(t_field),POINTER, SAVE :: f_T_guided(:),  f_ulon_guided(:), f_ulat_guided(:),f_T_guided_interp(:), &
159                                   f_ulon_guided_interp(:),f_ulat_guided_interp(:)
160    REAL(rstd), POINTER :: target_ue(:,:), ue(:,:), coef_e(:)
161    REAL(rstd), POINTER :: target_theta_rhodz(:,:,:), theta_rhodz(:,:,:), coef_i(:) 
162    INTEGER :: ind
163   
164    IF (abs(MOD(tt,REAL(nudging_time))-dt) < 1.0D-2) THEN
165       
166      CALL allocate_field(f_T_guided, field_t, type_real, llm, name='nudging_T')
167      CALL allocate_field(f_T_guided_interp, field_t, type_real, llm, name='nudging_T')
168      CALL allocate_field(f_ulon_guided_interp, field_t, type_real, llm, name='nudging_T')
169      CALL allocate_field(f_ulat_guided_interp, field_t, type_real, llm, name='nudging_T')
170      CALL allocate_field(f_ulon_guided, field_t, type_real, llm, name='nudging_ulon')
171      CALL allocate_field(f_ulat_guided, field_t, type_real, llm, name='nudging_ulat')
172      CALL allocate_field(f_pmid_target,field_t,type_real,llm,name='target_pressure') 
173      CALL xios_read_field("T_guided_read",f_T_guided)
174      CALL pression_mid(f_ps, f_pmid_target)
175      CALL vertical_remap_extdata(f_T_guided,f_pmid_target,f_T_guided_interp)
176      CALL transfert_request(f_T_guided,req_i0)
177      CALL transfert_request(f_T_guided_interp,req_i0)
178      CALL transfert_request(f_T_guided,req_i1)
179      CALL transfert_request(f_T_guided_interp,req_i1)
180      CALL temperature2theta_rhodz(f_ps,f_T_guided_interp,f_target_theta_rhodz)
181      CALL xios_read_field("ulon_guided_read",f_ulon_guided)
182      CALL xios_read_field("ulat_guided_read",f_ulat_guided)
183      CALL vertical_remap_extdata(f_ulon_guided,f_pmid_target,f_ulon_guided_interp)
184      CALL vertical_remap_extdata(f_ulat_guided,f_pmid_target,f_ulat_guided_interp)
185      CALL transfert_request(f_ulon_guided,req_i0)
186      CALL transfert_request(f_ulon_guided_interp,req_i0)
187      CALL transfert_request(f_ulat_guided,req_i0)
188      CALL transfert_request(f_ulat_guided_interp,req_i0)
189      CALL transfert_request(f_ulon_guided,req_i1)
190      CALL transfert_request(f_ulon_guided_interp,req_i1)
191      CALL transfert_request(f_ulat_guided,req_i1)
192      CALL transfert_request(f_ulat_guided_interp,req_i1)
193      CALL xios_write_field("ulat_guided_out",f_ulat_guided)
194      CALL xios_write_field("ulon_guided_out",f_ulon_guided)
195      CALL xios_write_field("T_guided_out",f_T_guided)
196      CALL ulonlat2un(f_ulon_guided_interp, f_ulat_guided,f_target_ue) 
197      CALL deallocate_field(f_T_guided)
198      CALL deallocate_field(f_T_guided_interp)
199      CALL deallocate_field(f_ulon_guided)
200      CALL deallocate_field(f_ulon_guided_interp)
201      CALL deallocate_field(f_ulat_guided)
202      CALL deallocate_field(f_ulat_guided_interp)
203      CALL deallocate_field(f_pmid_target)
204   
205    ENDIF
206
207
208    DO ind = 1 , ndomain
209       IF (.NOT. assigned_domain(ind)) CYCLE
210       CALL swap_dimensions(ind)
211       CALL swap_geometry(ind) 
212       coef_e = f_relax_coef_e(ind)
213       target_ue = f_target_ue(ind)
214       ue = f_u(ind)
215       CALL compute_guided_u(coef_e, target_ue, ue)
216       coef_i = f_relax_coef_i(ind)
217       target_theta_rhodz = f_target_theta_rhodz(ind)
218       theta_rhodz = f_theta_rhodz(ind)
219       CALL compute_guided_center(coef_i, target_theta_rhodz, theta_rhodz)
220    END DO
221
222  END SUBROUTINE guided
223
224  SUBROUTINE compute_guided_u(coef_e, target_ue, ue)
225    REAL(rstd), INTENT(IN) :: coef_e(iim*3*jjm)
226    REAL(rstd), INTENT(IN) :: target_ue(iim*3*jjm,llm)
227    REAL(rstd), INTENT(INOUT) :: ue(iim*3*jjm,llm)
228    INTEGER :: l, ij
229    DO l = ll_begin, ll_end
230       DO ij=ij_begin_ext, ij_end_ext
231          ue(ij+u_right,l) = ue(ij+u_right,l)*coef_e(ij+u_right) + &
232               target_ue(ij+u_right,l)*(1.-coef_e(ij+u_right))
233          ue(ij+u_lup,l) = ue(ij+u_lup,l)*coef_e(ij+u_lup) + &
234               target_ue(ij+u_lup,l)*(1.-coef_e(ij+u_lup))
235          ue(ij+u_ldown,l) = ue(ij+u_ldown,l)*coef_e(ij+u_ldown) + &
236               target_ue(ij+u_ldown,l)*(1.-coef_e(ij+u_ldown))
237       END DO
238    END DO
239  END SUBROUTINE compute_guided_u
240
241  SUBROUTINE compute_guided_center(coef_i, target_theta_rhodz, theta_rhodz)
242    REAL(rstd), INTENT(IN) :: coef_i(iim*jjm)
243    REAL(rstd), INTENT(IN) :: target_theta_rhodz(iim*jjm,llm,nqdyn)
244    REAL(rstd), INTENT(INOUT) :: theta_rhodz(iim*jjm,llm,nqdyn)
245    INTEGER :: l, ij, iq
246    DO iq=1, nqdyn
247       DO l = ll_begin, ll_end
248          DO ij=ij_begin_ext, ij_end_ext
249             theta_rhodz(ij,l,iq) = theta_rhodz(ij,l,iq)*coef_i(ij) + &
250                  target_theta_rhodz(ij,l,iq)*(1.-coef_i(ij))
251          END DO
252       END DO
253    END DO
254  END SUBROUTINE compute_guided_center
255 
256END MODULE nudging_mod
Note: See TracBrowser for help on using the repository browser.