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

Last change on this file since 57 was 50, checked in by lmdz, 25 years ago

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