source: codes/icosagcm/devel/src/diagnostics/observable.f90 @ 600

Last change on this file since 600 was 600, checked in by dubos, 7 years ago

devel : remove unused buffers in diagnostics/observable.f90
File size: 13.5 KB
Line 
1MODULE observable_mod
2  USE icosa
3  IMPLICIT NONE
4  PRIVATE
5
6  TYPE(t_field),POINTER, SAVE :: f_buf_i(:), &
7       f_buf_Fel(:), f_buf_uh(:), & ! horizontal velocity, different from prognostic velocity if NH
8       f_buf_ulon(:), f_buf_ulat(:)
9  TYPE(t_field),POINTER, SAVE :: f_buf_v(:), f_buf_s(:), f_buf_p(:)
10  TYPE(t_field),POINTER, SAVE :: f_pmid(:)
11
12! temporary shared variable for caldyn
13  TYPE(t_field),POINTER, SAVE :: f_theta(:)
14
15  PUBLIC init_observable, write_output_fields_basic, f_theta, f_buf_i
16
17CONTAINS
18 
19  SUBROUTINE init_observable
20    CALL allocate_field(f_buf_i,   field_t,type_real,llm,name="buffer_i")
21    CALL allocate_field(f_buf_p,   field_t,type_real,llm+1) 
22    CALL allocate_field(f_buf_ulon,field_t,type_real,llm, name="buf_ulon")
23    CALL allocate_field(f_buf_ulat,field_t,type_real,llm, name="buf_ulat")
24    CALL allocate_field(f_buf_uh,  field_u,type_real,llm, name="buf_uh")
25    CALL allocate_field(f_buf_Fel, field_u,type_real,llm+1, name="buf_F_el")
26    CALL allocate_field(f_buf_v,   field_z,type_real,llm, name="buf_v")
27    CALL allocate_field(f_buf_s,   field_t,type_real, name="buf_s")
28
29    CALL allocate_field(f_theta, field_t,type_real,llm,nqdyn,  name='theta') ! potential temperature
30    CALL allocate_field(f_pmid,  field_t,type_real,llm,  name='pmid')       ! mid layer pressure
31  END SUBROUTINE init_observable
32
33  SUBROUTINE write_output_fields_basic(init, f_phis, f_ps, f_mass, f_geopot, f_theta_rhodz, f_u, f_W, f_q)
34    USE xios_mod
35    USE disvert_mod
36    USE wind_mod
37    USE output_field_mod
38    USE omp_para
39    USE time_mod
40    USE xios_mod
41    USE earth_const
42    USE pression_mod
43    USE vertical_interp_mod
44    USE theta2theta_rhodz_mod
45    USE omega_mod
46    USE diagflux_mod
47    LOGICAL, INTENT(IN) :: init
48    INTEGER :: l
49    REAL :: scalar(1)
50    REAL :: mid_ap(llm)
51    REAL :: mid_bp(llm)
52
53    TYPE(t_field),POINTER :: f_phis(:), f_ps(:), f_mass(:), f_geopot(:), f_theta_rhodz(:), f_u(:), f_W(:), f_q(:)
54!    IF (is_master) PRINT *,'CALL write_output_fields_basic'
55
56    CALL transfert_request(f_ps,req_i1)
57   
58    IF(init) THEN
59       IF(is_master) PRINT *, 'Creating output files ...'
60       scalar(1)=dt
61       IF (is_omp_master) CALL xios_send_field("timestep", scalar)
62       scalar(1)=preff
63       IF (is_omp_master) CALL xios_send_field("preff", scalar)
64       IF (is_omp_master) CALL xios_send_field("ap",ap)
65       IF (is_omp_master) CALL xios_send_field("bp",bp)
66       DO l=1,llm
67          mid_ap(l)=(ap(l)+ap(l+1))/2
68          mid_bp(l)=(bp(l)+bp(l+1))/2
69       ENDDO
70       IF (is_omp_master) CALL xios_send_field("mid_ap",mid_ap)
71       IF (is_omp_master) CALL xios_send_field("mid_bp",mid_bp)
72
73       CALL output_field("phis",f_phis)
74       CALL output_field("Ai",geom%Ai) 
75       IF(is_master) PRINT *, '... done creating output files. Writing initial condition ...'
76    END IF
77
78    IF(nqdyn>1) THEN
79       CALL divide_by_mass(2, f_mass, f_theta_rhodz, f_buf_i)
80       IF(init) THEN
81          CALL output_field("dyn_q_init",f_buf_i)
82       ELSE
83          CALL output_field("dyn_q",f_buf_i)
84       END IF
85    END IF
86
87    CALL divide_by_mass(1, f_mass, f_theta_rhodz, f_buf_i)
88    IF(init) THEN
89       CALL output_field("theta_init",f_buf_i)
90    ELSE
91       CALL output_field("theta",f_buf_i)
92    END IF
93
94    CALL pression_mid(f_ps, f_pmid)
95    CALL diagnose_temperature(f_pmid, f_q, f_buf_i) ! f_buf_i : IN = theta, out = T
96
97    IF(init) THEN
98       CALL output_field("temp_init",f_buf_i)
99    ELSE
100       CALL output_field("temp",f_buf_i)
101       CALL vertical_interp(f_pmid,f_buf_i,f_buf_s,85000.)
102       CALL output_field("t850",f_buf_s)
103       CALL vertical_interp(f_pmid,f_buf_i,f_buf_s,50000.)
104       CALL output_field("t500",f_buf_s)
105       CALL vertical_interp(f_pmid,f_buf_i,f_buf_s,preff)
106       CALL output_field("SST",f_buf_s)       
107    END IF
108   
109    CALL progonostic_vel_to_horiz(f_geopot, f_ps, f_mass, f_u, f_W, f_buf_uh, f_buf_i)
110    CALL transfert_request(f_buf_uh,req_e1_vect) 
111    CALL un2ulonlat(f_buf_uh, f_buf_ulon, f_buf_ulat)
112    IF(init) THEN
113       CALL output_field("uz_init",f_buf_i)
114       CALL output_field("ulon_init",f_buf_ulon)
115       CALL output_field("ulat_init",f_buf_ulat)
116       CALL output_field("p_init",f_pmid)
117       CALL output_field("ps_init",f_ps)
118       CALL output_field("mass_init",f_mass)
119       CALL output_field("geopot_init",f_geopot)
120       CALL output_field("q_init",f_q)
121       IF(is_master) PRINT *, 'Done writing initial condition ...'
122    ELSE
123       CALL output_field("uz",f_buf_i)
124       CALL output_field("ulon",f_buf_ulon)
125       CALL output_field("ulat",f_buf_ulat)
126       CALL output_field("p",f_pmid)
127       CALL output_field("ps",f_ps)
128       CALL output_field("mass",f_mass)
129       CALL output_field("geopot",f_geopot)
130       CALL output_field("q",f_q)
131
132       !       CALL output_field("exner",f_pk)
133       !       CALL output_field("pv",f_qv)
134       
135       CALL vertical_interp(f_pmid,f_buf_ulon,f_buf_s,85000.)
136       CALL output_field("u850",f_buf_s)
137       CALL vertical_interp(f_pmid,f_buf_ulon,f_buf_s,50000.)
138       CALL output_field("u500",f_buf_s)
139       
140       CALL vertical_interp(f_pmid,f_buf_ulat,f_buf_s,85000.)
141       CALL output_field("v850",f_buf_s)
142       CALL vertical_interp(f_pmid,f_buf_ulat,f_buf_s,50000.)
143       CALL output_field("v500",f_buf_s)
144
145       CALL vertical_interp(f_pmid,f_buf_i,f_buf_s,85000.)
146       CALL output_field("w850",f_buf_s)
147       CALL vertical_interp(f_pmid,f_buf_i,f_buf_s,50000.)
148       CALL output_field("w500",f_buf_s)   
149
150       CALL w_omega(f_ps, f_u, f_buf_i)
151       CALL output_field("omega",f_buf_i)
152       CALL vertical_interp(f_pmid,f_buf_i,f_buf_s,85000.)
153       CALL output_field("omega850",f_buf_s)
154       CALL vertical_interp(f_pmid,f_buf_i,f_buf_s,50000.)
155       CALL output_field("omega500",f_buf_s)
156    END IF
157
158    IF(.NOT. init) THEN
159       IF(diagflux_on) THEN
160          CALL flux_centered_lonlat(1./(itau_out*dt) , f_massfluxt, f_buf_ulon, f_buf_ulat)
161          CALL output_field("mass_t", f_masst)
162          CALL output_field("massflux_lon",f_buf_ulon)
163          CALL output_field("massflux_lat",f_buf_ulat)
164
165          CALL transfert_request(f_epotfluxt,req_e1_vect) 
166          CALL flux_centered_lonlat(1./(itau_out*dt) , f_epotfluxt, f_buf_ulon, f_buf_ulat)
167          CALL output_field("epot_t", f_epot)
168          CALL output_field("epotflux_lon",f_buf_ulon)
169          CALL output_field("epotflux_lat",f_buf_ulat)
170
171          CALL transfert_request(f_ekinfluxt,req_e1_vect) 
172          CALL flux_centered_lonlat(1./(itau_out*dt) , f_ekinfluxt, f_buf_ulon, f_buf_ulat)
173          CALL output_field("ekin_t", f_ekin)
174          CALL output_field("ekinflux_lon",f_buf_ulon)
175          CALL output_field("ekinflux_lat",f_buf_ulat)
176
177          CALL transfert_request(f_enthalpyfluxt,req_e1_vect) 
178          CALL flux_centered_lonlat(1./(itau_out*dt) , f_enthalpyfluxt, f_buf_ulon, f_buf_ulat)
179          CALL output_field("enthalpy_t", f_enthalpy)
180          CALL output_field("enthalpyflux_lon",f_buf_ulon)
181          CALL output_field("enthalpyflux_lat",f_buf_ulat)
182
183          CALL qflux_centered_lonlat(1./(itau_out*dt) , f_qfluxt, f_qfluxt_lon, f_qfluxt_lat)
184          CALL output_field("qmass_t", f_qmasst)
185          CALL output_field("qflux_lon",f_qfluxt_lon)
186          CALL output_field("qflux_lat",f_qfluxt_lat)
187          CALL zero_qfluxt  ! restart accumulating fluxes
188       END IF
189    END IF
190  END SUBROUTINE write_output_fields_basic
191 
192 !------------------- Conversion from prognostic to observable variables ------------------
193
194  SUBROUTINE progonostic_vel_to_horiz(f_geopot, f_ps, f_rhodz, f_u, f_W, f_uh, f_uz)
195    USE disvert_mod
196    TYPE(t_field), POINTER :: f_geopot(:), f_ps(:), f_rhodz(:), &
197         f_u(:), f_W(:), f_uz(:), &  ! IN
198         f_uh(:)                         ! OUT
199    REAL(rstd),POINTER :: geopot(:,:), ps(:), rhodz(:,:), u(:,:), W(:,:), uh(:,:), uz(:,:), F_el(:,:)
200    INTEGER :: ind
201   
202    DO ind=1,ndomain
203       IF (.NOT. assigned_domain(ind)) CYCLE
204       CALL swap_dimensions(ind)
205       CALL swap_geometry(ind)
206       geopot = f_geopot(ind)
207       rhodz = f_rhodz(ind)
208       u = f_u(ind)
209       W = f_W(ind)
210       uh  = f_uh(ind)
211       F_el  = f_buf_Fel(ind)
212       IF(caldyn_eta==eta_mass) THEN
213          ps=f_ps(ind)
214          CALL compute_rhodz(.TRUE., ps, rhodz)
215       END IF
216       uz = f_uz(ind)
217       !$OMP BARRIER
218       CALL compute_prognostic_vel_to_horiz(geopot,rhodz,u,W, F_el, uh,uz)
219       !$OMP BARRIER
220    END DO
221  END SUBROUTINE progonostic_vel_to_horiz
222 
223  SUBROUTINE compute_prognostic_vel_to_horiz(Phi, rhodz, u, W, F_el, uh, uz)
224    USE omp_para
225    REAL(rstd), INTENT(IN) :: Phi(iim*jjm,llm+1)
226    REAL(rstd), INTENT(IN) :: rhodz(iim*jjm,llm)
227    REAL(rstd), INTENT(IN) :: u(3*iim*jjm,llm)
228    REAL(rstd), INTENT(IN) :: W(iim*jjm,llm+1)
229    REAL(rstd), INTENT(OUT) :: uh(3*iim*jjm,llm)
230    REAL(rstd), INTENT(OUT) :: uz(iim*jjm,llm)
231    INTEGER :: ij,l
232    REAL(rstd) :: F_el(3*iim*jjm,llm+1)
233    REAL(rstd) :: uu_right, uu_lup, uu_ldown, W_el, DePhil
234    ! NB : u and uh are not in DEC form, they are normal components   
235    ! => we must divide by de
236    IF(hydrostatic) THEN
237       uh(:,:)=u(:,:)
238       uz(:,:)=0.
239    ELSE
240       DO l=ll_begin, ll_endp1 ! compute on l levels (interfaces)
241          DO ij=ij_begin_ext, ij_end_ext
242             ! Compute on edge 'right'
243             W_el   = .5*( W(ij,l)+W(ij+t_right,l) )
244             DePhil = ne_right*(Phi(ij+t_right,l)-Phi(ij,l))
245             F_el(ij+u_right,l)   = DePhil*W_el/de(ij+u_right)
246             ! Compute on edge 'lup'
247             W_el   = .5*( W(ij,l)+W(ij+t_lup,l) )
248             DePhil = ne_lup*(Phi(ij+t_lup,l)-Phi(ij,l))
249             F_el(ij+u_lup,l)   = DePhil*W_el/de(ij+u_lup)
250             ! Compute on edge 'ldown'
251             W_el   = .5*( W(ij,l)+W(ij+t_ldown,l) )
252             DePhil = ne_ldown*(Phi(ij+t_ldown,l)-Phi(ij,l))
253             F_el(ij+u_ldown,l) = DePhil*W_el/de(ij+u_ldown)
254          END DO
255       END DO
256
257       ! We need a barrier here because we compute F_el above and do a vertical average below
258       !$OMP BARRIER
259
260       DO l=ll_begin, ll_end ! compute on k levels (full levels)
261          DO ij=ij_begin_ext, ij_end_ext
262             ! w = vertical momentum = g^-2*dPhi/dt = uz/g
263             uz(ij,l) = (.5*g)*(W(ij,l)+W(ij,l+1))/rhodz(ij,l)
264             ! uh = u-w.grad(Phi) = u - uz.grad(z)
265             uh(ij+u_right,l) = u(ij+u_right,l) - (F_el(ij+u_right,l)+F_el(ij+u_right,l+1)) / (rhodz(ij,l)+rhodz(ij+t_right,l))
266             uh(ij+u_lup,l)   = u(ij+u_lup,l)   - (F_el(ij+u_lup,l)+F_el(ij+u_lup,l+1))     / (rhodz(ij,l)+rhodz(ij+t_lup,l))
267             uh(ij+u_ldown,l) = u(ij+u_ldown,l) - (F_el(ij+u_ldown,l)+F_el(ij+u_ldown,l+1)) / (rhodz(ij,l)+rhodz(ij+t_ldown,l))
268          END DO
269       END DO
270
271    END IF
272  END SUBROUTINE compute_prognostic_vel_to_horiz
273
274  SUBROUTINE diagnose_temperature(f_pmid,f_q,f_temp)
275    USE icosa
276    USE pression_mod
277    IMPLICIT NONE
278    TYPE(t_field), POINTER :: f_pmid(:)           ! IN
279    TYPE(t_field), POINTER :: f_q(:)            ! IN
280    TYPE(t_field), POINTER :: f_temp(:)         ! INOUT
281   
282    REAL(rstd), POINTER :: pmid(:,:)
283    REAL(rstd), POINTER :: q(:,:,:)
284    REAL(rstd), POINTER :: temp(:,:)
285    INTEGER :: ind
286   
287    DO ind=1,ndomain
288       IF (.NOT. assigned_domain(ind)) CYCLE
289       CALL swap_dimensions(ind)
290       CALL swap_geometry(ind)
291       pmid=f_pmid(ind)
292       q=f_q(ind)
293       temp=f_temp(ind)
294       CALL compute_diagnose_temp(pmid,q,temp)
295    END DO
296  END SUBROUTINE diagnose_temperature
297 
298  SUBROUTINE compute_diagnose_temp(pmid,q,temp)
299    USE omp_para
300    USE pression_mod
301    REAL(rstd),INTENT(IN)    :: pmid(iim*jjm,llm)
302    REAL(rstd),INTENT(IN)    :: q(iim*jjm,llm,nqtot)
303    REAL(rstd),INTENT(INOUT) :: temp(iim*jjm,llm)
304
305    REAL(rstd) :: Rd, p_ik, theta_ik, temp_ik, qv, chi, Rmix
306    INTEGER :: ij,l
307
308    Rd = kappa*cpp
309    DO l=ll_begin,ll_end
310       DO ij=ij_begin,ij_end
311          p_ik = pmid(ij,l)
312          theta_ik = temp(ij,l)
313          qv = q(ij,l,1) ! water vaper mixing ratio = mv/md
314          SELECT CASE(caldyn_thermo)
315          CASE(thermo_theta)
316             temp_ik = theta_ik*((p_ik/preff)**kappa)
317          CASE(thermo_entropy)
318             temp_ik = Treff*exp((theta_ik + Rd*log(p_ik/preff))/cpp)
319          CASE(thermo_moist)
320             Rmix = Rd+qv*Rv
321             chi = ( theta_ik + Rmix*log(p_ik/preff) ) / (cpp + qv*cppv) ! log(T/Treff)
322             temp_ik = Treff*exp(chi)
323          END SELECT
324          IF(physics_thermo==thermo_fake_moist) temp_ik=temp_ik/(1+0.608*qv)
325          temp(ij,l)=temp_ik
326       END DO
327    END DO
328  END SUBROUTINE compute_diagnose_temp
329
330  SUBROUTINE Tv2T(f_Tv, f_q, f_T)
331    TYPE(t_field), POINTER :: f_TV(:)
332    TYPE(t_field), POINTER :: f_q(:)
333    TYPE(t_field), POINTER :: f_T(:)
334   
335    REAL(rstd),POINTER :: Tv(:,:), q(:,:,:), T(:,:)
336    INTEGER :: ind
337   
338    DO ind=1,ndomain
339       IF (.NOT. assigned_domain(ind)) CYCLE
340       CALL swap_dimensions(ind)
341       CALL swap_geometry(ind)
342       Tv=f_Tv(ind)
343       T=f_T(ind)
344       SELECT CASE(physics_thermo)
345       CASE(thermo_dry)
346          T=Tv
347       CASE(thermo_fake_moist)
348          q=f_q(ind)
349          T=Tv/(1+0.608*q(:,:,1))
350       END SELECT
351    END DO
352  END SUBROUTINE Tv2T
353
354  SUBROUTINE divide_by_mass(iq, f_mass, f_theta_rhodz, f_theta)
355    INTEGER, INTENT(IN) :: iq
356    TYPE(t_field), POINTER :: f_mass(:), f_theta_rhodz(:), f_theta(:)
357    REAL(rstd), POINTER :: mass(:,:), theta_rhodz(:,:,:), theta(:,:)
358    INTEGER :: ind
359    DO ind=1,ndomain
360       IF (.NOT. assigned_domain(ind)) CYCLE
361       CALL swap_dimensions(ind)
362       CALL swap_geometry(ind)
363       mass=f_mass(ind)
364       theta_rhodz=f_theta_rhodz(ind)
365       theta=f_theta(ind)
366       theta(:,:) = theta_rhodz(:,:,iq) / mass(:,:)
367    END DO
368  END SUBROUTINE divide_by_mass
369   
370END MODULE observable_mod
Note: See TracBrowser for help on using the repository browser.