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

Last change on this file since 3981 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.