source: LMDZ4/trunk/libf/dyn3d/etat0_netcdf.F @ 724

Last change on this file since 724 was 724, checked in by Laurent Fairhead, 18 years ago

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