source: LMDZ4/trunk/libf/phytherm/read_pstoke.F @ 1098

Last change on this file since 1098 was 814, checked in by Laurent Fairhead, 17 years ago

Rajout de la physique utilisant les thermiques FH
LF

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 14.2 KB
Line 
1!
2! $Header$
3!
4c
5c
6        subroutine read_pstoke(irec,
7     .   zrec,zklono,zklevo,airefi,phisfi,
8     .   t,mfu,mfd,en_u,de_u,en_d,de_d,coefh,
9     .   fm_therm,en_therm,
10     .   frac_impa,frac_nucl,pyu1,pyv1,ftsol,psrf)
11
12C******************************************************************************
13C  Frederic HOURDIN, Abderrahmane IDELKADI
14C Lecture des parametres physique stockes online necessaires pour
15C recalculer offline le transport de traceurs sur une grille 2x plus fine que
16C celle online
17C A FAIRE : une seule routine au lieu de 2 (lectflux, redecoupe)!
18C******************************************************************************
19
20
21       IMPLICIT NONE
22
23#include "netcdf.inc"
24#include "dimensions.h"
25#include "paramet.h"
26#include "comconst.h"
27#include "comgeom.h"
28#include "temps.h"
29#include "ener.h"
30#include "logic.h"
31#include "description.h"
32#include "serre.h"
33#include "indicesol.h"
34#include "control.h"
35#include "dimphy.h"
36       
37          integer klono,klevo,imo,jmo
38          parameter (imo=iim/2,jmo=(jjm+1)/2)
39          parameter(klono=(jmo-1)*imo+2,klevo=llm)
40          REAL phisfi(klono)
41          REAL phisfi2(imo,jmo+1),airefi2(imo,jmo+1)
42
43          REAL mfu(klono,klevo), mfd(klono,klevo)
44          REAL en_u(klono,klevo), de_u(klono,klevo)
45          REAL en_d(klono,klevo), de_d(klono,klevo)
46          REAL coefh(klono,klevo)
47           REAL fm_therm(klono,klevo),en_therm(klono,klevo)
48
49          REAL mfu2(imo,jmo+1,klevo), mfd2(imo,jmo+1,klevo)
50          REAL en_u2(imo,jmo+1,klevo), de_u2(imo,jmo+1,klevo)
51          REAL en_d2(imo,jmo+1,klevo), de_d2(imo,jmo+1,klevo)
52          REAL coefh2(imo,jmo+1,klevo)
53           REAL fm_therm2(imo,jmo+1,klevo)
54           REAL en_therm2(imo,jmo+1,klevo)
55
56          REAL pl(klevo)
57          integer irec
58          integer xid,yid,zid,tid
59          real zrec,zklono,zklevo,zim,zjm
60          integer ncrec,ncklono,ncklevo,ncim,ncjm
61
62          real airefi(klono)
63          character*20 namedim
64
65c  !! attention !!
66c attention il y a aussi le pb de def klono
67c dim de phis??
68         
69         
70          REAL frac_impa(klono,klevo), frac_nucl(klono,klevo)
71          REAL frac_impa2(imo,jmo+1,klevo),
72     .     frac_nucl2(imo,jmo+1,klevo)
73          REAL pyu1(klono), pyv1(klono)
74          REAL pyu12(imo,jmo+1), pyv12(imo,jmo+1)
75          REAL ftsol(klono,nbsrf)
76          REAL psrf(klono,nbsrf)
77          REAL ftsol1(klono),ftsol2(klono),ftsol3(klono),ftsol4(klono)
78          REAL psrf1(klono),psrf2(klono),psrf3(klono),psrf4(klono)
79          REAL ftsol12(imo,jmo+1),ftsol22(imo,jmo+1),
80     .     ftsol32(imo,jmo+1),
81     .     ftsol42(imo,jmo+1)
82          REAL psrf12(imo,jmo+1),psrf22(imo,jmo+1),psrf32(imo,jmo+1),
83     .     psrf42(imo,jmo+1)
84                REAL t(klono,klevo)
85                REAL t2(imo,jmo+1,klevo)       
86          integer ncidp
87          save ncidp
88                integer varidt
89          integer varidmfu, varidmfd, varidps, varidenu, variddeu       
90          integer varidend,varidded,varidch,varidfi,varidfn
91           integer varidfmth,varidenth
92          integer varidyu1,varidyv1,varidpl,varidai,varididvt
93          integer varidfts1,varidfts2,varidfts3,varidfts4
94          integer varidpsr1,varidpsr2,varidpsr3,varidpsr4
95          save varidmfu, varidmfd, varidps, varidenu, variddeu
96          save varidend,varidded,varidch,varidfi,varidfn
97           save varidfmth,varidenth
98          save varidyu1,varidyv1,varidpl,varidai,varididvt
99          save varidfts1,varidfts2,varidfts3,varidfts4
100          save varidpsr1,varidpsr2,varidpsr3,varidpsr4
101                save varidt
102
103          integer l, i
104          integer start(4),count(4),status
105          real rcode
106          logical first
107          save first
108          data first/.true./
109
110
111
112c ---------------------------------------------
113c   Initialisation de la lecture des fichiers
114c ---------------------------------------------
115
116      if (irec .eq. 0) then
117
118            ncidp=NCOPN('phystoke.nc',NCNOWRIT,rcode)
119
120            varidps=NCVID(ncidp,'phis',rcode)
121            print*,'ncidp,varidps',ncidp,varidps
122
123            varidpl=NCVID(ncidp,'sig_s',rcode)
124            print*,'ncidp,varidpl',ncidp,varidpl
125
126            varidai=NCVID(ncidp,'aire',rcode)
127            print*,'ncidp,varidai',ncidp,varidai
128
129c A FAIRE: Es-il necessaire de stocke t?
130                varidt=NCVID(ncidp,'t',rcode)
131                print*,'ncidp,varidt',ncidp,varidt
132
133            varidmfu=NCVID(ncidp,'mfu',rcode)
134            print*,'ncidp,varidmfu',ncidp,varidmfu
135
136            varidmfd=NCVID(ncidp,'mfd',rcode)
137            print*,'ncidp,varidmfd',ncidp,varidmfd
138
139            varidenu=NCVID(ncidp,'en_u',rcode)
140            print*,'ncidp,varidenu',ncidp,varidenu
141
142            variddeu=NCVID(ncidp,'de_u',rcode)
143            print*,'ncidp,variddeu',ncidp,variddeu
144
145            varidend=NCVID(ncidp,'en_d',rcode)
146            print*,'ncidp,varidend',ncidp,varidend
147       
148            varidded=NCVID(ncidp,'de_d',rcode)
149            print*,'ncidp,varidded',ncidp,varidded
150       
151            varidch=NCVID(ncidp,'coefh',rcode)
152            print*,'ncidp,varidch',ncidp,varidch
153       
154c abder (pour thermiques)
155             varidfmth=NCVID(ncidp,'fm_th',rcode)
156             print*,'ncidp,varidfmth',ncidp,varidfmth
157
158             varidenth=NCVID(ncidp,'en_th',rcode)
159             print*,'ncidp,varidenth',ncidp,varidenth
160
161            varidfi=NCVID(ncidp,'frac_impa',rcode)
162            print*,'ncidp,varidfi',ncidp,varidfi
163       
164            varidfn=NCVID(ncidp,'frac_nucl',rcode)
165            print*,'ncidp,varidfn',ncidp,varidfn
166       
167            varidyu1=NCVID(ncidp,'pyu1',rcode)
168            print*,'ncidp,varidyu1',ncidp,varidyu1
169       
170            varidyv1=NCVID(ncidp,'pyv1',rcode)
171            print*,'ncidp,varidyv1',ncidp,varidyv1
172       
173            varidfts1=NCVID(ncidp,'ftsol1',rcode)
174            print*,'ncidp,varidfts1',ncidp,varidfts1
175       
176            varidfts2=NCVID(ncidp,'ftsol2',rcode)
177            print*,'ncidp,varidfts2',ncidp,varidfts2
178         
179            varidfts3=NCVID(ncidp,'ftsol3',rcode)
180            print*,'ncidp,varidfts3',ncidp,varidfts3
181 
182            varidfts4=NCVID(ncidp,'ftsol4',rcode)
183            print*,'ncidp,varidfts4',ncidp,varidfts4
184       
185            varidpsr1=NCVID(ncidp,'psrf1',rcode)
186            print*,'ncidp,varidpsr1',ncidp,varidpsr1
187       
188            varidpsr2=NCVID(ncidp,'psrf2',rcode)
189            print*,'ncidp,varidpsr2',ncidp,varidpsr2
190       
191            varidpsr3=NCVID(ncidp,'psrf3',rcode)
192            print*,'ncidp,varidpsr3',ncidp,varidpsr3
193
194            varidpsr4=NCVID(ncidp,'psrf4',rcode)
195            print*,'ncidp,varidpsr4',ncidp,varidpsr4
196       
197c ID pour les dimensions
198
199            status = nf_inq_dimid(ncidp,'y',yid)
200            status = nf_inq_dimid(ncidp,'x',xid)
201            status = nf_inq_dimid(ncidp,'sig_s',zid)
202            status = nf_inq_dimid(ncidp,'time_counter',tid)
203
204c lecture des dimensions
205
206            status = nf_inq_dim(ncidp,yid,namedim,ncjm)
207            status = nf_inq_dim(ncidp,xid,namedim,ncim)
208            status = nf_inq_dim(ncidp,zid,namedim,ncklevo)
209            status = nf_inq_dim(ncidp,tid,namedim,ncrec)
210       
211            zrec=ncrec
212            zklevo=ncklevo
213            zim=ncim
214            zjm=ncjm
215       
216            zklono=zim*(zjm-2)+2
217       
218            write(*,*) 'read_pstoke : zrec = ', zrec
219            write(*,*) 'read_pstoke : zklevo = ', zklevo
220            write(*,*) 'read_pstoke : zim = ', zim
221            write(*,*) 'read_pstoke : zjm = ', zjm
222            write(*,*) 'read_pstoke : zklono = ', zklono
223
224c niveaux de pression
225#ifdef NC_DOUBLE
226      status=NF_GET_VARA_DOUBLE(ncidp,varidpl,1,zklevo,pl)
227#else
228      status=NF_GET_VARA_REAL(ncidp,varidpl,1,zklevo,pl)
229#endif
230
231c lecture de aire et phis
232       
233      start(1)=1
234      start(2)=1
235      start(3)=1
236      start(4)=0
237
238      count(1)=zim
239      count(2)=zjm
240      count(3)=1
241      count(4)=0
242
243c phis
244#ifdef NC_DOUBLE
245      status=NF_GET_VARA_DOUBLE(ncidp,varidps,start,count,phisfi2)
246#else
247      status=NF_GET_VARA_REAL(ncidp,varidps,start,count,phisfi2)
248#endif
249      call gr_ecrit_fi(1,klono,imo,jmo+1,phisfi2,phisfi)
250
251c aire
252#ifdef NC_DOUBLE
253      status=NF_GET_VARA_DOUBLE(ncidp,varidai,start,count,airefi2)
254#else
255      status=NF_GET_VARA_REAL(ncidp,varidai,start,count,airefi2)
256#endif
257       call gr_ecrit_fi(1,klono,imo,jmo+1,airefi2,airefi)
258      else
259
260      print*,'ok1'
261
262c ---------------------
263c   lecture des champs
264c ---------------------
265       
266        print*,'WARNING!!! Il n y a pas de test de coherence'
267        print*,'sur le nombre de niveaux verticaux dans le fichier nc'
268
269      start(1)=1
270      start(2)=1
271      start(3)=1
272      start(4)=irec
273
274      count(1)=zim
275      count(2)=zjm
276      count(3)=zklevo
277      count(4)=1
278
279
280C *** Lessivage******************************************************
281c frac_impa
282#ifdef NC_DOUBLE
283      status=NF_GET_VARA_DOUBLE(ncidp,varidfi,start,count,frac_impa2)
284#else
285      status=NF_GET_VARA_REAL(ncidp,varidfi,start,count,frac_impa2)
286#endif
287      call gr_ecrit_fi(klevo,klono,imo,jmo+1,frac_impa2,frac_impa)
288
289c frac_nucl
290#ifdef NC_DOUBLE
291      status=NF_GET_VARA_DOUBLE(ncidp,varidfn,start,count,frac_nucl2)
292#else
293      status=NF_GET_VARA_REAL(ncidp,varidfn,start,count,frac_nucl2)
294#endif
295      call gr_ecrit_fi(klevo,klono,imo,jmo+1,frac_nucl2,frac_nucl)
296
297C*** Temperature ******************************************************
298c abder t
299#ifdef NC_DOUBLE
300      status=NF_GET_VARA_DOUBLE(ncidp,varidt,start,count,t2)
301#else
302      status=NF_GET_VARA_REAL(ncidp,varidt,start,count,t2)
303#endif
304      call gr_ecrit_fi(klevo,klono,imo,jmo+1,t2,t)
305
306C*** Flux pour le calcul de la convection TIEDTK ***********************
307c mfu
308#ifdef NC_DOUBLE
309      status=NF_GET_VARA_DOUBLE(ncidp,varidmfu,start,count,mfu2)
310#else
311      status=NF_GET_VARA_REAL(ncidp,varidmfu,start,count,mfu2)
312#endif
313      call gr_ecrit_fi(klevo,klono,imo,jmo+1,mfu2,mfu)
314
315c mfd
316#ifdef NC_DOUBLE
317      status=NF_GET_VARA_DOUBLE(ncidp,varidmfd,start,count,mfd2)
318#else
319      status=NF_GET_VARA_REAL(ncidp,varidmfd,start,count,mfd2)
320#endif
321      call gr_ecrit_fi(klevo,klono,imo,jmo+1,mfd2,mfd)
322
323c en_u
324#ifdef NC_DOUBLE
325      status=NF_GET_VARA_DOUBLE(ncidp,varidenu,start,count,en_u2)
326#else
327      status=NF_GET_VARA_REAL(ncidp,varidenu,start,count,en_u2)
328#endif
329      call gr_ecrit_fi(klevo,klono,imo,jmo+1,en_u2,en_u)
330
331c de_u
332#ifdef NC_DOUBLE
333      status=NF_GET_VARA_DOUBLE(ncidp,variddeu,start,count,de_u2)
334#else
335      status=NF_GET_VARA_REAL(ncidp,variddeu,start,count,de_u2)
336#endif
337      call gr_ecrit_fi(klevo,klono,imo,jmo+1,de_u2,de_u)
338
339c en_d
340#ifdef NC_DOUBLE
341      status=NF_GET_VARA_DOUBLE(ncidp,varidend,start,count,en_d2)
342#else
343      status=NF_GET_VARA_REAL(ncidp,varidend,start,count,en_d2)
344#endif
345      call gr_ecrit_fi(klevo,klono,imo,jmo+1,en_d2,en_d)
346
347c de_d
348#ifdef NC_DOUBLE
349      status=NF_GET_VARA_DOUBLE(ncidp,varidded,start,count,de_d2)
350#else
351      status=NF_GET_VARA_REAL(ncidp,varidded,start,count,de_d2)
352#endif
353      call gr_ecrit_fi(klevo,klono,imo,jmo+1,de_d2,de_d)
354
355C **** Coeffecient du mellange turbulent**********************************
356c coefh
357#ifdef NC_DOUBLE
358      status=NF_GET_VARA_DOUBLE(ncidp,varidch,start,count,coefh2)
359#else
360      status=NF_GET_VARA_REAL(ncidp,varidch,start,count,coefh2)
361#endif
362       call gr_ecrit_fi(klevo,klono,imo,jmo+1,coefh2,coefh)
363
364C*** Flux ascendant et entrant pour les Thermiques************************
365cabder thermiques
366#ifdef NC_DOUBLE
367      status=NF_GET_VARA_DOUBLE(ncidp,varidfmth,start,count,fm_therm2)
368#else
369      status=NF_GET_VARA_REAL(ncidp,varidfmth,start,count,fm_therm2)
370#endif
371      call gr_ecrit_fi(klevo,klono,imo,jmo+1,fm_therm2,fm_therm)
372
373#ifdef NC_DOUBLE
374      status=NF_GET_VARA_DOUBLE(ncidp,varidenth,start,count,en_therm2)
375#else
376      status=NF_GET_VARA_REAL(ncidp,varidenth,start,count,en_therm2)
377#endif
378      call gr_ecrit_fi(klevo,klono,imo,jmo+1,en_therm2,en_therm)
379
380C*** Vitesses aux sol ******************************************************
381      start(3)=irec
382      start(4)=0
383      count(3)=1
384      count(4)=0
385c pyu1
386#ifdef NC_DOUBLE
387      status=NF_GET_VARA_DOUBLE(ncidp,varidyu1,start,count,pyu12)
388#else
389      status=NF_GET_VARA_REAL(ncidp,varidyu1,start,count,pyu12)
390#endif
391      call gr_ecrit_fi(1,klono,imo,jmo+1,pyu12,pyu1)
392
393c pyv1
394#ifdef NC_DOUBLE
395      status=NF_GET_VARA_DOUBLE(ncidp,varidyv1,start,count,pyv12)
396#else
397      status=NF_GET_VARA_REAL(ncidp,varidyv1,start,count,pyv12)
398#endif
399      call gr_ecrit_fi(1,klono,imo,jmo+1,pyv12,pyv1)
400
401C*** Temperature au sol ********************************************
402c ftsol1
403#ifdef NC_DOUBLE
404      status=NF_GET_VARA_DOUBLE(ncidp,varidfts1,start,count,ftsol12)
405#else
406      status=NF_GET_VARA_REAL(ncidp,varidfts1,start,count,ftsol12)
407#endif
408       call gr_ecrit_fi(1,klono,imo,jmo+1,ftsol12,ftsol1)
409
410c ftsol2
411#ifdef NC_DOUBLE
412      status=NF_GET_VARA_DOUBLE(ncidp,varidfts2,start,count,ftsol22)
413#else
414      status=NF_GET_VARA_REAL(ncidp,varidfts2,start,count,ftsol22)
415#endif
416      call gr_ecrit_fi(1,klono,imo,jmo+1,ftsol22,ftsol2)
417
418c ftsol3
419#ifdef NC_DOUBLE
420      status=NF_GET_VARA_DOUBLE(ncidp,varidfts3,start,count,ftsol32)
421#else
422      status=NF_GET_VARA_REAL(ncidp,varidfts3,start,count,ftsol32)
423#endif
424      call gr_ecrit_fi(1,klono,imo,jmo+1,ftsol32,ftsol3)
425
426c ftsol4
427#ifdef NC_DOUBLE
428      status=NF_GET_VARA_DOUBLE(ncidp,varidfts4,start,count,ftsol42)
429#else
430      status=NF_GET_VARA_REAL(ncidp,varidfts4,start,count,ftsol42)
431#endif
432      call gr_ecrit_fi(1,klono,imo,jmo+1,ftsol42,ftsol4)
433
434C*** Nature du sol **************************************************
435c psrf1
436#ifdef NC_DOUBLE
437      status=NF_GET_VARA_DOUBLE(ncidp,varidpsr1,start,count,psrf12)
438#else
439      status=NF_GET_VARA_REAL(ncidp,varidpsr1,start,count,psrf12)
440#endif
441      call gr_ecrit_fi(1,klono,imo,jmo+1,psrf12,psrf1)
442
443c psrf2
444#ifdef NC_DOUBLE
445      status=NF_GET_VARA_DOUBLE(ncidp,varidpsr2,start,count,psrf22)
446#else
447      status=NF_GET_VARA_REAL(ncidp,varidpsr2,start,count,psrf22)
448#endif
449      call gr_ecrit_fi(1,klono,imo,jmo+1,psrf22,psrf2)
450
451c psrf3
452#ifdef NC_DOUBLE
453      status=NF_GET_VARA_DOUBLE(ncidp,varidpsr3,start,count,psrf32)
454#else
455      status=NF_GET_VARA_REAL(ncidp,varidpsr3,start,count,psrf32)
456#endif
457      call gr_ecrit_fi(1,klono,imo,jmo+1,psrf32,psrf3)
458
459c psrf4
460#ifdef NC_DOUBLE
461      status=NF_GET_VARA_DOUBLE(ncidp,varidpsr4,start,count,psrf42)
462#else
463      status=NF_GET_VARA_REAL(ncidp,varidpsr4,start,count,psrf42)
464#endif
465      call gr_ecrit_fi(1,klono,imo,jmo+1,psrf42,psrf4)
466       
467          do i = 1,klono
468       
469        psrf(i,1) = psrf1(i)
470        psrf(i,2) = psrf2(i)
471        psrf(i,3) = psrf3(i)
472        psrf(i,4) = psrf4(i)
473 
474        ftsol(i,1) = ftsol1(i)
475        ftsol(i,2) = ftsol2(i)
476        ftsol(i,3) = ftsol3(i)
477        ftsol(i,4) = ftsol4(i)
478       
479          enddo
480       
481        endif
482       
483        return
484       
485        end
486
Note: See TracBrowser for help on using the repository browser.