source: codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/physics_lmdz_generic.f90 @ 239

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

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