source: codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/physics_lmdz5.f90 @ 245

Last change on this file since 245 was 221, checked in by ymipsl, 10 years ago

Creating temporary dynamico/lmdz/saturn branche

YM
File size: 14.5 KB
Line 
1MODULE physics_lmdz5_mod
2  USE field_mod
3  USE transfert_mod
4 
5  INTEGER,SAVE :: nbp_phys
6  TYPE(t_message) :: req_u
7
8  TYPE(t_field),POINTER :: f_p(:) 
9  TYPE(t_field),POINTER :: f_pks(:) 
10  TYPE(t_field),POINTER :: f_pk(:) 
11  TYPE(t_field),POINTER :: f_p_layer(:)   
12  TYPE(t_field),POINTER :: f_theta(:)   
13  TYPE(t_field),POINTER :: f_phi(:)   
14  TYPE(t_field),POINTER :: f_Temp(:)   
15  TYPE(t_field),POINTER :: f_ulon(:)   
16  TYPE(t_field),POINTER :: f_ulat(:)   
17  TYPE(t_field),POINTER :: f_dulon(:)
18  TYPE(t_field),POINTER :: f_dulat(:)
19  TYPE(t_field),POINTER :: f_dTemp(:)
20  TYPE(t_field),POINTER :: f_dq(:)
21  TYPE(t_field),POINTER :: f_dps(:)
22  TYPE(t_field),POINTER :: f_duc(:)
23
24  INTEGER :: start_clock
25  INTEGER :: stop_clock
26  INTEGER :: count_clock=0
27
28CONTAINS
29
30  SUBROUTINE init_physics
31  USE icosa
32  USE domain_mod
33  USE dimensions
34  USE mpi_mod
35  USE mpipara
36  IMPLICIT NONE
37  INTEGER  :: distrib(0:mpi_size-1)
38  INTEGER  :: ind,i,j,ij,pos
39
40  REAL(rstd),ALLOCATABLE :: latfi(:)
41  REAL(rstd),ALLOCATABLE :: lonfi(:)
42  REAL(rstd),ALLOCATABLE :: airefi(:)
43   
44
45    CALL allocate_field(f_p,field_t,type_real,llm+1)
46    CALL allocate_field(f_pks,field_t,type_real)
47    CALL allocate_field(f_pk,field_t,type_real,llm)
48    CALL allocate_field(f_p_layer,field_t,type_real,llm)
49    CALL allocate_field(f_theta,field_t,type_real,llm)
50    CALL allocate_field(f_phi,field_t,type_real,llm)
51    CALL allocate_field(f_Temp,field_t,type_real,llm)
52    CALL allocate_field(f_ulon,field_t,type_real,llm)
53    CALL allocate_field(f_ulat,field_t,type_real,llm)
54    CALL allocate_field(f_dulon,field_t,type_real,llm)
55    CALL allocate_field(f_dulat,field_t,type_real,llm)
56    CALL allocate_field(f_dTemp,field_t,type_real,llm)
57    CALL allocate_field(f_dq,field_t,type_real,llm,nqtot)
58    CALL allocate_field(f_dps,field_t,type_real)
59    CALL allocate_field(f_duc,field_t,type_real,3,llm)   
60   
61   
62    nbp_phys=0
63    DO ind=1,ndomain
64      CALL swap_dimensions(ind)
65      nbp_phys=nbp_phys+ii_nb*jj_nb
66    ENDDO
67   
68    CALL MPI_ALLGATHER(nbp_phys,1,MPI_INTEGER,distrib,1,MPI_INTEGER,comm_icosa,ierr)
69   
70    ALLOCATE(latfi(nbp_phys))
71    ALLOCATE(lonfi(nbp_phys))
72    ALLOCATE(airefi(nbp_phys))
73   
74    pos=0
75    DO ind=1,ndomain
76      CALL swap_dimensions(ind)
77      CALL swap_geometry(ind)
78      DO j=jj_begin,jj_end
79        DO i=ii_begin,ii_end
80          ij=(j-1)*iim+i
81          pos=pos+1
82          CALL xyz2lonlat(xyz_i(ij,:),lonfi(pos),latfi(pos))
83          airefi(pos)=Ai(ij) 
84        ENDDO
85      ENDDO
86    ENDDO
87   
88    CALL init_gcm_lmdz(nbp_phys,mpi_size,distrib,latfi,lonfi,airefi)
89   
90  END SUBROUTINE init_physics
91 
92  SUBROUTINE physics(it,f_phis, f_ps, f_theta_rhodz, f_u, f_wflux, f_q)
93  USE ICOSA
94  USE time_mod
95  USE disvert_mod
96  USE transfert_mod
97  IMPLICIT NONE
98    INTEGER,INTENT(IN)    :: it
99    TYPE(t_field),POINTER :: f_phis(:)
100    TYPE(t_field),POINTER :: f_ps(:)
101    TYPE(t_field),POINTER :: f_theta_rhodz(:)
102    TYPE(t_field),POINTER :: f_u(:)
103    TYPE(t_field),POINTER :: f_wflux(:)
104    TYPE(t_field),POINTER :: f_q(:)
105 
106    REAL(rstd),POINTER :: phis(:)
107    REAL(rstd),POINTER :: ps(:)
108    REAL(rstd),POINTER :: theta_rhodz(:,:)
109    REAL(rstd),POINTER :: u(:,:)
110    REAL(rstd),POINTER :: wflux(:,:)
111    REAL(rstd),POINTER :: q(:,:,:)
112    REAL(rstd),POINTER :: p(:,:)
113    REAL(rstd),POINTER :: pks(:)
114    REAL(rstd),POINTER :: pk(:,:)
115    REAL(rstd),POINTER :: p_layer(:,:)
116    REAL(rstd),POINTER :: theta(:,:)
117    REAL(rstd),POINTER :: phi(:,:)
118    REAL(rstd),POINTER :: Temp(:,:)
119    REAL(rstd),POINTER :: ulon(:,:)
120    REAL(rstd),POINTER :: ulat(:,:)
121    REAL(rstd),POINTER :: dulon(:,:)
122    REAL(rstd),POINTER :: dulat(:,:)
123    REAL(rstd),POINTER :: dTemp(:,:)
124    REAL(rstd),POINTER :: dq(:,:,:)
125    REAL(rstd),POINTER :: dps(:)
126    REAL(rstd),POINTER :: duc(:,:,:)
127
128
129    INTEGER :: ind
130   
131    REAL(rstd),ALLOCATABLE,SAVE :: ps_phy(:)
132    REAL(rstd),ALLOCATABLE,SAVE :: p_phy(:,:)
133    REAL(rstd),ALLOCATABLE,SAVE :: p_layer_phy(:,:)
134    REAL(rstd),ALLOCATABLE,SAVE :: Temp_phy(:,:)
135    REAL(rstd),ALLOCATABLE,SAVE :: phis_phy(:)
136    REAL(rstd),ALLOCATABLE,SAVE :: phi_phy(:,:)
137    REAL(rstd),ALLOCATABLE,SAVE :: ulon_phy(:,:)
138    REAL(rstd),ALLOCATABLE,SAVE :: ulat_phy(:,:)
139    REAL(rstd),ALLOCATABLE,SAVE :: q_phy(:,:,:)
140    REAL(rstd),ALLOCATABLE,SAVE :: wflux_phy(:,:)
141    REAL(rstd),ALLOCATABLE,SAVE :: dulon_phy(:,:)
142    REAL(rstd),ALLOCATABLE,SAVE :: dulat_phy(:,:)
143    REAL(rstd),ALLOCATABLE,SAVE :: dTemp_phy(:,:)
144    REAL(rstd),ALLOCATABLE,SAVE :: dq_phy(:,:,:)
145    REAL(rstd),ALLOCATABLE,SAVE :: dps_phy(:)
146    REAL :: dtphy 
147    INTEGER,SAVE :: offset
148!$OMP THREADPRIVATE(offset)
149    LOGICAL :: lafin=.FALSE.
150    LOGICAL,SAVE :: first=.TRUE.
151!$OMP THREADPRIVATE(first)
152
153   
154    IF (first) THEN
155      first=.FALSE.
156      CALL init_message(f_u,req_e1_vect,req_u)
157!$OMP MASTER
158      ALLOCATE(ps_phy(nbp_phys))
159      ALLOCATE(p_phy(nbp_phys,llm+1))
160      ALLOCATE(p_layer_phy(nbp_phys,llm))
161      ALLOCATE(Temp_phy(nbp_phys,llm))
162      ALLOCATE(phis_phy(nbp_phys))
163      ALLOCATE(phi_phy(nbp_phys,llm))
164      ALLOCATE(ulon_phy(nbp_phys,llm))
165      ALLOCATE(ulat_phy(nbp_phys,llm))
166      ALLOCATE(q_phy(nbp_phys,llm,nqtot))
167      ALLOCATE(wflux_phy(nbp_phys,llm))
168      ALLOCATE(dulon_phy(nbp_phys,llm))
169      ALLOCATE(dulat_phy(nbp_phys,llm))
170      ALLOCATE(dTemp_phy(nbp_phys,llm))
171      ALLOCATE(dq_phy(nbp_phys,llm,nqtot))
172      ALLOCATE(dps_phy(nbp_phys))
173!$OMP END MASTER
174!$OMP BARRIER
175    ENDIF
176
177!$OMP MASTER       
178    CALL update_calendar(it)
179!$OMP END MASTER
180!$OMP BARRIER
181    dtphy=itau_physics*dt
182   
183    CALL transfert_message(f_u,req_u)
184    offset=0
185   
186    DO ind=1,ndomain
187      CALL swap_dimensions(ind)
188      CALL swap_geometry(ind)
189     
190      phis=f_phis(ind)
191      ps=f_ps(ind)
192      theta_rhodz=f_theta_rhodz(ind)
193      u=f_u(ind)
194      q=f_q(ind)
195      wflux=f_wflux(ind)
196      p=f_p(ind)
197      pks=f_pks(ind)
198      pk=f_pk(ind)
199      p_layer=f_p_layer(ind)
200      theta=f_theta(ind)
201      phi=f_phi(ind)
202      Temp=f_Temp(ind)
203      ulon=f_ulon(ind)
204      ulat=f_ulat(ind)
205           
206      CALL grid_icosa_to_physics
207     
208    ENDDO
209
210!$OMP BARRIER
211!$OMP MASTER
212    CALL SYSTEM_CLOCK(start_clock)
213!$OMP END MASTER
214    CALL calfis_icosa(current_time, dtphy, lafin, presnivs,      &
215                      p_phy, p_layer_phy, phi_phy, phis_phy, ulon_phy, ulat_phy, Temp_phy, q_phy, wflux_phy, &
216                      dulon_phy, dulat_phy, dTemp_phy, dq_phy, dps_phy  )
217
218
219!$OMP MASTER
220    CALL SYSTEM_CLOCK(stop_clock)
221    count_clock=count_clock+stop_clock-start_clock
222!$OMP END MASTER
223
224!$OMP BARRIER                       
225    offset=0 
226
227    DO ind=1,ndomain
228      CALL swap_dimensions(ind)
229      CALL swap_geometry(ind)
230     
231      theta_rhodz=f_theta_rhodz(ind)
232      u=f_u(ind)
233      q=f_q(ind)
234      ps=f_ps(ind)
235      dulon=f_dulon(ind)
236      dulat=f_dulat(ind)
237      Temp=f_temp(ind)
238      dTemp=f_dTemp(ind)
239      dq=f_dq(ind)
240      dps=f_dps(ind)
241      duc=f_duc(ind)
242      p=f_p(ind)
243      pks=f_pks(ind)
244      pk=f_pk(ind)
245     
246      CALL grid_physics_to_icosa
247    ENDDO
248
249!$OMP BARRIER
250     
251  CONTAINS
252   
253    SUBROUTINE grid_physics_to_icosa
254    USE theta2theta_rhodz_mod
255    USE omp_para
256    IMPLICIT NONE
257      INTEGER :: i,j,ij,l,iq
258         
259      DO j=jj_begin,jj_end
260        DO i=ii_begin,ii_end
261          ij=(j-1)*iim+i
262          offset=offset+1
263         
264          dulon(ij,ll_begin:ll_end)=dulon_phy(offset,ll_begin:ll_end)
265          dulat(ij,ll_begin:ll_end)=dulat_phy(offset,ll_begin:ll_end)
266          dTemp(ij,ll_begin:ll_end)=dTemp_phy(offset,ll_begin:ll_end)
267          Temp(ij,ll_begin:ll_end) = Temp_phy(offset,ll_begin:ll_end)
268          dq(ij,ll_begin:ll_end,:)=dq_phy(offset,ll_begin:ll_end,:)
269          if (omp_first) dps(ij)=dps_phy(offset)
270        ENDDO
271      ENDDO
272   
273      DO l=ll_begin,ll_end
274        DO j=jj_begin,jj_end
275          DO i=ii_begin,ii_end
276            ij=(j-1)*iim+i
277            duc(ij,:,l)=dulon(ij,l)*elon_i(ij,:)+dulat(ij,l)*elat_i(ij,:)
278          ENDDO
279        ENDDO
280      ENDDO
281
282      DO l=ll_begin,ll_end
283        DO j=jj_begin,jj_end
284          DO i=ii_begin,ii_end
285            ij=(j-1)*iim+i
286            u(ij+u_right,l) = u(ij+u_right,l) + dtphy * sum( 0.5*(duc(ij,:,l) + duc(ij+t_right,:,l))*ep_e(ij+u_right,:) )
287            u(ij+u_lup,l) = u(ij+u_lup,l) + dtphy * sum( 0.5*(duc(ij,:,l) + duc(ij+t_lup,:,l))*ep_e(ij+u_lup,:) )
288            u(ij+u_ldown,l) = u(ij+u_ldown,l) + dtphy*sum( 0.5*(duc(ij,:,l) + duc(ij+t_ldown,:,l))*ep_e(ij+u_ldown,:) )
289          ENDDO
290        ENDDO
291      ENDDO         
292
293      DO l=ll_begin,ll_end
294        DO j=jj_begin,jj_end
295          DO i=ii_begin,ii_end
296            ij=(j-1)*iim+i
297            Temp(ij,l)=Temp(ij,l)+ dtphy * dTemp(ij,l)
298          ENDDO
299        ENDDO
300      ENDDO         
301     
302      DO iq=1,nqtot
303        DO l=ll_begin,ll_end
304          DO j=jj_begin,jj_end
305            DO i=ii_begin,ii_end
306              ij=(j-1)*iim+i
307              q(ij,l,iq)=q(ij,l,iq)+ dtphy * dq(ij,l,iq)
308            ENDDO
309          ENDDO
310        ENDDO 
311      ENDDO
312
313!$OMP BARRIER
314     
315     IF (omp_first) THEN
316       DO j=jj_begin,jj_end
317         DO i=ii_begin,ii_end
318           ij=(j-1)*iim+i
319           ps(ij)=ps(ij)+ dtphy * dps(ij)
320          ENDDO
321       ENDDO
322     ENDIF
323
324!     CALL compute_temperature2theta_rhodz(ps,Temp,theta_rhodz,0)
325
326! compute pression
327!$OMP BARRIER
328      DO    l    = ll_begin,ll_endp1
329        DO j=jj_begin,jj_end
330          DO i=ii_begin,ii_end
331            ij=(j-1)*iim+i
332            p(ij,l) = ap(l) + bp(l) * ps(ij)
333          ENDDO
334        ENDDO
335      ENDDO
336
337!$OMP BARRIER
338
339! compute exner
340       
341       IF (omp_first) THEN
342         DO j=jj_begin,jj_end
343            DO i=ii_begin,ii_end
344               ij=(j-1)*iim+i
345               pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
346            ENDDO
347         ENDDO
348       ENDIF
349
350       ! 3D : pk
351       DO l = ll_begin,ll_end
352          DO j=jj_begin,jj_end
353             DO i=ii_begin,ii_end
354                ij=(j-1)*iim+i
355                pk(ij,l) = cpp * ((.5/preff)*(p(ij,l)+p(ij,l+1))) ** kappa
356             ENDDO
357          ENDDO
358       ENDDO
359
360!$OMP BARRIER
361
362!   compute theta, temperature and pression at layer
363    DO    l    = ll_begin, ll_end
364      DO j=jj_begin,jj_end
365        DO i=ii_begin,ii_end
366          ij=(j-1)*iim+i
367          theta_rhodz(ij,l) = temp(ij,l) * ((p(ij,l)-p(ij,l+1))/g) / (pk(ij,l) / cpp )
368        ENDDO
369      ENDDO
370    ENDDO
371   
372    END SUBROUTINE grid_physics_to_icosa
373   
374   
375    SUBROUTINE grid_icosa_to_physics
376    USE pression_mod
377    USE exner_mod
378    USE theta2theta_rhodz_mod
379    USE geopotential_mod
380    USE wind_mod
381    USE omp_para
382    IMPLICIT NONE
383   
384    REAL(rstd) :: uc(3)
385    INTEGER :: i,j,ij,l
386   
387
388! compute pression
389
390      DO    l    = ll_begin,ll_endp1
391        DO j=jj_begin,jj_end
392          DO i=ii_begin,ii_end
393            ij=(j-1)*iim+i
394            p(ij,l) = ap(l) + bp(l) * ps(ij)
395          ENDDO
396        ENDDO
397      ENDDO
398
399!$OMP BARRIER
400
401! compute exner
402       
403       IF (omp_first) THEN
404         DO j=jj_begin,jj_end
405            DO i=ii_begin,ii_end
406               ij=(j-1)*iim+i
407               pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
408            ENDDO
409         ENDDO
410       ENDIF
411
412       ! 3D : pk
413       DO l = ll_begin,ll_end
414          DO j=jj_begin,jj_end
415             DO i=ii_begin,ii_end
416                ij=(j-1)*iim+i
417                pk(ij,l) = cpp * ((.5/preff)*(p(ij,l)+p(ij,l+1))) ** kappa
418             ENDDO
419          ENDDO
420       ENDDO
421
422!$OMP BARRIER
423
424!   compute theta, temperature and pression at layer
425    DO    l    = ll_begin, ll_end
426      DO j=jj_begin,jj_end
427        DO i=ii_begin,ii_end
428          ij=(j-1)*iim+i
429          theta(ij,l) = theta_rhodz(ij,l) / ((p(ij,l)-p(ij,l+1))/g)
430          Temp(ij,l) = theta(ij,l) * pk(ij,l) / cpp
431          p_layer(ij,l)=preff*(pk(ij,l)/cpp)**(1./kappa) 
432        ENDDO
433      ENDDO
434    ENDDO
435
436
437!!! Compute geopotential
438       
439  ! for first layer
440  IF (omp_first) THEN
441    DO j=jj_begin,jj_end
442      DO i=ii_begin,ii_end
443        ij=(j-1)*iim+i
444        phi( ij,1 ) = phis( ij ) + theta(ij,1) * ( pks(ij) - pk(ij,1) )
445      ENDDO
446    ENDDO
447  ENDIF
448!!-> implicit flush on phi(:,1)
449         
450  ! for other layers
451   DO l = ll_beginp1, ll_end
452     DO j=jj_begin,jj_end
453       DO i=ii_begin,ii_end
454         ij=(j-1)*iim+i
455         phi(ij,l) =  0.5 * ( theta(ij,l)  + theta(ij,l-1) )  & 
456                          * (  pk(ij,l-1) -  pk(ij,l)    )
457       ENDDO
458     ENDDO
459   ENDDO       
460
461!$OMP BARRIER
462
463   DO l = 2, llm
464!$OMP DO
465     DO j=jj_begin,jj_end
466       DO i=ii_begin,ii_end
467         ij=(j-1)*iim+i
468         phi(ij,l) = phi(ij,l)+ phi(ij,l-1)
469       ENDDO
470     ENDDO
471   ENDDO
472! --> IMPLICIT FLUSH on phi
473
474
475
476! compute wind centered lon lat compound
477    DO l=ll_begin,ll_end
478      DO j=jj_begin,jj_end
479        DO i=ii_begin,ii_end
480          ij=(j-1)*iim+i
481          uc(:)=1/Ai(ij)*                                                                                                &
482                        ( ne(ij,right)*u(ij+u_right,l)*le(ij+u_right)*((xyz_v(ij+z_rdown,:)+xyz_v(ij+z_rup,:))/2-centroid(ij,:))  &
483                         + ne(ij,rup)*u(ij+u_rup,l)*le(ij+u_rup)*((xyz_v(ij+z_rup,:)+xyz_v(ij+z_up,:))/2-centroid(ij,:))          &
484                         + ne(ij,lup)*u(ij+u_lup,l)*le(ij+u_lup)*((xyz_v(ij+z_up,:)+xyz_v(ij+z_lup,:))/2-centroid(ij,:))          &
485                         + ne(ij,left)*u(ij+u_left,l)*le(ij+u_left)*((xyz_v(ij+z_lup,:)+xyz_v(ij+z_ldown,:))/2-centroid(ij,:))    &
486                         + ne(ij,ldown)*u(ij+u_ldown,l)*le(ij+u_ldown)*((xyz_v(ij+z_ldown,:)+xyz_v(ij+z_down,:))/2-centroid(ij,:))&
487                         + ne(ij,rdown)*u(ij+u_rdown,l)*le(ij+u_rdown)*((xyz_v(ij+z_down,:)+xyz_v(ij+z_rdown,:))/2-centroid(ij,:)))
488          ulon(ij,l)=sum(uc(:)*elon_i(ij,:))
489          ulat(ij,l)=sum(uc(:)*elat_i(ij,:)) 
490        ENDDO
491      ENDDO
492    ENDDO
493
494!$OMP BARRIER
495
496
497     
498      DO j=jj_begin,jj_end
499        DO i=ii_begin,ii_end
500          ij=(j-1)*iim+i
501          offset=offset+1
502
503          IF (omp_first) ps_phy(offset) = ps(ij)
504          p_phy(offset,ll_begin:ll_endp1) = p(ij,ll_begin:ll_endp1)
505          p_layer_phy(offset,ll_begin:ll_end) = p_layer(ij,ll_begin:ll_end)
506          Temp_phy(offset,ll_begin:ll_end) = Temp(ij,ll_begin:ll_end)
507          IF (omp_first) phis_phy(offset) = phis(ij)
508          phi_phy(offset,ll_begin:ll_end) = phi(ij,ll_begin:ll_end)-phis(ij)
509          ulon_phy(offset,ll_begin:ll_end) = ulon(ij,ll_begin:ll_end)
510          ulat_phy(offset,ll_begin:ll_end) = ulat(ij,ll_begin:ll_end)
511          q_phy(offset,ll_begin:ll_end,:) = q(ij,ll_begin:ll_end,:)
512          wflux_phy(offset,ll_begin:ll_end) = wflux(ij,ll_begin:ll_end)
513        ENDDO
514      ENDDO
515     
516     END SUBROUTINE grid_icosa_to_physics
517
518  END SUBROUTINE physics
519 
520 
521END MODULE physics_lmdz5_mod
522   
523 
Note: See TracBrowser for help on using the repository browser.