source: LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3d/etat0_netcdf.F @ 2482

Last change on this file since 2482 was 1379, checked in by lguez, 15 years ago

Added optional ozone tracer with chemistry parameterized by Daniel
Cariolle. This tracer is passive: it has no influence on the rest of
the simulation.

Added variable "zmasse" in file "histrac.nc". Corrected long name of
variable "pplay" in "histrac.nc". Changed name of variable "t" to "T"
in "histrac.nc", corrected long name and unit.

In "phytrac", moved definition of "zmasse" toward the beginning of the
procedure, so that "zmasse" can be given as argument to
"traclmdz". Also added arguments "julien", "gmtime" and "xlon" to
"traclmdz". The four additional arguments are required for the ozone
tracer.

In module "traclmdz_mod", factorized declaration "implicit none" that
was in each procedure. There is now an equivalent single declaration
at the module level.

In procedure "traclmdz", removed variable "delp". Use "zmasse * rg"
instead since we now have "zmasse" as an argument.

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