Subversion Repositories pvinversion.ecmwf

Rev

Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
3 michaesp 1
      PROGRAM coastline
2
 
3
c     *********************************************************************
4
c     * Add the coastline and land sea mask to the REF file               *
5
c     * Michael Sprenger / Winter 2006,07                                 *
6
c     *********************************************************************
7
 
8
      implicit none
9
 
10
c     --------------------------------------------------------------------
11
c     Declaration of parameters and variables
12
c     --------------------------------------------------------------------
13
 
14
c     Input parameters
15
      character*80 coastfile
16
      character*80 reffile
17
      real         clon,clat,crot
18
 
19
c     Output arrays
20
      integer      nmax
21
      parameter    (nmax=100000)
22
      real,allocatable, dimension (:,:)   :: landsea
23
      real,allocatable, dimension (:,:)   :: xcor,ycor
24
      real         xpos(nmax),ypos(nmax)
25
      real         lon(nmax), lat(nmax)
26
      real         rlon(nmax),rlat(nmax)
27
      integer      ncoast
28
      integer      rnx,rny
29
      real         rxmin,rymin,rxmax,rymax
30
      real         rdx,rdy
31
 
32
c     Auxiliary variables
33
      integer      i,j
34
      character*80 name
35
      integer      cdfid,cstid
36
      character*80 cfn
37
      integer      ierr,stat
38
      character*80 varname
39
      real         varmin(4),varmax(4),stag(4)
40
      integer      vardim(4),ndim  
41
      real         mdv
42
      integer      isok
43
      integer      nvars
44
      integer      ntimes
45
      real         time
46
      character*80 vnam(100)
47
 
48
c     --------------------------------------------------------------------
49
c     Preparations
50
c     --------------------------------------------------------------------
51
 
52
      print*,'*********************************************************'
53
      print*,'* coastline                                             *'
54
      print*,'*********************************************************'
55
 
56
 
57
c     Read parameter file
58
      open(10,file='fort.10')
59
       read(10,*) reffile
60
       read(10,*) coastfile
61
       read(10,*) name,clon
62
       if ( trim(name).ne.'CLON' ) stop
63
       read(10,*) name,clat
64
       if ( trim(name).ne.'CLAT' ) stop      
65
       read(10,*) name,crot
66
       if ( trim(name).ne.'CROT' ) stop    
67
      close(10)
68
 
69
      print*
70
      print*,trim(reffile)
71
      print*,trim(coastfile)
72
      print*,clon,clat,crot
73
 
74
c     Read grid parameters
75
      call cdfopn(reffile,cdfid,ierr)
76
      if (ierr.ne.0) goto 998
77
      call getvars(cdfid,nvars,vnam,ierr)
78
      if (ierr.ne.0) goto 998     
79
      varname='X'      
80
      call getdef(cdfid,varname,ndim,mdv,vardim,
81
     >            varmin,varmax,stag,ierr)
82
      if (ierr.ne.0) goto 998
83
      call gettimes(cdfid,time,ntimes,ierr)
84
      if (ierr.ne.0) goto 998
85
      call clscdf(cdfid,ierr)
86
      rnx=vardim(1)
87
      rny=vardim(2)
88
      rxmin=varmin(1)
89
      rymin=varmin(2)
90
      rxmax=varmax(1)
91
      rymax=varmax(2)
92
      rdx=(rxmax-rxmin)/real(rnx-1)
93
      rdy=(rymax-rymin)/real(rny-1)
94
 
95
c     Allocate memory
96
      allocate(landsea(rnx,rny),stat=stat)
97
      if (stat.ne.0) print*,'*** error allocating array landsea ***'
98
      allocate(xcor(rnx,rny),stat=stat)
99
      if (stat.ne.0) print*,'*** error allocating array xcor ***'
100
      allocate(ycor(rnx,rny),stat=stat)
101
      if (stat.ne.0) print*,'*** error allocating array ycor ***'
102
 
103
c     Read data
104
      call cdfopn(reffile,cdfid,ierr)
105
      if (ierr.ne.0) goto 998
106
      varname='X'
107
      call getdat(cdfid,varname,time,1,xcor,ierr)
108
      if (ierr.ne.0) goto 998
109
      varname='Y'
110
      call getdat(cdfid,varname,time,1,ycor,ierr)
111
      if (ierr.ne.0) goto 998
