File: C:\NOAA\NEMS_11731\src\atmos\gfs\dyn\twriteg_rst.f

1           subroutine twriteg_rst(fname,IOPROC,fhour,idate,
2          x           si,pdryini,global_lats_a,lonsperlat,lats_nodes_a,
3          x           psg,ttg,uug,vvg,rqg,zsg,
4          &           kdt,nfcstdate7)
5     !
6     !-------------------------------------------------------------------
7     !*** program log
8     !*** Dec, 2009 Jun Wang:  write spectral variables for restart
9     !-------------------------------------------------------------------
10     !
11           use gfs_dyn_resol_def
12           use gfs_dyn_layout1
13           use gfs_dyn_io_header
14           use gfs_dyn_coordinate_def	
15           use namelist_dynamics_def
16           use gfs_dyn_mpi_def
17           use module_nemsio
18     !
19           implicit none
20     !
21           character(*),intent(in) :: fname
22           integer,intent(in) :: ioproc
23           real(kind=kind_evod),intent(in) :: fhour
24           integer,intent(in) :: idate(4)
25           integer,intent(in) :: kdt,nfcstdate7(7)
26     !
27           real(kind=kind_evod) si(levp1)
28           REAL(KIND=KIND_IO8),intent(in) :: pdryini
29           integer,intent(in) ::          global_lats_a(latg)
30           integer,intent(in) ::          lonsperlat(latg)
31           integer,intent(in) ::          lats_nodes_a(nodes)
32     !
33           REAL(KIND=KIND_GRID),intent(in) :: zsg(lonf,lats_node_a_max)
34           REAL(KIND=KIND_GRID),intent(in) :: psg(lonf,lats_node_a_max)
35           REAL(KIND=KIND_GRID),intent(in) :: ttg(lonf,lats_node_a_max,levs)
36           REAL(KIND=KIND_GRID),intent(in) :: uug(lonf,lats_node_a_max,levs)
37           REAL(KIND=KIND_GRID),intent(in) :: vvg(lonf,lats_node_a_max,levs)
38           REAL(KIND=KIND_GRID),intent(in) :: rqg(lonf,lats_node_a_max,levh)
39     !
40     !local variables:
41           REAL(kind=8) t1,t2,t3,t4,t5,t6,ta,tb,rtc
42     !
43           integer              ierr,j,k,l,lenrec,locl,n,node,jrec,idate7(7)
44     !
45           real(kind=kind_ior),allocatable ::  tmp8(:,:),buf(:)
46           real(kind=kind_ior),allocatable ::  GZ(:)
47     !
48           type(nemsio_gfile) gfile
49           integer kps,ktt,kuu,kvv,krq
50           integer nfhour,nfminute,nfsecondn,nfsecondd,nrec,nmeta
51           integer nmetavari,nmetavarr8,nmetaaryi,fieldsize
52           character(16),allocatable :: recname(:),reclevtyp(:)
53           character(16),allocatable :: variname(:),varr8name(:),aryiname(:)
54           integer,allocatable :: reclev(:)
55           integer,allocatable :: varival(:),aryilen(:),aryival(:,:)
56           real(kind=kind_io4),allocatable :: vcoord4(:,:,:)
57           real(8),allocatable :: varr8val(:)
58           integer iret, ipt_lats
59           integer  il,ilen,i,msgtag,ls_diml,nodesl,ij
60           logical first
61           save first,nmetavari,nmetavarr8,nmetaaryi,recname,reclevtyp, 
62          &     reclev,variname,varr8name,aryiname,varival,varr8val,
63          &     aryilen,aryival,vcoord4,nmeta,nrec,GZ
64           data first /.true./
65     !
66     !
67           real(kind=kind_mpi_r),allocatable :: grid_node (:,:,:)
68           real(kind=kind_mpi_r),allocatable :: grid_nodes(:,:,:,:)
69     !
70           real(kind=kind_mpi_r),allocatable :: grid_gz_nodes(:,:,:)
71     !
72           integer      lan,lat,iblk,lons_lat,lon,NJEFF,nn,lv
73     !
74     !---------------------------------------------------------------------
75     !
76     !       print *,' enter twriteg_rst ' 
77     !       print *,'lonf=',lonf,'lats_node_a_max=',lats_node_a_max,
78     !     &  'total_levels=',3*levs+1*levh+1,'lonsperlat=',lonsperlat
79     
80           allocate ( grid_node ( lonf,lats_node_a_max,3*levs+1*levh+1 ) )
81     !
82           fieldsize=sum(lonsperlat)
83     !      print *,'fieldsize=',fieldsize
84     !
85     !collect data 
86           kps=1
87           ktt=kps+1
88           kuu=ktt+levs
89           kvv=kuu+levs
90           krq=kvv+levs
91     !
92           do j=1,lats_node_a_max
93             do i=1,lonf
94               grid_node(i,j,kps) = psg(i,j)
95             enddo
96           enddo
97     !
98           do k=1,levs
99             do j=1,lats_node_a_max
100               do i=1,lonf
101                 grid_node(i,j,ktt+k-1) = ttg(i,j,k)
102                 grid_node(i,j,kuu+k-1) = uug(i,j,k)
103                 grid_node(i,j,kvv+k-1) = vvg(i,j,k)
104               enddo
105             enddo
106           enddo
107     !
108           do k=1,levh
109             do j=1,lats_node_a_max
110               do i=1,lonf
111                 grid_node(i,j,krq+k-1) = rqg(i,j,k)
112               enddo
113             enddo
114           enddo
115     !
116     !      print *,'pdryini=',pdryini,'lats_nodes_a=',lats_nodes_a, 
117     !     &  'nodes=',nodes
118     !
119           if ( me .eq. ioproc ) then
120              allocate ( grid_nodes ( lonf,lats_node_a_max,
121          &                           3*levs+1*levh+1, nodes ),stat=ierr )
122              if(first) then
123                allocate ( grid_gz_nodes (lonf,lats_node_a_max,nodes),
124          &                stat=ierr )
125              endif
126              if (ierr .ne. 0) then
127                call mpi_abort(mpi_comm_all,ierr,i)
128              endif
129           endif
130     !
131           if(nodes>1) then
132             lenrec = lonf*lats_node_a_max * (3*levs+1*levh+1)
133     !
134             t1=rtc()
135     !        print *,'after allocate grid_nodes,lenrec=',lenrec
136             call mpi_gather( grid_node , lenrec, MPI_R_MPI_R,
137          x                 grid_nodes, lenrec, MPI_R_MPI_R,
138          x                 ioproc, MPI_COMM_ALL, ierr)
139           else
140             grid_nodes(:,:,:,1)=grid_node(:,:,:)
141           endif
142           deallocate(grid_node)
143           if(first) then
144             if(nodes>1) then
145               lenrec=lonf*lats_node_a_max 
146               call mpi_gather( zsg , lenrec, MPI_R_MPI_R,
147          x                 grid_gz_nodes, lenrec, MPI_R_MPI_R,
148          x                 ioproc, MPI_COMM_ALL, ierr)
149             else
150               grid_gz_nodes(:,:,1)=zsg(:,:)
151             endif
152             if ( me .eq. ioproc ) then
153     !
154             allocate ( tmp8 ( lonf,latg ) )
155             allocate ( gz(fieldsize) )
156             ipt_lats=1
157             do node=1,nodes
158                do j=1,lats_nodes_a(node)
159                 lat=global_lats_a(ipt_lats-1+j)
160                 do i=1,lonf
161                   tmp8(i,lat)=grid_gz_nodes(i,j,node)
162                 enddo
163               enddo
164               ipt_lats=ipt_lats+lats_nodes_a(node)
165             enddo
166             ij=0
167             do j=1,latg
168                 do i=1,lonsperlat(j)
169                   ij=ij+1
170                   GZ(ij) = tmp8(i,j)
171                 enddo
172             enddo
173             deallocate(tmp8)
174             deallocate ( grid_gz_nodes )
175     !         print *,'gz=',maxval(gz),minval(gz)
176            endif
177           endif
178     !
179           t2=rtc()
180     !
181           IF (me.eq.ioproc) THEN
182      
183     !        print *,' in TWRITEG fhour=',fhour,'ij=',ij,'fieldsize=',
184     !     &    fieldsize
185     !
186             if (first) then
187     !
188               nmeta=6
189               nrec=3*levs+1*levh+2
190               allocate(vcoord4(levp1,3,2))
191               vcoord4=0.
192               idvt    = (ntoz-1) + 10 * (ntcw-1)
193     
194               nmetavari=16
195               allocate(variname(nmetavari),varival(nmetavari))
196               variname(1:nmetavari)=(/"latb     ","lonb     ","itrun    ", 
197          &   "iorder   ","irealf   ","igen     ","latf     ","lonf     ",
198          &   "icen2    ","idpp     ","idvt     ","idrun    ","idusr    ",
199          &   "ivs      ","nvcoord  ","NTIMESTEP"/)
200               varival(1:nmetavari-1)=(/latb,lonb,itrun,2,2,igen,
201          &    latg,lonf,icen2,0,idvt,0,0,200509,nvcoord /)
202     
203               nmetaaryi=2
204               allocate(aryiname(nmetaaryi),aryilen(nmetaaryi))
205               aryiname(1)='iens'
206               aryilen(1)=2
207               aryiname(2)='FCSTDATE'
208               aryilen(2)=7
209               allocate(aryival(maxval(aryilen),nmetaaryi))
210               aryival(1,1)=ienst
211               aryival(2,1)=iensi
212     !!
213               if (gen_coord_hybrid) then					! hmhj
214     
215     !            idvc    = vctype						! hmhj
216                 idvc    = 3      					! hmhj
217                 idvm    = 32    ! 1: ln(ps) 2:ps				! hmhj
218                 idsl    = 2    ! idsl=2 for middle of layer		! hmhj
219                 do k=1,levp1							! hmhj
220                   vcoord4(k,1,1)=ak5(k)*1000.				! hmhj
221                   vcoord4(k,2,1)=bk5(k)					! hmhj
222                   vcoord4(k,3,1)=ck5(k)*1000.				! hmhj
223                 enddo								! hmhj
224     
225               else if (hybrid) then						! hmhj
226                 idvc    = 2    ! for hybrid vertical coord.
227                 do k=1,levp1
228                   vcoord4(k,1,1)=ak5(levp1+1-k)*1000.
229                   vcoord4(k,2,1)=bk5(levp1+1-k)
230     !           if(me.eq.0)print 190,k,head%vcoord(k,1),head%vcoord(k,2)
231     190         format('in twrite k=',i2,'  ak5r4=',f13.6,'  bk5r4=',e13.5)
232                 enddo
233               else
234                 idvc    = 1    ! for sigma vertical coord. (default)
235                 vcoord4(:,1,1) = si (:)
236               endif
237     !
238     !-- field infomation
239               allocate(recname(nrec),reclevtyp(nrec),reclev(nrec))
240               recname(1)='hgt'
241               recname(2)='pres'
242               recname(3:levs+2)='tmp'
243               recname(levs+3:2*levs+2)='ugrd'
244               recname(2*levs+3:3*levs+2)='vgrd'
245               recname(3*levs+3:3*levs+levh+2)='tracer'
246               reclevtyp(1)='sfc'
247               reclevtyp(2)='sfc'
248               reclevtyp(3:3*levs+2)='mid layer'
249               reclevtyp(3*levs+3:3*levs+levh+2)='tracer layer'
250               reclev(1)=1
251               reclev(2)=1
252               do i=1,levs
253                 reclev(i+2)=i
254                 reclev(i+2+levs)=i
255                 reclev(i+2+2*levs)=i
256               enddo
257               do i=1,levh
258                 reclev(i+2+3*levs)=i
259               enddo
260     !
261               nmetavarr8=2
262               allocate(varr8name(nmetavarr8),varr8val(nmetavarr8))
263               varr8name(1:nmetavarr8)=(/'pdryini','fhour  '/)
264     !
265     !endof first
266             endif
267     !
268             idate7=0;idate7(7)=1
269             idate7(1)=idate(4)
270             idate7(2:3)=idate(2:3)
271             idate7(4)=idate(1)
272     !
273             nfhour=int(fhour)
274             nfminute=int((fhour-nfhour)*60)
275             nfsecondn=int(((fhour-nfhour)*60-nfminute)*60)
276             nfsecondd=1
277     !
278             varr8val(1)=pdryini
279             varr8val(2)=fhour
280             varival(nmetavari)=kdt
281             aryival(1:7,2)=nfcstdate7(1:7)
282     !
283             call nemsio_open(gfile,fname,'write',iret,modelname='GFS',   
284          &    gdatatype='bin8',idate=idate7,nfhour=nfhour,nfminute=nfminute,
285          &    nfsecondn=nfsecondn,nfsecondd=nfsecondd,dimx=fieldsize,dimy=1,
286          &    dimz=levs,nmeta=nmeta,jcap=jcap,idsl=idsl,
287          &    idvm=idvm,idvc=idvc,ntrac=ntrac,nrec=nrec,ncldt=ncld,
288          &    recname=recname,reclevtyp=reclevtyp,reclev=reclev,
289          &    vcoord=vcoord4,extrameta=.true.,nmetavari=nmetavari,
290          &    nmetavarr8=nmetavarr8,nmetaaryi=nmetaaryi,
291          &    variname=variname,varival=varival,varr8name=varr8name,
292          &    varr8val=varr8val,aryiname=aryiname,aryilen=aryilen,
293          &    aryival=aryival)
294     !
295     !--- write out data
296     !
297            jrec=1
298            call nemsio_writerec(gfile,1,GZ,iret=iret)
299     !
300            allocate ( tmp8 ( lonf,latg ) )
301            allocate ( buf(fieldsize) )
302            do k=1,3*levs+1*levh+1
303              jrec=k+1
304              ipt_lats=1
305              do node=1,nodes
306                do j=1,lats_nodes_a(node)
307                 lat=global_lats_a(ipt_lats-1+j)
308                 do i=1,lonf
309                   tmp8(i,lat)=grid_nodes(i,j,k,node)
310                 enddo
311                enddo
312                ipt_lats=ipt_lats+lats_nodes_a(node)
313              enddo
314              ij=0
315              do j=1,latg
316                 do i=1,lonsperlat(j)
317                   ij=ij+1
318                   buf(ij) = tmp8(i,j)
319                 enddo
320              enddo
321              call nemsio_writerec(gfile,jrec,buf,iret=iret)
322            end do
323            deallocate(buf,tmp8)
324            deallocate(grid_nodes)
325            call nemsio_close(gfile,iret=iret)
326     !
327             t4=rtc  ()
328     !sela print *, ' DISK TIME FOR SIG TWRITEO WRT ',t4-t3
329     !
330           endif   !me.eq.ioproc
331     !!
332           if(first) then
333               first = .false.
334           endif
335     !
336           call mpi_barrier(MPI_COMM_ALL,ierr)
337     !      print *,' leave twrites_nemsio ' 
338     
339           return
340           end
341