source: trunk/INTERP2/make_weight.pro @ 2

Last change on this file since 2 was 2, checked in by pinsard, 17 years ago

initial import from /usr/work/fvi/OPA/geomag/

File size: 15.3 KB
Line 
1;------------------------------------------------------------
2;------------------------------------------------------------
3;------------------------------------------------------------
4;
5; NAME: make_weight.pro
6;
7; PURPOSE: Creates for a given ORCA grid, and its ORCA_GEO equivalent
8;          grid  the required
9;
10; CATEGORY : Subroutine
11;
12; CALLING SEQUENCE : make_weight
13;
14; INPUTS :
15;         All the variables of a given ORCA meshmask file and its ORCA_GEO equivalent
16; OUTPUTS :
17;         Weights and pointers required to make interpolation between
18;         a field on the ORCA T grid and a field on the ORCA_GEO T grid
19;
20; COMMON BLOCKS:
21;         None
22;
23; SIDE EFFECTS:
24;
25; RESTRICTIONS:
26;
27; EXAMPLE:
28;
29; MODIFICATION HISTORY: 08/2002 Robinson Hordoir
30;
31;------------------------------------------------------------
32;------------------------------------------------------------
33;------------------------------------------------------------
34
35pro make_weight
36
37
38@initorca
39@naminterp2
40@init_path2
41
42
43;*****************************************
44; Error file opening
45;*****************************************
46
47openu, lun_err,errorfile,/get_lun
48
49;************************
50; Compute Geogrid
51;************************
52; Let's get our reference latitude
53
54northernmost=where(gphit EQ max(gphit))
55
56jnorth=northernmost/jpi
57
58inorth=min(northernmost-jnorth*jpi)
59
60ref_latitude=gphit(inorth,*)
61
62; Let's get the reference longitude
63
64equador=where(abs(gphit) EQ  min(abs(gphit)) )
65
66ref_longitude=glamt(equador)
67
68; Now we have the new grid
69
70
71glamt_reg=ref_longitude#replicate(1,n_elements(ref_latitude))
72
73
74gphit_reg=replicate(1,n_elements(ref_longitude))#ref_latitude
75
76; Get the output mask
77
78mask_reg=read_ncdf('tmask',FILENAME=outputgrid,/NOSTRUCT,/cont_nofill)
79
80; Get the latitude and longitude with extra boundary lines
81
82get_fullmesh, glamt_full,gphit_full
83
84sph_coord_full=replicate(9999.,3,jpiglo*jpjglo)
85
86sph_coord_full(0,*)=reform(glamt_full,jpiglo*jpjglo)
87sph_coord_full(1,*)=reform(gphit_full,jpiglo*jpjglo)
88sph_coord_full(2,*)=replicate(1.,jpiglo*jpjglo)
89
90rec_grid_full=cv_coord(/DEGREES,FROM_SPHERE=sph_coord_full,/TO_RECT)
91
92xx_grid_full=rec_grid_full(0,*)
93yy_grid_full=rec_grid_full(1,*)   
94zz_grid_full=rec_grid_full(2,*)
95
96
97sph_coord=replicate(9999.,3,jpi*jpj)
98
99radius=replicate(1.,jpi,jpj)
100radius(0:jpi/2,jpj-1)=1000.
101
102sph_coord(0,*)=reform(glamt,jpi*jpj)
103sph_coord(1,*)=reform(gphit,jpi*jpj)
104sph_coord(2,*)=reform(radius,jpi*jpj)
105
106rec_grid=cv_coord(/DEGREES,FROM_SPHERE=sph_coord,/TO_RECT)
107
108xx_grid=rec_grid(0,*)
109yy_grid=rec_grid(1,*)   
110zz_grid=rec_grid(2,*)
111
112
113
114; Init array to stock pointers
115
116  pointer_Tx=replicate(999,jpi,jpj,4)
117  pointer_Ty=replicate(999,jpi,jpj,4)
118; Init weight matric
119
120  weight_T=double(replicate(0.,jpi,jpj,4))
121
122; Get the latitude to start from
123
124not_found=1
125for jj=0, jpj-1 do begin
126
127   lat_max_ORCA=max( gphit(*,jj)     )
128   lat_max_reg =max( gphit_reg(*,jj) )
129
130
131   if ((lat_max_ORCA NE lat_max_reg) AND (not_found EQ 1)) then begin
132   jpj_start=jj-2
133   not_found=0
134   endif
135
136end
137
138
139    for jj=0,jpj-1 do begin
140 
141    for ji=0,jpi-1 do begin
142
143
144    if mask_reg(ji,jj) NE 0 then begin
145
146; Coordinates of our point
147
148    xx=glamt_reg(ji,jj)
149    yy=gphit_reg(ji,jj)
150   
151    sph_coord=replicate(1.,3)
152    sph_coord(0)=xx
153    sph_coord(1)=yy
154
155; Convert them into cartesian
156   
157    rec_coord=cv_coord(/DEGREES,FROM_SPHERE=sph_coord,/TO_RECT)
158
159    xx_bis=replicate(rec_coord(0),jpi*jpj)
160    yy_bis=replicate(rec_coord(1),jpi*jpj)
161    zz_bis=replicate(rec_coord(2),jpi*jpj)
162
163
164   
165
166;****************************************************************
167; Finds in which grid cell point stand   
168;****************************************************************
169   
170    xx_bis=reform(xx_bis,jpi,jpj)
171    yy_bis=reform(yy_bis,jpi,jpj)
172    zz_bis=reform(zz_bis,jpi,jpj)
173
174    xx_grid_full=reform(xx_grid_full,jpiglo,jpjglo)
175    yy_grid_full=reform(yy_grid_full,jpiglo,jpjglo)
176    zz_grid_full=reform(zz_grid_full,jpiglo,jpjglo)
177
178   
179    check_inside,xx_bis ,yy_bis ,zz_bis,$
180                 xx_grid_full(1:jpiglo-2,0:jpjglo-2),yy_grid_full(1:jpiglo-2,0:jpjglo-2),zz_grid_full(1:jpiglo-2,0:jpjglo-2),$
181                 xx_grid_full(2:jpiglo-1,0:jpjglo-2),yy_grid_full(2:jpiglo-1,0:jpjglo-2),zz_grid_full(2:jpiglo-1,0:jpjglo-2),$
182                 xx_grid_full(2:jpiglo-1,1:jpjglo-1),yy_grid_full(2:jpiglo-1,1:jpjglo-1),zz_grid_full(2:jpiglo-1,1:jpjglo-1),$
183                 xx_grid_full(1:jpiglo-2,1:jpjglo-1),yy_grid_full(1:jpiglo-2,1:jpjglo-1),zz_grid_full(1:jpiglo-2,1:jpjglo-1), answer
184
185   
186     
187    if answer(0) EQ -1 then begin
188
189    print, 'ERROR IN MAKE WEIGHT : POSITION NOT LOCATED ON GRID'
190    stop
191
192    endif
193   
194   
195    answer_j=answer/jpi
196    answer_i=answer-answer_j*jpi
197   
198   
199    answer_j=answer_j(where(answer_i EQ min(answer_i)))
200    answer_i=answer_i(where(answer_i EQ min(answer_i)))
201   
202    answer_j=answer_j(0)
203    answer_i=answer_i(0)   
204
205    pointer_Tx(ji,jj,0)=answer_i
206    pointer_Ty(ji,jj,0)=answer_j
207    point_x1=answer_i
208    point_y1=answer_j
209 
210    pointer_Tx(ji,jj,1)=answer_i+1
211    pointer_Ty(ji,jj,1)=answer_j
212    point_x2=answer_i+1
213    point_y2=answer_j
214
215    pointer_Tx(ji,jj,2)=answer_i+1
216    pointer_Ty(ji,jj,2)=answer_j+1
217    point_x3=answer_i+1
218    point_y3=answer_j+1
219
220    pointer_Tx(ji,jj,3)=answer_i
221    pointer_Ty(ji,jj,3)=answer_j+1
222    point_x4=answer_i
223    point_y4=answer_j+1
224
225   
226
227; Now we know which 4 points are located nearby our output T point.
228
229; Let's get the coordinates of this points in long and lat
230
231    glamt_ref=replicate(0.,4)
232    gphit_ref=replicate(0.,4)
233 
234 
235   
236    glamt_ref(0)=glamt_full(point_x1+1,point_y1)             
237    gphit_ref(0)=gphit_full(point_x1+1,point_y1)
238
239    glamt_ref(1)=glamt_full(point_x2+1,point_y2)             
240    gphit_ref(1)=gphit_full(point_x2+1,point_y2)
241   
242    glamt_ref(2)=glamt_full(point_x3+1,point_y3)             
243    gphit_ref(2)=gphit_full(point_x3+1,point_y3)
244
245   
246    glamt_ref(3)=glamt_full(point_x4+1,point_y4)             
247    gphit_ref(3)=gphit_full(point_x4+1,point_y4)
248 
249   
250
251; Let's make things simple for coordinates
252 
253   sph_coord_ref=replicate(0.,3,5)
254
255   sph_coord_ref(0,0)=glamt_reg(ji,jj)
256   sph_coord_ref(1,0)=gphit_reg(ji,jj)
257   sph_coord_ref(2,0)=1.
258
259   sph_coord_ref(0,1:4)=glamt_ref
260   sph_coord_ref(1,1:4)=gphit_ref
261   sph_coord_ref(2,1:4)=1.
262
263   rec_coord_ref=cv_coord(/DEGREES,FROM_SPHERE=sph_coord_ref,/TO_RECT)
264 
265   xx=rec_coord_ref(0,0)
266   yy=rec_coord_ref(1,0)
267   zz=rec_coord_ref(2,0)
268
269   x1=rec_coord_ref(0,1)
270   y1=rec_coord_ref(1,1)
271   z1=rec_coord_ref(2,1)
272 
273   x2=rec_coord_ref(0,2)
274   y2=rec_coord_ref(1,2)
275   z2=rec_coord_ref(2,2)
276
277   x3=rec_coord_ref(0,3)
278   y3=rec_coord_ref(1,3)
279   z3=rec_coord_ref(2,3)
280
281   x4=rec_coord_ref(0,4)
282   y4=rec_coord_ref(1,4)
283   z4=rec_coord_ref(2,4)
284
285   xG=total(rec_coord_ref(0,1:4))/4.
286   yG=total(rec_coord_ref(1,1:4))/4.
287   zG=total(rec_coord_ref(2,1:4))/4.
288   
289   rec_coord_ref(0,*)=rec_coord_ref(0,*)-xG
290   rec_coord_ref(1,*)=rec_coord_ref(1,*)-yG
291   rec_coord_ref(2,*)=rec_coord_ref(2,*)-zG
292
293   XT=rec_coord_ref(*,1:4)
294 
295   XM=MATRIX_MULTIPLY(XT,XT,/BTRANSPOSE)
296   EVAL=HQR(ELMHES(XM),/DOUBLE)
297   EVEC=EIGENVEC(XM,EVAL)
298
299   EVAL=FLOAT(EVAL)
300   EVEC=FLOAT(EVEC)
301   ORDER=SORT(EVAL) 
302   
303
304   xx=total(rec_coord_ref(*,0)*EVEC(*,ORDER(1)))
305   yy=total(rec_coord_ref(*,0)*EVEC(*,ORDER(2)))
306
307   x1=total(rec_coord_ref(*,1)*EVEC(*,ORDER(1)))
308   y1=total(rec_coord_ref(*,1)*EVEC(*,ORDER(2)))
309
310   x2=total(rec_coord_ref(*,2)*EVEC(*,ORDER(1)))
311   y2=total(rec_coord_ref(*,2)*EVEC(*,ORDER(2)))
312
313   x3=total(rec_coord_ref(*,3)*EVEC(*,ORDER(1)))
314   y3=total(rec_coord_ref(*,3)*EVEC(*,ORDER(2)))
315
316   x4=total(rec_coord_ref(*,4)*EVEC(*,ORDER(1)))
317   y4=total(rec_coord_ref(*,4)*EVEC(*,ORDER(2)))
318   
319
320;**************************************************
321; Computes i coordinate of point in the grid cell *
322;**************************************************
323
324; Calculate coordinates of j vector
325
326; Vector j0 is vector AD/(AD)
327
328xj0=x4-x1
329
330yj0=y4-y1
331
332; Normalised vector
333
334norm_j0=sqrt(xj0*xj0+yj0*yj0)
335
336xj0=xj0/norm_j0
337
338yj0=yj0/norm_j0
339
340
341; Vector j1 is vector BC/(BC)
342
343xj1=x3-x2
344
345yj1=y3-y2
346
347; Normalised vector
348
349norm_j1=sqrt(xj1*xj1+yj1*yj1)
350
351xj1=xj1/norm_j1
352
353yj1=yj1/norm_j1
354
355; Caculate coefficient for solving 2nd degree equation in i
356
357a_coeff= (x2-x1)*(yj0-yj1)-(y2-y1)*(xj0-xj1)
358
359b_coeff= -(x2-x1)*yj0-(xx-x1)*(yj0-yj1)+(yy-y1)*(xj0-xj1)+(y2-y1)*xj0
360   
361c_coeff= (xx-x1)*yj0-(yy-y1)*xj0
362
363
364delta=max([b_coeff*b_coeff-4*a_coeff*c_coeff,0])
365
366
367
368iP=999.
369
370if abs(a_coeff) GT 1.E-5 then begin
371
372sol1=(-b_coeff+sqrt(delta))/(2.*a_coeff)
373
374sol2=(-b_coeff-sqrt(delta))/(2.*a_coeff)
375
376dist_sol1=sol1*sol1+(sol1-1)*(sol1-1)
377dist_sol2=sol2*sol2+(sol2-1)*(sol2-1)
378
379iP=sol1
380
381if dist_sol1 LT dist_sol2 then iP=sol1
382if dist_sol2 LT dist_sol1 then iP=sol2
383
384
385endif else begin
386
387iP=-c_coeff/b_coeff
388
389endelse
390
391
392
393
394;**************************************************
395; Computes j coordinate of point in the grid cell *
396;**************************************************
397
398; Calculate coordinates of i vector
399
400; Vector i0 is vector AB/(AB)
401
402xi0=x2-x1
403
404yi0=y2-y1
405
406; Normalised vector
407
408norm_i0=sqrt(xi0*xi0+yi0*yi0)
409if norm_i0 LT 0. then stio
410xi0=xi0/norm_i0
411
412yi0=yi0/norm_i0
413
414
415; Vector i1 is vector DC/(DC)
416
417xi1=x3-x4
418
419yi1=y3-y4
420
421; Normalised vector
422
423norm_i1=sqrt(xi1*xi1+yi1*yi1)
424if norm_i1 LT 0. then stop
425xi1=xi1/norm_i1
426
427yi1=yi1/norm_i1
428
429; Caculate coefficient for solving 2nd degree equation in j
430
431a_coeff= (x4-x1)*(yi0-yi1)-(y4-y1)*(xi0-xi1)
432
433b_coeff= -(x4-x1)*yi0-(xx-x1)*(yi0-yi1)+(yy-y1)*(xi0-xi1)+(y4-y1)*xi0
434   
435c_coeff= (xx-x1)*yi0-(yy-y1)*xi0
436
437
438delta=max([b_coeff*b_coeff-4*a_coeff*c_coeff,0])
439
440
441jP=999.
442
443if abs(a_coeff) GT 1.E-5  then begin
444
445sol1=(-b_coeff+sqrt(delta))/(2.*a_coeff)
446
447sol2=(-b_coeff-sqrt(delta))/(2.*a_coeff)
448
449dist_sol1=sol1*sol1+(sol1-1)*(sol1-1)
450dist_sol2=sol2*sol2+(sol2-1)*(sol2-1)
451
452jP=sol1
453
454if dist_sol1 LT dist_sol2 then jP=sol1
455if dist_sol2 LT dist_sol1 then jP=sol2
456
457
458endif else begin
459
460jP=-c_coeff/b_coeff
461
462endelse
463   
464   
465   weight_T(ji,jj,0) = (1.-iP)*(1.-jP)
466
467   weight_T(ji,jj,1) = iP*(1.-jP)
468
469   weight_T(ji,jj,2) = iP*jP
470
471   weight_T(ji,jj,3) = (1.-iP)*jP
472
473   
474   if min(weight_T(ji,jj,*)) LT 0 OR max(weight_T(ji,jj,*)) GT 1 then begin
475
476      printf, lun_err,'JI=',ji,' JJ=',jj
477      printf, lun_err,'IP=',iP,' JP=',jP
478
479     
480      printf, lun_err,'XP=',xx,' YP=',yy
481      printf, lun_err,'XA=',x1,' YA=',y1
482      printf, lun_err,'XB=',x2,' YB=',y2
483      printf, lun_err,'XC=',x3,' YC=',y3
484      printf, lun_err,'XD=',x4,' YD=',y4
485     
486   endif
487 
488
489   endif
490
491 
492  end
493
494end
495
496close, lun_err
497free_lun,lun_err
498
499; Convert to full format
500
501
502  weight_full=replicate(0.,jpiglo,jpjglo,4)
503  pointer_i_full=replicate(0,jpiglo,jpjglo,4)
504  pointer_j_full=replicate(0,jpiglo,jpjglo,4)
505
506  weight_full(1:jpiglo-2,0:jpjglo-2,*)=weight_T
507  pointer_i_full(1:jpiglo-2,0:jpjglo-2,*)=pointer_Tx
508  pointer_j_full(1:jpiglo-2,0:jpjglo-2,*)=pointer_Ty
509
510
511
512
513; Let's save Weight functions and the pointer arrays
514
515
516; Creates Netcdf File to save this
517
518
519   vargrid = 'T'
520
521; Name
522   idout = NCDF_CREATE(weightfile,/clobber)
523 
524   NCDF_CONTROL, idout, /nofill
525; Dimension
526   xidout = NCDF_DIMDEF(idout, 'x', jpiglo)
527   yidout = NCDF_DIMDEF(idout, 'y', jpjglo)
528   didout = NCDF_DIMDEF(idout, 'deptht', 1)
529   tidout = NCDF_DIMDEF(idout, 'time_counter', /unlimited)
530
531
532
533; Attributes   
534   id_0  = NCDF_VARDEF(idout, 'nav_lon'     , [xidout, yidout                ], /FLOAT)
535   id_1  = NCDF_VARDEF(idout, 'nav_lat'     , [xidout, yidout                ], /FLOAT)
536   id_2  = NCDF_VARDEF(idout, 'deptht'      , [                didout        ], /FLOAT)
537   id_3  = NCDF_VARDEF(idout, 'time_counter', [                        tidout], /FLOAT)
538 
539   id0   = NCDF_VARDEF(idout, 'weight_1'      , [xidout, yidout,didout,tidout ], /FLOAT)
540   id1   = NCDF_VARDEF(idout, 'weight_2'      , [xidout, yidout,didout,tidout  ], /FLOAT)
541   id2   = NCDF_VARDEF(idout, 'weight_3'      , [xidout, yidout,didout,tidout  ], /FLOAT)
542   id3   = NCDF_VARDEF(idout, 'weight_4'      , [xidout, yidout,didout,tidout  ], /FLOAT)
543   id4   = NCDF_VARDEF(idout, 'pointer_1_i'     , [xidout, yidout ,didout,tidout ], /LONG)
544   id5   = NCDF_VARDEF(idout, 'pointer_2_i'     , [xidout, yidout ,didout,tidout ], /LONG)
545   id6   = NCDF_VARDEF(idout, 'pointer_3_i'     , [xidout, yidout,didout,tidout  ], /LONG)
546   id7   = NCDF_VARDEF(idout, 'pointer_4_i'     , [xidout, yidout,didout,tidout  ], /LONG)
547   id8   = NCDF_VARDEF(idout, 'pointer_1_j'     , [xidout, yidout,didout,tidout  ], /LONG)
548   id9   = NCDF_VARDEF(idout, 'pointer_2_j'     , [xidout, yidout,didout,tidout  ], /LONG)
549   id10  = NCDF_VARDEF(idout, 'pointer_3_j'     , [xidout, yidout,didout,tidout  ], /LONG)
550   id11  = NCDF_VARDEF(idout, 'pointer_4_j'     , [xidout, yidout,didout,tidout  ], /LONG)
551
552; Variable 0
553   NCDF_ATTPUT, idout, id_0, 'units', 'degrees_east'
554   NCDF_ATTPUT, idout, id_0, 'long_name', 'Longitude'
555   NCDF_ATTPUT, idout, id_0, 'nav_model', 'Default grid'
556; Variable 1
557   NCDF_ATTPUT, idout, id_1, 'units', 'degrees_north'
558   NCDF_ATTPUT, idout, id_1, 'long_name', 'Latitude'
559   NCDF_ATTPUT, idout, id_1, 'nav_model', 'Default grid'
560; Variable 2
561   NCDF_ATTPUT, idout, id_2, 'units','meters'
562   NCDF_ATTPUT, idout, id_2, 'long_name','Depth'
563   NCDF_ATTPUT, idout, id_2, 'nav_model','Default grid'
564; Variable3
565   NCDF_ATTPUT, idout, id_3, 'units', 'seconds since 0001-01-15 12:00:00 '
566   NCDF_ATTPUT, idout, id_3, 'calendar','noleap'
567   NCDF_ATTPUT, idout, id_3, 'title', 'Time'
568   NCDF_ATTPUT, idout, id_3, 'long_name', 'Time axis'
569   NCDF_ATTPUT, idout, id_3, 'time_origin','0000-DEC-15 00:00:00'
570
571
572; Variable 0
573   NCDF_ATTPUT, idout, id0, 'units','no unit'
574   NCDF_ATTPUT, idout, id0, 'long_name','W1'
575 
576; Variable 1
577   NCDF_ATTPUT, idout, id1, 'units','no unit'
578   NCDF_ATTPUT, idout, id1, 'long_name','W2'
579; Variable 0
580   NCDF_ATTPUT, idout, id2, 'units','no unit'
581   NCDF_ATTPUT, idout, id2, 'long_name','W3'
582
583; Variable 0
584   NCDF_ATTPUT, idout, id3, 'units','no unit'
585   NCDF_ATTPUT, idout, id3, 'long_name','W4'
586
587; Variable 0
588   NCDF_ATTPUT, idout, id4, 'units','no unit'
589   NCDF_ATTPUT, idout, id4, 'long_name','i1'
590 
591; Variable 1
592   NCDF_ATTPUT, idout, id5, 'units','no unit'
593   NCDF_ATTPUT, idout, id5, 'long_name','i2'
594; Variable 0
595   NCDF_ATTPUT, idout, id6, 'units','no unit'
596   NCDF_ATTPUT, idout, id6, 'long_name','i3'
597
598; Variable 0
599   NCDF_ATTPUT, idout, id7, 'units','no unit'
600   NCDF_ATTPUT, idout, id7, 'long_name','i4'
601
602
603; Variable 0
604   NCDF_ATTPUT, idout, id8, 'units','no unit'
605   NCDF_ATTPUT, idout, id8, 'long_name','j1'
606 
607; Variable 1
608   NCDF_ATTPUT, idout, id9, 'units','no unit'
609   NCDF_ATTPUT, idout, id9, 'long_name','j2'
610; Variable 0
611   NCDF_ATTPUT, idout, id10, 'units','no unit'
612   NCDF_ATTPUT, idout, id10, 'long_name','j3'
613
614; Variable 0
615   NCDF_ATTPUT, idout, id11, 'units','no unit'
616   NCDF_ATTPUT, idout, id11, 'long_name','j4'
617
618
619
620
621
622
623
624   NCDF_CONTROL, idout, /ENDEF
625
626; Writting
627
628   NCDF_VARPUT, idout, id_0, glamt
629
630   NCDF_VARPUT, idout, id_1, gphit
631
632   prof=5.
633   NCDF_VARPUT, idout, id_2, prof
634   temps=lonarr(12)
635   for i=0,0 do temps[i] = 86400l*(julday(1+i,15,1)-julday(1,15,1))
636   NCDF_VARPUT, idout, id_3, temps(0)
637
638 
639   NCDF_VARPUT, idout, id0,weight_full(*,*,0)
640   NCDF_VARPUT, idout, id1,weight_full(*,*,1)
641   NCDF_VARPUT, idout, id2,weight_full(*,*,2)
642   NCDF_VARPUT, idout, id3,weight_full(*,*,3)
643   NCDF_VARPUT, idout, id4,pointer_i_full(*,*,0)
644   NCDF_VARPUT, idout, id5,pointer_i_full(*,*,1)
645   NCDF_VARPUT, idout, id6,pointer_i_full(*,*,2)
646   NCDF_VARPUT, idout, id7,pointer_i_full(*,*,3)
647   NCDF_VARPUT, idout, id8,pointer_j_full(*,*,0)
648   NCDF_VARPUT, idout, id9,pointer_j_full(*,*,1)
649   NCDF_VARPUT, idout, id10,pointer_j_full(*,*,2)
650   NCDF_VARPUT, idout, id11,pointer_j_full(*,*,3)
651 
652
653
654   NCDF_CLOSE, idout
655
656
657
658end
Note: See TracBrowser for help on using the repository browser.