112
      call clscdf(cdfid,ierr)
113
 
114
c     Read the coastfile
115
      open(10,file=coastfile)
116
 
117
       ncoast=1
118
 100   read(10,*,end=200) lon(ncoast),lat(ncoast)
119
       ncoast=ncoast+1
120
       goto 100
121
 
122
 200  close(10)
123
      ncoast=ncoast-1
124
 
125
c     Handle missing data values
126
      mdv=-100000.
127
      do i=ncoast,nmax
128
         lon(i)=mdv
129
         lat(i)=mdv
130
      enddo
131
 
132
c     --------------------------------------------------------------------
133
c     Transform coast(lat/lon) to coast(x/y)
134
c     --------------------------------------------------------------------
135
 
136
      call  getenvir (clon,clat,crot,
137
     >                xcor,ycor,rnx,rny,rxmin,rymin,rdx,rdy,mdv,
138
     >                lon,lat,rlon,rlat,xpos,ypos,nmax)
139
 
140
c     --------------------------------------------------------------------
141
c     Write output 
142
c     --------------------------------------------------------------------
143
 
144
c     Open output file and set some general parameters
145
      call cdfwopn(reffile,cdfid,ierr)
146
      if (ierr.ne.0) goto 998
147
      vardim(1)=1
148
      vardim(2)=1
149
      vardim(3)=nmax
150
      ndim=3
151
 
152
c     Write LAT
153
      isok=0
154
      varname='COAST_LAT'
155
      call check_varok(isok,varname,vnam,nvars)
156
      if (isok.eq.0) then
157
         call putdef(cdfid,varname,ndim,mdv,vardim,
158
     >               varmin,varmax,stag,ierr)
159
         if (ierr.ne.0) goto 998
160
      endif
161
      call putdat(cdfid,varname,time,0,lat,ierr)
162
      if (ierr.ne.0) goto 998
163
      print*,'W ',trim(varname),' ',trim(reffile)
164
 
165
c     Write LON
166
      isok=0
167
      varname='COAST_LON'
168
      call check_varok(isok,varname,vnam,nvars)
169
      if (isok.eq.0) then
170
         call putdef(cdfid,varname,ndim,mdv,vardim,
171
     >               varmin,varmax,stag,ierr)
172
         if (ierr.ne.0) goto 998
173
      endif
174
      call putdat(cdfid,varname,time,0,lon,ierr)
175
      if (ierr.ne.0) goto 998
176
      print*,'W ',trim(varname),' ',trim(reffile)
177
 
178
c     Write RLAT
179
      isok=0
180
      varname='COAST_RLAT'
181
      call check_varok(isok,varname,vnam,nvars)
182
      if (isok.eq.0) then
183
         call putdef(cdfid,varname,ndim,mdv,vardim,
184
     >               varmin,varmax,stag,ierr)
185
         if (ierr.ne.0) goto 998
186
      endif
187
      call putdat(cdfid,varname,time,0,rlat,ierr)
188
      if (ierr.ne.0) goto 998
189
      print*,'W ',trim(varname),' ',trim(reffile)
190
 
191
c     Write RLON
192
      isok=0
193
      varname='COAST_RLON'
194
      call check_varok(isok,varname,vnam,nvars)
195
      if (isok.eq.0) then
196
         call putdef(cdfid,varname,ndim,mdv,vardim,
197
     >               varmin,varmax,stag,ierr)
198
         if (ierr.ne.0) goto 998
199
      endif
200
      call putdat(cdfid,varname,time,0,rlon,ierr)
201
      if (ierr.ne.0) goto 998
202
      print*,'W ',trim(varname),' ',trim(reffile)
203
 
204
c     Write Y
205
      isok=0
206
      varname='COAST_Y'
207
      call check_varok(isok,varname,vnam,nvars)
208
      if (isok.eq.0) then
209
         call putdef(cdfid,varname,ndim,mdv,vardim,
210
     >               varmin,varmax,stag,ierr)
211
         if (ierr.ne.0) goto 998
212
      endif
213
      call putdat(cdfid,varname,time,0,ypos,ierr)
214
      if (ierr.ne.0) goto 998
215
      print*,'W ',trim(varname),' ',trim(reffile)
216
 
217
c     Write X
218
      isok=0
219
      varname='COAST_X'
220
      call check_varok(isok,varname,vnam,nvars)
221
      if (isok.eq.0) then
