source: trunk/LMDZ.VENUS/libf/phyvenus/iophy.F90 @ 3556

Last change on this file since 3556 was 1719, checked in by emillour, 7 years ago

Venus GCM:
Stop using "abort_gcm" in the physics (it is a routine from the dynamics); use "abort_physic" instead.
EM

File size: 13.1 KB
Line 
1!
2! $Header$
3!
4module iophy
5 
6! abd  REAL,private,allocatable,dimension(:),save :: io_lat
7! abd  REAL,private,allocatable,dimension(:),save :: io_lon
8  REAL,allocatable,dimension(:),save :: io_lat
9  REAL,allocatable,dimension(:),save :: io_lon
10  INTEGER, save :: phys_domain_id
11  INTEGER, save :: npstn
12  INTEGER, allocatable, dimension(:), save :: nptabij
13 
14  INTERFACE histwrite_phy
15    MODULE PROCEDURE histwrite2d_phy,histwrite3d_phy
16  END INTERFACE
17
18  INTERFACE histbeg_phy_all
19    MODULE PROCEDURE histbeg_phy,histbeg_phy_points
20  END INTERFACE
21
22
23contains
24
25  subroutine init_iophy_new(rlat,rlon)
26  USE dimphy, only: klon
27  USE mod_phys_lmdz_para
28  USE mod_grid_phy_lmdz, only: nbp_lon, nbp_lat, klon_glo
29  USE ioipsl
30  implicit none
31    real,dimension(klon),intent(in) :: rlon ! longitudes, in degrees
32    real,dimension(klon),intent(in) :: rlat ! latitudes, in degrees
33
34    REAL,dimension(klon_glo)        :: rlat_glo
35    REAL,dimension(klon_glo)        :: rlon_glo
36   
37    INTEGER,DIMENSION(2) :: ddid
38    INTEGER,DIMENSION(2) :: dsg
39    INTEGER,DIMENSION(2) :: dsl
40    INTEGER,DIMENSION(2) :: dpf
41    INTEGER,DIMENSION(2) :: dpl
42    INTEGER,DIMENSION(2) :: dhs
43    INTEGER,DIMENSION(2) :: dhe
44    INTEGER :: i   
45
46    CALL gather(rlat,rlat_glo)
47    CALL bcast(rlat_glo)
48    CALL gather(rlon,rlon_glo)
49    CALL bcast(rlon_glo)
50   
51!$OMP MASTER 
52    ALLOCATE(io_lat(nbp_lat))
53    IF (klon_glo == 1) THEN
54      io_lat(1)=rlat_glo(1)
55    ELSE
56      io_lat(1)=rlat_glo(1)
57      io_lat(nbp_lat)=rlat_glo(klon_glo)
58      DO i=2,nbp_lat-1
59        io_lat(i)=rlat_glo(2+(i-2)*nbp_lon)
60      ENDDO
61    ENDIF
62
63    ALLOCATE(io_lon(nbp_lon))
64    IF (klon_glo == 1) THEN
65      io_lon(1)=rlon_glo(1)
66    ELSE
67      io_lon(1:nbp_lon)=rlon_glo(2:nbp_lon+1)
68    ENDIF
69
70! POUR VENUS, LA CARTE EST A L ENVERS !!
71    io_lat = -1.*io_lat
72    io_lon = -1.*io_lon
73
74    ddid=(/ 1,2 /)
75    dsg=(/ nbp_lon, nbp_lat /)
76    dsl=(/ nbp_lon, jj_nb /)
77    dpf=(/ 1,jj_begin /)
78    dpl=(/ nbp_lon, jj_end /)
79    dhs=(/ ii_begin-1,0 /)
80    IF (mpi_rank==mpi_size-1) THEN
81      dhe=(/0,0/)
82    ELSE
83      dhe=(/ nbp_lon-ii_end,0 /) 
84    ENDIF
85   
86    call flio_dom_set(mpi_size,mpi_rank,ddid,dsg,dsl,dpf,dpl,dhs,dhe, &
87                      'APPLE',phys_domain_id)
88
89!$OMP END MASTER
90     
91  end subroutine init_iophy_new
92
93  subroutine init_iophy(lat,lon)
94  USE dimphy
95  USE mod_phys_lmdz_para
96  use ioipsl, only: flio_dom_set
97  USE mod_grid_phy_lmdz, ONLY : nbp_lon, nbp_lat
98  implicit none
99    real,dimension(nbp_lon),intent(in) :: lon
100    real,dimension(nbp_lat),intent(in) :: lat
101
102    INTEGER,DIMENSION(2) :: ddid
103    INTEGER,DIMENSION(2) :: dsg
104    INTEGER,DIMENSION(2) :: dsl
105    INTEGER,DIMENSION(2) :: dpf
106    INTEGER,DIMENSION(2) :: dpl
107    INTEGER,DIMENSION(2) :: dhs
108    INTEGER,DIMENSION(2) :: dhe
109
110!$OMP MASTER 
111    allocate(io_lat(nbp_lat))
112    io_lat(:)=lat(:)
113    allocate(io_lon(nbp_lon))
114    io_lon(:)=lon(:)
115   
116    ddid=(/ 1,2 /)
117    dsg=(/ nbp_lon, nbp_lat /)
118    dsl=(/ nbp_lon, jj_nb /)
119    dpf=(/ 1,jj_begin /)
120    dpl=(/ nbp_lon, jj_end /)
121    dhs=(/ ii_begin-1,0 /)
122    if (mpi_rank==mpi_size-1) then
123      dhe=(/0,0/)
124    else
125      dhe=(/ nbp_lon-ii_end,0 /) 
126    endif
127   
128    call flio_dom_set(mpi_size,mpi_rank,ddid,dsg,dsl,dpf,dpl,dhs,dhe, &
129                      'APPLE',phys_domain_id)
130
131!$OMP END MASTER
132     
133  end subroutine init_iophy
134 
135  subroutine histbeg_phy(name,itau0,zjulian,dtime,nhori,nid_day)
136  USE mod_phys_lmdz_para, only: jj_begin, jj_end, jj_nb, is_sequential
137  USE mod_grid_phy_lmdz, ONLY : nbp_lon, nbp_lat
138  use ioipsl, only: histbeg
139  use write_field
140  implicit none
141   
142    character*(*), intent(IN) :: name
143    integer, intent(in) :: itau0
144    real,intent(in) :: zjulian
145    real,intent(in) :: dtime
146    integer,intent(out) :: nhori
147    integer,intent(out) :: nid_day
148
149!$OMP MASTER   
150    if (is_sequential) then
151      call histbeg(name,nbp_lon,io_lon, jj_nb,io_lat(jj_begin:jj_end), &
152                   1,nbp_lon,1,jj_nb,itau0, zjulian, dtime, nhori, nid_day)
153    else
154      call histbeg(name,nbp_lon,io_lon, jj_nb,io_lat(jj_begin:jj_end), &
155                   1,nbp_lon,1,jj_nb,itau0, zjulian, dtime, nhori, nid_day,phys_domain_id)
156    endif
157!$OMP END MASTER
158 
159  end subroutine histbeg_phy
160
161  subroutine histbeg_phy_points(rlon,rlat,pim,tabij,ipt,jpt, &
162             plon,plat,plon_bounds,plat_bounds, &
163             nname,itau0,zjulian,dtime,nnhori,nnid_day)
164  USE dimphy, only: klon
165  USE mod_phys_lmdz_para
166  USE mod_grid_phy_lmdz, only: klon_glo, nbp_lon, nbp_lat
167  use ioipsl, only: histbeg
168
169  implicit none
170
171    real,dimension(klon),intent(in) :: rlon
172    real,dimension(klon),intent(in) :: rlat
173    integer, intent(in) :: itau0
174    real,intent(in) :: zjulian
175    real,intent(in) :: dtime
176    integer, intent(in) :: pim
177    integer, intent(out) :: nnhori
178    character(len=20), intent(in) :: nname
179    INTEGER, intent(out) :: nnid_day
180    integer :: i
181    REAL,dimension(klon_glo)        :: rlat_glo
182    REAL,dimension(klon_glo)        :: rlon_glo
183    INTEGER, DIMENSION(pim), intent(in)  :: tabij
184    REAL,dimension(pim), intent(in) :: plat, plon
185    INTEGER,dimension(pim), intent(in) :: ipt, jpt
186    REAL,dimension(pim,2), intent(out) :: plat_bounds, plon_bounds
187
188    INTEGER, SAVE :: tabprocbeg, tabprocend
189!$OMP THREADPRIVATE(tabprocbeg, tabprocend)
190    INTEGER :: ip
191    INTEGER, PARAMETER :: nip=1
192    INTEGER :: npproc
193    REAL, allocatable, dimension(:) :: npplat, npplon
194    REAL, allocatable, dimension(:,:) :: npplat_bounds, npplon_bounds
195    REAL, dimension(nbp_lon,nbp_lat) :: zx_lon, zx_lat
196
197    CALL gather(rlat,rlat_glo)
198    CALL bcast(rlat_glo)
199    CALL gather(rlon,rlon_glo)
200    CALL bcast(rlon_glo)
201
202!$OMP MASTER
203    DO i=1,pim
204
205!    print*,'CFMIP_iophy i tabij lon lat',i,tabij(i),plon(i),plat(i)
206
207     plon_bounds(i,1)=rlon_glo(tabij(i)-1)
208     plon_bounds(i,2)=rlon_glo(tabij(i)+1)
209     if(plon_bounds(i,2).LE.0..AND.plon_bounds(i,1).GE.0.) THEN
210      if(rlon_glo(tabij(i)).GE.0.) THEN
211       plon_bounds(i,2)=-1*plon_bounds(i,2)
212      endif
213     endif
214     if(plon_bounds(i,2).GE.0..AND.plon_bounds(i,1).LE.0.) THEN
215      if(rlon_glo(tabij(i)).LE.0.) THEN
216       plon_bounds(i,2)=-1*plon_bounds(i,2)
217      endif
218     endif
219!
220     IF ( tabij(i).LE.nbp_lon) THEN
221      plat_bounds(i,1)=rlat_glo(tabij(i))
222     ELSE
223      plat_bounds(i,1)=rlat_glo(tabij(i)-nbp_lon)
224     ENDIF
225     plat_bounds(i,2)=rlat_glo(tabij(i)+nbp_lon)
226!
227!    print*,'CFMIP_iophy point i lon lon_bds',i,plon_bounds(i,1),rlon_glo(tabij(i)),plon_bounds(i,2)
228!    print*,'CFMIP_iophy point i lat lat_bds',i,plat_bounds(i,1),rlat_glo(tabij(i)),plat_bounds(i,2)
229!
230    ENDDO
231    if (is_sequential) then
232
233     npstn=pim
234     IF(.NOT. ALLOCATED(nptabij)) THEN
235      ALLOCATE(nptabij(pim))
236     ENDIF
237     DO i=1,pim
238      nptabij(i)=tabij(i)
239     ENDDO
240
241       CALL gr_fi_ecrit(1,klon,nbp_lon,nbp_lat,rlon_glo,zx_lon)
242       if ((nbp_lon*nbp_lat).gt.1) then
243       DO i = 1, nbp_lon
244         zx_lon(i,1) = rlon_glo(i+1)
245         zx_lon(i,nbp_lat) = rlon_glo(i+1)
246       ENDDO
247       endif
248       CALL gr_fi_ecrit(1,klon,nbp_lon,nbp_lat,rlat_glo,zx_lat)
249
250    DO i=1,pim
251!    print*,'CFMIP_iophy i tabij lon lat',i,tabij(i),plon(i),plat(i)
252
253     plon_bounds(i,1)=zx_lon(ipt(i)-1,jpt(i))
254     plon_bounds(i,2)=zx_lon(ipt(i)+1,jpt(i))
255
256     if (ipt(i).EQ.1) then
257      plon_bounds(i,1)=zx_lon(nbp_lon,jpt(i))
258      plon_bounds(i,2)=360.+zx_lon(ipt(i)+1,jpt(i))
259     endif
260 
261     if (ipt(i).EQ.nbp_lon) then
262      plon_bounds(i,2)=360.+zx_lon(1,jpt(i))
263     endif
264
265     plat_bounds(i,1)=zx_lat(ipt(i),jpt(i)-1)
266     plat_bounds(i,2)=zx_lat(ipt(i),jpt(i)+1)
267
268     if (jpt(i).EQ.1) then
269      plat_bounds(i,1)=zx_lat(ipt(i),1)+0.001
270      plat_bounds(i,2)=zx_lat(ipt(i),1)-0.001
271     endif
272 
273     if (jpt(i).EQ.nbp_lat) then
274      plat_bounds(i,1)=zx_lat(ipt(i),nbp_lat)+0.001
275      plat_bounds(i,2)=zx_lat(ipt(i),nbp_lat)-0.001
276     endif
277!
278!    print*,'CFMIP_iophy point i lon lon_bds',i,plon_bounds(i,1),rlon(tabij(i)),plon_bounds(i,2)
279!    print*,'CFMIP_iophy point i lat lat_bds',i,plat_bounds(i,1),rlat(tabij(i)),plat_bounds(i,2)
280!
281    ENDDO
282!    print*,'iophy is_sequential nname, nnhori, nnid_day=',trim(nname),nnhori,nnid_day
283     call histbeg(nname,pim,plon,plon_bounds, &
284                           plat,plat_bounds, &
285                           itau0, zjulian, dtime, nnhori, nnid_day)
286    else
287     npproc=0
288     DO ip=1, pim
289      tabprocbeg=klon_mpi_begin
290      tabprocend=klon_mpi_end
291      IF(tabij(ip).GE.tabprocbeg.AND.tabij(ip).LE.tabprocend) THEN
292       npproc=npproc+1
293       npstn=npproc
294      ENDIF
295     ENDDO
296!    print*,'CFMIP_iophy mpi_rank npstn',mpi_rank,npstn
297     IF(.NOT. ALLOCATED(nptabij)) THEN
298      ALLOCATE(nptabij(npstn))
299      ALLOCATE(npplon(npstn), npplat(npstn))
300      ALLOCATE(npplon_bounds(npstn,2), npplat_bounds(npstn,2))
301     ENDIF
302     npproc=0
303     DO ip=1, pim
304      IF(tabij(ip).GE.tabprocbeg.AND.tabij(ip).LE.tabprocend) THEN
305       npproc=npproc+1
306       nptabij(npproc)=tabij(ip)
307!      print*,'mpi_rank npproc ip plon plat tabij=',mpi_rank,npproc,ip, &
308!      plon(ip),plat(ip),tabij(ip)
309       npplon(npproc)=plon(ip)
310       npplat(npproc)=plat(ip)
311       npplon_bounds(npproc,1)=plon_bounds(ip,1)
312       npplon_bounds(npproc,2)=plon_bounds(ip,2)
313       npplat_bounds(npproc,1)=plat_bounds(ip,1)
314       npplat_bounds(npproc,2)=plat_bounds(ip,2)
315!!!
316!!! print qui sert a reordonner les points stations selon l'ordre CFMIP
317!!! ne pas enlever
318        print*,'iophy_mpi rank ip lon lat',mpi_rank,ip,plon(ip),plat(ip)
319!!!
320      ENDIF
321     ENDDO
322     call histbeg(nname,npstn,npplon,npplon_bounds, &
323                            npplat,npplat_bounds, &
324                            itau0,zjulian,dtime,nnhori,nnid_day,phys_domain_id)
325    endif
326!$OMP END MASTER
327
328  end subroutine histbeg_phy_points
329 
330  subroutine histwrite2d_phy(nid,lpoint,name,itau,field)
331  USE dimphy, only: klon
332  USE mod_phys_lmdz_para
333  USE ioipsl, only: histwrite
334  USE mod_grid_phy_lmdz, ONLY : nbp_lon, nbp_lat
335  implicit none
336   
337    integer,intent(in) :: nid
338    logical,intent(in) :: lpoint
339    character*(*), intent(IN) :: name
340    integer, intent(in) :: itau
341    real,dimension(:),intent(in) :: field
342    REAL,dimension(klon_mpi) :: buffer_omp
343    INTEGER, allocatable, dimension(:) :: index2d
344    REAL :: Field2d(nbp_lon,jj_nb)
345
346    integer :: ip
347    real,allocatable,dimension(:) :: fieldok
348
349    IF (size(field)/=klon) CALL abort_physic('iophy::histwrite2d','Field first dimension not equal to klon',1)
350   
351    CALL Gather_omp(field,buffer_omp)   
352!$OMP MASTER
353    CALL grid1Dto2D_mpi(buffer_omp,Field2d)
354    if(.NOT.lpoint) THEN
355     ALLOCATE(index2d(nbp_lon*jj_nb))
356     ALLOCATE(fieldok(nbp_lon*jj_nb))
357     CALL histwrite(nid,name,itau,Field2d,nbp_lon*jj_nb,index2d)
358    else
359     ALLOCATE(fieldok(npstn))
360     ALLOCATE(index2d(npstn))
361
362     if(is_sequential) then
363!     klon_mpi_begin=1
364!     klon_mpi_end=klon
365      DO ip=1, npstn
366       fieldok(ip)=buffer_omp(nptabij(ip))
367      ENDDO
368     else
369      DO ip=1, npstn
370!     print*,'histwrite2d is_sequential npstn ip name nptabij',npstn,ip,name,nptabij(ip)
371       IF(nptabij(ip).GE.klon_mpi_begin.AND. &
372          nptabij(ip).LE.klon_mpi_end) THEN
373         fieldok(ip)=buffer_omp(nptabij(ip)-klon_mpi_begin+1)
374       ENDIF
375      ENDDO
376     endif
377     CALL histwrite(nid,name,itau,fieldok,npstn,index2d)
378!
379    endif
380    deallocate(index2d)
381    deallocate(fieldok)
382!$OMP END MASTER   
383  end subroutine histwrite2d_phy
384
385  subroutine histwrite3d_phy(nid,lpoint,name,itau,field)
386  USE dimphy, only: klon
387  USE mod_phys_lmdz_para
388  USE mod_grid_phy_lmdz, ONLY : nbp_lon, nbp_lat
389  use ioipsl, only: histwrite
390  implicit none
391   
392    integer,intent(in) :: nid
393    logical,intent(in) :: lpoint
394    character*(*), intent(IN) :: name
395    integer, intent(in) :: itau
396    real,dimension(:,:),intent(in) :: field  ! --> field(klon,:)
397    REAL,dimension(klon_mpi,size(field,2)) :: buffer_omp
398    REAL :: Field3d(nbp_lon,jj_nb,size(field,2))
399    INTEGER :: ip, n, nlev
400    INTEGER, ALLOCATABLE, dimension(:) :: index3d
401    real,allocatable, dimension(:,:) :: fieldok
402
403    IF (size(field,1)/=klon) CALL abort_physic('iophy::histwrite3d','Field first dimension not equal to klon',1)
404    nlev=size(field,2)
405
406!   print*,'hist3d_phy mpi_rank npstn=',mpi_rank,npstn
407
408!   DO ip=1, npstn
409!    print*,'hist3d_phy mpi_rank nptabij',mpi_rank,nptabij(ip)
410!   ENDDO
411
412    CALL Gather_omp(field,buffer_omp)
413!$OMP MASTER
414    CALL grid1Dto2D_mpi(buffer_omp,field3d)
415    if(.NOT.lpoint) THEN
416     ALLOCATE(index3d(nbp_lon*jj_nb*nlev))
417     ALLOCATE(fieldok(nbp_lon*jj_nb,nlev))
418     CALL histwrite(nid,name,itau,Field3d,nbp_lon*jj_nb*nlev,index3d)
419    else
420      nlev=size(field,2)
421      ALLOCATE(index3d(npstn*nlev))
422      ALLOCATE(fieldok(npstn,nlev))
423
424      if(is_sequential) then
425!      klon_mpi_begin=1
426!      klon_mpi_end=klon
427       DO n=1, nlev
428       DO ip=1, npstn
429        fieldok(ip,n)=buffer_omp(nptabij(ip),n)
430       ENDDO
431       ENDDO
432      else
433       DO n=1, nlev
434       DO ip=1, npstn
435        IF(nptabij(ip).GE.klon_mpi_begin.AND. &
436         nptabij(ip).LE.klon_mpi_end) THEN
437         fieldok(ip,n)=buffer_omp(nptabij(ip)-klon_mpi_begin+1,n)
438        ENDIF
439       ENDDO
440       ENDDO
441      endif
442      CALL histwrite(nid,name,itau,fieldok,npstn*nlev,index3d)
443    endif
444  deallocate(index3d)
445  deallocate(fieldok)
446!$OMP END MASTER   
447  end subroutine histwrite3d_phy
448 
449end module iophy
Note: See TracBrowser for help on using the repository browser.