source: trunk/MESOSCALE/LMD_MM_MARS/SRC/POSTPROC/w2g.f @ 134

Last change on this file since 134 was 19, checked in by aslmd, 14 years ago

spiga:mineur

File size: 41.1 KB
Line 
1!=================================
2!! Aymeric Spiga - 10/2009
3!! NB: this is a slightly modified version of WRFnc2ctl.f called w2g
4!! NB: it is able to read Martian files
5!=================================
6
7
8!! This program read a WRF netcdf file and create the necessary .ctl
9!! files to display WRF netcdf directly with GrADS.
10!!
11!! Note, due to staggering, you get 4 .ctl files, one for mass points, one
12!! for U points, one for V points and one for W points.
13!! They need to be open all at once in GrADS.
14!!
15!! When using this method to display WRF in GrADS, one cannot inpertolate
16!! to pressure/height levels, and no extra diagnostics are available
17!
18!=================================Make Executable============================
19!  Make executable:
20!    DEC Alpha
21!      f90 WRFnc2ctl.f -L/usr/local/netcdf/lib -lnetcdf -lm  \
22!      -I/usr/local/netcdf/include  -free  -o WRFnc2ctl
23!
24!   linux flags
25!      pgf90 WRFnc2ctl.f -L/usr/local/netcdf/lib -lnetcdf -lm  \
26!      -I/usr/local/netcdf/include  -Mfree  -o WRFnc2ctl
27!
28!   Sun flags
29!      f90 WRFnc2ctl.f -L/usr/local/netcdf/lib -lnetcdf -lm  \
30!      -I/usr/local/netcdf/include  -free  -o WRFnc2ctl
31!
32!   SGI flags
33!      f90 WRFnc2ctl.f -L/usr/local/netcdf/lib -lnetcdf -lm  \
34!      -I/usr/local/netcdf/include  -freeform  -o WRFnc2ctl
35!
36!   IBM flags
37!      xlf WRFnc2ctl.f -L/usr/local/lib32/r4i4 -lnetcdf -lm  \
38!      -I/usr/local/netcdf/include  -qfree=f90  -o WRFnc2ctl
39!
40!   Mac flags (with xlf compiler)
41!      xlf WRFnc2ctl.f -L/usr/local/netcdf-xlf/lib -lnetcdf -lm  \
42!      -I/usr/local/netcdf-xlf/include  -qfree=f90  -o WRFnc2ctl
43!
44!   If you extra compile flags for other computers - please send along
45!
46!=================================Run Program================================
47!  Run program:
48!
49!          WRFnc2ctl   -i wrf_netcdf_input_file   -o ctl_output_name
50!
51!=================================Options====================================
52!
53! -help     : Print help information                           
54! -h        : Print help information                           
55! -i        : WRF input file (netcdf format)                   
56! -o        : ctl output name                                 
57! -debug    : Print some debug information
58!
59!============================================================================
60!
61!  December 2006
62!  Cindy Bruyere
63!
64     program WRFnc2ctl
65
66     implicit none
67     include 'netcdf.inc'
68
69  character (len=200)                                 :: input_file                  ! netcdf input file
70  logical                                             :: have_vert_stag=.FALSE.      ! also need W.ctl file
71  integer                                             :: cdfid                       ! input file unit number
72  real                                                :: rcode                       ! return status code for netcdf
73  integer                                             :: ndims, nvars, natts         ! number of dims, var and att in file
74  integer                                             :: unlimdimid                  ! unlimited dims ID - not used
75
76  character (len=80)                                  :: case                        ! ctl file name
77  character (len=80)                                  :: ctl_U, ctl_V, ctl_W, ctl_M  ! ctl file name
78  character (len=200), allocatable, dimension(:)      :: M_vars                      ! variables to ctl
79  character (len=200),              dimension(10)     :: U_vars, V_vars, W_vars      ! variables to ctl
80  integer                                             :: i_M, i_U, i_V, i_W          ! number of t/u/v/w variables
81  character (len=200)                                 :: var_string                  ! tmp var string
82  character (len=40)                                  :: tdef                        ! ctl tdef
83  integer                                             :: ntimes                      ! number of times in input file
84  integer                                             :: end_ctl_files               ! sometimes we will have 3, sometimes 4 .ctl files
85
86  character (len=31), allocatable, dimension(:)       :: dnam                        ! name of netcdf DIMS
87  integer,            allocatable, dimension(:)       :: dval                        ! value of netcdf DIMS
88  character (len=40)                                  :: title                       ! global att from netcdf file
89  character (len=1)                                   :: gridtype                    ! global att from netcdf file
90  real                                                :: cen_lon, cen_lat            ! global att from netcdf file
91  real                                                :: stand_lon                   ! global att from netcdf file
92  real                                                :: truelat1, truelat2          ! global att from netcdf file
93  real                                                :: dx, dy                      ! global att from netcdf file
94  integer                                             :: map_proj                    ! global att from netcdf file
95
96  character (len=80)                                  :: varnam                      ! var from netcdf file
97  integer                                             :: idvar                       ! varnam id number     
98  integer                                             :: ivtype                      ! varnam type         
99  integer                                             :: idm                         ! varnam - number of dims
100  integer                                             :: natt                        ! varnam attributes
101  character (len=80)                                  :: description                 ! varnam description
102  character (len=80)                                  :: units                       ! varnam units
103  character (len=80)                                  :: varnam_small                ! lowecase of varnam - for ctl file
104  character (len=80)                                  :: varout                      ! varnam + varnam_small - for ctl file
105
106  integer                                             :: len_unt, len_title          ! length of character strings 
107  integer                                             :: len_string, len_des         ! length of character strings 
108  integer                                             :: len_var                     ! length of character strings 
109
110  integer                                             :: dims(4), ishape(10)         ! netcdf array dims and shapes
111  character,          allocatable, dimension(:,:,:,:) :: times                       ! Times array from netcdf file
112  real,               allocatable, dimension(:,:,:,:) :: znu, znw                    ! arrays from netcdf file
113  real,               allocatable, dimension(:,:,:,:) :: xlat_M, xlon_M              ! u-staggered lat/lon     
114  real,               allocatable, dimension(:,:,:,:) :: xlat_U, xlon_U              ! u-staggered lat/lon     
115  real,               allocatable, dimension(:,:,:,:) :: xlat_V, xlon_V              ! v-staggered lat/lon     
116  real,                            dimension(4)       :: xlat_a, xlon_a              ! lat/lon corners
117  real                                                :: lat_ll_u, lat_ll_v          ! lat/lon corners
118  real                                                :: lat_ur_u, lat_ur_v          ! lat/lon corners
119  real                                                :: lon_ll_u, lon_ll_v          ! lat/lon corners
120  real                                                :: lon_ur_u, lon_ur_v          ! lat/lon corners
121  integer                                             :: iweg, isng, ibtg            ! staggered dims in input_file
122  real,                            dimension(16)      :: lats16, lons16              ! lat/lon corners
123 
124  real                                                :: xi, xj                      ! PS defs
125  real                                                :: r, re                       ! PS defs
126  real                                                :: ref_lat, ref_lon            ! PS defs
127  real                                                :: rlat, rlong                 ! PS defs
128  real                                                :: earthr, radpd               ! PS defs
129  real                                                :: dx_ps                       ! PS defs
130  integer                                             :: ipole                       ! PS defs
131
132  integer                                             :: ilon                        ! flag indicating dateline in domain
133  real                                                :: abslatmin, abslatmax        ! min/max lat - for Lambert def's
134  real                                                :: abslonmin, abslonmax        ! min/max lon - for Lambert def's
135  integer                                             :: ipoints, jpoints            ! points used for Lambert def's     
136  real                                                :: dxll, slack                 ! Lambert area setup (xdef/ydef)   
137
138  integer                                             :: itmp, i, k                  ! tmp values                   
139  real                                                :: temp                        ! tmp values                   
140  real,               allocatable, dimension(:,:,:,:) :: tmp_var                     ! tmp values           
141
142  logical                                             :: mercator_defs               ! for better mercator projections
143  logical                                             :: debug                       ! for debug
144
145  integer                                             :: test1, test2  !!ajout
146
147
148   write(*,*) "  "
149   write(*,*) "================================================ "
150   write(*,*) "WRFnc2ctl for ARW WRF files      "
151   write(*,*) " Only works correctly for gradsnc from GrADS 1.9 "
152
153  tdef = '1 linear 00z01jan2000  1hr'
154
155! GET INPUT FILE AND OUTPUT CASE NAME FOR COMMAND LINE
156  if (debug) write(*,*) "READ arguments from command line"
157  call read_args(input_file,case,mercator_defs,debug)
158
159
160!! OPEN THE INPUT FILE
161   if (debug) write(*,*) "OPENING netcdf input file"           
162   rcode = nf_open (input_file, NF_NOWRITE, cdfid)
163   call handle_err ("Opening input netCDF file", rcode)
164   rcode = nf_get_att_text(cdfid, nf_global, 'GRIDTYPE', gridtype)
165   if ( trim(gridtype) .ne. 'C' ) then
166      write(*,*) 'Can only deal with ARW data for now'
167      STOP
168   endif
169
170
171!! GET BASIC INFORMATION ABOUT THE FILE
172   if (debug) write(*,*) "READ global attributes from netcdf file"
173   stand_lon = -99999.99   !!! backward compatibility
174   rcode = nf_get_att_text(cdfid, nf_global, 'TITLE',     title)
175   rcode = nf_get_att_real(cdfid, nf_global, 'DX',        dx)
176   rcode = nf_get_att_real(cdfid, nf_global, 'DY',        dy)
177   rcode = nf_get_att_real(cdfid, nf_global, 'CEN_LAT',   cen_lat)
178   rcode = nf_get_att_real(cdfid, nf_global, 'CEN_LON',   cen_lon)
179   rcode = nf_get_att_real(cdfid, nf_global, 'STAND_LON', stand_lon)
180   rcode = nf_get_att_real(cdfid, nf_global, 'TRUELAT1',  truelat1)
181   rcode = nf_get_att_real(cdfid, nf_global, 'TRUELAT2',  truelat2)
182   rcode = nf_get_att_int (cdfid, nf_global, 'MAP_PROJ',  map_proj)
183   rcode = nf_get_att_int (cdfid, nf_global, 'WEST-EAST_GRID_DIMENSION',  test1)    !!ajout
184   rcode = nf_get_att_int (cdfid, nf_global, 'SOUTH-NORTH_GRID_DIMENSION',  test2)  !!ajout
185   if (stand_lon == -99999.99 ) stand_lon = cen_lat    !!! dealing with an old WRF file
186
187   if (debug) write(*,*) "READ dims from netcdf file"
188   rcode = nf_inq(cdfid, nDims, nVars, nAtts, unlimDimID)
189
190   print *, 'BEWARE if 64 bits NETCDF...'
191   print *, cdfid, nDims, nVars, nAtts, unlimDimID
192
193
194   allocate (M_vars(nVars))
195   allocate (dval(nDims))
196   allocate (dnam(nDims))
197   iweg = 1 
198   isng = 1 
199   ibtg = 1 
200   do i = 1,nDims
201     rcode = nf_inq_dim(cdfid, i, dnam(i), dval(i))
202     if ( dnam(i)(1:9)  == 'west_east'          ) iweg = max(iweg,dval(i))
203     if ( dnam(i)(1:11) == 'south_north'        ) isng = max(isng,dval(i))
204
205     if ( dnam(i)(1:10) == 'bottom_top'         ) ibtg = max(ibtg,dval(i))
206     if ( trim(dnam(i)) == 'num_metgrid_levels' ) ibtg = max(ibtg,dval(i))
207     if ( dnam(i)(1:4)  == 'land'               ) ibtg = max(ibtg,dval(i))
208     if ( dnam(i)(1:4)  == 'soil'               ) ibtg = max(ibtg,dval(i))
209     if ( dnam(i)(1:5)  == 'month'              ) ibtg = max(ibtg,dval(i))
210     if ( dnam(i)(1:5)  == 'z-dim'              ) ibtg = max(ibtg,dval(i))
211
212     if ( dnam(i)(1:15) == 'bottom_top_stag'    ) have_vert_stag = .TRUE.   
213   enddo
214
215!! OPEN THE OUTPUT FILES
216   write(*,*) 
217   write(*,*) "CREATING .ctl files for you case: ", case
218   write(*,*) 
219!   ctl_M = trim(case)//"_M.ctl"
220   ctl_M = trim(case)//".ctl"
221   ctl_U = trim(case)//"_U.ctl"
222   ctl_V = trim(case)//"_V.ctl"
223   ctl_W = trim(case)//"_W.ctl"
224   open (10, file=ctl_M)
225!   open (11, file=ctl_U)
226!   open (12, file=ctl_V)
227   if (have_vert_stag) open (13, file=ctl_W)
228
229!!-------------------------------------------------
230!! AS: ici fix pour les fichiers sortie de api qui sont deja de-stagg
231!horizontalement
232if (iweg .ne. test1) then
233    iweg=iweg+1
234    test1=-1
235    PRINT *, 'IF YOU HAVE A BOXED FILE w STAGG VARIABLES it is a PROBLEM'
236endif
237if (isng .ne. test2) then
238    isng=isng+1
239    test2=-1
240endif
241if (test1 .ne. -1) open (11, file=ctl_U) 
242if (test2 .ne. -1) open (12, file=ctl_V) 
243!!-------------------------------------------------
244
245  lat_ll_u = -999.99
246  lat_ll_v = -999.99
247  lat_ur_u = -999.99
248  lat_ur_v = -999.99
249  lon_ll_u = -999.99
250  lon_ll_v = -999.99
251  lon_ur_u = -999.99
252  lon_ur_v = -999.99
253
254!===========================================================================================
255
256      i_U = 0
257      i_V = 0
258      i_W = 0
259      i_M = 0
260
261      DO idvar = 1,nVars        !! LOOP THROUGH ALL VARIABLES IN FILE
262        rcode = nf_inq_var(cdfid, idvar, varnam, ivtype, idm, ishape, natt)
263        call handle_err ("ERROR reading variable", rcode)
264        if (debug) write(*,*) 
265        if (debug) write(*,*) "Variable:",idvar,"out of",nVars
266        if (debug) write(*,*) "DEALING with variable: ", trim(varnam)
267        if (.not. debug ) write(*,'("   Variable ",i2," :  ",A10,$)') idvar,varnam
268        if (.not. debug .and. idm .gt. 2 ) write(*,*) " (*)"
269        if (.not. debug .and. idm .le. 2 ) write(*,*) "    "
270
271!! GET THE DIMS FOR INPUT AND OUTPUT FROM THE SHAPE
272        dims  = 1
273        do i = 1,idm
274          dims(i)  = dval(ishape(i))
275        enddo
276        if (debug) write(*,*) " DIMS: ", dims
277
278
279!! GET THE UNITS AND DESCRIPTION OF THE FIELD
280        units = "                                                                                "
281        description = "                                                                                "
282        rcode = NF_GET_ATT_TEXT(cdfid, idvar, "units", units )
283        rcode = NF_GET_ATT_TEXT(cdfid, idvar, "description", description )
284        len_var = len_trim(varnam)
285        len_des = len_trim(description)
286        len_unt = len_trim(units)
287        if (debug) write(*,*) " DESCRIPTION: ", description(1:len_des)
288        if (debug) write(*,*) " UNITS: ", units(1:len_unt)
289
290!! GET THE LOWER CASE OF EACH varnam
291        varnam_small = " " 
292        do i = 1,len_var
293          if (varnam(i:i) .ge. 'A' .and. varnam(i:i) .le. 'Z') then
294            itmp = ichar(varnam(i:i)) + 32
295            varnam_small(i:i) = achar(itmp)
296          else
297            varnam_small(i:i) = varnam(i:i)
298          endif
299        enddo
300        write(varout,*) varnam(1:len_var)//"=>"//varnam_small(1:len_var)
301        len_var = max(22,len_trim(varout))
302
303!! BUILD THE CTL var LIST
304        IF     ( idm == 4 ) THEN
305          write(var_string,'(A,i4,"  t,z,y,x  ",A,"  (",A,")")') varout(1:len_var), dims(3), &
306description(1:len_des), units(1:len_unt)
307        ELSEIF ( idm == 3 ) THEN
308          write(var_string,'(A,i4,"  t,y,x    ",A,"  (",A,")")') varout(1:len_var), 0, &
309description(1:len_des), units(1:len_unt)
310        ENDIF
311        IF     ( idm .gt. 2 ) THEN
312          len_string = len_trim(var_string)
313          if     ( dims(1) == iweg ) then
314               i_U = i_U + 1
315               U_vars(i_U) = var_string(1:len_string)
316          elseif ( dims(2) == isng ) then
317               i_V = i_V + 1
318               V_vars(i_V) = var_string(1:len_string)
319          elseif ( dims(3) == ibtg .and. have_vert_stag) then
320               i_W = i_W + 1
321               W_vars(i_W) = var_string(1:len_string)
322          else
323               i_M = i_M + 1
324               M_vars(i_M) = var_string(1:len_string)
325          endif
326          if (debug) write(*,*) " CTL: ",var_string(1:len_string)
327        ENDIF
328
329
330!!!!
331!!!! WE NEED KEEP SOME EXTRA INFO BEFORE READING THE NEXT VARIABLE
332!!!!
333
334
335!! GET THE NUMBER OF TIMES IN FILE AND BUILD tdef
336        IF (varnam == 'Times' ) THEN
337           allocate (times(dims(1), dims(2), dims(3), dims(4)))
338           rcode = nf_get_var_text(cdfid, idvar, times)
339           ntimes = dims(2)
340           call time_calc( times, ntimes, tdef, debug )
341        ENDIF
342
343!! GET VERTICAL INFORMATION FOR zdef, AND xlon/xlat FOR pdef
344        IF (varnam == 'ZNU' ) THEN
345           allocate (znu(dims(1), dims(2), dims(3), dims(4)))
346           rcode = nf_get_var_real(cdfid, idvar, znu)
347        ENDIF
348        IF (varnam == 'ZNW' ) THEN
349           allocate (znw(dims(1), dims(2), dims(3), dims(4)))
350           rcode = nf_get_var_real(cdfid, idvar, znw)
351        ENDIF
352        IF (varnam == 'XLAT'  ) THEN
353           allocate (xlat_M(dims(1), dims(2), dims(3), dims(4)))
354           rcode = nf_get_var_real(cdfid, idvar, xlat_M)
355        ENDIF
356        IF (varnam == 'XLAT_M'  ) THEN
357           allocate (xlat_M(dims(1), dims(2), dims(3), dims(4)))
358           rcode = nf_get_var_real(cdfid, idvar, xlat_M)
359        ENDIF
360        IF (varnam == 'XLAT_U'  ) THEN
361           allocate (xlat_U(dims(1), dims(2), dims(3), dims(4)))
362           rcode = nf_get_var_real(cdfid, idvar, xlat_U)
363        ENDIF
364        IF (varnam == 'XLAT_V'  ) THEN
365           allocate (xlat_V(dims(1), dims(2), dims(3), dims(4)))
366           rcode = nf_get_var_real(cdfid, idvar, xlat_V)
367        ENDIF
368        IF (varnam == 'XLONG' ) THEN
369           allocate (xlon_M(dims(1), dims(2), dims(3), dims(4)))
370           rcode = nf_get_var_real(cdfid, idvar, xlon_M)
371        ENDIF
372        IF (varnam == 'XLONG_M' ) THEN
373           allocate (xlon_M(dims(1), dims(2), dims(3), dims(4)))
374           rcode = nf_get_var_real(cdfid, idvar, xlon_M)
375        ENDIF
376        IF (varnam == 'XLONG_U'  ) THEN
377           allocate (xlon_U(dims(1), dims(2), dims(3), dims(4)))
378           rcode = nf_get_var_real(cdfid, idvar, xlon_U)
379        ENDIF
380        IF (varnam == 'XLONG_V' ) THEN
381           allocate (xlon_V(dims(1), dims(2), dims(3), dims(4)))
382           rcode = nf_get_var_real(cdfid, idvar, xlon_V)
383        ENDIF
384
385!! GET SOME CORNER INFORMATION
386        IF (varnam == 'LAT_LL_U' ) THEN
387           allocate (tmp_var(dims(1), dims(2), dims(3), dims(4)))
388           rcode = nf_get_var_real(cdfid, idvar, tmp_var)
389           lat_ll_u = tmp_var(1,1,1,1)
390           deallocate (tmp_var)
391        ENDIF
392        IF (varnam == 'LAT_UR_U' ) THEN
393           allocate (tmp_var(dims(1), dims(2), dims(3), dims(4)))
394           rcode = nf_get_var_real(cdfid, idvar, tmp_var)
395           lat_ur_u = tmp_var(1,1,1,1)
396           deallocate (tmp_var)
397        ENDIF
398        IF (varnam == 'LAT_LL_V' ) THEN
399           allocate (tmp_var(dims(1), dims(2), dims(3), dims(4)))
400           rcode = nf_get_var_real(cdfid, idvar, tmp_var)
401           lat_ll_v = tmp_var(1,1,1,1)
402           deallocate (tmp_var)
403        ENDIF
404        IF (varnam == 'LAT_UR_V' ) THEN
405           allocate (tmp_var(dims(1), dims(2), dims(3), dims(4)))
406           rcode = nf_get_var_real(cdfid, idvar, tmp_var)
407           lat_ur_v = tmp_var(1,1,1,1)
408           deallocate (tmp_var)
409        ENDIF
410
411        IF (varnam == 'LON_LL_U' ) THEN
412           allocate (tmp_var(dims(1), dims(2), dims(3), dims(4)))
413           rcode = nf_get_var_real(cdfid, idvar, tmp_var)
414           lon_ll_u = tmp_var(1,1,1,1)
415           deallocate (tmp_var)
416        ENDIF
417        IF (varnam == 'LON_UR_U' ) THEN
418           allocate (tmp_var(dims(1), dims(2), dims(3), dims(4)))
419           rcode = nf_get_var_real(cdfid, idvar, tmp_var)
420           lon_ur_u = tmp_var(1,1,1,1)
421           deallocate (tmp_var)
422        ENDIF
423        IF (varnam == 'LON_LL_V' ) THEN
424           allocate (tmp_var(dims(1), dims(2), dims(3), dims(4)))
425           rcode = nf_get_var_real(cdfid, idvar, tmp_var)
426           lon_ll_v = tmp_var(1,1,1,1)
427           deallocate (tmp_var)
428        ENDIF
429        IF (varnam == 'LON_UR_V' ) THEN
430           allocate (tmp_var(dims(1), dims(2), dims(3), dims(4)))
431           rcode = nf_get_var_real(cdfid, idvar, tmp_var)
432           lon_ur_v = tmp_var(1,1,1,1)
433           deallocate (tmp_var)
434        ENDIF
435
436
437     
438      ENDDO       !! END OF VAR LOOP
439      write(*,*) " "
440      write(*,*) "      (*) = variables that will be available in CTL files"
441      write(*,*) " "
442      write(*,*) "DONE reading variables"
443      if (debug) write(*,*) " "
444!===========================================================================================
445 
446!! CREATE THE CTL FILES
447
448   
449   if (debug) write(*,*) "START writing to .ctl files"
450   len_title = len_trim(title)
451   end_ctl_files = 12
452   if (have_vert_stag) end_ctl_files = 13
453   do i = 10,end_ctl_files
454      write(i,'("dset ^",A)') trim(input_file)
455      write(i,'("dtype netcdf")') 
456      write(i, '("undef 1.e37")')
457      write(i,'("title ",A)') title(1:len_title)
458   enddo
459   
460   IF (map_proj == 1) THEN                                     !! Lambert Projection
461     if (debug) write(*,*) "We have Lambert Projection data"
462      !PDEF for M
463      write(10,'("pdef ",i4," ",i3," lcc ",f7.3," ",f8.3," ",     &
464&          f8.3," ",f8.3," ",f9.5," ",f9.5," ",f10.5," ",         &
465&          f10.3," ",f10.3)')                                     &
466           iweg-1,isng-1,cen_lat,cen_lon,iweg/2.,isng/2.,         &
467           truelat1,truelat2,stand_lon,dx,dy
468      !PDEF for U
469      if (test1 .ne. -1) write(11,'("pdef ",i4," ",i3," lcc ",f7.3," ",f8.3," ",     &
470&          f8.3," ",f8.3," ",f9.5," ",f9.5," ",f10.5," ",         &
471&          f10.3," ",f10.3)')                                     &
472           iweg,isng-1,cen_lat,cen_lon,(iweg+1.)/2.,isng/2.,      &
473           truelat1,truelat2,stand_lon,dx,dy
474      !PDEF for V
475      if (test2 .ne. -1) write(12,'("pdef ",i4," ",i3," lcc ",f7.3," ",f8.3," ",     &
476&          f8.3," ",f8.3," ",f9.5," ",f9.5," ",f10.5," ",         &
477&          f10.3," ",f10.3)')                                     &
478           iweg-1,isng,cen_lat,cen_lon,iweg/2.,(isng+1.)/2.,      &
479           truelat1,truelat2,stand_lon,dx,dy
480      if (have_vert_stag) then
481       !PDEF for W
482       write(13,'("pdef ",i4," ",i3," lcc ",f7.3," ",f8.3," ",    &
483&           f8.3," ",f8.3," ",f9.5," ",f9.5," ",f10.5," ",        &
484&           f10.3," ",f10.3)')                                    &
485            iweg-1,isng-1,cen_lat,cen_lon,iweg/2.,isng/2.,        &
486            truelat1,truelat2,stand_lon,dx,dy
487      endif
488
489      !make sure truelat1 is the larger number
490       if (truelat1 < truelat2) then
491          temp = truelat1
492          truelat1 = truelat2
493          truelat2 = temp
494       endif
495
496       xlon_a(1) = xlon_M(1,1,1,1)
497       xlon_a(2) = xlon_M(1,isng-1,1,1)
498       xlon_a(3) = xlon_M(iweg-1,1,1,1)
499       xlon_a(4) = xlon_M(iweg-1,isng-1,1,1)
500
501       !check for dateline
502       ilon = 0
503       if ( abs(xlon_a(1) - xlon_a(3)) .GT. 180. ) ilon = 1
504       if ( abs(xlon_a(2) - xlon_a(4)) .GT. 180. ) ilon = 1
505
506       abslatmin = minval(xlat_M)
507       abslatmax = maxval(xlat_M)
508       abslonmin =  99999.
509       abslonmax = -99999.
510
511       do i=1,4
512         IF ( xlon_a(i) .lt. 0.0 .AND. ilon .eq. 1 ) THEN
513           abslonmin=min(abslonmin,360.+xlon_a(i))
514           abslonmax=max(abslonmax,360.+xlon_a(i))
515         ELSE
516           abslonmin=min(abslonmin,xlon_a(i))
517           abslonmax=max(abslonmax,xlon_a(i))
518         ENDIF
519       enddo
520
521       dxll = (dx/1000.)/59./2.  !!MARS MARS MARS
522       slack = ((dx/1000.)*3.)/100.
523       ipoints = int(((abslonmax-abslonmin)+(3.*slack))/dxll)
524       jpoints = int(((abslatmax-abslatmin)+(3.*slack))/dxll)
525
526       do i = 10,end_ctl_files
527         write(i,'("xdef ",i4," linear ",f10.5," ",f12.8)') ipoints, &
528                     (abslonmin-1.5*slack),dxll
529         write(i,'("ydef ",i4," linear ",f10.5," ",f12.8)') jpoints, &
530                     (abslatmin-1.5*slack),dxll
531       enddo
532   ENDIF   
533
534   IF (map_proj == 2) THEN                                     !! Polar Stereo Projection
535     if (debug) write(*,*) "We have Polar Stereo Projection data"
536
537       ipole = 1
538       if (truelat1 .lt. 0.0) ipole = -1
539       radpd  = .01745329
540       !earthr = 6.3712E6
541       earthr = 3.3972E6       
542 
543       dx_ps = dx/( (1.0+sin(ipole*TRUELAT1*radpd))/(1.0+sin(60*radpd)) ) !! true dx (meter) at 60degrees
544       re    = ( earthr * (1.0 + sin(60*radpd)) ) / dx_ps          !!
545
546
547       IF ( ipole == 1 ) THEN
548
549         rlong  = (xlon_M(1,1,1,1) - stand_lon) * radpd
550         rlat   = xlat_M(1,1,1,1) * radpd
551         r    =  (re * cos(rlat)) / (1.0 + sin(rlat))
552         xi   =  (1. - r * sin(rlong))
553         xj   =  (1. + r * cos(rlong))
554
555         !PDEF for M
556         write(10,'("pdef ",i3," ",i3," nps ",f8.3," ",f8.3," ", &
557&              f12.4," ",f12.7)') iweg-1, isng-1,                &
558               xi, xj, stand_lon, dx_ps*0.001
559         !PDEF for U
560         if (test1 .ne. -1) write(11,'("pdef ",i3," ",i3," nps ",f8.3," ",f8.3," ", &
561&              f12.4," ",f12.7)') iweg, isng-1,                  &
562               xi, xj, stand_lon, dx_ps*0.001
563         !PDEF for V
564         if (test2 .ne. -1) write(12,'("pdef ",i3," ",i3," nps ",f8.3," ",f8.3," ", &
565&              f12.4," ",f12.7)') iweg-1, isng,                  &
566               xi, xj, stand_lon, dx_ps*0.001
567         if (have_vert_stag) then
568           !PDEF for W
569           write(13,'("pdef ",i3," ",i3," nps ",f8.3," ",f8.3," ", &
570&                f12.4," ",f12.7)') iweg-1, isng-1,                &
571                 xi, xj, stand_lon, dx_ps*0.001
572         endif 
573
574       ELSE IF ( ipole == -1 ) THEN
575
576         rlong  = (180.0 + xlon_M(1,1,1,1) - stand_lon) * radpd
577         rlat   = ipole*xlat_M(1,1,1,1) * radpd
578         r    =  (re * cos(rlat)) / (1.0 + sin(rlat))
579         xi   =  (1. - ipole * r * sin(rlong))
580         xj   =  (1. + r * cos(rlong))
581
582         !PDEF for M
583         write(10,'("pdef ",i3," ",i3," sps ",f8.3," ",f8.3," ", &
584&              f12.4," ",f12.7)') iweg-1, isng-1,                &
585               xi, xj, (180+stand_lon), -0.001*dx_ps
586         !PDEF for U
587         if (test1 .ne. -1) write(11,'("pdef ",i3," ",i3," sps ",f8.3," ",f8.3," ", &
588&              f12.4," ",f12.7)') iweg, isng-1,                  &
589               xi, xj, (180+stand_lon), -0.001*dx_ps
590         !PDEF for V
591         if (test2 .ne. -1) write(12,'("pdef ",i3," ",i3," sps ",f8.3," ",f8.3," ", &
592&              f12.4," ",f12.7)') iweg-1, isng,                  &
593               xi, xj, (180+stand_lon), -0.001*dx_ps
594         if (have_vert_stag) then
595         !PDEF for W
596           write(13,'("pdef ",i3," ",i3," sps ",f8.3," ",f8.3," ", &
597&                f12.4," ",f12.7)') iweg-1, isng-1,                &
598                 xi, xj, (180+stand_lon), -0.001*dx_ps
599         endif
600 
601       END IF
602
603
604       abslonmin = minval(xlon_M)
605       abslonmax = maxval(xlon_M)
606       abslatmin = minval(xlat_M)
607       abslatmax = maxval(xlat_M)
608
609       dxll = (dx_ps/1000.)/59./2.   !! MARS MARS
610       ipoints = int(((abslonmax-abslonmin)+2.)/dxll)-2
611       jpoints = int(((abslatmax-abslatmin)+2.)/dxll)-3
612
613       ref_lon = abslonmin-1.
614       IF ( abslonmax .ge. 175. .AND. abslonmax .le. 180.) THEN
615         IF ( abslonmin .ge. -180. .AND. abslonmin .le. -175.) THEN
616           ref_lon = 0.0
617         END IF
618       END IF
619
620       ref_lat = abslatmin-1.
621       IF ( ipole == -1 ) ref_lat = abslatmin
622
623
624       do i = 10,end_ctl_files
625         write(i,'("xdef ",i4," linear ",f6.1," ",f12.8)') ipoints, &
626                     ref_lon,dxll
627         write(i,'("ydef ",i4," linear ",f6.1," ",f12.8)') jpoints, &
628                     ref_lat,dxll
629       enddo
630
631   ENDIF   
632
633   IF (map_proj == 3) THEN                                     !! Mercator Projection
634       if (debug) write(*,*) "We have Mercator Projection data"
635
636       !check to see if we have the corner lat/lon
637       if ( allocated(xlat_U) ) then
638         lat_ll_u = xlat_U(1,1,1,1)
639         lat_ur_u = xlat_U(iweg,isng-1,1,1)
640       endif
641       if ( allocated(xlat_V) ) then
642         lat_ll_v = xlat_V(1,1,1,1)
643         lat_ur_v = xlat_V(iweg-1,isng,1,1)
644       endif
645       if ( allocated(xlon_U) ) then
646         lon_ll_u = xlon_U(1,1,1,1)
647         lon_ur_u = xlon_U(iweg,isng-1,1,1)
648       endif
649       if ( allocated(xlon_V) ) then
650         lon_ll_v = xlon_V(1,1,1,1)
651         lon_ur_v = xlon_V(iweg-1,isng,1,1)
652       endif
653
654       if ( lat_ll_u==-999.99 .or. lat_ur_u==-999.99 .or. lat_ll_v==-999.99 .or. lat_ur_v==-999.99 .or. &
655            lon_ll_u==-999.99 .or. lon_ur_u==-999.99 .or. lon_ll_v==-999.99 .or. lon_ur_v==-999.99 ) then ! try getting them from the global att
656          if (debug) write(*,*) ' NO STAGGERED LAT/LON DATA AVAILBLE - try and get it from meta data '
657          rcode = nf_get_att_real(cdfid, nf_global, 'corner_lons', lons16)
658          if ( rcode == 0 ) then     
659            rcode = nf_get_att_real(cdfid, nf_global, 'corner_lats', lats16)
660            lat_ll_u = lats16( 5)
661            lat_ur_u = lats16( 7)
662            lat_ll_v = lats16( 9)
663            lat_ur_v = lats16(11) 
664            lon_ll_u = lons16( 5)
665            lon_ur_u = lons16( 7)
666            lon_ll_v = lons16( 9)
667            lon_ur_v = lons16(11)
668          else
669             write(*,*) ' NO STAGGERED LAT/LON DATA AVAILBLE - FAKE IT'
670            lat_ll_u = xlat_M(1,1,1,1)           - abs(xlat_M(2,1,1,1)-xlat_M(1,1,1,1))/2.0
671            lat_ur_u = xlat_M(iweg-1,isng-1,1,1) + abs(xlat_M(2,1,1,1)-xlat_M(1,1,1,1))/2.0
672            lat_ll_v = xlat_M(1,1,1,1)           - abs(xlat_M(1,2,1,1)-xlat_M(1,1,1,1))/2.0
673            lat_ur_v = xlat_M(iweg-1,isng-1,1,1) + abs(xlat_M(1,2,1,1)-xlat_M(1,1,1,1))/2.0
674            lon_ll_u = xlon_M(1,1,1,1)           - abs(xlon_M(2,1,1,1)-xlon_M(1,1,1,1))/2.0
675            lon_ur_u = xlon_M(iweg-1,isng-1,1,1) + abs(xlon_M(2,1,1,1)-xlon_M(1,1,1,1))/2.0
676            lon_ll_v = xlon_M(1,1,1,1)           - abs(xlon_M(1,2,1,1)-xlon_M(1,1,1,1))/2.0
677            lon_ur_v = xlon_M(iweg-1,isng-1,1,1) + abs(xlon_M(1,2,1,1)-xlon_M(1,1,1,1))/2.0
678          endif
679       endif
680
681       !check for dateline
682       ilon = 0
683       if ( abs(xlon_M(1,1,1,1)      - xlon_M(iweg-1,1,1,1))      .GT. 180. ) ilon = 1
684       if ( abs(xlon_M(1,isng-1,1,1) - xlon_M(iweg-1,isng-1,1,1)) .GT. 180. ) ilon = 1
685
686       IF ( ilon == 1 ) THEN
687
688         !! For M
689         WRITE(10,'("xdef ",i4," linear ",f9.4," ",f8.4)')   &
690            iweg-1,xlon_M(1,1,1,1), &
691            (abs(xlon_M(1,1,1,1)-(360.+xlon_M(iweg-1,isng-1,1,1)))/(iweg-2))
692         !! For U
693         if (test1 .ne. -1) WRITE(11,'("xdef ",i4," linear ",f9.4," ",f8.4)')   &
694            iweg,lon_ll_u, &
695            (abs(lon_ll_u-(360.+lon_ur_u))/(iweg-1))
696         !! For V
697         if (test2 .ne. -1) WRITE(12,'("xdef ",i4," linear ",f9.4," ",f8.4)')   &
698            iweg-1,lon_ll_v, &
699            (abs(lon_ll_v-(360.+lon_ur_v))/(iweg-2))
700         if (have_vert_stag) then
701           !! For W
702           WRITE(13,'("xdef ",i4," linear ",f9.4," ",f8.4)')   &
703              iweg-1,xlon_M(1,1,1,1), &
704              (abs(xlon_M(1,1,1,1)-(360.+xlon_M(iweg-1,isng-1,1,1)))/(iweg-2))
705         endif
706
707       ELSE
708
709         IF ( mercator_defs .AND. allocated(xlon_U) .AND. allocated(xlon_V) ) THEN
710           WRITE(10,'("xdef ",i4," levels ")') iweg-1  !! M
711           if (test1 .ne. -1) WRITE(11,'("xdef ",i4," levels ")') iweg    !! U
712           if (test2 .ne. -1) WRITE(12,'("xdef ",i4," levels ")') iweg-1  !! V
713           if (have_vert_stag) WRITE(13,'("xdef ",i4," levels ")') iweg-1  !! W
714           DO i = 1,iweg-1
715              WRITE(10,*) xlon_M(i,1,1,1)
716           END DO
717           DO i = 1,iweg
718              if (test1 .ne. -1) WRITE(11,*) xlon_U(i,1,1,1)
719           END DO
720           DO i = 1,iweg-1
721              if (test2 .ne. -1) WRITE(12,*) xlon_V(i,1,1,1)
722           END DO
723           if (have_vert_stag) then
724             DO i = 1,iweg-1
725                WRITE(13,*) xlon_M(i,1,1,1)
726             END DO
727           endif 
728         ELSE
729           !! For M
730           WRITE(10,'("xdef ",i4," linear ",f9.4," ",f8.4)')   &
731              iweg-1,xlon_M(1,1,1,1), &
732              (abs(xlon_M(1,1,1,1)-(xlon_M(iweg-1,isng-1,1,1)))/(iweg-2))
733           !! For U
734           if (test1 .ne. -1) WRITE(11,'("xdef ",i4," linear ",f9.4," ",f8.4)')   &
735              iweg,lon_ll_u, &
736              (abs(lon_ll_u-(lon_ur_u))/(iweg-1))
737           !! For V
738           if (test2 .ne. -1) WRITE(12,'("xdef ",i4," linear ",f9.4," ",f8.4)')   &
739              iweg-1,lon_ll_v, &
740              (abs(lon_ll_v-(lon_ur_v))/(iweg-2))
741           if (have_vert_stag) then
742             !! For W
743             WRITE(13,'("xdef ",i4," linear ",f9.4," ",f8.4)')   &
744                iweg-1,xlon_M(1,1,1,1), &
745                (abs(xlon_M(1,1,1,1)-(xlon_M(iweg-1,isng-1,1,1)))/(iweg-2))
746           endif
747         END IF
748
749       ENDIF
750
751       IF ( mercator_defs .AND. allocated(xlat_U) .AND. allocated(xlat_V) ) THEN
752          WRITE(10,'("ydef ",i4," levels ")') isng-1  !! M
753          if (test1 .ne. -1) WRITE(11,'("ydef ",i4," levels ")') isng-1  !! U
754          if (test2 .ne. -1) WRITE(12,'("ydef ",i4," levels ")') isng    !! V
755          if (have_vert_stag) WRITE(13,'("ydef ",i4," levels ")') isng-1  !! W
756          DO i = 1,isng-1
757            WRITE(10,*) xlat_M(1,i,1,1)
758          END DO
759          DO i = 1,isng-1
760            if (test1 .ne. -1) WRITE(11,*) xlat_U(1,i,1,1)
761          END DO
762          DO i = 1,isng
763            if (test2 .ne. -1) WRITE(12,*) xlat_V(1,i,1,1)
764          END DO
765          if (have_vert_stag) then
766            DO i = 1,isng-1
767              WRITE(13,*) xlat_M(1,i,1,1)
768            END DO
769          endif 
770       ELSE
771          !! For M
772          WRITE(10,'("ydef ",i4," linear ",f9.4," ",f8.4)')   &
773                    isng-1,xlat_M(1,1,1,1),(abs(xlat_M(1,1,1,1)-xlat_M(iweg-1,isng-1,1,1))/(isng-2))
774          !! For U
775          if (test1 .ne. -1) WRITE(11,'("ydef ",i4," linear ",f9.4," ",f8.4)')   &
776                    isng-1,lat_ll_u,(abs(lat_ll_u-lat_ur_u)/(isng-2))
777          !! For V
778          if (test2 .ne. -1) WRITE(12,'("ydef ",i4," linear ",f9.4," ",f8.4)')   &
779                    isng-1,lat_ll_v,(abs(lat_ll_v-lat_ur_v)/(isng-2))
780          if (have_vert_stag) then
781            !! For W
782            WRITE(13,'("ydef ",i4," linear ",f9.4," ",f8.4)')   &
783                      isng-1,xlat_M(1,1,1,1),(abs(xlat_M(1,1,1,1)-xlat_M(iweg-1,isng-1,1,1))/(isng-2))
784          endif
785       ENDIF
786
787   ENDIF
788
789
790
791   if (have_vert_stag) then
792     !FOR M
793       write(10,'("zdef  ",i3, " levels  ")') ibtg-1
794       do k = 1,ibtg-1
795          write(10,'(f10.5)') znu(k,1,1,1)
796       enddo
797     !FOR U
798       if (test1 .ne. -1) write(11,'("zdef  ",i3, " levels  ")') ibtg-1
799       do k = 1,ibtg-1
800          if (test1 .ne. -1) write(11,'(f10.5)') znu(k,1,1,1)
801       enddo
802     !FOR V
803       if (test2 .ne. -1) write(12,'("zdef  ",i3, " levels  ")') ibtg-1
804       do k = 1,ibtg-1
805          if (test2 .ne. -1) write(12,'(f10.5)') znu(k,1,1,1)
806       enddo
807     !FOR W
808       write(13,'("zdef  ",i3, " levels  ")') ibtg
809       do k = 1,ibtg
810          write(13,'(f10.5)') znw(k,1,1,1)
811       enddo
812   else
813     do i = 10,end_ctl_files
814       write(i,'("zdef  ",i3," linear 1 1")') ibtg
815     enddo
816   endif
817
818
819   !TDEF
820   do i = 10,end_ctl_files
821     write(i,'("tdef  ",A)') tdef
822   enddo
823
824
825   if (debug) write(*,*) "WRITING variables to .ctl files" 
826   !VARS
827   write(10,'("vars  ",i3)') i_M
828   do i = 1,i_M
829      write(10,'(A)') trim(M_vars(i))
830   enddo
831   if (test1 .ne. -1) write(11,'("vars  ",i3)') i_U
832   do i = 1,i_U
833      if (test1 .ne. -1) write(11,'(A)') trim(U_vars(i))
834   enddo
835   if (test2 .ne. -1) write(12,'("vars  ",i3)') i_V
836   do i = 1,i_V
837      if (test2 .ne. -1) write(12,'(A)') trim(V_vars(i))
838   enddo
839   if (have_vert_stag) then
840     write(13,'("vars  ",i3)') i_W
841     do i = 1,i_W
842        write(13,'(A)') trim(W_vars(i))
843     enddo
844   endif
845
846
847   do i = 10,end_ctl_files
848      write(i,'("endvars")') 
849   enddo
850   CALL SYSTEM ( '( rm -rf fort.11 fort.12 )' )
851
852!===========================================================================================
853
854  rcode = nf_close(cdfid) 
855
856  write(*,*) "Successful "
857  write(*,*) "================================================ "
858  print*,"  "
859
860  end program WRFnc2ctl
861!------------------------------------------------------------------------------
862
863  subroutine time_calc( times, ntimes, tdef, debug )
864
865  character (len=19)             :: times(ntimes), time
866  character (len=3)              :: mth, t_ind
867  character (len=40)             :: tdef
868  integer                        :: year,month,day,hh1,hh2,mn1,mn2,seconds,hourint,minsint,tdiff
869  integer                        :: ntimes
870  logical                        :: debug
871
872
873!! Time comes in as YYYY-MM-DD_HH:MM:SS
874!!                  1234 67 90 23 45 89
875   time = times(1)
876   read(time(1:4),*)   year
877   read(time(6:7),*)   month
878   read(time(9:10),*)  day
879   read(time(12:13),*) hh1
880   read(time(15:16),*) mn1
881
882    if (month == 0) return     
883    if (month == 1) mth = 'jan'
884    if (month == 2) mth = 'feb'
885    if (month == 3) mth = 'mar'
886    if (month == 4) mth = 'apr'
887    if (month == 5) mth = 'may'
888    if (month == 6) mth = 'jun'
889    if (month == 7) mth = 'jul'
890    if (month == 8) mth = 'aug'
891    if (month == 9) mth = 'sep'
892    if (month ==10) mth = 'oct'
893    if (month ==11) mth = 'nov'
894    if (month ==12) mth = 'dec'
895
896   !if (debug)
897     write(*,*) "Start date is: ",year,"-",month,"-",day,"_",hh1,":",mn1
898
899   if (ntimes .ge. 2) then
900     time = times(2)
901     read(time(12:13),*) hh2
902     read(time(15:16),*) mn2
903     hourint = abs(hh2-hh1)
904     minsint = abs(mn2-mn1)
905     if (hourint == 0 ) then
906       tdiff = minsint                 
907       t_ind = 'mn'                       
908     else
909       tdiff = hourint                 
910       t_ind = 'hr'                       
911     endif
912   else
913     tdiff = 1                 
914     t_ind = 'hr'                       
915   endif
916
917    tdef = '                               '
918   
919    !!!! TEMPORARY STUFF
920    !'00z01jan2000  1hr'
921    !hh1 = 00
922    !day = 01
923    !mth = 'jan'
924    !year = 2000
925    !tdiff = 1
926    !t_ind = 'hr'   
927    !!!! TEMPORARY STUFF
928
929    write (tdef,'(i4," linear ",i2,"z",i2,A,i4,"  ",i3,A)') ntimes, hh1, day, mth, year, tdiff, t_ind
930    !write (tdef,'(i4," linear ",A,"Z",A,A,i4,"  ",i3,A)') ntimes, time(12:13), time(9:10), mth, year, tdiff, t_ind
931    !! PB de +1 dans la date
932
933    !write (tdef,'(i4, " linear 00Z01jan2000  1hr")') ntimes
934    !if (debug)
935    write(*,*) "TDEF is set to: ",tdef
936
937  end subroutine time_calc
938
939!------------------------------------------------------------------------------
940
941  subroutine help_info
942
943  print*," "
944  print*," WRFnc2ctl   -i wrf_netcdf_input_file   -o ctl_output_name"
945  print*," "
946  print*," Current options available are:"
947  print*," -help     : Print this information"                           
948  print*," -h        : Print this information"                           
949  print*," -i        : WRF input file (netcdf format)"                 
950  print*," -o        : ctl output name - default is ARWout"                                 
951  print*," -mercator : Use if mercator projection is stretched"         
952  print*," -debug    : Print some debug information"                   
953
954  STOP
955
956  end subroutine help_info
957
958!------------------------------------------------------------------------------
959  subroutine read_args(input_file,case,mercator_defs,debug) 
960
961  implicit none
962  character (len=200)   :: input_file                       
963  character (len=80)    :: case                       
964  logical               :: mercator_defs, debug
965
966  integer, external     :: iargc
967  integer               :: numarg, i
968  character (len=80)    :: dummy
969
970
971! set up some defaults first
972  input_file = " "
973  case       = "ARWout"
974  numarg     = iargc()
975  i          = 1
976  mercator_defs = .FALSE.
977  debug = .FALSE.
978
979
980  if (numarg == 0) call help_info
981
982  do while (i <= numarg)
983    call getarg(i,dummy)
984
985    if (dummy(1:1) == "-") then    ! We have an option, else it is the filename
986
987      SELECTCASE (trim(dummy))
988          CASE ("-help")
989               call help_info
990          CASE ("-h")
991               call help_info
992          CASE ("-i")
993               i = i+1
994               call getarg(i,dummy)         ! read input_file name
995               input_file = dummy
996          CASE ("-o")
997               i = i+1
998               call getarg(i,dummy)         ! read case name
999               case = dummy
1000          CASE ("-mercator")
1001               mercator_defs = .TRUE.
1002          CASE ("-debug")
1003               debug = .TRUE.
1004          CASE DEFAULT
1005               call help_info
1006      END SELECT
1007    else
1008      input_file = dummy                     ! if option -i was not used for input file
1009    endif
1010
1011      i = i+1
1012
1013  enddo
1014
1015  if (input_file == " ") call help_info
1016
1017
1018  end subroutine read_args
1019
1020!------------------------------------------------------------------------------
1021      subroutine handle_err(message,nf_status)
1022      include "netcdf.inc"
1023      integer                 :: nf_status
1024      character (len=80)      :: message
1025      if (nf_status .ne. nf_noerr) then
1026         write(*,*)  'ERROR: ', trim(message)
1027         STOP
1028      endif
1029      end subroutine handle_err
1030!------------------------------------------------------------------------------
1031
Note: See TracBrowser for help on using the repository browser.