Subversion Repositories lagranto.um

Rev

Rev 16 | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
16 michaesp 1
C     ********************************************************************
2
 
3
      program ptos
4
 
5
C     ********************************************************************
6
 
7
C     Purpose
8
C     -------
9
C
10
C     	Calculates secondary data files from primary data files
11
C     	(based upon IVE-routines).
12
C     	The program is invoked by the shell-script p2s.
13
C
14
C     Author
15
C     ------
16
C
17
C     	H. Wernli	April 96
18
C
19
C     Modifications
20
C     -------------
21
C
22
CSue Adpated to calculate f*geostrophic wind components and geostrophic 
23
C momentum components. 
24
C lat, lon, rlat, rlon modified to arrays so that can be precalculated
25
C lat in radians, rlat in degrees (similarly for lon)
26
C     ********************************************************************
27
 
28
 
29
      include "um_dims.inc"
30
 
31
      real	time(ntmax)
32
      real	sp(nxmax*nymax),cl(nxmax*nymax),f(nxmax*nymax)
33
      real      oro(nxmax*nymax),ps(nxmax*nymax)
34
      real	var(nxmax*nymax*nzmax),th(nxmax*nymax*nzmax),
35
     >		pv(nxmax*nymax*nzmax),the(nxmax*nymax*nzmax),
36
     >		thw(nxmax*nymax*nzmax),
37
     >		rh(nxmax*nymax*nzmax),dhr(nxmax*nymax*nzmax),
38
     >		tt(nxmax*nymax*nzmax),qq(nxmax*nymax*nzmax),
39
     >		uu(nxmax*nymax*nzmax),vv(nxmax*nymax*nzmax),
40
     >		ww(nxmax*nymax*nzmax),fug(nxmax*nymax*nzmax),
41
     >          fvg(nxmax*nymax*nzmax),dspdx(nxmax*nymax),
42
     >          dspdy(nxmax*nymax),dvardx(nxmax*nymax*nzmax),
43
     >          dvardy(nxmax*nymax*nzmax),Mg(nxmax*nymax*nzmax),
44
     >          Ng(nxmax*nymax*nzmax)
45
      character*80 cdfnam,cstnam,cdf_file
46
      integer	cdfid,cdfid1,cstid,ierr,ndim,vardim(4)
47
      real	dx,dy,mdv,varmin(4),varmax(4),stag(4)
48
      real      aklev(nzmax),bklev(nzmax),aklay(nzmax),bklay(nzmax),
49
     >		ak(nzmax),bk(nzmax)
50
      real lon_eq,lat_eq,u_ll,v_ll,Mg_ll,Ng_ll,lon1,lat1
51
      integer	nx,ny,nz,ntimes,i,j,k,n,kk,jj
52
      integer	stdate(5)
53
      integer	mode,qmode,zdef
54
      real	rlat(nxmax*nymax),rlon(nxmax*nymax),
55
     >          lat(nxmax*nymax),lon(nxmax*nymax)
56
      real	pollon,pollat,yphys,xphys
57
      real	phstoph,lmstolm
58
      logical	prelev
59
 
60
      real      pi
61
      data      pi /3.141592654/
62
 
63
      integer iw, jw, kkw, w
64
      parameter (w=7)
65
      real weight(w,w), sumweight, tmpu(nxmax*nymax*nzmax), 
66
     > tmpv(nxmax*nymax*nzmax), ug(nxmax*nymax*nzmax), 
67
     > vg(nxmax*nymax*nzmax)
68
      write(*,*)'*** start of program ptos ***'
69
 
70
C     Read filename
71
 
72
      read(9,10)cdfnam
73
      write(*,*) 'cdfnam is ',cdfnam
74
   10 format(a13)
75
 
76
C     Read mode and qmode
77
 
78
      read(9,*)mode
79
      read(9,*)qmode
80
      if (mode.eq.10) read(9,*)zdef
81
 
82
C     Open files and get infos about data domain
83
 
84
      if (mode.eq.10) then
