source: LMDZ.3.3/trunk/libf/dyn3d/etat0_netcdf.F

Last change on this file was 259, checked in by lmdz, 23 years ago

Nouveaux programmes pour la creation des etats initiaux et des conditions aux limites. Levan
LF
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 11.7 KB
Line 
1C
2C $Header$
3C
4      SUBROUTINE etat0_netcdf ( interbar )
5   
6      USE startvar
7      USE ioipsl
8      !
9      IMPLICIT NONE
10      !
11#include "dimensions.h"
12#include "paramet.h"
13      !
14      !
15      INTEGER, PARAMETER :: KIDIA=1, KFDIA=iim*(jjm-1)+2,
16     .KLON=KFDIA-KIDIA+1,KLEV=llm
17      !
18#include "comgeom2.h"
19#include "comvert.h"
20#include "comconst.h"
21      !
22      LOGICAL interbar
23      REAL :: latfi(klon), lonfi(klon)
24      REAL :: orog(iip1,jjp1), rugo(iip1,jjp1), masque(iip1,jjp1),
25     . psol(iip1, jjp1), phis(iip1, jjp1)
26      REAL :: p3d(iip1, jjp1, llm+1)
27      REAL :: uvent(iip1, jjp1, llm)
28      REAL :: vvent(iip1, jjm, llm)
29      REAL :: t3d(iip1, jjp1, llm), tpot(iip1, jjp1, llm)
30      REAL :: q3d(iip1, jjp1, llm,nqmx), qsat(iip1, jjp1, llm)
31      REAL :: tsol(klon), qsol(klon), sn(klon), radsol(klon)
32      REAL :: deltat(klon), rugmer(klon), agesno(klon)
33      REAL :: zmea(iip1*jjp1), zstd(iip1*jjp1)
34      REAL :: zsig(iip1*jjp1), zgam(iip1*jjp1), zthe(iip1*jjp1)
35      REAL :: zpic(iip1*jjp1), zval(iip1*jjp1), rugsrel(iip1*jjp1)
36      REAL :: qd(iip1, jjp1, llm)
37      !
38      CHARACTER*80 :: varname
39      !
40      INTEGER :: i,j, ig, l,ii1,ii2
41      REAL :: xpi
42      !
43      REAL :: alpha(iip1,jjp1,llm),beta(iip1,jjp1,llm)
44      REAL :: pk(iip1,jjp1,llm), pls(iip1,jjp1,llm), pks(ip1jmp1)
45      REAL :: workvar(iip1,jjp1,llm)
46      !
47      REAL ::  prefkap, unskap
48      !
49      REAL :: q_sat
50      EXTERNAL q_sat
51      real :: time_step,t_ops,t_wrt
52
53#include "comdissnew.h"
54#include "control.h"
55#include "serre.h"
56#include "clesph0.h"
57
58      INTEGER  ::        longcles
59      PARAMETER      ( longcles  = 20 )
60      REAL :: clesphy0 ( longcles       )
61      REAL :: p(iip1,jjp1,llm)
62      INTEGER :: itau, iday
63      REAL :: masse(iip1,jjp1,llm)
64      REAL :: xpn,xps,xppn(iim),xpps(iim)
65      real :: time
66      REAL :: phi(ip1jmp1,llm)
67      REAL :: pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
68      REAL :: w(ip1jmp1,llm)
69      REAL ::phystep,co2_ppm,solaire
70      INTEGER :: radpas
71       real zrel(iip1*jjp1),chmin,chmax
72
73      CHARACTER*80 :: visu_file
74      INTEGER :: visuid
75      !
76      !   Constantes
77      !
78      pi     = 4. * ATAN(1.)
79      rad    = 6371229.
80      omeg   = 4.* ASIN(1.)/(24.*3600.)
81      g      = 9.8
82      daysec = 86400.
83      kappa  = 0.2857143
84      cpp    = 1004.70885
85      !
86      preff     = 101325.
87      unskap = 1./kappa
88      !
89      jmp1    = jjm + 1
90      !
91      !    Construct a grid
92      !
93
94      CALL defrun_new(99,.TRUE.,clesphy0)
95
96      dtvr   = daysec/FLOAT(day_step)
97      print*,'dtvr',dtvr
98
99      CALL inicons0()
100      CALL inigeom()
101      !
102      CALL inifilr()
103      !
104      latfi(1) = ASIN(1.0)
105      DO j = 2, jjm
106        DO i = 1, iim
107          latfi((j-2)*iim+1+i)=  rlatu(j)
108        ENDDO
109      ENDDO
110      latfi(klon) = - ASIN(1.0)
111      !
112      lonfi(1) = 0.0
113      DO j = 2, jjm
114        DO i = 1, iim
115          lonfi((j-2)*iim+1+i) =  rlonv(i)
116        ENDDO
117      ENDDO
118      lonfi(klon) = 0.0
119      !
120      xpi = 2.0 * ASIN(1.0)
121      DO ig = 1, klon
122        latfi(ig) = latfi(ig) * 180.0 / xpi
123        lonfi(ig) = lonfi(ig) * 180.0 / xpi
124      ENDDO
125      !
126      varname = 'relief'
127      ! This line needs to be replaced by a call to restget to get the values in the restart file
128      orog(:,:) = 0.0
129       CALL startget(varname, iip1, jjp1, rlonv, rlatu, orog, 0.0 ,
130     , jjm ,rlonu,rlatv , interbar )
131      !
132      WRITE(*,*) 'OUT OF GET VARIABLE : Relief'
133      WRITE(*,'(49I1)') INT(orog(:,:))
134      !
135      varname = 'rugosite'
136      ! This line needs to be replaced by a call to restget to get the values in the restart file
137      rugo(:,:) = 0.0
138       CALL startget(varname, iip1, jjp1, rlonv, rlatu, rugo, 0.0 ,
139     , jjm, rlonu,rlatv , interbar )
140      !
141      WRITE(*,*) 'OUT OF GET VARIABLE : Rugosite'
142      WRITE(*,'(49I1)') INT(rugo(:,:)*10)
143      !
144      varname = 'masque'
145      ! This line needs to be replaced by a call to restget to get the values in the restart file
146      masque(:,:) = 0.0
147       CALL startget(varname, iip1, jjp1, rlonv, rlatu, masque, 0.0,
148     , jjm ,rlonu,rlatv , interbar )
149      !
150c      WRITE(*,*) 'MASQUE construit : Masque'
151c      WRITE(*,'(49I1)') INT(masque(:,:))
152      !
153      !
154      varname = 'psol'
155      psol(:,:) = 0.0
156      CALL startget(varname, iip1, jjp1, rlonv, rlatu, psol, 0.0 ,
157     , jjm ,rlonu,rlatv , interbar )
158      !
159      !  Compute here the pressure on the intermediate levels. One would expect that this is available in the GCM
160      !  anyway.
161      !
162      WRITE(*,*) 'PSOL :', psol(10,20)
163      WRITE(*,*) ap(:), bp(:)
164      CALL pression(ip1jmp1, ap, bp, psol, p3d)
165      WRITE(*,*) 'P3D :', p3d(10,20,:)
166      CALL exner_hyb(ip1jmp1, psol, p3d, alpha, beta, pks, pk, workvar)
167      WRITE(*,*) 'PK:', pk(10,20,:)
168      !
169      !
170      !
171      prefkap =  preff  ** kappa
172      WRITE(*,*) 'unskap, cpp,  preff :', unskap, cpp,  preff
173      DO l = 1, llm
174        DO j=1,jjp1
175          DO i =1, iip1
176            pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap
177           ENDDO
178        ENDDO
179      ENDDO
180      !
181      WRITE(*,*) 'PLS :', pls(10,20,:)
182      !
183      varname = 'surfgeo'
184      phis(:,:) = 0.0
185      CALL startget(varname, iip1, jjp1, rlonv, rlatu, phis, 0.0 ,
186     , jjm ,rlonu,rlatv, interbar )
187      !
188      varname = 'u'
189      uvent(:,:,:) = 0.0
190      CALL startget(varname, iip1, jjp1, rlonu, rlatu, llm, pls,
191     . workvar, uvent, 0.0, jjm ,rlonv, rlatv, interbar )
192      ! 
193      varname = 'v'
194      vvent(:,:,:) = 0.0
195      CALL startget(varname, iip1, jjm, rlonv, rlatv, llm, pls,
196     . workvar, vvent, 0.0, jjp1, rlonu, rlatu, interbar )
197      !
198      varname = 't'
199      t3d(:,:,:) = 0.0
200      CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
201     . workvar, t3d, 0.0 , jjm, rlonu, rlatv , interbar )
202      !
203      WRITE(*,*) 'T3D min,max:',minval(t3d(:,:,:)),
204     .                          maxval(t3d(:,:,:))
205      varname = 'tpot'
206      tpot(:,:,:) = 0.0
207      CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
208     . pk, tpot, 0.0 , jjm, rlonu, rlatv , interbar )
209      !
210      WRITE(*,*) 'T3D min,max:',minval(t3d(:,:,:)),
211     .                          maxval(t3d(:,:,:))
212      WRITE(*,*) 'PLS min,max:',minval(pls(:,:,:)),
213     .                          maxval(pls(:,:,:))
214      DO l = 1, llm
215        DO j=1,jjp1
216          DO i =1, iip1-1
217           qsat(i,j,l) = q_sat(t3d(i,j,l),pls(i,j,l)/100. )
218          ENDDO
219          qsat(iip1,j,l) = qsat(1,j,l)
220        ENDDO
221      ENDDO
222      WRITE(*,*) 'QSAT min,max:',minval(qsat(:,:,:)),
223     .                           maxval(qsat(:,:,:))
224      !
225      WRITE(*,*) 'QSAT :', qsat(10,20,:)
226      !
227      varname = 'q'
228      qd(:,:,:) = 0.0
229      q3d(:,:,:,:) = 0.0
230      WRITE(*,*) 'QSAT min,max:',minval(qsat(:,:,:)),
231     .                           maxval(qsat(:,:,:))
232      CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
233     . qsat, qd, 0.0, jjm, rlonu, rlatv , interbar )
234      q3d(:,:,:,1) = qd(:,:,:)
235      !
236      varname = 'tsol'
237      ! This line needs to be replaced by a call to restget to get the values in the restart file
238      tsol(:) = 0.0
239      CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, tsol, 0.0,
240     .    jjm, rlonu, rlatv , interbar )
241      !
242      WRITE(*,*) 'TSOL construit :'
243      WRITE(*,'(48I3)') INT(TSOL(2:klon)-273)
244      !
245      varname = 'qsol'
246      qsol(:) = 0.0
247      CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, qsol, 0.0,
248     .   jjm, rlonu, rlatv , interbar )
249      !
250      varname = 'snow'
251      sn(:) = 0.0
252      CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, sn, 0.0,
253     .    jjm, rlonu, rlatv , interbar )
254      !
255      varname = 'rads'
256      radsol(:) = 0.0
257      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,radsol,0.0,
258     .    jjm, rlonu, rlatv , interbar )
259      !
260      varname = 'deltat'
261      deltat(:) = 0.0
262      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,deltat,0.0,
263     .     jjm, rlonu, rlatv , interbar )
264      !
265      varname = 'rugmer'
266      rugmer(:) = 0.0
267      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,rugmer,0.0,
268     .     jjm, rlonu, rlatv , interbar )
269      !
270      varname = 'agsno'
271      agesno(:) = 0.0
272      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,agesno,0.0,
273     .     jjm, rlonu, rlatv , interbar )
274
275      varname = 'zmea'
276      zmea(:) = 0.0
277      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zmea,0.0,
278     .     jjm, rlonu, rlatv , interbar )
279
280      varname = 'zstd'
281      zstd(:) = 0.0
282      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zstd,0.0,
283     .     jjm, rlonu, rlatv , interbar )
284      varname = 'zsig'
285      zsig(:) = 0.0
286      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zsig,0.0,
287     .     jjm, rlonu, rlatv , interbar )
288      varname = 'zgam'
289      zgam(:) = 0.0
290      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zgam,0.0,
291     .     jjm, rlonu, rlatv , interbar )
292      varname = 'zthe'
293      zthe(:) = 0.0
294      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zthe,0.0,
295     .     jjm, rlonu, rlatv , interbar )
296      varname = 'zpic'
297      zpic(:) = 0.0
298      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zpic,0.0,
299     .     jjm, rlonu, rlatv , interbar )
300      varname = 'zval'
301      zval(:) = 0.0
302      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zval,0.0,
303     .     jjm, rlonu, rlatv , interbar )
304c
305      rugsrel(:) = 0.0
306      IF(ok_orodr)  THEN
307        DO i = 1, iip1* jjp1
308         rugsrel(i) = MAX( 1.e-05, zstd(i)* zsig(i) /2. )
309        ENDDO
310
311      WRITE(*,*) 'OUT OF GET VARIABLE : Rugosite  rugsrel'
312      WRITE(*,'(49I1)') INT(rugsrel(:)*10)
313
314      ENDIF
315c
316c  Calcul intermediaire
317c
318      CALL massdair( p3d, masse  )
319c
320
321      DO  l = 1, llm
322        DO  i    = 1, iim
323          xppn(i) = aire( i, 1   ) * masse(  i     ,  1   , l )
324          xpps(i) = aire( i,jjp1 ) * masse(  i     , jjp1 , l )
325        ENDDO
326          xpn      = SUM(xppn)/apoln
327          xps      = SUM(xpps)/apols
328        DO i   = 1, iip1
329          masse(   i   ,   1     ,  l )   = xpn
330          masse(   i   ,   jjp1  ,  l )   = xps
331        ENDDO
332      ENDDO
333      q3d(iip1,:,:,:) = q3d(1,:,:,:)
334      phis(iip1,:) = phis(1,:)
335
336C Ecriture
337
338
339      CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   ,
340     *                tetagdiv, tetagrot , tetatemp              )
341      itau = 0
342      iday = dayref +itau/day_step
343      time = FLOAT(itau-(iday-dayref)*day_step)/day_step
344c     
345      IF(time.GT.1)  THEN
346       time = time - 1
347       iday = iday + 1
348      ENDIF
349      CALL geopot  ( ip1jmp1, tpot  , pk , pks,  phis  , phi   )
350     
351      CALL caldyn0 ( itau,uvent,vvent,tpot,psol,masse,pk,phis ,
352     *                phi,w, pbaru,pbarv,time+iday-dayref   )
353      print*,'sortie caldyn0'     
354      CALL dynredem0("start.nc",dayref,anneeref,phis,nqmx)
355      print*,'sortie dynredem0'
356      CALL dynredem1("start.nc",0.0,vvent,uvent,tpot,q3d,nqmx,masse ,
357     .                            psol)
358      print*,'sortie dynredem1'
359C
360C Ecriture etat initial physique
361C
362      phystep   = dtvr * FLOAT(iphysiq)
363      radpas    = NINT (86400./phystep/ FLOAT(nbapp_rad) )
364      co2_ppm   = 330.0
365      solaire   = 1370.0
366
367      call physdem(lonfi, latfi, phystep,radpas,co2_ppm,
368     .                   solaire,tsol, qsol,
369     .                   sn, radsol, deltat, rugmer,
370     .                   agesno, zmea, zstd, zsig,
371     .                   zgam, zthe, zpic, zval,
372     .                   rugsrel)
373cc
374
375C     Sortie Visu pour les champs dynamiques
376      print*,'sortie visu'
377      time_step = 1.
378      t_ops = 2.
379      t_wrt = 2.
380      itau = 2.
381      visu_file='Etat0_visu.nc'
382      CALL initdynav(visu_file,dayref,anneeref,time_step,
383     .              t_ops, t_wrt, nqmx, visuid)
384      CALL writedynav(visuid, nqmx, itau,vvent ,
385     .                uvent,tpot,pk,phi,q3d,masse,psol,phis)
386      print*,'entree histclo'
387      CALL histclo
388      RETURN
389      !
390      END SUBROUTINE etat0_netcdf
Note: See TracBrowser for help on using the repository browser.