WebSVN – pvinversion.ecmwf – Blame – /trunk/pvin/z2s.f

Rev	Author	Line No.	Line
3	michaesp	1	`PROGRAM z2s`
		2
		3	`c Calculate secondary fields on z levels`
		4	`c Michael Sprenger / Summer 2006`
		5
		6	`implicit none`
		7
		8	`c ---------------------------------------------------------------`
		9	`c Declaration of variables`
		10	`c ---------------------------------------------------------------`
		11
		12	`c Variables for output Z file : height level`
		13	`character*80 cfn`
		14	`real varmin(4),varmax(4),stag(4)`
		15	`integer vardim(4)`
		16	`real mdv`
		17	`integer ndim`
		18	`integer nx,ny,nz`
		19	`real xmin,xmax,ymin,ymax,dx,dy`
		20	`integer ntimes`
		21	`real aklev(1000),bklev(1000)`
		22	`real aklay(1000),bklay(1000)`
		23	`real time`
		24	`real pollon,pollat`
		25	`integer idate(5)`
		26	`integer nfields`
		27	`character*80 field_nam(100)`
		28	`real,allocatable, dimension (:,:,:,:) :: field_dat`
		29	`real,allocatable, dimension (:,:,:) :: z3`
		30	`real,allocatable, dimension (:,:) :: x2,y2,f2`
		31	`real,allocatable, dimension (:,:,:) :: out`
		32	`real,allocatable, dimension (:,:,:) :: inp`
		33	`integer nvars`
		34	`character*80 vnam(100),varname`
		35	`integer isok`
		36	`integer cdfid,cstid`
		37
		38	`c Parameter file`
		39	`character*80 fieldname`
		40	`integer nfilt`
		41	`character*80 ofn,gri`
		42
		43	`c Auxiliary variables`
		44	`integer ierr,stat`
		45	`integer i,j,k,n`
		46	`real,allocatable, dimension (:,:) :: out2`
		47	`character*80 vnam2(100)`
		48	`integer nvars2`
		49
		50	`c ---------------------------------------------------------------`
		51	`c Preparations`
		52	`c ---------------------------------------------------------------`
		53
		54	`print,'********************************************************'`
		55	`print,' z2s *'`
		56	`print,'********************************************************'`
		57
		58	`c Read parameter file`
		59	`open(10,file='fort.10')`
		60	`read(10,*) fieldname`
		61	`read(10,*) ofn`
		62	`read(10,*) gri`
		63	`read(10,*) nfilt`
		64	`close(10)`
		65
		66	`c Get grid description for Z file : height level`
		67	`call cdfopn(ofn,cdfid,ierr)`
		68	`if (ierr.ne.0) goto 998`
		69	`call getcfn(cdfid,cfn,ierr)`
		70	`if (ierr.ne.0) goto 998`
		71	`call cdfopn(cfn,cstid,ierr)`
		72	`if (ierr.ne.0) goto 998`
		73	`call getdef(cdfid,'T',ndim,mdv,vardim,`
		74	`> varmin,varmax,stag,ierr)`
		75	`if (ierr.ne.0) goto 998`
		76	`nx =vardim(1)`
		77	`ny =vardim(2)`
		78	`nz =vardim(3)`
		79	`xmin=varmin(1)`
		80	`ymin=varmin(2)`
		81	`call getlevs(cstid,nz,aklev,bklev,aklay,bklay,ierr)`
		82	`if (ierr.ne.0) goto 998`
		83	`call getgrid(cstid,dx,dy,ierr)`
		84	`if (ierr.ne.0) goto 998`
		85	`xmax=xmin+real(nx-1)*dx`
		86	`ymax=ymin+real(ny-1)*dy`
		87	`call gettimes(cdfid,time,ntimes,ierr)`
		88	`if (ierr.ne.0) goto 998`
		89	`call getstart(cstid,idate,ierr)`
		90	`if (ierr.ne.0) goto 998`
		91	`call getpole(cstid,pollon,pollat,ierr)`
		92	`if (ierr.ne.0) goto 998`
		93	`call getvars(cdfid,nvars,vnam,ierr)`
		94	`if (ierr.ne.0) goto 998`
		95	`call clscdf(cstid,ierr)`
		96	`if (ierr.ne.0) goto 998`
		97	`call clscdf(cdfid,ierr)`
		98	`if (ierr.ne.0) goto 998`
		99
		100	`c Get a list of all variables on the GRID file`
		101	`if (trim(gri).ne.trim(ofn)) then`
		102	`call cdfopn(gri,cdfid,ierr)`
		103	`if (ierr.ne.0) goto 998`
		104	`call getvars(cdfid,nvars2,vnam2,ierr)`
		105	`if (ierr.ne.0) goto 998`
		106	`do i=1,nvars2`
		107	`vnam(nvars+i)=vnam2(i)`
		108	`enddo`
		109	`nvars=nvars+nvars2`
		110	`endif`
		111
		112	`c Check whether calculation of <fieldname> is supported`
		113	`if ( (fieldname.ne.'TH' ).and.`
		114	`> (fieldname.ne.'NSQ') .and.`
		115	`> (fieldname.ne.'PV' ) .and.`
		116	`> (fieldname.ne.'RHO')) then`
		117	`print*,'Variable not supported ',trim(fieldname)`
		118	`stop`
		119	`endif`
		120
		121	`c Set dependencies`
		122	`if (fieldname.eq.'TH') then`
		123	`nfields=2`
		124	`field_nam(1)='T'`
		125	`field_nam(2)='P'`
		126	`else if (fieldname.eq.'RHO') then`
		127	`nfields=2`
		128	`field_nam(1)='T'`
		129	`field_nam(2)='P'`
		130	`else if (fieldname.eq.'NSQ') then`
		131	`nfields=2`
		132	`field_nam(1)='T'`
		133	`field_nam(2)='Q'`
		134	`else if (fieldname.eq.'PV') then`
		135	`nfields=4`
		136	`field_nam(1)='T'`
		137	`field_nam(2)='P'`
		138	`field_nam(3)='U'`
		139	`field_nam(4)='V'`
		140	`endif`
		141
		142	`c Allocate memory`
		143	`allocate(field_dat(nfields,nx,ny,nz),stat=stat)`
		144	`if (stat.ne.0) print*,'error allocating field_dat'`
		145	`allocate(out(nx,ny,nz),stat=stat)`
		146	`if (stat.ne.0) print*,'error allocating out'`
		147	`allocate(inp(nx,ny,nz),stat=stat)`
		148	`if (stat.ne.0) print*,'error allocating inp'`
		149	`allocate(z3(nx,ny,nz),stat=stat)`
		150	`if (stat.ne.0) print*,'error allocating z3'`
		151	`allocate(x2(nx,ny),stat=stat)`
		152	`if (stat.ne.0) print*,'error allocating x2'`
		153	`allocate(y2(nx,ny),stat=stat)`
		154	`if (stat.ne.0) print*,'error allocating y2'`
		155	`allocate(f2(nx,ny),stat=stat)`
		156	`if (stat.ne.0) print*,'error allocating f2'`
		157
		158	`c Allocate auxiliary fields`
		159	`allocate(out2(nx,ny),stat=stat)`
		160	`if (stat.ne.0) print*,'error allocating out2'`
		161
		162	`c Read X grid`
		163	`varname='X'`
		164	`isok=0`
		165	`call check_varok (isok,varname,vnam,nvars)`
		166	`if (isok.eq.0) then`
		167	`print*,'Variable ',trim(varname),' missing... Stop'`
		168	`stop`
		169	`endif`
		170	`call cdfopn(gri,cdfid,ierr)`
		171	`if (ierr.ne.0) goto 998`
		172	`call getdat(cdfid,varname,time,0,x2,ierr)`
		173	`if (ierr.ne.0) goto 998`
		174	`print*,'R ',trim(varname),' ',trim(gri)`
		175	`call clscdf(cdfid,ierr)`
		176	`if (ierr.ne.0) goto 998`
		177
		178	`c Read Y grid`
		179	`varname='Y'`
		180	`isok=0`
		181	`call check_varok (isok,varname,vnam,nvars)`
		182	`if (isok.eq.0) then`
		183	`print*,'Variable ',trim(varname),' missing... Stop'`
		184	`stop`
		185	`endif`
		186	`call cdfopn(gri,cdfid,ierr)`
		187	`if (ierr.ne.0) goto 998`
		188	`call getdat(cdfid,varname,time,0,y2,ierr)`
		189	`if (ierr.ne.0) goto 998`
		190	`print*,'R ',trim(varname),' ',trim(gri)`
		191	`call clscdf(cdfid,ierr)`
		192	`if (ierr.ne.0) goto 998`
		193
		194	`c Read Coriolis parametzer`
		195	`varname='CORIOL'`
		196	`isok=0`
		197	`call check_varok (isok,varname,vnam,nvars)`
		198	`if (isok.eq.0) then`
		199	`print*,'Variable ',trim(varname),' missing... Stop'`
		200	`stop`
		201	`endif`
		202	`call cdfopn(gri,cdfid,ierr)`
		203	`if (ierr.ne.0) goto 998`
		204	`call getdat(cdfid,varname,time,0,f2,ierr)`
		205	`if (ierr.ne.0) goto 998`
		206	`print*,'R ',trim(varname),' ',trim(gri)`
		207	`call clscdf(cdfid,ierr)`
		208	`if (ierr.ne.0) goto 998`
		209
		210	`c Init height levels`
		211	`do i=1,nx`
		212	`do j=1,ny`
		213	`do k=1,nz`
		214	`z3(i,j,k)=aklay(k)`
		215	`enddo`
		216	`enddo`
		217	`enddo`
		218
		219	`c Load needed variables`
		220	`do n=1,nfields`
		221
		222	`c Check whether variable is available on file`
		223	`varname=field_nam(n)`
		224	`isok=0`
		225	`call check_varok (isok,varname,vnam,nvars)`
		226
		227	`c Variable is available on file`
		228	`if (isok.eq.1) then`
		229
		230	`call cdfopn(ofn,cdfid,ierr)`
		231	`if (ierr.ne.0) goto 998`
		232	`call getdat(cdfid,varname,time,0,inp,ierr)`
		233	`if (ierr.ne.0) goto 998`
		234	`print*,'R ',trim(varname),' ',trim(ofn)`
		235	`call clscdf(cdfid,ierr)`
		236	`if (ierr.ne.0) goto 998`
		237
		238	`do i=1,nx`
		239	`do j=1,ny`
		240	`do k=1,nz`
		241	`field_dat(n,i,j,k)=inp(i,j,k)`
		242	`enddo`
		243	`enddo`
		244	`enddo`
		245
		246	`else`
		247	`print*,'Variable missing : ',trim(varname)`
		248	`stop`
		249	`endif`
		250
		251	`enddo`
		252
		253	`c Change unit of pressure from hPa to Pa`
		254	`n=0`
		255	`do i=1,nfields`
		256	`if (field_nam(i).eq.'P') n=i`
		257	`enddo`
		258	`if (n.ne.0) then`
		259	`do i=1,nx`
		260	`do j=1,ny`
		261	`do k=1,nz`
		262	`field_dat(n,i,j,k)=100.*field_dat(n,i,j,k)`
		263	`enddo`
		264	`enddo`
		265	`enddo`
		266	`endif`
		267
		268	`c ----------------------------------------------------------------`
		269	`c Calculations`
		270	`c ----------------------------------------------------------------`
		271
		272	`c Call to the defining routines`
		273	`if (fieldname.eq.'RHO') then`
		274
		275	`print*,'C ',trim(fieldname)`
		276	`call calc_rho (out, ! RHO`
		277	`> field_dat(1,:,:,:), ! T`
		278	`> field_dat(2,:,:,:), ! P`
		279	`> nx,ny,nz,mdv)`
		280
		281	`else if (fieldname.eq.'TH') then`
		282
		283	`print*,'C ',trim(fieldname)`
		284	`call calc_theta (out, ! TH`
		285	`> field_dat(1,:,:,:), ! T`
		286	`> field_dat(2,:,:,:), ! P`
		287	`> nx,ny,nz,mdv)`
		288
		289	`else if (fieldname.eq.'NSQ') then`
		290
		291	`print*,'C ',trim(fieldname)`
		292	`call calc_nsq (out, ! NSQ`
		293	`> field_dat(1,:,:,:), ! T`
		294	`> field_dat(2,:,:,:), ! Q`
		295	`> z3, ! Z`
		296	`> nx,ny,nz,mdv)`
		297
		298	`else if (fieldname.eq.'PV') then`
		299
		300	`print*,'C ',trim(fieldname)`
		301	`call calc_pv (out, ! PV`
		302	`> field_dat(1,:,:,:), ! T`
		303	`> field_dat(2,:,:,:), ! P`
		304	`> field_dat(3,:,:,:), ! U`
		305	`> field_dat(4,:,:,:), ! V`
		306	`> z3,x2,y2,f2, ! Z,X,Y,CORIOL`
		307	`> nx,ny,nz,mdv)`
		308
		309	`endif`
		310
		311	`c Horizontal filtering of the resulting fields`
		312	`print*,'F ',trim(fieldname)`
		313	`do k=1,nz`
		314
		315	`do i=1,nx`
		316	`do j=1,ny`
		317	`out2(i,j)=out(i,j,k)`
		318	`enddo`
		319	`enddo`
		320	`do n=1,nfilt`
		321	`call filt2d (out2,out2,nx,ny,1.,mdv,0,0,0,0)`
		322	`enddo`
		323
		324	`do i=1,nx`
		325	`do j=1,ny`
		326	`out(i,j,k)=out2(i,j)`
		327	`enddo`
		328	`enddo`
		329
		330	`enddo`
		331
		332	`c ----------------------------------------------------------------`
		333	`c Save result onto netcdf file`
		334	`c ----------------------------------------------------------------`
		335
		336	`call cdfwopn(ofn,cdfid,ierr)`
		337	`if (ierr.ne.0) goto 998`
		338	`isok=0`
		339	`varname=fieldname`
		340	`call check_varok(isok,varname,vnam,nvars)`
		341	`if (isok.eq.0) then`
		342	`call putdef(cdfid,varname,ndim,mdv,vardim,`
		343	`> varmin,varmax,stag,ierr)`
		344	`if (ierr.ne.0) goto 998`
		345	`endif`
		346	`call putdat(cdfid,varname,time,0,out,ierr)`
		347	`if (ierr.ne.0) goto 998`
		348	`print*,'W ',trim(varname),' ',trim(ofn)`
		349	`call clscdf(cdfid,ierr)`
		350	`if (ierr.ne.0) goto 998`
		351
		352
		353	`c ----------------------------------------------------------------`
		354	`c Exception handling`
		355	`c ----------------------------------------------------------------`
		356
		357	`stop`
		358
		359	`998 print*,'Problem with netcdf file'`
		360	`stop`
		361
		362	`end`
		363
		364
		365
		366	`c ****************************************************************`
		367	`c * SUBROUTINE SECTION: AUXILIARY ROUTINES *`
		368	`c ****************************************************************`
		369
		370	`c ----------------------------------------------------------------`
		371	`c Check whether variable is found on netcdf file`
		372	`c ----------------------------------------------------------------`
		373
		374	`subroutine check_varok (isok,varname,varlist,nvars)`
		375
		376	`c Check whether the variable <varname> is in the list <varlist(nvars)>.`
		377	`c If this is the case, <isok> is incremented by 1. Otherwise <isok>`
		378	`c keeps its value.`
		379
		380	`implicit none`
		381
		382	`c Declaraion of subroutine parameters`
		383	`integer isok`
		384	`integer nvars`
		385	`character*80 varname`
		386	`character*80 varlist(nvars)`
		387
		388	`c Auxiliary variables`
		389	`integer i`
		390
		391	`c Main`
		392	`do i=1,nvars`
		393	`if (trim(varname).eq.trim(varlist(i))) isok=isok+1`
		394	`enddo`
		395
		396	`end`
		397
		398
		399	`c ****************************************************************`
		400	`c * SUBROUTINE SECTION: CALCULATE SECONDARY FIELDS *`
		401	`c ****************************************************************`
		402
		403	`c -----------------------------------------------------------------`
		404	`c Calculate density`
		405	`c -----------------------------------------------------------------`
		406
		407	`subroutine calc_rho (rho3,t3,p3,`
		408	`> nx,ny,nz,mdv)`
		409
		410	`c Calculate the density <rho3> (in kg/m^3) from temperature <t3>`
		411	`c (in deg C) and pressure <p3> (in hPa).`
		412
		413	`implicit none`
		414
		415	`c Declaration of subroutine parameters`
		416	`integer nx,ny,nz`
		417	`real mdv`
		418	`real rho3(nx,ny,nz)`
		419	`real t3 (nx,ny,nz)`
		420	`real p3 (nx,ny,nz)`
		421
		422	`c Declaration of physical constants`
		423	`real eps`
		424	`parameter (eps=0.01)`
		425	`real rd`
		426	`parameter (rd=287.)`
		427	`real tzero`
		428	`parameter (tzero=273.15)`
		429
		430	`c Auxiliary variables`
		431	`integer i,j,k`
		432
		433	`c Calculation`
		434	`do i=1,nx`
		435	`do j=1,ny`
		436	`do k=1,nz`
		437
		438	`if ((abs(t3(i,j,k)-mdv).gt.eps).and.`
		439	`> (abs(p3(i,j,k)-mdv).gt.eps)) then`
		440
		441	`rho3(i,j,k)=1./rd*p3(i,j,k)/(t3(i,j,k)+tzero)`
		442
		443	`else`
		444	`rho3(i,j,k)=mdv`
		445	`endif`
		446
		447	`enddo`
		448	`enddo`
		449	`enddo`
		450
		451	`end`
		452
		453
		454	`c -----------------------------------------------------------------`
		455	`c Calculate potential temperature`
		456	`c -----------------------------------------------------------------`
		457
		458	`subroutine calc_theta (th3,t3,p3,`
		459	`> nx,ny,nz,mdv)`
		460
		461	`c Calculate potential temperature <th3> (in K) from temperature <t3>`
		462	`c (in deg C) and pressure <p3> (in hPa).`
		463
		464	`implicit none`
		465
		466	`c Declaration of subroutine parameters`
		467	`integer nx,ny,nz`
		468	`real mdv`
		469	`real th3(nx,ny,nz)`
		470	`real t3 (nx,ny,nz)`
		471	`real p3 (nx,ny,nz)`
		472
		473	`c Declaration of physical constants`
		474	`real rdcp,tzero,p0`
		475	`parameter (rdcp=0.286)`
		476	`parameter (tzero=273.15)`
		477	`parameter (p0=100000.)`
		478	`real eps`
		479	`parameter (eps=0.01)`
		480
		481	`c Auxiliary variables`
		482	`integer i,j,k`
		483
		484	`c Calculation`
		485	`do i=1,nx`
		486	`do j=1,ny`
		487	`do k=1,nz`
		488
		489	`if ((abs(t3(i,j,k)-mdv).gt.eps).and.`
		490	`> (abs(p3(i,j,k)-mdv).gt.eps)) then`
		491
		492	`th3(i,j,k)=(t3(i,j,k)+tzero)( (p0/p3(i,j,k))*rdcp )`
		493
		494	`else`
		495	`th3(i,j,k)=mdv`
		496	`endif`
		497
		498	`enddo`
		499	`enddo`
		500	`enddo`
		501
		502	`end`
		503
		504	`c -----------------------------------------------------------------`
		505	`c Calculate stratification`
		506	`c -----------------------------------------------------------------`
		507
		508	`subroutine calc_nsq (nsq3,t3,q3,`
		509	`> z3,nx,ny,nz,mdv)`
		510
		511	`c Calculate stratification <nsq3> on the model grid. The grid is`
		512	`c specified in the horizontal by <xmin,ymin,dx,dy,nx,ny>. The number`
		513	`c of vertical levels is <nz>. The input field are: temperature <t3>,`
		514	`c specific humidity <q3>, height <z3> and horizontal grid <x2>, <y2>.`
		515
		516	`implicit none`
		517
		518	`c Declaration of subroutine parameters`
		519	`integer nx,ny,nz`
		520	`real nsq3(nx,ny,nz)`
		521	`real t3 (nx,ny,nz)`
		522	`real q3 (nx,ny,nz)`
		523	`real z3 (nx,ny,nz)`
		524	`real x2 (nx,ny)`
		525	`real y2 (nx,ny)`
		526	`real mdv`
		527
		528	`c Physical and numerical parameters`
		529	`real scale`
		530	`parameter (scale=1.)`
		531	`real g`
		532	`parameter (g=9.80665)`
		533	`real eps`
		534	`parameter (eps=0.01)`
		535	`real tzero`
		536	`parameter (tzero=273.15)`
		537	`real kappa`
		538	`parameter (kappa=0.6078)`
		539	`real cp`
		540	`parameter (cp=1005.)`
		541	`real zerodiv`
		542	`parameter (zerodiv=0.0000000001)`
		543	`real R`
		544	`parameter (R=287.)`
		545
		546	`c Auxiliary variables`
		547	`real tv (nx,ny,nz)`
		548	`real dtvdz(nx,ny,nz)`
		549	`integer i,j,k`
		550	`real scaledz`
		551
		552	`c Calculate 3d virtual temperature`
		553	`do k=1,nz`
		554	`do i=1,nx`
		555	`do j=1,ny`
		556	`if ((abs(t3(i,j,k)-mdv).gt.eps).and.`
		557	`> (abs(q3(i,j,k)-mdv).gt.eps))`
		558	`> then`
		559	`tv(i,j,k) = (t3(i,j,k)+tzero)(1.+kappaq3(i,j,k))`
		560	`else`
		561	`tv(i,j,k) = mdv`
		562	`endif`
		563	`enddo`
		564	`enddo`
		565	`enddo`
		566
		567	`c Vertical derivative of virtual temperature`
		568	`call deriv (dtvdz,tv,'z',z3,x2,y2,nx,ny,nz,mdv)`
		569
		570	`c Stratification`
		571	`do i=1,nx`
		572	`do j=1,ny`
		573	`do k=1,nz`
		574	`if ((abs(dtvdz(i,j,k)-mdv).gt.eps).and.`
		575	`> (abs(tv (i,j,k)-mdv).gt.eps))`
		576	`> then`
		577	`nsq3(i,j,k) = g/tv(i,j,k) * (dtvdz(i,j,k) + g/cp)`
		578	`else`
		579	`nsq3(i,j,k) = mdv`
		580	`endif`
		581	`enddo`
		582	`enddo`
		583	`enddo`
		584
		585	`end`
		586
		587	`c -----------------------------------------------------------------`
		588	`c Calculate potential vorticity`
		589	`c -----------------------------------------------------------------`
		590
		591	`subroutine calc_pv (pv3,t3,p3,u3,v3,`
		592	`> z3,x2,y2,f2,nx,ny,nz,mdv)`
		593
		594	`c Calculate potential vorticity <pv3> on the model grid. The grid is`
		595	`c specified in the horizontal by <xmin,ymin,dx,dy,nx,ny>. The number`
		596	`c of vertical levels is <nz>. The input field are: potential temperature`
		597	`c <th3>, horizontal wind <u3> and <v3>, density <rho3>, height <z3>.`
		598	`c The terms including the vertical velocity are neglected in the calculation`
		599	`c of the PV.`
		600
		601	`implicit none`
		602
		603	`c Declaration of subroutine parameters`
		604	`integer nx,ny,nz`
		605	`real pv3 (nx,ny,nz)`
		606	`real t3 (nx,ny,nz)`
		607	`real p3 (nx,ny,nz)`
		608	`real u3 (nx,ny,nz)`
		609	`real v3 (nx,ny,nz)`
		610	`real z3 (nx,ny,nz)`
		611	`real x2 (nx,ny)`
		612	`real y2 (nx,ny)`
		613	`real f2 (nx,ny)`
		614
		615	`real mdv`
		616
		617	`c Physical and numerical parameters`
		618	`real scale`
		619	`parameter (scale=1.E6)`
		620	`real omega`
		621	`parameter (omega=7.292E-5)`
		622	`real pi180`
		623	`parameter (pi180=3.141592654/180.)`
		624	`real eps`
		625	`parameter (eps=0.01)`
		626
		627	`c Auxiliary variables`
		628	`real dtdz(nx,ny,nz)`
		629	`real dtdx(nx,ny,nz)`
		630	`real dtdy(nx,ny,nz)`
		631	`real dudz(nx,ny,nz)`
		632	`real dvdz(nx,ny,nz)`
		633	`real dvdx(nx,ny,nz)`
		634	`real dudy(nx,ny,nz)`
		635	`real rho3(nx,ny,nz)`
		636	`real th3 (nx,ny,nz)`
		637
		638	`integer i,j,k`
		639
		640	`c Calculate density and potential temperature`
		641	`call calc_rho (rho3,t3,p3,nx,ny,nz,mdv)`
		642	`call calc_theta (th3 ,t3,p3,nx,ny,nz,mdv)`
		643
		644	`c Calculate needed derivatives`
		645	`call deriv (dudz, u3,'z',z3,x2,y2,nx,ny,nz,mdv)`
		646	`call deriv (dvdz, v3,'z',z3,x2,y2,nx,ny,nz,mdv)`
		647	`call deriv (dtdz,th3,'z',z3,x2,y2,nx,ny,nz,mdv)`
		648	`call deriv (dtdx,th3,'x',z3,x2,y2,nx,ny,nz,mdv)`
		649	`call deriv (dtdy,th3,'y',z3,x2,y2,nx,ny,nz,mdv)`
		650	`call deriv (dvdx, v3,'x',z3,x2,y2,nx,ny,nz,mdv)`
		651	`call deriv (dudy, u3,'y',z3,x2,y2,nx,ny,nz,mdv)`
		652
		653	`c Calculate potential vorticity`
		654	`do j=1,ny`
		655	`do i=1,nx`
		656	`do k=1,nz`
		657
		658	`c Evaluate PV formula with missing data check`
		659	`if ((abs(dtdz(i,j,k)-mdv).gt.eps).and.`
		660	`> (abs(dudz(i,j,k)-mdv).gt.eps).and.`
		661	`> (abs(dtdy(i,j,k)-mdv).gt.eps).and.`
		662	`> (abs(dvdz(i,j,k)-mdv).gt.eps).and.`
		663	`> (abs(rho3(i,j,k)-mdv).gt.eps).and.`
		664	`> (abs(dtdx(i,j,k)-mdv).gt.eps).and.`
		665	`> (abs(dvdx(i,j,k)-mdv).gt.eps).and.`
		666	`> (abs(dudy(i,j,k)-mdv).gt.eps)) then`
		667
		668	`pv3(i,j,k)=scale1./rho3(i,j,k)`
		669	`> ((dvdx(i,j,k)-dudy(i,j,k)+f2(i,j))*dtdz(i,j,k)`
		670	`> +dudz(i,j,k)*dtdy(i,j,k)`
		671	`> -dvdz(i,j,k)*dtdx(i,j,k))`
		672	`else`
		673	`pv3(i,j,k)=mdv`
		674	`endif`
		675
		676	`enddo`
		677	`enddo`
		678	`enddo`
		679
		680	`end`
		681
		682
		683	`c ****************************************************************`
		684	`c * SUBROUTINE SECTION: GRID HANDLING *`
		685	`c ****************************************************************`
		686
		687	`c -----------------------------------------------------------------`
		688	`c Horizontal and vertical derivatives for 3d fields`
		689	`c -----------------------------------------------------------------`
		690
		691	`subroutine deriv (df,f,direction,`
		692	`> z3,x2,y2,nx,ny,nz,mdv)`
		693
		694	`c Calculate horizontal and vertical derivatives of the 3d field <f>.`
		695	`c The direction of the derivative is specified in <direction>`
		696	`c 'x','y' : Horizontal derivative in x and y direction`
		697	`c 'p','z','t','m' : Vertical derivative (pressure, height, theta, model)`
		698	`c The 3d field <z3> specifies the isosurfaces along which the horizontal`
		699	`c derivatives are calculated or the levels for the vertical derivatives.`
		700
		701	`implicit none`
		702
		703	`c Input and output parameters`
		704	`integer nx,ny,nz`
		705	`real df (nx,ny,nz)`
		706	`real f (nx,ny,nz)`
		707	`real z3 (nx,ny,nz)`
		708	`real x2 (nx,ny)`
		709	`real y2 (nx,ny)`
		710	`character direction`
		711	`real mdv`
		712
		713	`c Numerical and physical parameters`
		714	`real pi180`
		715	`parameter (pi180=3.141592654/180.)`
		716	`real zerodiv`
		717	`parameter (zerodiv=0.00000001)`
		718	`real eps`
		719	`parameter (eps=0.01)`
		720
		721	`c Auxiliary variables`
		722	`integer i,j,k`
		723	`real vmin,vmax`
		724	`real scale,lat`
		725	`real vu,vl,vuvl,vlvu`
		726	`integer o,u,w,e,n,s`
		727
		728	`c Vertical derivative`
		729	`if ((direction.eq.'z').or.`
		730	`> (direction.eq.'th').or.`
		731	`> (direction.eq.'p').or.`
		732	`> (direction.eq.'m').and.`
		733	`> (nz.gt.1)) then`
		734
		735	`do i=1,nx`
		736	`do j=1,ny`
		737	`do k=1,nz`
		738
		739	`o=k+1`
		740	`if (o.gt.nz) o=nz`
		741	`u=k-1`
		742	`if (u.lt.1) u=1`
		743
		744	`if ((abs(f(i,j,o)-mdv).gt.eps).and.`
		745	`> (abs(f(i,j,u)-mdv).gt.eps).and.`
		746	`> (abs(f(i,j,k)-mdv).gt.eps)) then`
		747
		748	`vu = z3(i,j,k)-z3(i,j,o)`
		749	`vl = z3(i,j,u)-z3(i,j,k)`
		750	`vuvl = vu/(vl+zerodiv)`
		751	`vlvu = 1./(vuvl+zerodiv)`
		752
		753	`df(i,j,k) = 1./(vu+vl)`
		754	`> * (vuvl*(f(i,j,u)-f(i,j,k))`
		755	`> + vlvu*(f(i,j,k)-f(i,j,o)))`
		756
		757	`else`
		758	`df(i,j,k) = mdv`
		759	`endif`
		760
		761	`enddo`
		762	`enddo`
		763	`enddo`
		764
		765	`c Horizontal derivative in the y direction: 3d`
		766	`elseif (direction.eq.'y') then`
		767
		768	`do i=1,nx`
		769	`do j=1,ny`
		770	`do k=1,nz`
		771
		772	`s=j-1`
		773	`if (s.lt.1) s=1`
		774	`n=j+1`
		775	`if (n.gt.ny) n=ny`
		776
		777	`if ((abs(f(i,n,k)-mdv).gt.eps).and.`
		778	`> (abs(f(i,j,k)-mdv).gt.eps).and.`
		779	`> (abs(f(i,s,k)-mdv).gt.eps)) then`
		780
		781	`vu = 1000.*(y2(i,j)-y2(i,n))`
		782	`vl = 1000.*(y2(i,s)-y2(i,j))`
		783	`vuvl = vu/(vl+zerodiv)`
		784	`vlvu = 1./(vuvl+zerodiv)`
		785
		786	`df(i,j,k) = 1./(vu+vl)`
		787	`> * (vuvl*(f(i,s,k)-f(i,j,k))`
		788	`> + vlvu*(f(i,j,k)-f(i,n,k)))`
		789
		790	`else`
		791	`df(i,j,k) = mdv`
		792	`endif`
		793
		794	`enddo`
		795	`enddo`
		796	`enddo`
		797
		798	`c Horizontal derivative in the x direction: 3d`
		799	`elseif (direction.eq.'x') then`
		800
		801	`do i=1,nx`
		802	`do j=1,ny`
		803	`do k=1,nz`
		804
		805	`w=i-1`
		806	`if (w.lt.1) w=1`
		807	`e=i+1`
		808	`if (e.gt.nx) e=nx`
		809
		810	`if ((abs(f(w,j,k)-mdv).gt.eps).and.`
		811	`> (abs(f(i,j,k)-mdv).gt.eps).and.`
		812	`> (abs(f(e,j,k)-mdv).gt.eps)) then`
		813
		814	`vu = 1000.*(x2(i,j)-x2(e,j))`
		815	`vl = 1000.*(x2(w,j)-x2(i,j))`
		816	`vuvl = vu/(vl+zerodiv)`
		817	`vlvu = 1./(vuvl+zerodiv)`
		818
		819	`df(i,j,k) = 1./(vu+vl)`
		820	`> * (vuvl*(f(w,j,k)-f(i,j,k))`
		821	`> + vlvu*(f(i,j,k)-f(e,j,k)))`
		822
		823	`else`
		824	`df(i,j,k) = mdv`
		825	`endif`
		826
		827	`enddo`
		828	`enddo`
		829	`enddo`
		830
		831	`c Undefined direction for derivative`
		832	`else`
		833
		834	`print*,'Invalid direction of derivative... Stop'`
		835	`stop`
		836
		837	`endif`
		838
		839	`end`
		840
		841	`c -----------------------------------------------------------------`
		842	`c Horizontal filter`
		843	`c -----------------------------------------------------------------`
		844
		845	`subroutine filt2d (a,af,nx,ny,fil,misdat,`
		846	`& iperx,ipery,ispol,inpol)`
		847
		848	`c Apply a conservative diffusion operator onto the 2d field a,`
		849	`c with full missing data checking.`
		850	`c`
		851	`c a real inp array to be filtered, dimensioned (nx,ny)`
		852	`c af real out filtered array, dimensioned (nx,ny), can be`
		853	`c equivalenced with array a in the calling routine`
		854	`c f1 real workarray, dimensioned (nx+1,ny)`
		855	`c f2 real workarray, dimensioned (nx,ny+1)`
		856	`c fil real inp filter-coeff., 0<afil<=1. Maximum filtering with afil=1`
		857	`c corresponds to one application of the linear filter.`
		858	`c misdat real inp missing-data value, a(i,j)=misdat indicates that`
		859	`c the corresponding value is not available. The`
		860	`c misdat-checking can be switched off with with misdat=0.`
		861	`c iperx int inp periodic boundaries in the x-direction (1=yes,0=no)`
		862	`c ipery int inp periodic boundaries in the y-direction (1=yes,0=no)`
		863	`c inpol int inp northpole at j=ny (1=yes,0=no)`
		864	`c ispol int inp southpole at j=1 (1=yes,0=no)`
		865	`c`
		866	`c Christoph Schaer, 1993`
		867
		868	`c argument declaration`
		869	`integer nx,ny`
		870	`real a(nx,ny),af(nx,ny),fil,misdat`
		871	`integer iperx,ipery,inpol,ispol`
		872
		873	`c local variable declaration`
		874	`integer i,j,is`
		875	`real fh`
		876	`real f1(nx+1,ny),f2(nx,ny+1)`
		877
		878	`c compute constant fh`
		879	`fh=0.125*fil`
		880
		881	`c compute fluxes in x-direction`
		882	`if (misdat.eq.0.) then`
		883	`do j=1,ny`
		884	`do i=2,nx`
		885	`f1(i,j)=a(i-1,j)-a(i,j)`
		886	`enddo`
		887	`enddo`
		888	`else`
		889	`do j=1,ny`
		890	`do i=2,nx`
		891	`if ((a(i,j).eq.misdat).or.(a(i-1,j).eq.misdat)) then`
		892	`f1(i,j)=0.`
		893	`else`
		894	`f1(i,j)=a(i-1,j)-a(i,j)`
		895	`endif`
		896	`enddo`
		897	`enddo`
		898	`endif`
		899	`if (iperx.eq.1) then`
		900	`c do periodic boundaries in the x-direction`
		901	`do j=1,ny`
		902	`f1(1,j)=f1(nx,j)`
		903	`f1(nx+1,j)=f1(2,j)`
		904	`enddo`
		905	`else`
		906	`c set boundary-fluxes to zero`
		907	`do j=1,ny`
		908	`f1(1,j)=0.`
		909	`f1(nx+1,j)=0.`
		910	`enddo`
		911	`endif`
		912
		913	`c compute fluxes in y-direction`
		914	`if (misdat.eq.0.) then`
		915	`do j=2,ny`
		916	`do i=1,nx`
		917	`f2(i,j)=a(i,j-1)-a(i,j)`
		918	`enddo`
		919	`enddo`
		920	`else`
		921	`do j=2,ny`
		922	`do i=1,nx`
		923	`if ((a(i,j).eq.misdat).or.(a(i,j-1).eq.misdat)) then`
		924	`f2(i,j)=0.`
		925	`else`
		926	`f2(i,j)=a(i,j-1)-a(i,j)`
		927	`endif`
		928	`enddo`
		929	`enddo`
		930	`endif`
		931	`c set boundary-fluxes to zero`
		932	`do i=1,nx`
		933	`f2(i,1)=0.`
		934	`f2(i,ny+1)=0.`
		935	`enddo`
		936	`if (ipery.eq.1) then`
		937	`c do periodic boundaries in the x-direction`
		938	`do i=1,nx`
		939	`f2(i,1)=f2(i,ny)`
		940	`f2(i,ny+1)=f2(i,2)`
		941	`enddo`
		942	`endif`
		943	`if (iperx.eq.1) then`
		944	`if (ispol.eq.1) then`
		945	`c do south-pole`
		946	`is=(nx-1)/2`
		947	`do i=1,nx`
		948	`f2(i,1)=-f2(mod(i-1+is,nx)+1,2)`
		949	`enddo`
		950	`endif`
		951	`if (inpol.eq.1) then`
		952	`c do north-pole`
		953	`is=(nx-1)/2`
		954	`do i=1,nx`
		955	`f2(i,ny+1)=-f2(mod(i-1+is,nx)+1,ny)`
		956	`enddo`
		957	`endif`
		958	`endif`
		959
		960	`c compute flux-convergence -> filter`
		961	`if (misdat.eq.0.) then`
		962	`do j=1,ny`
		963	`do i=1,nx`
		964	`af(i,j)=a(i,j)+fh*(f1(i,j)-f1(i+1,j)+f2(i,j)-f2(i,j+1))`
		965	`enddo`
		966	`enddo`
		967	`else`
		968	`do j=1,ny`
		969	`do i=1,nx`
		970	`if (a(i,j).eq.misdat) then`
		971	`af(i,j)=misdat`
		972	`else`
		973	`af(i,j)=a(i,j)+fh*(f1(i,j)-f1(i+1,j)+f2(i,j)-f2(i,j+1))`
		974	`endif`
		975	`enddo`
		976	`enddo`
		977	`endif`
		978	`end`

Subversion Repositories pvinversion.ecmwf

(root)/trunk/pvin/z2s.f @ 4 – Rev 3