85
        call cdfwopn('P'//trim(cdfnam),cdfid1,ierr)
86
      else
87
        call cdfopn('P'//trim(cdfnam),cdfid1,ierr)
88
      endif
89
      call getcfn(cdfid1,cstnam,ierr)
90
      call cdfopn(trim(cstnam),cstid,ierr)
91
 
92
      call getdef(cdfid1,'T',ndim,mdv,vardim,varmin,varmax,stag,ierr)
93
      if (ierr.ne.0) goto 920
94
      mdv=-999.98999
95
 
96
C     Get the levels, pole, etc.
97
 
98
      nx=vardim(1)
99
      ny=vardim(2)
100
      nz=vardim(3)
101
 
102
      call getgrid(cstid,dx,dy,ierr)
103
      call getlevs(cstid,nz,aklev,bklev,aklay,bklay,ierr)
104
      call getpole(cstid,pollon,pollat,ierr)
105
      call getstart(cstid,stdate,ierr)
106
C      write(*,*) 'dx ',dx,' dy ',dy
107
C      write(*,*) 'pollon ', pollon,' pollat ',pollat
108
 
109
C     Determine if data is on pressure or model levels
110
 
111
      prelev=.true.
112
      do k=1,nz
113
        if (bklev(k).ne.0.) prelev=.false.
114
      enddo
115
C      write(*,*) ' prelev ',prelev,'aklev ',aklev,' bklev ',bklev
116
 
117
CSue Calculate real lats and lons
118
      if ((pollon.ne.0.).or.(pollat.ne.90.)) then
119
        do j=1,ny
120
          do i=1,nx
121
            jj=i+(j-1)*nx
122
            rlat(jj)=varmin(2)+(j-1)*dy
123
            rlon(jj)=varmin(1)+(i-1)*dx      
124
            yphys=phstoph(rlat(jj),rlon(jj),pollat,pollon)
125
C           if I use sind(lat in deg): troubles at the N-pole
126
            lat(jj)=2.*pi*yphys/360.
127
            xphys=lmstolm(rlat(jj),rlon(jj),pollat,pollon)
128
            lon(jj)=2.*pi*xphys/360.
129
          enddo
130
        enddo
131
        else
132
        do j=1,ny
133
          do i=1,nx
134
            jj=i+(j-1)*nx
135
            lon(jj)=varmin(1)+(i-1)*dx 
136
            lat(jj)=varmin(2)+(j-1)*dy
137
          enddo
138
        enddo
139
      endif
140
 
141
C ---------------------------------
142
 
143
C     Calculate cos(latitude) array and the coriolis parameter
144
 
145
      if ((pollon.ne.0.).or.(pollat.ne.90.)) then
146
        do j=1,ny
147
          do i=1,nx
148
            jj=i+(j-1)*nx
149
            cl(i+(j-1)*nx)=cos(pi*rlat(jj)/180.)
150
            f(i+(j-1)*nx)=0.000145444*sin(lat(jj))
151
          enddo
152
        enddo
153
      else
154
        do j=1,ny
155
          do i=1,nx
156
            cl(i+(j-1)*nx)=cos(lat(jj))
157
            f(i+(j-1)*nx)=0.000145444*sin(lat(jj))
158
          enddo
159
        enddo
160
      endif
161
 
162
C     Determine if data is on levels or layers
163
 
164
      if (stag(3).eq.-0.5) then
165
        do k=1,nz
166
          ak(k)=aklay(k)
167
          bk(k)=bklay(k)
168
        enddo
169
      else
170
        do k=1,nz
171
          ak(k)=aklev(k)
172
          bk(k)=bklev(k)
173
        enddo
174
      endif
175
 
176
C     Get all the fields
177
 
178
      call gettimes(cdfid1,time,ntimes,ierr)
179
 
180
C     Loop over all times
181
      write(*,*) 'ntimes = ',ntimes
182
 
183
      do n=1,ntimes
184
 
185
        if (.not.prelev) then
186
          call getdat(cdfid1,'PS',time(n),0,sp,ierr)
187
          if (ierr.ne.0) goto 921
188
        else
189
          call getdat(cdfid1,'PS',time(n),0,ps,ierr)
190
          if (ierr.ne.0) goto 921
191
        endif
192
        call getdat(cdfid1,'T',time(n),0,tt,ierr)
193
        if (ierr.ne.0) goto 920
194
        if ((mode.eq.0).or.(mode.eq.2).or.(mode.eq.4).or.
195
     >      (mode.eq.5).or.(mode.eq.6).or.(mode.eq.1).or.
196
     >	    (mode.eq.9).or.(mode.eq.10)) then
197
          if (qmode.eq.1) then
198
            call getdat(cdfid1,'Q',time(n),0,qq,ierr)
199
          else
200
            call getdat(cdfid1,'QD',time(n),0,qq,ierr)
201
          endif
202
          if (ierr.ne.0) goto 922
203
        endif
204
        if (mode.lt.9 .or. mode.eq.10) then
205
          call getdat(cdfid1,'U',time(n),0,uu,ierr)
206
          if (ierr.ne.0) goto 923
207
          call getdat(cdfid1,'V',time(n),0,vv,ierr)
208
          if (ierr.ne.0) goto 924
209
          if ((mode.eq.0).or.(mode.eq.2).or.(mode.eq.8)) then
210
            call getdat(cdfid1,'OMEGA',time(n),0,ww,ierr)
211
            if (ierr.ne.0) goto 925
212
          endif
213
        endif
214
 
215
        if (mode.eq.10) then 
216
 
217
C         Calculation of the geopotential
218
 
219
C         first get the orography
220
          call getoro(oro,dx,dy,stdate(1),
221
     >                varmin(1),varmin(2),nx,ny,ierr)
222
          if (ierr.eq.2) then
223
            stop '*** error in subrountine getoro ***'
224
          endif
225
 
226
          call geopot(var,qq,tt,oro,sp,nx,ny,nz,ak,bk)
227
 
228
          if (zdef.eq.0) then
229
            if (n.eq.1) then
230
	      call putdef(cdfid1,'Z',4,mdv,
231
     >                    vardim,varmin,varmax,stag,ierr)
232
              write(*,*)
233
              write(*,*)'*** variable Z created on P-file'
234
            endif
235
          endif
236
 
237
          call putdat(cdfid1,'Z',time(n),0,var,ierr)
238
c          goto 900
239
        endif
240
 
241
C       Create the secondary data file
242
 
243
        if (n.eq.1) then
244
          cdf_file = 'S'//trim(cdfnam)
245
          call crecdf(cdf_file,cdfid,varmin,varmax,
246
     >                3,trim(cstnam),ierr)
247
          if (ierr.ne.0) goto 996
248
          write(*,*)
249
          write(*,*)'*** NetCDF file S',trim(cdfnam),
250
     >              ' created'
251
        endif
252
 
253
C       Put surface pressure on S-file. 
254
 
255
        if (.not.prelev) then
256
          vardim(3)=1
257
          if (n.eq.1) then
258
            call putdef(cdfid,'PS',4,mdv,vardim,varmin,varmax,stag,ierr)
259
            write(*,*)
260
            write(*,*)'*** variable PS created on S-file'
261
          endif
262
          call putdat(cdfid,'PS',time(n),0,sp,ierr)
263
          vardim(3)=nz
264
        else
265
          vardim(3)=1
266
          if (n.eq.1) then
267
            call putdef(cdfid,'PS',4,mdv,vardim,varmin,varmax,stag,ierr)
268
            write(*,*)
269
            write(*,*)'*** variable PS created on S-file'
270
          endif
271
          call putdat(cdfid,'PS',time(n),0,ps,ierr)
272
          vardim(3)=nz
273
        endif
274
 
275
C       Calculate the secondary data variables
276
 
277
C       Calculation of potential temperature
278
 
279
        if (mode.lt.9 .or. mode.eq.10) then
280
          call pottemp(th,tt,sp,nx,ny,nz,ak,bk)
281
 
282
          if (n.eq.1) then
283
	    call putdef(cdfid,'TH',4,mdv,
284
     >                  vardim,varmin,varmax,stag,ierr)
285
            write(*,*)
286
            write(*,*)'*** variable TH created on S-file'
287
          endif
288
 
289
          call putdat(cdfid,'TH',time(n),0,th,ierr)
290
        endif
291
 
292
C       Calculation of relative humidity
293
 
294
        if ((mode.eq.0).or.(mode.eq.2).or.(mode.eq.4).or.
295
     >      (mode.eq.5).or.(mode.eq.6).or.(mode.eq.1).or.
296
     >      (mode.eq.10)) then
297
 
298
          call relhum(rh,qq,tt,sp,nx,ny,nz,ak,bk)
299
 
300
          if (n.eq.1) then
301
	    call putdef(cdfid,'RH',4,mdv,
302
     >                  vardim,varmin,varmax,stag,ierr)
303
            write(*,*)
304
            write(*,*)'*** variable RH created on S-file'
305
          endif
306
 
307
          call putdat(cdfid,'RH',time(n),0,rh,ierr)
308
        endif
309
 
310
C       Calculation of relative vorticity
311
 
312
        if (mode.eq.999) then
313
          call relvort(var,uu,vv,sp,cl,f,nx,ny,nz,ak,bk,
314
     >                 varmin,varmax)
315
 
316
          if (n.eq.1) then
317
            call putdef(cdfid,'VO',4,mdv,
318
     >                  vardim,varmin,varmax,stag,ierr)
319
            write(*,*)
320
            write(*,*)'*** variable VO created on S-file'
321
          endif
322
 
323
          call putdat(cdfid,'VO',time(n),0,var,ierr)
324
        endif
325
 
326
C       Calculation of potential vorticity
327
 
328
        if (mode.eq.5) then
329
          call potvort(pv,uu,vv,th,sp,cl,f,nx,ny,nz,ak,bk,
330
     >    	     varmin,varmax)
331
 
332
          if (n.eq.1) then
333
	    call putdef(cdfid,'PV',4,mdv,
334
     >                  vardim,varmin,varmax,stag,ierr)
335
            write(*,*)
336
            write(*,*)'*** variable PV created on S-file'
337
          endif
338
 
339
          call putdat(cdfid,'PV',time(n),0,pv,ierr)
340
        endif
341
 
342
C       Calculation of equivalent potential temperature
343
 
344
        if ((mode.eq.0).or.(mode.eq.2).or.(mode.eq.1)
345
     >        .or.(mode.eq.10)) then
346
          call equpot(var,tt,qq,sp,nx,ny,nz,ak,bk)
347
 
348
          if (n.eq.1) then
349
	    call putdef(cdfid,'THE',4,mdv,
350
     >                  vardim,varmin,varmax,stag,ierr)
351
            write(*,*)
352
            write(*,*)'*** variable THE created on S-file'
353
          endif
354
           call putdat(cdfid,'THE',time(n),0,var,ierr)
355
        endif
356
 
357
C       Calculation of wet-bulb potential temperature
358
 
359
        if ((mode.eq.0).or.(mode.eq.2).or.(mode.eq.4)
360
     >       .or.(mode.eq.1).or.(mode.eq.5)) then
361
          call equpot(the,tt,qq,sp,nx,ny,nz,ak,bk)
362
          call wetbpt(var,the,sp,nx,ny,nz,ak,bk)
363
 
364
          if (n.eq.1) then
365
	    call putdef(cdfid,'THW',4,mdv,
366
     >                  vardim,varmin,varmax,stag,ierr)
367
            write(*,*)
368
            write(*,*)'*** variable THW created on S-file'
369
          endif
370
 
371
          call putdat(cdfid,'THW',time(n),0,var,ierr)
372
        endif
373
 
374
C       Calculation of the vertical gradient of the wet-bulb potential
375
C       temperature
376
 
377
        if ((mode.eq.4).or.(mode.eq.5)) then
378
c          if (mode.eq.5) then
379
c must calculate thw
380
c          call equpot(the,tt,qq,sp,nx,ny,nz,ak,bk)
381
c          call wetbpt(thw,the,sp,nx,ny,nz,ak,bk)
382
c          endif
383
 
384
          call dwetbptdp(var,thw,sp,nx,ny,nz,ak,bk)
385
 
386
          if (n.eq.1) then
387
            call putdef(cdfid,'DTHWDP',4,mdv,
388
     >                  vardim,varmin,varmax,stag,ierr)
389
            write(*,*)
390
            write(*,*)'*** variable DTHWDP created on S-file'
391
          endif
392
 
393
          call putdat(cdfid,'DTHWDP',time(n),0,var,ierr)
394
        endif
395
 
396
C       Calculation of the vertical gradient of potential temperature
397
 
398
        if ((mode.eq.999).or.(mode.eq.5)) then
399
          call dpottdp(var,th,sp,nx,ny,nz,ak,bk)
400
 
401
          if (n.eq.1) then
402
            call putdef(cdfid,'DTHDP',4,mdv,
403
     >                  vardim,varmin,varmax,stag,ierr)
404
            write(*,*)
405
            write(*,*)'*** variable DTHDP created on S-file'
406
          endif
407
 
408
          call putdat(cdfid,'DTHDP',time(n),0,var,ierr)
409
        endif
410
 
411
C       Calculation of the Richardson number
412
 
413
        if (mode.eq.4) then
414
          call richardson(var,uu,vv,th,sp,nx,ny,nz,ak,bk,mdv)
415
 
416
          if (n.eq.1) then
417
            call putdef(cdfid,'RI',4,mdv,
418
     >                  vardim,varmin,varmax,stag,ierr)
419
            write(*,*)
420
            write(*,*)'*** variable RI created on S-file'
421
          endif
422
 
423
          call putdat(cdfid,'RI',time(n),0,var,ierr)
424
        endif
425
 
426
C       Calculation of diabatic heating rate
427
 
428
        if ((mode.eq.0).or.(mode.eq.2).or.(mode.eq.8)) then
429
          call diabheat(dhr,th,ww,rh,sp,nx,ny,nz,ak,bk)
430
 
431
          if (n.eq.1) then
432
	    call putdef(cdfid,'CH',4,mdv,
433
     >                  vardim,varmin,varmax,stag,ierr)
434
            write(*,*)
435
            write(*,*)'*** variable CH created on S-file'
436
          endif
437
 
438
          call putdat(cdfid,'CH',time(n),0,dhr,ierr)
439
        endif
440
 
441
C       Calculation of diabatic PV rate
442
 
443
        if ((mode.eq.0).or.(mode.eq.2).or.(mode.eq.7)) then
444
          call diabpvr(var,uu,vv,dhr,sp,cl,f,nx,ny,nz,ak,bk,
445
     >                 varmin,varmax)
446
 
447
          if (n.eq.1) then
448
	    call putdef(cdfid,'PVR',4,mdv,
449
     >                  vardim,varmin,varmax,stag,ierr)
450
            write(*,*)
451
            write(*,*)'*** variable PVR created on S-file'
452
          endif
453
 
454
          call putdat(cdfid,'PVR',time(n),0,var,ierr)
455
        endif
456
 
457
C       Calculation of the geopotential
458
 
459
        if (mode.eq.9) then
460
 
461
C         first get the orography
462
          call getoro(oro,dx,dy,stdate(1),
463
     >                varmin(1),varmin(2),nx,ny,ierr)
464
          if (ierr.eq.2) then
465
            stop '*** error in subrountine getoro ***'
466
          endif
467
 
468
          call geopot(var,qq,tt,oro,sp,nx,ny,nz,ak,bk)
469
 
470
          if (n.eq.1) then
471
	    call putdef(cdfid,'Z',4,mdv,
472
     >                  vardim,varmin,varmax,stag,ierr)
473
            write(*,*)
474
            write(*,*)'*** variable Z created on S-file'
475
          endif
476
 
477
          call putdat(cdfid,'Z',time(n),0,var,ierr)
478
        endif
479
 
480
C       Calculation of the theta gradient
17 michaesp 481
 
16 michaesp 482
        if (mode.eq.11) then
483
 
484
          call pottemp(th,tt,sp,nx,ny,nz,ak,bk)
485
          call gradth(var,th,sp,prelev,cl,nx,ny,nz,ak,bk,varmin,varmax)
486
 
487
          if (n.eq.1) then
488
	    call putdef(cdfid,'GRADTH',4,mdv,
489
     >                  vardim,varmin,varmax,stag,ierr)
490
            write(*,*)
491
            write(*,*)'*** variable GRADTH created on S-file'
492
          endif
493
 
494
          call putdat(cdfid,'GRADTH',time(n),0,var,ierr)
495
        endif
496
 
497
C     Calculate Real W-E and N-S winds
498
 
499
      if (mode.eq.4) then
500
      if ((pollon.ne.0.).or.(pollat.ne.90.)) then
501
      do k=1,nz
502
        do j=1,ny
503
          lat_eq=varmin(2)+(j-1)*dy
504
          do i=1,nx
505
            lon_eq=varmin(1)+(i-1)*dx
506
            kk=i+(j-1)*nx+(k-1)*nx*ny
507
            u_ll=uvtougm(uu(kk),vv(kk),lon_eq,lat_eq,pollon,pollat)
508
            v_ll=uvtovgm(uu(kk),vv(kk),lon_eq,lat_eq,pollon,pollat)
509
            uu(kk)=u_ll
510
            vv(kk)=v_ll
511
          enddo
512
        enddo
513
      enddo
514
          if (n.eq.1) then
515
           call putdef(cdfid,'UREAL',4,mdv,
516
     >                  vardim,varmin,varmax,stag,ierr)
517
            write(*,*)
518
            write(*,*)'*** variable UREAL created on S-file'
519
          endif
520
 
521
          call putdat(cdfid,'UREAL',time(n),0,uu,ierr)
522
 
523
          if (n.eq.1) then
524
           call putdef(cdfid,'VREAL',4,mdv,
525
     >                  vardim,varmin,varmax,stag,ierr)
526
            write(*,*)
527
            write(*,*)'*** variable VREAL created on S-file'
528
          endif
529
 
530
          call putdat(cdfid,'VREAL',time(n),0,vv,ierr)
531
      endif
532
      endif
533
 
534
CSue  Calculate f*geostrophic winds
535
 
536
      if (mode.eq.999) then
537
 
538
C         first get the orography
539
          call getoro(oro,dx,dy,stdate(1),
540
     >                varmin(1),varmin(2),nx,ny,ierr)
541
          if (ierr.eq.2) then
542
            stop '*** error in subroutine getoro ***'
543
          endif
544
 
545
          call geopot(var,qq,tt,oro,sp,nx,ny,nz,ak,bk)
546
 
547
 
548
C         Then calculate the horizontal derivatives 
549
 
550
      if (prelev) then
551
        call ddh2(var,dvardx,cl,'X',nx,ny,1,varmin,varmax)
552
        call ddh2(var,dvardy,cl,'Y',nx,ny,1,varmin,varmax)
553
      else
554
        call ddh2(sp,dspdx,cl,'X',nx,ny,1,varmin,varmax)
555
        call ddh2(sp,dspdy,cl,'Y',nx,ny,1,varmin,varmax)
556
        call ddh3(var,dvardx,sp,dspdx,cl,'X',nx,ny,nz,
557
     >       varmin,varmax,ak,bk)
558
        call ddh3(var,dvardy,sp,dspdy,cl,'Y',nx,ny,nz,
559
     >       varmin,varmax,ak,bk)
560
      endif
561
 
562
 
563
C          Finally calculate Real W-E and N-S geostrophic winds
564
 
565
      if ((pollon.ne.0.).or.(pollat.ne.90.)) then
566
      do k=1,nz
567
        do j=1,ny
568
          lat_eq=varmin(2)+(j-1)*dy
569
          do i=1,nx
570
            lon_eq=varmin(1)+(i-1)*dx
571
            kk=i+(j-1)*nx+(k-1)*nx*ny
572
            fug(kk) = -9.806665*dvardy(kk)
573
            fvg(kk) = 9.806665*dvardx(kk)
574
            u_ll=uvtougm(fug(kk),fvg(kk),
575
     >           lon_eq,lat_eq,pollon,pollat)
576
            v_ll=uvtovgm(fug(kk),fvg(kk),
577
     >           lon_eq,lat_eq,pollon,pollat)
578
            fug(kk)=u_ll
579
            fvg(kk)=v_ll
580
          enddo
581
        enddo
582
      enddo
583
      else
584
      do k=1,nz
585
        do j=1,ny
586
          do i=1,nx          
587
            kk=i+(j-1)*nx+(k-1)*nx*ny
588
            fug(kk) = -9.806665*dvardy(kk)
589
            fvg(kk) = 9.806665*dvardx(kk)
590
          enddo
591
        enddo
592
       enddo
593
       endif
594
c       write(21) fug
595
c       write(22) fvg
596
c       write(23) f
597
 
598
          if (n.eq.1) then
599
           call putdef(cdfid,'FUGREAL',4,mdv,
600
     >                  vardim,varmin,varmax,stag,ierr)
601
            write(*,*)
602
            write(*,*)'*** variable FUGREAL created on S-file'
603
          endif
604
 
605
          call putdat(cdfid,'FUGREAL',time(n),0,fug,ierr)
606
 
607
          if (n.eq.1) then
608
           call putdef(cdfid,'FVGREAL',4,mdv,
609
     >                  vardim,varmin,varmax,stag,ierr)
610
            write(*,*)
611
            write(*,*)'*** variable FVGREAL created on S-file'
612
          endif
613
 
614
          call putdat(cdfid,'FVGREAL',time(n),0,fvg,ierr)
615
 
616
      endif ! if mode
617
 
618
CSue  Calculate smoothed geostrophic winds
619
 
620
      if (mode.eq.4) then
621
 
622
C         first get the orography
623
          call getoro(oro,dx,dy,stdate(1),
624
     >                varmin(1),varmin(2),nx,ny,ierr)
625
          if (ierr.eq.2) then
626
            stop '*** error in subroutine getoro ***'
627
          endif
628
 
629
          call geopot(var,qq,tt,oro,sp,nx,ny,nz,ak,bk)
630
 
631
 
632
C         Then calculate the horizontal derivatives 
633
 
634
      if (prelev) then
635
        call ddh2(var,dvardx,cl,'X',nx,ny,1,varmin,varmax)
636
        call ddh2(var,dvardy,cl,'Y',nx,ny,1,varmin,varmax)
637
      else
638
        call ddh2(sp,dspdx,cl,'X',nx,ny,1,varmin,varmax)
639
        call ddh2(sp,dspdy,cl,'Y',nx,ny,1,varmin,varmax)
640
        call ddh3(var,dvardx,sp,dspdx,cl,'X',nx,ny,nz,
641
     >       varmin,varmax,ak,bk)
642
        call ddh3(var,dvardy,sp,dspdy,cl,'Y',nx,ny,nz,
643
     >       varmin,varmax,ak,bk)
644
      endif
645
 
646
 
647
c Calculated model grid orientated ug and vg
648
      do k=1,nz
649
        do j=1,ny
650
          do i=1,nx          
651
            kk=i+(j-1)*nx+(k-1)*nx*ny
652
            jj=i+(j-1)*nx
653
            ug(kk) = -9.806665*dvardy(kk)/f(jj)
654
            vg(kk) = 9.806665*dvardx(kk)/f(jj)
655
          enddo
656
        enddo
657
       enddo
658
 
659
 
660
c Smooth winds
661
c      data weight/1, 2, 1, 2, 3, 2, 1, 2, 1/ 
662
      data weight/49*1/
663
      do k=1,nz
664
        do j=1,ny
665
          do i=1,nx
666
            kk=i+(j-1)*nx+(k-1)*nx*ny             
667
            sumweight = 0
668
            do jw = 1, w
669
               do iw = 1, w
670
                  iiw = i+iw-(int(w/2.)+1)
671
                  jjw = j+jw-(int(w/2.)+1)
672
                  if (iiw.lt.1) iiw = 1
673
                  if (jjw.lt.1) jjw = 1 
674
 
675
c            if (iiw.gt.nx) then 
676
c               write(*,*) 'nx ',nx,' ny ',ny,' i ',i,' j ',j
677
c               write(*,*) 'jw ',jw,' iw ',iw
678
c               write(*,*) 'ug ', ug(iiw+(jjw-1)*nx+(k-1)*nx*ny),
679
c     >                    'vg ', vg(iiw+(jjw-1)*nx+(k-1)*nx*ny)
680
c            do tmpjw = 1, w
681
c               do tmpiw = 1, w
682
c                  tmpiiw = i+tmpiw-(int(w/2.)+1)
683
c                  tmpjjw = j+tmpjw-(int(w/2.)+1)
684
c                  if (tmpiiw.lt.1) tmpiiw = 1
685
c                  if (tmpjjw.lt.1) tmpjjw = 1                
686
c                  if (tmpiiw.gt.nx) tmpiiw = nx
687
c                  if (tmpjjw.gt.ny) tmpjjw = ny
688
c                  kkw = tmpiiw+(tmpjjw-1)*nx+(k-1)*nx*ny   
689
c                  write(*,*) 'tmpiiw ', tmpiiw,' tmpjjw ', tmpjjw,
690
c     >                       ' vg ',vg(kkw),' vg ',vg(kkw)
691
c                enddo
692
c            enddo
693
c            stop
694
c            endif
695
 
696
                  if (iiw.gt.nx) iiw = nx
697
                  if (jjw.gt.ny) jjw = ny                  
698
                  kkw = iiw+(jjw-1)*nx+(k-1)*nx*ny   
699
                  tmpu(kk) = tmpu(kk) + ug(kkw)*weight(iw,jw)
700
                  tmpv(kk) = tmpv(kk) + vg(kkw)*weight(iw,jw) 
701
                  sumweight = sumweight+weight(iw,jw)
702
                enddo
703
            enddo
704
            tmpu(kk) = tmpu(kk)/sumweight
705
            tmpv(kk) = tmpv(kk)/sumweight
706
         enddo
707
        enddo
708
      enddo
709
 
710
      do k=1,nz
711
        do j=1,ny
712
          do i=1,nx
713
            kk=i+(j-1)*nx+(k-1)*nx*ny
714
       	    ug(kk) = tmpu(kk)
715
       	    vg(kk) = tmpv(kk)
716
         enddo
717
        enddo
718
      enddo
719
 
720
C          Finally calculate Real W-E and N-S geostrophic winds
721
 
722
      if ((pollon.ne.0.).or.(pollat.ne.90.)) then
723
      do k=1,nz
724
        do j=1,ny
725
          lat_eq=varmin(2)+(j-1)*dy
726
          do i=1,nx
727
            lon_eq=varmin(1)+(i-1)*dx
728
            kk=i+(j-1)*nx+(k-1)*nx*ny
729
            jj=i+(j-1)*nx
730
            u_ll=uvtougm(ug(kk),vg(kk),
731
     >           lon_eq,lat_eq,pollon,pollat)
732
            v_ll=uvtovgm(ug(kk),vg(kk),
733
     >           lon_eq,lat_eq,pollon,pollat)
734
            ug(kk)=u_ll
735
            vg(kk)=v_ll
736
          enddo
737
        enddo
738
      enddo
739
      endif
740
c       write(21) fug
741
c       write(22) fvg
742
c       write(23) f
743
 
744
c      open(111,file='ug.dat',status = 'unknown',form = 'unformatted')
745
c      write(111) ug      
746
c      open(112,file='vg.dat',status = 'unknown',form = 'unformatted')
747
c      write(112) vg
748
c      open(113,file='geo.dat',status = 'unknown',form = 'unformatted')
749
c      write(113) var
750
c      stop
751
 
752
          if (n.eq.1) then
753
           call putdef(cdfid,'UGREAL',4,mdv,
754
     >                  vardim,varmin,varmax,stag,ierr)
755
            write(*,*)
756
            write(*,*)'*** variable UGREAL created on S-file'
757
          endif
758
 
759
          call putdat(cdfid,'UGREAL',time(n),0,ug,ierr)
760
 
761
          if (n.eq.1) then
762
           call putdef(cdfid,'VGREAL',4,mdv,
763
     >                  vardim,varmin,varmax,stag,ierr)
764
            write(*,*)
765
            write(*,*)'*** variable VGREAL created on S-file'
766
          endif
767
 
768
          call putdat(cdfid,'VGREAL',time(n),0,vg,ierr)
769
 
770
 
771
cSue Calculate geostropic momentum components
772
 
773
      do k=1,nz
774
        do j=1,ny
775
          do i=1,nx    
776
            jj=i+(j-1)*nx
777
            kk=i+(j-1)*nx+(k-1)*nx*ny
778
            Mg(kk) = 0.000145444*6.378e6*sin(lat(jj))*
779
     >           (lon(jj) - lon(1))*cos(lat(jj)) 
780
     >           + vg(kk)
781
            Ng(kk) = 0.000145444*6.378e6*
782
     >           (cos(lat(1)) - cos(lat(jj))) - 
783
     >           ug(kk)
784
          enddo
785
        enddo
786
      enddo
787
 
788
          if (n.eq.1) then
789
           call putdef(cdfid,'MGREAL',4,mdv,
790
     >                  vardim,varmin,varmax,stag,ierr)
791
            write(*,*)
792
            write(*,*)'*** variable MGREAL created on S-file'
793
          endif
794
 
795
          call putdat(cdfid,'MGREAL',time(n),0,Mg,ierr)
796
 
797
          if (n.eq.1) then
798
           call putdef(cdfid,'NGREAL',4,mdv,
799
     >                  vardim,varmin,varmax,stag,ierr)
800
            write(*,*)
801
            write(*,*)'*** variable NGREAL created on S-file'
802
          endif
803
 
804
          call putdat(cdfid,'NGREAL',time(n),0,Ng,ierr)
805
 
806
       endif
807
 
808
      enddo !loop over n
809
 
810
 
811
C     Close the NetCDF files
812
 
813
      call clscdf(cdfid,ierr)
814
 900  continue
815
      call clscdf(cdfid1,ierr)
816
      call clscdf(cstid,ierr)
817
 
818
      goto 999
819
 
820
  920 stop '*** error: variable T not found on P-file ***'
821
  921 stop '*** error: variable PS not found on P-file ***'
822
  922 stop '*** error: variable Q not found on P-file ***'
823
  923 stop '*** error: variable U not found on P-file ***'
824
  924 stop '*** error: variable V not found on P-file ***'
825
  925 stop '*** error: variable OMEGA not found on P-file ***'
826
 
827
  996 stop '*** error: could not create S-file ***'
828
  999 continue
829
      end
830
 
831
      subroutine pottemp(pt,t,sp,ie,je,ke,ak,bk)
832
c     ==========================================
833
 
834
c     argument declaration
835
      integer	ie,je,ke
836
      real	pt(ie,je,ke),t(ie,je,ke),sp(ie,je),
837
     >		ak(ke),bk(ke)
838
 
839
c     variable declaration
840
      integer	i,j,k
841
      real	rdcp,tzero
842
      data	rdcp,tzero /0.286,273.15/
843
 
844
c     statement-functions for the computation of pressure
845
      real	pp,psrf
846
      integer	is
847
      pp(is)=ak(is)+bk(is)*psrf
848
 
849
c     computation of potential temperature
850
      do i=1,ie
851
      do j=1,je
852
        psrf=sp(i,j)
853
        do k=1,ke
854
c         distinction of temperature in K and deg C
855
          if (t(i,j,k).lt.100.) then
856
            pt(i,j,k)=(t(i,j,k)+tzero)*( (1000./pp(k))**rdcp )
857
          else
858
            pt(i,j,k)=t(i,j,k)*( (1000./pp(k))**rdcp )
859
          endif
860
        enddo
861
      enddo
862
      enddo
863
      end
864
 
865
      subroutine gradth(gth,th,sp,prelev,cl,ie,je,ke,ak,bk,vmin,vmax)
866
C     ===============================================================
867
 
868
c     argument declaration
869
      integer   ie,je,ke
870
      real      gth(ie,je,ke),th(ie,je,ke),sp(ie,je),cl(ie,je)
871
      real      ak(ke),bk(ke),vmin(4),vmax(4)
872
      logical	prelev
873
 
874
c     variable declaration
875
      include "um_dims.inc"
876
      real      dthdx(nxmax*nymax*nzmax),dthdy(nxmax*nymax*nzmax)
877
      real      dspdx(nxmax*nymax),dspdy(nxmax*nymax)
878
      integer   i,j,k,ind,ind2
879
 
880
      if (prelev) then
881
        call ddh2(th,dthdx,cl,'X',ie,je,1,vmin,vmax)
882
        call ddh2(th,dthdy,cl,'Y',ie,je,1,vmin,vmax)
883
      else
884
        call ddh2(sp,dspdx,cl,'X',ie,je,1,vmin,vmax)
885
        call ddh2(sp,dspdy,cl,'Y',ie,je,1,vmin,vmax)
886
        call ddh3(th,dthdx,sp,dspdx,cl,'X',ie,je,ke,vmin,vmax,ak,bk)
887
        call ddh3(th,dthdy,sp,dspdy,cl,'Y',ie,je,ke,vmin,vmax,ak,bk)
888
      endif
889
 
890
      do j=1,je
891
      do i=1,ie
892
        ind2=i+(j-1)*ie
893
        do k=1,ke
894
          ind=ind2+(k-1)*ie*je
895
          gth(i,j,k)=sqrt(dthdx(ind)**2.+dthdy(ind)**2.)
896
        enddo
897
      enddo
898
      enddo
899
      end
900
 
901
      subroutine geopot(psi,q,t,oro,sp,ie,je,ke,ak,bk)
902
c     ================================================
903
 
904
c     argument declaration
905
      integer   ie,je,ke
906
      real      psi(ie,je,ke),t(ie,je,ke),q(ie,je,ke),oro(ie,je),
907
     >          sp(ie,je),ak(ke),bk(ke)
908
 
909
c     variable declaration
910
      integer   i,j,k
911
      real      r,c,g
912
      data	r,c,g /287.,0.608,9.8/
913
 
914
c     statement-functions for the computation of pressure
915
      real      pp,psrf
916
      integer   is
917
      pp(is)=ak(is)+bk(is)*psrf
918
 
919
c     integration of geopotential height(special for first layer)
920
      do i=1,ie
921
      do j=1,je
922
        psrf=sp(i,j)
923
        psi(i,j,1)=oro(i,j)+r/g*
924
     >                 ((t(i,j,1)+273.15)*(1.+c*q(i,j,1))*
925
     >                (psrf-pp(1))/(0.5*(psrf+pp(1))))
926
c        psi(i,j,1)=1./g*(oro(i,j)
927
c     >                 +r*(t(i,j,1)+273.15)*(1.+c*q(i,j,1))*
928
c     >                (psrf-pp(1))/(0.5*(psrf+pp(1))))
929
      enddo
930
      enddo
931
      do j=1,je
932
      do i=1,ie
933
        psrf=sp(i,j)
934
        do k=2,ke
935
          psi(i,j,k)=psi(i,j,k-1)+r/g*
936
     >                  ((t(i,j,k-1)+273.15)*(1.+c*q(i,j,k-1))+
937
     >                   (t(i,j,k)+273.15)*(1.+c*q(i,j,k)))*
938
     >                  (pp(k-1)-pp(k))/(pp(k-1)+pp(k))
939
        enddo
940
      enddo
941
      enddo
942
      end
943
 
944
      subroutine getoro(oro,dx,dy,starty,lonmin,latmin,ie,je,ierr)
945
c     ===========================================================
946
c     reads the orography for the actual data domain from file
947
 
948
      include 'netcdf.inc'
949
      include "um_dims.inc"
950
 
951
c     argument declaration
952
      integer  ie,je
953
      real     or(nxmax*nymax)
954
      real     oro(ie,je), fld(nxmax,nymax)
955
 
956
c     variable declaration
957
      integer look(45)
958
      real rook(19)
959
      character*52 name
960
      character*52 filename
961
      integer  starty
962
 
963
c     open file with orography values and get them
964
c      write(*,*) 'Inset name of orography file'
965
c      read(*,*) name
966
c      filename = '/export/wombat/wombat-01/sws98slg/'//name
967
      write(*,*) ' Using orography file MESO17_orog_high_res.pp ' 
968
      open(10,file = 
969
     >'/export/cloud/stingjet/rb904381/umdata/orog_NAE.pp32',
970
     >  form = 'unformatted', status = 'old')
971
cNH     >'/export/wombat/wombat-01/sws98slg/MESO17_orog_high_res.pp',
972
cNH     >  form = 'unformatted', status = 'old')
973
c      open(10,file = '/export/wombat/wombat-01/sws98slg/'//name,
974
c     >  form = 'unformatted', status = 'old')
975
      read(10) look, rook
976
      nx = look(19)
977
      ny = look(18)
978
      read(10) (or(i), i = 1, nx*ny)
979
      close(10)
980
 
981
c data written starting at N-W corner, following loops assign data from S-W 
982
c corner and write out.
983
      do j = ny-1,0,-1
984
         i = j*nx
985
         do n = 1, nx
986
            fld(n,ny-j) = or(n+i)
987
         enddo
988
      enddo
989
        do i=1,ie
990
          do j=1,je
991
            oro(i,j)=fld(i+1,j+1)
992
          enddo
993
        enddo
994
 
995
 
996
      end
997
 
998
      subroutine wetbpt(thw,the,sp,ie,je,ke,ak,bk)
999
c     ============================================
1000
 
1001
c     argument declaration
1002
      integer   ie,je,ke
1003
      real      thw(ie,je,ke),the(ie,je,ke),
1004
     >          sp(ie,je),ak(ke),bk(ke)
1005
 
1006
c     variable declaration
1007
      real	tsa
1008
      integer	i,j,k
1009
 
1010
      do k=1,ke
1011
      do j=1,je
1012
      do i=1,ie
1013
        thw(i,j,k)=tsa(the(i,j,k)-273.15,1000.)+273.15
1014
      enddo
1015
      enddo
1016
      enddo
1017
      end
1018
 
1019
      subroutine dwetbptdp(dthwdp,thw,sp,ie,je,ke,ak,bk)
1020
c     ==================================================
1021
 
1022
c     argument declaration
1023
      integer   ie,je,ke
1024
      real      dthwdp(ie,je,ke),thw(ie,je,ke),
1025
     >          sp(ie,je),ak(ke),bk(ke)
1026
 
1027
      call ddp(thw,dthwdp,sp,ie,je,ke,ak,bk)
1028
 
1029
      end
1030
 
1031
      subroutine dpottdp(dthdp,th,sp,ie,je,ke,ak,bk)
1032
c     ==============================================
1033
 
1034
c     argument declaration
1035
      integer   ie,je,ke
1036
      real      dthdp(ie,je,ke),th(ie,je,ke),
1037
     >          sp(ie,je),ak(ke),bk(ke)
1038
 
1039
      call ddp(th,dthdp,sp,ie,je,ke,ak,bk)
1040
 
1041
      end
1042
 
1043
      subroutine richardson(ri,uu,vv,th,sp,ie,je,ke,ak,bk,mdv)
1044
c     ========================================================
1045
 
1046
c     argument declaration
1047
      integer   ie,je,ke
1048
      real      ri(ie,je,ke),uu(ie,je,ke),vv(ie,je,ke),
1049
     >          th(ie,je,ke),sp(ie,je),ak(ke),bk(ke),mdv
1050
c     variable declaration
1051
      include "um_dims.inc"
1052
      real      vel(nxmax,nymax,nzmax),dthdp(nxmax,nymax,nzmax),
1053
     &          dveldp(nxmax,nymax,nzmax)
1054
 
1055
c     statement-functions for the computation of pressure
1056
      real      pp,psrf
1057
      integer   is
1058
      pp(is)=ak(is)+bk(is)*psrf
1059
 
1060
      do k=1,ke
1061
      do j=1,je
1062
      do i=1,ie
1063
        vel(i,j,k)=(uu(i,j,k)**2.+vv(i,j,k)**2.)**0.5
1064
      enddo
1065
      enddo
1066
      enddo
1067
 
1068
      call ddp(th,dthdp,sp,ie,je,ke,ak,bk)
1069
      call ddp(vel,dveldp,sp,ie,je,ke,ak,bk)
1070
 
1071
c     computation of the richardson number
1072
c     use: kappa-1 = -0.714
1073
c          p0**(-kappa) = 0.0372
1074
      do i=1,ie
1075
      do j=1,je
1076
        psrf=sp(i,j)
1077
        do k=1,ke
1078
          if (abs(dveldp(i,j,k)).lt.0.0001) then
1079
            ri(i,j,k)=1000.
1080
          else
1081
            ri(i,j,k)=-287.*((100.*pp(k))**(-.714))*.0372*
1082
     >                dthdp(i,j,k)/(dveldp(i,j,k)**2.)
1083
          endif
1084
          if (ri(i,j,k).gt.1000.) ri(i,j,k)=1000.
1085
        enddo
1086
      enddo
1087
      enddo
1088
 
1089
      end
1090
 
1091
      real function tsa(os,p)
1092
c     =======================
1093
 
1094
C     This function returns the temperature tsa (celsius) on a saturation
1095
C     adiabat at pressure p (millibars). os is the equivalent potential
1096
C     temperature of the parcel (celsius). sign(a,b) replaces the
1097
C     algebraic sign of a with that of b.
1098
C     b is an empirical constant approximately equal to 0.001 of the latent
1099
C     heat of vaporization for water divided by the specific heat at constant
1100
C     pressure for dry air.
1101
 
1102
      real      a,b,os,p,tq,d,tqk,x,w
1103
      integer   i
1104
 
1105
      data b/2.6518986/
1106
      a=os+273.16
1107
 
1108
C     tq is the first guess for tsa
1109
 
1110
      tq=253.16
1111
 
1112
C     d is an initial value used in the iteration below
1113
 
1114
      d=120.
1115
 
1116
C     Iterate to obtain sufficient accuracy....see table 1, p.8
1117
C     of Stipanuk (1973) for equation used in iteration
1118
      do 1 i=1,12
1119
        tqk=tq-273.16
1120
        d=d/2.
1121
        x=a*exp(-b*w(tqk,p)/tq)-tq*((1000./p)**.286)
1122
        if (abs(x).lt.1e-7) go to 2
1123
        tq=tq+sign(d,x)
1124
 1    continue
1125
 2    tsa=tq-273.16
1126
      end
1127
 
1128
      real function w(t,p)
1129
c     ====================
1130
 
1131
C     This function returns the mixing ratio (grams of water vapor per
1132
C     kilogram of dry air) given the dew point (celsius) and pressure
1133
C     (millibars). If the temperature is input instead of the
1134
C     dew point, then saturation mixing ratio (same units) is returned.
1135
C     The formula is found in most meteorological texts.
1136
 
1137
      real      t,p,tkel,x,esat
1138
 
1139
      tkel=t+273.16
1140
      x=esat(tkel)      ! our function esat requires t in Kelvin
1141
      w=622.*x/(p-x)
1142
      end
1143
 
1144
      real function esat(t)
1145
c     =====================
1146
 
1147
C     This function returns the saturation vapor pressure over water (mb)
1148
C     given the temperature (Kelvin).
1149
C     The algorithm is due to Nordquist, W. S. ,1973: "Numerical
1150
C     Approximations of Selected Meteorological Parameters for Cloud
1151
C     Physics Problems" ECOM-5475, Atmospheric Sciences Laboratory, U. S.
1152
C     Army Electronics Command, White Sands Missile Range, New Mexico 88002.
1153
 
1154
      real p1,p2,c1,t
1155
 
1156
      p1=11.344-0.0303998*t
1157
      p2=3.49149-1302.8844/t
1158
      c1=23.832241-5.02808*log10(t)
1159
      esat=10.**(c1-1.3816e-7*10.**p1+8.1328e-3*10.**p2-2949.076/t)
1160
      end
1161
 
1162
      subroutine equpot(ap,t,q,sp,ie,je,ke,ak,bk)
1163
c     ===========================================
1164
 
1165
c     argument declaration
1166
      integer  ie,je,ke
1167
      real     ap(ie,je,ke),t(ie,je,ke),sp(ie,je)
1168
      real     q(ie,je,ke),ak(ke),bk(ke)
1169
 
1170
c     variable declaration
1171
      integer  i,j,k
1172
      real     rdcp,tzero
1173
      data     rdcp,tzero /0.286,273.15/
1174
 
1175
c     statement-functions for the computation of pressure
1176
      real      pp,psrf
1177
      integer   is
1178
      pp(is)=ak(is)+bk(is)*psrf
1179
 
1180
c     computation of potential temperature
1181
      do i=1,ie
1182
      do j=1,je
1183
        psrf=sp(i,j)
1184
        do k=1,ke
1185
          ap(i,j,k) = (t(i,j,k)+tzero)*(1000./pp(k))
1186
     +       **(0.2854*(1.0-0.28*q(i,j,k)))*exp(
1187
     +       (3.376/(2840.0/(3.5*alog(t(i,j,k)+tzero)-alog(
1188
     +       100.*pp(k)*max(1.0E-10,q(i,j,k))/(0.622+0.378*
1189
     +       q(i,j,k)))-0.1998)+55.0)-0.00254)*1.0E3*
1190
     +       max(1.0E-10,q(i,j,k))*(1.0+0.81*q(i,j,k)))
1191
        enddo
1192
      enddo
1193
      enddo
1194
      end
1195
 
1196
      subroutine relhum(rh,q,t,sp,ie,je,ke,ak,bk)
1197
c     ===========================================
1198
 
1199
c     argument declaration
1200
      integer   ie,je,ke
1201
      real      rh(ie,je,ke),t(ie,je,ke),q(ie,je,ke),
1202
     >		sp(ie,je),ak(ke),bk(ke)
1203
 
1204
c     variable declaration
1205
      integer   i,j,k
1206
      real	rdcp,tzero
1207
      real	b1,b2w,b3,b4w,r,rd,gqd,ge
1208
      data	rdcp,tzero /0.286,273.15/
1209
      data	b1,b2w,b3,b4w,r,rd /6.1078, 17.2693882, 273.16, 35.86,
1210
     &		287.05, 461.51/
1211
 
1212
c     statement-functions for the computation of pressure
1213
      real      pp,psrf
1214
      integer   is
1215
      pp(is)=ak(is)+bk(is)*psrf
1216
 
1217
      do i=1,ie
1218
      do j=1,je
1219
        psrf=sp(i,j)
1220
        do k=1,ke
1221
          ge = b1*exp(b2w*(t(i,j,k))/(t(i,j,k)+b3-b4w))
1222
          gqd= r/rd*ge/(pp(k)-(1.-r/rd)*ge)
1223
          rh(i,j,k)=100.*q(i,j,k)/gqd
1224
        enddo
1225
      enddo
1226
      enddo
1227
      end
1228
 
1229
      subroutine diabheat(dhr,t,w,rh,sp,ie,je,ke,ak,bk)
1230
c     =================================================
1231
 
1232
c     argument declaration
1233
      integer   ie,je,ke
1234
      real      dhr(ie,je,ke),t(ie,je,ke),w(ie,je,ke),
1235
     &          rh(ie,je,ke),sp(ie,je),ak(ke),bk(ke)
1236
 
1237
c     variable declaration
1238
      integer   i,j,k
1239
      real      p0,kappa,tzero
1240
      data      p0,kappa,tzero /1000.,0.286,273.15/
1241
      real      blog10,cp,r,lw,eps
1242
      data      blog10,cp,r,lw,eps /.08006,1004.,287.,2.5e+6,0.622/
1243
      real      esat,c,tt
1244
 
1245
c     statement-functions for the computation of pressure
1246
      real      pp,psrf
1247
      integer   is
1248
      pp(is)=ak(is)+bk(is)*psrf
1249
 
1250
c     computation of diabatic heating rate
1251
      do i=1,ie
1252
      do j=1,je
1253
        psrf=sp(i,j)
1254
        do k=1,ke
1255
          if (rh(i,j,k).lt.80.) then    ! only moist air of interest
1256
            dhr(i,j,k)=0.               ! cond. heating rate set to zero
1257
          else if (w(i,j,k).gt.0.) then ! cond. heating only for ascent
1258
            dhr(i,j,k)=0.
1259
          else
1260
            tt=t(i,j,k)*((pp(k)/p0)**kappa)   ! temp. from pot.temp.
1261
            c=lw/cp*eps*blog10*esat(tt)/pp(k)
1262
            dhr(i,j,k)=21600.*		      ! in units K per 6 hours
1263
     >        (1.-exp(.2*(80.-rh(i,j,k))))    ! weighting fun. for 80<RH<100
1264
     >        *(-c*kappa*t(i,j,k)*w(i,j,k)/(100.*pp(k)))/(1.+c)
1265
          endif
1266
        enddo
1267
      enddo
1268
      enddo
1269
      end
1270
 
1271
      subroutine diabpvr(dpvr,uu,vv,dhr,sp,cl,f,ie,je,ke,ak,bk,
1272
     >           vmin,vmax)
1273
C     =========================================================
1274
 
1275
c     argument declaration
1276
      integer   ie,je,ke
1277
      real      dpvr(ie,je,ke),uu(ie,je,ke),vv(ie,je,ke),
1278
     >          dhr(ie,je,ke),sp(ie,je),cl(ie,je),f(ie,je)
1279
      real      ak(ke),bk(ke),vmin(4),vmax(4)
1280
 
1281
c     variable declaration
1282
      include "um_dims.inc"
1283
      real      dhrdp(nxmax*nymax*nzmax),dudp(nxmax*nymax*nzmax),
1284
     >          dvdp(nxmax*nymax*nzmax),dvdx(nxmax*nymax*nzmax),
1285
     >          dhrdx(nxmax*nymax*nzmax),dudy(nxmax*nymax*nzmax),
1286
     >          dhrdy(nxmax*nymax*nzmax)
1287
      real      dspdx(nxmax*nymax),dspdy(nxmax*nymax)
1288
      integer   i,j,k,ind,ind2
1289
 
1290
      call ddp(dhr,dhrdp,sp,ie,je,ke,ak,bk)
1291
      call ddp(uu,dudp,sp,ie,je,ke,ak,bk)
1292
      call ddp(vv,dvdp,sp,ie,je,ke,ak,bk)
1293
      call ddh2(sp,dspdx,cl,'X',ie,je,1,vmin,vmax)
1294
      call ddh2(sp,dspdy,cl,'Y',ie,je,1,vmin,vmax)
1295
      call ddh3(dhr,dhrdx,sp,dspdx,cl,'X',ie,je,ke,vmin,vmax,ak,bk)
1296
      call ddh3(vv,dvdx,sp,dspdx,cl,'X',ie,je,ke,vmin,vmax,ak,bk)
1297
      call ddh3(dhr,dhrdy,sp,dspdy,cl,'Y',ie,je,ke,vmin,vmax,ak,bk)
1298
      call ddh3(uu,dudy,sp,dspdy,cl,'Y',ie,je,ke,vmin,vmax,ak,bk)
1299
 
1300
      do j=1,je
1301
      do i=1,ie
1302
        ind2=i+(j-1)*ie
1303
        do k=1,ke
1304
          ind=ind2+(k-1)*ie*je
1305
          dpvr(i,j,k)=-1.E6*9.80665*(
1306
     >		    -dvdp(ind)*dhrdx(ind)+dudp(ind)*dhrdy(ind)
1307
     >		    +(-dudy(ind)+dvdx(ind)+f(i,j))*dhrdp(ind))
1308
        enddo
1309
      enddo
1310
      enddo
1311
      end
1312
 
1313
      subroutine potvort(pv,uu,vv,th,sp,cl,f,ie,je,ke,ak,bk,
1314
     >		 vmin,vmax)
1315
C     ======================================================
1316
 
1317
c     argument declaration
1318
      integer   ie,je,ke
1319
      real      pv(ie,je,ke),uu(ie,je,ke),vv(ie,je,ke),
1320
     >          th(ie,je,ke),sp(ie,je),
1321
     >		cl(ie,je),f(ie,je)
1322
      real	ak(ke),bk(ke),vmin(4),vmax(4)
1323
 
1324
c     variable declaration
1325
      include "um_dims.inc"
1326
      real      dthdp(nxmax*nymax*nzmax),dudp(nxmax*nymax*nzmax),
1327
     >		dvdp(nxmax*nymax*nzmax),dvdx(nxmax*nymax*nzmax),
1328
     >          dthdx(nxmax*nymax*nzmax),dudy(nxmax*nymax*nzmax),
1329
     >		dthdy(nxmax*nymax*nzmax)
1330
      real      dspdx(nxmax*nymax),dspdy(nxmax*nymax)
1331
      integer	i,j,k,ind,ind2
1332
 
1333
      call ddp(th,dthdp,sp,ie,je,ke,ak,bk)
1334
      call ddp(uu,dudp,sp,ie,je,ke,ak,bk)
1335
      call ddp(vv,dvdp,sp,ie,je,ke,ak,bk)
1336
      call ddh2(sp,dspdx,cl,'X',ie,je,1,vmin,vmax)
1337
      call ddh2(sp,dspdy,cl,'Y',ie,je,1,vmin,vmax)
1338
      call ddh3(th,dthdx,sp,dspdx,cl,'X',ie,je,ke,vmin,vmax,ak,bk)
1339
      call ddh3(vv,dvdx,sp,dspdx,cl,'X',ie,je,ke,vmin,vmax,ak,bk)
1340
      call ddh3(th,dthdy,sp,dspdy,cl,'Y',ie,je,ke,vmin,vmax,ak,bk)
1341
      call ddh3(uu,dudy,sp,dspdy,cl,'Y',ie,je,ke,vmin,vmax,ak,bk)
1342
 
1343
      do j=1,je
1344
      do i=1,ie
1345
        ind2=i+(j-1)*ie
1346
        do k=1,ke
1347
          ind=ind2+(k-1)*ie*je
1348
          pv(i,j,k)=1.E6*9.80665*(
1349
     >		    -(-dudy(ind)+dvdx(ind)+f(i,j))*dthdp(ind)
1350
     >		    -(dudp(ind)*dthdy(ind)-dvdp(ind)*dthdx(ind)))
1351
        enddo
1352
      enddo
1353
      enddo
1354
      end
1355
 
1356
      subroutine relvort(vo,uu,vv,sp,cl,f,ie,je,ke,ak,bk,
1357
     >           vmin,vmax)
1358
C     ===================================================
1359
 
1360
c     argument declaration
1361
      integer   ie,je,ke
1362
      real      vo(ie,je,ke),uu(ie,je,ke),vv(ie,je,ke),
1363
     >          sp(ie,je),cl(ie,je),f(ie,je)
1364
      real      ak(ke),bk(ke),vmin(4),vmax(4)
1365
 
1366
c     variable declaration
1367
      include "um_dims.inc"
1368
      real      dvdx(nxmax*nymax*nzmax),dudy(nxmax*nymax*nzmax)
1369
      real      dspdx(nxmax*nymax),dspdy(nxmax*nymax)
1370
      integer   i,j,k,ind,ind2
1371
 
1372
      call ddh2(sp,dspdx,cl,'X',ie,je,1,vmin,vmax)
1373
      call ddh2(sp,dspdy,cl,'Y',ie,je,1,vmin,vmax)
1374
      call ddh3(vv,dvdx,sp,dspdx,cl,'X',ie,je,ke,vmin,vmax,ak,bk)
1375
      call ddh3(uu,dudy,sp,dspdy,cl,'Y',ie,je,ke,vmin,vmax,ak,bk)
1376
 
1377
      do j=1,je
1378
      do i=1,ie
1379
        ind2=i+(j-1)*ie
1380
        do k=1,ke
1381
          ind=ind2+(k-1)*ie*je
1382
          vo(i,j,k)=1.E4*(-dudy(ind)+dvdx(ind))
1383
        enddo
1384
      enddo
1385
      enddo
1386
      end
1387
 
1388
 
1389
      subroutine ddp(a,d,sp,ie,je,ke,ak,bk)
1390
c--------------------------------------------------------------------------
1391
c     Purpose: VERTICAL DERIVATIVE
1392
c       Compute the vertical derivative without missing data checking.
1393
c       The derivative is taken from array 'a' in the direction of 'P'.
1394
c       The result is stored in array 'd'.
1395
c       3 point weighted centered differencing is used.
1396
c       The vertical level-structure of the data is of the form
1397
c       p=ak(k)+bk(k)*ps.
1398
c--------------------------------------------------------------------------
1399
 
1400
c     declaration of arguments
1401
      integer       ie,je,ke
1402
      real          a(ie,je,ke),d(ie,je,ke),sp(ie,je)
1403
 
1404
c     variable declaration
1405
      integer       i,j,k
1406
      real          dpu,dpl,quot,fac,psrf
1407
      real          ak(ke),bk(ke)
1408
 
1409
c     specify scaling factor associated with derivation with respect
1410
c     to pressure
1411
      fac=0.01
1412
 
1413
c     compute vertical 3 point derivative
1414
c     ---------------------------
1415
c     3-point vertical derivative
1416
c     ---------------------------
1417
      do j=1,je
1418
      do i=1,ie
1419
c       get surface pressure at current grid-point
1420
        psrf=sp(i,j)
1421
c       points at k=1
1422
        dpu=(ak(1)+bk(1)*psrf)-(ak(2)+bk(2)*psrf)
1423
        d(i,j,1)=(a(i,j,1)-a(i,j,2))*fac/dpu
1424
c       points at 1<k<ke
1425
        do k=2,ke-1
1426
          dpu=(ak(k)+bk(k)*psrf)-(ak(k+1)+bk(k+1)*psrf)
1427
          dpl=(ak(k-1)+bk(k-1)*psrf)-(ak(k)+bk(k)*psrf)
1428
          quot=dpu/dpl
1429
          d(i,j,k)=(quot*(a(i,j,k-1)-a(i,j,k))
1430
     &              +1./quot*(a(i,j,k)-a(i,j,k+1)))*fac/(dpu+dpl)
1431
        enddo
1432
c       points at k=ke
1433
        dpl=(ak(ke-1)+bk(ke-1)*psrf)-(ak(ke)+bk(ke)*psrf)
1434
        d(i,j,ke)=(a(i,j,ke-1)-a(i,j,ke))*fac/dpl
1435
      enddo
1436
      enddo
1437
 
1438
      end
1439
 
1440
      subroutine ddh3(a,d,ps,dps,cl,dir,ie,je,ke,datmin,datmax,ak,bk)
1441
c-----------------------------------------------------------------------
1442
c     Purpose: HORIZONTAL DERIVATIVE ON PRESSURE-SURFACES WITHOUT MISSING DATA
1443
c       The derivative is taken from array 'a' in the direction of 'dir',
1444
c       where 'dir' is either 'X','Y'. The result is stored in array 'd'.
1445
c       The routine accounts for derivatives at the pole and periodic
1446
c       boundaries in the longitudinal direction (depending on
1447
c       the value of datmin, datmax). If the data-set does not reach to
1448
c       the pole, a one-sided derivative is taken. Pole-treatment is only
1449
c       carried out if the data-set covers 360 deg in longitude, and it
1450
c       requires that ie=4*ii+1, where ii is an integer.
1451
c     History:
1452
c       Daniel Luethi
1453
c-----------------------------------------------------------------------
1454
 
1455
c     declaration of arguments
1456
      integer       ie,je,ke
1457
      real          a(ie,je,ke),d(ie,je,ke),cl(ie,je)
1458
      real          ps(ie,je),dps(ie,je)
1459
      real          datmin(4),datmax(4)
1460
      character*(*) dir
1461
 
1462
c     variable declaration
1463
      integer       i,j,k
1464
      real          ak(ke),bk(ke),as(ke),bs(ke)
1465
 
1466
c     compute vertical derivatives of ak's and bk's
1467
      do k=2,ke-1
1468
        as(k)=(ak(k-1)-ak(k+1))/2.
1469
        bs(k)=(bk(k-1)-bk(k+1))/2.
1470
      enddo
1471
      as(1 )=ak(1)-ak(2)
1472
      bs(1 )=bk(1)-bk(2)
1473
      as(ke)=ak(ke-1)-ak(ke)
1474
      bs(ke)=bk(ke-1)-bk(ke)
1475
 
1476
c     compute horizontal derivatives on sigma surfaces
1477
      call ddh2(a,d,cl,dir,ie,je,ke,datmin,datmax)
1478
 
1479
c     apply correction for horizontal derivative on p-surfaces
1480
      do j=1,je
1481
      do i=1,ie
1482
        do k=2,ke-1
1483
          d(i,j,k)=d(i,j,k)+bk(k)*dps(i,j)/2./(as(k)+
1484
     &               bs(k)*ps(i,j))*(a(i,j,k+1)-a(i,j,k-1))
1485
        enddo
1486
        k=1
1487
        d(i,j,k)=d(i,j,k)+bk(k)*dps(i,j)/(as(k)+
1488
     &             bs(k)*ps(i,j))*(a(i,j,k+1)-a(i,j,k))
1489
        k=ke
1490
        d(i,j,k)=d(i,j,k)+bk(k)*dps(i,j)/(as(k)+
1491
     &             bs(k)*ps(i,j))*(a(i,j,k)-a(i,j,k-1))
1492
      enddo
1493
      enddo
1494
      end
1495
 
1496
      subroutine ddh2(a,d,cl,dir,ie,je,ke,datmin,datmax)
1497
c-----------------------------------------------------------------------
1498
c     Purpose: HORIZONTAL DERIVATIVE ON DATA-SURFACES WITHOUT MISSING DATA
1499
c       Compute the horizontal derivative without missing data checking.
1500
c       The derivative is taken from array 'a' in the direction of 'dir',
1501
c       where 'dir' is either 'X','Y'. The result is stored in array 'd'.
1502
c          The routine accounts for derivatives at the pole and periodic
1503
c       boundaries in the longitudinal direction (depending on
1504
c       the value of datmin, datmax). If the data-set does not reach to
1505
c       the pole, a one-sided derivative is taken. Pole-treatment is only
1506
c       carried out if the data-set covers 360 deg in longitude, and it
1507
c       requires that ie=4*ii+1, where ii is an integer.
1508
c-----------------------------------------------------------------------
1509
 
1510
c     declaration of arguments
1511
      integer       ie,je,ke
1512
      real          a(ie,je,ke),d(ie,je,ke),cl(ie,je)
1513
      real          datmin(4),datmax(4)
1514
      character*(*) dir
1515
 
1516
c     local variable declaration
1517
      integer       i,j,k,ip1,im1,jp1,jm1,ip,im,j1,j2
1518
      real          dlat,dlon,coslat,dx,dy,dxr,dyr
1519
      integer       northpl,southpl,lonper
1520
 
1521
c     rerd and circ are the mean radius and diameter of the earth in meter
1522
      real          rerd,circ,pi
1523
      data          rerd,circ,pi /6.37e6,4.e7,3.141592654/
1524
 
1525
c     compute flags for pole and periodic treatment
1526
      southpl=0
1527
      northpl=0
1528
      lonper =0
1529
      j1=1
1530
      j2=je
1531
      if (abs(datmax(1)-datmin(1)-360.).lt.1.e-3) then
1532
        lonper=1
1533
        if (abs(datmin(2)+90.).lt.1.e-3) then
1534
          southpl=1
1535
          j1=2
1536
        endif
1537
        if (abs(datmax(2)-90.).lt.1.e-3) then
1538
          northpl=1
1539
          j2=je-1
1540
        endif
1541
      endif
1542
 
1543
      dlon=((datmax(1)-datmin(1))/float(ie-1)) *pi/180.
1544
      dlat=((datmax(2)-datmin(2))/float(je-1)) *pi/180.
1545
 
1546
c      print *,'Computing derivative ',dir(1:1),
1547
c     &        ' of an array dimensioned ',ie,je,ke
1548
 
1549
      if (dir(1:1).eq.'X') then
1550
c       -----------------------------
1551
c       derivation in the x-direction
1552
c       -----------------------------
1553
        do k=1,ke
1554
 
1555
c         do gridpoints at j1<=j<=j2
1556
          do j=j1,j2
1557
           coslat=cl(1,j)
1558
 
1559
c          do regular gridpoints at 1<i<ie, 1<j<je
1560
           dx =rerd*coslat*dlon
1561
           dxr=1./(2.*dx)
1562
           do i=2,ie-1
1563
            ip1=i+1
1564
            im1=i-1
1565
            d(i,j,k)=dxr*(a(ip1,j,k)-a(im1,j,k))
1566
           enddo ! i-loop
1567
c          completed regular gridpoints at 1<i<ie, 1<j<je
1568
 
1569
c          do gridpoints at i=1, i=ie, 1<j<je
1570
           if (lonper.eq.1) then
1571
c           use periodic boundaries
1572
            i=1
1573
            ip1=2
1574
            im1=ie-1
1575
            d(i,j,k)=dxr*(a(ip1,j,k)-a(im1,j,k))
1576
            d(ie,j,k)=d(1,j,k)
1577
           else
1578
c           use one-sided derivatives
1579
            dxr=1./dx
1580
            i=1
1581
            ip1=2
1582
            im1=1
1583
            d(i,j,k)=dxr*(a(ip1,j,k)-a(im1,j,k))
1584
            i=ie
1585
            ip1=ie
1586
            im1=ie-1
1587
            d(i,j,k)=dxr*(a(ip1,j,k)-a(im1,j,k))
1588
           endif
1589
c          completed gridpoints at i=1, i=ie, j1<=j<=j2
1590
 
1591
          enddo ! j-loop
1592
c         completed gridpoints at 1<j<je
1593
 
1594
c         do gridpoints at j=je
1595
          if (northpl.eq.1) then
1596
c          for these gridpoints, the derivative in the x-direction is a
1597
c          derivative in the y-direction at another pole-gridpoint
1598
           dy =rerd*dlat
1599
           dyr=1./(2.*dy)
1600
           j=je
1601
           jp1=je-1
1602
           jm1=je-1
1603
           do i=1,ie
1604
            ip=mod(i-1+  (ie-1)/4,ie)+1
1605
            im=mod(i-1+3*(ie-1)/4,ie)+1
1606
            d(i,j,k)=dyr*(a(ip,jp1,k)-a(im,jm1,k))
1607
           enddo    ! i-loop
1608
c          completed gridpoints at j=je
1609
          endif
1610
c         do gridpoints at j=1
1611
          if (southpl.eq.1) then
1612
           dy =rerd*dlat
1613
           dyr=1./(2.*dy)
1614
           j=1
1615
           jp1=2
1616
           jm1=2
1617
           do i=1,ie
1618
            ip=mod(i-1+  (ie-1)/4,ie)+1
1619
            im=mod(i-1+3*(ie-1)/4,ie)+1
1620
            d(i,j,k)=dyr*(a(ip,jp1,k)-a(im,jm1,k))
1621
           enddo    ! i-loop
1622
          endif
1623
c         completed gridpoints at j=1
1624
 
1625
        enddo    ! k-loop
1626
 
1627
      else if (dir(1:1).eq.'Y') then
1628
c       -----------------------------
1629
c       derivation in the y-direction
1630
c       -----------------------------
1631
        dy =dlat*rerd
1632
        dyr=1./(2.*dy)
1633
        do k=1,ke
1634
         do i=1,ie
1635
 
1636
c          do regular gridpoints
1637
           do j=2,je-1
1638
            jp1=j+1
1639
            jm1=j-1
1640
            d(i,j,k)=dyr*(a(i,jp1,k)-a(i,jm1,k))
1641
           enddo
1642
 
1643
c          do gridpoints at j=je
1644
           if (northpl.eq.1) then
1645
c           pole-treatment
1646
            j=je
1647
            jm1=j-1
1648
            jp1=j-1
1649
            ip=mod(i-1+(ie-1)/2,ie)+1
1650
            im=i
1651
            d(i,j,k)=dyr*(a(ip,jp1,k)-a(im,jm1,k))
1652
           else
1653
c           one-sided derivative
1654
            j=je
1655
            jm1=j-1
1656
            jp1=j
1657
            d(i,j,k)=2.*dyr*(a(i,jp1,k)-a(i,jm1,k))
1658
           endif
1659
c          completed gridpoints at j=je
1660
 
1661
c          do gridpoints at j=1
1662
           if (southpl.eq.1) then
1663
c           pole-treatment
1664
            j=1
1665
            jm1=2
1666
            jp1=2
1667
            ip=i
1668
            im=mod(i-1+(ie-1)/2,ie)+1
1669
            d(i,j,k)=dyr*(a(ip,jp1,k)-a(im,jm1,k))
1670
           else
1671
c           one-sided derivative
1672
            j=1
1673
            jm1=1
1674
            jp1=2
1675
            d(i,j,k)=2.*dyr*(a(i,jp1,k)-a(i,jm1,k))
1676
           endif
1677
c          completed gridpoints at j=1
1678
 
1679
         enddo
1680
        enddo
1681
 
1682
 
1683
      endif
1684
      end
1685
c--------------------------------------------------------------------------
1686
c     transformation of coordinates with a rotated pole
1687
c--------------------------------------------------------------------------
1688
c     the following routines allow for the transformation of
1689
c     spherical coordinates to a spherical coordinate system with a rotated
1690
c     pole, and vice versa.
1691
c
1692
c     the arguments of the routines use the following conventions
1693
c        PHI,LAM:       latitude and longitude in regular spherical coordinates
1694
c        lat,lon:       dito
1695
c        PHIS,LAMS:     latitude and longitude in rotated spherical coordinates
1696
c        latp,lonp:     dito
1697
c        polphi,pollat: latitude and longitude of the rotated pole as
1698
c                       expressed in regular spherical coordinates
1699
c        latpol,lonpol: dito
1700
c        u,v:           windvector in regular spherical coordinates
1701
c        up,vp:         windvector in rotated spherical coordinates
1702
c     all angles are expressed in degrees.
1703
c
1704
c     For a given rotated pole location (pollat,polphi), the routines	
1705
c     provide the following transformations:
1706
c         REAL FUNCTION PHTOPHS: (lam,phi) -> phis
1707
c         REAL FUNCTION PHSTOPH: (lams,phis) -> phi
1708
c         REAL FUNCTION LMTOLMS: (lam,phi) -> lams
1709
c         REAL FUNCTION LMSTOLM: (lams,phis) -> lams
1710
c         real function uvtouem: (u,v,lon,lat) -> up
1711
c         real function uvtovem: (u,v,lon,lat) -> vp
1712
c         real function uvtougm: (up,vp,lonp,latp) -> u
1713
c         real function uvtovgm: (up,vp,lonp,latp) -> v
1714
c--------------------------------------------------------------------------
1715
 
1716
      REAL FUNCTION PHTOPHS (PHI, LAM, POLPHI, POLLAM)
1717
C
1718
C%Z% Modul %M%, V%I% vom %G%, extrahiert am %H%
1719
C
1720
C**** PHTOPHS  -   FC:UMRECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE PHI
1721
C****                 AUF EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)
1722
C****                 IM ROTIERTEN SYSTEM. DER NORDPOL DES SYSTEMS HAT
1723
C****                 DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
1724
C**   AUFRUF   :   PHI = PHTOPHS (PHI, LAM, POLPHI, POLLAM)
1725
C**   ENTRIES  :   KEINE
1726
C**   ZWECK    :   UMRECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE PHI AUF
1727
C**                EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM
1728
C**                ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
1729
C**                DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
1730
C**   VERSIONS-
1731
C**   DATUM    :   03.05.90
1732
C**
1733
C**   EXTERNALS:   KEINE
1734
C**   EINGABE-
1735
C**   PARAMETER:   PHI    REAL BREITE DES PUNKTES IM GEOGR. SYSTEM
1736
C**                LAM    REAL LAENGE DES PUNKTES IM GEOGR. SYSTEM
1737
C**                POLPHI REAL GEOGR.BREITE DES N-POLS DES ROT. SYSTEMS
1738
C**                POLLAM REAL GEOGR.LAENGE DES N-POLS DES ROT. SYSTEMS
1739
C**   AUSGABE-
1740
C**   PARAMETER:   ROTIERTE BREITE PHIS ALS WERT DER FUNKTION
1741
C**                ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)
1742
C**
1743
C**   COMMON-
1744
C**   BLOECKE  :   KEINE
1745
C**
1746
C**   FEHLERBE-
1747
C**   HANDLUNG :   KEINE
1748
C**   VERFASSER:   G. DE MORSIER
1749
 