222
         call putdef(cdfid,varname,ndim,mdv,vardim,
223
     >               varmin,varmax,stag,ierr)
224
         if (ierr.ne.0) goto 998
225
      endif
226
      call putdat(cdfid,varname,time,0,xpos,ierr)
227
      if (ierr.ne.0) goto 998
228
      print*,'W ',trim(varname),' ',trim(reffile)
229
 
230
c     Close output file
231
      call clscdf(cdfid,ierr)
232
      if (ierr.ne.0) goto 998
233
 
234
 
235
c     --------------------------------------------------------------------
236
c     Exception handling
237
c     --------------------------------------------------------------------
238
 
239
      stop
240
 
241
 998  print*,'Problems with input netcdf file'
242
      stop
243
 
244
      END
245
 
246
 
247
c     --------------------------------------------------------------------------------
248
c     Subroutine to get environment of strcof
249
c     --------------------------------------------------------------------------------
250
 
251
      SUBROUTINE getenvir (clon,clat,rotation,
252
     >                     xcor,ycor,rnx,rny,rxmin,rymin,rdx,rdy,mdv,
253
     >                     lon,lat,rlon,rlat,xpos,ypos,n)
254
 
255
c     Rotate from a local quasi-cartesian coordiante system into lat/lon coordinate 
256
c     system.
257
 
258
      implicit none
259
 
260
c     Declaration of input and output parameters
261
      integer     n
262
      real        clon,clat,rotation
263
      real        lon(n), lat(n)
264
      real        rlon(n),rlat(n)
265
      real        xpos(n),ypos(n)
266
      integer     rnx,rny
267
      real        xcor(rnx,rny),ycor(rnx,rny)
268
      real        rxmin,rymin,rdx,rdy
269
      real        mdv
270
 
271
c     Set numerical and physical constants
272
      real	  g2r
273
      parameter   (g2r=0.0174533)
274
      real        pi180
275
      parameter   (pi180=3.14159265359/180.)
276
      real        eps
277
      parameter   (eps=0.0001)
278
 
279
c     Auxiliary variables 
280
      real         pollon,pollat
281
      integer      i,j
282
      real         xind,yind
283
      real         rindx,rindy
284
      integer      indx,indy,indr,indu
285
      real         frac0i,frac0j,frac1i,frac1j
286
      real         rlon1(n),rlat1(n)
287
 
288
c     Externals
289
      real     lmtolms,phtophs
290
      external lmtolms,phtophs
291
 
292
c     First rotation
293
      pollon=clon-180.
294
      if (pollon.lt.-180.) pollon=pollon+360.
295
      pollat=90.-clat
296
      do i=1,n
297
 
298
         if ( (abs(lon(i)-mdv).gt.eps).and.
299
     >        (abs(lat(i)-mdv).gt.eps) ) then
300
 
301
c           First rotation
302
            pollon=clon-180.
303
            if (pollon.lt.-180.) pollon=pollon+360.
304
            pollat=90.-clat
305
            rlon1(i)=lmtolms(lat(i),lon(i),pollat,pollon)
306
            rlat1(i)=phtophs(lat(i),lon(i),pollat,pollon)            
307
 
308
c           Second coordinate transformation
309
            pollon=-180.
310
            pollat=90.+rotation
311
            rlon(i)=90.+lmtolms(rlat1(i),rlon1(i)-90.,pollat,pollon)
312
            rlat(i)=phtophs(rlat1(i),rlon1(i)-90.,pollat,pollon)   
313
 
314
c           Get rotated latitude and longitude
315
 100        if (rlon(i).lt.rxmin) then
316
               rlon(i)=rlon(i)+360.
317
               goto 100
318
            endif
319
 102        if (rlon(i).gt.(rxmin+real(rnx-1)*rdx)) then
320
               rlon(i)=rlon(i)-360.
321
               goto 102
322
            endif
323
            if ( (rlon(i).lt.rxmin).or.
324
     >           (rlon(i).gt.(rxmin+real(rnx-1)*rdx)).or.
325
     >           (rlat(i).lt.rymin).or.
326
     >           (rlat(i).gt.(rymin+real(rny-1)*rdy)) ) then
327
 
328
               rlon(i)=mdv
329
               rlat(i)=mdv
330
            endif
331
 
332
c           Get x and y coordinates
333
            rindx=(rlon(i)-rxmin)/rdx+1.
