NCEPLIBS-nemsio  2.5.3
•All Data Structures Files
mkgfsnemsioctl.f90
Go to the documentation of this file.
1 
6  program main
7  use nemsio_module
8  implicit none
9 !
10 !---------------------------------------------------------------------------
11  type(nemsio_gfile) :: gfile
12  character(255) cin,cindx
13  real,allocatable :: data(:,:)
14 !---------------------------------------------------------------------------
15 !--- nemsio meta data
16  integer nrec,im,jm,lm,idate(7),nfhour,tlmeta,nsoil,fieldsize
17  real(4),allocatable :: lat1(:),lat(:),lon1(:)
18  real(8),allocatable :: slat(:)
19  character(16),allocatable:: recname(:),reclevtyp(:)
20  integer,allocatable:: reclev(:)
21 !---------------------------------------------------------------------------
22 !--- local vars
23  character(3) cmon
24  character(16) reclevtyp_sht
25  character(32) ctldate,varname
26  character(35) sweep_blanks
27  real(8) radi
28  real lon_stt,lat_stt,lon_intl,lat_intl
29  integer i,n,j,krec,iret,io_unit,idrt,nmeta
30 !---------------------------------------------------------------------------
31 !---------------------------------------------------------------------------
32 !
33  call nemsio_init(iret=iret)
34  if(iret/=0) print *,'ERROR: nemsio_init '
35 !
36 !---------------------------------------------------------------------------
37 !*** read nemsio grd header info
38 !---------------------------------------------------------------------------
39 !--- open gfile for reading
40  call getarg(1,cin)
41  print *,'filename is cin=',trim(cin)
42  if(trim(cin)=='') then
43  print *,'usage: mknemsioctl input_nemsio_file_name'
44  stop
45  endif
46  call getarg(2,cindx)
47  if(trim(cindx)=='') then
48  cindx=trim(cin)//'.ctl'
49  endif
50 
51  call nemsio_open(gfile,trim(cin),'READ',iret=iret)
52  if(iret/=0) print *,'Error: open nemsio file,',trim(cin),' iret=',iret
53 
54  call nemsio_getfilehead(gfile,iret=iret,nrec=nrec,dimx=im,dimy=jm, &
55  dimz=lm,idate=idate,nfhour=nfhour,nmeta=nmeta, &
56  nsoil=nsoil,idrt=idrt,tlmeta=tlmeta)
57 !
58  fieldsize=im*jm
59  allocate(recname(nrec),reclevtyp(nrec),reclev(nrec))
60  allocate(lat1(fieldsize),lat(jm),lon1(fieldsize))
61  call nemsio_getfilehead(gfile,iret=iret,recname=recname, &
62  reclevtyp=reclevtyp,reclev=reclev,lat=lat1,lon=lon1)
63 !
64  call nemsio_close(gfile,iret=iret)
65 !
66  call nemsio_finalize()
67 !
68 !---------------------------------------------------------------------------
69 !****** write .ctl file
70 !---------------------------------------------------------------------------
71 !
72 !-- get date
73  call cmonth(idate(2),cmon)
74  write(ctldate,'(i2.2,a,i2.2,a3,i4.4)')idate(4),'Z',idate(3) &
75  ,cmon,idate(1)
76 !-- get Gaussian grid
77  if(idrt==4) then
78 !-- data has lat/lon info
79  if (nmeta>=8) then
80  lon_stt=lon1(1)
81  lon_intl=lon1(2)-lon1(1)
82  do j=1,jm
83  lat(j) = lat1(1+(j-1)*im)
84  enddo
85  else
86  call splat8(idrt,jm,slat)
87  radi=180.0/(4.*atan(1.0))
88  do j=1,jm
89  lat(j) = asin(slat(j)) * radi
90  enddo
91  lon_stt=0.
92  lon_intl=360./real(jm)
93  endif
94  elseif(idrt==0) then
95  if(nmeta>=8) then
96  lon_stt=lon1(1)
97  lon_intl=lon1(2)-lon1(1)
98  lat_stt=-abs(lat1(1))
99  lat_intl=abs(lat1(1+im)-lat1(1))
100 !! print*, "im=",im,' jm=',jm
101 !! print*, "lon1=",lon1
102 !! print*, "lat_intl=",lat_intl
103 !! print*, "lat1=",lat1
104  else
105  if (mod(jm,2) == 0) then
106  lat_intl=180.0/real(jm)
107  lat_stt=-90.0+0.5*lat_intl
108  else
109  lat_intl=180.0/real(jm-1)
110  lat_stt=-90.0
111  endif
112  lon_intl=360./real(im)
113  lon_stt=0.5*lon_intl
114  endif
115  endif
116 
117  io_unit=650
118  open(io_unit,file=trim(cindx),form='formatted')
119 !
120  write(io_unit,105)trim(cin)
121  write(io_unit,106)
122 
123  if (idrt == 4 )then
124  write(io_unit,107)
125  elseif(idrt==0) then
126  if(nmeta>=8) then
127  if (lat1(1)>lat1(im+1)) then
128  write(io_unit,107)
129  else
130  write(io_unit,1077)
131  endif
132  else
133  write(io_unit,1077)
134  endif
135  endif
136 
137  write(io_unit,108)tlmeta
138  write(io_unit,109)
139  write(io_unit,111)im,lon_stt,lon_intl
140  if (idrt==4) then
141  write(io_unit,112)jm
142  write(io_unit,113)(lat(i),i=jm,1,-1)
143  else
144  write(io_unit,120)jm,lat_stt,lat_intl
145  endif
146  write(io_unit,114)lm
147  if(nfhour/=0) then
148  write(io_unit,115)1,trim(ctldate),nfhour
149  else
150  write(io_unit,116)1,trim(ctldate)
151  endif
152 
153 !
154  105 FORMAT('dset ^',a)
155  106 FORMAT('undef 9.99E+20')
156  107 FORMAT('options big_endian sequential yrev')
157  1077 FORMAT('options big_endian sequential')
158  108 FORMAT('fileheader',i12.0)
159  109 FORMAT('title gfs nemsioi file')
160 
161  111 FORMAT('xdef ',i6,' linear ',f9.6,' ',f9.6)
162  112 FORMAT('ydef ',i6,' levels')
163  113 FORMAT(10f11.6)
164 
165  114 FORMAT('zdef ',i6,' linear 1 1')
166  115 FORMAT('tdef ',i6,' linear ',a12,' ',i6,'hr')
167  116 FORMAT('tdef ',i6,' linear ',a12,' 1yr')
168 !
169  120 FORMAT('ydef ',i6,' linear ',f13.6,' ',f9.6)
170 !
171  krec=0
172  do n=1,nrec
173  if(reclev(n)==1) then
174  krec=krec+1
175  endif
176  enddo
177 
178  WRITE(io_unit,'(A,I6)')'VARS ',krec
179 
180  n=1
181  do while (n<=nrec)
182  reclevtyp_sht=reclevtyp(n)
183  if(trim(reclevtyp_sht) == "convect-cld bot") then
184  reclevtyp_sht="cvb"
185  elseif (trim(reclevtyp_sht) == "convect-cld top") then
186  reclevtyp_sht="cvt"
187  elseif (trim(reclevtyp_sht) == "high cld bot") then
188  reclevtyp_sht="hcb"
189  elseif (trim(reclevtyp_sht) == "high cld top") then
190  reclevtyp_sht="hct"
191  elseif (trim(reclevtyp_sht) == "mid cld bot") then
192  reclevtyp_sht="mcb"
193  elseif (trim(reclevtyp_sht) == "mid cld top") then
194  reclevtyp_sht="mct"
195  elseif (trim(reclevtyp_sht) == "low cld bot") then
196  reclevtyp_sht="lcb"
197  elseif (trim(reclevtyp_sht) == "low cld top") then
198  reclevtyp_sht="lct"
199  elseif (trim(reclevtyp_sht) == "convect-cld laye") then
200  reclevtyp_sht="cvcl"
201  elseif (trim(reclevtyp_sht) == "bndary-layer cld") then
202  reclevtyp_sht="bdrlc"
203  elseif (trim(reclevtyp_sht) == "atmos col") then
204  reclevtyp_sht="acol"
205  elseif (trim(reclevtyp_sht) == "high cld lay") then
206  reclevtyp_sht="hcl"
207  elseif (trim(reclevtyp_sht) == "mid cld lay") then
208  reclevtyp_sht="mcl"
209  elseif (trim(reclevtyp_sht) == "low cld lay") then
210  reclevtyp_sht="lcl"
211  elseif (trim(reclevtyp_sht) == "2 m above gnd") then
212  reclevtyp_sht="2m"
213  elseif (trim(reclevtyp_sht) == "10 m above gnd") then
214  reclevtyp_sht="10m"
215  endif
216  varname=sweep_blanks(trim(recname(n))//trim(reclevtyp_sht))
217  if(trim(reclevtyp(n))=='mid layer') then
218  write(io_unit,'(a16,i3,a)')varname,lm,' 99 model layer'
219  n=n+lm
220  elseif(trim(reclevtyp(n))=='soil layer') then
221  write(io_unit,'(a16,i3,a)')varname,nsoil,' 99 soil layer'
222  n=n+nsoil
223  elseif(trim(reclevtyp(n))=='2 m above gnd') then
224  recname(n)=trim(recname(n))//'2m'
225  write(io_unit,'(a16,a7,a)')varname,' 0 99 ',trim(reclevtyp(n))
226  n=n+1
227  elseif(trim(reclevtyp(n))=='10 m above gnd') then
228  recname(n)=trim(recname(n))//'10m'
229  write(io_unit,'(a16,a7,a)')varname,' 0 99 ',trim(reclevtyp(n))
230  n=n+1
231  else
232  write(io_unit,'(a16,a7,a)')varname,' 0 99 ',trim(reclevtyp(n))
233  n=n+1
234  endif
235  enddo
236 
237  write(io_unit,'(A8)')'endvars'
238  close(io_unit)
239 
240 !---------------------------------------------------------------------------
241 !****** clean up
242 !---------------------------------------------------------------------------
243  deallocate(recname,reclevtyp,reclev,lat,lat1,lon1)
244 !---------------------------------------------------------------------------
245 !
246 ! - - - - -- - -- - -- - -- - - -- - -- -- - -- - -- - - -- - - - -- - --
247  stop
248 
249  end program
250 ! - - - - -- - -- - -- - -- - - -- - -- -- - -- - -- - - -- - - - -- - --
251  SUBROUTINE cmonth(IMON,CMON)
252 !
253 !-----------------------------------------------------------------------
254 !*** Convert month
255 !-----------------------------------------------------------------------
256 !
257  INTEGER,INTENT(IN) :: IMON
258  CHARACTER(LEN=3) :: CMON
259 !
260 !-----------------------------------------------------------------------
261 !
262  SELECT CASE (imon)
263  CASE(1)
264  cmon='Jan'
265  CASE(2)
266  cmon='Feb'
267  CASE(3)
268  cmon='Mar'
269  CASE(4)
270  cmon='Apr'
271  CASE(5)
272  cmon='May'
273  CASE(6)
274  cmon='Jun'
275  CASE(7)
276  cmon='Jul'
277  CASE(8)
278  cmon='Aug'
279  CASE(9)
280  cmon='Sep'
281  CASE(10)
282  cmon='Oct'
283  CASE(11)
284  cmon='Nov'
285  CASE(12)
286  cmon='Dec'
287  END SELECT
288 !
289 !-----------------------------------------------------------------------
290 !
291  END SUBROUTINE cmonth
292 
293 !
294  sUBROUTINE splat4(IDRT,JMAX,ASLAT)
295 !
296  implicit none
297  integer,intent(in) :: idrt,jmax
298  real(4),intent(out) :: ASLAT(JMAX)
299  INTEGER,PARAMETER:: KD=selected_real_kind(15,45)
300  REAL(KIND=kd):: pk(jmax/2),pkm1(jmax/2),pkm2(jmax/2)
301  REAL(KIND=kd):: aslatd(jmax/2),sp,spmax,eps=10.d0*epsilon(sp)
302  integer,PARAMETER:: JZ=50
303  REAL(8) BZ(JZ)
304  DATA bz / 2.4048255577d0, 5.5200781103d0, &
305  8.6537279129d0, 11.7915344391d0, 14.9309177086d0, 18.0710639679d0, &
306  21.2116366299d0, 24.3524715308d0, 27.4934791320d0, 30.6346064684d0, &
307  33.7758202136d0, 36.9170983537d0, 40.0584257646d0, 43.1997917132d0, &
308  46.3411883717d0, 49.4826098974d0, 52.6240518411d0, 55.7655107550d0, &
309  58.9069839261d0, 62.0484691902d0, 65.1899648002d0, 68.3314693299d0, &
310  71.4729816036d0, 74.6145006437d0, 77.7560256304d0, 80.8975558711d0, &
311  84.0390907769d0, 87.1806298436d0, 90.3221726372d0, 93.4637187819d0, &
312  96.6052679510d0, 99.7468198587d0, 102.888374254d0, 106.029930916d0, &
313  109.171489649d0, 112.313050280d0, 115.454612653d0, 118.596176630d0, &
314  121.737742088d0, 124.879308913d0, 128.020877005d0, 131.162446275d0, &
315  134.304016638d0, 137.445588020d0, 140.587160352d0, 143.728733573d0, &
316  146.870307625d0, 150.011882457d0, 153.153458019d0, 156.295034268d0 /
317  REAL(8):: DLT,D1=1.d0
318  INTEGER(4):: JHE,JHO,J0=0
319  real(8),PARAMETER :: PI=3.14159265358979d0,c=(1.d0-(2.d0/pi)**2)*0.25d0
320  real(8) r
321  integer jh,js,n,j
322 !C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
323 !C GAUSSIAN LATITUDES
324  IF(idrt.EQ.4) THEN
325  jh=jmax/2
326  jhe=(jmax+1)/2
327  r=1.d0/sqrt((jmax+0.5d0)**2+c)
328  DO j=1,min(jh,jz)
329  aslatd(j)=cos(bz(j)*r)
330  ENDDO
331  DO j=jz+1,jh
332  aslatd(j)=cos((bz(jz)+(j-jz)*pi)*r)
333  ENDDO
334  spmax=1.d0
335  DO WHILE(spmax.GT.eps)
336  spmax=0.d0
337  DO j=1,jh
338  pkm1(j)=1.d0
339  pk(j)=aslatd(j)
340  ENDDO
341  DO n=2,jmax
342  DO j=1,jh
343  pkm2(j)=pkm1(j)
344  pkm1(j)=pk(j)
345  pk(j)=((2*n-1)*aslatd(j)*pkm1(j)-(n-1)*pkm2(j))/n
346  ENDDO
347  ENDDO
348  DO j=1,jh
349  sp=pk(j)*(1.d0-aslatd(j)**2)/(jmax*(pkm1(j)-aslatd(j)*pk(j)))
350  aslatd(j)=aslatd(j)-sp
351  spmax=max(spmax,abs(sp))
352  ENDDO
353  ENDDO
354 !CDIR$ IVDEP
355  DO j=1,jh
356  aslat(j)=aslatd(j)
357  aslat(jmax+1-j)=-aslat(j)
358  ENDDO
359  IF(jhe.GT.jh) THEN
360  aslat(jhe)=0.d0
361  ENDIF
362 !C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
363 !C EQUALLY-SPACED LATITUDES INCLUDING POLES
364  ELSEIF(idrt.EQ.0) THEN
365  jh=jmax/2
366  jhe=(jmax+1)/2
367  jho=jhe-1
368  dlt=pi/(jmax-1)
369  aslat(1)=1.d0
370  DO j=2,jh
371  aslat(j)=cos((j-1)*dlt)
372  ENDDO
373 !CDIR$ IVDEP
374  DO j=1,jh
375  aslat(jmax+1-j)=-aslat(j)
376  ENDDO
377  IF(jhe.GT.jh) THEN
378  aslat(jhe)=0.d0
379  ENDIF
380 !C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
381 !C EQUALLY-SPACED LATITUDES EXCLUDING POLES
382  ELSEIF(idrt.EQ.256) THEN
383  jh=jmax/2
384  jhe=(jmax+1)/2
385  jho=jhe
386  dlt=pi/jmax
387  aslat(1)=1.d0
388  DO j=1,jh
389  aslat(j)=cos((j-0.5)*dlt)
390  ENDDO
391 !CDIR$ IVDEP
392  DO j=1,jh
393  aslat(jmax+1-j)=-aslat(j)
394  ENDDO
395  IF(jhe.GT.jh) THEN
396  aslat(jhe)=0.d0
397  ENDIF
398  ENDIF
399 !C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
400  end subroutine splat4
401 !
402 !----------------------------------------------------------------------
403  SUBROUTINE splat8(IDRT,JMAX,ASLAT)
404 !$$$
405  implicit none
406  integer,intent(in) :: idrt,jmax
407  real(8),intent(out) :: ASLAT(JMAX)
408  INTEGER,PARAMETER:: KD=selected_real_kind(15,45)
409  REAL(KIND=kd):: pk(jmax/2),pkm1(jmax/2),pkm2(jmax/2)
410  REAL(KIND=kd):: aslatd(jmax/2),sp,spmax,eps=10.d0*epsilon(sp)
411  integer,PARAMETER:: JZ=50
412  REAL(8) BZ(JZ)
413  DATA bz / 2.4048255577d0, 5.5200781103d0, &
414  8.6537279129d0, 11.7915344391d0, 14.9309177086d0, 18.0710639679d0, &
415  21.2116366299d0, 24.3524715308d0, 27.4934791320d0, 30.6346064684d0, &
416  33.7758202136d0, 36.9170983537d0, 40.0584257646d0, 43.1997917132d0, &
417  46.3411883717d0, 49.4826098974d0, 52.6240518411d0, 55.7655107550d0, &
418  58.9069839261d0, 62.0484691902d0, 65.1899648002d0, 68.3314693299d0, &
419  71.4729816036d0, 74.6145006437d0, 77.7560256304d0, 80.8975558711d0, &
420  84.0390907769d0, 87.1806298436d0, 90.3221726372d0, 93.4637187819d0, &
421  96.6052679510d0, 99.7468198587d0, 102.888374254d0, 106.029930916d0, &
422  109.171489649d0, 112.313050280d0, 115.454612653d0, 118.596176630d0, &
423  121.737742088d0, 124.879308913d0, 128.020877005d0, 131.162446275d0, &
424  134.304016638d0, 137.445588020d0, 140.587160352d0, 143.728733573d0, &
425  146.870307625d0, 150.011882457d0, 153.153458019d0, 156.295034268d0 /
426  REAL(8):: DLT,D1=1.d0
427  INTEGER(4):: JHE,JHO,J0=0
428  real(8),PARAMETER :: PI=3.14159265358979d0,c=(1.d0-(2.d0/pi)**2)*0.25d0
429  real(8) r
430  integer jh,js,n,j
431 !C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
432 !C GAUSSIAN LATITUDES
433  IF(idrt.EQ.4) THEN
434  jh=jmax/2
435  jhe=(jmax+1)/2
436  r=1.d0/sqrt((jmax+0.5d0)**2+c)
437  DO j=1,min(jh,jz)
438  aslatd(j)=cos(bz(j)*r)
439  ENDDO
440  DO j=jz+1,jh
441  aslatd(j)=cos((bz(jz)+(j-jz)*pi)*r)
442  ENDDO
443  spmax=1.d0
444  DO WHILE(spmax.GT.eps)
445  spmax=0.d0
446  DO j=1,jh
447  pkm1(j)=1.d0
448  pk(j)=aslatd(j)
449  ENDDO
450  DO n=2,jmax
451  DO j=1,jh
452  pkm2(j)=pkm1(j)
453  pkm1(j)=pk(j)
454  pk(j)=((2*n-1)*aslatd(j)*pkm1(j)-(n-1)*pkm2(j))/n
455  ENDDO
456  ENDDO
457  DO j=1,jh
458  sp=pk(j)*(1.d0-aslatd(j)**2)/(jmax*(pkm1(j)-aslatd(j)*pk(j)))
459  aslatd(j)=aslatd(j)-sp
460  spmax=max(spmax,abs(sp))
461  ENDDO
462  ENDDO
463 !CDIR$ IVDEP
464  DO j=1,jh
465  aslat(j)=aslatd(j)
466  aslat(jmax+1-j)=-aslat(j)
467  ENDDO
468  IF(jhe.GT.jh) THEN
469  aslat(jhe)=0.d0
470  ENDIF
471 !C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
472 !C EQUALLY-SPACED LATITUDES INCLUDING POLES
473  ELSEIF(idrt.EQ.0) THEN
474  jh=jmax/2
475  jhe=(jmax+1)/2
476  jho=jhe-1
477  dlt=pi/(jmax-1)
478  aslat(1)=1.d0
479  DO j=2,jh
480  aslat(j)=cos((j-1)*dlt)
481  ENDDO
482 !CDIR$ IVDEP
483  DO j=1,jh
484  aslat(jmax+1-j)=-aslat(j)
485  ENDDO
486  IF(jhe.GT.jh) THEN
487  aslat(jhe)=0.d0
488  ENDIF
489 !C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
490 !C EQUALLY-SPACED LATITUDES EXCLUDING POLES
491  ELSEIF(idrt.EQ.256) THEN
492  jh=jmax/2
493  jhe=(jmax+1)/2
494  jho=jhe
495  dlt=pi/jmax
496  aslat(1)=1.d0
497  DO j=1,jh
498  aslat(j)=cos((j-0.5d0)*dlt)
499  ENDDO
500 !DIR$ IVDEP
501  DO j=1,jh
502  aslat(jmax+1-j)=-aslat(j)
503  ENDDO
504  IF(jhe.GT.jh) THEN
505  aslat(jhe)=0.d0
506  ENDIF
507  ENDIF
508 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
509  end subroutine splat8
510 !-----------------------------------------------------------------------
511 
512  character(35) function sweep_blanks(in_str)
513 
514 !
515  implicit none
516 !
517  character(*), intent(in) :: in_str
518  character(35) :: out_str
519  character :: ch
520  integer :: j
521 
522  out_str = " "
523  do j=1, len_trim(in_str)
524  ! get j-th char
525  ch = in_str(j:j)
526  if (ch .eq. "-") then
527  out_str = trim(out_str) // "_"
528  else if (ch .ne. " ") then
529  out_str = trim(out_str) // ch
530  endif
531  sweep_blanks = out_str
532  end do
533  end function sweep_blanks
534