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

Last change on this file since 332 was 332, checked in by lmdzadmin, 22 years ago

Menage
LF
  • 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      endif
421
422
423C
424C lecture du fichier glace de terre pour fixer la fraction de terre
425C et de glace de terre
426C
427      CALL flininfo("landiceref.nc", iml_lic, jml_lic,llm_tmp, ttm_tmp
428     $    , fid)
429      ALLOCATE(lat_lic(iml_lic, jml_lic), stat=iret)
430      ALLOCATE(lon_lic(iml_lic, jml_lic), stat=iret)
431      ALLOCATE(dlon_lic(iml_lic), stat=iret)
432      ALLOCATE(dlat_lic(jml_lic), stat=iret)
433      ALLOCATE(fraclic(iml_lic, jml_lic), stat=iret)
434      CALL flinopen("landiceref.nc", .FALSE., iml_lic, jml_lic, llm_tmp
435     $    , lon_lic, lat_lic, lev, ttm_tmp, itaul, date, dt, fid)
436      CALL flinget(fid, 'landice', iml_lic, jml_lic, llm_tmp, ttm_tmp
437     $    , 1, 1, fraclic)
438      CALL flinclo(fid)
439C
440C interpolation sur la grille T du modele
441C
442      WRITE(*,*) 'dimensions de landice iml_lic, jml_lic : ',
443     $    iml_lic, jml_lic
444c
445C sil les coordonnees sont en degres, on les transforme
446C
447      IF( MAXVAL( lon_lic(:,:) ) .GT. 2.0 * asin(1.0) )  THEN
448          lon_lic(:,:) = lon_lic(:,:) * 2.0* ASIN(1.0) / 180.
449      ENDIF
450      IF( maxval( lat_lic(:,:) ) .GT. 2.0 * asin(1.0)) THEN
451          lat_lic(:,:) = lat_lic(:,:) * 2.0 * asin(1.0) / 180.
452      ENDIF
453
454      dlon_lic(1 : iml_lic) = lon_lic(1 : iml_lic, 1)
455      dlat_lic(1 : jml_lic) = lat_lic(1 , 1 : jml_lic)
456C
457      CALL grille_m(iml_lic, jml_lic, dlon_lic, dlat_lic, fraclic
458     $    ,iim, jjp1,
459     $    rlonv, rlatu, flic_tmp(1 : iim, 1 : jjp1))
460c$$$      flic_tmp(1 : iim, 1 : jjp1) = champint(1: iim, 1 : jjp1)
461      flic_tmp(iip1, 1 : jjp1) = flic_tmp(1 , 1 : jjp1)
462C
463C passage sur la grille physique
464C
465      CALL gr_dyn_fi(1, iip1, jjp1, klon, flic_tmp,
466     $    pctsrf(1:klon, is_lic))
467C adequation avec le maque terre/mer
468c      zmasq(157) = 0.
469      WHERE (pctsrf(1 : klon, is_lic) .LT. EPSFRA )
470          pctsrf(1 : klon, is_lic) = 0.
471      END WHERE
472      WHERE (zmasq( 1 : klon) .LT. EPSFRA)
473          pctsrf(1 : klon, is_lic) = 0.
474      END WHERE
475      pctsrf(1 : klon, is_ter) = zmasq(1 : klon)
476      DO ji = 1, klon
477        IF (zmasq(ji) .GT. EPSFRA) THEN
478            IF ( pctsrf(ji, is_lic) .GE. zmasq(ji)) THEN
479                pctsrf(ji, is_lic) = zmasq(ji)
480                pctsrf(ji, is_ter) = 0.
481            ELSE
482                pctsrf(ji,is_ter) = zmasq(ji) - pctsrf(ji, is_lic)
483                IF (pctsrf(ji,is_ter) .LT. EPSFRA) THEN
484                    pctsrf(ji,is_ter) = 0.
485                    pctsrf(ji, is_lic) = zmasq(ji)
486                ENDIF
487            ENDIF
488        ENDIF
489      END DO
490C
491C sous surface ocean et glace de mer (pour demarrer on met glace de mer a 0)
492C
493      pctsrf(1 : klon, is_oce) = (1. - zmasq(1 : klon))
494
495
496      WHERE (pctsrf(1 : klon, is_oce) .LT. EPSFRA)
497          pctsrf(1 : klon, is_oce) = 0.
498      END WHERE
499
500      if (couple) pctsrf(1 : klon, is_oce) = ocemask_fi(1 : klon)
501
502      isst = 0
503      where (pctsrf(2:klon-1,is_oce) >0.) isst = 1
504C
505C verif que somme des sous surface = 1
506C
507      ji=count( (abs( sum(pctsrf(1 : klon, 1 : nbsrf), dim = 2)) - 1.0 )
508     $    .GT. EPSFRA)
509      IF (ji .NE. 0) THEN
510          WRITE(*,*) 'pb repartition sous maille pour ',ji,' points'
511      ENDIF
512
513
514
515
516
517C Calcul intermediaire
518c
519      CALL massdair( p3d, masse  )
520c
521
522      print *,' ALPHAX ',alphax
523
524      DO  l = 1, llm
525        DO  i    = 1, iim
526          xppn(i) = aire( i, 1   ) * masse(  i     ,  1   , l )
527          xpps(i) = aire( i,jjp1 ) * masse(  i     , jjp1 , l )
528        ENDDO
529          xpn      = SUM(xppn)/apoln
530          xps      = SUM(xpps)/apols
531        DO i   = 1, iip1
532          masse(   i   ,   1     ,  l )   = xpn
533          masse(   i   ,   jjp1  ,  l )   = xps
534        ENDDO
535      ENDDO
536      q3d(iip1,:,:,:) = q3d(1,:,:,:)
537      phis(iip1,:) = phis(1,:)
538
539C Ecriture
540
541
542      CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   ,
543     *                tetagdiv, tetagrot , tetatemp              )
544      print*,'sortie inidissip'
545      itau = 0
546      iday = dayref +itau/day_step
547      time = FLOAT(itau-(iday-dayref)*day_step)/day_step
548c     
549      IF(time.GT.1)  THEN
550       time = time - 1
551       iday = iday + 1
552      ENDIF
553      CALL geopot  ( ip1jmp1, tpot  , pk , pks,  phis  , phi   )
554      print*,'sortie geopot'
555     
556      CALL caldyn0 ( itau,uvent,vvent,tpot,psol,masse,pk,phis ,
557     *                phi,w, pbaru,pbarv,time+iday-dayref   )
558       print*,'sortie caldyn0'     
559      CALL dynredem0("start.nc",dayref,anneeref,phis,nqmx)
560      print*,'sortie dynredem0'
561      CALL dynredem1("start.nc",0.0,vvent,uvent,tpot,q3d,nqmx,masse ,
562     .                            psol)
563      print*,'sortie dynredem1'
564C
565C Ecriture etat initial physique
566C
567      phystep   = dtvr * FLOAT(iphysiq)
568      radpas    = NINT (86400./phystep/ FLOAT(nbapp_rad) )
569      co2_ppm   = 330.0
570      solaire   = 1370.0
571
572c      call physdem(lonfi, latfi, phystep,radpas,co2_ppm,
573c     .                   solaire,tsol, qsol,
574c     .                   sn, radsol, deltat, rugmer,
575c     .                   agesno, zmea, zstd, zsig,
576c     .                   zgam, zthe, zpic, zval,
577c     .                   rugsrel)
578
579c
580c Initialisation
581c tsol, qsol, sn,albe, evap,tsoil,rain_fall, snow_fall,solsw, sollw,frugs
582c
583      tsolsrf(:,is_ter) = tsol
584      tsolsrf(:,is_lic) = tsol
585      tsolsrf(:,is_oce) = tsol
586      tsolsrf(:,is_sic) = tsol
587      snsrf(:,is_ter) = sn
588      snsrf(:,is_lic) = sn
589      snsrf(:,is_oce) = sn
590      snsrf(:,is_sic) = sn
591      albe(:,is_ter) = 0.08
592      albe(:,is_lic) = 0.6
593      albe(:,is_oce) = 0.5
594      albe(:,is_sic) = 0.6
595      evap(:,:) = 0.
596      qsolsrf(:,is_ter) = 150
597      qsolsrf(:,is_lic) = 150
598      qsolsrf(:,is_oce) = 150.
599      qsolsrf(:,is_sic) = 150.
600      do i = 1, nbsrf
601        do j = 1, nsoilmx
602          tsoil(:,j,i) = tsol
603        enddo
604      enddo
605      rain_fall = 0.; snow_fall = 0.
606      solsw = 165.
607      sollw = -53.
608      t_ancien = 273.15
609      q_ancien = 0.
610      agesno = 0.
611      deltat = 0.
612      frugs(1:klon,is_oce) = rugmer(1:klon)
613      frugs(1:klon,is_ter) = MAX(1.0e-05, zstd(1:klon)*zsig(1:klon)/2.0)
614      frugs(1:klon,is_lic) = MAX(1.0e-05, zstd(1:klon)*zsig(1:klon)/2.0)
615      frugs(1:klon,is_sic) = 0.001
616
617      call physdem("startphy.nc",phystep,radpas, co2_ppm, solaire,
618     $    latfi, lonfi, pctsrf, tsolsrf, tsoil, deltat, qsolsrf, snsrf,
619     $    albe, evap, rain_fall, snow_fall, solsw, sollw,
620     $    radsol, frugs,  agesno,
621     $    zmea, zstd, zsig, zgam, zthe, zpic, zval, rugsrel,
622     $    t_ancien, q_ancien)
623
624C     Sortie Visu pour les champs dynamiques
625      print*,'sortie visu'
626      time_step = 1.
627      t_ops = 2.
628      t_wrt = 2.
629      itau = 2.
630      visu_file='Etat0_visu.nc'
631      CALL initdynav(visu_file,dayref,anneeref,time_step,
632     .              t_ops, t_wrt, nqmx, visuid)
633      CALL writedynav(visuid, nqmx, itau,vvent ,
634     .                uvent,tpot,pk,phi,q3d,masse,psol,phis)
635      print*,'entree histclo'
636      CALL histclo
637      RETURN
638      !
639      END SUBROUTINE etat0_netcdf
Note: See TracBrowser for help on using the repository browser.