source: codes/icosagcm/trunk/src/physics/physics_external.F90 @ 871

Last change on this file since 871 was 871, checked in by ymipsl, 5 years ago

experimental : add smooth physics tendency for external physics.

Activated with param :

phys_smooth_tendency=y

To be tested...

YM
File size: 4.9 KB
Line 
1MODULE physics_external_mod
2  USE field_mod
3 
4  INTEGER,SAVE :: it
5!$OMP THREADPRIVATE(it)
6
7  TYPE(t_field),POINTER,SAVE :: f_phis(:)
8  TYPE(t_field),POINTER,SAVE :: f_ps(:)
9  TYPE(t_field),POINTER,SAVE :: f_theta_rhodz(:)
10  TYPE(t_field),POINTER,SAVE :: f_u(:)
11  TYPE(t_field),POINTER,SAVE :: f_wflux(:)
12  TYPE(t_field),POINTER,SAVE :: f_q(:)
13
14  TYPE(t_field),POINTER,SAVE :: f_theta_rhodz0(:)
15  TYPE(t_field),POINTER,SAVE :: f_u0(:)
16  TYPE(t_field),POINTER,SAVE :: f_q0(:)
17 
18  TYPE(t_field),POINTER,SAVE :: f_dtheta_rhodz(:)
19  TYPE(t_field),POINTER,SAVE :: f_du(:)
20  TYPE(t_field),POINTER,SAVE :: f_dq(:)
21
22  TYPE(t_field),POINTER,SAVE :: f_rhodz(:)
23  TYPE(t_field),POINTER,SAVE :: f_rhodz0(:)
24 
25  LOGICAL,SAVE :: phys_smooth_tendency
26!$OMP THREADPRIVATE(phys_smooth_tendency) 
27
28
29CONTAINS
30
31  SUBROUTINE init_physics
32  USE icosa
33  IMPLICIT NONE
34
35    CALL initialize_external_physics
36    CALL allocate_field(f_theta_rhodz0, field_t, type_real, llm, nqdyn, name='theta_rhodz0')
37    CALL allocate_field(f_u0,field_u,type_real,llm,name='u0')
38    CALL allocate_field(f_q0,field_t,type_real,llm,nqtot,'q0')
39
40    CALL allocate_field(f_dtheta_rhodz, field_t, type_real, llm, nqdyn, name='theta_rhodz0')
41    CALL allocate_field(f_du,field_u,type_real,llm,name='u0')
42    CALL allocate_field(f_dq,field_t,type_real,llm,nqtot,'q0')
43
44    CALL allocate_field(f_rhodz, field_t, type_real, llm, name='rhodz')
45   
46    phys_smooth_tendency=.FALSE.
47    CALL getin("phys_smooth_tendency",phys_smooth_tendency)
48   
49       
50  END SUBROUTINE init_physics
51 
52  SUBROUTINE physics(it_,f_phis_, f_ps_, f_theta_rhodz_, f_u_, f_wflux_, f_q_)
53  USE icosa
54  USE field_mod
55  USE mpipara
56  USE omp_para
57  USE xios
58  USE domain_mod
59  USE time_mod
60  USE disvert_mod
61  IMPLICIT NONE
62    INTEGER,INTENT(IN)    :: it_
63    TYPE(t_field),POINTER :: f_phis_(:)
64    TYPE(t_field),POINTER :: f_ps_(:)
65    TYPE(t_field),POINTER :: f_theta_rhodz_(:)
66    TYPE(t_field),POINTER :: f_u_(:)
67    TYPE(t_field),POINTER :: f_wflux_(:)
68    TYPE(t_field),POINTER :: f_q_(:)
69
70    REAL(rstd),POINTER    :: theta_rhodz(:,:,:), theta_rhodz0(:,:,:), dtheta_rhodz(:,:,:) 
71    REAL(rstd),POINTER    :: u(:,:), u0(:,:), du(:,:)
72    REAL(rstd),POINTER    :: q(:,:,:),q0(:,:,:),dq(:,:,:)
73    REAL(rstd),POINTER    :: ps(:)
74    REAL(rstd),POINTER    :: rhodz(:,:)
75    INTEGER :: ind, iq
76   
77   
78!$OMP BARRIER
79!$OMP MASTER
80    f_phis=>f_phis_
81    f_ps=>f_ps_
82    f_theta_rhodz=>f_theta_rhodz_
83    f_u=>f_u_
84    f_wflux=>f_wflux_
85    f_q=>f_q_
86!$OMP END MASTER
87!$OMP BARRIER
88
89    IF (phys_smooth_tendency) THEN
90   
91      IF (MOD(it_,itau_physics)==1) THEN
92        DO ind=1, ndomain
93          IF (.NOT. assigned_domain(ind)) CYCLE
94          CALL swap_dimensions(ind)
95          CALL swap_geometry(ind)
96          theta_rhodz=f_theta_rhodz(ind)
97          theta_rhodz0=f_theta_rhodz0(ind)
98          u=f_u(ind)
99          u0=f_u0(ind)
100          q=f_q(ind)
101          q0=f_q0(ind)
102          ps=f_ps(ind)
103          rhodz=f_rhodz(ind)
104       
105          theta_rhodz0(:,:,1)=theta_rhodz(:,:,1)
106          u0=u
107          q0=q
108          CALL compute_rhodz(.TRUE., ps, rhodz)
109        ENDDO
110     
111        IF (is_omp_master) CALL xios_timer_suspend("dynamico")
112        it = it_-1 + itau_physics
113        CALL external_physics
114        IF (is_omp_master) CALL xios_timer_resume("dynamico")
115
116        DO ind=1, ndomain
117          IF (.NOT. assigned_domain(ind)) CYCLE
118          CALL swap_dimensions(ind) 
119          CALL swap_geometry(ind) 
120          theta_rhodz=f_theta_rhodz(ind)
121          theta_rhodz0=f_theta_rhodz0(ind)
122          u=f_u(ind)
123          u0=f_u0(ind)
124          q=f_q(ind)
125          q0=f_q0(ind)
126          dtheta_rhodz=f_dtheta_rhodz(ind)
127          du=f_du(ind)
128          dq=f_dq(ind)
129          rhodz=f_rhodz(ind)
130       
131          dtheta_rhodz(:,:,1)=(theta_rhodz(:,:,1)-theta_rhodz0(:,:,1))/itau_physics
132          du=(u-u0)/itau_physics
133   
134          DO iq=1, nqtot
135            dq(:,:,iq)=((q(:,:,iq)-q0(:,:,iq))/itau_physics)/rhodz(:,:)
136          ENDDO
137         
138          theta_rhodz(:,:,1)=theta_rhodz0(:,:,1)
139          u=u0
140          q=q0
141        ENDDO
142     ENDIF
143   
144     DO ind=1, ndomain
145       IF (.NOT. assigned_domain(ind)) CYCLE
146       CALL swap_dimensions(ind)
147       CALL swap_geometry(ind)
148         
149       theta_rhodz=f_theta_rhodz(ind)
150       u=f_u(ind)
151       q=f_q(ind)
152       dtheta_rhodz=f_dtheta_rhodz(ind)
153       du=f_du(ind)
154       dq=f_dq(ind)
155       rhodz=f_rhodz(ind)
156       ps=f_ps(ind)
157
158       u=u+du
159       theta_rhodz=theta_rhodz+dtheta_rhodz
160       CALL compute_rhodz(.TRUE., ps, rhodz)
161       DO iq=1, nqtot
162         q(:,:,iq)=q(:,:,iq)+dq(:,:,iq)*rhodz(:,:)
163       ENDDO
164     ENDDO
165   
166   ELSE
167       
168     IF (MOD(it_,itau_physics)==0) THEN
169       it=it_
170       IF (is_omp_master) CALL xios_timer_suspend("dynamico")
171       CALL external_physics
172       IF (is_omp_master) CALL xios_timer_resume("dynamico")
173     ENDIF
174   
175   ENDIF     
176
177
178  END SUBROUTINE physics
179 
180 
181END MODULE physics_external_mod
182   
183 
Note: See TracBrowser for help on using the repository browser.