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

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

devel: Fix to Changeset 880
File size: 9.3 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
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    REAL(rstd), INTENT(IN):: tt
151    TYPE(t_field),POINTER :: f_ps(:)
152    TYPE(t_field),POINTER :: f_phis(:)
153    TYPE(t_field),POINTER :: f_theta_rhodz(:)
154    TYPE(t_field),POINTER :: f_u(:) 
155    TYPE(t_field),POINTER :: f_q(:)
156    TYPE(t_field),POINTER, SAVE :: f_T_guided(:),  f_ulon_guided(:), f_ulat_guided(:),f_pressure_mid(:) 
157    REAL(rstd), POINTER :: target_ue(:,:), ue(:,:), coef_e(:)
158    REAL(rstd), POINTER :: target_theta_rhodz(:,:,:), theta_rhodz(:,:,:), coef_i(:) 
159    INTEGER :: ind
160
161    IF (MOD(INT(tt),nudging_time)==dt) THEN
162   
163      CALL allocate_field(f_T_guided, field_t, type_real, llm, name='nudging_T')
164      CALL allocate_field(f_ulon_guided, field_t, type_real, llm, name='nudging_ulon')
165      CALL allocate_field(f_ulat_guided, field_t, type_real, llm, name='nudging_ulat')
166
167      CALL xios_read_field("T_guided_read",f_T_guided)
168      CALL transfert_request(f_T_guided,req_i0)
169      CALL temperature2theta_rhodz(f_ps,f_T_guided,f_target_theta_rhodz)
170      CALL xios_read_field("ulon_guided_read",f_ulon_guided)
171      CALL xios_read_field("ulat_guided_read",f_ulat_guided)
172      CALL transfert_request(f_ulon_guided,req_i0)
173      CALL transfert_request(f_ulat_guided,req_i0)
174      CALL transfert_request(f_ulon_guided,req_i1)
175      CALL transfert_request(f_ulat_guided,req_i1)
176      CALL xios_write_field("ulat_guided_out",f_ulat_guided)
177      CALL xios_write_field("ulon_guided_out",f_ulon_guided)
178      CALL ulonlat2un(f_ulon_guided, f_ulat_guided,f_target_ue)
179
180      CALL deallocate_field(f_T_guided)
181      CALL deallocate_field(f_ulon_guided)
182      CALL deallocate_field(f_ulat_guided)
183   
184    ENDIF
185
186
187    DO ind = 1 , ndomain
188       IF (.NOT. assigned_domain(ind)) CYCLE
189       CALL swap_dimensions(ind)
190       CALL swap_geometry(ind) 
191       coef_e = f_relax_coef_e(ind)
192       target_ue = f_target_ue(ind)
193       ue = f_u(ind)
194       CALL compute_guided_u(coef_e, target_ue, ue)
195       coef_i = f_relax_coef_i(ind)
196       target_theta_rhodz = f_target_theta_rhodz(ind)
197       theta_rhodz = f_theta_rhodz(ind)
198       CALL compute_guided_center(coef_i, target_theta_rhodz, theta_rhodz)
199    END DO
200
201  END SUBROUTINE guided
202
203  SUBROUTINE compute_guided_u(coef_e, target_ue, ue)
204    REAL(rstd), INTENT(IN) :: coef_e(iim*3*jjm)
205    REAL(rstd), INTENT(IN) :: target_ue(iim*3*jjm,llm)
206    REAL(rstd), INTENT(INOUT) :: ue(iim*3*jjm,llm)
207    INTEGER :: l, ij
208    DO l = ll_begin, ll_end
209       DO ij=ij_begin_ext, ij_end_ext
210          ue(ij+u_right,l) = ue(ij+u_right,l)*coef_e(ij+u_right) + &
211               target_ue(ij+u_right,l)*(1.-coef_e(ij+u_right))
212          ue(ij+u_lup,l) = ue(ij+u_lup,l)*coef_e(ij+u_lup) + &
213               target_ue(ij+u_lup,l)*(1.-coef_e(ij+u_lup))
214          ue(ij+u_ldown,l) = ue(ij+u_ldown,l)*coef_e(ij+u_ldown) + &
215               target_ue(ij+u_ldown,l)*(1.-coef_e(ij+u_ldown))
216       END DO
217    END DO
218  END SUBROUTINE compute_guided_u
219
220  SUBROUTINE compute_guided_center(coef_i, target_theta_rhodz, theta_rhodz)
221    REAL(rstd), INTENT(IN) :: coef_i(iim*jjm)
222    REAL(rstd), INTENT(IN) :: target_theta_rhodz(iim*jjm,llm,nqdyn)
223    REAL(rstd), INTENT(INOUT) :: theta_rhodz(iim*jjm,llm,nqdyn)
224    INTEGER :: l, ij, iq
225    DO iq=1, nqdyn
226       DO l = ll_begin, ll_end
227          DO ij=ij_begin_ext, ij_end_ext
228             theta_rhodz(ij,l,iq) = theta_rhodz(ij,l,iq)*coef_i(ij) + &
229                  target_theta_rhodz(ij,l,iq)*(1.-coef_i(ij))
230          END DO
231       END DO
232    END DO
233  END SUBROUTINE compute_guided_center
234 
235END MODULE nudging_mod
Note: See TracBrowser for help on using the repository browser.