source: LMDZ.3.3/branches/rel-1-0-patch/libf/dyn3d/etat0_netcdf.F @ 4855

Last change on this file since 4855 was 253, checked in by (none), 23 years ago

This commit was manufactured by cvs2svn to create branch
'rel-1-0-patch'.
  • 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      q3d(:,:,:,:) = 0.0
220      WRITE(*,*) 'QSAT min,max:',minval(qsat(:,:,:)),
221     .                           maxval(qsat(:,:,:))
222      CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
223     . qsat, qd, 0.0)
224      q3d(:,:,:,1) = qd(:,:,:)
225      !
226      varname = 'tsol'
227      ! This line needs to be replaced by a call to restget to get the values in the restart file
228      tsol(:) = 0.0
229      CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, tsol, 0.0)
230      !
231      WRITE(*,*) 'TSOL construit :'
232      WRITE(*,'(48I3)') INT(TSOL(2:klon)-273)
233      !
234      varname = 'qsol'
235      qsol(:) = 0.0
236      CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, qsol, 0.0)
237      !
238      varname = 'snow'
239      sn(:) = 0.0
240      CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, sn, 0.0)
241      !
242      varname = 'rads'
243      radsol(:) = 0.0
244      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,radsol,0.0)
245      !
246      varname = 'deltat'
247      deltat(:) = 0.0
248      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,deltat,0.0)
249      !
250      varname = 'rugmer'
251      rugmer(:) = 0.0
252      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,rugmer,0.0)
253      !
254      varname = 'agsno'
255      agesno(:) = 0.0
256      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,agesno,0.0)
257
258      varname = 'zmea'
259      zmea(:) = 0.0
260      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zmea,0.0)
261      varname = 'zstd'
262      zstd(:) = 0.0
263      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zstd,0.0)
264      varname = 'zsig'
265      zsig(:) = 0.0
266      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zsig,0.0)
267      varname = 'zgam'
268      zgam(:) = 0.0
269      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zgam,0.0)
270      varname = 'zthe'
271      zthe(:) = 0.0
272      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zthe,0.0)
273      varname = 'zpic'
274      zpic(:) = 0.0
275      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zpic,0.0)
276      varname = 'zval'
277      zval(:) = 0.0
278      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zval,0.0)
279      rugsrel(:) = 0.0
280C Calcul intermediaire
281c
282      CALL massdair( p3d, masse  )
283c
284
285      print *,' ALPHAX ',alphax
286
287      DO  l = 1, llm
288        DO  i    = 1, iim
289          xppn(i) = aire( i, 1   ) * masse(  i     ,  1   , l )
290          xpps(i) = aire( i,jjp1 ) * masse(  i     , jjp1 , l )
291        ENDDO
292          xpn      = SUM(xppn)/apoln
293          xps      = SUM(xpps)/apols
294        DO i   = 1, iip1
295          masse(   i   ,   1     ,  l )   = xpn
296          masse(   i   ,   jjp1  ,  l )   = xps
297        ENDDO
298      ENDDO
299      q3d(iip1,:,:,:) = q3d(1,:,:,:)
300      phis(iip1,:) = phis(1,:)
301
302C Ecriture
303
304
305      CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   ,
306     *                tetagdiv, tetagrot , tetatemp              )
307      print*,'sortie inidissip'
308      itau = 0
309      iday = dayref +itau/day_step
310      time = FLOAT(itau-(iday-dayref)*day_step)/day_step
311c     
312      IF(time.GT.1)  THEN
313       time = time - 1
314       iday = iday + 1
315      ENDIF
316      CALL geopot  ( ip1jmp1, tpot  , pk , pks,  phis  , phi   )
317      print*,'sortie geopot'
318     
319      CALL caldyn0 ( itau,uvent,vvent,tpot,psol,masse,pk,phis ,
320     *                phi,w, pbaru,pbarv,time+iday-dayref   )
321       print*,'sortie caldyn0'     
322      CALL dynredem0("start.nc",dayref,anneeref,phis,nqmx)
323      print*,'sortie dynredem0'
324      CALL dynredem1("start.nc",0.0,vvent,uvent,tpot,q3d,nqmx,masse ,
325     .                            psol)
326      print*,'sortie dynredem1'
327C
328C Ecriture etat initial physique
329C
330      phystep   = dtvr * FLOAT(iphysiq)
331      radpas    = NINT (86400./phystep/ FLOAT(nbapp_rad) )
332      co2_ppm   = 330.0
333      solaire   = 1370.0
334
335      call physdem(lonfi, latfi, phystep,radpas,co2_ppm,
336     .                   solaire,tsol, qsol,
337     .                   sn, radsol, deltat, rugmer,
338     .                   agesno, zmea, zstd, zsig,
339     .                   zgam, zthe, zpic, zval,
340     .                   rugsrel)
341
342C     Sortie Visu pour les champs dynamiques
343      print*,'sortie visu'
344      time_step = 1.
345      t_ops = 2.
346      t_wrt = 2.
347      itau = 2.
348      visu_file='Etat0_visu.nc'
349      CALL initdynav(visu_file,dayref,anneeref,time_step,
350     .              t_ops, t_wrt, nqmx, visuid)
351      CALL writedynav(visuid, nqmx, itau,vvent ,
352     .                uvent,tpot,pk,phi,q3d,masse,psol,phis)
353      print*,'entree histclo'
354      CALL histclo
355      RETURN
356      !
357      END SUBROUTINE etat0_netcdf
Note: See TracBrowser for help on using the repository browser.