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

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

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