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

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

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