source: LMDZ6/trunk/libf/phylmd/read_pstoke.f90 @ 5408

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