source: LMDZ4/branches/pre_V3/libf/dyn3dpar/etat0_netcdf.F @ 3837

Last change on this file since 3837 was 630, checked in by Laurent Fairhead, 20 years ago

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