source: LMDZ6/branches/contrails/libf/phylmd/read_pstoke0.f90 @ 5452

Last change on this file since 5452 was 5270, checked in by abarral, 2 months ago

Replace F77 netcdf library by F90 netcdf library
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 14.9 KB
Line 
1
2! $Id: read_pstoke0.f90 5270 2024-10-24 11:55:38Z aborella $
3
4
5
6SUBROUTINE read_pstoke0(irec, zrec, zkon, zkev, airefi, phisfi, t, mfu, mfd, &
7    en_u, de_u, en_d, de_d, coefh, fm_therm, en_therm, frac_impa, frac_nucl, &
8    pyu1, pyv1, ftsol, psrf)
9
10  ! ******************************************************************************
11  ! Frederic HOURDIN, Abderrahmane IDELKADI
12  ! Lecture des parametres physique stockes online necessaires pour
13  ! recalculer offline le transport des traceurs sur la meme grille que
14  ! online
15  ! A FAIRE : une seule routine au lieu de 2 (lectflux, redecoupe)!
16  ! ******************************************************************************
17
18  USE netcdf
19  USE dimphy
20  USE indice_sol_mod
21  USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, nbp_lev
22
23  IMPLICIT NONE
24
25  INTEGER kon, kev, zkon, zkev
26!  PARAMETER (kon=iim*(jjm-1)+2, kev=llm)
27  REAL :: phisfi(nbp_lon*(nbp_lat-2)+2) !phisfi(kon)
28  REAL,ALLOCATABLE :: phisfi2(:,:) !phisfi2(nbp_lon, nbp_lat)
29  REAL,ALLOCATABLE :: airefi2(:,:) !airefi2(nbp_lon, nbp_lat)
30
31  REAL :: mfu(nbp_lon*(nbp_lat-2)+2,nbp_lev) !mfu(kon, kev)
32  REAL :: mfd(nbp_lon*(nbp_lat-2)+2,nbp_lev) !mfd(kon, kev)
33  REAL :: en_u(nbp_lon*(nbp_lat-2)+2,nbp_lev) !en_u(kon, kev)
34  REAL :: de_u(nbp_lon*(nbp_lat-2)+2,nbp_lev) !de_u(kon, kev)
35  REAL :: en_d(nbp_lon*(nbp_lat-2)+2,nbp_lev) !en_d(kon, kev)
36  REAL :: de_d(nbp_lon*(nbp_lat-2)+2,nbp_lev) !de_d(kon, kev)
37  REAL :: coefh(nbp_lon*(nbp_lat-2)+2,nbp_lev) !coefh(kon, kev)
38
39  ! abd 25 11 02
40  ! Thermiques
41  REAL :: fm_therm(nbp_lon*(nbp_lat-2)+2,nbp_lev) !fm_therm(kon, kev)
42  REAL :: en_therm(nbp_lon*(nbp_lat-2)+2,nbp_lev) !en_therm(kon, kev)
43  REAL :: t(nbp_lon*(nbp_lat-2)+2,nbp_lev) !t(kon, kev)
44
45  REAL,ALLOCATABLE :: mfu2(:,:,:) !mfu2(nbp_lon, nbp_lat, kev)
46  REAL,ALLOCATABLE :: mfd2(:,:,:) !mfd2(nbp_lon, nbp_lat, kev)
47  REAL,ALLOCATABLE :: en_u2(:,:,:) !en_u2(nbp_lon, nbp_lat, kev)
48  REAL,ALLOCATABLE :: de_u2(:,:,:) !de_u2(nbp_lon, nbp_lat, kev)
49  REAL,ALLOCATABLE :: en_d2(:,:,:) !en_d2(nbp_lon, nbp_lat, kev)
50  REAL,ALLOCATABLE :: de_d2(:,:,:) !de_d2(nbp_lon, nbp_lat, kev)
51  REAL,ALLOCATABLE :: coefh2(:,:,:) !coefh2(nbp_lon, nbp_lat, kev)
52  REAL,ALLOCATABLE :: t2(:,:,:) !t2(nbp_lon, nbp_lat, kev)
53  ! Thermiques
54  REAL,ALLOCATABLE :: fm_therm2(:,:,:) !fm_therm2(nbp_lon, nbp_lat, kev)
55  REAL,ALLOCATABLE :: en_therm2(:,:,:) !en_therm2(nbp_lon, nbp_lat, kev)
56
57  REAL,ALLOCATABLE :: pl(:) !pl(kev)
58  INTEGER irec
59  INTEGER xid, yid, zid, tid
60  INTEGER zrec, zim, zjm
61  INTEGER ncrec, nckon, nckev, ncim, ncjm
62
63  REAL :: airefi(nbp_lon*(nbp_lat-2)+2) !airefi(kon)
64  CHARACTER *20 namedim
65
66  ! !! attention !!
67  ! attention il y a aussi le pb de def kon
68  ! dim de phis??
69
70  REAL :: frac_impa(nbp_lon*(nbp_lat-2)+2,nbp_lev) !frac_impa(kon, kev)
71  REAL :: frac_nucl(nbp_lon*(nbp_lat-2)+2,nbp_lev) !frac_nucl(kon, kev)
72  REAL,ALLOCATABLE :: frac_impa2(:,:,:) !frac_impa2(nbp_lon, nbp_lat, kev)
73  REAL,ALLOCATABLE :: frac_nucl2(:,:,:) !frac_nucl2(nbp_lon, nbp_lat, kev)
74  REAL :: pyu1(nbp_lon*(nbp_lat-2)+2) !pyu1(kon)
75  REAL :: pyv1(nbp_lon*(nbp_lat-2)+2) !pyv1(kon)
76  REAL,ALLOCATABLE :: pyu12(:,:), pyv12(:,:) !pyu12(nbp_lon, nbp_lat), pyv12(nbp_lon, nbp_lat)
77  REAL :: ftsol(nbp_lon*(nbp_lat-2)+2,nbp_lev) !ftsol(kon, nbsrf)
78  REAL :: psrf(nbp_lon*(nbp_lat-2)+2,nbp_lev) !psrf(kon, nbsrf)
79  REAL,ALLOCATABLE :: ftsol1(:),ftsol2(:) !ftsol1(kon), ftsol2(kon)
80  REAL,ALLOCATABLE :: ftsol3(:),ftsol4(:) !ftsol3(kon), ftsol4(kon)
81  REAL,ALLOCATABLE :: psrf1(:), psrf2(:) !psrf1(kon), psrf2(kon)
82  REAL,ALLOCATABLE :: psrf3(:), psrf4(:) !psrf3(kon), psrf4(kon)
83  REAL,ALLOCATABLE :: ftsol12(:,:) !ftsol12(nbp_lon, nbp_lat)
84  REAL,ALLOCATABLE :: ftsol22(:,:) !ftsol22(nbp_lon, nbp_lat)
85  REAL,ALLOCATABLE :: ftsol32(:,:) !ftsol32(nbp_lon, nbp_lat)
86  REAL,ALLOCATABLE :: ftsol42(:,:) !ftsol42(nbp_lon, nbp_lat)
87  REAL,ALLOCATABLE :: psrf12(:,:) !psrf12(nbp_lon, nbp_lat)
88  REAL,ALLOCATABLE :: psrf22(:,:) !psrf22(nbp_lon, nbp_lat)
89  REAL,ALLOCATABLE :: psrf32(:,:) !psrf32(nbp_lon, nbp_lat)
90  REAL,ALLOCATABLE :: psrf42(:,:) !psrf42(nbp_lon, nbp_lat)
91
92  INTEGER,SAVE :: ncidp
93  INTEGER,SAVE :: varidmfu, varidmfd, varidps, varidenu, variddeu
94  INTEGER,SAVE :: varidt
95  INTEGER,SAVE :: varidend, varidded, varidch, varidfi, varidfn
96  ! therm
97  INTEGER,SAVE :: varidfmth, varidenth
98  INTEGER,SAVE :: varidyu1, varidyv1, varidpl, varidai, varididvt
99  INTEGER,SAVE :: varidfts1, varidfts2, varidfts3, varidfts4
100  INTEGER,SAVE :: varidpsr1, varidpsr2, varidpsr3, varidpsr4
101
102  INTEGER l, i
103  INTEGER start(4), count(4), status
104  REAL rcode
105  LOGICAL,SAVE :: first=.TRUE.
106
107  ! Allocate arrays
108  kon=nbp_lon*(nbp_lat-2)+2
109  kev=nbp_lev
110
111  ALLOCATE(phisfi2(nbp_lon, nbp_lat))
112  ALLOCATE(airefi2(nbp_lon, nbp_lat))
113  ALLOCATE(mfu2(nbp_lon, nbp_lat, kev))
114  ALLOCATE(mfd2(nbp_lon, nbp_lat, kev))
115  ALLOCATE(en_u2(nbp_lon, nbp_lat, kev))
116  ALLOCATE(de_u2(nbp_lon, nbp_lat, kev))
117  ALLOCATE(en_d2(nbp_lon, nbp_lat, kev))
118  ALLOCATE(de_d2(nbp_lon, nbp_lat, kev))
119  ALLOCATE(coefh2(nbp_lon, nbp_lat, kev))
120  ALLOCATE(t2(nbp_lon, nbp_lat, kev))
121  ALLOCATE(fm_therm2(nbp_lon, nbp_lat, kev))
122  ALLOCATE(en_therm2(nbp_lon, nbp_lat, kev))
123  ALLOCATE(pl(kev))
124  ALLOCATE(frac_impa2(nbp_lon, nbp_lat, kev))
125  ALLOCATE(frac_nucl2(nbp_lon, nbp_lat, kev))
126  ALLOCATE(pyu12(nbp_lon, nbp_lat), pyv12(nbp_lon, nbp_lat))
127  ALLOCATE(ftsol1(kon), ftsol2(kon))
128  ALLOCATE(ftsol3(kon), ftsol4(kon))
129  ALLOCATE(psrf1(kon), psrf2(kon))
130  ALLOCATE(psrf3(kon), psrf4(kon))
131  ALLOCATE(ftsol12(nbp_lon, nbp_lat))
132  ALLOCATE(ftsol22(nbp_lon, nbp_lat))
133  ALLOCATE(ftsol32(nbp_lon, nbp_lat))
134  ALLOCATE(ftsol42(nbp_lon, nbp_lat))
135  ALLOCATE(psrf12(nbp_lon, nbp_lat))
136  ALLOCATE(psrf22(nbp_lon, nbp_lat))
137  ALLOCATE(psrf32(nbp_lon, nbp_lat))
138  ALLOCATE(psrf42(nbp_lon, nbp_lat))
139
140  ! ---------------------------------------------
141  ! Initialisation de la lecture des fichiers
142  ! ---------------------------------------------
143
144  IF (irec==0) THEN
145
146    rcode = nf90_open('phystoke.nc', nf90_nowrite, ncidp)
147
148    rcode = nf90_inq_varid(ncidp, 'phis', varidps)
149    PRINT *, 'ncidp,varidps', ncidp, varidps
150
151    rcode = nf90_inq_varid(ncidp, 'sig_s', varidpl)
152    PRINT *, 'ncidp,varidpl', ncidp, varidpl
153
154    rcode = nf90_inq_varid(ncidp, 'aire', varidai)
155    PRINT *, 'ncidp,varidai', ncidp, varidai
156
157    rcode = nf90_inq_varid(ncidp, 't', varidt)
158    PRINT *, 'ncidp,varidt', ncidp, varidt
159
160    rcode = nf90_inq_varid(ncidp, 'mfu', varidmfu)
161    PRINT *, 'ncidp,varidmfu', ncidp, varidmfu
162
163    rcode = nf90_inq_varid(ncidp, 'mfd', varidmfd)
164    PRINT *, 'ncidp,varidmfd', ncidp, varidmfd
165
166    rcode = nf90_inq_varid(ncidp, 'en_u', varidenu)
167    PRINT *, 'ncidp,varidenu', ncidp, varidenu
168
169    rcode = nf90_inq_varid(ncidp, 'de_u', variddeu)
170    PRINT *, 'ncidp,variddeu', ncidp, variddeu
171
172    rcode = nf90_inq_varid(ncidp, 'en_d', varidend)
173    PRINT *, 'ncidp,varidend', ncidp, varidend
174
175    rcode = nf90_inq_varid(ncidp, 'de_d', varidded)
176    PRINT *, 'ncidp,varidded', ncidp, varidded
177
178    rcode = nf90_inq_varid(ncidp, 'coefh', varidch)
179    PRINT *, 'ncidp,varidch', ncidp, varidch
180
181    ! Thermiques
182    rcode = nf90_inq_varid(ncidp, 'fm_th', varidfmth)
183    PRINT *, 'ncidp,varidfmth', ncidp, varidfmth
184
185    rcode = nf90_inq_varid(ncidp, 'en_th', varidenth)
186    PRINT *, 'ncidp,varidenth', ncidp, varidenth
187
188    rcode = nf90_inq_varid(ncidp, 'frac_impa', varidfi)
189    PRINT *, 'ncidp,varidfi', ncidp, varidfi
190
191    rcode = nf90_inq_varid(ncidp, 'frac_nucl', varidfn)
192    PRINT *, 'ncidp,varidfn', ncidp, varidfn
193
194    rcode = nf90_inq_varid(ncidp, 'pyu1', varidyu1)
195    PRINT *, 'ncidp,varidyu1', ncidp, varidyu1
196
197    rcode = nf90_inq_varid(ncidp, 'pyv1', varidyv1)
198    PRINT *, 'ncidp,varidyv1', ncidp, varidyv1
199
200    rcode = nf90_inq_varid(ncidp, 'ftsol1', varidfts1)
201    PRINT *, 'ncidp,varidfts1', ncidp, varidfts1
202
203    rcode = nf90_inq_varid(ncidp, 'ftsol2', varidfts2)
204    PRINT *, 'ncidp,varidfts2', ncidp, varidfts2
205
206    rcode = nf90_inq_varid(ncidp, 'ftsol3', varidfts3)
207    PRINT *, 'ncidp,varidfts3', ncidp, varidfts3
208
209    rcode = nf90_inq_varid(ncidp, 'ftsol4', varidfts4)
210    PRINT *, 'ncidp,varidfts4', ncidp, varidfts4
211
212    rcode = nf90_inq_varid(ncidp, 'psrf1', varidpsr1)
213    PRINT *, 'ncidp,varidpsr1', ncidp, varidpsr1
214
215    rcode = nf90_inq_varid(ncidp, 'psrf2', varidpsr2)
216    PRINT *, 'ncidp,varidpsr2', ncidp, varidpsr2
217
218    rcode = nf90_inq_varid(ncidp, 'psrf3', varidpsr3)
219    PRINT *, 'ncidp,varidpsr3', ncidp, varidpsr3
220
221    rcode = nf90_inq_varid(ncidp, 'psrf4', varidpsr4)
222    PRINT *, 'ncidp,varidpsr4', ncidp, varidpsr4
223
224    ! ID pour les dimensions
225
226    status = nf90_inq_dimid(ncidp, 'y', yid)
227    status = nf90_inq_dimid(ncidp, 'x', xid)
228    status = nf90_inq_dimid(ncidp, 'sig_s', zid)
229    status = nf90_inq_dimid(ncidp, 'time_counter', tid)
230
231    ! lecture des dimensions
232
233    status = nf90_inquire_dimension(ncidp, yid, namedim, ncjm)
234    status = nf90_inquire_dimension(ncidp, xid, namedim, ncim)
235    status = nf90_inquire_dimension(ncidp, zid, namedim, nckev)
236    status = nf90_inquire_dimension(ncidp, tid, namedim, ncrec)
237
238    zrec = ncrec
239    zkev = nckev
240    zim = ncim
241    zjm = ncjm
242
243    zkon = zim*(zjm-2) + 2
244
245    WRITE (*, *) 'read_pstoke : zrec = ', zrec
246    WRITE (*, *) 'read_pstoke : kev = ', zkev
247    WRITE (*, *) 'read_pstoke : zim = ', zim
248    WRITE (*, *) 'read_pstoke : zjm = ', zjm
249    WRITE (*, *) 'read_pstoke : kon = ', zkon
250
251    ! niveaux de pression
252
253    status = nf90_get_var(ncidp, varidpl, pl, [1], [kev])
254
255    ! lecture de aire et phis
256
257    start(1) = 1
258    start(2) = 1
259    start(3) = 1
260    start(4) = 0
261
262    count(1) = zim
263    count(2) = zjm
264    count(3) = 1
265    count(4) = 0
266
267
268    ! **** Geopotentiel au sol ***************************************
269    ! phis
270    status = nf90_get_var(ncidp, varidps, phisfi2, start, count)
271    CALL gr_ecrit_fi(1, kon, nbp_lon, nbp_lat, phisfi2, phisfi)
272
273    ! **** Aires des mails aux sol ************************************
274    ! aire
275    status = nf90_get_var(ncidp, varidai, airefi2, start, count)
276    CALL gr_ecrit_fi(1, kon, nbp_lon, nbp_lat, airefi2, airefi)
277  ELSE
278
279    PRINT *, 'ok1'
280
281    ! ---------------------
282    ! lecture des champs
283    ! ---------------------
284
285    PRINT *, 'WARNING!!! Il n y a pas de test de coherence'
286    PRINT *, 'sur le nombre de niveaux verticaux dans le fichier nc'
287
288    start(1) = 1
289    start(2) = 1
290    start(3) = 1
291    start(4) = irec
292
293    count(1) = zim
294    count(2) = zjm
295    count(3) = kev
296    count(4) = 1
297
298    ! **** Temperature ********************************************
299    ! A FAIRE : Es-ce necessaire ?
300
301    ! abder t
302    status = nf90_get_var(ncidp, varidt, t2, start, count)
303    CALL gr_ecrit_fi(kev, kon, nbp_lon, nbp_lat, t2, t)
304
305    ! **** Flux pour la convection (Tiedtk)
306    ! ********************************************
307    ! mfu
308    status = nf90_get_var(ncidp, varidmfu, mfu2, start, count)
309    CALL gr_ecrit_fi(kev, kon, nbp_lon, nbp_lat, mfu2, mfu)
310
311    ! mfd
312    status = nf90_get_var(ncidp, varidmfd, mfd2, start, count)
313    CALL gr_ecrit_fi(kev, kon, nbp_lon, nbp_lat, mfd2, mfd)
314
315    ! en_u
316    status = nf90_get_var(ncidp, varidenu, en_u2, start, count)
317    CALL gr_ecrit_fi(kev, kon, nbp_lon, nbp_lat, en_u2, en_u)
318
319    ! de_u
320    status = nf90_get_var(ncidp, variddeu, de_u2, start, count)
321    CALL gr_ecrit_fi(kev, kon, nbp_lon, nbp_lat, de_u2, de_u)
322
323    ! en_d
324    status = nf90_get_var(ncidp, varidend, en_d2, start, count)
325    CALL gr_ecrit_fi(kev, kon, nbp_lon, nbp_lat, en_d2, en_d)
326
327    ! de_d
328    status = nf90_get_var(ncidp, varidded, de_d2, start, count)
329    CALL gr_ecrit_fi(kev, kon, nbp_lon, nbp_lat, de_d2, de_d)
330
331    ! **** Coefficient de mellange turbulent
332    ! *******************************************
333    ! coefh
334    PRINT *, 'LECTURE de coefh a irec =', irec
335    status = nf90_get_var(ncidp, varidch, coefh2, start, count)
336    CALL gr_ecrit_fi(kev, kon, nbp_lon, nbp_lat, coefh2, coefh)
337    ! call dump2d(iip1,jjp1,coefh2(1,2),'COEFH2READ   ')
338    ! call dump2d(iim ,jjm ,coefh (2,2),'COEFH2READ   ')
339
340    ! **** Flux ascendants et entrant dans le thermique
341    ! **********************************
342    ! Thermiques
343    PRINT *, 'LECTURE de fm_therm a irec =', irec
344    status = nf90_get_var(ncidp, varidfmth, fm_therm2, start, count)
345    CALL gr_ecrit_fi(kev, kon, nbp_lon, nbp_lat, fm_therm2, fm_therm)
346    PRINT *, 'LECTURE de en_therm a irec =', irec
347    status = nf90_get_var(ncidp, varidenth, en_therm2, start, count)
348    CALL gr_ecrit_fi(kev, kon, nbp_lon, nbp_lat, en_therm2, en_therm)
349
350    ! **** Coefficients de lessivage
351    ! *******************************************
352    ! frac_impa
353    status = nf90_get_var(ncidp, varidfi, frac_impa2, start, count)
354    CALL gr_ecrit_fi(kev, kon, nbp_lon, nbp_lat, frac_impa2, frac_impa)
355
356    ! frac_nucl
357
358    status = nf90_get_var(ncidp, varidfn, frac_nucl2, start, count)
359    CALL gr_ecrit_fi(kev, kon, nbp_lon, nbp_lat, frac_nucl2, frac_nucl)
360
361    ! **** Vents aux sol ********************************************
362
363    start(3) = irec
364    start(4) = 0
365    count(3) = 1
366    count(4) = 0
367
368    ! pyu1
369    PRINT *, 'LECTURE de yu1 a irec =', irec
370    status = nf90_get_var(ncidp, varidyu1, pyu12, start, count)
371    CALL gr_ecrit_fi(1, kon, nbp_lon, nbp_lat, pyu12, pyu1)
372
373    ! pyv1
374    PRINT *, 'LECTURE de yv1 a irec =', irec
375    status = nf90_get_var(ncidp, varidyv1, pyv12, start, count)
376    CALL gr_ecrit_fi(1, kon, nbp_lon, nbp_lat, pyv12, pyv1)
377
378    ! **** Temerature au sol ********************************************
379    ! ftsol1
380    PRINT *, 'LECTURE de ftsol1 a irec =', irec
381    status = nf90_get_var(ncidp, varidfts1, ftsol12, start, count)
382    CALL gr_ecrit_fi(1, kon, nbp_lon, nbp_lat, ftsol12, ftsol1)
383
384    ! ftsol2
385    PRINT *, 'LECTURE de ftsol2 a irec =', irec
386    status = nf90_get_var(ncidp, varidfts2, ftsol22, start, count)
387    CALL gr_ecrit_fi(1, kon, nbp_lon, nbp_lat, ftsol22, ftsol2)
388
389    ! ftsol3
390    PRINT *, 'LECTURE de ftsol3 a irec =', irec
391    status = nf90_get_var(ncidp, varidfts3, ftsol32, start, count)
392    CALL gr_ecrit_fi(1, kon, nbp_lon, nbp_lat, ftsol32, ftsol3)
393
394    ! ftsol4
395    status = nf90_get_var(ncidp, varidfts4, ftsol42, start, count)
396    CALL gr_ecrit_fi(1, kon, nbp_lon, nbp_lat, ftsol42, ftsol4)
397
398    ! **** Nature sol ********************************************
399    ! psrf1
400    status = nf90_get_var(ncidp, varidpsr1, psrf12, start, count)
401    ! call dump2d(iip1-1,jjm+1,psrf12,'PSRF1NC')
402    CALL gr_ecrit_fi(1, kon, nbp_lon, nbp_lat, psrf12, psrf1)
403
404    ! psrf2
405    status = nf90_get_var(ncidp, varidpsr2, psrf22, start, count)
406    ! call dump2d(iip1-1,jjm+1,psrf22,'PSRF2NC')
407    CALL gr_ecrit_fi(1, kon, nbp_lon, nbp_lat, psrf22, psrf2)
408
409    ! psrf3
410    status = nf90_get_var(ncidp, varidpsr3, psrf32, start, count)
411    CALL gr_ecrit_fi(1, kon, nbp_lon, nbp_lat, psrf32, psrf3)
412
413    ! psrf4
414    status = nf90_get_var(ncidp, varidpsr4, psrf42, start, count)
415    CALL gr_ecrit_fi(1, kon, nbp_lon, nbp_lat, psrf42, psrf4)
416
417    DO i = 1, kon
418
419      psrf(i, 1) = psrf1(i)
420      psrf(i, 2) = psrf2(i)
421      psrf(i, 3) = psrf3(i)
422      ! test abderr
423      ! print*,'Dans read_pstoke psrf3 =',psrf3(i),i
424      psrf(i, 4) = psrf4(i)
425
426      ftsol(i, 1) = ftsol1(i)
427      ftsol(i, 2) = ftsol2(i)
428      ftsol(i, 3) = ftsol3(i)
429      ftsol(i, 4) = ftsol4(i)
430
431    END DO
432
433  END IF
434
435  RETURN
436
437END SUBROUTINE read_pstoke0
438
Note: See TracBrowser for help on using the repository browser.