source: LMDZ.3.3/branches/rel-LF/libf/dyn3d/create_limit.F @ 279

Last change on this file since 279 was 278, checked in by lmdzadmin, 23 years ago

Remplacement defrun_new part conf_gcm
LF

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