1750
      REAL        LAM,PHI,POLPHI,POLLAM
1751
 
1752
      DATA        ZRPI18 , ZPIR18  / 57.2957795 , 0.0174532925 /
1753
 
1754
      ZSINPOL = SIN(ZPIR18*POLPHI)
1755
      ZCOSPOL = COS(ZPIR18*POLPHI)
1756
      ZLAMPOL = ZPIR18*POLLAM
1757
      ZPHI    = ZPIR18*PHI
1758
      ZLAM    = LAM
1759
      IF(ZLAM.GT.180.0) ZLAM = ZLAM - 360.0
1760
      ZLAM    = ZPIR18*ZLAM
1761
      ZARG    = ZCOSPOL*COS(ZPHI)*COS(ZLAM-ZLAMPOL) + ZSINPOL*SIN(ZPHI)
1762
 
1763
      PHTOPHS = ZRPI18*ASIN(ZARG)
1764
 
1765
      RETURN
1766
      END
1767
      REAL FUNCTION PHSTOPH (PHIS, LAMS, POLPHI, POLLAM)
1768
C
1769
C%Z% Modul %M%, V%I% vom %G%, extrahiert am %H%
1770
C
1771
C**** PHSTOPH  -   FC:BERECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE FUER
1772
C****                 EINEN PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM
1773
C****                 ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
1774
C****                 DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
1775
C**   AUFRUF   :   PHI = PHSTOPH (PHIS, LAMS, POLPHI, POLLAM)
1776
C**   ENTRIES  :   KEINE
1777
C**   ZWECK    :   BERECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE FUER
1778
C**                EINEN PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM
1779
C**                ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
1780
C**                DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
1781
C**   VERSIONS-
1782
C**   DATUM    :   03.05.90
1783
C**
1784
C**   EXTERNALS:   KEINE
1785
C**   EINGABE-
1786
C**   PARAMETER:   PHIS     REAL   GEOGR. BREITE DES PUNKTES IM ROT.SYS.
1787
C**                LAMS     REAL   GEOGR. LAENGE DES PUNKTES IM ROT.SYS.
1788
C**                POLPHI   REAL   WAHRE GEOGR. BREITE DES NORDPOLS
1789
C**                POLLAM   REAL   WAHRE GEOGR. LAENGE DES NORDPOLS
1790
C**   AUSGABE-
1791
C**   PARAMETER:   WAHRE GEOGRAPHISCHE BREITE ALS WERT DER FUNKTION
1792
C**                ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)
1793
C**
1794
C**   COMMON-
1795
C**   BLOECKE  :   KEINE
1796
C**
1797
C**   FEHLERBE-
1798
C**   HANDLUNG :   KEINE
1799
C**   VERFASSER:   D.MAJEWSKI
1800
 