334
            rindy=(rlat(i)-rymin)/rdy+1.
335
            indx=int(rindx)
336
            indy=int(rindy)
337
            if ((indx.gt.1).and.(indx.lt.rnx).and.
338
     >          (indy.gt.1).and.(indy.lt.rny)) then
339
 
340
               indr=indx+1
341
               if (indr.gt.rnx) indr=1
342
               indu=indy+1
343
               if (indu.gt.rny) indu=rny
344
 
345
               frac0i=rindx-float(indx)
346
               frac0j=rindy-float(indy)
347
               frac1i=1.-frac0i
348
               frac1j=1.-frac0j
349
               xpos(i) = xcor(indx ,indy) * frac1i * frac1j
350
     &                 + xcor(indx ,indu) * frac1i * frac0j
351
     &                 + xcor(indr ,indy) * frac0i * frac1j
352
     &                 + xcor(indr ,indu) * frac0i * frac0j
353
               ypos(i) = ycor(indx ,indy) * frac1i * frac1j
354
     &                 + ycor(indx ,indu) * frac1i * frac0j
355
     &                 + ycor(indr ,indy) * frac0i * frac1j
356
     &                 + ycor(indr ,indu) * frac0i * frac0j
357
 
358
            else
359
 
360
               xpos(i)=mdv
361
               ypos(i)=mdv
362
 
363
            endif
364
 
365
         else
366
 
367
            rlon(i)=mdv
368
            rlat(i)=mdv
369
            xpos(i)=mdv
370
            ypos(i)=mdv
371
 
372
         endif
373
 
374
      enddo
375
 
376
      END    
377
 
378
 
379
c     --------------------------------------------------------------------------------
380
c     Transformation routine: LMSTOLM and PHSTOPH from library gm2em
381
c     --------------------------------------------------------------------------------
382
 
383
      REAL FUNCTION LMTOLMS (PHI, LAM, POLPHI, POLLAM)
384
C
385
C%Z% Modul %M%, V%I% vom %G%, extrahiert am %H%
386
C
387
C**** LMTOLMS  -   FC:UMRECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE LAM
388
C****                 AUF EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)
389
C****                 IM ROTIERTEN SYSTEM. DER NORDPOL DES SYSTEMS HAT
390
C****                 DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
391
C**   AUFRUF   :   LAM = LMTOLMS (PHI, LAM, POLPHI, POLLAM)
392
C**   ENTRIES  :   KEINE
393
C**   ZWECK    :   UMRECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE LAM AUF
394
C**                EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM
395
C**                ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
396
C**                DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
397
C**   VERSIONS-
398
C**   DATUM    :   03.05.90
399
C**
400
C**   EXTERNALS:   KEINE
401
C**   EINGABE-
402
C**   PARAMETER:   PHI    REAL BREITE DES PUNKTES IM GEOGR. SYSTEM
403
C**                LAM    REAL LAENGE DES PUNKTES IM GEOGR. SYSTEM
404
C**                POLPHI REAL GEOGR.BREITE DES N-POLS DES ROT. SYSTEMS
405
C**                POLLAM REAL GEOGR.LAENGE DES N-POLS DES ROT. SYSTEMS
406
C**   AUSGABE-
407
C**   PARAMETER:   WAHRE GEOGRAPHISCHE LAENGE ALS WERT DER FUNKTION
408
C**                ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)
409
C**
410
C**   COMMON-
411
C**   BLOECKE  :   KEINE
412
C**
413
C**   FEHLERBE-
414
C**   HANDLUNG :   KEINE
415
C**   VERFASSER:   G. DE MORSIER
416
 
417
      REAL        LAM,PHI,POLPHI,POLLAM
418
 
419
      DATA        ZRPI18 , ZPIR18  / 57.2957795 , 0.0174532925 /
420
 
421
      ZSINPOL = SIN(ZPIR18*POLPHI)
422
      ZCOSPOL = COS(ZPIR18*POLPHI)
423
      ZLAMPOL =     ZPIR18*POLLAM
424
      ZPHI    =     ZPIR18*PHI
425
      ZLAM    = LAM
426
      IF(ZLAM.GT.180.0) ZLAM = ZLAM - 360.0
427
      ZLAM    = ZPIR18*ZLAM
428
 
429
      ZARG1   = - SIN(ZLAM-ZLAMPOL)*COS(ZPHI)
