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

1           SUBROUTINE TREADG_nemsio(gfilename,FHOUR,IDATE,
2          &                   zsg,psg,ttg,uug,vvg,rqg,
3          &                   pdryini,IPRINT,
4          &                   global_lats_a,lats_nodes_a,lonsperlat)
5     !
6     !-------------------------------------------------------------------
7     !*** program log
8     !*** Dec, 2009 Jun Wang:  read grid point variables for restart
9     !-------------------------------------------------------------------
10     !
11     !!
12           use gfs_dyn_resol_def
13           use gfs_dyn_layout1
14           use gfs_dyn_coordinate_def					
15           use gfs_dyn_io_header
16           use namelist_dynamics_def
17           use gfs_dyn_vert_def
18           use gfs_dyn_mpi_def
19           use gfs_dyn_physcons, rerth => con_rerth
20          &,             grav => con_g, rkap => con_rocp
21          &,             cpd => con_cp
22           use module_nemsio
23           use gfs_dyn_tracer_config, ONLY : gfs_dyn_tracer   ! generalized tracer
24     !
25           IMPLICIT NONE
26           character(*) gfilename
27           REAL(KIND=KIND_EVOD) FHOUR
28           INTEGER              IDATE(4),NTRACI, ntozi, ntcwi, ncldi
29          &,                    latbi, lonbi, levsi, jcapi,
30          &                     latgi, lonfi, latri, lonri,idate7(7)
31     !
32           integer              lats_nodes_a(nodes)
33           INTEGER              IPRINT
34           INTEGER              J,K,L,LOCL,N,lv,kk
35           integer              i,lan,lat,iblk,lons_lat,il,lon,njeff,nn
36           REAL(KIND=KIND_EVOD) TRUN,WAVES,XLAYERS
37           REAL(KIND=KIND_EVOD) XI(LEVP1),XL(LEVS)
38           REAL(KIND=KIND_EVOD) sikp1(levp1)
39           REAL(KIND=KIND_IO4)   VTID,RUNID4,fhour4,XNCLD,xgf
40           REAL(KIND=KIND_grid)  PDRYINI
41           real(kind=kind_io4), allocatable ::  vcoord4(:,:,:)
42     ! for generalized tracers
43           integer                      nreci
44           character*16               :: vname
45     !
46           type (nemsio_gfile) gfile_in
47     !
48           integer              iret, num_dta
49           integer              im,jm,fieldsize
50           real(kind=kind_evod) psurfff
51           real(kind=kind_evod) pressk, tem
52           real(kind=kind_evod), parameter :: rkapi=1.0/rkap,
53          &                                   rkapp1=1.0+rkap
54     !
55           integer kmsk(lonf,latg), global_lats_a(latg), lonsperlat(latg)
56           real(kind=kind_io8) buffo(lonf,lats_node_a)
57          &,                   buff2(lonf,lats_node_a)
58           real(kind=kind_evod) teref(levp1),ck5p(levp1)			! hmhj
59     !
60           real (kind=kind_io8), allocatable ::  nemsio_data(:)
61     !!
62           real(kind=kind_grid) zsg(lonf,lats_node_a)
63           real(kind=kind_grid) psg(lonf,lats_node_a)
64           real(kind=kind_grid) uug(lonf,lats_node_a,levs)
65           real(kind=kind_grid) vvg(lonf,lats_node_a,levs)
66           real(kind=kind_grid) ttg(lonf,lats_node_a,levs)
67           real(kind=kind_grid) rqg(lonf,lats_node_a,levh)
68     !
69     !------------------------------------------------------------------
70           print *,'in treadg_nemsio,gfilename=',trim(gfilename)
71     !
72     !--- Input file is in grid-point space - use nems_io package
73     !
74           call nemsio_init()
75     !
76           call nemsio_open(gfile_in,trim(gfilename),'read',iret)
77           print *,'after nemsio_open,iret=',iret
78     !
79           call nemsio_getfilehead(gfile_in,iret=iret,
80          &  idate=idate7,
81          &  dimy=im,dimx=jm,dimz=levsi,jcap=jcapi,
82          &  ntrac=ntraci,idsl=idsl,idvc=idvc,idvm=idvm,
83          &  ncldt=ncldt)
84     !      print *,'after nemsio_getfilehead,iret=',iret
85     !!
86           call nemsio_getheadvar(gfile_in,'fhour',fhour,iret=iret)
87           call nemsio_getheadvar(gfile_in,'iorder',iorder,iret=iret)
88           call nemsio_getheadvar(gfile_in,'irealf',irealf,iret=iret)
89           call nemsio_getheadvar(gfile_in,'igen',igen,iret=iret)
90           call nemsio_getheadvar(gfile_in,'latf',latgi,iret=iret)
91           call nemsio_getheadvar(gfile_in,'lonf',lonfi,iret=iret)
92           call nemsio_getheadvar(gfile_in,'icen2',icen2,iret=iret)
93           call nemsio_getheadvar(gfile_in,'iens',iens,iret=iret)
94           call nemsio_getheadvar(gfile_in,'idpp',idpp,iret=iret)
95           call nemsio_getheadvar(gfile_in,'idvt',idvt,iret=iret)
96           call nemsio_getheadvar(gfile_in,'idrun',idrun,iret=iret)
97           call nemsio_getheadvar(gfile_in,'itrun',itrun,iret=iret)
98           call nemsio_getheadvar(gfile_in,'nvcoord',nvcoord,iret=iret)
99           call nemsio_getheadvar(gfile_in,'pdryini',pdryini,iret=iret)
100     !
101           if (me == 0) then
102             print *,'iret=',iret,'idvt=',idvt,' nvcoord=',nvcoord,
103          &     ' levsi=',levsi,'ntoz=',ntoz,'idate7=',idate7,
104          &   'lonf=',lonf,'lonfi=',lonfi,'latg=',latg,'latgi=',latgi,
105          &   'jcap=',jcap,'jcapi=',jcapi,'levs=',levs,'levsi=',levsi,
106          &   'idvc=',idvc,'idvm=',idvm,'idsl=',idsl,'itrun=',itrun,
107          &   'gen_coord_hybrid=',gen_coord_hybrid,'pdryini=',pdryini
108             if(lonf .ne. lonfi .or. latg .ne. latgi .or.
109          &     jcap .ne. jcapi .or. levs .ne. levsi) then
110               print *,' Input resolution and the model resolutions are'
111          &,  ' different- run aborted'
112               call mpi_quit(777)
113             endif
114           endif
115     !
116           allocate (vcoord4(levsi+1,3,2))
117           if(.not.allocated(vcoord)) allocate (vcoord(levsi+1,nvcoord))
118           call nemsio_getfilehead(gfile_in,iret=iret,vcoord=vcoord4)
119     !
120           vcoord(:,1:nvcoord) = vcoord4(:,1:nvcoord,1)
121     !      if (me .eq. 0) print *,' vcoord=',vcoord(1:10,1:nvcoord)
122           deallocate (vcoord4)
123     !
124     !---  for generalized tracers
125     ! retrieve nreci, recnamei, reclevtypi, and reclevi
126           call nemsio_getfilehead(gfile_in,iret=iret,nrec=nreci)
127           if (me == 0) then
128             print *, 'LU_TRC: nreci =', nreci, iret
129           endif
130     !
131     !---
132           if (gen_coord_hybrid) then                                        ! hmhj
133     
134             sfcpress_id  = mod(idvm , 10)
135             thermodyn_id = mod(idvm/10 , 10)
136     !   ak bk ck in file have the same order as model                       ! hmhj
137             do k=1,levp1                                                    ! hmhj
138               ak5(k) = vcoord(k,1)/1000.                                    ! hmhj
139               bk5(k) = vcoord(k,2)                                          ! hmhj
140               ck5(k) = vcoord(k,3)/1000.                                    ! hmhj
141             enddo                                                           ! hmhj
142             vertcoord_id=0                                                  ! hmhj
143             do k=1,levp1                                                    ! hmhj
144               if( ck5(k).ne.0.0 ) vertcoord_id=3                            ! hmhj
145             enddo
146     ! provide better estimated press                                        ! hmhj
147             psurfff = 101.3d0
148             if( thermodyn_id.eq.3 ) then                                    ! hmhj
149               do k=1,levs                                                   ! hmhj
150                 thref(k) = 300.0*cpd                                        ! hmhj
151                 teref(k) = 255.0*cpd                                        ! hmhj
152               enddo                                                         ! hmhj
153             else                                                            ! hmhj
154              do k=1,levp1                                                   ! hmhj
155               thref(k) = 300.0                                              ! hmhj
156               teref(k) = 255.0                                              ! hmhj
157              enddo                                                          ! hmhj
158             endif
159             ck5p(levp1) = ck5(levp1)                                        ! hmhj
160             do k=1,levs                                                     ! hmhj
161               ck5p(k) = ck5(k)*(teref(k)/thref(k))**rkapi                   ! hmhj
162             enddo
163             if( me.eq.0 ) then                                              ! hmhj
164               do k=1,levp1                                                  ! hmhj
165                 pressk=ak5(k)+bk5(k)*psurfff+ck5p(k)                        ! hmhj
166                 print 180,k,ak5(k),bk5(k),ck5(k),pressk                     ! hmhj
167      180     format('k=',i2,'  ak5=',f13.6,'  bk5=',e13.5,                  ! hmhj
168          &            '   ck5=',f13.6,'  closed pressk=',f10.6)             ! hmhj
169               enddo                                                         ! hmhj
170             endif                                                           ! hmhj
171             do k=1,levp1                                                    ! hmhj
172               si(k) = ak5(k)/psurfff + bk5(k) + ck5p(k)/psurfff             ! hmhj
173             enddo                                                           ! hmhj
174             do k=1,levs                                                     ! hmhj
175               sl(k) = 0.5*(si(k)+si(k+1))                                   ! hmhj
176             enddo                                                           ! hmhj
177     
178           else if (hybrid .and. idvc .eq. 2) then
179     !       idsl=slid  !=2,pk=0.5*(p(k+1/2)+p(k-1/2)) check alfa(1)  am_bm
180     !   ak bk order in "sigma" file is bottom to top !!!!!!!!!!!!!!!!!!
181             psurfff = 101.3
182             do k=1,levp1
183               ak5(k) = vcoord(levp1+1-k,1)/1000.
184               bk5(k) = vcoord(levp1+1-k,2)
185               pressk = ak5(k) + bk5(k)*psurfff
186     
187               if(me.eq.0)print 190,k,ak5(k),bk5(k),pressk
188     190       format('k=',i2,'  ak5=',E14.6,'  bk5=',e14.6,
189          &           '  pressk=',E14.6)
190     
191             enddo
192             do k=1,levs
193               dbk(k) = bk5(k+1)-bk5(k)
194               bkl(k) = (bk5(k+1)+bk5(k))*0.5
195               ck(k)  = ak5(k+1)*bk5(k)-ak5(k)*bk5(k+1)
196               if(me.eq.0)print 200,k,dbk(k),ck(k)
197     200       format('k=',i2,'  dbk=',f8.6,'  ck=',e13.5)
198             enddo
199     !
200     ! hmhj give an estimated si and sl for dynamics
201             do k=1,levs+1
202               si(levs+2-k) = ak5(k)/psurfff + bk5(k) !ak(k) bk(k) go top to bottom
203             enddo
204             do k=1,levs
205               sl(k) = 0.5*(si(k)+si(k+1))
206             enddo
207     !
208           elseif (idvc .le. 1) then
209             si(:)    = vcoord(:,1)
210             sik(:)   = si(:) ** rkap
211             sikp1(:) = si(:) ** rkapp1
212             do k=1,levs
213               tem      = rkapp1 * (si(k) - si(k+1))
214               slk(k)   = (sikp1(k)-sikp1(k+1))/tem
215               sl(k)    = slk(k) ** rkapi
216     !         sl(k)    = ((sikp1(k)-sikp1(k+1))/tem)**rkapi
217               if (me .eq. 0) print 250, k, si(k), sl(k)
218     250       format('k=',i2,'  si=',f8.6,'  sl=',e13.5)
219             enddo
220           else
221             print *,' Non compatible Initial state IDVC=',idvc
222          &,' iret=',iret
223             call MPI_QUIT(333)
224           endif
225     !
226           idate(1)       = idate7(4)
227           idate(2:3)       = idate7(2:3)
228           idate(4)       = idate7(1)
229           WAVES       = jcap
230           XLAYERS     = levs
231           itrun       = itrun
232           icen        = 7
233           ienst       = iens(1)
234           iensi       = iens(2)
235           ntraci = ntrac
236           if (idvt .gt. 0.0) then
237             ntcwi = idvt / 10
238             ntozi = idvt - ntcwi * 10 + 1
239             ntcwi = ntcwi + 1
240             ncldi = ncldt
241           elseif(ntraci .eq. 2) then
242             ntozi = 2
243             ntcwi = 0
244             ncldi = 0
245           elseif(ntraci .eq. 3) then
246             ntozi = 2
247             ntcwi = 3
248             ncldi = 1
249           else
250             ntozi = 0
251             ntcwi = 0
252             ncldi = 0
253           endif
254     
255     !
256     !
257           IF (me.eq.0) THEN
258             write(0,*)'gfile,in treadeo fhour,idate=',gfilename,fhour,idate
259          &, ' ntozi=',ntozi,' ntcwi=',ntcwi,' ncldi=',ncldi
260          &, ' ntraci=',ntraci,' tracers=',ntrac,' idvt=',idvt
261          &, ' ncldt=',ncldt,' idvc=',idvc,' jcap=',jcap
262          &, ' pdryini=',pdryini,'ntoz=',ntoz,'idate=',idate
263           ENDIF
264     !
265           fieldsize=im*jm
266           allocate (nemsio_data(fieldsize))
267     !  Read orog
268           call nemsio_readrecv(gfile_in,'hgt','sfc',1,nemsio_data,
269          &     iret=iret)
270           call split2d_rdGRD(nemsio_data,zsg,fieldsize,global_lats_a,
271          &     lonsperlat)
272     !
273     !  Read ps
274           call nemsio_readrecv(gfile_in,'pres','sfc',1,nemsio_data,
275          &     iret=iret)
276     !      print *,'aftr read in pressfc=',maxval(nemsio_data),
277     !     &    minval(nemsio_data)
278           call split2d_rdGRD(nemsio_data,psg,fieldsize,global_lats_a,
279          &  lonsperlat)
280     !      print *,'in treadg_nemsio,psg=',maxval(psg(1:lonf,1:lats_node_a)),
281     !     &    minval(psg(1:lonf,1:lats_node_a)),'psg(1:5,1)=',psg(1:5,1),
282     !     & 'psg(1:5,2)=',psg(1:5,2),'psg(1:5,3)=',psg(1:5,3)
283     !
284     !  Read u
285           do k=1,levs
286             call nemsio_readrecv(gfile_in,'ugrd','mid layer',k,nemsio_data,
287          &       iret=iret)
288             call split2d_rdGRD(nemsio_data,uug(:,:,k),fieldsize,
289          &    global_lats_a,lonsperlat)
290           enddo
291     !  Read v
292           do k=1,levs
293             call nemsio_readrecv(gfile_in,'vgrd','mid layer',k,nemsio_data,
294          &    iret=iret)
295             call split2d_rdGRD(nemsio_data,vvg(:,:,k),fieldsize,
296          &    global_lats_a,lonsperlat)
297           enddo
298     !  Read T   -- this is real temperature
299           do k=1,levs
300             call nemsio_readrecv(gfile_in,'tmp','mid layer',k,nemsio_data,
301          &    iret=iret)
302             call split2d_rdGRD(nemsio_data,ttg(:,:,k),fieldsize,
303          &    global_lats_a,lonsperlat)
304           enddo
305     !
306     !  Initial Tracers with zero
307     !
308           rqg(:,:,:) = 0.0
309     
310     !! Generalized tracers: 
311     !! Loop through ntrac to read in met + chem tracers
312     !*
313           do k=1,levh
314               call nemsio_readrecv(gfile_in,'tracer',
315          &                   'tracer layer',k,nemsio_data,iret=iret)
316               if(iret == 0) then
317     !            if(me==0) print *,'TRACER read in ok -',
318     !     &                gfs_dyn_tracer%vname(n),k
319     !      print *,'aftr read in tmp=',maxval(nemsio_data),
320     !     &    minval(nemsio_data),'k=',k
321                 call split2d_rdGRD(nemsio_data,rqg(:,:,k),fieldsize,
322          &      global_lats_a,lonsperlat)
323     !      print *,'in treadg_nemsio,k=',k,'rqg=',
324     !     &    maxval(rqg(1:lonf,1:lats_node_a,k)),
325     !     &    minval(rqg(1:lonf,1:lats_node_a,k))
326               else
327                 if(me==0) print *,'TRACER: tracer not found in input; ',
328          &         'set chem tracer to default values',me,k
329               endif
330            enddo
331     !
332           call nemsio_close(gfile_in,iret)
333           call nemsio_finalize()
334     !
335           deallocate(nemsio_data)
336           iprint=0
337      
338     !!!!
339           RETURN
340           END
341