1801
      REAL        LAMS,PHIS,POLPHI,POLLAM
1802
 
1803
      DATA        ZRPI18 , ZPIR18  / 57.2957795 , 0.0174532925 /
1804
 
1805
      SINPOL = SIN(ZPIR18*POLPHI)
1806
      COSPOL = COS(ZPIR18*POLPHI)
1807
      ZPHIS  = ZPIR18*PHIS
1808
      ZLAMS  = LAMS
1809
      IF(ZLAMS.GT.180.0) ZLAMS = ZLAMS - 360.0
1810
      ZLAMS  = ZPIR18*ZLAMS
1811
      ARG     = COSPOL*COS(ZPHIS)*COS(ZLAMS) + SINPOL*SIN(ZPHIS)
1812
 
1813
      PHSTOPH = ZRPI18*ASIN(ARG)
1814
 
1815
      RETURN
1816
      END
1817
      REAL FUNCTION LMTOLMS (PHI, LAM, POLPHI, POLLAM)
1818
C
1819
C%Z% Modul %M%, V%I% vom %G%, extrahiert am %H%
1820
C
1821
C**** LMTOLMS  -   FC:UMRECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE LAM
1822
C****                 AUF EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)
1823
C****                 IM ROTIERTEN SYSTEM. DER NORDPOL DES SYSTEMS HAT
1824
C****                 DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
1825
C**   AUFRUF   :   LAM = LMTOLMS (PHI, LAM, POLPHI, POLLAM)
1826
C**   ENTRIES  :   KEINE
1827
C**   ZWECK    :   UMRECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE LAM AUF
1828
C**                EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM
1829
C**                ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
1830
C**                DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
1831
C**   VERSIONS-
1832
C**   DATUM    :   03.05.90
1833
C**
1834
C**   EXTERNALS:   KEINE
1835
C**   EINGABE-
1836
C**   PARAMETER:   PHI    REAL BREITE DES PUNKTES IM GEOGR. SYSTEM
1837
C**                LAM    REAL LAENGE DES PUNKTES IM GEOGR. SYSTEM
1838
C**                POLPHI REAL GEOGR.BREITE DES N-POLS DES ROT. SYSTEMS
1839
C**                POLLAM REAL GEOGR.LAENGE DES N-POLS DES ROT. SYSTEMS
1840
C**   AUSGABE-
1841
C**   PARAMETER:   WAHRE GEOGRAPHISCHE LAENGE ALS WERT DER FUNKTION
1842
C**                ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)
1843
C**
1844
C**   COMMON-
1845
C**   BLOECKE  :   KEINE
1846
C**
1847
C**   FEHLERBE-
1848
C**   HANDLUNG :   KEINE
1849
C**   VERFASSER:   G. DE MORSIER
1850
 