430
      ZARG2   = - ZSINPOL*COS(ZPHI)*COS(ZLAM-ZLAMPOL)+ZCOSPOL*SIN(ZPHI)
431
      IF (ABS(ZARG2).LT.1.E-30) THEN
432
        IF (ABS(ZARG1).LT.1.E-30) THEN
433
          LMTOLMS =   0.0
434
        ELSEIF (ZARG1.GT.0.) THEN
435
              LMTOLMS =  90.0
436
            ELSE
437
              LMTOLMS = -90.0
438
            ENDIF
439
      ELSE
440
        LMTOLMS = ZRPI18*ATAN2(ZARG1,ZARG2)
441
      ENDIF
442
 
443
      RETURN
444
      END
445
 
446
 
447
      REAL FUNCTION PHTOPHS (PHI, LAM, POLPHI, POLLAM)
448
C
449
C%Z% Modul %M%, V%I% vom %G%, extrahiert am %H%
450
C
451
C**** PHTOPHS  -   FC:UMRECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE PHI
452
C****                 AUF EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)
453
C****                 IM ROTIERTEN SYSTEM. DER NORDPOL DES SYSTEMS HAT
454
C****                 DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
455
C**   AUFRUF   :   PHI = PHTOPHS (PHI, LAM, POLPHI, POLLAM)
456
C**   ENTRIES  :   KEINE
457
C**   ZWECK    :   UMRECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE PHI AUF
458
C**                EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM
459
C**                ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
460
C**                DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
461
C**   VERSIONS-
462
C**   DATUM    :   03.05.90
463
C**
464
C**   EXTERNALS:   KEINE
465
C**   EINGABE-
466
C**   PARAMETER:   PHI    REAL BREITE DES PUNKTES IM GEOGR. SYSTEM
467
C**                LAM    REAL LAENGE DES PUNKTES IM GEOGR. SYSTEM
468
C**                POLPHI REAL GEOGR.BREITE DES N-POLS DES ROT. SYSTEMS
469
C**                POLLAM REAL GEOGR.LAENGE DES N-POLS DES ROT. SYSTEMS
470
C**   AUSGABE-
471
C**   PARAMETER:   ROTIERTE BREITE PHIS ALS WERT DER FUNKTION
472
C**                ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)
473
C**
474
C**   COMMON-
475
C**   BLOECKE  :   KEINE
476
C**
477
C**   FEHLERBE-
478
C**   HANDLUNG :   KEINE
479
C**   VERFASSER:   G. DE MORSIER
480
 
481
      REAL        LAM,PHI,POLPHI,POLLAM
482
 
483
      DATA        ZRPI18 , ZPIR18  / 57.2957795 , 0.0174532925 /
484
 
485
      ZSINPOL = SIN(ZPIR18*POLPHI)
486
      ZCOSPOL = COS(ZPIR18*POLPHI)
487
      ZLAMPOL = ZPIR18*POLLAM
488
      ZPHI    = ZPIR18*PHI
489
      ZLAM    = LAM
490
      IF(ZLAM.GT.180.0) ZLAM = ZLAM - 360.0
491
      ZLAM    = ZPIR18*ZLAM
492
      ZARG    = ZCOSPOL*COS(ZPHI)*COS(ZLAM-ZLAMPOL) + ZSINPOL*SIN(ZPHI)
493
 
494
      PHTOPHS = ZRPI18*ASIN(ZARG)
495
 
496
      RETURN
497
      END
498
 
499
 
500
c     ----------------------------------------------------------------
501
c     Check whether variable is found on netcdf file
502
c     ----------------------------------------------------------------
503
 
504
      subroutine check_varok (isok,varname,varlist,nvars)
505
 
506
c     Check whether the variable <varname> is in the list <varlist(nvars)>.
507
c     If this is the case, <isok> is incremented by 1. Otherwise <isok>
508
c     keeps its value.
509
 
510
      implicit none
511
 
512
c     Declaraion of subroutine parameters
513
      integer      isok
514
      integer      nvars
515
      character*80 varname
516
      character*80 varlist(nvars)
517
 
518
c     Auxiliary variables
519
      integer      i
520
 
521
c     Main
522
      do i=1,nvars
523
         if (trim(varname).eq.trim(varlist(i))) isok=isok+1
524
      enddo
525
 
526
      end
527
 
528