source: LMDZ.3.3/branches/rel-LF/libf/dyn3d/etat0_netcdf.F @ 316

Last change on this file since 316 was 278, checked in by lmdzadmin, 23 years ago

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