1851
      REAL        LAM,PHI,POLPHI,POLLAM
1852
 
1853
      DATA        ZRPI18 , ZPIR18  / 57.2957795 , 0.0174532925 /
1854
 
1855
      ZSINPOL = SIN(ZPIR18*POLPHI)
1856
      ZCOSPOL = COS(ZPIR18*POLPHI)
1857
      ZLAMPOL =     ZPIR18*POLLAM
1858
      ZPHI    =     ZPIR18*PHI
1859
      ZLAM    = LAM
1860
      IF(ZLAM.GT.180.0) ZLAM = ZLAM - 360.0
1861
      ZLAM    = ZPIR18*ZLAM
1862
 
1863
      ZARG1   = - SIN(ZLAM-ZLAMPOL)*COS(ZPHI)
1864
      ZARG2   = - ZSINPOL*COS(ZPHI)*COS(ZLAM-ZLAMPOL)+ZCOSPOL*SIN(ZPHI)
1865
      IF (ABS(ZARG2).LT.1.E-30) THEN
1866
        IF (ABS(ZARG1).LT.1.E-30) THEN
1867
          LMTOLMS =   0.0
1868
        ELSEIF (ZARG1.GT.0.) THEN
1869
              LMTOLMS =  90.0
1870
            ELSE
