source: LMDZ.3.3/tags/IPSL-CM4_LJ29/libf/dyn3d/create_limit.F @ 494

Last change on this file since 494 was 494, checked in by (none), 20 years ago

This commit was manufactured by cvs2svn to create tag
'IPSL-CM4_LJ29'.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 32.7 KB
Line 
1c $Header$
2      PROGRAM create_limit
3      USE startvar
4      USE ioipsl
5      IMPLICIT none
6c
7c-------------------------------------------------------------
8C Author : L. Fairhead
9C Date   : 27/01/94
10C Objet  : Construction des fichiers de conditions aux limites
11C          pour le nouveau
12C          modele a partir de fichiers de climatologie. Les deux
13C          grilles doivent etre regulieres
14c
15c Modifie par z.x.li (le23mars1994)
16c Modifie par L. Fairhead (fairhead@lmd.jussieu.fr) septembre 1999
17c                         pour lecture netcdf dans LMDZ.3.3
18c modifie par P. Braconnot pour utiliser la version sous-surfaces
19c-------------------------------------------------------------
20c
21#include "dimensions.h"
22#include "paramet.h"
23#include "control.h"
24#include "logic.h"
25#include "netcdf.inc"
26#include "comvert.h"
27#include "comgeom2.h"
28#include "comconst.h"
29#include "dimphy.h"
30#include "indicesol.h"
31c-----------------------------------------------------------------------
32      REAL phy_nat(klon,360)
33      real phy_nat0(klon)
34      REAL phy_alb(klon,360)
35      REAL phy_sst(klon,360)
36      REAL phy_bil(klon,360)
37      REAL phy_rug(klon,360)
38      REAL phy_ice(klon)
39      logical interbar
40      parameter (interbar=.false.)
41CPB
42c      REAL phy_icet(klon,360)
43c      REAL phy_oce(klon,360)
44      real pctsrf_t(klon,nbsrf,360)
45      real pctsrf(klon,nbsrf)
46      REAL verif
47c
48      REAL masque(iip1,jjp1)
49      REAL mask(iim,jjp1)
50CPB
51C newlmt indique l'utilisation de la sous-maille fractionnelle
52C tandis que l'ancien codage utilise l'indicateur du sol (0,1,2,3)
53      LOGICAL newlmt, fracterre
54      PARAMETER(newlmt=.TRUE.)
55      PARAMETER(fracterre = .TRUE.)
56CPB
57C Declarations pour le champ de depart
58      INTEGER imdep, jmdep,lmdep
59      INTEGER ibid, jbid, tbid
60      PARAMETER (ibid = 400,       ! >360 pts
61     .           jbid = 200,       ! >181 pts
62     .           tbid = 60)        ! >52 semaines
63      REAL champ(ibid*jbid)
64      REAL dlon(ibid), dlat(jbid), timecoord(tbid)
65c
66      INTEGER ibid_msk, jbid_msk
67      PARAMETER(ibid_msk=2200,jbid_msk=1100)
68      REAL champ_msk(ibid_msk*jbid_msk)
69      REAL dlon_msk(ibid_msk), dlat_msk(jbid_msk)
70      REAL*4 zbidon(ibid_msk*jbid_msk)
71C Declarations pour le champ interpole 2D
72      REAL champint(iim,jjp1)
73
74C Declarations pour le champ interpole 3D
75      REAL champtime(iim,jjp1,tbid)
76      REAL timeyear(tbid)
77      REAL champan(iip1,jjp1,366)
78
79C Declarations pour l'inteprolation verticale
80      REAL ax(tbid), ay(tbid)
81      REAL by
82      REAL yder(tbid)
83
84
85      INTEGER ierr
86      INTEGER dimfirst(3)
87      INTEGER dimlast(3)
88c
89      INTEGER nid, ndim, ntim
90      INTEGER dims(2), debut(2), epais(2)
91      INTEGER id_tim
92      INTEGER id_NAT, id_SST, id_BILS, id_RUG, id_ALB
93CPB
94      INTEGER id_FOCE, id_FSIC, id_FTER, id_FLIC
95
96      INTEGER i, j, k, l, ji
97c declarations pour lecture glace de mer
98      INTEGER :: iml_lic, jml_lic, llm_tmp, ttm_tmp, iret
99      INTEGER :: itaul(1), fid
100      REAL :: lev(1), date, dt
101      REAL, ALLOCATABLE, DIMENSION(:,:) :: lon_lic, lat_lic
102      REAL, ALLOCATABLE, DIMENSION(:)  :: dlon_lic, dlat_lic
103      REAL, ALLOCATABLE, DIMENSION (:,:) :: fraclic
104      REAL :: flic_tmp(iip1, jjp1)
105c Diverses variables locales
106      REAL time
107! pour la lecture du fichier masque ocean
108      integer :: nid_o2a
109      logical :: couple = .false.
110      INTEGER :: iml_omask, jml_omask
111      REAL, ALLOCATABLE, DIMENSION(:,:) :: lon_omask, lat_omask
112      REAL, ALLOCATABLE, DIMENSION(:)  :: dlon_omask, dlat_omask
113      REAL, ALLOCATABLE, DIMENSION (:,:) :: ocemask, ocetmp
114      real, dimension(klon) :: ocemask_fi
115
116
117      INTEGER          longcles
118      PARAMETER      ( longcles = 20 )
119      REAL  clesphy0 ( longcles      )
120#include "serre.h"
121      INTEGER ncid,varid,ndimid(4),dimid
122      character*30 namedim
123      CHARACTER*80 :: varname
124
125c initialisations:
126!      OPEN (8,file='run.def',form='formatted')
127!      CALL defrun_new(8,.TRUE.,clesphy0)
128!      CLOSE(8)
129      CALL conf_gcm( 99, .TRUE. , clesphy0 )
130
131      pi     = 4. * ATAN(1.)
132      rad    = 6 371 229.
133      omeg   = 4.* ASIN(1.)/(24.*3600.)
134      g      = 9.8
135      daysec = 86400.
136      kappa  = 0.2857143
137      cpp    = 1004.70885
138      dtvr    = daysec/FLOAT(day_step)
139
140c
141ccc      CALL iniconst ( non  indispensable )
142
143      CALL inigeom
144c
145c
146C Traitement du relief au sol
147c
148      write(*,*) 'Fabrication masque'
149
150      varname = 'masque'
151      masque(:,:) = 0.0
152      CALL startget(varname, iip1, jjp1, rlonv, rlatu, masque, 0.0, jjm
153     ,     ,rlonu,rlatv , interbar )
154      pctsrf=0.
155      varname = 'zmasq'
156      zmasq(:) = 0.
157      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zmasq,0.0,
158     .jjm ,rlonu,rlatv , interbar)
159      WHERE (zmasq(1 : klon) .LT. EPSFRA)
160          zmasq(1 : klon) = 0.
161      END WHERE
162      WHERE (1 - zmasq(1 : klon) .LT. EPSFRA)
163          zmasq(1 : klon) = 1.
164      END WHERE
165!      WRITE(*,*)zmasq
166
167      IF ( fracterre ) THEN
168          DO i = 1, iim
169            masque(i,1) = masque(i,1)
170            masque(i,jjp1) = masque(i,jjp1)
171          END DO
172      ELSE
173          DO i = 1, iim
174            masque(i,1) = FLOAT(NINT(masque(i,1)))
175            masque(i,jjp1) = FLOAT(NINT(masque(i,jjp1)))
176          END DO
177      ENDIF
178c$$$      DO i = 1, iim
179c$$$      DO j = 1, jjp1
180c$$$         mask(i,j) = masque(i,j)
181c$$$      ENDDO
182c$$$      ENDDO
183c$$$      CALL gr_dyn_fi(1, iip1, jjp1, klon, masque, phy_nat0)
184      phy_nat0(1:klon) = zmasq(1:klon)
185      mask = 0.
186      DO j = 1, jjp1
187        DO i = 1, iim
188          IF ( masque(i,j) .GE. EPSFRA) mask (i,j) = 1
189        END DO
190      END DO 
191C
192C En cas de simulation couplee, lecture du masque ocean issu du modele ocean
193C utilise pour calculer les poids et pour assurer l'adequation entre les
194C fractions d'ocean vu par l'atmosphere et l'ocean
195C
196
197      write(*,*)'Essai de lecture masque ocean'
198      iret = nf_open("o2a.nc", NF_NOWRITE, nid_o2a)
199      if (iret .ne. 0) then
200        write(*,*)'ATTENTION!! pas de fichier o2a.nc trouve'
201        write(*,*)'Run force'
202      else
203        couple = .true.
204        iret = nf_close(nid_o2a)
205        call flininfo("o2a.nc", iml_omask, jml_omask, llm_tmp, ttm_tmp
206     $    , nid_o2a)
207        if (iml_omask /= iim .or. jml_omask /= jjp1) then
208          write(*,*)'Dimensions non compatibles pour masque ocean'
209          write(*,*)'iim = ',iim,' iml_omask = ',iml_omask
210          write(*,*)'jjp1 = ',jjp1,' jml_omask = ',jml_omask
211          stop
212        endif
213        ALLOCATE(lat_omask(iml_omask, jml_omask), stat=iret)
214        ALLOCATE(lon_omask(iml_omask, jml_omask), stat=iret)
215        ALLOCATE(dlon_omask(iml_omask), stat=iret)
216        ALLOCATE(dlat_omask(jml_omask), stat=iret)
217        ALLOCATE(ocemask(iml_omask, jml_omask), stat=iret)
218        ALLOCATE(ocetmp(iml_omask, jml_omask), stat=iret)
219        CALL flinopen("o2a.nc", .FALSE., iml_omask, jml_omask, llm_tmp
220     $    , lon_omask, lat_omask, lev, ttm_tmp, itaul, date, dt, fid)
221        CALL flinget(fid, 'OceMask', iml_omask, jml_omask, llm_tmp,
222     $      ttm_tmp, 1, 1, ocetmp)
223        CALL flinclo(fid)
224        dlon_omask(1 : iml_omask) = lon_omask(1 : iml_omask, 1)
225        dlat_omask(1 : jml_omask) = lat_omask(1 , 1 : jml_omask)
226        ocemask = ocetmp
227        if (dlat_omask(1) < dlat_omask(jml_omask)) then
228          do j = 1, jml_omask
229            ocemask(:,j) = ocetmp(:,jml_omask-j+1)
230          enddo
231        endif
232C
233C passage masque ocean a la grille physique
234C
235        ocemask_fi(1) = ocemask(1,1)
236        do j = 2, jjm
237          do i = 1, iim
238            ocemask_fi((j-2)*iim + i + 1) = ocemask(i,j)
239          enddo
240        enddo
241        ocemask_fi(klon) = ocemask(1,jjp1)
242        zmasq = 1. - ocemask_fi
243      endif
244
245
246C
247C lecture du fichier glace de terre pour fixer la fraction de terre
248C et de glace de terre
249C
250      CALL flininfo("landiceref.nc", iml_lic, jml_lic,llm_tmp, ttm_tmp
251     $    , fid)
252      ALLOCATE(lat_lic(iml_lic, jml_lic), stat=iret)
253      ALLOCATE(lon_lic(iml_lic, jml_lic), stat=iret)
254      ALLOCATE(dlon_lic(iml_lic), stat=iret)
255      ALLOCATE(dlat_lic(jml_lic), stat=iret)
256      ALLOCATE(fraclic(iml_lic, jml_lic), stat=iret)
257      CALL flinopen("landiceref.nc", .FALSE., iml_lic, jml_lic, llm_tmp
258     $    , lon_lic, lat_lic, lev, ttm_tmp, itaul, date, dt, fid)
259      CALL flinget(fid, 'landice', iml_lic, jml_lic, llm_tmp, ttm_tmp
260     $    , 1, 1, fraclic)
261      CALL flinclo(fid)
262C
263C interpolation sur la grille T du modele
264C
265      WRITE(*,*) 'dimensions de landice iml_lic, jml_lic : ',
266     $    iml_lic, jml_lic
267c
268C sil les coordonnees sont en degres, on les transforme
269C
270      IF( MAXVAL( lon_lic(:,:) ) .GT. 2.0 * asin(1.0) )  THEN
271          lon_lic(:,:) = lon_lic(:,:) * 2.0* ASIN(1.0) / 180.
272      ENDIF
273      IF( maxval( lat_lic(:,:) ) .GT. 2.0 * asin(1.0)) THEN
274          lat_lic(:,:) = lat_lic(:,:) * 2.0 * asin(1.0) / 180.
275      ENDIF
276
277      dlon_lic(1 : iml_lic) = lon_lic(1 : iml_lic, 1)
278      dlat_lic(1 : jml_lic) = lat_lic(1 , 1 : jml_lic)
279C
280      CALL grille_m(iml_lic, jml_lic, dlon_lic, dlat_lic, fraclic
281     $    ,iim, jjp1,
282     $    rlonv, rlatu, flic_tmp(1 : iim, 1 : jjp1))
283c$$$      flic_tmp(1 : iim, 1 : jjp1) = champint(1: iim, 1 : jjp1)
284      flic_tmp(iip1, 1 : jjp1) = flic_tmp(1 , 1 : jjp1)
285C
286C passage sur la grille physique
287C
288      CALL gr_dyn_fi(1, iip1, jjp1, klon, flic_tmp,
289     $    pctsrf(1:klon, is_lic))
290C adequation avec le maque terre/mer
291      WHERE (pctsrf(1 : klon, is_lic) .LT. EPSFRA )
292          pctsrf(1 : klon, is_lic) = 0.
293      END WHERE
294      WHERE (zmasq( 1 : klon) .LT. EPSFRA)
295          pctsrf(1 : klon, is_lic) = 0.
296      END WHERE
297      pctsrf(1 : klon, is_ter) = zmasq(1 : klon)
298      DO ji = 1, klon
299        IF (zmasq(ji) .GT. EPSFRA) THEN
300            IF ( pctsrf(ji, is_lic) .GE. zmasq(ji)) THEN
301                pctsrf(ji, is_lic) = zmasq(ji)
302                pctsrf(ji, is_ter) = 0.
303            ELSE
304                pctsrf(ji,is_ter) = zmasq(ji) - pctsrf(ji, is_lic)
305                IF (pctsrf(ji,is_ter) .LT. EPSFRA) THEN
306                    pctsrf(ji,is_ter) = 0.
307                    pctsrf(ji, is_lic) = zmasq(ji)
308                ENDIF
309            ENDIF
310        ENDIF
311      END DO
312c
313c
314C Traitement de la rugosite
315c
316      PRINT*, 'Traitement de la rugosite'
317      ierr = NF_OPEN('Rugos.nc', NF_NOWRITE, ncid)
318      if (ierr.ne.0) then
319        print *, NF_STRERROR(ierr)
320        STOP
321      ENDIF
322
323      ierr = NF_INQ_VARID(ncid,'RUGOS',varid)
324      if (ierr.ne.0) then
325        print *, NF_STRERROR(ierr)
326        STOP
327      ENDIF
328      ierr = NF_INQ_VARDIMID (ncid,varid,ndimid)
329      if (ierr.ne.0) then
330        print *, NF_STRERROR(ierr)
331        STOP
332      ENDIF
333      ierr = NF_INQ_DIM(ncid,ndimid(1), namedim, imdep)
334      if (ierr.ne.0) then
335        print *, NF_STRERROR(ierr)
336        STOP
337      ENDIF
338      print*,'variable ', namedim, 'dimension ', imdep
339      ierr = NF_INQ_VARID(ncid,namedim,dimid)
340      if (ierr.ne.0) then
341        print *, NF_STRERROR(ierr)
342        STOP
343      ENDIF
344#ifdef NC_DOUBLE
345      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon)
346#else
347      ierr = NF_GET_VAR_REAL(ncid,dimid,dlon)
348#endif
349      if (ierr.ne.0) then
350        print *, NF_STRERROR(ierr)
351        STOP
352      ENDIF
353      ierr = NF_INQ_DIM(ncid,ndimid(2), namedim, jmdep)
354      if (ierr.ne.0) then
355        print *, NF_STRERROR(ierr)
356        STOP
357      ENDIF
358      print*,'variable ', namedim, 'dimension ', jmdep
359      ierr = NF_INQ_VARID(ncid,namedim,dimid)
360      if (ierr.ne.0) then
361        print *, NF_STRERROR(ierr)
362        STOP
363      ENDIF
364#ifdef NC_DOUBLE
365      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat)
366#else
367      ierr = NF_GET_VAR_REAL(ncid,dimid,dlat)
368#endif
369      if (ierr.ne.0) then
370        print *, NF_STRERROR(ierr)
371        STOP
372      ENDIF
373      ierr = NF_INQ_DIM(ncid,ndimid(3), namedim, lmdep)
374      if (ierr.ne.0) then
375        print *, NF_STRERROR(ierr)
376        STOP
377      ENDIF
378      print*,'variable ', namedim, 'dimension ', lmdep
379      ierr = NF_INQ_VARID(ncid,namedim,dimid)
380      if (ierr.ne.0) then
381        print *, NF_STRERROR(ierr)
382        STOP
383      ENDIF
384#ifdef NC_DOUBLE
385      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord)
386#else
387      ierr = NF_GET_VAR_REAL(ncid,dimid,timecoord)
388#endif
389      if (ierr.ne.0) then
390        print *, NF_STRERROR(ierr)
391        STOP
392      ENDIF
393c
394      DO l = 1, lmdep
395         dimfirst(1) = 1
396         dimfirst(2) = 1
397         dimfirst(3) = l
398c
399         dimlast(1) = imdep
400         dimlast(2) = jmdep
401         dimlast(3) = 1
402c
403         PRINT*,'Lecture temporelle et int. horizontale ',l,timecoord(l)
404         print*,dimfirst,dimlast
405#ifdef NC_DOUBLE
406         ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ)
407#else
408         ierr = NF_GET_VARA_REAL(ncid,varid,dimfirst,dimlast,champ)
409#endif
410         if (ierr.ne.0) then
411           print *, NF_STRERROR(ierr)
412           STOP
413         ENDIF
414   
415         CALL rugosite(imdep, jmdep, dlon, dlat, champ,
416     .             iim, jjp1, rlonv, rlatu, champint, mask)
417         DO j = 1,jjp1
418         DO i = 1, iim
419            champtime (i,j,l) = champint(i,j)
420         ENDDO
421         ENDDO
422      ENDDO
423c      write(70,*)champtime
424c
425      DO l = 1, lmdep
426         timeyear(l) = timecoord(l)
427      ENDDO
428
429      PRINT 222, timeyear
430222   FORMAT(2x,' Time year ',10f6.1)
431c
432       
433      PRINT*, 'Interpolation temporelle dans l annee'
434
435
436      DO j = 1, jjp1
437      DO i = 1, iim
438          DO l = 1, lmdep
439            ax(l) = timeyear(l)
440            ay(l) = champtime (i,j,l)
441          ENDDO
442          CALL SPLINE(ax,ay,lmdep,1.e30,1.e30,yder)
443          DO k = 1, 360
444            time = FLOAT(k-1)
445            CALL SPLINT(ax,ay,yder,lmdep,time,by)
446            champan(i,j,k) = by
447          ENDDO
448      ENDDO
449      ENDDO
450      DO k = 1, 360
451      DO j = 1, jjp1
452         champan(iip1,j,k) = champan(1,j,k)
453      ENDDO
454      ENDDO
455c
456      DO k = 1, 360
457         CALL gr_dyn_fi(1,iip1,jjp1,klon,champan(1,1,k), phy_rug(1,k))
458      ENDDO
459c
460      ierr = NF_CLOSE(ncid)
461c
462c
463C Traitement de la glace oceanique
464c
465      PRINT*, 'Traitement de la glace oceanique'
466      ierr = NF_OPEN('AMIP.nc', NF_NOWRITE, ncid)
467      if (ierr.ne.0) then
468        print *, NF_STRERROR(ierr)
469        STOP
470      ENDIF
471
472      ierr = NF_INQ_VARID(ncid,'SEA_ICE',varid)
473      if (ierr.ne.0) then
474        print *, NF_STRERROR(ierr)
475        STOP
476      ENDIF
477      ierr = NF_INQ_VARDIMID (ncid,varid,ndimid)
478      if (ierr.ne.0) then
479        print *, NF_STRERROR(ierr)
480        STOP
481      ENDIF
482      ierr = NF_INQ_DIM(ncid,ndimid(1), namedim, imdep)
483      if (ierr.ne.0) then
484        print *, NF_STRERROR(ierr)
485        STOP
486      ENDIF
487      print*,'variable ', namedim, 'dimension ', imdep
488      ierr = NF_INQ_VARID(ncid,namedim,dimid)
489      if (ierr.ne.0) then
490        print *, NF_STRERROR(ierr)
491        STOP
492      ENDIF
493#ifdef NC_DOUBLE
494      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon)
495#else
496      ierr = NF_GET_VAR_REAL(ncid,dimid,dlon)
497#endif
498      if (ierr.ne.0) then
499        print *, NF_STRERROR(ierr)
500        STOP
501      ENDIF
502      ierr = NF_INQ_DIM(ncid,ndimid(2), namedim, jmdep)
503      if (ierr.ne.0) then
504        print *, NF_STRERROR(ierr)
505        STOP
506      ENDIF
507      print*,'variable ', namedim, jmdep
508      ierr = NF_INQ_VARID(ncid,namedim,dimid)
509      if (ierr.ne.0) then
510        print *, NF_STRERROR(ierr)
511        STOP
512      ENDIF
513#ifdef NC_DOUBLE
514      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat)
515#else
516      ierr = NF_GET_VAR_REAL(ncid,dimid,dlat)
517#endif
518      if (ierr.ne.0) then
519        print *, NF_STRERROR(ierr)
520        STOP
521      ENDIF
522      ierr = NF_INQ_DIM(ncid,ndimid(3), namedim, lmdep)
523      if (ierr.ne.0) then
524        print *, NF_STRERROR(ierr)
525        STOP
526      ENDIF
527      print*,'variable ', namedim, lmdep
528      ierr = NF_INQ_VARID(ncid,namedim,dimid)
529      if (ierr.ne.0) then
530        print *, NF_STRERROR(ierr)
531        STOP
532      ENDIF
533#ifdef NC_DOUBLE
534      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord)
535#else
536      ierr = NF_GET_VAR_REAL(ncid,dimid,timecoord)
537#endif
538      if (ierr.ne.0) then
539        print *, NF_STRERROR(ierr)
540        STOP
541      ENDIF
542c
543      DO l = 1, lmdep
544         dimfirst(1) = 1
545         dimfirst(2) = 1
546         dimfirst(3) = l
547c
548         dimlast(1) = imdep
549         dimlast(2) = jmdep
550         dimlast(3) = 1
551c
552         PRINT*,'Lecture temporelle et int. horizontale ',l,timecoord(l)
553#ifdef NC_DOUBLE
554         ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ)
555#else
556         ierr = NF_GET_VARA_REAL(ncid,varid,dimfirst,dimlast,champ)
557#endif
558         if (ierr.ne.0) then
559           print *, NF_STRERROR(ierr)
560           STOP
561         ENDIF
562
563         CALL sea_ice(imdep, jmdep, dlon, dlat, champ,
564     .             iim, jjp1, rlonv, rlatu, champint)
565         DO j = 1,jjp1
566         DO i = 1, iim
567            champtime (i,j,l) = champint(i,j)
568         ENDDO
569         ENDDO
570      ENDDO
571c
572      DO l = 1, lmdep
573         timeyear(l) = timecoord(l)
574      ENDDO
575      PRINT 222,  timeyear
576c
577      PRINT*, 'Interpolation temporelle'
578      DO j = 1, jjp1
579      DO i = 1, iim
580          DO l = 1, lmdep
581            ax(l) = timeyear(l)
582            ay(l) = champtime (i,j,l)
583          ENDDO
584          CALL SPLINE(ax,ay,lmdep,1.e30,1.e30,yder)
585          DO k = 1, 360
586            time = FLOAT(k-1)
587            CALL SPLINT(ax,ay,yder,lmdep,time,by)
588            champan(i,j,k) = by
589          ENDDO
590      ENDDO
591      ENDDO
592      DO k = 1, 360
593      DO j = 1, jjp1
594         champan(iip1, j, k) = champan(1, j, k)
595      ENDDO
596      ENDDO
597c
598c      WRITE(*,*) 'phy_nat'
599c     WRITE(*,'(72f4.1)') phy_nat0(1:klon)
600c
601      DO k = 1, 360
602        CALL gr_dyn_fi(1, iip1, jjp1, klon,
603     .      champan(1,1,k), phy_ice)
604        IF ( newlmt) THEN
605
606CPB  en attendant de mettre fraction de terre
607c
608          WHERE(phy_ice(1:klon) .GE. 1.) phy_ice(1 : klon) = 1.
609          WHERE(phy_ice(1:klon) .LT. EPSFRA) phy_ice(1 : klon) = 0.
610c
611          IF (fracterre ) THEN
612c            WRITE(*,*) 'passe dans cas fracterre'
613            pctsrf_t(:,is_ter,k) = pctsrf(:,is_ter)
614            pctsrf_t(:,is_lic,k) = pctsrf(:,is_lic)
615            pctsrf_t(1:klon,is_sic,k) =   phy_ice(1:klon)
616     $            - pctsrf_t(1:klon,is_lic,k)
617c§§ Il y a des cas ou il y a de la glace dans landiceref et pas dans AMIP
618            WHERE (pctsrf_t(1:klon,is_sic,k) .LE. 0)
619                pctsrf_t(1:klon,is_sic,k) = 0.
620            END WHERE
621            WHERE( 1. - zmasq(1:klon) .LT. EPSFRA)
622                pctsrf_t(1:klon,is_sic,k) = 0.
623                pctsrf_t(1:klon,is_oce,k) = 0.
624            END WHERE
625            DO i = 1, klon
626c$$              pctsrf_t(i,is_sic,k) = (1. - pctsrf_t(i,is_lic,k) -
627c$$     .                               pctsrf_t(i,is_ter,k)) * phy_ice(i)
628c$$              pctsrf_t(i,is_oce,k) = 1. - pctsrf_t(i,is_lic,k) -
629c$$     .                      pctsrf_t(i,is_ter,k) - pctsrf_t(i,is_sic,k)
630              IF ( 1. - zmasq(i) .GT. EPSFRA) THEN
631                  IF ( pctsrf_t(i,is_sic,k) .GE. 1 - zmasq(i)) THEN
632                      pctsrf_t(i,is_sic,k) = 1 - zmasq(i)
633                      pctsrf_t(i,is_oce,k) = 0.
634                  ELSE
635                      pctsrf_t(i,is_oce,k) = 1 - zmasq(i)
636     $                    - pctsrf_t(i,is_sic,k)
637                      IF (pctsrf_t(i,is_oce,k) .LT. EPSFRA) THEN
638                          pctsrf_t(i,is_oce,k) = 0.
639                          pctsrf_t(i,is_sic,k) = 1 - zmasq(i)
640                      ENDIF
641                  ENDIF
642              ENDIF 
643              if (pctsrf_t(i,is_oce,k) .lt. 0.) then
644                  WRITE(*,*) 'pb sous maille au point : i,k '
645     $                , i,k,pctsrf_t(:,is_oce,k)
646              ENDIF
647              IF ( abs( pctsrf_t(i, is_ter,k) + pctsrf_t(i, is_lic,k) +
648     $            pctsrf_t(i, is_oce,k) + pctsrf_t(i, is_sic,k)  - 1.)
649     $            .GT. EPSFRA) THEN
650                  WRITE(*,*) 'physiq : pb sous surface au point ', i,
651     $                pctsrf_t(i, 1 : nbsrf,k), phy_ice(i)
652              ENDIF
653            END DO
654        ELSE
655            DO i = 1, klon
656              pctsrf_t(i,is_ter,k) = pctsrf(i,is_ter)
657              IF (NINT(pctsrf(i,is_ter)).EQ.1 ) THEN
658                pctsrf_t(i,is_sic,k) = 0.
659                pctsrf_t(i,is_oce,k) = 0.                 
660                IF(phy_ice(i) .GE. 1.e-5) THEN
661                  pctsrf_t(i,is_lic,k) = phy_ice(i)
662                  pctsrf_t(i,is_ter,k) = pctsrf_t(i,is_ter,k)
663     .                                   - pctsrf_t(i,is_lic,k)
664                ELSE
665                  pctsrf_t(i,is_lic,k) = 0.
666                ENDIF
667              ELSE
668                pctsrf_t(i,is_lic,k) = 0.
669                IF(phy_ice(i) .GE. 1.e-5) THEN
670                  pctsrf_t(i,is_ter,k) = 0.
671                  pctsrf_t(i,is_sic,k) = phy_ice(i)
672                  pctsrf_t(i,is_oce,k) = 1. - pctsrf_t(i,is_sic,k)
673                ELSE
674                  pctsrf_t(i,is_sic,k) = 0.
675                  pctsrf_t(i,is_oce,k) = 1.                     
676                ENDIF
677              ENDIF
678              verif = pctsrf_t(i,is_ter,k) +
679     .                pctsrf_t(i,is_oce,k) +
680     .                pctsrf_t(i,is_sic,k) +
681     .                pctsrf_t(i,is_lic,k)
682              IF ( verif .LT. 1. - 1.e-5 .OR.
683     $             verif .GT. 1 + 1.e-5) THEN 
684                WRITE(*,*) 'pb sous maille au point : i,k,verif '
685     $                    , i,k,verif
686              ENDIF
687            END DO
688          ENDIF
689        ELSE 
690          DO i = 1, klon
691            phy_nat(i,k) = phy_nat0(i)
692            IF ( (phy_ice(i) - 0.5).GE.1.e-5 ) THEN
693              IF (NINT(phy_nat0(i)).EQ.0) THEN
694                phy_nat(i,k) = 3.0
695              ELSE
696                phy_nat(i,k) = 2.0
697              ENDIF
698            ENDIF
699          END DO
700        ENDIF
701      ENDDO
702c
703      ierr = NF_CLOSE(ncid)
704c
705c
706C Traitement de la sst
707c
708      PRINT*, 'Traitement de la sst'
709      ierr = NF_OPEN('AMIP.nc', NF_NOWRITE, ncid)
710      if (ierr.ne.0) then
711        print *, NF_STRERROR(ierr)
712        STOP
713      ENDIF
714
715      ierr = NF_INQ_VARID(ncid,'SST',varid)
716      if (ierr.ne.0) then
717        print *, NF_STRERROR(ierr)
718        STOP
719      ENDIF
720      ierr = NF_INQ_VARDIMID (ncid,varid,ndimid)
721      if (ierr.ne.0) then
722        print *, NF_STRERROR(ierr)
723        STOP
724      ENDIF
725      ierr = NF_INQ_DIM(ncid,ndimid(1), namedim, imdep)
726      if (ierr.ne.0) then
727        print *, NF_STRERROR(ierr)
728        STOP
729      ENDIF
730      print*,'variable ', namedim,'dimension ', imdep
731      ierr = NF_INQ_VARID(ncid,namedim,dimid)
732      if (ierr.ne.0) then
733        print *, NF_STRERROR(ierr)
734        STOP
735      ENDIF
736#ifdef NC_DOUBLE
737      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon)
738#else
739      ierr = NF_GET_VAR_REAL(ncid,dimid,dlon)
740#endif
741
742      if (ierr.ne.0) then
743        print *, NF_STRERROR(ierr)
744        STOP
745      ENDIF
746      ierr = NF_INQ_DIM(ncid,ndimid(2), namedim, jmdep)
747      if (ierr.ne.0) then
748        print *, NF_STRERROR(ierr)
749        STOP
750      ENDIF
751      print*,'variable ', namedim, 'dimension ', jmdep
752      ierr = NF_INQ_VARID(ncid,namedim,dimid)
753      if (ierr.ne.0) then
754        print *, NF_STRERROR(ierr)
755        STOP
756      ENDIF
757#ifdef NC_DOUBLE
758      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat)
759#else
760      ierr = NF_GET_VAR_REAL(ncid,dimid,dlat)
761#endif
762      if (ierr.ne.0) then
763        print *, NF_STRERROR(ierr)
764        STOP
765      ENDIF
766      ierr = NF_INQ_DIM(ncid,ndimid(3), namedim, lmdep)
767      if (ierr.ne.0) then
768        print *, NF_STRERROR(ierr)
769        STOP
770      ENDIF
771      print*,'variable ', namedim, 'dimension ', lmdep
772      ierr = NF_INQ_VARID(ncid,namedim,dimid)
773      if (ierr.ne.0) then
774        print *, NF_STRERROR(ierr)
775        STOP
776      ENDIF
777#ifdef NC_DOUBLE
778      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord)
779#else
780      ierr = NF_GET_VAR_REAL(ncid,dimid,timecoord)
781#endif
782      if (ierr.ne.0) then
783        print *, NF_STRERROR(ierr)
784        STOP
785      ENDIF
786c
787      DO l = 1, lmdep
788         dimfirst(1) = 1
789         dimfirst(2) = 1
790         dimfirst(3) = l
791c
792         dimlast(1) = imdep
793         dimlast(2) = jmdep
794         dimlast(3) = 1
795c
796         PRINT*,'Lecture temporelle et int. horizontale ',l,timecoord(l)
797#ifdef NC_DOUBLE
798         ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ)
799#else
800         ierr = NF_GET_VARA_REAL(ncid,varid,dimfirst,dimlast,champ)
801#endif
802         if (ierr.ne.0) then
803           print *, NF_STRERROR(ierr)
804           STOP
805         ENDIF
806         CALL grille_m(imdep, jmdep, dlon, dlat, champ,
807     .             iim, jjp1, rlonv, rlatu, champint)
808
809         DO j = 1,jjp1
810         DO i = 1, iim
811            champtime (i,j,l) = champint(i,j)
812         ENDDO
813         ENDDO
814      ENDDO
815c
816      DO l = 1, lmdep
817         timeyear(l) = timecoord(l)
818      ENDDO
819      print 222,  timeyear
820c
821C interpolation temporelle
822      DO j = 1, jjp1
823      DO i = 1, iim
824          DO l = 1, lmdep
825            ax(l) = timeyear(l)
826            ay(l) = champtime (i,j,l)
827          ENDDO
828          CALL SPLINE(ax,ay,lmdep,1.e30,1.e30,yder)
829          DO k = 1, 360
830            time = FLOAT(k-1)
831            CALL SPLINT(ax,ay,yder,lmdep,time,by)
832            champan(i,j,k) = by
833          ENDDO
834      ENDDO
835      ENDDO
836      DO k = 1, 360
837      DO j = 1, jjp1
838         champan(iip1,j,k) = champan(1,j,k)
839      ENDDO
840      ENDDO
841c
842      DO k = 1, 360
843         CALL gr_dyn_fi(1, iip1, jjp1, klon,
844     .                  champan(1,1,k), phy_sst(1,k))
845      ENDDO
846c
847      WHERE(phy_sst .LT. 271.35) phy_sst = 271.35
848      ierr = NF_CLOSE(ncid)
849c
850c
851C Traitement de l'albedo
852c
853      PRINT*, 'Traitement de l albedo'
854      ierr = NF_OPEN('Albedo.nc', NF_NOWRITE, ncid)
855      if (ierr.ne.0) then
856        print *, NF_STRERROR(ierr)
857        STOP
858      ENDIF
859      ierr = NF_INQ_VARID(ncid,'ALBEDO',varid)
860      if (ierr.ne.0) then
861        print *, NF_STRERROR(ierr)
862        STOP
863      ENDIF
864      ierr = NF_INQ_VARDIMID (ncid,varid,ndimid)
865      if (ierr.ne.0) then
866        print *, NF_STRERROR(ierr)
867        STOP
868      ENDIF
869      ierr = NF_INQ_DIM(ncid,ndimid(1), namedim, imdep)
870      if (ierr.ne.0) then
871        print *, NF_STRERROR(ierr)
872        STOP
873      ENDIF
874      print*,'variable ', namedim, 'dimension ', imdep
875      ierr = NF_INQ_VARID(ncid,namedim,dimid)
876      if (ierr.ne.0) then
877        print *, NF_STRERROR(ierr)
878        STOP
879      ENDIF
880#ifdef NC_DOUBLE
881      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon)
882#else
883      ierr = NF_GET_VAR_REAL(ncid,dimid,dlon)
884#endif
885      if (ierr.ne.0) then
886        print *, NF_STRERROR(ierr)
887        STOP
888      ENDIF
889      ierr = NF_INQ_DIM(ncid,ndimid(2), namedim, jmdep)
890      if (ierr.ne.0) then
891        print *, NF_STRERROR(ierr)
892        STOP
893      ENDIF
894      print*,'variable ', namedim, 'dimension ', jmdep
895      ierr = NF_INQ_VARID(ncid,namedim,dimid)
896      if (ierr.ne.0) then
897        print *, NF_STRERROR(ierr)
898        STOP
899      ENDIF
900#ifdef NC_DOUBLE
901      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat)
902#else
903      ierr = NF_GET_VAR_REAL(ncid,dimid,dlat)
904#endif
905      if (ierr.ne.0) then
906        print *, NF_STRERROR(ierr)
907        STOP
908      ENDIF
909      ierr = NF_INQ_DIM(ncid,ndimid(3), namedim, lmdep)
910      if (ierr.ne.0) then
911        print *, NF_STRERROR(ierr)
912        STOP
913      ENDIF
914      print*,'variable ', namedim, 'dimension ', lmdep
915      ierr = NF_INQ_VARID(ncid,namedim,dimid)
916      if (ierr.ne.0) then
917        print *, NF_STRERROR(ierr)
918        STOP
919      ENDIF
920#ifdef NC_DOUBLE
921      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord)
922#else
923      ierr = NF_GET_VAR_REAL(ncid,dimid,timecoord)
924#endif
925      if (ierr.ne.0) then
926        print *, NF_STRERROR(ierr)
927        STOP
928      ENDIF
929c
930      DO l = 1, lmdep
931         dimfirst(1) = 1
932         dimfirst(2) = 1
933         dimfirst(3) = l
934c
935         dimlast(1) = imdep
936         dimlast(2) = jmdep
937         dimlast(3) = 1
938c
939         PRINT*,'Lecture temporelle et int. horizontale ',l,timecoord(l)
940#ifdef NC_DOUBLE
941         ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ)
942#else
943         ierr = NF_GET_VARA_REAL(ncid,varid,dimfirst,dimlast,champ)
944#endif
945         if (ierr.ne.0) then
946           print *, NF_STRERROR(ierr)
947           STOP
948         ENDIF
949         CALL grille_m(imdep, jmdep, dlon, dlat, champ,
950     .             iim, jjp1, rlonv, rlatu, champint)
951c
952         DO j = 1,jjp1
953         DO i = 1, iim
954            champtime (i, j, l) = champint(i, j)
955         ENDDO
956         ENDDO
957      ENDDO
958c
959      DO l = 1, lmdep
960         timeyear(l) = timecoord(l)
961      ENDDO
962      print 222,  timeyear
963c
964C interpolation temporelle
965      DO j = 1, jjp1
966      DO i = 1, iim
967          DO l = 1, lmdep
968            ax(l) = timeyear(l)
969            ay(l) = champtime (i, j, l)
970          ENDDO
971          CALL SPLINE(ax,ay,lmdep,1.e30,1.e30,yder)
972          DO k = 1, 360
973            time = FLOAT(k-1)
974            CALL SPLINT(ax,ay,yder,lmdep,time,by)
975            champan(i,j,k) = by
976          ENDDO
977      ENDDO
978      ENDDO
979      DO k = 1, 360
980      DO j = 1, jjp1
981         champan(iip1, j, k) = champan(1, j, k)
982      ENDDO
983      ENDDO
984c
985      DO k = 1, 360
986         CALL gr_dyn_fi(1, iip1, jjp1, klon,
987     .                  champan(1,1,k), phy_alb(1,k))
988      ENDDO
989c
990      ierr = NF_CLOSE(ncid)
991c
992c
993      DO k = 1, 360
994      DO i = 1, klon
995         phy_bil(i,k) = 0.0
996      ENDDO
997      ENDDO
998c
999      PRINT*, 'Ecriture du fichier limit'
1000c
1001      ierr = NF_CREATE ("limit.nc", NF_CLOBBER, nid)
1002c
1003      ierr = NF_PUT_ATT_TEXT (nid, NF_GLOBAL, "title", 30,
1004     .                       "Fichier conditions aux limites")
1005      ierr = NF_DEF_DIM (nid, "points_physiques", klon, ndim)
1006      ierr = NF_DEF_DIM (nid, "time", NF_UNLIMITED, ntim)
1007c
1008      dims(1) = ndim
1009      dims(2) = ntim
1010c
1011      ierr = NF_DEF_VAR (nid, "TEMPS", NF_FLOAT, 1,ntim, id_tim)
1012      ierr = NF_PUT_ATT_TEXT (nid, id_tim, "title", 17,
1013     .                        "Jour dans l annee")
1014      IF (newlmt) THEN
1015c
1016          ierr = NF_DEF_VAR (nid, "FOCE", NF_FLOAT, 2,dims, id_FOCE)
1017          ierr = NF_PUT_ATT_TEXT (nid, id_FOCE, "title", 14,
1018     .                        "Fraction ocean")
1019c
1020          ierr = NF_DEF_VAR (nid, "FSIC", NF_FLOAT, 2,dims, id_FSIC)
1021          ierr = NF_PUT_ATT_TEXT (nid, id_FSIC, "title", 21,
1022     .                        "Fraction glace de mer")
1023c
1024          ierr = NF_DEF_VAR (nid, "FTER", NF_FLOAT, 2,dims, id_FTER)
1025          ierr = NF_PUT_ATT_TEXT (nid, id_FTER, "title", 14,
1026     .                        "Fraction terre")
1027c
1028          ierr = NF_DEF_VAR (nid, "FLIC", NF_FLOAT, 2,dims, id_FLIC)
1029          ierr = NF_PUT_ATT_TEXT (nid, id_FLIC, "title", 17,
1030     .                        "Fraction land ice")
1031c
1032      ELSE
1033          ierr = NF_DEF_VAR (nid, "NAT", NF_FLOAT, 2,dims, id_NAT)
1034          ierr = NF_PUT_ATT_TEXT (nid, id_NAT, "title", 23,
1035     .                        "Nature du sol (0,1,2,3)")
1036      ENDIF
1037C
1038      ierr = NF_DEF_VAR (nid, "SST", NF_FLOAT, 2,dims, id_SST)
1039      ierr = NF_PUT_ATT_TEXT (nid, id_SST, "title", 35,
1040     .                        "Temperature superficielle de la mer")
1041      ierr = NF_DEF_VAR (nid, "BILS", NF_FLOAT, 2,dims, id_BILS)
1042      ierr = NF_PUT_ATT_TEXT (nid, id_BILS, "title", 32,
1043     .                        "Reference flux de chaleur au sol")
1044      ierr = NF_DEF_VAR (nid, "ALB", NF_FLOAT, 2,dims, id_ALB)
1045      ierr = NF_PUT_ATT_TEXT (nid, id_ALB, "title", 19,
1046     .                        "Albedo a la surface")
1047      ierr = NF_DEF_VAR (nid, "RUG", NF_FLOAT, 2,dims, id_RUG)
1048      ierr = NF_PUT_ATT_TEXT (nid, id_RUG, "title", 8,
1049     .                        "Rugosite")
1050c
1051      ierr = NF_ENDDEF(nid)
1052c
1053      DO k = 1, 360
1054c
1055      debut(1) = 1
1056      debut(2) = k
1057      epais(1) = klon
1058      epais(2) = 1
1059c
1060#ifdef NC_DOUBLE
1061      ierr = NF_PUT_VAR1_DOUBLE (nid,id_tim,k,DBLE(k))
1062c
1063      IF (newlmt ) THEN
1064          ierr = NF_PUT_VARA_DOUBLE (nid,id_FOCE,debut,epais
1065     $        ,pctsrf_t(1,is_oce,k))
1066          ierr = NF_PUT_VARA_DOUBLE (nid,id_FSIC,debut,epais
1067     $        ,pctsrf_t(1,is_sic,k))
1068          ierr = NF_PUT_VARA_DOUBLE (nid,id_FTER,debut,epais
1069     $        ,pctsrf_t(1,is_ter,k))
1070          ierr = NF_PUT_VARA_DOUBLE (nid,id_FLIC,debut,epais
1071     $        ,pctsrf_t(1,is_lic,k))
1072      ELSE
1073          ierr = NF_PUT_VARA_DOUBLE (nid,id_NAT,debut,epais
1074     $        ,phy_nat(1,k))
1075      ENDIF
1076c
1077      ierr = NF_PUT_VARA_DOUBLE (nid,id_SST,debut,epais,phy_sst(1,k))
1078      ierr = NF_PUT_VARA_DOUBLE (nid,id_BILS,debut,epais,phy_bil(1,k))
1079      ierr = NF_PUT_VARA_DOUBLE (nid,id_ALB,debut,epais,phy_alb(1,k))
1080      ierr = NF_PUT_VARA_DOUBLE (nid,id_RUG,debut,epais,phy_rug(1,k))
1081#else
1082      ierr = NF_PUT_VAR1_REAL (nid,id_tim,k,FLOAT(k))
1083      IF (newlmt ) THEN
1084          ierr = NF_PUT_VARA_REAL (nid,id_FOCE,debut,epais
1085     $        ,pctsrf_t(1,is_oce,k))
1086          ierr = NF_PUT_VARA_REAL (nid,id_FSIC,debut,epais
1087     $        ,pctsrf_t(1,is_sic,k))
1088          ierr = NF_PUT_VARA_REAL (nid,id_FTER,debut,epais
1089     $        ,pctsrf_t(1,is_ter,k))
1090          ierr = NF_PUT_VARA_REAL (nid,id_FLIC,debut,epais
1091     $        ,pctsrf_t(1,is_lic,k))
1092      ELSE
1093          ierr = NF_PUT_VARA_REAL (nid,id_NAT,debut,epais
1094     $        ,phy_nat(1,k))
1095      ENDIF
1096      ierr = NF_PUT_VARA_REAL (nid,id_SST,debut,epais,phy_sst(1,k))
1097      ierr = NF_PUT_VARA_REAL (nid,id_BILS,debut,epais,phy_bil(1,k))
1098      ierr = NF_PUT_VARA_REAL (nid,id_ALB,debut,epais,phy_alb(1,k))
1099      ierr = NF_PUT_VARA_REAL (nid,id_RUG,debut,epais,phy_rug(1,k))
1100#endif
1101c
1102      ENDDO
1103c
1104      ierr = NF_CLOSE(nid)
1105c
1106      STOP
1107      END
1108
Note: See TracBrowser for help on using the repository browser.