source: codes/icosagcm/trunk/src/observable.f90 @ 417

Last change on this file since 417 was 417, checked in by dubos, 8 years ago

Fixed computation of omega
File size: 13.7 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_uh(:), & ! horizontal velocity, different from prognostic velocity if NH
8       f_buf_ulon(:), f_buf_ulat(:), &
9       f_buf_u3d(:) ! unused, remove ?
10  TYPE(t_field),POINTER, SAVE :: f_buf1_i(:), f_buf2_i(:)
11  TYPE(t_field),POINTER, SAVE :: f_buf_v(:), f_buf_s(:), f_buf_p(:)
12  TYPE(t_field),POINTER, SAVE :: f_pmid(:)
13
14! temporary shared variable for caldyn
15  TYPE(t_field),POINTER, SAVE :: f_theta(:)
16
17  PUBLIC init_observable, write_output_fields_basic, f_theta
18
19!$OMP THREADPRIVATE(first_output)
20 
21CONTAINS
22 
23  SUBROUTINE init_observable
24    CALL allocate_field(f_buf_i,   field_t,type_real,llm,name="buffer_i")
25    CALL allocate_field(f_buf1_i,   field_t,type_real,llm,name="buffer1_i")
26    CALL allocate_field(f_buf2_i,   field_t,type_real,llm,name="buffer2_i")
27    CALL allocate_field(f_buf_p,   field_t,type_real,llm+1) 
28    CALL allocate_field(f_buf_u3d, field_t,type_real,3,llm)  ! 3D vel at cell centers
29    CALL allocate_field(f_buf_ulon,field_t,type_real,llm, name="buf_ulon")
30    CALL allocate_field(f_buf_ulat,field_t,type_real,llm, name="buf_ulat")
31    CALL allocate_field(f_buf_uh,  field_u,type_real,llm, name="buf_uh")
32    CALL allocate_field(f_buf_v,   field_z,type_real,llm, name="buf_v")
33    CALL allocate_field(f_buf_s,   field_t,type_real, name="buf_s")
34
35    CALL allocate_field(f_theta, field_t,type_real,llm,nqdyn,  name='theta') ! potential temperature
36    CALL allocate_field(f_pmid,  field_t,type_real,llm,  name='pmid')       ! mid layer pressure
37  END SUBROUTINE init_observable
38
39  SUBROUTINE write_output_fields_basic(init, f_phis, f_ps, f_mass, f_geopot, f_theta_rhodz, f_u, f_W, f_q)
40    USE xios
41    USE disvert_mod
42    USE wind_mod
43    USE output_field_mod
44    USE omp_para
45    USE time_mod
46    USE xios
47    USE earth_const
48    USE pression_mod
49    USE vertical_interp_mod
50    USE theta2theta_rhodz_mod
51    USE omega_mod
52    LOGICAL, INTENT(IN) :: init
53    INTEGER :: l
54    REAL :: scalar(1)
55    REAL :: mid_ap(llm)
56    REAL :: mid_bp(llm)
57
58    TYPE(t_field),POINTER :: f_phis(:), f_ps(:), f_mass(:), f_geopot(:), f_theta_rhodz(:), f_u(:), f_W(:), f_q(:)
59!    IF (is_master) PRINT *,'CALL write_output_fields_basic'
60
61    CALL transfert_request(f_ps,req_i1)
62   
63    IF(init) THEN
64       scalar(1)=dt
65       CALL xios_send_field("timestep", scalar)
66       scalar(1)=preff
67       CALL xios_send_field("preff", scalar)
68       CALL xios_send_field("ap",ap)
69       CALL xios_send_field("bp",bp)
70       DO l=1,llm
71          mid_ap(l)=(ap(l)+ap(l+1))/2
72          mid_bp(l)=(bp(l)+bp(l+1))/2
73       ENDDO
74       CALL xios_send_field("mid_ap",mid_ap)
75       CALL xios_send_field("mid_bp",mid_bp)
76
77       CALL output_field("phis",f_phis)
78       CALL output_field("Ai",geom%Ai)       
79    END IF
80
81    CALL divide_by_mass(1, f_mass, f_theta_rhodz, f_buf_i)
82    IF(init) THEN
83       CALL output_field("theta_init",f_buf_i)
84    ELSE
85       CALL output_field("theta",f_buf_i)
86    END IF
87
88    IF(nqdyn>1) THEN
89       CALL divide_by_mass(2, f_mass, f_theta_rhodz, f_buf_i)
90       IF(init) THEN
91          CALL output_field("dyn_q_init",f_buf_i)
92       ELSE
93          CALL output_field("dyn_q",f_buf_i)
94       END IF
95    END IF
96
97    CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_buf_i) ; 
98    CALL Tv2T(f_buf_i,f_q,f_buf1_i)
99    IF(init) THEN
100       CALL output_field("temp_init",f_buf_i)
101    ELSE
102       CALL output_field("temp",f_buf_i)
103       CALL vertical_interp(f_ps,f_buf_i,f_buf_s,85000.)
104       CALL output_field("t850",f_buf_s)
105       CALL vertical_interp(f_ps,f_buf_i,f_buf_s,50000.)
106       CALL output_field("t500",f_buf_s)
107       CALL vertical_interp(f_ps,f_buf_i,f_buf_s,preff)
108       CALL output_field("SST",f_buf_s)       
109    END IF
110   
111    CALL progonostic_vel_to_horiz(f_geopot, f_ps, f_mass, f_u, f_W, f_buf_uh, f_buf_i)
112    CALL transfert_request(f_buf_uh,req_e1_vect) 
113    CALL un2ulonlat(f_buf_uh, f_buf_ulon, f_buf_ulat)
114    CALL pression_mid(f_ps, f_pmid)
115    IF(init) THEN
116       CALL output_field("uz_init",f_buf_i)
117       CALL output_field("ulon_init",f_buf_ulon)
118       CALL output_field("ulat_init",f_buf_ulat)
119       CALL output_field("p_init",f_pmid)
120       CALL output_field("ps_init",f_ps)
121       CALL output_field("mass_init",f_mass)
122       CALL output_field("geopot_init",f_geopot)
123       CALL output_field("q_init",f_q)
124    ELSE
125       CALL output_field("uz",f_buf_i)
126       CALL output_field("ulon",f_buf_ulon)
127       CALL output_field("ulat",f_buf_ulat)
128       CALL output_field("p",f_pmid)
129       CALL output_field("ps",f_ps)
130       CALL output_field("mass",f_mass)
131       CALL output_field("geopot",f_geopot)
132       CALL output_field("q",f_q)
133
134       !       CALL output_field("exner",f_pk)
135       !       CALL output_field("pv",f_qv)
136       
137       CALL vertical_interp(f_ps,f_buf_ulon,f_buf_s,85000.)
138       CALL output_field("u850",f_buf_s)
139       CALL vertical_interp(f_ps,f_buf_ulon,f_buf_s,50000.)
140       CALL output_field("u500",f_buf_s)
141       
142       CALL vertical_interp(f_ps,f_buf_ulat,f_buf_s,85000.)
143       CALL output_field("v850",f_buf_s)
144       CALL vertical_interp(f_ps,f_buf_ulat,f_buf_s,50000.)
145       CALL output_field("v500",f_buf_s)
146
147       CALL vertical_interp(f_ps,f_buf_i,f_buf_s,85000.)
148       CALL output_field("w850",f_buf_s)
149       CALL vertical_interp(f_ps,f_buf_i,f_buf_s,50000.)
150       CALL output_field("w500",f_buf_s)   
151
152       CALL w_omega(f_ps, f_u, f_buf_i)
153       CALL output_field("omega",f_buf_i)
154       CALL vertical_interp(f_ps,f_buf_i,f_buf_s,85000.)
155       CALL output_field("omega850",f_buf_s)
156       CALL vertical_interp(f_ps,f_buf_i,f_buf_s,50000.)
157       CALL output_field("omega500",f_buf_s)
158    END IF
159
160  END SUBROUTINE write_output_fields_basic
161 
162  SUBROUTINE write_output_fields(f_ps, f_phis, f_dps, f_u, f_theta_rhodz, f_q, &
163       f_buf_i, f_buf_v, f_buf_i3, f_buf1_i, f_buf2_i, f_buf_s, f_buf_p)
164    USE vorticity_mod
165    USE theta2theta_rhodz_mod
166    USE pression_mod
167    USE omega_mod
168    USE write_field_mod
169    USE vertical_interp_mod
170    USE wind_mod
171    TYPE(t_field),POINTER :: f_ps(:), f_phis(:), f_u(:), f_theta_rhodz(:), f_q(:), f_dps(:), &
172         f_buf_i(:), f_buf_v(:), f_buf_i3(:), f_buf1_i(:), f_buf2_i(:), f_buf_s(:), f_buf_p(:)
173   
174    REAL(rstd) :: out_pression_level
175    CHARACTER(LEN=255) :: str_pression
176    CHARACTER(LEN=255) :: physics_type
177   
178    out_pression_level=0.
179    CALL getin("out_pression_level",out_pression_level) 
180    WRITE(str_pression,*) INT(out_pression_level/100)
181    str_pression=ADJUSTL(str_pression)
182   
183    CALL writefield("ps",f_ps)
184    CALL writefield("dps",f_dps)
185    CALL writefield("phis",f_phis)
186    CALL vorticity(f_u,f_buf_v)
187    CALL writefield("vort",f_buf_v)
188
189    CALL w_omega(f_ps, f_u, f_buf_i)
190    CALL writefield('omega', f_buf_i)
191    IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
192      CALL vertical_interp(f_ps,f_buf_i,f_buf_s,out_pression_level)
193      CALL writefield("omega"//TRIM(str_pression),f_buf_s)
194    ENDIF
195   
196    ! Temperature
197!    CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_buf_i) ; ! FIXME
198     
199    CALL getin('physics',physics_type)
200    IF (TRIM(physics_type)=='dcmip') THEN
201      CALL Tv2T(f_buf_i,f_q,f_buf1_i) 
202      CALL writefield("T",f_buf1_i)
203      IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
204        CALL vertical_interp(f_ps,f_buf1_i,f_buf_s,out_pression_level)
205        CALL writefield("T"//TRIM(str_pression),f_buf_s)
206      ENDIF
207    ELSE
208      CALL writefield("T",f_buf_i)
209      IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
210        CALL vertical_interp(f_ps,f_buf_i,f_buf_s,out_pression_level)
211        CALL writefield("T"//TRIM(str_pression),f_buf_s)
212      ENDIF
213    ENDIF
214   
215    ! velocity components
216    CALL un2ulonlat(f_u, f_buf1_i, f_buf2_i)
217    CALL writefield("ulon",f_buf1_i)
218    CALL writefield("ulat",f_buf2_i)
219
220    IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
221      CALL vertical_interp(f_ps,f_buf1_i,f_buf_s,out_pression_level)
222      CALL writefield("ulon"//TRIM(str_pression),f_buf_s)
223      CALL vertical_interp(f_ps,f_buf2_i,f_buf_s,out_pression_level)
224      CALL writefield("ulat"//TRIM(str_pression),f_buf_s)
225    ENDIF
226   
227    ! geopotential ! FIXME
228    CALL thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_buf_s,f_buf_p,f_buf1_i,f_buf2_i,f_buf_i)
229    CALL writefield("p",f_buf_p)
230!    CALL writefield("phi",f_geopot)   ! geopotential
231    CALL writefield("theta",f_buf1_i) ! potential temperature
232    CALL writefield("pk",f_buf2_i)    ! Exner pressure
233 
234  END SUBROUTINE write_output_fields
235
236!------------------- Conversion from prognostic to observable variables ------------------
237
238  SUBROUTINE progonostic_vel_to_horiz(f_geopot, f_ps, f_rhodz, f_u, f_W, f_uh, f_uz)
239    USE disvert_mod
240    TYPE(t_field), POINTER :: f_geopot(:), f_ps(:), f_rhodz(:), &
241         f_u(:), f_W(:), f_uz(:), &  ! IN
242         f_uh(:)                         ! OUT
243    REAL(rstd),POINTER :: geopot(:,:), ps(:), rhodz(:,:), u(:,:), W(:,:), uh(:,:), uz(:,:)
244    INTEGER :: ind
245   
246    DO ind=1,ndomain
247       IF (.NOT. assigned_domain(ind)) CYCLE
248       CALL swap_dimensions(ind)
249       CALL swap_geometry(ind)
250       geopot = f_geopot(ind)
251       rhodz = f_rhodz(ind)
252       u = f_u(ind)
253       W = f_W(ind)
254       uh  = f_uh(ind)
255       IF(caldyn_eta==eta_mass) THEN
256          ps=f_ps(ind)
257          CALL compute_rhodz(.TRUE., ps, rhodz)
258       END IF
259       uz = f_uz(ind)
260       CALL compute_prognostic_vel_to_horiz(geopot,rhodz,u,W,uh,uz)
261    END DO
262  END SUBROUTINE progonostic_vel_to_horiz
263 
264  SUBROUTINE compute_prognostic_vel_to_horiz(Phi, rhodz, u, W, uh, uz)
265    USE omp_para
266    REAL(rstd), INTENT(IN) :: Phi(iim*jjm,llm+1)
267    REAL(rstd), INTENT(IN) :: rhodz(iim*jjm,llm)
268    REAL(rstd), INTENT(IN) :: u(3*iim*jjm,llm)
269    REAL(rstd), INTENT(IN) :: W(iim*jjm,llm+1)
270    REAL(rstd), INTENT(OUT) :: uh(3*iim*jjm,llm)
271    REAL(rstd), INTENT(OUT) :: uz(iim*jjm,llm)
272    INTEGER :: ij,l
273    REAL(rstd) :: F_el(3*iim*jjm,llm+1)
274    REAL(rstd) :: uu_right, uu_lup, uu_ldown, W_el, DePhil
275    ! NB : u and uh are not in DEC form, they are normal components   
276    ! => we must divide by de
277    IF(hydrostatic) THEN
278       uh(:,:)=u(:,:)
279       uz(:,:)=0.
280    ELSE
281       DO l=ll_begin, ll_endp1 ! compute on l levels (interfaces)
282          DO ij=ij_begin_ext, ij_end_ext
283             ! Compute on edge 'right'
284             W_el   = .5*( W(ij,l)+W(ij+t_right,l) )
285             DePhil = ne_right*(Phi(ij+t_right,l)-Phi(ij,l))
286             F_el(ij+u_right,l)   = DePhil*W_el/de(ij+u_right)
287             ! Compute on edge 'lup'
288             W_el   = .5*( W(ij,l)+W(ij+t_lup,l) )
289             DePhil = ne_lup*(Phi(ij+t_lup,l)-Phi(ij,l))
290             F_el(ij+u_lup,l)   = DePhil*W_el/de(ij+u_lup)
291             ! Compute on edge 'ldown'
292             W_el   = .5*( W(ij,l)+W(ij+t_ldown,l) )
293             DePhil = ne_ldown*(Phi(ij+t_ldown,l)-Phi(ij,l))
294             F_el(ij+u_ldown,l) = DePhil*W_el/de(ij+u_ldown)
295          END DO
296       END DO
297
298       DO l=ll_begin, ll_end ! compute on k levels (full levels)
299          DO ij=ij_begin_ext, ij_end_ext
300             ! w = vertical momentum = g^-2*dPhi/dt = uz/g
301             uz(ij,l) = (.5*g)*(W(ij,l)+W(ij,l+1))/rhodz(ij,l)
302             ! uh = u-w.grad(Phi) = u - uz.grad(z)
303             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))
304             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))
305             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))
306          END DO
307       END DO
308
309    END IF
310  END SUBROUTINE compute_prognostic_vel_to_horiz
311
312  SUBROUTINE thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_pks,f_p,f_theta,f_pk,f_phi) 
313    USE field_mod
314    USE pression_mod
315    USE exner_mod
316    USE geopotential_mod
317    USE theta2theta_rhodz_mod
318    TYPE(t_field), POINTER :: f_ps(:), f_phis(:), f_theta_rhodz(:), &  ! IN
319         f_pks(:), f_p(:), f_theta(:), f_pk(:), f_phi(:)               ! OUT
320    REAL(rstd),POINTER :: pk(:,:), p(:,:), theta(:,:), theta_rhodz(:,:,:), &
321         phi(:,:), phis(:), ps(:), pks(:)
322    INTEGER :: ind
323   
324    DO ind=1,ndomain
325       IF (.NOT. assigned_domain(ind)) CYCLE
326       CALL swap_dimensions(ind)
327       CALL swap_geometry(ind)
328       ps = f_ps(ind)
329       p  = f_p(ind)
330!$OMP BARRIER
331       CALL compute_pression(ps,p,0)
332       pk = f_pk(ind)
333       pks = f_pks(ind)
334!$OMP BARRIER
335       CALL compute_exner(ps,p,pks,pk,0)
336!$OMP BARRIER
337       theta_rhodz = f_theta_rhodz(ind)
338       theta = f_theta(ind)
339       CALL compute_theta_rhodz2theta(ps, theta_rhodz(:,:,1),theta,0)
340       phis = f_phis(ind)
341       phi = f_phi(ind)
342       CALL compute_geopotential(phis,pks,pk,theta,phi,0)
343    END DO
344
345  END SUBROUTINE thetarhodz2geopot
346 
347  SUBROUTINE Tv2T(f_Tv, f_q, f_T)
348    TYPE(t_field), POINTER :: f_TV(:)
349    TYPE(t_field), POINTER :: f_q(:)
350    TYPE(t_field), POINTER :: f_T(:)
351   
352    REAL(rstd),POINTER :: Tv(:,:), q(:,:,:), T(:,:)
353    INTEGER :: ind
354   
355    DO ind=1,ndomain
356       IF (.NOT. assigned_domain(ind)) CYCLE
357       CALL swap_dimensions(ind)
358       CALL swap_geometry(ind)
359       Tv=f_Tv(ind)
360       q=f_q(ind)
361       T=f_T(ind)
362       T=Tv/(1+0.608*q(:,:,1))
363    END DO
364   
365  END SUBROUTINE Tv2T
366
367  SUBROUTINE divide_by_mass(iq, f_mass, f_theta_rhodz, f_theta)
368    INTEGER, INTENT(IN) :: iq
369    TYPE(t_field), POINTER :: f_mass(:), f_theta_rhodz(:), f_theta(:)
370    REAL(rstd), POINTER :: mass(:,:), theta_rhodz(:,:,:), theta(:,:)
371    INTEGER :: ind
372    DO ind=1,ndomain
373       IF (.NOT. assigned_domain(ind)) CYCLE
374       CALL swap_dimensions(ind)
375       CALL swap_geometry(ind)
376       mass=f_mass(ind)
377       theta_rhodz=f_theta_rhodz(ind)
378       theta=f_theta(ind)
379       theta(:,:) = theta_rhodz(:,:,iq) / mass(:,:)
380    END DO
381  END SUBROUTINE divide_by_mass
382   
383END MODULE observable_mod
Note: See TracBrowser for help on using the repository browser.