1871
              LMTOLMS = -90.0
1872
            ENDIF
1873
      ELSE
1874
        LMTOLMS = ZRPI18*ATAN2(ZARG1,ZARG2)
1875
      ENDIF
1876
 
1877
      RETURN
1878
      END
1879
      REAL FUNCTION LMSTOLM (PHIS, LAMS, POLPHI, POLLAM)
1880
C
1881
C%Z% Modul %M%, V%I% vom %G%, extrahiert am %H%
1882
C
1883
C**** LMSTOLM  -   FC:BERECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE FUER
1884
C****                 EINEN PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)
1885
C****                 IM ROTIERTEN SYSTEM. DER NORDPOL DES SYSTEMS HAT
1886
C****                 DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
1887
C**   AUFRUF   :   LAM = LMSTOLM (PHIS, LAMS, POLPHI, POLLAM)
1888
C**   ENTRIES  :   KEINE
1889
C**   ZWECK    :   BERECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE FUER
1890
C**                EINEN PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)
1891
C**                IM ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
1892
C**                DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
1893
C**   VERSIONS-
1894
C**   DATUM    :   03.05.90
1895
C**
1896
C**   EXTERNALS:   KEINE
1897
C**   EINGABE-
1898
C**   PARAMETER:   PHIS     REAL   GEOGR. BREITE DES PUNKTES IM ROT.SYS.
1899
C**                LAMS     REAL   GEOGR. LAENGE DES PUNKTES IM ROT.SYS.
1900
C**                POLPHI   REAL   WAHRE GEOGR. BREITE DES NORDPOLS
1901
C**                POLLAM   REAL   WAHRE GEOGR. LAENGE DES NORDPOLS
1902
C**   AUSGABE-
1903
C**   PARAMETER:   WAHRE GEOGRAPHISCHE LAENGE ALS WERT DER FUNKTION
1904
C**                ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)
1905
C**
1906
C**   COMMON-
1907
C**   BLOECKE  :   KEINE
1908
C**
1909
C**   FEHLERBE-
1910
C**   HANDLUNG :   KEINE
1911
C**   VERFASSER:   D.MAJEWSKI
1912
 
