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

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

Creating temporary dynamico/lmdz/saturn branche

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