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

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

New : positive advection option for theta
File size: 9.8 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
13! temporary shared variable for caldyn
14  TYPE(t_field),POINTER, SAVE :: f_theta(:)
15
16  PUBLIC init_observable, write_output_fields_basic, f_theta
17
18CONTAINS
19 
20  SUBROUTINE init_observable
21    CALL allocate_field(f_buf_i,   field_t,type_real,llm,name="buffer_i")
22    CALL allocate_field(f_buf_p,   field_t,type_real,llm+1) 
23    CALL allocate_field(f_buf_u3d, field_t,type_real,3,llm)  ! 3D vel at cell centers
24    CALL allocate_field(f_buf_ulon,field_t,type_real,llm, name="buf_ulon")
25    CALL allocate_field(f_buf_ulat,field_t,type_real,llm, name="buf_ulat")
26    CALL allocate_field(f_buf_uh,  field_u,type_real,llm, name="buf_uh")
27    CALL allocate_field(f_buf_v,   field_z,type_real,llm, name="buf_v")
28    CALL allocate_field(f_buf_s,   field_t,type_real, name="buf_s")
29
30    CALL allocate_field(f_theta, field_t,type_real,llm,  name='theta')   ! potential temperature
31  END SUBROUTINE init_observable
32 
33  SUBROUTINE write_output_fields_basic(f_ps, f_mass, f_geopot, f_u, f_W, f_q)
34    USE wind_mod
35    USE output_field_mod
36    USE omp_para
37    TYPE(t_field),POINTER :: f_ps(:), f_mass(:), f_geopot(:), f_u(:), f_W(:), f_q(:)
38!    IF (is_master) PRINT *,'CALL write_output_fields_basic'
39    CALL progonostic_vel_to_horiz(f_geopot, f_ps, f_mass, f_u, f_W, f_buf_uh, f_buf_i)
40    CALL transfert_request(f_buf_uh,req_e1_vect) 
41    CALL output_field("uz",f_buf_i)
42    CALL un2ulonlat(f_buf_uh, f_buf_ulon, f_buf_ulat)
43    CALL output_field("ulon",f_buf_ulon)
44    CALL output_field("ulat",f_buf_ulat)
45    CALL output_field("ps",f_ps)
46    CALL output_field("Ai",geom%Ai)
47    !       CALL output_field("dps",f_dps)
48    CALL output_field("mass",f_mass)
49    CALL output_field("geopot",f_geopot)
50    !       CALL output_field("dmass",f_dmass)
51    !       CALL output_field("vort",f_qv)
52    CALL output_field("theta",f_theta)
53    !       CALL output_field("exner",f_pk)
54    !       CALL output_field("pv",f_qv)
55    CALL output_field("q",f_q)
56  END SUBROUTINE write_output_fields_basic
57 
58  SUBROUTINE write_output_fields(f_ps, f_phis, f_dps, f_u, f_theta_rhodz, f_q, &
59       f_buf_i, f_buf_v, f_buf_i3, f_buf1_i, f_buf2_i, f_buf_s, f_buf_p)
60    USE vorticity_mod
61    USE theta2theta_rhodz_mod
62    USE pression_mod
63    USE omega_mod
64    USE write_field_mod
65    USE vertical_interp_mod
66    USE wind_mod
67    TYPE(t_field),POINTER :: f_ps(:), f_phis(:), f_u(:), f_theta_rhodz(:), f_q(:), f_dps(:), &
68         f_buf_i(:), f_buf_v(:), f_buf_i3(:), f_buf1_i(:), f_buf2_i(:), f_buf_s(:), f_buf_p(:)
69   
70    REAL(rstd) :: out_pression_level
71    CHARACTER(LEN=255) :: str_pression
72    CHARACTER(LEN=255) :: physics_type
73   
74    out_pression_level=0.
75    CALL getin("out_pression_level",out_pression_level) 
76    WRITE(str_pression,*) INT(out_pression_level/100)
77    str_pression=ADJUSTL(str_pression)
78   
79    CALL writefield("ps",f_ps)
80    CALL writefield("dps",f_dps)
81    CALL writefield("phis",f_phis)
82    CALL vorticity(f_u,f_buf_v)
83    CALL writefield("vort",f_buf_v)
84
85    CALL w_omega(f_ps, f_u, f_buf_i)
86    CALL writefield('omega', f_buf_i)
87    IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
88      CALL vertical_interp(f_ps,f_buf_i,f_buf_s,out_pression_level)
89      CALL writefield("omega"//TRIM(str_pression),f_buf_s)
90    ENDIF
91   
92    ! Temperature
93!    CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_buf_i) ; ! FIXME
94     
95    CALL getin('physics',physics_type)
96    IF (TRIM(physics_type)=='dcmip') THEN
97      CALL Tv2T(f_buf_i,f_q,f_buf1_i) 
98      CALL writefield("T",f_buf1_i)
99      IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
100        CALL vertical_interp(f_ps,f_buf1_i,f_buf_s,out_pression_level)
101        CALL writefield("T"//TRIM(str_pression),f_buf_s)
102      ENDIF
103    ELSE
104      CALL writefield("T",f_buf_i)
105      IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
106        CALL vertical_interp(f_ps,f_buf_i,f_buf_s,out_pression_level)
107        CALL writefield("T"//TRIM(str_pression),f_buf_s)
108      ENDIF
109    ENDIF
110   
111    ! velocity components
112    CALL un2ulonlat(f_u, f_buf1_i, f_buf2_i)
113    CALL writefield("ulon",f_buf1_i)
114    CALL writefield("ulat",f_buf2_i)
115
116    IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
117      CALL vertical_interp(f_ps,f_buf1_i,f_buf_s,out_pression_level)
118      CALL writefield("ulon"//TRIM(str_pression),f_buf_s)
119      CALL vertical_interp(f_ps,f_buf2_i,f_buf_s,out_pression_level)
120      CALL writefield("ulat"//TRIM(str_pression),f_buf_s)
121    ENDIF
122   
123    ! geopotential ! FIXME
124    CALL thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_buf_s,f_buf_p,f_buf1_i,f_buf2_i,f_buf_i)
125    CALL writefield("p",f_buf_p)
126!    CALL writefield("phi",f_geopot)   ! geopotential
127    CALL writefield("theta",f_buf1_i) ! potential temperature
128    CALL writefield("pk",f_buf2_i)    ! Exner pressure
129 
130  END SUBROUTINE write_output_fields
131
132!------------------- Conversion from prognostic to observable variables ------------------
133
134  SUBROUTINE progonostic_vel_to_horiz(f_geopot, f_ps, f_rhodz, f_u, f_W, f_uh, f_uz)
135    USE disvert_mod
136    TYPE(t_field), POINTER :: f_geopot(:), f_ps(:), f_rhodz(:), &
137         f_u(:), f_W(:), f_uz(:), &  ! IN
138         f_uh(:)                         ! OUT
139    REAL(rstd),POINTER :: geopot(:,:), ps(:), rhodz(:,:), u(:,:), W(:,:), uh(:,:), uz(:,:)
140    INTEGER :: ind
141   
142    DO ind=1,ndomain
143       IF (.NOT. assigned_domain(ind)) CYCLE
144       CALL swap_dimensions(ind)
145       CALL swap_geometry(ind)
146       geopot = f_geopot(ind)
147       rhodz = f_rhodz(ind)
148       u = f_u(ind)
149       W = f_W(ind)
150       uh  = f_uh(ind)
151       IF(caldyn_eta==eta_mass) THEN
152          ps=f_ps(ind)
153          CALL compute_rhodz(.TRUE., ps, rhodz)
154       END IF
155       uz = f_uz(ind)
156       CALL compute_prognostic_vel_to_horiz(geopot,rhodz,u,W,uh,uz)
157    END DO
158  END SUBROUTINE progonostic_vel_to_horiz
159 
160  SUBROUTINE compute_prognostic_vel_to_horiz(Phi, rhodz, u, W, uh, uz)
161    USE omp_para
162    REAL(rstd), INTENT(IN) :: Phi(iim*jjm,llm+1)
163    REAL(rstd), INTENT(IN) :: rhodz(iim*jjm,llm)
164    REAL(rstd), INTENT(IN) :: u(3*iim*jjm,llm)
165    REAL(rstd), INTENT(IN) :: W(iim*jjm,llm+1)
166    REAL(rstd), INTENT(OUT) :: uh(3*iim*jjm,llm)
167    REAL(rstd), INTENT(OUT) :: uz(iim*jjm,llm)
168    INTEGER :: ij,l
169    REAL(rstd) :: F_el(3*iim*jjm,llm+1)
170    REAL(rstd) :: uu_right, uu_lup, uu_ldown, W_el, DePhil
171    ! NB : u and uh are not in DEC form, they are normal components   
172    ! => we must divide by de
173    IF(hydrostatic) THEN
174       uh(:,:)=u(:,:)
175       uz(:,:)=0.
176    ELSE
177       DO l=ll_begin, ll_endp1 ! compute on l levels (interfaces)
178          DO ij=ij_begin_ext, ij_end_ext
179             ! Compute on edge 'right'
180             W_el   = .5*( W(ij,l)+W(ij+t_right,l) )
181             DePhil = ne_right*(Phi(ij+t_right,l)-Phi(ij,l))
182             F_el(ij+u_right,l)   = DePhil*W_el/de(ij+u_right)
183             ! Compute on edge 'lup'
184             W_el   = .5*( W(ij,l)+W(ij+t_lup,l) )
185             DePhil = ne_lup*(Phi(ij+t_lup,l)-Phi(ij,l))
186             F_el(ij+u_lup,l)   = DePhil*W_el/de(ij+u_lup)
187             ! Compute on edge 'ldown'
188             W_el   = .5*( W(ij,l)+W(ij+t_ldown,l) )
189             DePhil = ne_ldown*(Phi(ij+t_ldown,l)-Phi(ij,l))
190             F_el(ij+u_ldown,l) = DePhil*W_el/de(ij+u_ldown)
191          END DO
192       END DO
193
194       DO l=ll_begin, ll_end ! compute on k levels (full levels)
195          DO ij=ij_begin_ext, ij_end_ext
196             ! w = vertical momentum = g^-2*dPhi/dt = uz/g
197             uz(ij,l) = (.5*g)*(W(ij,l)+W(ij,l+1))/rhodz(ij,l)
198             ! uh = u-w.grad(Phi) = u - uz.grad(z)
199             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))
200             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))
201             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))
202          END DO
203       END DO
204
205    END IF
206  END SUBROUTINE compute_prognostic_vel_to_horiz
207
208  SUBROUTINE thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_pks,f_p,f_theta,f_pk,f_phi) 
209    USE field_mod
210    USE pression_mod
211    USE exner_mod
212    USE geopotential_mod
213    USE theta2theta_rhodz_mod
214    TYPE(t_field), POINTER :: f_ps(:), f_phis(:), f_theta_rhodz(:), &  ! IN
215         f_pks(:), f_p(:), f_theta(:), f_pk(:), f_phi(:)               ! OUT
216    REAL(rstd),POINTER :: pk(:,:), p(:,:), theta(:,:), theta_rhodz(:,:), &
217         phi(:,:), phis(:), ps(:), pks(:)
218    INTEGER :: ind
219   
220    DO ind=1,ndomain
221       IF (.NOT. assigned_domain(ind)) CYCLE
222       CALL swap_dimensions(ind)
223       CALL swap_geometry(ind)
224       ps = f_ps(ind)
225       p  = f_p(ind)
226!$OMP BARRIER
227       CALL compute_pression(ps,p,0)
228       pk = f_pk(ind)
229       pks = f_pks(ind)
230!$OMP BARRIER
231       CALL compute_exner(ps,p,pks,pk,0)
232!$OMP BARRIER
233       theta_rhodz = f_theta_rhodz(ind)
234       theta = f_theta(ind)
235       CALL compute_theta_rhodz2theta(ps, theta_rhodz,theta,0)
236       phis = f_phis(ind)
237       phi = f_phi(ind)
238       CALL compute_geopotential(phis,pks,pk,theta,phi,0)
239    END DO
240
241  END SUBROUTINE thetarhodz2geopot
242 
243  SUBROUTINE Tv2T(f_Tv, f_q, f_T)
244    TYPE(t_field), POINTER :: f_TV(:)
245    TYPE(t_field), POINTER :: f_q(:)
246    TYPE(t_field), POINTER :: f_T(:)
247   
248    REAL(rstd),POINTER :: Tv(:,:), q(:,:,:), T(:,:)
249    INTEGER :: ind
250   
251    DO ind=1,ndomain
252       IF (.NOT. assigned_domain(ind)) CYCLE
253       CALL swap_dimensions(ind)
254       CALL swap_geometry(ind)
255       Tv=f_Tv(ind)
256       q=f_q(ind)
257       T=f_T(ind)
258       T=Tv/(1+0.608*q(:,:,1))
259    END DO
260   
261  END SUBROUTINE Tv2T
262 
263END MODULE observable_mod
Note: See TracBrowser for help on using the repository browser.