1913
      REAL        LAMS,PHIS,POLPHI,POLLAM
1914
 
1915
      DATA        ZRPI18 , ZPIR18  / 57.2957795 , 0.0174532925 /
1916
 
1917
      ZSINPOL = SIN(ZPIR18*POLPHI)
1918
      ZCOSPOL = COS(ZPIR18*POLPHI)
1919
      ZLAMPOL = ZPIR18*POLLAM
1920
      ZPHIS   = ZPIR18*PHIS
1921
      ZLAMS   = LAMS
1922
      IF(ZLAMS.GT.180.0) ZLAMS = ZLAMS - 360.0
1923
      ZLAMS   = ZPIR18*ZLAMS
1924
 
1925
      ZARG1   = SIN(ZLAMPOL)*(- ZSINPOL*COS(ZLAMS)*COS(ZPHIS)  +
1926
     1                          ZCOSPOL*           SIN(ZPHIS)) -
1927
     2          COS(ZLAMPOL)*           SIN(ZLAMS)*COS(ZPHIS)
1928
      ZARG2   = COS(ZLAMPOL)*(- ZSINPOL*COS(ZLAMS)*COS(ZPHIS)  +
1929
     1                          ZCOSPOL*           SIN(ZPHIS)) +
1930
     2          SIN(ZLAMPOL)*           SIN(ZLAMS)*COS(ZPHIS)
