source: trunk/LMDZ.TITAN.old/libf/phytitan/iophy.F90 @ 3557

Last change on this file since 3557 was 1545, checked in by emillour, 9 years ago

Venus and Titan GCMs:

Adaptation wrt previous changes for Titan and Venus where
longitude and latitude arrays (in phycommon/geometry_mod) were overwritten
with values from startphy.nc files, where values are given in degrees.
For the sake of homegeneity with other physics package, revert to "default"
behaviour: longitude/latitude are in radians and longitude_deg/latitude_deg
are in degrees.
Also added checking in phyetat0 that the longitude/latitude read in the
restartphy.nc files match the ones provided by the dynamics.

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