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

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

Implement restartability for dynamico

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    IF(it==itaumax) THEN
183      lafin=.TRUE.
184    ELSE
185      lafin=.FALSE.
186    ENDIF
187    CALL transfert_message(f_u,req_u)
188    offset=0
189   
190    DO ind=1,ndomain
191      CALL swap_dimensions(ind)
192      CALL swap_geometry(ind)
193     
194      phis=f_phis(ind)
195      ps=f_ps(ind)
196      theta_rhodz=f_theta_rhodz(ind)
197      u=f_u(ind)
198      q=f_q(ind)
199      wflux=f_wflux(ind)
200      p=f_p(ind)
201      pks=f_pks(ind)
202      pk=f_pk(ind)
203      p_layer=f_p_layer(ind)
204      theta=f_theta(ind)
205      phi=f_phi(ind)
206      Temp=f_Temp(ind)
207      ulon=f_ulon(ind)
208      ulat=f_ulat(ind)
209           
210      CALL grid_icosa_to_physics
211     
212    ENDDO
213
214!$OMP BARRIER
215!$OMP MASTER
216    CALL SYSTEM_CLOCK(start_clock)
217!$OMP END MASTER
218    CALL calfis_icosa(current_time, dtphy, lafin, presnivs,      &
219                      p_phy, p_layer_phy, phi_phy, phis_phy, ulon_phy, ulat_phy, Temp_phy, q_phy, wflux_phy, &
220                      dulon_phy, dulat_phy, dTemp_phy, dq_phy, dps_phy  )
221
222
223!$OMP MASTER
224    CALL SYSTEM_CLOCK(stop_clock)
225    count_clock=count_clock+stop_clock-start_clock
226!$OMP END MASTER
227
228!$OMP BARRIER                       
229    offset=0 
230
231    DO ind=1,ndomain
232      CALL swap_dimensions(ind)
233      CALL swap_geometry(ind)
234     
235      theta_rhodz=f_theta_rhodz(ind)
236      u=f_u(ind)
237      q=f_q(ind)
238      ps=f_ps(ind)
239      dulon=f_dulon(ind)
240      dulat=f_dulat(ind)
241      Temp=f_temp(ind)
242      dTemp=f_dTemp(ind)
243      dq=f_dq(ind)
244      dps=f_dps(ind)
245      duc=f_duc(ind)
246      p=f_p(ind)
247      pks=f_pks(ind)
248      pk=f_pk(ind)
249     
250      CALL grid_physics_to_icosa
251    ENDDO
252
253!$OMP BARRIER
254     
255  CONTAINS
256   
257    SUBROUTINE grid_physics_to_icosa
258    USE theta2theta_rhodz_mod
259    USE omp_para
260    IMPLICIT NONE
261      INTEGER :: i,j,ij,l,iq
262         
263      DO j=jj_begin,jj_end
264        DO i=ii_begin,ii_end
265          ij=(j-1)*iim+i
266          offset=offset+1
267         
268          dulon(ij,ll_begin:ll_end)=dulon_phy(offset,ll_begin:ll_end)
269          dulat(ij,ll_begin:ll_end)=dulat_phy(offset,ll_begin:ll_end)
270          dTemp(ij,ll_begin:ll_end)=dTemp_phy(offset,ll_begin:ll_end)
271          Temp(ij,ll_begin:ll_end) = Temp_phy(offset,ll_begin:ll_end)
272          dq(ij,ll_begin:ll_end,:)=dq_phy(offset,ll_begin:ll_end,:)
273          if (omp_first) dps(ij)=dps_phy(offset)
274        ENDDO
275      ENDDO
276   
277      DO l=ll_begin,ll_end
278        DO j=jj_begin,jj_end
279          DO i=ii_begin,ii_end
280            ij=(j-1)*iim+i
281            duc(ij,:,l)=dulon(ij,l)*elon_i(ij,:)+dulat(ij,l)*elat_i(ij,:)
282          ENDDO
283        ENDDO
284      ENDDO
285
286      DO l=ll_begin,ll_end
287        DO j=jj_begin,jj_end
288          DO i=ii_begin,ii_end
289            ij=(j-1)*iim+i
290            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,:) )
291            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,:) )
292            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,:) )
293          ENDDO
294        ENDDO
295      ENDDO         
296
297      DO l=ll_begin,ll_end
298        DO j=jj_begin,jj_end
299          DO i=ii_begin,ii_end
300            ij=(j-1)*iim+i
301            Temp(ij,l)=Temp(ij,l)+ dtphy * dTemp(ij,l)
302          ENDDO
303        ENDDO
304      ENDDO         
305     
306      DO iq=1,nqtot
307        DO l=ll_begin,ll_end
308          DO j=jj_begin,jj_end
309            DO i=ii_begin,ii_end
310              ij=(j-1)*iim+i
311              q(ij,l,iq)=q(ij,l,iq)+ dtphy * dq(ij,l,iq)
312            ENDDO
313          ENDDO
314        ENDDO 
315      ENDDO
316
317!$OMP BARRIER
318     
319     IF (omp_first) THEN
320       DO j=jj_begin,jj_end
321         DO i=ii_begin,ii_end
322           ij=(j-1)*iim+i
323           ps(ij)=ps(ij)+ dtphy * dps(ij)
324          ENDDO
325       ENDDO
326     ENDIF
327
328!     CALL compute_temperature2theta_rhodz(ps,Temp,theta_rhodz,0)
329
330! compute pression
331!$OMP BARRIER
332      DO    l    = ll_begin,ll_endp1
333        DO j=jj_begin,jj_end
334          DO i=ii_begin,ii_end
335            ij=(j-1)*iim+i
336            p(ij,l) = ap(l) + bp(l) * ps(ij)
337          ENDDO
338        ENDDO
339      ENDDO
340
341!$OMP BARRIER
342
343! compute exner
344       
345       IF (omp_first) THEN
346         DO j=jj_begin,jj_end
347            DO i=ii_begin,ii_end
348               ij=(j-1)*iim+i
349               pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
350            ENDDO
351         ENDDO
352       ENDIF
353
354       ! 3D : pk
355       DO l = ll_begin,ll_end
356          DO j=jj_begin,jj_end
357             DO i=ii_begin,ii_end
358                ij=(j-1)*iim+i
359                pk(ij,l) = cpp * ((.5/preff)*(p(ij,l)+p(ij,l+1))) ** kappa
360             ENDDO
361          ENDDO
362       ENDDO
363
364!$OMP BARRIER
365
366!   compute theta, temperature and pression at layer
367    DO    l    = ll_begin, ll_end
368      DO j=jj_begin,jj_end
369        DO i=ii_begin,ii_end
370          ij=(j-1)*iim+i
371          theta_rhodz(ij,l) = temp(ij,l) * ((p(ij,l)-p(ij,l+1))/g) / (pk(ij,l) / cpp )
372        ENDDO
373      ENDDO
374    ENDDO
375   
376    END SUBROUTINE grid_physics_to_icosa
377   
378   
379    SUBROUTINE grid_icosa_to_physics
380    USE pression_mod
381    USE exner_mod
382    USE theta2theta_rhodz_mod
383    USE geopotential_mod
384    USE wind_mod
385    USE omp_para
386    IMPLICIT NONE
387   
388    REAL(rstd) :: uc(3)
389    INTEGER :: i,j,ij,l
390   
391
392! compute pression
393
394      DO    l    = ll_begin,ll_endp1
395        DO j=jj_begin,jj_end
396          DO i=ii_begin,ii_end
397            ij=(j-1)*iim+i
398            p(ij,l) = ap(l) + bp(l) * ps(ij)
399          ENDDO
400        ENDDO
401      ENDDO
402
403!$OMP BARRIER
404
405! compute exner
406       
407       IF (omp_first) THEN
408         DO j=jj_begin,jj_end
409            DO i=ii_begin,ii_end
410               ij=(j-1)*iim+i
411               pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
412            ENDDO
413         ENDDO
414       ENDIF
415
416       ! 3D : pk
417       DO l = ll_begin,ll_end
418          DO j=jj_begin,jj_end
419             DO i=ii_begin,ii_end
420                ij=(j-1)*iim+i
421                pk(ij,l) = cpp * ((.5/preff)*(p(ij,l)+p(ij,l+1))) ** kappa
422             ENDDO
423          ENDDO
424       ENDDO
425
426!$OMP BARRIER
427
428!   compute theta, temperature and pression at layer
429    DO    l    = ll_begin, ll_end
430      DO j=jj_begin,jj_end
431        DO i=ii_begin,ii_end
432          ij=(j-1)*iim+i
433          theta(ij,l) = theta_rhodz(ij,l) / ((p(ij,l)-p(ij,l+1))/g)
434          Temp(ij,l) = theta(ij,l) * pk(ij,l) / cpp
435          p_layer(ij,l)=preff*(pk(ij,l)/cpp)**(1./kappa) 
436        ENDDO
437      ENDDO
438    ENDDO
439
440
441!!! Compute geopotential
442       
443  ! for first layer
444  IF (omp_first) THEN
445    DO j=jj_begin,jj_end
446      DO i=ii_begin,ii_end
447        ij=(j-1)*iim+i
448        phi( ij,1 ) = phis( ij ) + theta(ij,1) * ( pks(ij) - pk(ij,1) )
449      ENDDO
450    ENDDO
451  ENDIF
452!!-> implicit flush on phi(:,1)
453         
454  ! for other layers
455   DO l = ll_beginp1, ll_end
456     DO j=jj_begin,jj_end
457       DO i=ii_begin,ii_end
458         ij=(j-1)*iim+i
459         phi(ij,l) =  0.5 * ( theta(ij,l)  + theta(ij,l-1) )  & 
460                          * (  pk(ij,l-1) -  pk(ij,l)    )
461       ENDDO
462     ENDDO
463   ENDDO       
464
465!$OMP BARRIER
466
467   DO l = 2, llm
468!$OMP DO
469     DO j=jj_begin,jj_end
470       DO i=ii_begin,ii_end
471         ij=(j-1)*iim+i
472         phi(ij,l) = phi(ij,l)+ phi(ij,l-1)
473       ENDDO
474     ENDDO
475   ENDDO
476! --> IMPLICIT FLUSH on phi
477
478
479
480! compute wind centered lon lat compound
481    DO l=ll_begin,ll_end
482      DO j=jj_begin,jj_end
483        DO i=ii_begin,ii_end
484          ij=(j-1)*iim+i
485          uc(:)=1/Ai(ij)*                                                                                                &
486                        ( 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,:))  &
487                         + 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,:))          &
488                         + 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,:))          &
489                         + 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,:))    &
490                         + 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,:))&
491                         + 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,:)))
492          ulon(ij,l)=sum(uc(:)*elon_i(ij,:))
493          ulat(ij,l)=sum(uc(:)*elat_i(ij,:)) 
494        ENDDO
495      ENDDO
496    ENDDO
497
498!$OMP BARRIER
499
500
501     
502      DO j=jj_begin,jj_end
503        DO i=ii_begin,ii_end
504          ij=(j-1)*iim+i
505          offset=offset+1
506
507          IF (omp_first) ps_phy(offset) = ps(ij)
508          p_phy(offset,ll_begin:ll_endp1) = p(ij,ll_begin:ll_endp1)
509          p_layer_phy(offset,ll_begin:ll_end) = p_layer(ij,ll_begin:ll_end)
510          Temp_phy(offset,ll_begin:ll_end) = Temp(ij,ll_begin:ll_end)
511          IF (omp_first) phis_phy(offset) = phis(ij)
512          phi_phy(offset,ll_begin:ll_end) = phi(ij,ll_begin:ll_end)-phis(ij)
513          ulon_phy(offset,ll_begin:ll_end) = ulon(ij,ll_begin:ll_end)
514          ulat_phy(offset,ll_begin:ll_end) = ulat(ij,ll_begin:ll_end)
515          q_phy(offset,ll_begin:ll_end,:) = q(ij,ll_begin:ll_end,:)
516          wflux_phy(offset,ll_begin:ll_end) = wflux(ij,ll_begin:ll_end)
517        ENDDO
518      ENDDO
519     
520     END SUBROUTINE grid_icosa_to_physics
521
522  END SUBROUTINE physics
523 
524 
525END MODULE physics_lmdz5_mod
526   
527 
Note: See TracBrowser for help on using the repository browser.