source: LMDZ.3.3/tags/IPSL-CM4_0/libf/dyn3d/etat0_netcdf.F @ 331

Last change on this file since 331 was 331, checked in by (none), 22 years ago

This commit was manufactured by cvs2svn to create tag 'IPSL-CM4_0'.
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 20.0 KB
Line 
1c
2c $Header$
3c
4      SUBROUTINE etat0_netcdf (interbar, masque, pctsrf)
5   
6      USE startvar
7      USE ioipsl
8      !
9      IMPLICIT NONE
10      !
11#include "netcdf.inc"
12#include "dimensions.h"
13#include "paramet.h"
14      !
15      !
16!      INTEGER, PARAMETER :: KIDIA=1, KFDIA=iim*(jjm-1)+2,
17!     .KLON=KFDIA-KIDIA+1,KLEV=llm
18      !
19#include "comgeom2.h"
20#include "comvert.h"
21#include "comconst.h"
22#include "indicesol.h"
23#include "dimphy.h"
24#include "dimsoil.h"
25      !
26      LOGICAL interbar
27      REAL :: latfi(klon), lonfi(klon)
28      REAL :: orog(iip1,jjp1), rugo(iip1,jjp1), masque(iip1,jjp1),
29     . psol(iip1, jjp1), phis(iip1, jjp1)
30      REAL :: p3d(iip1, jjp1, llm+1)
31      REAL :: uvent(iip1, jjp1, llm)
32      REAL :: vvent(iip1, jjm, llm)
33      REAL :: t3d(iip1, jjp1, llm), tpot(iip1, jjp1, llm)
34      REAL :: q3d(iip1, jjp1, llm,nqmx), qsat(iip1, jjp1, llm)
35      REAL :: tsol(klon), qsol(klon), sn(klon)
36      REAL :: tsolsrf(klon,nbsrf), qsolsrf(klon,nbsrf),snsrf(klon,nbsrf)
37      REAL :: albe(klon,nbsrf), evap(klon,nbsrf)
38      REAL :: tsoil(klon,nsoilmx,nbsrf)
39      REAL :: radsol(klon),rain_fall(klon), snow_fall(klon)
40      REAL :: solsw(klon), sollw(klon)
41      REAL :: deltat(klon), frugs(klon,nbsrf), agesno(klon),rugmer(klon)
42      REAL :: zmea(iip1*jjp1), zstd(iip1*jjp1)
43      REAL :: zsig(iip1*jjp1), zgam(iip1*jjp1), zthe(iip1*jjp1)
44      REAL :: zpic(iip1*jjp1), zval(iip1*jjp1), rugsrel(iip1*jjp1)
45      REAL :: qd(iip1, jjp1, llm)
46      REAL :: pctsrf(klon, nbsrf)
47      REAL :: t_ancien(klon,klev), q_ancien(klon,klev)      !
48      ! declarations pour lecture glace de mer
49      REAL :: rugv(klon)
50      INTEGER :: iml_lic, jml_lic, llm_tmp, ttm_tmp, iret
51      INTEGER :: itaul(1), fid
52      REAL :: lev(1), date, dt
53      REAL, ALLOCATABLE, DIMENSION(:,:) :: lon_lic, lat_lic
54      REAL, ALLOCATABLE, DIMENSION(:)  :: dlon_lic, dlat_lic
55      REAL, ALLOCATABLE, DIMENSION (:,:) :: fraclic
56      REAL :: flic_tmp(iip1, jjp1)
57      REAL :: champint(iim, jjp1)
58      !
59
60      CHARACTER*80 :: varname
61      !
62      INTEGER :: i,j, ig, l, ji,ii1,ii2
63      REAL :: xpi
64      !
65      REAL :: alpha(iip1,jjp1,llm),beta(iip1,jjp1,llm)
66      REAL :: pk(iip1,jjp1,llm), pls(iip1,jjp1,llm), pks(ip1jmp1)
67      REAL :: workvar(iip1,jjp1,llm)
68      !
69      REAL ::  prefkap, unskap
70      !
71      REAL :: q_sat
72      EXTERNAL q_sat
73      real :: time_step,t_ops,t_wrt
74
75#include "comdissnew.h"
76#include "control.h"
77#include "serre.h"
78#include "clesph0.h"
79
80      INTEGER  ::        longcles
81      PARAMETER      ( longcles  = 20 )
82      REAL :: clesphy0 ( longcles       )
83      REAL :: p(iip1,jjp1,llm)
84      INTEGER :: itau, iday
85      REAL :: masse(iip1,jjp1,llm)
86      REAL :: xpn,xps,xppn(iim),xpps(iim)
87      real :: time
88      REAL :: phi(ip1jmp1,llm)
89      REAL :: pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
90      REAL :: w(ip1jmp1,llm)
91      REAL ::phystep,co2_ppm,solaire
92      INTEGER :: radpas
93       real zrel(iip1*jjp1),chmin,chmax
94
95      CHARACTER*80 :: visu_file
96      INTEGER :: visuid
97
98! pour la lecture du fichier masque ocean
99      integer :: nid_o2a
100      logical :: couple = .false.
101      INTEGER :: iml_omask, jml_omask
102      REAL, ALLOCATABLE, DIMENSION(:,:) :: lon_omask, lat_omask
103      REAL, ALLOCATABLE, DIMENSION(:)  :: dlon_omask, dlat_omask
104      REAL, ALLOCATABLE, DIMENSION (:,:) :: ocemask, ocetmp
105      real, dimension(klon) :: ocemask_fi
106      integer :: isst(klon-2)
107
108      !
109      !   Constantes
110      !
111      pi     = 4. * ATAN(1.)
112      rad    = 6371229.
113      omeg   = 4.* ASIN(1.)/(24.*3600.)
114      g      = 9.8
115      daysec = 86400.
116      kappa  = 0.2857143
117      cpp    = 1004.70885
118      !
119      preff     = 101325.
120      unskap = 1./kappa
121      !
122      jmp1    = jjm + 1
123      !
124      !    Construct a grid
125      !
126
127!      CALL defrun_new(99,.TRUE.,clesphy0)
128      CALL conf_gcm( 99, .TRUE. , clesphy0 )
129
130      dtvr   = daysec/FLOAT(day_step)
131      print*,'dtvr',dtvr
132
133      CALL inicons0()
134      CALL inigeom()
135      !
136      CALL inifilr()
137      !
138      latfi(1) = ASIN(1.0)
139      DO j = 2, jjm
140        DO i = 1, iim
141          latfi((j-2)*iim+1+i)=  rlatu(j)
142        ENDDO
143      ENDDO
144      latfi(klon) = - ASIN(1.0)
145      !
146      lonfi(1) = 0.0
147      DO j = 2, jjm
148        DO i = 1, iim
149          lonfi((j-2)*iim+1+i) =  rlonv(i)
150        ENDDO
151      ENDDO
152      lonfi(klon) = 0.0
153      !
154      xpi = 2.0 * ASIN(1.0)
155      DO ig = 1, klon
156        latfi(ig) = latfi(ig) * 180.0 / xpi
157        lonfi(ig) = lonfi(ig) * 180.0 / xpi
158      ENDDO
159      !
160      varname = 'relief'
161      ! This line needs to be replaced by a call to restget to get the values in the restart file
162      orog(:,:) = 0.0
163       CALL startget(varname, iip1, jjp1, rlonv, rlatu, orog, 0.0 ,
164     , jjm ,rlonu,rlatv , interbar )
165      !
166      WRITE(*,*) 'OUT OF GET VARIABLE : Relief'
167!      WRITE(*,'(49I1)') INT(orog(:,:))
168      !
169      varname = 'rugosite'
170      ! This line needs to be replaced by a call to restget to get the values in the restart file
171      rugo(:,:) = 0.0
172       CALL startget(varname, iip1, jjp1, rlonv, rlatu, rugo, 0.0 ,
173     , jjm, rlonu,rlatv , interbar )
174      !
175      WRITE(*,*) 'OUT OF GET VARIABLE : Rugosite'
176!      WRITE(*,'(49I1)') INT(rugo(:,:)*10)
177      !
178      varname = 'masque'
179      ! This line needs to be replaced by a call to restget to get the values in the restart file
180      masque(:,:) = 0.0
181       CALL startget(varname, iip1, jjp1, rlonv, rlatu, masque, 0.0,
182     , jjm ,rlonu,rlatv , interbar )
183      !
184      WRITE(*,*) 'MASQUE construit : Masque'
185!      WRITE(*,'(49I1)') INT(masque(:,:))
186      !
187      !
188
189
190C
191C on initialise les sous surfaces
192C
193      pctsrf=0.
194      !cree le masque a partir du fichier relief
195      varname = 'zmasq'
196      zmasq(:) = 0.
197      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zmasq,0.0,
198     , jjm ,rlonu,rlatv , interbar )
199      WHERE (zmasq(1 : klon) .LT. EPSFRA)
200          zmasq(1 : klon) = 0.
201      END WHERE
202      WHERE (1. - zmasq(1 : klon) .LT. EPSFRA)
203          zmasq(1 : klon) = 1.
204      END WHERE
205
206
207
208
209      varname = 'psol'
210      psol(:,:) = 0.0
211      CALL startget(varname, iip1, jjp1, rlonv, rlatu, psol, 0.0 ,
212     , jjm ,rlonu,rlatv , interbar )
213      !
214      !  Compute here the pressure on the intermediate levels. One would expect that this is available in the GCM
215      !  anyway.
216      !
217!      WRITE(*,*) 'PSOL :', psol(10,20)
218!      WRITE(*,*) ap(:), bp(:)
219      CALL pression(ip1jmp1, ap, bp, psol, p3d)
220!      WRITE(*,*) 'P3D :', p3d(10,20,:)
221      CALL exner_hyb(ip1jmp1, psol, p3d, alpha, beta, pks, pk, workvar)
222!      WRITE(*,*) 'PK:', pk(10,20,:)
223      !
224      !
225      !
226      prefkap =  preff  ** kappa
227!      WRITE(*,*) 'unskap, cpp,  preff :', unskap, cpp,  preff
228      DO l = 1, llm
229        DO j=1,jjp1
230          DO i =1, iip1
231            pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap
232           ENDDO
233        ENDDO
234      ENDDO
235      !
236!      WRITE(*,*) 'PLS :', pls(10,20,:)
237      !
238      varname = 'surfgeo'
239      phis(:,:) = 0.0
240      CALL startget(varname, iip1, jjp1, rlonv, rlatu, phis, 0.0 ,
241     , jjm ,rlonu,rlatv, interbar )
242      !
243      varname = 'u'
244      uvent(:,:,:) = 0.0
245      CALL startget(varname, iip1, jjp1, rlonu, rlatu, llm, pls,
246     . workvar, uvent, 0.0, jjm ,rlonv, rlatv, interbar )
247      ! 
248      varname = 'v'
249      vvent(:,:,:) = 0.0
250      CALL startget(varname, iip1, jjm, rlonv, rlatv, llm, pls,
251     . workvar, vvent, 0.0, jjp1, rlonu, rlatu, interbar )
252      !
253      varname = 't'
254      t3d(:,:,:) = 0.0
255      CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
256     . workvar, t3d, 0.0 , jjm, rlonu, rlatv , interbar )
257      !
258      WRITE(*,*) 'T3D min,max:',minval(t3d(:,:,:)),
259     .                          maxval(t3d(:,:,:))
260      varname = 'tpot'
261      tpot(:,:,:) = 0.0
262      CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
263     . pk, tpot, 0.0 , jjm, rlonu, rlatv , interbar )
264      !
265      WRITE(*,*) 'T3D min,max:',minval(t3d(:,:,:)),
266     .                          maxval(t3d(:,:,:))
267      WRITE(*,*) 'PLS min,max:',minval(pls(:,:,:)),
268     .                          maxval(pls(:,:,:))
269      DO l = 1, llm
270        DO j=1,jjp1
271          DO i =1, iip1-1
272           qsat(i,j,l) = q_sat(t3d(i,j,l),pls(i,j,l)/100. )
273          ENDDO
274          qsat(iip1,j,l) = qsat(1,j,l)
275        ENDDO
276      ENDDO
277      WRITE(*,*) 'QSAT min,max:',minval(qsat(:,:,:)),
278     .                           maxval(qsat(:,:,:))
279      !
280      WRITE(*,*) 'QSAT :', qsat(10,20,:)
281      !
282      varname = 'q'
283      qd(:,:,:) = 0.0
284      q3d(:,:,:,:) = 0.0
285      WRITE(*,*) 'QSAT min,max:',minval(qsat(:,:,:)),
286     .                           maxval(qsat(:,:,:))
287      CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
288     . qsat, qd, 0.0, jjm, rlonu, rlatv , interbar )
289      q3d(:,:,:,1) = qd(:,:,:)
290      !
291      varname = 'tsol'
292      ! This line needs to be replaced by a call to restget to get the values in the restart file
293      tsol(:) = 0.0
294      CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, tsol, 0.0,
295     .    jjm, rlonu, rlatv , interbar )
296      !
297      WRITE(*,*) 'TSOL construit :'
298!      WRITE(*,'(48I3)') INT(TSOL(2:klon)-273)
299      !
300      varname = 'qsol'
301      qsol(:) = 0.0
302      CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, qsol, 0.0,
303     .   jjm, rlonu, rlatv , interbar )
304      !
305      varname = 'snow'
306      sn(:) = 0.0
307      CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, sn, 0.0,
308     .    jjm, rlonu, rlatv , interbar )
309      !
310      varname = 'rads'
311      radsol(:) = 0.0
312      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,radsol,0.0,
313     .    jjm, rlonu, rlatv , interbar )
314      !
315      varname = 'deltat'
316      deltat(:) = 0.0
317      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,deltat,0.0,
318     .     jjm, rlonu, rlatv , interbar )
319      !
320      varname = 'rugmer'
321      rugmer(:) = 0.0
322      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,rugmer,0.0,
323     .     jjm, rlonu, rlatv , interbar )
324      !
325      varname = 'agsno'
326      agesno(:) = 0.0
327      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,agesno,0.0,
328     .     jjm, rlonu, rlatv , interbar )
329
330      varname = 'zmea'
331      zmea(:) = 0.0
332      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zmea,0.0,
333     .     jjm, rlonu, rlatv , interbar )
334
335      varname = 'zstd'
336      zstd(:) = 0.0
337      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zstd,0.0,
338     .     jjm, rlonu, rlatv , interbar )
339      varname = 'zsig'
340      zsig(:) = 0.0
341      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zsig,0.0,
342     .     jjm, rlonu, rlatv , interbar )
343      varname = 'zgam'
344      zgam(:) = 0.0
345      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zgam,0.0,
346     .     jjm, rlonu, rlatv , interbar )
347      varname = 'zthe'
348      zthe(:) = 0.0
349      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zthe,0.0,
350     .     jjm, rlonu, rlatv , interbar )
351      varname = 'zpic'
352      zpic(:) = 0.0
353      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zpic,0.0,
354     .     jjm, rlonu, rlatv , interbar )
355      varname = 'zval'
356      zval(:) = 0.0
357      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zval,0.0,
358     .     jjm, rlonu, rlatv , interbar )
359c
360      rugsrel(:) = 0.0
361      IF(ok_orodr)  THEN
362        DO i = 1, iip1* jjp1
363         rugsrel(i) = MAX( 1.e-05, zstd(i)* zsig(i) /2. )
364        ENDDO
365      ENDIF
366C
367C En cas de simulation couplee, lecture du masque ocean issu du modele ocean
368C utilise pour calculer les poids et pour assurer l'adequation entre les
369C fractions d'ocean vu par l'atmosphere et l'ocean
370C
371
372      write(*,*)'Essai de lecture masque ocean'
373      iret = nf_open("o2a.nc", NF_NOWRITE, nid_o2a)
374      if (iret .ne. 0) then
375        write(*,*)'ATTENTION!! pas de fichier o2a.nc trouve'
376        write(*,*)'Run force'
377      else
378        couple = .true.
379        iret = nf_close(nid_o2a)
380        call flininfo("o2a.nc", iml_omask, jml_omask, llm_tmp, ttm_tmp
381     $    , nid_o2a)
382        if (iml_omask /= iim .or. jml_omask /= jjp1) then
383          write(*,*)'Dimensions non compatibles pour masque ocean'
384          write(*,*)'iim = ',iim,' iml_omask = ',iml_omask
385          write(*,*)'jjp1 = ',jjp1,' jml_omask = ',jml_omask
386          stop
387        endif
388        ALLOCATE(lat_omask(iml_omask, jml_omask), stat=iret)
389        ALLOCATE(lon_omask(iml_omask, jml_omask), stat=iret)
390        ALLOCATE(dlon_omask(iml_omask), stat=iret)
391        ALLOCATE(dlat_omask(jml_omask), stat=iret)
392        ALLOCATE(ocemask(iml_omask, jml_omask), stat=iret)
393        ALLOCATE(ocetmp(iml_omask, jml_omask), stat=iret)
394        CALL flinopen("o2a.nc", .FALSE., iml_omask, jml_omask, llm_tmp
395     $    , lon_omask, lat_omask, lev, ttm_tmp, itaul, date, dt, fid)
396        CALL flinget(fid, 'OceMask', iml_omask, jml_omask, llm_tmp,
397     $      ttm_tmp, 1, 1, ocetmp)
398        CALL flinclo(fid)
399        dlon_omask(1 : iml_omask) = lon_omask(1 : iml_omask, 1)
400        dlat_omask(1 : jml_omask) = lat_omask(1 , 1 : jml_omask)
401        ocemask = ocetmp
402        if (dlat_omask(1) < dlat_omask(jml_omask)) then
403          do j = 1, jml_omask
404            ocemask(:,j) = ocetmp(:,jml_omask-j+1)
405          enddo
406        endif
407C
408C passage masque ocean a la grille physique
409C
410        ocemask_fi(1) = ocemask(1,1)
411        do j = 2, jjm
412          do i = 1, iim
413            ocemask_fi((j-2)*iim + i + 1) = ocemask(i,j)
414          enddo
415        enddo
416        ocemask_fi(klon) = ocemask(1,jjp1)
417        zmasq = 1. - ocemask_fi
418        isst = 0
419        where (ocemask_fi(2:klon-1) >0.) isst = 1
420        write(45,'(72i1)')isst
421      endif
422
423
424C
425C lecture du fichier glace de terre pour fixer la fraction de terre
426C et de glace de terre
427C
428      CALL flininfo("landiceref.nc", iml_lic, jml_lic,llm_tmp, ttm_tmp
429     $    , fid)
430      ALLOCATE(lat_lic(iml_lic, jml_lic), stat=iret)
431      ALLOCATE(lon_lic(iml_lic, jml_lic), stat=iret)
432      ALLOCATE(dlon_lic(iml_lic), stat=iret)
433      ALLOCATE(dlat_lic(jml_lic), stat=iret)
434      ALLOCATE(fraclic(iml_lic, jml_lic), stat=iret)
435      CALL flinopen("landiceref.nc", .FALSE., iml_lic, jml_lic, llm_tmp
436     $    , lon_lic, lat_lic, lev, ttm_tmp, itaul, date, dt, fid)
437      CALL flinget(fid, 'landice', iml_lic, jml_lic, llm_tmp, ttm_tmp
438     $    , 1, 1, fraclic)
439      CALL flinclo(fid)
440C
441C interpolation sur la grille T du modele
442C
443      WRITE(*,*) 'dimensions de landice iml_lic, jml_lic : ',
444     $    iml_lic, jml_lic
445c
446C sil les coordonnees sont en degres, on les transforme
447C
448      IF( MAXVAL( lon_lic(:,:) ) .GT. 2.0 * asin(1.0) )  THEN
449          lon_lic(:,:) = lon_lic(:,:) * 2.0* ASIN(1.0) / 180.
450      ENDIF
451      IF( maxval( lat_lic(:,:) ) .GT. 2.0 * asin(1.0)) THEN
452          lat_lic(:,:) = lat_lic(:,:) * 2.0 * asin(1.0) / 180.
453      ENDIF
454
455      dlon_lic(1 : iml_lic) = lon_lic(1 : iml_lic, 1)
456      dlat_lic(1 : jml_lic) = lat_lic(1 , 1 : jml_lic)
457C
458      CALL grille_m(iml_lic, jml_lic, dlon_lic, dlat_lic, fraclic
459     $    ,iim, jjp1,
460     $    rlonv, rlatu, flic_tmp(1 : iim, 1 : jjp1))
461c$$$      flic_tmp(1 : iim, 1 : jjp1) = champint(1: iim, 1 : jjp1)
462      flic_tmp(iip1, 1 : jjp1) = flic_tmp(1 , 1 : jjp1)
463C
464C passage sur la grille physique
465C
466      CALL gr_dyn_fi(1, iip1, jjp1, klon, flic_tmp,
467     $    pctsrf(1:klon, is_lic))
468C adequation avec le maque terre/mer
469c      zmasq(157) = 0.
470      WHERE (pctsrf(1 : klon, is_lic) .LT. EPSFRA )
471          pctsrf(1 : klon, is_lic) = 0.
472      END WHERE
473      WHERE (zmasq( 1 : klon) .LT. EPSFRA)
474          pctsrf(1 : klon, is_lic) = 0.
475      END WHERE
476      pctsrf(1 : klon, is_ter) = zmasq(1 : klon)
477      DO ji = 1, klon
478        IF (zmasq(ji) .GT. EPSFRA) THEN
479            IF ( pctsrf(ji, is_lic) .GE. zmasq(ji)) THEN
480                pctsrf(ji, is_lic) = zmasq(ji)
481                pctsrf(ji, is_ter) = 0.
482            ELSE
483                pctsrf(ji,is_ter) = zmasq(ji) - pctsrf(ji, is_lic)
484                IF (pctsrf(ji,is_ter) .LT. EPSFRA) THEN
485                    pctsrf(ji,is_ter) = 0.
486                    pctsrf(ji, is_lic) = zmasq(ji)
487                ENDIF
488            ENDIF
489        ENDIF
490      END DO
491C
492C sous surface ocean et glace de mer (pour demarrer on met glace de mer a 0)
493C
494      pctsrf(1 : klon, is_oce) = (1. - zmasq(1 : klon))
495
496
497      WHERE (pctsrf(1 : klon, is_oce) .LT. EPSFRA)
498          pctsrf(1 : klon, is_oce) = 0.
499      END WHERE
500
501      if (couple) pctsrf(1 : klon, is_oce) = ocemask_fi(1 : klon)
502
503      isst = 0
504      where (pctsrf(2:klon-1,is_oce) >0.) isst = 1
505      write(46,'(72i1)')isst
506C
507C verif que somme des sous surface = 1
508C
509      ji=count( (abs( sum(pctsrf(1 : klon, 1 : nbsrf), dim = 2)) - 1.0 )
510     $    .GT. EPSFRA)
511      IF (ji .NE. 0) THEN
512          WRITE(*,*) 'pb repartition sous maille pour ',ji,' points'
513      ENDIF
514
515
516
517
518
519C Calcul intermediaire
520c
521      CALL massdair( p3d, masse  )
522c
523
524      print *,' ALPHAX ',alphax
525
526      DO  l = 1, llm
527        DO  i    = 1, iim
528          xppn(i) = aire( i, 1   ) * masse(  i     ,  1   , l )
529          xpps(i) = aire( i,jjp1 ) * masse(  i     , jjp1 , l )
530        ENDDO
531          xpn      = SUM(xppn)/apoln
532          xps      = SUM(xpps)/apols
533        DO i   = 1, iip1
534          masse(   i   ,   1     ,  l )   = xpn
535          masse(   i   ,   jjp1  ,  l )   = xps
536        ENDDO
537      ENDDO
538      q3d(iip1,:,:,:) = q3d(1,:,:,:)
539      phis(iip1,:) = phis(1,:)
540
541C Ecriture
542
543
544      CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   ,
545     *                tetagdiv, tetagrot , tetatemp              )
546      print*,'sortie inidissip'
547      itau = 0
548      iday = dayref +itau/day_step
549      time = FLOAT(itau-(iday-dayref)*day_step)/day_step
550c     
551      IF(time.GT.1)  THEN
552       time = time - 1
553       iday = iday + 1
554      ENDIF
555      CALL geopot  ( ip1jmp1, tpot  , pk , pks,  phis  , phi   )
556      print*,'sortie geopot'
557     
558      CALL caldyn0 ( itau,uvent,vvent,tpot,psol,masse,pk,phis ,
559     *                phi,w, pbaru,pbarv,time+iday-dayref   )
560       print*,'sortie caldyn0'     
561      CALL dynredem0("start.nc",dayref,anneeref,phis,nqmx)
562      print*,'sortie dynredem0'
563      CALL dynredem1("start.nc",0.0,vvent,uvent,tpot,q3d,nqmx,masse ,
564     .                            psol)
565      print*,'sortie dynredem1'
566C
567C Ecriture etat initial physique
568C
569      phystep   = dtvr * FLOAT(iphysiq)
570      radpas    = NINT (86400./phystep/ FLOAT(nbapp_rad) )
571      co2_ppm   = 330.0
572      solaire   = 1370.0
573
574c      call physdem(lonfi, latfi, phystep,radpas,co2_ppm,
575c     .                   solaire,tsol, qsol,
576c     .                   sn, radsol, deltat, rugmer,
577c     .                   agesno, zmea, zstd, zsig,
578c     .                   zgam, zthe, zpic, zval,
579c     .                   rugsrel)
580
581c
582c Initialisation
583c tsol, qsol, sn,albe, evap,tsoil,rain_fall, snow_fall,solsw, sollw,frugs
584c
585      tsolsrf(:,is_ter) = tsol
586      tsolsrf(:,is_lic) = tsol
587      tsolsrf(:,is_oce) = tsol
588      tsolsrf(:,is_sic) = tsol
589      snsrf(:,is_ter) = sn
590      snsrf(:,is_lic) = sn
591      snsrf(:,is_oce) = sn
592      snsrf(:,is_sic) = sn
593      albe(:,is_ter) = 0.08
594      albe(:,is_lic) = 0.6
595      albe(:,is_oce) = 0.5
596      albe(:,is_sic) = 0.6
597      evap(:,:) = 0.
598      qsolsrf(:,is_ter) = 150
599      qsolsrf(:,is_lic) = 150
600      qsolsrf(:,is_oce) = 150.
601      qsolsrf(:,is_sic) = 150.
602      do i = 1, nbsrf
603        do j = 1, nsoilmx
604          tsoil(:,j,i) = tsol
605        enddo
606      enddo
607      rain_fall = 0.; snow_fall = 0.
608      solsw = 165.
609      sollw = -53.
610      t_ancien = 273.15
611      q_ancien = 0.
612      agesno = 0.
613      deltat = 0.
614      frugs(1:klon,is_oce) = rugmer(1:klon)
615      frugs(1:klon,is_ter) = MAX(1.0e-05, zstd(1:klon)*zsig(1:klon)/2.0)
616      frugs(1:klon,is_lic) = MAX(1.0e-05, zstd(1:klon)*zsig(1:klon)/2.0)
617      frugs(1:klon,is_sic) = 0.001
618
619      call physdem("startphy.nc",phystep,radpas, co2_ppm, solaire,
620     $    latfi, lonfi, pctsrf, tsolsrf, tsoil, deltat, qsolsrf, snsrf,
621     $    albe, evap, rain_fall, snow_fall, solsw, sollw,
622     $    radsol, frugs,  agesno,
623     $    zmea, zstd, zsig, zgam, zthe, zpic, zval, rugsrel,
624     $    t_ancien, q_ancien)
625
626C     Sortie Visu pour les champs dynamiques
627      print*,'sortie visu'
628      time_step = 1.
629      t_ops = 2.
630      t_wrt = 2.
631      itau = 2.
632      visu_file='Etat0_visu.nc'
633      CALL initdynav(visu_file,dayref,anneeref,time_step,
634     .              t_ops, t_wrt, nqmx, visuid)
635      CALL writedynav(visuid, nqmx, itau,vvent ,
636     .                uvent,tpot,pk,phi,q3d,masse,psol,phis)
637      print*,'entree histclo'
638      CALL histclo
639      RETURN
640      !
641      END SUBROUTINE etat0_netcdf
Note: See TracBrowser for help on using the repository browser.