source: LMDZ4/trunk/libf/dyn3dpar/etat0_netcdf.F @ 1076

Last change on this file since 1076 was 1058, checked in by lmdzadmin, 16 years ago

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