source: LMDZ6/trunk/libf/phylmd/read_pstoke.F90 @ 5202

Last change on this file since 5202 was 5084, checked in by Laurent Fairhead, 12 months ago

Reverting to r4065. Updating fortran standard broke too much stuff. Will do it by smaller chunks
AB, LF

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