source: LMDZ.3.3/branches/LF/libf/dyn3d/etat0_netcdf.F @ 4190

Last change on this file since 4190 was 2, checked in by lmdz, 25 years ago

Initial revision
  • 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, llmp1)
23      REAL :: uvent(iip1, jjp1, llmp1)
24      REAL :: vvent(iip1, jjm, llmp1)
25      REAL :: t3d(iip1, jjp1, llmp1), tpot(iip1, jjp1, llmp1)
26      REAL :: q3d(iip1, jjp1, llmp1,nqmx), qsat(iip1, jjp1, llmp1)
27      REAL :: tsol(klon), qsol(klon), sn(klon), radsol(klon)
28      REAL :: deltat(klon), rugmer(klon), agesno(klon), zmea(klon)
29      REAL :: zstd(klon), zsig(klon), zgam(klon), zthe(klon)
30      REAL :: zpic(klon), zval(klon), rugsrel(klon)
31      REAL :: qd(iip1, jjp1, llmp1)
32      !
33      CHARACTER*80 :: varname
34      !
35      INTEGER :: i,j, ig, l
36      REAL :: xpi
37      !
38      REAL :: alpha(iip1,jjp1,llm),beta(iip1,jjp1,llm)
39      REAL :: pk(iip1,jjp1,llm), pls(iip1,jjp1,llm), pks(ip1jmp1)
40      REAL :: workvar(iip1,jjp1,llm)
41      !
42      REAL ::  prefkap, unskap
43      !
44      REAL :: q_sat
45      EXTERNAL q_sat
46      real :: time_step,t_ops,t_wrt
47
48#include "comdissnew.h"
49#include "control.h"
50#include "serre.h"
51#include "clesph0.h"
52
53      INTEGER  ::        longcles
54      PARAMETER      ( longcles  = 20 )
55      REAL :: clesphy0 ( longcles       )
56      REAL :: p(iip1,jjp1,llmp1)
57      INTEGER :: itau, iday
58      REAL :: masse(iip1,jjp1,llm)
59      REAL :: xpn,xps,xppn(iim),xpps(iim)
60      real :: time
61      REAL :: phi(ip1jmp1,llm)
62      REAL :: pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
63      REAL :: w(ip1jmp1,llm)
64      REAL ::phystep,radpas,co2_ppm,solaire
65
66      CHARACTER*80 :: visu_file
67      INTEGER :: visuid
68      !
69      !   Constantes
70      !
71      pi     = 4. * ATAN(1.)
72      rad    = 6371229.
73      omeg   = 4.* ASIN(1.)/(24.*3600.)
74      g      = 9.8
75      daysec = 86400.
76      kappa  = 0.2857143
77      cpp    = 1004.70885
78      !
79      preff     = 101325.
80      unskap = 1./kappa
81      !
82      jmp1    = jjm + 1
83      !
84      !    Construct a grid
85      !
86
87      CALL defrun_new(99,.TRUE.,clesphy0)
88
89      dtvr   = daysec/FLOAT(day_step)
90      print*,'dtvr',dtvr
91
92      CALL inicons0()
93      CALL inigeom()
94      !
95      CALL inifilr()
96      !
97      latfi(1) = ASIN(1.0)
98      DO j = 2, jjm
99        DO i = 1, iim
100          latfi((j-2)*iim+1+i)=  rlatu(j)
101        ENDDO
102      ENDDO
103      latfi(klon) = - ASIN(1.0)
104      !
105      lonfi(1) = 0.0
106      DO j = 2, jjm
107        DO i = 1, iim
108          lonfi((j-2)*iim+1+i) =  rlonv(i)
109        ENDDO
110      ENDDO
111      lonfi(klon) = 0.0
112      !
113      xpi = 2.0 * ASIN(1.0)
114      DO ig = 1, klon
115        latfi(ig) = latfi(ig) * 180.0 / xpi
116        lonfi(ig) = lonfi(ig) * 180.0 / xpi
117      ENDDO
118      !
119      varname = 'relief'
120      ! This line needs to be replaced by a call to restget to get the values in the restart file
121      orog(:,:) = 0.0
122       CALL startget(varname, iip1, jjp1, rlonv, rlatu, orog, 0.0)
123      !
124      WRITE(*,*) 'OUT OF GET VARIABLE : Relief'
125      WRITE(*,'(49I1)') INT(orog(:,:))
126      !
127      varname = 'rugosite'
128      ! This line needs to be replaced by a call to restget to get the values in the restart file
129      rugo(:,:) = 0.0
130       CALL startget(varname, iip1, jjp1, rlonv, rlatu, rugo, 0.0)
131      !
132      WRITE(*,*) 'OUT OF GET VARIABLE : Rugosite'
133      WRITE(*,'(49I1)') INT(rugo(:,:)*10)
134      !
135      varname = 'masque'
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, masque, 0.0)
139      !
140      WRITE(*,*) 'MASQUE construit : Masque'
141      WRITE(*,'(49I1)') INT(masque(:,:))
142      !
143      !
144      varname = 'psol'
145      psol(:,:) = 0.0
146      CALL startget(varname, iip1, jjp1, rlonv, rlatu, psol, 0.0)
147      !
148      !  Compute here the pressure on the intermediate levels. One would expect that this is available in the GCM
149      !  anyway.
150      !
151      WRITE(*,*) 'PSOL :', psol(10,20)
152      WRITE(*,*) ap(:), bp(:)
153      CALL pression(ip1jmp1, ap, bp, psol, p3d)
154      WRITE(*,*) 'P3D :', p3d(10,20,:)
155      CALL exner_hyb(ip1jmp1, psol, p3d, alpha, beta, pks, pk, workvar)
156      WRITE(*,*) 'PK:', pk(10,20,:)
157      !
158      !
159      !
160      prefkap =  preff  ** kappa
161      WRITE(*,*) 'unskap, cpp,  preff :', unskap, cpp,  preff
162      DO l = 1, llm
163        DO j=1,jjp1
164          DO i =1, iip1
165            pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap
166           ENDDO
167        ENDDO
168      ENDDO
169      !
170      WRITE(*,*) 'PLS :', pls(10,20,:)
171      !
172      varname = 'surfgeo'
173      phis(:,:) = 0.0
174      CALL startget(varname, iip1, jjp1, rlonv, rlatu, phis, 0.0)
175      phis = phis * 9.81
176      !
177      varname = 'u'
178      uvent(:,:,:) = 0.0
179      CALL startget(varname, iip1, jjp1, rlonu, rlatu, llm, pls,
180     . workvar, uvent, 0.0)
181      ! 
182      varname = 'v'
183      vvent(:,:,:) = 0.0
184      CALL startget(varname, iip1, jjm, rlonv, rlatv, llm, pls,
185     . workvar, vvent, 0.0)
186      !
187      varname = 't'
188      t3d(:,:,:) = 0.0
189      CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
190     . workvar, t3d, 0.0)
191      !
192      WRITE(*,*) 'T3D min,max:',minval(t3d(:,:,:)),
193     .                          maxval(t3d(:,:,:))
194      varname = 'tpot'
195      tpot(:,:,:) = 0.0
196      CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
197     . pk, tpot, 0.0)
198      !
199      WRITE(*,*) 'T3D min,max:',minval(t3d(:,:,:)),
200     .                          maxval(t3d(:,:,:))
201      WRITE(*,*) 'PLS min,max:',minval(pls(:,:,:)),
202     .                          maxval(pls(:,:,:))
203      DO l = 1, llm
204        DO j=1,jjp1
205          DO i =1, iip1-1
206           qsat(i,j,l) = q_sat(t3d(i,j,l),pls(i,j,l)/100. )
207          ENDDO
208          qsat(iip1,j,l) = qsat(1,j,l)
209        ENDDO
210      ENDDO
211      WRITE(*,*) 'QSAT min,max:',minval(qsat(:,:,:)),
212     .                           maxval(qsat(:,:,:))
213      !
214      WRITE(*,*) 'QSAT :', qsat(10,20,:)
215      !
216      varname = 'q'
217      qd(:,:,:) = 0.0
218      WRITE(*,*) 'QSAT min,max:',minval(qsat(:,:,:)),
219     .                           maxval(qsat(:,:,:))
220      CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
221     . qsat, qd, 0.0)
222      q3d(:,:,:,1) = qd(:,:,:)
223      !
224      varname = 'tsol'
225      ! This line needs to be replaced by a call to restget to get the values in the restart file
226      tsol(:) = 0.0
227      CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, tsol, 0.0)
228      !
229      WRITE(*,*) 'TSOL construit :'
230      WRITE(*,'(48I3)') INT(TSOL(2:klon)-273)
231      !
232      varname = 'qsol'
233      qsol(:) = 0.0
234      CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, qsol, 0.0)
235      !
236      varname = 'snow'
237      sn(:) = 0.0
238      CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, sn, 0.0)
239      !
240      varname = 'rads'
241      radsol(:) = 0.0
242      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,radsol,0.0)
243      !
244      varname = 'deltat'
245      deltat(:) = 0.0
246      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,deltat,0.0)
247      !
248      varname = 'rugmer'
249      rugmer(:) = 0.0
250      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,rugmer,0.0)
251      !
252      varname = 'agsno'
253      agesno(:) = 0.0
254      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,agesno,0.0)
255
256      varname = 'zmea'
257      zmea(:) = 0.0
258      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zmea,0.0)
259      varname = 'zstd'
260      zstd(:) = 0.0
261      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zstd,0.0)
262      varname = 'zsig'
263      zsig(:) = 0.0
264      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zsig,0.0)
265      varname = 'zgam'
266      zgam(:) = 0.0
267      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zgam,0.0)
268      varname = 'zthe'
269      zthe(:) = 0.0
270      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zthe,0.0)
271      varname = 'zpic'
272      zpic(:) = 0.0
273      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zpic,0.0)
274      varname = 'zval'
275      zval(:) = 0.0
276      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zval,0.0)
277      rugsrel(:) = 0.0
278C Calcul intermediaire
279c
280      CALL massdair( p3d, masse  )
281c
282
283      print *,' ALPHAX ',alphax
284
285      IF( alphax.NE.0. )   THEN
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      ENDIF
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.