1931
      IF (ABS(ZARG2).LT.1.E-30) THEN
1932
        IF (ABS(ZARG1).LT.1.E-30) THEN
1933
          LMSTOLM =   0.0
1934
        ELSEIF (ZARG1.GT.0.) THEN
1935
              LMSTOLAM =  90.0
1936
            ELSE
1937
              LMSTOLAM = -90.0
1938
            ENDIF
1939
      ELSE
1940
        LMSTOLM = ZRPI18*ATAN2(ZARG1,ZARG2)
1941
      ENDIF
1942
 
1943
      RETURN
1944
      END
1945
c-----------------------------------------------------------------
1946
      real function uvtouem(u,v,lon,lat,lonpol,latpol)
1947
c-----------------------------------------------------------------
1948
c     Rene Fehlmann, Wed Nov  1 13:31:28 MET 1995
1949
c     INPUT:
1950
c     - winds in geographical coordinated system
1951
c     - geographical  coordinates
1952
c     - coordinates of the rotated pole
1953
c     OUTPUT
1954
c     - u-component of the rotated wind
1955
c-----------------------------------------------------------------
1956
 
1957
      real               u,v,lon,lat,lonpol,latpol
1958
      real               lampol,phipol,lam,phi
1959
      real               cosdlam,cosphipol,cosphi
1960
      real               sindlam,sinphipol,sinphi
1961
 
1962
      data               pid180 / 0.0174532925199 /
1963
 
1964
c-----------------------------------------------------------------
1965
 
1966
      lampol=lonpol*pid180
1967
      phipol=latpol*pid180
1968
      lam=lon*pid180
1969
      phi=lat*pid180
1970
      cosdlam=cos(lampol-lam)
1971
      cosphipol=cos(phipol)
1972
      cosphi=cos(phi)
1973
      sinphipol=sin(phipol)
1974
      sinphi=sin(phi)
1975
      cosdlam=cosdlam
1976
      sindlam=sin(lampol-lam)
1977
 
1978
      uvtouem=-(sqrt(1.-(cosdlam*cosphipol*cosphi+ 
1979
     >     sinphipol*sinphi)**2.)*(-(u*cosdlam/
1980
     >     ((cosdlam*cosphi*sinphipol- 
1981
     >     cosphipol*sinphi))) + sindlam*
1982
     >     (-(v*sinphi/
1983
     >     ((cosdlam*cosphi*sinphipol - 
1984
     >     cosphipol*sinphi))) - 
1985
     >     cosphi*(-(v*cosphipol*cosphi) + 
1986
     >     sinphipol*(u*sindlam - v*cosdlam*sinphi))/
1987
     >     (cos(lampol - lam)*cosphi*sinphipol - 
1988
     >     cosphipol*sinphi)**2))/
1989
     >     (1. + cosphi**2. *sindlam**2./
1990
     >     (cosdlam*cosphi*sinphipol - 
1991
     >     cosphipol*sinphi)**2.))
1992
 
1993
      return
1994
      end
1995
c-----------------------------------------------------------------
1996
      real function uvtovem(u,v,lon,lat,lonpol,latpol)
1997
c-----------------------------------------------------------------
1998
c     Rene Fehlmann, Wed Nov  1 13:31:28 MET 1995
1999
c     INPUT:
2000
c     - winds in geographical coordinated system
2001
c     - geographical  coordinates
2002
c     - coordinates of the rotated pole
2003
c     OUTPUT
2004
c     - v-component of the rotated wind
2005
c-----------------------------------------------------------------
2006
 
2007
      real               u,v,lon,lat,lonpol,latpol
2008
      real               lampol,phipol,lam,phi
2009
 
2010
      real               cosdlam,cosphipol,cosphi
2011
      real               sindlam,sinphipol,sinphi
2012
 
2013
      data               pid180 / 0.0174532925199 /
2014
 
2015
c-----------------------------------------------------------------
2016
 
2017
      lampol=lonpol*pid180
2018
      phipol=latpol*pid180
2019
      lam=lon*pid180
2020
      phi=lat*pid180
2021
      cosdlam=cos(lampol-lam)
2022
      cosphipol=cos(phipol)
2023
      cosphi=cos(phi)
2024
      sinphipol=sin(phipol)
2025
      sinphi=sin(phi)
2026
      cosdlam=cosdlam
2027
      sindlam=sin(lampol-lam)
2028
 
2029
      uvtovem=(v*cosphi*sinphipol + cosphipol*
2030
     >     (u*sindlam-v*cosdlam*sinphi))/
2031
     >     sqrt(1.-(cosdlam*cosphipol*cosphi + 
2032
     >     sinphipol*sinphi)**2.)
2033
 
2034
      return
2035
      end
2036
c-----------------------------------------------------------------
2037
      real function uvtougm(up,vp,lonp,latp,lonpol,latpol)
2038
c-----------------------------------------------------------------
2039
c     Rene Fehlmann, Wed Nov  1 13:31:28 MET 1995
2040
c     INPUT:
2041
c     - rotated winds (i.e. wrt equatorial grid)
2042
c     - rotated coordinates (i.e. on equatorial grid)
2043
c     - coordinates of the rotated pole
2044
c     OUTPUT
2045
c     - u-component of the wind (i.e. real W-E wind)
2046
c-----------------------------------------------------------------
2047
 
2048
      real               up,vp,lonp,latp,lonpol,latpol
2049
      real               lampol,phipol,lamp,phip,pid180
2050
 
2051
      data               pid180 / 0.0174532925199 /
2052
 
2053
c-----------------------------------------------------------------
2054
 
2055
      phip=latp*pid180
2056
      lamp=lonp*pid180
2057
      phipol=latpol*pid180
2058
      lampol=lonpol*pid180
2059
 
2060
      uvtougm=8*Cos(ASin(Cos(lamp)*Cos(phipol)*Cos(phip) +
2061
     -     Sin(phipol)*Sin(phip)))*
2062
     -   (-2*vp*Sin(lamp - phipol) - 2*vp*Sin(lamp + phipol) - 
2063
     -     up*Sin(lamp - phipol - phip) - 
2064
     >     2*up*Sin(phipol - phip) - 
2065
     -     up*Sin(lamp + phipol - phip) + up*Sin(lamp - 
2066
     >     phipol + phip) - 
2067
     -     2*up*Sin(phipol + phip) + up*Sin(lamp + phipol + phip))/
2068
     -  (-20 + 4*Cos(2*lamp) + 2*Cos(2*(lamp - phipol)) - 
2069
     >     4*Cos(2*phipol) + 
2070
     -    2*Cos(2*(lamp + phipol)) - 4*Cos(lamp - 2*phipol - 2*phip) + 
2071
     -    4*Cos(lamp + 2*phipol - 2*phip) + 2*Cos(2*(lamp - phip)) + 
2072
     -    Cos(2*(lamp - phipol - phip)) + 6*Cos(2*(phipol - phip)) + 
2073
     -    Cos(2*(lamp + phipol - phip)) - 4*Cos(2*phip) + 
2074
     -    2*Cos(2*(lamp + phip)) + Cos(2*(lamp - phipol + phip)) + 
2075
     -    6*Cos(2*(phipol + phip)) + Cos(2*(lamp + phipol + phip)) + 
2076
     -    4*Cos(lamp - 2*phipol + 2*phip) - 
2077
     -    4*Cos(lamp + 2*phipol + 2*phip))
2078
 
2079
      return
2080
      end
2081
c-----------------------------------------------------------------
2082
      real function uvtovgm(up,vp,lonp,latp,lonpol,latpol)
2083
c-----------------------------------------------------------------
2084
c     Rene Fehlmann, Wed Nov  1 13:31:28 MET 1995
2085
c     INPUT:
2086
c     - rotated winds (i.e. wrt equatorial grid)
2087
c     - rotated coordinates (i.e. on equatorial grid)
2088
c     - coordinates of the rotated pole
2089
c     OUTPUT
2090
c     - v-component of the wind (i.e. real N-S wind)
2091
c-----------------------------------------------------------------
2092
 
2093
      real               up,vp,lonp,latp,lonpol,latpol
2094
      real               lampol,phipol,lamp,phip
2095
 
2096
      data               pid180 / 0.0174532925199 /
2097
 
2098
c-----------------------------------------------------------------
2099
 
2100
      phip=latp*pid180
2101
      lamp=lonp*pid180
2102
      phipol=latpol*pid180
2103
      lampol=lonpol*pid180
2104
 
2105
      uvtovgm=(vp*Cos(phip)*Sin(phipol) + 
2106
     -     Cos(phipol)*(-(up*Sin(lamp)) -
2107
     -     vp*Cos(lamp)*Sin(phip)))/
2108
     -     Sqrt(1 - (Cos(lamp)*Cos(phipol)*Cos(phip) + 
2109
     -       Sin(phipol)*Sin(phip))**2)
2110
 
2111
      return
2112
      end