source: LMDZ.3.3/trunk/libf/dyn3d/limit_netcdf.F @ 2793

Last change on this file since 2793 was 425, checked in by lmdz, 22 years ago

Dans le cas "oldice" la routine "sea_ice" utilise seulement l'interpol.
grille_m Levan

LF

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