Subversion Repositories tropofold.echam

Compare Revisions

No changes between revisions

Ignore whitespace Rev 4 → Rev 5

/tags/1.0/LICENSE.pdf
Cannot display: file marked as a binary type.
svn:mime-type = application/pdf
Property changes:
Added: svn:mime-type
+application/pdf
\ No newline at end of property
/tags/1.0/3d_labelling_and_fold_id.f90
0,0 → 1,1492
PROGRAM clustering
 
! *********************************************************************************
! * Interpolate fields to tropopause and assign labels with clustering algorithm *
! * Michael Sprenger / Spring 2006: initial version *
! * Michael Sprenger / 2012: horizontal propagation of LABEL 2 prohibited in *
! * lower troposphere *
! * Bojan Skerlak / Spring 2012 / Fall 2012: advanced LABEL 1-5 handling *
! * Jan 2013: - change horiz. prop. decision method to height(log(pressure)) *
! * - perform check of label for exotic cases where the pressure *
! * is higher than 800 hPa (deep folds) *
! * - additional cutoff-special-treatment prohibits identification *
! * of strat. air as cutoffs if below horizontal prop. limit *
! * (e.g. deep 'funnels' with 'bumps' on the side) *
! * - additional check in strat. (TH>380) for 'wrong' sign of PV *
! * that can lead to errors in calculating TROPO *
! * - additional check for PV towers (whole column gets label 2) *
! * July 2013: - replace fixed tropo_pv,tropo_th by user-defined values *
! * Jan 2014: - v2.0 Add specific humidity criterion for diabatic PV *
! * Feb 2014: - v2.1 Correct 'filling' when labels 2/3 or 2/5 meet *
! * Mar 2014: - v2.2 Filling of 2/3 and 2/5 also propagates downward via *
! * 'connect' subroutine. Thus, now also dird/diru have to be *
! * specified. *
! * - v2.3 Identify and write out tropopause folds (FOLD) 1/0 *
! * Nov 2014: - remove some commented sections and delete temporary stuff *
! * Dec 2014: - simply the code to only contain the core algorithm and no *
! * input/output section (which depends on data format and *
! * model setup) *
! * Andrea Pozzer (MPIC):
! * Dec 2014: - Implementation for netcdf I/O, Makefile and namelist reading *
! *********************************************************************************
!????????????????????????????????????????????????????????????????????????????????
!????????????????????????????????????????????????????????????????????????????????
!????????????????????????????????????????????????????????????????????????????????
!??????????????? ??????????????????????????
!??????????????? 1) 3D-LABELLING ??????????????????????????
!??????????????? 2) TROPOPAUSE FOLD IDENTIFICATION ??????????????????????????
!??????????????? ??????????????????????????
!????????????????????????????????????????????????????????????????????????????????
!????????????????????????????????????????????????????????????????????????????????
!MMM?????????????????????????????????????????????????????????????????????????????
!.....DMMI???????????????????????????????????????????MMMMMMMMMMM?I???????????????
!.........IMO????????????????????????????????????IMM ........... 8MM8???????????
!............MM?????????????????????????????????M$ ..............MMMI??????
!............. MD?????????????????????????????IM. PMIN .................OMMI??
!................MI???????????????????????????M . ^ ......................NM
!.................MO?????????????????????????M... | .............................
!.................. M???????????????????????M.... | .............................
!....................M8?????????????????????M.... | .........................
!.....................:M????????????????????M ... | DP .........................
!.......................MO???????????????????M... | .........................
!........................ M??????????????????M~.. | .............................
!..........................MM????????????????IM.. | .............................
!............................MI????????????????M. v .............................
!.............................:MI???????????????M. .............................
!...............................DMI??????????????M$ .............................
!................................ MM??????????????DM ............................
!.................................. MMI????????????IM ...........................
!.....................................MMI????????????M ..........................
!........................................MM??????????IM..........................
!......................................... $MMI???????M .........................
!........................................... MMMMMMM7. ..........................
!............................................. .............................
!............................................. PMAX .............................
!............................................. .............................
!................................................................................
 
USE mo_f2kcli ! command line interface
USE netcdf_tools
 
 
implicit none
 
! VERSION
character(len=*), parameter :: VERSION = '1.4'
! ... INTERNAL parameter
integer, parameter :: str_short = 30 ! length of short strings
integer, parameter :: str_long = 100 ! length of long strings
integer, parameter :: str_vlong = 500 ! length of very long strings
integer, parameter :: iou = 21 ! I/O unit
 
! FOR COMMAND LINE
CHARACTER(LEN=256) :: EXE ! program name
CHARACTER(LEN=80) :: CMD ! argument
INTEGER :: NARG ! number of arguments
 
! VARIABLES
INTEGER :: status ! status flag
LOGICAL :: l_invert = .FALSE. ! invert axis
LOGICAL :: l_Pa = .FALSE. ! convert units Pa->hPa
! Set tropopause thresholds
real :: tropo_pv,tropo_th
! diabatically produced PV (Q-threshold), Feb 2014:
real :: q_th
! masking number
real, parameter :: mdv =-999.0
 
! DATA/VARIABLE
! - input data
real, allocatable, dimension(:) :: x_data, y_data,z_data,t_data
real, allocatable, dimension(:) :: ak,bk
real, allocatable, dimension(:,:,:) :: aps
real, allocatable, dimension(:,:,:,:) :: press
real, allocatable, dimension(:,:,:,:) :: q_in,pv_in,pt_in
real, allocatable, dimension(:,:,:,:) :: q,pv,pt
character(len=str_long) :: x_units,y_units, z_units,t_units
character(len=str_long) :: aps_units
! - final results
real, allocatable, dimension(:,:,:) :: fold ! fold yes/no
real, allocatable, dimension(:,:,:) :: sfold ! shallow fold yes/no
real, allocatable, dimension(:,:,:) :: mfold ! medium fold yes/no
real, allocatable, dimension(:,:,:) :: dfold ! deep fold yes/no
real, allocatable, dimension(:,:,:) :: tp ! tropopause
real, allocatable, dimension(:,:,:) :: dp ! folding depth (Pa)
real, allocatable, dimension(:,:,:) :: pmin ! folding depth (Pa)
real, allocatable, dimension(:,:,:) :: pmax ! folding depth (Pa)
real, allocatable, dimension(:,:,:,:) :: label ! folding classification
real, allocatable, dimension(:,:,:,:) :: label_out
 
! Auxiliary variables
integer :: stat
integer :: i,j,k,t
real :: lev(100)
 
 
!namelist...
character(len=str_long) file_input
character(len=str_long) file_output
character(len=str_short) x_name,y_name,z_name,t_name
character(len=str_short) aps_name,hyam_name, hybm_name
character(len=str_short) q_name, pv_name, pt_name
 
! DATA DIMENSION
integer :: nx,ny,nz,nt
 
 
! -----------------------------------------------------------------
! Initialisation paramaters
! -----------------------------------------------------------------
! Humidity criterion for diabatically produced PV, BOJAN: this is the threshold Q=0.1g/kg
q_th=0.0001
tropo_pv = 2.0
tropo_th = 380.0
 
! -----------------------------------------------------------------
! Read command line
! -----------------------------------------------------------------
narg = command_argument_count() ! number of arguments
call get_command_argument(0,exe) ! program name
!
if (narg > 1) then
write(*,*) 'command-line error: too many arguments !'
call usage(trim(exe))
stop
end if
!
if (narg == 0) then
call usage(trim(exe))
stop
end if
!
call get_command_argument(1,cmd)
 
! -----------------------------------------------------------------
! Read namelist file
! -----------------------------------------------------------------
 
CALL read_nml(status, iou, TRIM(CMD))
IF (status /= 0) STOP
 
! -----------------------------------------------------------------
! Inquire netcdf file
! -----------------------------------------------------------------
CALL inquire_file(TRIM(file_input), &
x_name, y_name, z_name,t_name, &
nx, ny, nz, nt )
! IF (status > 0) STOP ! ERROR
 
 
WRITE(*,*) '==========================================================='
WRITE(*,*) " FILE DIMENSION: "
WRITE(*,*) " X :", nx
WRITE(*,*) " Y :", ny
WRITE(*,*) " Z :", nz
WRITE(*,*) " T :", nt
WRITE(*,*) '==========================================================='
 
! -----------------------------------------------------------------
! allocate data file
! -----------------------------------------------------------------
ALLOCATE(x_data(nx))
ALLOCATE(y_data(ny))
ALLOCATE(z_data(nz))
ALLOCATE(t_data(nt))
 
ALLOCATE(aps(nx,ny,nt))
ALLOCATE(ak(nz))
ALLOCATE(bk(nz))
ALLOCATE(press(nx,ny,nz,nt))
!
ALLOCATE(q(nx,ny,nz,nt))
ALLOCATE(pt(nx,ny,nz,nt))
ALLOCATE(pv(nx,ny,nz,nt))
ALLOCATE(q_in(nx,ny,nz,nt))
ALLOCATE(pt_in(nx,ny,nz,nt))
ALLOCATE(pv_in(nx,ny,nz,nt))
 
ALLOCATE(label(nx,ny,nz,nt))
ALLOCATE(fold(nx,ny,nt))
ALLOCATE(sfold(nx,ny,nt))
ALLOCATE(mfold(nx,ny,nt))
ALLOCATE(dfold(nx,ny,nt))
ALLOCATE(tp(nx,ny,nt))
ALLOCATE(dp(nx,ny,nt))
ALLOCATE(pmin(nx,ny,nt))
ALLOCATE(pmax(nx,ny,nt))
ALLOCATE(label_out(nx,ny,nz,nt))
 
! -----------------------------------------------------------------
! read netcdf file
! -----------------------------------------------------------------
call read_file(TRIM(file_input), TRIM(x_name), x_data)
call read_file(TRIM(file_input), TRIM(y_name), y_data)
call read_file(TRIM(file_input), TRIM(z_name), z_data)
call read_file(TRIM(file_input), TRIM(t_name), t_data)
 
call read_att(TRIM(file_input), TRIM(x_name), "units", x_units)
call read_att(TRIM(file_input), TRIM(y_name), "units", y_units)
call read_att(TRIM(file_input), TRIM(z_name), "units", z_units)
call read_att(TRIM(file_input), TRIM(t_name), "units", t_units)
 
! pressure data
call read_file(TRIM(file_input), TRIM(aps_name), aps)
call read_file(TRIM(file_input), TRIM(hyam_name), ak)
call read_file(TRIM(file_input), TRIM(hybm_name), bk)
! presure units
call read_att(TRIM(file_input), TRIM(aps_name), "units", aps_units)
 
! humidity
call read_file(TRIM(file_input), TRIM(q_name), q_in)
! potential temperature
call read_file(TRIM(file_input), TRIM(pt_name), pt_in)
! potential vorticity
call read_file(TRIM(file_input), TRIM(pv_name), pv_in)
 
! -----------------------------------------------------------------
! create pressure field and
! invert vertical axis (if needed)
! -----------------------------------------------------------------
! bottom => high pressure = 1
! top => low pressure = nz
if ((ak(1)+bk(1)*101325).lt.(ak(2)+bk(2)*101325)) then
l_invert=.TRUE.
write (*,*) 'Invert vertical axis'
endif
 
if (aps_units == "Pa") then
write (*,*) 'convert pressure from Pa to hPa'
l_Pa=.TRUE.
endif
if (l_invert) then
! to be inverted
do i=1,nx
do j=1,ny
do k=1,nz
do t=1,nt
if (l_Pa) then
press(i,j,nz-k+1,t)= (ak(k)+aps(i,j,t)*bk(k))/100.0 ! Pa -> hPa
else
press(i,j,nz-k+1,t)= (ak(k)+aps(i,j,t)*bk(k))
endif
q(i,j,nz-k+1,t) = q_in(i,j,k,t)
pt(i,j,nz-k+1,t)= pt_in(i,j,k,t)
pv(i,j,nz-k+1,t)= pv_in(i,j,k,t)
enddo
enddo
enddo
enddo
ELSE
! NOT to be inverted
do i=1,nx
do j=1,ny
do k=1,nz
do t=1,nt
if (l_Pa) then
press(i,j,k,t)= (ak(k)+aps(i,j,t)*bk(k))/100.0 ! Pa -> hPa
else
press(i,j,k,t)= ak(k)+aps(i,j,t)*bk(k)
endif
enddo
enddo
enddo
enddo
q(1:nx,1:ny,1:nz,1:nt) = q_in(1:nx,1:ny,1:nz,1:nt)
pt(1:nx,1:ny,1:nz,1:nt) = pt_in(1:nx,1:ny,1:nz,1:nt)
pv(1:nx,1:ny,1:nz,1:nt) = pv_in(1:nx,1:ny,1:nz,1:nt)
ENDIF
 
! -----------------------------------------------------------------
! perform calculation for each timestep
! -----------------------------------------------------------------
! TIME LOOP
do t=1,nt
 
write(*,*) "time step",t,"/", nt
! -------------------------------------------------------------------
! Part 1) Run clustering algorithm, get labels and interpolate to tropopause
! -------------------------------------------------------------------
 
! ------------------------------
! Running clustering algorithm
! ------------------------------
call tropopause (tp(:,:,t),label(:,:,:,t),press(:,:,:,t),pv(:,:,:,t),pt(:,:,:,t),q(:,:,:,t), &
mdv,y_data(:),nx,ny,nz,tropo_pv,tropo_th,q_th)
 
! -------------------------------------------------------------------
! Part 2) Identify tropopause folds from vertical cross-sections
! -----------------------------------------------------------------
! ------------------------------
! Identifying tropopause folds
! ------------------------------
call tropofold (label(:,:,:,t),press(:,:,:,t),pv(:,:,:,t),pt(:,:,:,t), &
fold(:,:,t), dp(:,:,t), pmin(:,:,t), pmax(:,:,t), sfold(:,:,t), mfold(:,:,t), dfold(:,:,t), &
mdv,nx,ny,nz,tropo_pv,tropo_th)
enddo
! END OF TIME LOOP
 
! -----------------------------------------------------------------
! re-invert vertical axis (if needed)
! -----------------------------------------------------------------
if (l_invert) then
! to be inverted
do i=1,nx
do j=1,ny
do k=1,nz
do t=1,nt
label_out(i,j,nz-k+1,t) = label(i,j,k,t)
enddo
enddo
enddo
enddo
ELSE
! NOT to be inverted
label_out(1:nx,1:ny,1:nz,1:nt) = label(1:nx,1:ny,1:nz,1:nt)
ENDIF
 
! -----------------------------------------------------------------
! write netcdf files
! -----------------------------------------------------------------
call nc_dump (TRIM(file_output), x_data, y_data, z_data, t_data, &
x_units,y_units, z_units,t_units, aps, ak, bk, label_out, fold, &
sfold, mfold, dfold, tp, dp, pmin, pmax)
 
! -----------------------------------------------------------------
! deallocate data file
! -----------------------------------------------------------------
DEALLOCATE(x_data)
DEALLOCATE(y_data)
DEALLOCATE(z_data)
DEALLOCATE(t_data)
 
DEALLOCATE(aps)
DEALLOCATE(ak)
DEALLOCATE(bk)
DEALLOCATE(press)
!
DEALLOCATE(q_in)
DEALLOCATE(pt_in)
DEALLOCATE(pv_in)
DEALLOCATE(q)
DEALLOCATE(pt)
DEALLOCATE(pv)
 
DEALLOCATE(label)
DEALLOCATE(label_out)
DEALLOCATE(fold)
DEALLOCATE(sfold)
DEALLOCATE(mfold)
DEALLOCATE(dfold)
DEALLOCATE(tp)
DEALLOCATE(dp)
DEALLOCATE(pmin)
DEALLOCATE(pmax)
 
CONTAINS
 
! ------------------------------------------------------------------------
SUBROUTINE USAGE(EXE)
CHARACTER (LEN=*) :: EXE
WRITE(*,*) '--------------------------------------------'
WRITE(*,*) '3d_labelling_and_fold_id, version: ',VERSION
WRITE(*,*) 'Authors: Andrea Pozzer, MPICH '
WRITE(*,*) ' Bojan Skerlak, ETH '
WRITE(*,*) ' Michael Sprenger, ETH '
WRITE(*,*) '--------------------------------------------'
WRITE(*,*) 'Usage: '//TRIM(EXE)//' <namelist-file>'
WRITE(*,*) 'See README.txt or EMAC.nml for example'
WRITE(*,*) '--------------------------------------------'
END SUBROUTINE USAGE
! ------------------------------------------------------------------------
 
! ------------------------------------------------------------------------
SUBROUTINE read_nml(status, iou, fname)
 
IMPLICIT NONE
 
! I/O
INTEGER, INTENT(OUT) :: status
INTEGER, INTENT(IN) :: iou
CHARACTER(LEN=*), INTENT(IN) :: fname
 
! LOCAL
CHARACTER(LEN=*), PARAMETER :: substr = 'read_nml'
LOGICAL :: lex ! file exists
INTEGER :: fstat ! file status
 
NAMELIST /CTRL/ file_input,file_output,X_name,Y_name,Z_name,T_name, &
APS_name,HYAM_name,HYBM_name,Q_name,PV_name,PT_name
 
status = 1 ! ERROR
 
WRITE(*,*) '==========================================================='
 
! CHECK IF FILE EXISTS
INQUIRE(file=TRIM(fname), exist=lex)
IF (.NOT.lex) THEN
WRITE(*,*) substr,': FILE DOES NOT EXIST (',TRIM(fname),')'
status = 1
RETURN
END IF
 
! OPEN FILE
OPEN(iou,file=TRIM(fname))
! READ NEMELIST
WRITE(*,*) 'READING NAMELIST ''CTRL'''//&
&' FROM '''//TRIM(fname),''' (unit ',iou,') ...'
!
READ(iou, NML=CTRL, IOSTAT=fstat)
!
IF (fstat /= 0) THEN
WRITE(*,*) substr,': READ ERROR IN NAMELIST ''CTRL'' (',TRIM(fname),')'
status = 3 ! READ ERROR IN NAMELIST
RETURN
END IF
 
WRITE(*,*) ' FILE INPUT : ', TRIM(file_input)
WRITE(*,*) ' FILE OUTPUT : ', TRIM(file_output)
WRITE(*,*) ' X_NAME : ', TRIM(x_name)
WRITE(*,*) ' Y_NAME : ', TRIM(y_name)
WRITE(*,*) ' Z_NAME : ', TRIM(z_name)
WRITE(*,*) ' T_NAME : ', TRIM(t_name)
WRITE(*,*) ' APS_NAME : ', TRIM(aps_name)
WRITE(*,*) ' HYAM_NAME : ', TRIM(hyam_name)
WRITE(*,*) ' HYBM_NAME : ', TRIM(hybm_name)
WRITE(*,*) ' Q_NAME : ', TRIM(q_name)
WRITE(*,*) ' PV_NAME : ', TRIM(pv_name)
WRITE(*,*) ' PT_NAME : ', TRIM(pt_name)
 
! CLOSE FILE
CLOSE(iou)
 
WRITE(*,*) '==========================================================='
 
status = 0
 
END SUBROUTINE read_nml
! ------------------------------------------------------------------------
 
 
end program clustering
 
 
! *******************************************************************************
! * Subroutine Section *
! *******************************************************************************
 
 
! --------------------------------------------------------
! Get grid points which are connected to a grid point: BOJAN: this is the 3D-connectedness criterion of the 3D-labelling algorithm
! --------------------------------------------------------
 
SUBROUTINE connect (outar,label,inx,iny,inz,dirh,dird,diru,nx,ny,nz)
 
! The input array <outar(nx,ny,nz)> is 0/1 labeled. Given a starting point
! <inx,iny,inz> where the label is 1, find all points with label 1 which are
! connected and attribute to all these points the label <label>.
! The 3D arrays <dir*(nx,ny,nz)> specifie whether an expansion is allowed
! dirh=1 where horizontal propagation is allowed,
! dird=1 where vertical propagation is allowed downward,
! diru=1 where vertical propagation is allowed upward, and
! a value of dir*=0 prohibits the respective propagation
 
! Declaration of subroutine parameters
integer :: nx,ny,nz
integer :: inx,iny,inz
integer :: outar(nx,ny,nz)
integer :: label
integer :: dirh(nx,ny,nz)
integer :: diru(nx,ny,nz)
integer :: dird(nx,ny,nz)
 
! Auxiliary variables
integer :: il,ir,jb,jf,ku,kd,im,jm,km
integer :: i,j,k
integer :: stack
integer :: indx(nx*ny*nz),indy(nx*ny*nz),indz(nx*ny*nz)
 
! Push the starting elements on the stack
stack=1
indx(stack)=inx
indy(stack)=iny
indz(stack)=inz
outar(inx,iny,inz)=label
 
! Define the indices of the neighbouring elements
100 continue
 
il=indx(stack)-1
if (il.lt.1) il=nx
ir=indx(stack)+1
if (ir.gt.nx) ir=1
jb=indy(stack)-1
if (jb.lt.1) jb=1
jf=indy(stack)+1
if (jf.gt.ny) jf=ny
ku=indz(stack)+1
if (ku.gt.nz) ku=nz
kd=indz(stack)-1
if (kd.lt.1) kd=1
im=indx(stack)
jm=indy(stack)
km=indz(stack)
stack=stack-1
 
! Check for index overflow
if (stack.ge.(nx*ny*nz-nx)) then
print*,'Stack overflow while clustering...'
stop
endif
 
! Mark all connected elements (build up the stack)
! Mark level below/above if dird/diru and dirh allowed
! BOJAN: U=up, D=down, R=right, L=left, M=mid, F=front, B=back (positions in a 27-element cube of nearest neighbours)
 
if ((dirh(im,jm,km).eq.1).and.(dird(im,jm,km).eq.1)) then
! below:
if ( outar(il,jf,kd).eq.1) then
outar(il,jf,kd)=label
stack=stack+1
indx(stack)=il
indy(stack)=jf
indz(stack)=kd
endif
if ( outar(im,jf,kd).eq.1) then
outar(im,jf,kd)=label
stack=stack+1
indx(stack)=im
indy(stack)=jf
indz(stack)=kd
endif
if ( outar(ir,jf,kd).eq.1) then
outar(ir,jf,kd)=label
stack=stack+1
indx(stack)=ir
indy(stack)=jf
indz(stack)=kd
endif
if (outar(il,jm,kd).eq.1) then
outar(il,jm,kd)=label
stack=stack+1
indx(stack)=il
indy(stack)=jm
indz(stack)=kd
endif
if (outar(ir,jm,kd).eq.1) then
outar(ir,jm,kd)=label
stack=stack+1
indx(stack)=ir
indy(stack)=jm
indz(stack)=kd
endif
if (outar(il,jb,kd).eq.1) then
outar(il,jb,kd)=label
stack=stack+1
indx(stack)=il
indy(stack)=jb
indz(stack)=kd
endif
if (outar(im,jb,kd).eq.1) then
outar(im,jb,kd)=label
stack=stack+1
indx(stack)=im
indy(stack)=jb
indz(stack)=kd
endif
if (outar(ir,jb,kd).eq.1) then
outar(ir,jb,kd)=label
stack=stack+1
indx(stack)=ir
indy(stack)=jb
indz(stack)=kd
endif
endif
if ((dirh(im,jm,km).eq.1).and.(diru(im,jm,km).eq.1)) then
if (outar(il,jf,ku).eq.1) then
outar(il,jf,ku)=label
stack=stack+1
indx(stack)=il
indy(stack)=jf
indz(stack)=ku
endif
if (outar(im,jf,ku).eq.1) then
outar(im,jf,ku)=label
stack=stack+1
indx(stack)=im
indy(stack)=jf
indz(stack)=ku
endif
if (outar(ir,jf,ku).eq.1) then
outar(ir,jf,ku)=label
stack=stack+1
indx(stack)=ir
indy(stack)=jf
indz(stack)=ku
endif
if (outar(il,jm,ku).eq.1) then
outar(il,jm,ku)=label
stack=stack+1
indx(stack)=il
indy(stack)=jm
indz(stack)=ku
endif
if (outar(ir,jm,ku).eq.1) then
outar(ir,jm,ku)=label
stack=stack+1
indx(stack)=ir
indy(stack)=jm
indz(stack)=ku
endif
if (outar(il,jb,ku).eq.1) then
outar(il,jb,ku)=label
stack=stack+1
indx(stack)=il
indy(stack)=jb
indz(stack)=ku
endif
if (outar(im,jb,ku).eq.1) then
outar(im,jb,ku)=label
stack=stack+1
indx(stack)=im
indy(stack)=jb
indz(stack)=ku
endif
if (outar(ir,jb,ku).eq.1) then
outar(ir,jb,ku)=label
stack=stack+1
indx(stack)=ir
indy(stack)=jb
indz(stack)=ku
endif
endif
! Mark level directly below and above if allowed (dird/u=1)
if (dird(im,jm,km).eq.1) then
if (outar(im,jm,kd).eq.1) then
outar(im,jm,kd)=label
stack=stack+1
indx(stack)=im
indy(stack)=jm
indz(stack)=kd
endif
endif
if (diru(im,jm,km).eq.1) then
if (outar(im,jm,ku).eq.1) then
outar(im,jm,ku)=label
stack=stack+1
indx(stack)=im
indy(stack)=jm
indz(stack)=ku
endif
endif
! Mark other points on same level if allowed (dirh=1)
if (dirh(im,jm,km).eq.1) then
if (outar(il,jf,km).eq.1) then
outar(il,jf,km)=label
stack=stack+1
indx(stack)=il
indy(stack)=jf
indz(stack)=km
endif
if (outar(im,jf,km).eq.1) then
outar(im,jf,km)=label
stack=stack+1
indx(stack)=im
indy(stack)=jf
indz(stack)=km
endif
if (outar(ir,jf,km).eq.1) then
outar(ir,jf,km)=label
stack=stack+1
indx(stack)=ir
indy(stack)=jf
indz(stack)=km
endif
if (outar(il,jm,km).eq.1) then
outar(il,jm,km)=label
stack=stack+1
indx(stack)=il
indy(stack)=jm
indz(stack)=km
endif
if (outar(ir,jm,km).eq.1) then
outar(ir,jm,km)=label
stack=stack+1
indx(stack)=ir
indy(stack)=jm
indz(stack)=km
endif
if (outar(il,jb,km).eq.1) then
outar(il,jb,km)=label
stack=stack+1
indx(stack)=il
indy(stack)=jb
indz(stack)=km
endif
if (outar(im,jb,km).eq.1) then
outar(im,jb,km)=label
stack=stack+1
indx(stack)=im
indy(stack)=jb
indz(stack)=km
endif
if (outar(ir,jb,km).eq.1) then
outar(ir,jb,km)=label
stack=stack+1
indx(stack)=ir
indy(stack)=jb
indz(stack)=km
endif
endif
 
if (stack.gt.0) goto 100
 
! Exit point
700 continue
 
return
 
end subroutine connect
 
! ---------------------------------------------------------
! Calculate the height of the tropopause: BOJAN: although this subroutine is called tropopause for historical reasons (long story), it actually is the core of the 3D-labelling algorithm :-)
! ---------------------------------------------------------
 
SUBROUTINE tropopause (f2,f3,p3,pv3,th3,q3,mdv,xlat,nx,ny,nz,tropo_pv,tropo_th,q_th)
 
! Get the pressure height of the tropopause (f2) and
! cluster the atmosphere into labels 1-5 (f3).
! Given: 3d field: Pressure <p3(nx,ny,nz)>,
! potential vorticity <pv3(nx,ny,nz)>, and potential
! temperature <th3(nx,ny,nz)>. The parameters <nx,ny,nz>
! characterize the grid. The missing data value is <mdv>.
! Bojan July 2013: tropo_pv and tropo_th have to be specified!
 
implicit none
 
! Declaration of subroutine parameters
integer,intent(in) :: nx,ny,nz
real,intent(out) :: f3(nx,ny,nz)
real,intent(out) :: f2(nx,ny)
real,intent(in) :: p3(nx,ny,nz)
real,intent(in) :: q3(nx,ny,nz)
real,intent(inout) :: pv3(nx,ny,nz)
real,intent(in) :: th3(nx,ny,nz)
real,intent(in) :: xlat(ny)
integer :: dirsh(nx,ny,nz),dirsd(nx,ny,nz),dirsu(nx,ny,nz)
integer :: dirth(nx,ny,nz),dirtd(nx,ny,nz),dirtu(nx,ny,nz)
integer :: dirdh(nx,ny,nz),dirdd(nx,ny,nz),dirdu(nx,ny,nz)
integer :: dirfh(nx,ny,nz),dirfd(nx,ny,nz),dirfu(nx,ny,nz)
real :: mdv
 
! Set tropopause thresholds
real :: tropo_pv,tropo_th
! diabatically produced PV (Q-threshold), Feb 2014:
real :: q_th
 
! Set permission for expansion of stratospheric label. A horizontal
! expansion is forbidden if less than <forbid_h> of the grid column
! below a grid point is tropospheric
real :: forbid_h
parameter (forbid_h=0.5)
 
! Internal variables
integer :: i,j,k,l
real :: pv2(nx,ny)
real :: th2(nx,ny)
real :: pvn(nx,ny),pvs(nx,ny)
real :: lat
integer :: st(nx,ny,nz),tr(nx,ny,nz),de(nx,ny,nz),fi(nx,ny,nz)
real :: sign
real :: lev(nx,ny)
integer :: tschk,kabove,kbelow,vertpvchk
real :: below, total
 
! ----- STEP 1a ----- BOJAN: these steps should be the same as described in the appendix of my PhD thesis
! Mark all points above 50 hPa as stratospheric
do i=1,nx
do j=1,ny
do k=1,nz
if (xlat(j).ge.0.) then
sign=1.
else
sign=-1.
endif
if (p3(i,j,k).lt.50.) then
pv3(i,j,k)=sign*1.5*tropo_pv
endif
enddo
enddo
enddo
 
! Mark points with PV > PV threshold (e.g. 2 pvu) OR TH > TH threshold (e.g. 380 K)
do i=1,nx
do j=1,ny
do k=1,nz
if ((abs(pv3(i,j,k)).ge.tropo_pv).or.(th3(i,j,k).ge.tropo_th)) then
st(i,j,k)=1 ! 'stratosphere'
tr(i,j,k)=0 ! 'tropopsphere'
de(i,j,k)=1 ! additional array for de-treatment
else
st(i,j,k)=0
tr(i,j,k)=1
de(i,j,k)=0
endif
enddo
enddo
enddo
 
! ----- STEP 1b -----
 
! Set the expansion permissions for the individual labels (1=allowed, 0=forbidden)
! h=horizontal, v=vertical
do i=1,nx
do j=1,ny
do k=1,nz
 
! Set permissions for stratospheric (st) label: always vertical, never horizontal
dirsh(i,j,k) = 0
dirsd(i,j,k) = 1
dirsu(i,j,k) = 1
 
! Set permissions for fill up thing label (fi): always horizontal, only downward, never upward
dirfh(i,j,k) = 1
dirfd(i,j,k) = 1
dirfu(i,j,k) = 0
 
! Set permissions for tropospheric label (can always propagate)
dirth(i,j,k) = 1
dirtd(i,j,k) = 1
dirtu(i,j,k) = 1
 
! Set permissions for deep fold (de) label: always horizontal, vertical only if air is dry
dirdh(i,j,k) = 1
! Allow vertical propagation of de label only for dry air (Q < Q-threshold)
if ( q3(i,j,k).le.q_th ) then
dirdd(i,j,k) = 1
dirdu(i,j,k) = 1
else
dirdd(i,j,k) = 0
dirdu(i,j,k) = 0
endif
!endif
 
enddo
enddo
enddo
 
! ----- STEP 2a -----
do i=1,nx
do j=1,ny
! Determine stratosphere by connectivity criterion: for all points (x/y) at highest
! model level with st=1, find all connected points with st=1 => they get label st=2
! this label can only propagate vertically and it is important that this happens before
! label st=5 can propagate from the bottom and overwrite st=1 with st=5
if ((st(i,j,nz).eq.1).and.(p3(i,j,nz).lt.500.)) then
call connect (st,2,i,j,nz,dirsh,dirsd,dirsu,nx,ny,nz)
endif
enddo
enddo
 
do i=1,nx
do j=1,ny
 
! ----- STEP 2b -----
 
! for de array, do the same thing. The big difference is that the horizontal propagation is
! always allowed. In the vertical direction, it is only allowed if RH<50%, ensuring that
! very deep tropopause folds, which are not connected by st are reached by de.
 
! This is needed because for some cases, stratospheric air (e.g. at the boundary of a deep-reaching
! streamer) is not identified as such (the label st=2 is not allowed to propagate
! horizontally). The de=2 label can always propagate horizontally and is in the end
! used to distinguish 'real' cutoffs (with st=1 and de=1) from
! deep reaching stratospheric air (with st=1 and de=2)
if ((de(i,j,nz).eq.1).and.(p3(i,j,nz).lt.500.)) then
call connect (de,2,i,j,nz,dirdh,dirdd,dirdu,nx,ny,nz)
endif
 
! ----- STEP 2c -----
 
! Determine troposphere analogeously: for all points (x/y) at lowest model
! level with tr=1 => all connected points get label tr=2
if ((tr(i,j,1).eq.1).and.(p3(i,j,1).gt.300.)) then
call connect (tr,2,i,j,1,dirth,dirtd,dirtu,nx,ny,nz)
endif
 
! ----- STEP 2d -----
! Determine surface-bound PV blobs (|PV|>2): for all points (x/y) at lowest
! model level with st=1 => all connected points get label st=5
! This label can always propagate! (use dirth,dirtv)
if ((st(i,j,1).eq.1).and.(p3(i,j,1).gt.300.)) then
call connect (st,5,i,j,1,dirth,dirtd,dirtu,nx,ny,nz)
endif
 
enddo
enddo
 
! ----- STEP 2e -----
! Remove tropospheric blobs in surface-bound PV blobs (tr=1 surrounded by st=5)
do i=1,nx
do j=1,ny
do k=1,nz
 
if (tr(i,j,k).eq.1) then
tschk=0
do l=k,nz
if (tschk.eq.0) then
! Check if air above is stratosphere (st=2) or top of atmosphere (nz)
if ( (st(i,j,l).eq.2) .or. (l.eq.nz) ) then
tschk=2
! Check if air above is surface-bound PV blob (st=5)
elseif (st(i,j,l).eq.5) then
tschk=5
endif
endif
enddo
! Mark this blob as tr=5
if (tschk.eq.5) then
call connect (tr,5,i,j,k,dirth,dirtd,dirtu,nx,ny,nz)
endif
endif
 
enddo
enddo
enddo
 
! ----- STEP 3 -----
! --------------------------------------------------------------------------------------------------
! Set the stratospheric and troposhperic labels (commented out version is 1=strat/0=everything else)
! --------------------------------------------------------------------------------------------------
do i=1,nx
do j=1,ny
do k=1,nz
f3(i,j,k)=0.
 
if (th3(i,j,k).gt.tropo_th) then
f3(i,j,k) = 2. ! stratosphere TH (2)
! f3(i,j,k) = 1. ! stratosphere TH (1)
else if (st(i,j,k).eq.2) then
f3(i,j,k) = 2. ! stratosphere PV (2)
! f3(i,j,k) = 1. ! stratosphere PV (1)
else if (tr(i,j,k).eq.2) then
f3(i,j,k) = 1. ! troposphere (1)
! f3(i,j,k) = 0. ! troposphere (0)
else if ((st(i,j,k).eq.1).and.(de(i,j,k).eq.1)) then
f3(i,j,k) = 3. ! strat cutoff (3)
! f3(i,j,k) = 0. ! strat cutoff (0)
else if ((st(i,j,k).eq.1).and.(de(i,j,k).eq.2)) then
f3(i,j,k) = 2. ! deep reaching stratospheric air (2)
! f3(i,j,k) = 1. ! deep reaching stratospheric air (1)
else if (tr(i,j,k).eq.1) then
f3(i,j,k) = 4. ! trop cutoff (4)
! f3(i,j,k) = 0. ! trop cutoff (0)
else if (st(i,j,k).eq.5) then
f3(i,j,k) = 5. ! surface-bound PV blob (5)
! f3(i,j,k) = 0. ! surface-bound PV blob (0)
else if (tr(i,j,k).eq.5) then
f3(i,j,k) = 1. ! tropospheric air within surface-bound PV blob (1)
! f3(i,j,k) = 0. ! tropospheric air within surface-bound PV blob (0)
endif
 
enddo
enddo
enddo
 
! --------------------------------------------------------------------------------------------------
! --------------------------------------------------------------------------------------------------
! ----- STEP 6 -----
do i=1,nx
do j=1,ny
do k=1,nz
fi(i,j,k)=0
if (( f3(i,j,k).eq.2 ) .or. ( f3(i,j,k).eq.3 ).or.( f3(i,j,k).eq.5 )) then
fi(i,j,k)=1 ! mark point as possible horizontal fill-up point
endif
end do
end do
end do
! Horizontally fill up 2/3 and 2/5 patches (only aestetically important ;-))
do i=1,nx
do j=1,ny
do k=1,nz
if ( f3(i,j,k).eq.3 ) then ! label 3
call connect (fi,3,i,j,k,dirfh,dirfd,dirfu,nx,ny,nz)
endif
if ( f3(i,j,k).eq.5 ) then ! label 3
call connect (fi,5,i,j,k,dirfh,dirfd,dirfu,nx,ny,nz)
endif
enddo
enddo
enddo
do i=1,nx
do j=1,ny
do k=1,nz
if ( fi(i,j,k).eq.3 ) then ! label 3 that was horizontally filled up
f3(i,j,k)=3.
do l=1,k-1 ! vertically fill up points below if label is 2
if ( f3(i,j,l).eq.2 ) then
f3(i,j,l)=3.
end if
end do
endif
if ( fi(i,j,k).eq.5 ) then ! label 5 that was horizontally filled up
f3(i,j,k)=5.
endif
enddo
enddo
enddo
 
! CHECK that all points obtained a label!
do i=1,nx
do j=1,ny
do k=1,nz
if ( f3(i,j,k).eq.0 ) then ! no label assigned
print*,'Error: no label assigned'
print*,'i,j,k = ',i,j,k
print*,'p3(i,j,k) = ',p3(i,j,k)
print*,'st(i,j,k) = ',st(i,j,k)
print*,'tr(i,j,k) = ',tr(i,j,k)
print*,'de(i,j,k) = ',de(i,j,k)
stop
endif
enddo
enddo
enddo
 
! --------------------------------------------------------------------------------------------------
! DONE with assigning labels, now interpolate field to tropopause pressure
! --------------------------------------------------------------------------------------------------
! BOJAN: here, we calculate the height of the tropopause (folds are ignored, see PhD thesis)
! Remove stratospheric blobs in troposphere and tropospheric blobs in stratosphere
! Brute force: Change PV values above and below tropo threshold
do i=1,nx
do j=1,ny
do k=1,nz
if (xlat(j).ge.0.) then
sign=1.
else
sign=-1.
endif
 
! Set PV of tropo. cutoff in stratosphere
if ( (st(i,j,k).eq.0).and.(tr(i,j,k).eq.1) ) then
pv3(i,j,k)=sign*1.5*tropo_pv
! Set PV of exotic cases where PV < -2 pvu (PV > +2 pvu) in NH (SH) => would otherwise be detected as TP
else if ( ((sign .gt. 0.).and.(pv3(i,j,k).lt.-tropo_pv).or.(sign .lt. 0.).and.(pv3(i,j,k).gt.tropo_pv)).and.(th3(i,j,k).gt.tropo_pv) ) then
pv3(i,j,k)=sign*1.5*tropo_pv
! Set PV of strat. cutoff in trop. or surface bound PV (only below 380K)
!else if ( (((st(i,j,k).eq.1).and.(de(i,j,k).eq.1)).or.(st(i,j,k).eq.5)).and.(th3(i,j,k).le.tropo_th) ) then
else if ( (f3(i,j,k).eq.3).or.(f3(i,j,k).eq.5) ) then
pv3(i,j,k)=sign*0.5*tropo_pv
endif
enddo
enddo
enddo
 
! Calculate the height of the dynamical tropopause (NH:2pvu, SH:-2pvu)
do i=1,nx
do j=1,ny
lev(i,j)=tropo_pv
enddo
enddo
call thipo_lev(p3,pv3,lev,pvn,nx,ny,nz,mdv,0) ! for NH
do i=1,nx
do j=1,ny
do k=1,nz
pv3(i,j,k)=-pv3(i,j,k) ! change sign of PV
enddo
enddo
enddo
call thipo_lev(p3,pv3,lev,pvs,nx,ny,nz,mdv,0) ! for SH
do i=1,nx
do j=1,ny
do k=1,nz
pv3(i,j,k)=-pv3(i,j,k) ! change sign of PV back
enddo
enddo
enddo
do j=1,ny
do i=1,nx
if (xlat(j).ge.0.) then
pv2(i,j)=pvn(i,j) ! NH
else
pv2(i,j)=pvs(i,j) ! SH
endif
enddo
enddo
 
! Special treatment for PV towers (in whole column, |PV|>2)
do i=1,nx
do j=1,ny
if (pv2(i,j).eq.mdv) then
vertpvchk=0
do k=1,nz
if (abs(pv3(i,j,k)).ge.tropo_pv) vertpvchk=vertpvchk+1
enddo
if(vertpvchk.eq.nz) then
!print*,'PV tower! Tropopause reaches ground: ',p3(i,j,1)
pv2(i,j)=p3(i,j,1)
endif
endif
enddo
enddo
 
! Calculate the height of the 380 K surface
do i=1,nx
do j=1,ny
lev(i,j)=tropo_th
enddo
enddo
call thipo_lev(p3,th3,lev,th2,nx,ny,nz,mdv,0)
 
! Set unreasonable values to mdv
do i=1,nx
do j=1,ny
if ((th2(i,j).lt.40.).or.(th2(i,j).gt.800.)) then ! Out of 'regular' limits => check label
if (th2(i,j).ne.mdv) then
!print*,'Warning: p(TP) (TH) out of limits:',th2(i,j)
!print*,'i,j = ',i,j
!print*,'LABEL profile:',f3(i,j,:)
do k=1,nz
if (p3(i,j,k).ge.th2(i,j)) then
kbelow=k
endif
if (p3(i,j,nz+1-k).le.th2(i,j)) then
kabove=nz+1-k
endif
enddo
if ((f3(i,j,kbelow).eq.1).and.(f3(i,j,kabove).eq.2)) then
!print*,'Point is okay!'
else
!print*,'Point is not okay => set to mdv'
th2(i,j)=mdv
endif
endif
endif
if ((pv2(i,j).lt.40.).or.(pv2(i,j).gt.800.)) then ! Out of 'regular' limits => check label
if (pv2(i,j).ne.mdv) then
!print*,'Warning: p(TP) (PV) out of limits:',pv2(i,j)
!print*,'i,j = ',i,j
!print*,'LABEL profile:',f3(i,j,:)
do k=1,nz
if (p3(i,j,k).ge.pv2(i,j)) then
kbelow=k
endif
if (p3(i,j,nz+1-k).le.pv2(i,j)) then
kabove=nz+1-k
endif
enddo
if ((f3(i,j,kbelow).eq.1).and.(f3(i,j,kabove).eq.2)) then
!print*,'Point is okay!'
elseif ((kabove.eq.1).and.(kbelow.eq.1)) then
!print*,'PV tower! Tropopause reaches ground: ',pv2(i,j)
else
!print*,'Point is not okay => set to mdv'
pv2(i,j)=mdv
endif
endif
 
endif
enddo
enddo
 
! --------------------------------------------------------------------------------------------------
! Set the resulting tropopause height (maximum pressure: 380 K, 2 PVU)
! --------------------------------------------------------------------------------------------------
do i=1,nx
do j=1,ny
if ((th2(i,j).ne.mdv).and.(pv2(i,j).eq.mdv)) then
f2(i,j)=th2(i,j)
else if ((th2(i,j).eq.mdv).and.(pv2(i,j).ne.mdv)) then
f2(i,j)=pv2(i,j)
else if (th2(i,j).gt.pv2(i,j)) then
f2(i,j)=th2(i,j)
else
f2(i,j)=pv2(i,j)
endif
enddo
enddo
! --------------------------------------------------------------------------------------------------
! DONE
! --------------------------------------------------------------------------------------------------
 
end subroutine tropopause
 
SUBROUTINE tropofold (label3,p3,pv3,th3,fold,dp, pmin, pmax, sfold,mfold,dfold,mdv,nx,ny,nz,tropo_pv,tropo_th)
 
! -------------------------------------------------------------------------------
! Part 2) Identify tropopause folds (multiple crossings of the tropopause in one column)
! -------------------------------------------------------------------------------
 
! Given the label field, determine whether a tropopause fold is present.
! The parameters <nx,ny,nx> characterize the grid. The missing data value is <mdv>.
 
implicit none
 
! Declaration of subroutine parameters
integer :: nx,ny,nz
real :: mdv
real :: label3(nx,ny,nz)
real :: p3(nx,ny,nz)
real :: pv3(nx,ny,nz)
real :: th3(nx,ny,nz)
real :: pre0,pre1,pre2,ok1,ok2
real :: tropo_pv,tropo_th
real :: fold(nx,ny)
real :: dp(nx,ny)
real :: pmin(nx,ny)
real :: pmax(nx,ny)
real :: sfold(nx,ny)
real :: mfold(nx,ny)
real :: dfold(nx,ny)
 
! Aux integers
integer :: i,j,k,l,k0,k1,k2,ncross
 
! Adjust fields (integer labels, absolute value of PV)
do i=1,nx
do j=1,ny
do k=1,nz
pv3(i,j,k) = abs(pv3(i,j,k))
enddo
enddo
enddo
 
! Definition of labels
! 1 : troposphere
! 2 : stratosphere
! 3 : stratosperic cutoff
! 4 : tropospheric cutoff
! 5 : surface-bound PV blobs
 
do i=1,nx
do j=1,ny
k0 = 0
k1 = 0
k2 = 0
ncross = 0
 
do k=nz-1,1,-1
 
! strat-trop transition
if ( (label3(i,j,k+1).eq.2).and.(label3(i,j,k).eq.1) ) then
 
ncross = ncross + 1
if ( k2.eq.0 ) then
k2 = k
elseif ( k1.ne.0 ) then
k0 = k
endif
endif
! special case for folds that are affected by the q-criterion (q<0.0001 only)
! example (top to bottom): 222211122231111 need to recognize the 2-3-1 transition
! transitions like 222222223111 should not be recognized, thus only if k2.ne.0
if ( (label3(i,j,k+1).eq.2).and.(label3(i,j,k).eq.3) ) then
 
ncross = ncross + 1
if (( k2.ne.0 ) .and. ( k1.ne.0 )) then
k0 = k
endif
endif
 
! trop-strat transition
if ( (label3(i,j,k+1).eq.1).and.(label3(i,j,k).eq.2) ) then
 
ncross = ncross + 1
if ( ( k2.ne.0 ).and.( k1.eq.0 ) ) then
k1 = k
endif
endif
enddo
 
! Skip further steps if we have a single tropopause
if ( ncross.le.2 ) goto 100
! Get the exact pressures at the crossings
pre0 = 0.
pre1 = 0.
pre2 = 0.
! lowest point (0) => pre0
ok1 = ( pv3(i,j,k0+1)-tropo_pv ) * ( pv3(i,j,k0)-tropo_pv )
ok2 = ( th3(i,j,k0+1)-tropo_th ) * ( th3(i,j,k0)-tropo_th )
 
if ( ok1.le.0. ) then
pre0 = p3(i,j,k0) + ( p3(i,j,k0+1) - p3(i,j,k0) ) * &
( tropo_pv - pv3(i,j,k0) )/( pv3(i,j,k0+1) - pv3(i,j,k0) )
elseif ( ok2.le.0. ) then
pre0 = p3(i,j,k0) + ( p3(i,j,k0+1) - p3(i,j,k0) ) * &
( tropo_th - th3(i,j,k0) )/( th3(i,j,k0+1) - th3(i,j,k0) )
endif
! middle point (1) => pre1
ok1 = ( pv3(i,j,k1+1)-tropo_pv ) * ( pv3(i,j,k1)-tropo_pv )
ok2 = ( th3(i,j,k1+1)-tropo_th ) * ( th3(i,j,k1)-tropo_th )
if ( ok1.le.0. ) then
pre1 = p3(i,j,k1) + ( p3(i,j,k1+1) - p3(i,j,k1) ) * &
( tropo_pv - pv3(i,j,k1) )/( pv3(i,j,k1+1) - pv3(i,j,k1) )
elseif ( ok2.le.0. ) then
pre1 = p3(i,j,k1) + ( p3(i,j,k1+1) - p3(i,j,k1) ) * &
( tropo_th - th3(i,j,k1) )/( th3(i,j,k1+1) - th3(i,j,k1) )
endif
! upper point (2) => pre2
ok1 = ( pv3(i,j,k2+1)-tropo_pv ) * ( pv3(i,j,k2)-tropo_pv )
ok2 = ( th3(i,j,k2+1)-tropo_th ) * ( th3(i,j,k2)-tropo_th )
if ( ok1.le.0. ) then
pre2 = p3(i,j,k2) + ( p3(i,j,k2+1) - p3(i,j,k2) ) * &
( tropo_pv - pv3(i,j,k2) )/( pv3(i,j,k2+1) - pv3(i,j,k2) )
elseif ( ok2.le.0. ) then
pre2 = p3(i,j,k2) + ( p3(i,j,k2+1) - p3(i,j,k2) ) * &
( tropo_th - th3(i,j,k2) )/( th3(i,j,k2+1) - th3(i,j,k2) )
endif
 
! Decide whether all pressure values are ok
if ( ( pre0.lt.p3(i,j,k0+1) ).or. &
( pre0.gt.p3(i,j,k0 ) ).or. &
( pre1.lt.p3(i,j,k1+1) ).or. &
( pre1.gt.p3(i,j,k1 ) ).or. &
( pre2.lt.p3(i,j,k2+1) ).or. &
( pre2.gt.p3(i,j,k2 ) ) ) goto 100
! Everything is fine: remember the fold
fold(i,j) = 1.
dp(i,j) = pre1 - pre2
pmin(i,j) = pre2
pmax(i,j) = pre0
 
! Exit point for loop
100 continue
 
enddo ! y
enddo ! x
 
where((fold .gt. 0.0) .and. (dp .ge. 50.0) .and. (dp .lt. 200.0))
sfold=1.0
elsewhere((fold .gt. 0.0) .and. (dp .ge. 200.0) .and. (dp .lt. 350.0))
mfold=1.0
elsewhere((fold .gt. 0.0) .and. (dp .ge. 350.0))
dfold=1.0
end where
end subroutine tropofold
 
! -----------------------------------------------------------
! Interpolation onto theta surface (top -> down): BOJAN: not sure you have all the interpolation subroutines used here, contact us or feel free to use any other interpolation routine that does the job :-)
! -----------------------------------------------------------
subroutine thipo_lev(var3d,th3d,lev,var,nx,ny,nz,mdv,mode)
 
! Interpolates the 3d variable var3d on the isentropic surface
! defined by lev. The interpolated field is returned as var.
! th3d denotes the 3d theta array.
! mode determines the way of vertical interpolation:
! mode=0 is for linear interpolation
! mode=1 is for logarithmic interpolation
 
integer :: nx,ny,nz,mode
real :: lev(nx,ny),mdv
real :: var3d(nx,ny,nz),th3d(nx,ny,nz),var(nx,ny)
 
integer :: i,j,k
real :: kind
real :: int3dm
 
do i=1,nx
do j=1,ny
 
if (lev(i,j).ne.mdv) then
 
kind=0.
do k=nz-1,1,-1
if ((th3d(i,j,k).le.lev(i,j)).and.(th3d(i,j,k+1).ge.lev(i,j))) then
kind=float(k)+(th3d(i,j,k)-lev(i,j))/(th3d(i,j,k)-th3d(i,j,k+1))
goto 100
endif
enddo
100 continue
 
if (kind.eq.0) then
var(i,j)=mdv
else
var(i,j)=int3dm(var3d,nx,ny,nz,float(i),float(j),kind,mdv)
endif
 
else
 
var(i,j)=mdv
 
endif
 
enddo
enddo
 
return
end subroutine thipo_lev
 
! -----------------------------------------------------------
! Interpolation onto theta surface (bottom -> up)
! -----------------------------------------------------------
 
subroutine pipo_lev(var3d,p3d,lev,var,nx,ny,nz,mdv,mode)
 
! Interpolates the 3d variable var3d on the pressure surface
! defined by lev. The interpolated field is returned as var.
! p3d denotes the 3d pressure array.
! mode determines the way of vertical interpolation:
! mode=0 is for linear interpolation
! mode=1 is for logarithmic interpolation
 
integer :: nx,ny,nz,mode
real :: lev(nx,ny),mdv
real :: var3d(nx,ny,nz),p3d(nx,ny,nz),var(nx,ny)
 
integer :: i,j,k
real :: kind
real :: int3dm
 
do i=1,nx
do j=1,ny
 
if (lev(i,j).ne.mdv) then
 
kind=0.
do k=1,nz-1
if ((p3d(i,j,k).ge.lev(i,j)).and.(p3d(i,j,k+1).le.lev(i,j))) then
kind=float(k)+(p3d(i,j,k)-lev(i,j))/(p3d(i,j,k)-p3d(i,j,k+1))
goto 100
endif
enddo
100 continue
 
if (kind.eq.0.) then
var(i,j)=mdv
else
var(i,j)=int3dm(var3d,nx,ny,nz,float(i),float(j),kind,mdv)
endif
 
else
 
var(i,j)=mdv
 
endif
 
enddo
enddo
 
return
end subroutine pipo_lev
 
 
! ------------------------------------------------------------------
! Auxiliary routines
! ------------------------------------------------------------------
 
/tags/1.0/3d_labelling_and_fold_id_20150126.zip
Cannot display: file marked as a binary type.
svn:mime-type = application/zip
Property changes:
Added: svn:mime-type
+application/zip
\ No newline at end of property
/tags/1.0/EMAC.nml
0,0 → 1,15
! -*- f90 -*-
&CTRL
file_input = '/ptmp/andrep/ZANIS/ZANIS__________20000801_0000_fold.nc'
file_output = 'test.nc'
X_name = 'lon'
Y_name = 'lat'
Z_name = 'lev'
T_name = 'time'
APS_name = 'aps'
HYAM_name = 'hyam'
HYBM_name = 'hybm'
Q_name = 'qm1'
PV_name = 'PV'
PT_name = 'tpot'
/
/tags/1.0/Makefile
0,0 → 1,26
# -*- Makefile -*-
# contact: Andrea Pozzer MPIC (andrea.pozzer@mpic.de)
# ----------------------------------------------
export
SHELL = sh
 
# hydra intel fortran
F90 = ifort
 
#debugging
#F90FLAGS = -autodouble -cpp -g -debug full -traceback -O3
#production
F90FLAGS = -autodouble -cpp -O3
 
FC = $(F90)
FFLAGS = $(F90FLAGS)
 
NETCDF_ROOT = /ptmp/mpcdata/software/x86_64-suse-linux/netcdf/v3.6.3_i
 
LIBS = -L$(NETCDF_ROOT)/lib -lnetcdf
INCLUDES = -I$(NETCDF_ROOT)/include
 
INSTALLDIR = .
# ----------------------------------------------
include main.mk
# ----------------------------------------------
/tags/1.0/README.txt
0,0 → 1,54
Tropopause folding calculator:
(*original reference to be added*)
 
original work by B.Skerlak (bojan.skerlak@env.ethz.ch)
extended for netcdf data by A.Pozzer (andrea.pozzer@mpic.de)
 
------------------------------------------------------------------
To use the algorithm:
 
1) Extract the zip file
2) Modify the Makefile to the current machine set-up
3) Compile to code with the command "make".
Other options possible with "make help".
4) Create a namelist file (filename.nml) with the specific of the
netcdf files to be analyzed (see EMAC.nml for example).
5) run the code with the command "3d_labelling_and_fold_id.exe filename.nml"
 
-------------------------------------------------------------------
NML format:
the namelist file should contain the following fields:
------------------------------------------------------------------
&CTRL
! input file (path included)
file_input = '/ptmp/andrep/ZANIS/ZANIS__________20000801_0000_fold.nc'
! output file (path included)
file_output = 'test.nc'
! name of the longitude dimension in the file
X_name = 'lon'
! name of the latitude dimension in the file
Y_name = 'lat'
! name of the level dimension in the file
Z_name = 'lev'
! name of the time dimension in the file
T_name = 'time'
! name of the surface pressure
APS_name = 'aps'
! name of the hybrid coefficient A
HYAM_name = 'hyam'
! name of the hybrid coefficient B
HYBM_name = 'hybm'
! name of the specific humidity field
Q_name = 'qm1'
! name of the potential vorticity field
PV_name = 'PV'
! name of the potential temperature field
PT_name = 'tpot'
/
------------------------------------------------------------------
WARNINGS:
1) The code works automatically only if all the needed data (i.e. Q,PV,PT)
are included in the same input file.
2) Only hybrid pressure coordinates included in the pressure
calculation.
------------------------------------------------------------------
/tags/1.0/libipo.f
0,0 → 1,1368
real function int2d(ar,n1,n2,rid,rjd)
c-----------------------------------------------------------------------
c Purpose:
c This subroutine interpolates a 2d-array to an arbitrary
c location within the grid.
c Arguments:
c ar real input surface pressure, define as ar(n1,n2)
c n1,n2 int input dimensions of ar
c ri,rj real input grid location to be interpolated to
c History:
c-----------------------------------------------------------------------
 
c argument declarations
integer n1,n2
real ar(n1,n2), rid,rjd
 
c local declarations
integer i,j,ip1,jp1,ih,jh
real frac0i,frac0j,frac1i,frac1j,ri,rj
 
c do linear interpolation
ri=amax1(1.,amin1(float(n1),rid))
rj=amax1(1.,amin1(float(n2),rjd))
ih=nint(ri)
jh=nint(rj)
 
c Check for interpolation in i
* if (abs(float(ih)-ri).lt.1.e-3) then
* i =ih
* ip1=ih
* else
i =min0(int(ri),n1-1)
ip1=i+1
* endif
 
c Check for interpolation in j
* if (abs(float(jh)-rj).lt.1.e-3) then
* j =jh
* jp1=jh
* else
j =min0(int(rj),n2-1)
jp1=j+1
* endif
 
if ((i.eq.ip1).and.(j.eq.jp1)) then
c no interpolation at all
int2d=ar(i,j)
else
frac0i=ri-float(i)
frac0j=rj-float(j)
frac1i=1.-frac0i
frac1j=1.-frac0j
int2d = ar(i ,j ) * frac1i * frac1j
& + ar(i ,jp1) * frac1i * frac0j
& + ar(ip1,j ) * frac0i * frac1j
& + ar(ip1,jp1) * frac0i * frac0j
endif
end
real function int2dm(ar,n1,n2,rid,rjd,misdat)
c-----------------------------------------------------------------------
c Purpose:
c This subroutine interpolates a 2d-array to an arbitrary
c location within the grid. The interpolation includes the
c testing of the missing data flag 'misdat'.
c Arguments:
c ar real input surface pressure, define as ar(n1,n2)
c n1,n2 int input dimensions of ar
c ri,rj real input grid location to be interpolated to
c misdat real input missing data flag (on if misdat<>0)
c Warning:
c This routine has not yet been seriously tested
c History:
c-----------------------------------------------------------------------
c argument declarations
integer n1,n2
real ar(n1,n2), rid,rjd, misdat
c local declarations
integer i,j,ip1,jp1,ih,jh
real frac0i,frac0j,frac1i,frac1j,ri,rj,int2d
c check if routine without missing data checking can be called instead
if (misdat.eq.0.) then
int2dm=int2d(ar,n1,n2,rid,rjd)
return
endif
c do linear interpolation
ri=amax1(1.,amin1(float(n1),rid))
rj=amax1(1.,amin1(float(n2),rjd))
ih=nint(ri)
jh=nint(rj)
c Check for interpolation in i
* if (abs(float(ih)-ri).lt.1.e-3) then
* i =ih
* ip1=ih
* else
i =min0(int(ri),n1-1)
ip1=i+1
* endif
c Check for interpolation in j
* if (abs(float(jh)-rj).lt.1.e-3) then
* j =jh
* jp1=jh
* else
j =min0(int(rj),n2-1)
jp1=j+1
* endif
if ((i.eq.ip1).and.(j.eq.jp1)) then
c no interpolation at all
int2dm=ar(i,j)
else
if ((misdat.eq.ar(i ,j )).or.
& (misdat.eq.ar(i ,jp1)).or.
& (misdat.eq.ar(ip1,j )).or.
& (misdat.eq.ar(ip1,jp1))) then
int2dm=misdat
else
frac0i=ri-float(i)
frac0j=rj-float(j)
frac1i=1.-frac0i
frac1j=1.-frac0j
int2dm = ar(i ,j ) * frac1i * frac1j
& + ar(i ,jp1) * frac1i * frac0j
& + ar(ip1,j ) * frac0i * frac1j
& + ar(ip1,jp1) * frac0i * frac0j
endif
endif
end
real function int2dp(ar,n1,n2,rid,rjd)
c-----------------------------------------------------------------------
c Purpose:
c This subroutine interpolates a 2d-array to an arbitrary
c location within the grid. The 2d-array is periodic in the
c i-direction: ar(0,.)=ar(n1,.) and ar(1,.)=ar(n1+1,.).
c Therefore rid can take values in the range 0,...,n1+1
c Arguments:
c ar real input surface pressure, define as ar(n1,n2)
c n1,n2 int input dimensions of ar
c ri,rj real input grid location to be interpolated to
c History:
c-----------------------------------------------------------------------
c argument declarations
integer n1,n2
real ar(0:n1+1,n2), rid,rjd
c local declarations
integer i,j,ip1,jp1,ih,jh
real frac0i,frac0j,frac1i,frac1j,ri,rj
c do linear interpolation
ri=amax1(0.,amin1(float(n1+1),rid))
rj=amax1(1.,amin1(float(n2),rjd))
ih=nint(ri)
jh=nint(rj)
c Check for interpolation in i
* if (abs(float(ih)-ri).lt.1.e-5) then
* i =ih
* ip1=ih
* else
i =min0(int(ri),n1-1)
ip1=i+1
* endif
c Check for interpolation in j
* if (abs(float(jh)-rj).lt.1.e-5) then
* j =jh
* jp1=jh
* else
j =min0(int(rj),n2-1)
jp1=j+1
* endif
if ((i.eq.ip1).and.(j.eq.jp1)) then
c no interpolation at all
int2dp=ar(i,j)
else
frac0i=ri-float(i)
frac0j=rj-float(j)
frac1i=1.-frac0i
frac1j=1.-frac0j
int2dp = ar(i ,j ) * frac1i * frac1j
& + ar(i ,jp1) * frac1i * frac0j
& + ar(ip1,j ) * frac0i * frac1j
& + ar(ip1,jp1) * frac0i * frac0j
endif
end
real function cint2d(ar,n1,n2,n3,rid,rjd,ikd)
c-----------------------------------------------------------------------
c Purpose:
c This subroutine interpolates a 3d-array to an arbitrary
c location in the horizontal. ikd specifies the level (must
c be an integer). A bicubic method is applied (following
c the Numerical Recipes).
c Arguments:
c ar real input field, define as ar(n1,n2,n3)
c n1,n2 int input dimensions of ar
c rid,rjd real input grid location to be interpolated to
c ikd int input level for interpolation
c History:
c-----------------------------------------------------------------------
 
c argument declarations
integer n1,n2,n3,ikd
real ar(n1,n2,n3), rid,rjd
 
c local declarations
integer i,j,k
real y(4),y1(4),y2(4),y12(4)
 
c indices of lower left corner of interpolation grid
 
i=int(rid)
j=int(rjd)
k=ikd
 
c do not interpolate if i or j are at the boundaries
 
if ((i.eq.1).or.(j.eq.1).or.(i.eq.n1).or.(j.eq.n2)) then
cint2d=ar(i,j,k)
return
endif
 
c define the arrays y,y1,y2,y12 necessary for the bicubic
c interpolation (following the Numerical Recipes).
 
y(1)=ar(i,j,k)
y(2)=ar(i+1,j,k)
y(3)=ar(i+1,j+1,k)
y(4)=ar(i,j+1,k)
 
y1(1)=-(ar(i-1,j,k)-ar(i+1,j,k))/2.
y1(2)=-(ar(i,j,k)-ar(i+2,j,k))/2.
y1(3)=-(ar(i,j+1,k)-ar(i+2,j+1,k))/2.
y1(4)=-(ar(i-1,j+1,k)-ar(i+1,j+1,k))/2.
 
y2(1)=-(ar(i,j-1,k)-ar(i,j+1,k))/2.
y2(2)=-(ar(i+1,j-1,k)-ar(i+1,j+1,k))/2.
y2(3)=-(ar(i+1,j,k)-ar(i+1,j+2,k))/2.
y2(4)=-(ar(i,j,k)-ar(i,j+2,k))/2.
 
y12(1)=(ar(i+1,j+1,k)-ar(i+1,j-1,k)-ar(i-1,j+1,k)+ar(i-1,j-1,k))/4.
y12(2)=(ar(i+2,j+1,k)-ar(i+2,j-1,k)-ar(i,j+1,k)+ar(i,j-1,k))/4.
y12(3)=(ar(i+2,j+2,k)-ar(i+2,j,k)-ar(i,j+2,k)+ar(i,j,k))/4.
y12(4)=(ar(i+1,j+2,k)-ar(i+1,j,k)-ar(i-1,j+2,k)+ar(i-1,j,k))/4.
 
call bcuint(y,y1,y2,y12,i,j,rid,rjd,cint2d)
return
end
real function int3d(ar,n1,n2,n3,rid,rjd,rkd)
c-----------------------------------------------------------------------
c Purpose:
c This subroutine interpolates a 3d-array to an arbitrary
c location within the grid.
c Arguments:
c ar real input surface pressure, define as ar(n1,n2,n3)
c n1,n2,n3 int input dimensions of ar
c ri,rj,rk real input grid location to be interpolated to
c History:
c-----------------------------------------------------------------------
 
c argument declarations
integer n1,n2,n3
real ar(n1,n2,n3), rid,rjd,rkd
 
c local declarations
integer i,j,k,ip1,jp1,kp1,ih,jh,kh
real frac0i,frac0j,frac0k,frac1i,frac1j,frac1k,ri,rj,rk
 
c do linear interpolation
ri=amax1(1.,amin1(float(n1),rid))
rj=amax1(1.,amin1(float(n2),rjd))
rk=amax1(1.,amin1(float(n3),rkd))
ih=nint(ri)
jh=nint(rj)
kh=nint(rk)
 
c Check for interpolation in i
* if (abs(float(ih)-ri).lt.1.e-3) then
* i =ih
* ip1=ih
* else
i =min0(int(ri),n1-1)
ip1=i+1
* endif
 
c Check for interpolation in j
if (abs(float(jh)-rj).lt.1.e-3) then
j =jh
jp1=jh
else
j =min0(int(rj),n2-1)
jp1=j+1
endif
 
c Check for interpolation in k
* if (abs(float(kh)-rk).lt.1.e-3) then
* k =kh
* kp1=kh
* else
k =min0(int(rk),n3-1)
kp1=k+1
* endif
 
if (k.eq.kp1) then
c no interpolation in k
if ((i.eq.ip1).and.(j.eq.jp1)) then
c no interpolation at all
int3d=ar(i,j,k)
c print *,'int3d 00: ',rid,rjd,rkd,int3d
else
c horizontal interpolation only
frac0i=ri-float(i)
frac0j=rj-float(j)
frac1i=1.-frac0i
frac1j=1.-frac0j
int3d = ar(i ,j ,k ) * frac1i * frac1j
& + ar(i ,jp1,k ) * frac1i * frac0j
& + ar(ip1,j ,k ) * frac0i * frac1j
& + ar(ip1,jp1,k ) * frac0i * frac0j
c print *,'int3d 10: ',rid,rjd,rkd,int3d
endif
else
frac0k=rk-float(k)
frac1k=1.-frac0k
if ((i.eq.ip1).and.(j.eq.jp1)) then
c vertical interpolation only
int3d = ar(i ,j ,k ) * frac1k
& + ar(i ,j ,kp1) * frac0k
c print *,'int3d 01: ',rid,rjd,rkd,int3d
else
c full 3d interpolation
frac0i=ri-float(i)
frac0j=rj-float(j)
frac1i=1.-frac0i
frac1j=1.-frac0j
int3d = ar(i ,j ,k ) * frac1i * frac1j * frac1k
& + ar(i ,jp1,k ) * frac1i * frac0j * frac1k
& + ar(ip1,j ,k ) * frac0i * frac1j * frac1k
& + ar(ip1,jp1,k ) * frac0i * frac0j * frac1k
& + ar(i ,j ,kp1) * frac1i * frac1j * frac0k
& + ar(i ,jp1,kp1) * frac1i * frac0j * frac0k
& + ar(ip1,j ,kp1) * frac0i * frac1j * frac0k
& + ar(ip1,jp1,kp1) * frac0i * frac0j * frac0k
c print *,'int3d 11: ',rid,rjd,rkd,int3d
endif
endif
end
real function int3dm(ar,n1,n2,n3,rid,rjd,rkd,misdat)
c-----------------------------------------------------------------------
c Purpose:
c This subroutine interpolates a 3d-array to an arbitrary
c location within the grid. The interpolation includes the
c testing of the missing data flag 'misdat'.
c Arguments:
c ar real input surface pressure, define as ar(n1,n2,n3)
c n1,n2,n3 int input dimensions of ar
c ri,rj,rk real input grid location to be interpolated to
c misdat real input missing data flag (on if misdat<>0)
c Warning:
c This routine has not yet been seriously tested
c History:
c-----------------------------------------------------------------------
 
c argument declarations
integer n1,n2,n3
real ar(n1,n2,n3), rid,rjd,rkd, misdat
 
c local declarations
integer i,j,k,ip1,jp1,kp1,ih,jh,kh
real frac0i,frac0j,frac0k,frac1i,frac1j,frac1k,ri,rj,rk,int3d
 
c check if routine without missing data checking can be called instead
if (misdat.eq.0.) then
int3dm=int3d(ar,n1,n2,n3,rid,rjd,rkd)
return
endif
 
c do linear interpolation
ri=amax1(1.,amin1(float(n1),rid))
rj=amax1(1.,amin1(float(n2),rjd))
rk=amax1(1.,amin1(float(n3),rkd))
ih=nint(ri)
jh=nint(rj)
kh=nint(rk)
 
c Check for interpolation in i
* if (abs(float(ih)-ri).lt.1.e-3) then
* i =ih
* ip1=ih
* else
i =min0(int(ri),n1-1)
ip1=i+1
* endif
 
c Check for interpolation in j
* if (abs(float(jh)-rj).lt.1.e-3) then
* j =jh
* jp1=jh
* else
j =min0(int(rj),n2-1)
jp1=j+1
* endif
 
c Check for interpolation in k
* if (abs(float(kh)-rk).lt.1.e-3) then
* k =kh
* kp1=kh
* else
k =min0(int(rk),n3-1)
kp1=k+1
* endif
 
if (k.eq.kp1) then
c no interpolation in k
if ((i.eq.ip1).and.(j.eq.jp1)) then
c no interpolation at all
if (misdat.eq.ar(i,j,k)) then
int3dm=misdat
else
int3dm=ar(i,j,k)
endif
c print *,'int3dm 00: ',rid,rjd,rkd,int3dm
else
c horizontal interpolation only
if ((misdat.eq.ar(i ,j ,k )).or.
& (misdat.eq.ar(i ,jp1,k )).or.
& (misdat.eq.ar(ip1,j ,k )).or.
& (misdat.eq.ar(ip1,jp1,k ))) then
int3dm=misdat
else
frac0i=ri-float(i)
frac0j=rj-float(j)
frac1i=1.-frac0i
frac1j=1.-frac0j
int3dm = ar(i ,j ,k ) * frac1i * frac1j
& + ar(i ,jp1,k ) * frac1i * frac0j
& + ar(ip1,j ,k ) * frac0i * frac1j
& + ar(ip1,jp1,k ) * frac0i * frac0j
c print *,'int3dm 10: ',rid,rjd,rkd,int3dm
endif
endif
else
frac0k=rk-float(k)
frac1k=1.-frac0k
if ((i.eq.ip1).and.(j.eq.jp1)) then
c vertical interpolation only
if ((misdat.eq.ar(i ,j ,k )).or.
& (misdat.eq.ar(i ,j ,kp1))) then
int3dm=misdat
else
int3dm = ar(i ,j ,k ) * frac1k
& + ar(i ,j ,kp1) * frac0k
c print *,'int3dm 01: ',rid,rjd,rkd,int3dm
endif
else
c full 3d interpolation
if ((misdat.eq.ar(i ,j ,k )).or.
& (misdat.eq.ar(i ,jp1,k )).or.
& (misdat.eq.ar(ip1,j ,k )).or.
& (misdat.eq.ar(ip1,jp1,k )).or.
& (misdat.eq.ar(i ,j ,kp1)).or.
& (misdat.eq.ar(i ,jp1,kp1)).or.
& (misdat.eq.ar(ip1,j ,kp1)).or.
& (misdat.eq.ar(ip1,jp1,kp1))) then
int3dm=misdat
else
frac0i=ri-float(i)
frac0j=rj-float(j)
frac1i=1.-frac0i
frac1j=1.-frac0j
int3dm = ar(i ,j ,k ) * frac1i * frac1j * frac1k
& + ar(i ,jp1,k ) * frac1i * frac0j * frac1k
& + ar(ip1,j ,k ) * frac0i * frac1j * frac1k
& + ar(ip1,jp1,k ) * frac0i * frac0j * frac1k
& + ar(i ,j ,kp1) * frac1i * frac1j * frac0k
& + ar(i ,jp1,kp1) * frac1i * frac0j * frac0k
& + ar(ip1,j ,kp1) * frac0i * frac1j * frac0k
& + ar(ip1,jp1,kp1) * frac0i * frac0j * frac0k
c print *,'int3dm 11: ',rid,rjd,rkd,int3dm
endif
endif
endif
end
real function int3dmlog(ar,n1,n2,n3,rid,rjd,rkd,misdat)
c-----------------------------------------------------------------------
c Purpose:
c This subroutine interpolates a 3d-array to an arbitrary
c location within the grid. The interpolation includes the
c testing of the missing data flag 'misdat'.
c Prior to vertical interpolations the log is taken from the array.
c Arguments:
c ar real input surface pressure, define as ar(n1,n2,n3)
c n1,n2,n3 int input dimensions of ar
c ri,rj,rk real input grid location to be interpolated to
c misdat real input missing data flag (on if misdat<>0)
c Warning:
c This routine has not yet been seriously tested
c History:
c-----------------------------------------------------------------------
 
c argument declarations
integer n1,n2,n3
real ar(n1,n2,n3), rid,rjd,rkd, misdat
 
c local declarations
integer i,j,k,ip1,jp1,kp1,ih,jh,kh
real frac0i,frac0j,frac0k,frac1i,frac1j,frac1k,ri,rj,rk,int3d
 
c print*,'hallo in SR int3dmlog'
c check if routine without missing data checking can be called instead
if (misdat.eq.0.) then
int3dmlog=int3d(ar,n1,n2,n3,rid,rjd,rkd)
return
endif
 
c do linear interpolation
ri=amax1(1.,amin1(float(n1),rid))
rj=amax1(1.,amin1(float(n2),rjd))
rk=amax1(1.,amin1(float(n3),rkd))
ih=nint(ri)
jh=nint(rj)
kh=nint(rk)
 
c Check for interpolation in i
* if (abs(float(ih)-ri).lt.1.e-3) then
* i =ih
* ip1=ih
* else
i =min0(int(ri),n1-1)
ip1=i+1
* endif
 
c Check for interpolation in j
* if (abs(float(jh)-rj).lt.1.e-3) then
* j =jh
* jp1=jh
* else
j =min0(int(rj),n2-1)
jp1=j+1
* endif
 
c Check for interpolation in k
* if (abs(float(kh)-rk).lt.1.e-3) then
* k =kh
* kp1=kh
* else
k =min0(int(rk),n3-1)
kp1=k+1
* endif
 
if (k.eq.kp1) then
c no interpolation in k
if ((i.eq.ip1).and.(j.eq.jp1)) then
c no interpolation at all
if (misdat.eq.ar(i,j,k)) then
int3dmlog=misdat
else
int3dmlog=ar(i,j,k)
endif
c print *,'int3dmlog 00: ',rid,rjd,rkd,int3dmlog
else
c horizontal interpolation only
if ((misdat.eq.ar(i ,j ,k )).or.
& (misdat.eq.ar(i ,jp1,k )).or.
& (misdat.eq.ar(ip1,j ,k )).or.
& (misdat.eq.ar(ip1,jp1,k ))) then
int3dmlog=misdat
else
frac0i=ri-float(i)
frac0j=rj-float(j)
frac1i=1.-frac0i
frac1j=1.-frac0j
int3dmlog = ar(i ,j ,k ) * frac1i * frac1j
& + ar(i ,jp1,k ) * frac1i * frac0j
& + ar(ip1,j ,k ) * frac0i * frac1j
& + ar(ip1,jp1,k ) * frac0i * frac0j
c print *,'int3dmlog 10: ',rid,rjd,rkd,int3dmlog
endif
endif
else
frac0k=rk-float(k)
frac1k=1.-frac0k
if ((i.eq.ip1).and.(j.eq.jp1)) then
c vertical interpolation only
if ((misdat.eq.ar(i ,j ,k )).or.
& (misdat.eq.ar(i ,j ,kp1))) then
int3dmlog=misdat
else
int3dmlog = log(ar(i ,j ,k )) * frac1k
& + log(ar(i ,j ,kp1)) * frac0k
int3dmlog = exp(int3dmlog)
c print *,'int3dmlog 01: ',rid,rjd,rkd,int3dmlog
endif
else
c full 3d interpolation
if ((misdat.eq.ar(i ,j ,k )).or.
& (misdat.eq.ar(i ,jp1,k )).or.
& (misdat.eq.ar(ip1,j ,k )).or.
& (misdat.eq.ar(ip1,jp1,k )).or.
& (misdat.eq.ar(i ,j ,kp1)).or.
& (misdat.eq.ar(i ,jp1,kp1)).or.
& (misdat.eq.ar(ip1,j ,kp1)).or.
& (misdat.eq.ar(ip1,jp1,kp1))) then
int3dmlog=misdat
else
frac0i=ri-float(i)
frac0j=rj-float(j)
frac1i=1.-frac0i
frac1j=1.-frac0j
int3dmlog = log(ar(i ,j ,k )) * frac1i * frac1j * frac1k
& + log(ar(i ,jp1,k )) * frac1i * frac0j * frac1k
& + log(ar(ip1,j ,k )) * frac0i * frac1j * frac1k
& + log(ar(ip1,jp1,k )) * frac0i * frac0j * frac1k
& + log(ar(i ,j ,kp1)) * frac1i * frac1j * frac0k
& + log(ar(i ,jp1,kp1)) * frac1i * frac0j * frac0k
& + log(ar(ip1,j ,kp1)) * frac0i * frac1j * frac0k
& + log(ar(ip1,jp1,kp1)) * frac0i * frac0j * frac0k
int3dmlog = exp(int3dmlog)
c print *,'int3dmlog 11: ',rid,rjd,rkd,int3dmlog
endif
endif
endif
end
 
real function int3dmvc(ar,n1,n2,n3,rid,rjd,rkd,misdat)
c-----------------------------------------------------------------------
c Purpose:
c This subroutine interpolates a 3d-array to an arbitrary
c location within the grid. The interpolation includes the
c testing of the missing data flag 'misdat'.
c In the vertical a Lagrangian cubic interpolation is
c performed.
c Arguments:
c ar real input surface pressure, define as ar(n1,n2,n3)
c n1,n2,n3 int input dimensions of ar
c ri,rj,rk real input grid location to be interpolated to
c misdat real input missing data flag (on if misdat<>0)
c Warning:
c This routine has not yet been seriously tested
c History:
c-----------------------------------------------------------------------
 
c argument declarations
integer n1,n2,n3
real ar(n1,n2,n3), rid,rjd,rkd, misdat
c local declarations
integer i,j,k,ip1,jp1,kp1,ih,jh,kh,klow,n
real frac0i,frac0j,frac1i,frac1j,ri,rj,rk
real int2d(4)
 
real int3dm
 
c if n3 < 4 then do call linear interpolation in the vertical
if (n3.lt.4) then
int3dmvc=int3dm(ar,n1,n2,n3,rid,rjd,rkd,misdat)
return
endif
c do linear interpolation in the horizontal
ri=amax1(1.,amin1(float(n1),rid))
rj=amax1(1.,amin1(float(n2),rjd))
rk=amax1(1.,amin1(float(n3),rkd))
ih=nint(ri)
jh=nint(rj)
kh=nint(rk)
c Check for interpolation in i
* if (abs(float(ih)-ri).lt.1.e-3) then
* i =ih
* ip1=ih
* else
i =min0(int(ri),n1-1)
ip1=i+1
* endif
c Check for interpolation in j
* if (abs(float(jh)-rj).lt.1.e-3) then
* j =jh
* jp1=jh
* else
j =min0(int(rj),n2-1)
jp1=j+1
* endif
c Check for interpolation in k
* if (abs(float(kh)-rk).lt.1.e-3) then
* k =kh
* kp1=kh
* else
k =min0(int(rk),n3-1)
kp1=k+1
* endif
if (k.eq.kp1) then
c no interpolation in k
if ((i.eq.ip1).and.(j.eq.jp1)) then
c no interpolation at all
int3dmvc=ar(i,j,k)
c print *,'int3dmvc 00: ',rid,rjd,rkd,int3dmvc
else
c horizontal interpolation only
if ((misdat.eq.ar(i ,j ,k )).or.
& (misdat.eq.ar(i ,jp1,k )).or.
& (misdat.eq.ar(ip1,j ,k )).or.
& (misdat.eq.ar(ip1,jp1,k ))) then
int3dmvc=misdat
else
frac0i=ri-float(i)
frac0j=rj-float(j)
frac1i=1.-frac0i
frac1j=1.-frac0j
int3dmvc = ar(i ,j ,k ) * frac1i * frac1j
& + ar(i ,jp1,k ) * frac1i * frac0j
& + ar(ip1,j ,k ) * frac0i * frac1j
& + ar(ip1,jp1,k ) * frac0i * frac0j
c print *,'int3dmvc 10: ',rid,rjd,rkd,int3dmvc
endif
endif
else
if ((i.eq.ip1).and.(j.eq.jp1)) then
c vertical interpolation only
if ((misdat.eq.ar(i ,j ,k )).or.
& (misdat.eq.ar(i ,j ,kp1))) then
int3dmvc=misdat
else
if (k-1.lt.1) then
klow=1
else if (k+2.gt.n3) then
klow=n3-3
else
klow=k-1
endif
call cubint(ar(i,j,klow),ar(i,j,klow+1),ar(i,j,klow+2),
& ar(i,j,klow+3),klow,rk,int3dmvc)
endif
else
c full 3d interpolation
if ((misdat.eq.ar(i ,j ,k )).or.
& (misdat.eq.ar(i ,jp1,k )).or.
& (misdat.eq.ar(ip1,j ,k )).or.
& (misdat.eq.ar(ip1,jp1,k )).or.
& (misdat.eq.ar(i ,j ,kp1)).or.
& (misdat.eq.ar(i ,jp1,kp1)).or.
& (misdat.eq.ar(ip1,j ,kp1)).or.
& (misdat.eq.ar(ip1,jp1,kp1))) then
int3dmvc=misdat
else
frac0i=ri-float(i)
frac0j=rj-float(j)
frac1i=1.-frac0i
frac1j=1.-frac0j
if (k-1.lt.1) then
klow=1
else if (k+2.gt.n3) then
klow=n3-3
else
klow=k-1
endif
do n=1,4
int2d(n) = ar(i ,j ,klow-1+n ) * frac1i * frac1j
& + ar(i ,jp1,klow-1+n ) * frac1i * frac0j
& + ar(ip1,j ,klow-1+n ) * frac0i * frac1j
& + ar(ip1,jp1,klow-1+n ) * frac0i * frac0j
enddo
call cubint(int2d(1),int2d(2),int2d(3),int2d(4),
& klow,rk,int3dmvc)
endif
endif
endif
end
real function int4d(ar,n1,n2,n3,n4,rid,rjd,rkd,rld)
c-----------------------------------------------------------------------
c Purpose:
c This subroutine interpolates a 4d-array to an arbitrary
c location within the grid.
c Arguments:
c ar real input surface pressure, define as ar(n1,n2,n3)
c n1,..,n4 int input dimensions of ar
c ri,..,rl real input grid location to be interpolated to
c History:
c-----------------------------------------------------------------------
 
c argument declarations
integer n1,n2,n3,n4
real ar(n1,n2,n3,n4), rid,rjd,rkd,rld
 
c local declarations
integer l,lp1,lh
real frac0l,frac1l,rl,int3d
 
c do linear interpolation in l-direction
rl=amax1(1.,amin1(float(n4),rld))
lh=nint(rl)
 
c Check for interpolation in l
* if (abs(float(lh)-rl).lt.1.e-3) then
* l =lh
* lp1=lh
* else
l =min0(int(rl),n4-1)
lp1=l+1
* endif
 
if (l.eq.lp1) then
c no interpolation in l
int4d=int3d(ar(1,1,1,l),n1,n2,n3,rid,rjd,rkd)
else
c interpolation in l
frac0l=rl-float(l)
frac1l=1.-frac0l
int4d = int3d(ar(1,1,1,l ),n1,n2,n3,rid,rjd,rkd) * frac1l
& + int3d(ar(1,1,1,lp1),n1,n2,n3,rid,rjd,rkd) * frac0l
endif
end
real function int4dm(ar,n1,n2,n3,n4,rid,rjd,rkd,rld,misdat)
c-----------------------------------------------------------------------
c Purpose:
c This subroutine interpolates a 4d-array to an arbitrary
c location within the grid. The interpolation includes the
c testing of the missing data flag 'misdat'.
c Arguments:
c ar real input surface pressure, define as ar(n1,n2,n3)
c n1,..,n4 int input dimensions of ar
c ri,..,rl real input grid location to be interpolated to
c misdat real input missing data flag (on if misdat<>0)
c Warning:
c This routine has not yet been seriously tested.
c History:
c-----------------------------------------------------------------------
 
c argument declarations
integer n1,n2,n3,n4
real ar(n1,n2,n3,n4), rid,rjd,rkd,rld, misdat
 
c local declarations
integer l,lp1,lh
real frac0l,frac1l,rl,rint0,rint1,int4d,int3dm
 
c check whether missing data checking is required
if (misdat.eq.0.) then
int4dm=int4d(ar,n1,n2,n3,n4,rid,rjd,rkd,rld)
return
endif
 
c do linear interpolation in l-direction
rl=amax1(1.,amin1(float(n4),rld))
lh=nint(rl)
 
c Check for interpolation in l
* if (abs(float(lh)-rl).lt.1.e-3) then
* l =lh
* lp1=lh
* else
l =min0(int(rl),n4-1)
lp1=l+1
* endif
 
if (l.eq.lp1) then
c no interpolation in l
int4dm = int3dm(ar(1,1,1,l),n1,n2,n3,rid,rjd,rkd,misdat)
else
c interpolation in l
frac0l=rl-float(l)
frac1l=1.-frac0l
rint0 = int3dm(ar(1,1,1,l ),n1,n2,n3,rid,rjd,rkd,misdat)
rint1 = int3dm(ar(1,1,1,lp1),n1,n2,n3,rid,rjd,rkd,misdat)
if ((rint0.eq.misdat).or.(rint1.eq.misdat)) then
int4dm = misdat
else
int4dm = rint0*frac1l + rint1*frac0l
endif
endif
end
real function int3dl(ar,n1,n2,n3,levels,rid,rjd,rkd)
c-----------------------------------------------------------------------
c Purpose:
c This subroutine interpolates a 3d-array to an arbitrary
c location within the grid. The vertical interpolation is linear
c in log(pressure).
c Arguments:
c ar real input surface pressure, define as ar(n1,n2,n3)
c n1,n2,n3 int input dimensions of ar
c levels real input array contains pressure levels for ar
c ri,rj,rk real input grid location to be interpolated to
c History:
c Based on int3d July 93
c-----------------------------------------------------------------------
 
c argument declarations
integer n1,n2,n3
real ar(n1,n2,n3), rid,rjd,rkd
real levels(n3)
 
c local declarations
real pval
integer i,j,k,ip1,jp1,kp1,ih,jh,kh
real frac0i,frac0j,frac0k,frac1i,frac1j,frac1k,ri,rj,rk
 
c do linear interpolation
ri=amax1(1.,amin1(float(n1),rid))
rj=amax1(1.,amin1(float(n2),rjd))
rk=amax1(1.,amin1(float(n3),rkd))
ih=nint(ri)
jh=nint(rj)
kh=nint(rk)
 
c Check for interpolation in i
* if (abs(float(ih)-ri).lt.1.e-3) then
* i =ih
* ip1=ih
* else
i =min0(int(ri),n1-1)
ip1=i+1
* endif
 
c Check for interpolation in j
* if (abs(float(jh)-rj).lt.1.e-3) then
* j =jh
* jp1=jh
* else
j =min0(int(rj),n2-1)
jp1=j+1
* endif
 
c Check for interpolation in k
* if (abs(float(kh)-rk).lt.1.e-3) then
* k =kh
* kp1=kh
* else
k =min0(int(rk),n3-1)
kp1=k+1
* endif
 
if (k.eq.kp1) then
c no interpolation in k
if ((i.eq.ip1).and.(j.eq.jp1)) then
c no interpolation at all
int3dl=ar(i,j,k)
c print *,'int3dl 00: ',rid,rjd,rkd,int3dl
else
c horizontal interpolation only
frac0i=ri-float(i)
frac0j=rj-float(j)
frac1i=1.-frac0i
frac1j=1.-frac0j
int3dl = ar(i ,j ,k ) * frac1i * frac1j
& + ar(i ,jp1,k ) * frac1i * frac0j
& + ar(ip1,j ,k ) * frac0i * frac1j
& + ar(ip1,jp1,k ) * frac0i * frac0j
c print *,'int3dl 10: ',rid,rjd,rkd,int3dl
endif
else
* frac0k=rk-float(k)
c calculate the pressure value to be interpolated to
pval=levels(int(rk))
> -(rk-aint(rk))*(levels(int(rk))-levels(int(rk)+1))
frac0k=log(levels(int(rk))/pval)
& /log(levels(int(rk))/levels(int(rk)+1))
frac1k=1.-frac0k
if ((i.eq.ip1).and.(j.eq.jp1)) then
c vertical interpolation only
int3dl = ar(i ,j ,k ) * frac1k
& + ar(i ,j ,kp1) * frac0k
c print *,'int3dl 01: ',rid,rjd,rkd,int3dl
else
c full 3d interpolation
frac0i=ri-float(i)
frac0j=rj-float(j)
frac1i=1.-frac0i
frac1j=1.-frac0j
int3dl = ar(i ,j ,k ) * frac1i * frac1j * frac1k
& + ar(i ,jp1,k ) * frac1i * frac0j * frac1k
& + ar(ip1,j ,k ) * frac0i * frac1j * frac1k
& + ar(ip1,jp1,k ) * frac0i * frac0j * frac1k
& + ar(i ,j ,kp1) * frac1i * frac1j * frac0k
& + ar(i ,jp1,kp1) * frac1i * frac0j * frac0k
& + ar(ip1,j ,kp1) * frac0i * frac1j * frac0k
& + ar(ip1,jp1,kp1) * frac0i * frac0j * frac0k
c print *,'int3dl 11: ',rid,rjd,rkd,int3dl
endif
endif
end
subroutine bcucof(y,y1,y2,y12,c)
c-----------------------------------------------------------------------
c Given arrays y,y1,y2 and y12, each of length 4, containing the
c function, gradients, and cross derivative at the four grid points
c of a rectangular grid cell (numbered counterclockwise from the
c lower left), and given d1 and d2, the length of the grid cell in
c the 1- and 2-directions, this routine returns the table c that is
c used by routine bcuint for biqubic interpolation.
c Source: Numerical Recipes, Fortran Version, p.99
c-----------------------------------------------------------------------
real c(4,4),y(4),y1(4),y2(4),y12(4),cl(16),x(16)
integer wt(16,16)
 
data wt/1,0,-3,2,4*0,-3,0,9,-6,2,0,-6,4,
> 8*0,3,0,-9,6,-2,0,6,-4,
> 10*0,9,-6,2*0,-6,4,2*0,3,-2,6*0,-9,6,2*0,6,-4,
> 4*0,1,0,-3,2,-2,0,6,-4,1,0,-3,2,8*0,-1,0,3,-2,1,0,-3,2,
> 10*0,-3,2,2*0,3,-2,6*0,3,-2,2*0,-6,4,2*0,3,-2,
> 0,1,-2,1,5*0,-3,6,-3,0,2,-4,2,9*0,3,-6,3,0,-2,4,-2,
> 10*0,-3,3,2*0,2,-2,2*0,-1,1,6*0,3,-3,2*0,-2,2,
> 5*0,1,-2,1,0,-2,4,-2,0,1,-2,1,9*0,-1,2,-1,0,1,-2,1,
> 10*0,1,-1,2*0,-1,1,6*0,-1,1,2*0,2,-2,2*0,-1,1/
 
real xx
integer i,j,k,l
 
do i=1,4 ! pack a temporary vector x
x(i)=y(i)
x(i+4)=y1(i)
x(i+8)=y2(i)
x(i+12)=y12(i)
enddo
 
do i=1,16 ! matrix multiply by the stored table
xx=0.
do k=1,16
xx=xx+wt(i,k)*x(k)
enddo
cl(i)=xx
enddo
 
l=0
do i=1,4 ! unpack the result into the output table
do j=1,4
l=l+1
c(i,j)=cl(l)
enddo
enddo
 
return
end
subroutine bcuint(y,y1,y2,y12,x1l,x2l,x1,x2,ansy)
c-----------------------------------------------------------------------
c Bicubic interpolation within a grid square. Input quantities are
c y,y1,y2,y12 (as described in bcucof); x1l and x1u, the lower and
c upper coordinates of the grid square in the 1-direction; x2l and
c x2u likewise for the 2-direction; and x1,x2, the coordinates of
c the desired point for the interpolation. The interplated function
c value is returned as ansy. This routine calls bcucof.
c Source: Numerical Recipes, Fortran Version, p.99/100
c !!! changed the proposed code !!!
c-----------------------------------------------------------------------
 
real y(4),y1(4),y2(4),y12(4),c(4,4)
real ansy,x1,x2,t,u
integer i,x1l,x2l
 
call bcucof(y,y1,y2,y12,c)
 
t=x1-real(x1l)
u=x2-real(x2l)
 
ansy=0.
 
do i=4,1,-1
ansy=t*ansy+((c(i,4)*u+c(i,3))*u+c(i,2))*u+c(i,1)
enddo
 
return
end
subroutine cubint(ya1,ya2,ya3,ya4,k,x,y)
c-----------------------------------------------------------------------
c Interface routine for SR polint for special case of cubic
c interpolation in index space, with xa=k,k+1,k+2,k+3
c-----------------------------------------------------------------------
integer k
real ya1,ya2,ya3,ya4,x,y
 
integer n
real xa(4),ya(4),dy
 
do n=1,4
xa(1)=real(k)
xa(2)=real(k+1)
xa(3)=real(k+2)
xa(4)=real(k+3)
ya(1)=ya1
ya(2)=ya2
ya(3)=ya3
ya(4)=ya4
enddo
 
call polint(xa,ya,4,x,y,dy)
 
return
end
subroutine polint(xa,ya,n,x,y,dy)
c-----------------------------------------------------------------------
c Given arrays xa and ya, each of length n, and given a value x, this
c routine returns a value y, and an error estimate dy. If P(x) is the
c polynomial of degree n-1 such that p(xa(i))=ya(i),i=1,...,n, then
c the returned value y=p(x)
c Source: Numerical Recipes, Fortran Version, p.82
c-----------------------------------------------------------------------
 
integer nmax,n
parameter(nmax=10)
real xa(n),ya(n),x,y,dy
integer i,m,ns
real c(nmax),d(nmax),den,dif,dift,ho,hp,w
 
ns=1
dif=abs(x-xa(1))
 
do i=1,n
dift=abs(x-xa(i))
if (dift.lt.dif) then
ns=i
dif=dift
endif
c(i)=ya(i)
d(i)=ya(i)
enddo
 
y=ya(ns)
ns=ns-1
do m=1,n-1
do i=1,n-m
ho=xa(i)-x
hp=xa(i+m)-x
w=c(i+1)-d(i)
den=ho-hp
den=w/den
d(i)=hp*den
c(i)=ho*den
enddo
if (2*ns.lt.n-m) then
dy=c(ns+1)
else
dy=d(ns)
ns=ns-1
endif
y=y+dy
enddo
 
return
end
subroutine filt2d (a,af,f1,f2,nx,ny,fil,misdat,
& iperx,ipery,ispol,inpol)
c =============================================================
c Apply a conservative diffusion operator onto the 2d field a,
c with full missing data checking.
c
c a real inp array to be filtered, dimensioned (nx,ny)
c af real out filtered array, dimensioned (nx,ny), can be
c equivalenced with array a in the calling routine
c f1 real workarray, dimensioned (nx+1,ny)
c f2 real workarray, dimensioned (nx,ny+1)
c fil real inp filter-coeff., 0<afil<=1. Maximum filtering with afil=1
c corresponds to one application of the linear filter.
c misdat real inp missing-data value, a(i,j)=misdat indicates that
c the corresponding value is not available. The
c misdat-checking can be switched off with with misdat=0.
c iperx int inp periodic boundaries in the x-direction (1=yes,0=no)
c ipery int inp periodic boundaries in the y-direction (1=yes,0=no)
c inpol int inp northpole at j=ny (1=yes,0=no)
c ispol int inp southpole at j=1 (1=yes,0=no)
c
c Christoph Schaer, 1993
c argument declaration
integer nx,ny
real a(nx,ny),af(nx,ny),f1(nx+1,ny),f2(nx,ny+1),fil,misdat
integer iperx,ipery,inpol,ispol
c local variable declaration
integer i,j,is
real fh
c compute constant fh
fh=0.125*fil
c compute fluxes in x-direction
if (misdat.eq.0.) then
do j=1,ny
do i=2,nx
f1(i,j)=a(i-1,j)-a(i,j)
enddo
enddo
else
do j=1,ny
do i=2,nx
if ((a(i,j).eq.misdat).or.(a(i-1,j).eq.misdat)) then
f1(i,j)=0.
else
f1(i,j)=a(i-1,j)-a(i,j)
endif
enddo
enddo
endif
if (iperx.eq.1) then
c do periodic boundaries in the x-direction
do j=1,ny
f1(1,j)=f1(nx,j)
f1(nx+1,j)=f1(2,j)
enddo
else
c set boundary-fluxes to zero
do j=1,ny
f1(1,j)=0.
f1(nx+1,j)=0.
enddo
endif
c compute fluxes in y-direction
if (misdat.eq.0.) then
do j=2,ny
do i=1,nx
f2(i,j)=a(i,j-1)-a(i,j)
enddo
enddo
else
do j=2,ny
do i=1,nx
if ((a(i,j).eq.misdat).or.(a(i,j-1).eq.misdat)) then
f2(i,j)=0.
else
f2(i,j)=a(i,j-1)-a(i,j)
endif
enddo
enddo
endif
c set boundary-fluxes to zero
do i=1,nx
f2(i,1)=0.
f2(i,ny+1)=0.
enddo
if (ipery.eq.1) then
c do periodic boundaries in the x-direction
do i=1,nx
f2(i,1)=f2(i,ny)
f2(i,ny+1)=f2(i,2)
enddo
endif
if (iperx.eq.1) then
if (ispol.eq.1) then
c do south-pole
is=(nx-1)/2
do i=1,nx
f2(i,1)=-f2(mod(i-1+is,nx)+1,2)
enddo
endif
if (inpol.eq.1) then
c do north-pole
is=(nx-1)/2
do i=1,nx
f2(i,ny+1)=-f2(mod(i-1+is,nx)+1,ny)
enddo
endif
endif
c compute flux-convergence -> filter
if (misdat.eq.0.) then
do j=1,ny
do i=1,nx
af(i,j)=a(i,j)+fh*(f1(i,j)-f1(i+1,j)+f2(i,j)-f2(i,j+1))
enddo
enddo
else
do j=1,ny
do i=1,nx
if (a(i,j).eq.misdat) then
af(i,j)=misdat
else
af(i,j)=a(i,j)+fh*(f1(i,j)-f1(i+1,j)+f2(i,j)-f2(i,j+1))
endif
enddo
enddo
endif
end
 
subroutine pipo(var3d,p3d,lev,var,nx,ny,nz,mdv,mode)
C ====================================================
 
C Interpolates the 3d variable var3d on the pressure surface
C defined by lev. The interpolated field is returned as var.
C p3d denotes the 3d pressure array.
C mode determines the way of vertical interpolation:
C mode=0 is for linear interpolation
C mode=1 is for logarithmic interpolation
 
integer nx,ny,nz,mode
real lev,mdv
real var3d(nx,ny,nz),p3d(nx,ny,nz),var(nx,ny)
 
integer i,j,k
real kind
real int3dm
 
do i=1,nx
do j=1,ny
 
kind=0.
do k=1,nz-1
if ((p3d(i,j,k).ge.lev).and.(p3d(i,j,k+1).le.lev)) then
kind=float(k)+(p3d(i,j,k)-lev)/
> (p3d(i,j,k)-p3d(i,j,k+1))
goto 100
endif
enddo
100 continue
 
if (kind.eq.0.) then
var(i,j)=mdv
else
var(i,j)=int3dm(var3d,nx,ny,nz,float(i),float(j),kind,mdv)
endif
 
enddo
enddo
 
return
end
 
subroutine thipo(var3d,th3d,lev,var,nx,ny,nz,mdv,mode)
C ======================================================
 
C Interpolates the 3d variable var3d on the isentropic surface
C defined by lev. The interpolated field is returned as var.
C th3d denotes the 3d theta array.
C mode determines the way of vertical interpolation:
C mode=0 is for linear interpolation
C mode=1 is for logarithmic interpolation
 
integer nx,ny,nz,mode
real lev,mdv
real var3d(nx,ny,nz),th3d(nx,ny,nz),var(nx,ny)
 
integer i,j,k
real kind
real int3dm
 
do i=1,nx
do j=1,ny
 
kind=0
do k=1,nz-1
if ((th3d(i,j,k).le.lev).and.(th3d(i,j,k+1).ge.lev)) then
kind=float(k)+(th3d(i,j,k)-lev)/
> (th3d(i,j,k)-th3d(i,j,k+1))
goto 100
endif
enddo
100 continue
 
if (kind.eq.0) then
var(i,j)=mdv
else
var(i,j)=int3dm(var3d,nx,ny,nz,float(i),float(j),kind,mdv)
endif
 
enddo
enddo
 
return
end
/tags/1.0/main.mk
0,0 → 1,56
### name of the executable that will be produced
PROG = $(INSTALLDIR)/3d_labelling_and_fold_id.exe
 
DATE=$(shell date +"%Y%m%d")
 
# complete list of all f90 source files
SRCS = mo_f2kcli.f90 netcdf_tools.f90 libipo.f 3d_labelling_and_fold_id.f90
 
# --------------------------------------------------------------------
OBJS = mo_f2kcli.o netcdf_tools.o libipo.o 3d_labelling_and_fold_id.o
 
all: $(PROG)
 
$(PROG): $(OBJS)
$(F90) $(F90FLAGS) $(OBJS) $(LIBS) -o $@
 
list:
@echo "------------------------------------------------"
@echo "SRCS = $(SRCS)"
@echo "------------------------------------------------"
@echo
@echo "------------------------------------------------"
@echo "OBJS = $(OBJS)"
@echo "------------------------------------------------"
 
zip:
zip -r 3d_labelling_and_fold_id_$(DATE).zip *
clean:
rm -f *.mod
rm -f *.o
 
distclean: clean
rm -f $(PROG)
 
help:
@echo ''
@echo ' posible make targets:'
@echo ' -------------------------------------------------------'
@echo ' make : build all'
@echo ' make all : (= make)'
@echo ''
@echo ' make clean : delete libs/mod/obj'
@echo ' gmake distclean : clean up distribution'
@echo ''
@echo ' gmake zip : zip source-code'
@echo ' -------------------------------------------------------'
@echo ''
 
 
%.o: %.f90
$(F90) $(F90FLAGS) $(LIBS) $(INCLUDES) -c $<
%.o: %.f
$(FC) $(FFLAGS) $(LIBS) $(INCLUDES) -c $<
 
## ------------------------------------------------------------------
/tags/1.0/mo_f2kcli.f90
0,0 → 1,205
! F2KCLI : Fortran 200x Command Line Interface
! copyright Interactive Software Services Ltd. 2001
! For conditions of use see manual.txt
!
! Platform : Unix/Linux
! Compiler : Any Fortran 9x compiler supporting IARGC/GETARG
! which counts the first true command line argument
! after the program name as argument number one.
! (Excludes compilers which require a special USE
! statement to make IARGC/GETARG available).
! To compile : f90 -c f2kcli.f90
! (exact compiler name will vary)
! Implementer : Lawson B. Wakefield, I.S.S. Ltd.
! Date : February 2001
!
MODULE MO_F2KCLI
#if defined (NAG)
USE f90_unix
!#else
! EXTERNAL :: GETARG, IARG, IARGC
#endif
! INTRINSIC :: LEN_TRIM, PRESENT, LEN
! PRIVATE :: LEN_TRIM, PRESENT, LEN
! INTEGER :: IARGC
!
CONTAINS
!
SUBROUTINE GET_COMMAND(COMMAND,LENGTH,STATUS)
!
! Description. Returns the entire command by which the program was
! invoked.
!
! Class. Subroutine.
!
! Arguments.
! COMMAND (optional) shall be scalar and of type default character.
! It is an INTENT(OUT) argument. It is assigned the entire command
! by which the program was invoked. If the command cannot be
! determined, COMMAND is assigned all blanks.
! LENGTH (optional) shall be scalar and of type default integer. It is
! an INTENT(OUT) argument. It is assigned the significant length
! of the command by which the program was invoked. The significant
! length may include trailing blanks if the processor allows commands
! with significant trailing blanks. This length does not consider any
! possible truncation or padding in assigning the command to the
! COMMAND argument; in fact the COMMAND argument need not even be
! present. If the command length cannot be determined, a length of
! 0 is assigned.
! STATUS (optional) shall be scalar and of type default integer. It is
! an INTENT(OUT) argument. It is assigned the value 0 if the
! command retrieval is sucessful. It is assigned a processor-dependent
! non-zero value if the command retrieval fails.
!
CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: COMMAND
INTEGER , INTENT(OUT), OPTIONAL :: LENGTH
INTEGER , INTENT(OUT), OPTIONAL :: STATUS
!
INTEGER :: IARG,NARG,IPOS
INTEGER , SAVE :: LENARG
CHARACTER(LEN=2000), SAVE :: ARGSTR
LOGICAL , SAVE :: GETCMD = .TRUE.
!
! Under Unix we must reconstruct the command line from its constituent
! parts. This will not be the original command line. Rather it will be
! the expanded command line as generated by the shell.
!
IF (GETCMD) THEN
NARG = IARGC()
IF (NARG > 0) THEN
IPOS = 1
DO IARG = 1,NARG
CALL GETARG(IARG,ARGSTR(IPOS:))
LENARG = LEN_TRIM(ARGSTR)
IPOS = LENARG + 2
IF (IPOS > LEN(ARGSTR)) EXIT
END DO
ELSE
ARGSTR = ' '
LENARG = 0
ENDIF
GETCMD = .FALSE.
ENDIF
IF (PRESENT(COMMAND)) COMMAND = ARGSTR
IF (PRESENT(LENGTH)) LENGTH = LENARG
IF (PRESENT(STATUS)) STATUS = 0
RETURN
END SUBROUTINE GET_COMMAND
!
INTEGER FUNCTION COMMAND_ARGUMENT_COUNT()
!
! Description. Returns the number of command arguments.
!
! Class. Inquiry function
!
! Arguments. None.
!
! Result Characteristics. Scalar default integer.
!
! Result Value. The result value is equal to the number of command
! arguments available. If there are no command arguments available
! or if the processor does not support command arguments, then
! the result value is 0. If the processor has a concept of a command
! name, the command name does not count as one of the command
! arguments.
!
COMMAND_ARGUMENT_COUNT = IARGC()
RETURN
END FUNCTION COMMAND_ARGUMENT_COUNT
!
SUBROUTINE GET_COMMAND_ARGUMENT(NUMBER,VALUE,LENGTH,STATUS)
!
! Description. Returns a command argument.
!
! Class. Subroutine.
!
! Arguments.
! NUMBER shall be scalar and of type default integer. It is an
! INTENT(IN) argument. It specifies the number of the command
! argument that the other arguments give information about. Useful
! values of NUMBER are those between 0 and the argument count
! returned by the COMMAND_ARGUMENT_COUNT intrinsic.
! Other values are allowed, but will result in error status return
! (see below). Command argument 0 is defined to be the command
! name by which the program was invoked if the processor has such
! a concept. It is allowed to call the GET_COMMAND_ARGUMENT
! procedure for command argument number 0, even if the processor
! does not define command names or other command arguments.
! The remaining command arguments are numbered consecutively from
! 1 to the argument count in an order determined by the processor.
! VALUE (optional) shall be scalar and of type default character.
! It is an INTENT(OUT) argument. It is assigned the value of the
! command argument specified by NUMBER. If the command argument value
! cannot be determined, VALUE is assigned all blanks.
! LENGTH (optional) shall be scalar and of type default integer.
! It is an INTENT(OUT) argument. It is assigned the significant length
! of the command argument specified by NUMBER. The significant
! length may include trailing blanks if the processor allows command
! arguments with significant trailing blanks. This length does not
! consider any possible truncation or padding in assigning the
! command argument value to the VALUE argument; in fact the
! VALUE argument need not even be present. If the command
! argument length cannot be determined, a length of 0 is assigned.
! STATUS (optional) shall be scalar and of type default integer.
! It is an INTENT(OUT) argument. It is assigned the value 0 if
! the argument retrieval is sucessful. It is assigned a
! processor-dependent non-zero value if the argument retrieval fails.
!
! NOTE
! One possible reason for failure is that NUMBER is negative or
! greater than COMMAND_ARGUMENT_COUNT().
!
INTEGER , INTENT(IN) :: NUMBER
CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: VALUE
INTEGER , INTENT(OUT), OPTIONAL :: LENGTH
INTEGER , INTENT(OUT), OPTIONAL :: STATUS
!
! A temporary variable for the rare case case where LENGTH is
! specified but VALUE is not. An arbitrary maximum argument length
! of 1000 characters should cover virtually all situations.
!
CHARACTER(LEN=1000) :: TMPVAL
!
! Possible error codes:
! 1 = Argument number is less than minimum
! 2 = Argument number exceeds maximum
!
IF (NUMBER < 0) THEN
IF (PRESENT(VALUE )) VALUE = ' '
IF (PRESENT(LENGTH)) LENGTH = 0
IF (PRESENT(STATUS)) STATUS = 1
RETURN
ELSE IF (NUMBER > IARGC()) THEN
IF (PRESENT(VALUE )) VALUE = ' '
IF (PRESENT(LENGTH)) LENGTH = 0
IF (PRESENT(STATUS)) STATUS = 2
RETURN
END IF
!
! Get the argument if VALUE is present
!
IF (PRESENT(VALUE)) CALL GETARG(NUMBER,VALUE)
!
! The LENGTH option is fairly pointless under Unix.
! Trailing spaces can only be specified using quotes.
! Since the command line has already been processed by the
! shell before the application sees it, we have no way of
! knowing the true length of any quoted arguments. LEN_TRIM
! is used to ensure at least some sort of meaningful result.
!
IF (PRESENT(LENGTH)) THEN
IF (PRESENT(VALUE)) THEN
LENGTH = LEN_TRIM(VALUE)
ELSE
CALL GETARG(NUMBER,TMPVAL)
LENGTH = LEN_TRIM(TMPVAL)
END IF
END IF
!
! Since GETARG does not return a result code, assume success
!
IF (PRESENT(STATUS)) STATUS = 0
RETURN
END SUBROUTINE GET_COMMAND_ARGUMENT
!
END MODULE MO_F2KCLI
/tags/1.0/netcdf_tools.f90
0,0 → 1,654
module netcdf_tools
 
USE netcdf
implicit none
 
INTRINSIC :: TRIM, ADJUSTL, NINT, SIZE
 
interface read_file
module procedure read_file_1d
module procedure read_file_2d
module procedure read_file_3d
module procedure read_file_4d
end interface
 
contains
 
SUBROUTINE inquire_file(fname, &
x_name, y_name, z_name,t_name, &
nx, ny, nz, nt)
 
 
IMPLICIT NONE
 
 
! I/O
CHARACTER(LEN=*), INTENT(IN) :: fname ! filename
CHARACTER(LEN=*), INTENT(IN) :: x_name ! name of dimension in file dimension
CHARACTER(LEN=*), INTENT(IN) :: y_name ! name of dimension in file dimension
CHARACTER(LEN=*), INTENT(IN) :: z_name ! name of dimension in file dimension
CHARACTER(LEN=*), INTENT(IN) :: t_name ! name of dimension in file dimension
INTEGER, INTENT(OUT) :: nx ! file dimension lenght
INTEGER, INTENT(OUT) :: ny ! file dimension length
INTEGER, INTENT(OUT) :: nz ! file dimension length
INTEGER, INTENT(OUT) :: nt ! file dimension length
 
! LOCAL
INTEGER :: status
CHARACTER(LEN=*), PARAMETER :: substr = 'inquire_file'
INTEGER,SAVE :: ncid ! netCDF-ID
INTEGER :: dimid, varid
LOGICAL :: file_exists ! checking existence of file
 
 
INTEGER, DIMENSION(:), ALLOCATABLE :: date_file
 
CHARACTER(LEN=30) :: name_dim ! line
 
 
INQUIRE(FILE=TRIM(fname), EXIST=file_exists)
IF (.not.file_exists) THEN
WRITE(*,*) 'File ',TRIM(fname),' does NOT exist! skipping....'
STOP
ENDIF
 
! OPEN FILE
CALL NFERR( status, &
nf90_open(TRIM(fname), NF90_NOWRITE, ncid) &
,1)
! latitude dimension check
CALL NFERR( status, &
nf90_inq_dimid(ncid, x_name, dimid ) &
,2)
CALL NFERR( status, &
nf90_Inquire_Dimension(ncid, dimid, name_dim, nx ) &
,3)
! longitude dimension check
CALL NFERR( status, &
nf90_inq_dimid(ncid, y_name, dimid ) &
,4)
CALL NFERR( status, &
nf90_Inquire_Dimension(ncid, dimid, name_dim, ny ) &
,5)
! vertical dimension check
CALL NFERR( status, &
nf90_inq_dimid(ncid, z_name, dimid ) &
,6)
CALL NFERR( status, &
nf90_Inquire_Dimension(ncid, dimid, name_dim, nz ) &
,7)
! time dimension check
CALL NFERR( status, &
nf90_inq_dimid(ncid, t_name, dimid ) &
,8)
CALL NFERR( status, &
nf90_Inquire_Dimension(ncid, dimid, name_dim, nt ) &
,9)
 
 
!CLOSE FILE
CALL NFERR( status, &
nf90_close(ncid) &
,14)
 
! RETURN
status = 0
 
END SUBROUTINE inquire_file
 
! ------------------------------------------------------------------
 
! ------------------------------------------------------------------------
SUBROUTINE read_file_1D(fname, varname, data_file)
 
 
IMPLICIT NONE
 
! I/O
CHARACTER(LEN=*), INTENT(IN) :: fname ! filename
CHARACTER(LEN=*), INTENT(IN) :: varname ! variable name
REAL, DIMENSION(:), INTENT(OUT):: data_file ! INTENT(OUT)
! LOCAL
INTEGER :: status
CHARACTER(LEN=*), PARAMETER :: substr = 'read_file_1d'
INTEGER,SAVE :: ncid ! netCDF-ID
INTEGER :: dimid, varid
 
! OPEN FILE
CALL NFERR( status, &
nf90_open(TRIM(fname), NF90_NOWRITE, ncid) &
,21)
 
CALL NFERR( status, &
nf90_inq_varid(ncid, TRIM(varname), varid ) &
,22)
IF (status.ne.0) then
write(*,*) "variable not found in NetCDF file, skipping"
STOP
ENDIF
 
CALL NFERR( status, &
nf90_get_var(ncid, varid, data_file ) &
,23)
 
!CLOSE FILE
CALL NFERR( status, &
nf90_close(ncid) &
,24)
 
! RETURN
status = 0
END SUBROUTINE read_file_1D
! ------------------------------------------------------------------------
! ------------------------------------------------------------------------
SUBROUTINE read_file_2D(fname, varname, data_file)
 
 
IMPLICIT NONE
 
! I/O
CHARACTER(LEN=*), INTENT(IN) :: fname ! filename
CHARACTER(LEN=*), INTENT(IN) :: varname ! variable name
REAL, DIMENSION(:,:), INTENT(OUT):: data_file ! INTENT(OUT)
! LOCAL
INTEGER :: status
CHARACTER(LEN=*), PARAMETER :: substr = 'read_file_2d'
INTEGER,SAVE :: ncid ! netCDF-ID
INTEGER :: dimid, varid
 
! OPEN FILE
CALL NFERR( status, &
nf90_open(TRIM(fname), NF90_NOWRITE, ncid) &
,21)
 
CALL NFERR( status, &
nf90_inq_varid(ncid, TRIM(varname), varid ) &
,22)
IF (status.ne.0) then
write(*,*) "variable not found in NetCDF file, skipping"
STOP
ENDIF
 
CALL NFERR( status, &
nf90_get_var(ncid, varid, data_file ) &
,23)
 
!CLOSE FILE
CALL NFERR( status, &
nf90_close(ncid) &
,24)
 
! RETURN
status = 0
END SUBROUTINE read_file_2D
! ------------------------------------------------------------------------
! ------------------------------------------------------------------------
SUBROUTINE read_file_3D(fname, varname, data_file)
 
 
IMPLICIT NONE
 
! I/O
CHARACTER(LEN=*), INTENT(IN) :: fname ! filename
CHARACTER(LEN=*), INTENT(IN) :: varname ! variable name
REAL, DIMENSION(:,:,:), INTENT(OUT):: data_file ! INTENT(OUT)
! LOCAL
INTEGER :: status
CHARACTER(LEN=*), PARAMETER :: substr = 'read_file_3d'
INTEGER,SAVE :: ncid ! netCDF-ID
INTEGER :: dimid, varid
 
! OPEN FILE
CALL NFERR( status, &
nf90_open(TRIM(fname), NF90_NOWRITE, ncid) &
,21)
 
CALL NFERR( status, &
nf90_inq_varid(ncid, TRIM(varname), varid ) &
,22)
IF (status.ne.0) then
write(*,*) "variable not found in NetCDF file, skipping"
STOP
ENDIF
 
CALL NFERR( status, &
nf90_get_var(ncid, varid, data_file ) &
,23)
 
!CLOSE FILE
CALL NFERR( status, &
nf90_close(ncid) &
,24)
 
! RETURN
status = 0
END SUBROUTINE read_file_3D
! ------------------------------------------------------------------------
! ------------------------------------------------------------------------
SUBROUTINE read_file_4D(fname, varname, data_file)
 
 
IMPLICIT NONE
 
! I/O
CHARACTER(LEN=*), INTENT(IN) :: fname ! filename
CHARACTER(LEN=*), INTENT(IN) :: varname ! variable name
REAL, DIMENSION(:,:,:,:), INTENT(OUT):: data_file ! INTENT(OUT)
! LOCAL
INTEGER :: status
CHARACTER(LEN=*), PARAMETER :: substr = 'read_file_4d'
INTEGER,SAVE :: ncid ! netCDF-ID
INTEGER :: dimid, varid
 
! OPEN FILE
CALL NFERR( status, &
nf90_open(TRIM(fname), NF90_NOWRITE, ncid) &
,21)
 
CALL NFERR( status, &
nf90_inq_varid(ncid, TRIM(varname), varid ) &
,22)
IF (status.ne.0) then
write(*,*) "variable not found in NetCDF file, skipping"
STOP
ENDIF
 
CALL NFERR( status, &
nf90_get_var(ncid, varid, data_file ) &
,23)
 
!CLOSE FILE
CALL NFERR( status, &
nf90_close(ncid) &
,24)
 
! RETURN
status = 0
END SUBROUTINE read_file_4D
! ------------------------------------------------------------------------
 
! ------------------------------------------------------------------------
SUBROUTINE read_att(fname, varname, attname, str)
 
 
IMPLICIT NONE
 
! I/O
CHARACTER(LEN=*), INTENT(IN) :: fname ! filename
CHARACTER(LEN=*), INTENT(IN) :: varname ! variable name
CHARACTER(LEN=*), INTENT(IN) :: attname ! attribute name
CHARACTER(LEN=*), INTENT(OUT) :: str ! date string
! LOCAL
INTEGER :: status
CHARACTER(LEN=*), PARAMETER :: substr = 'read_startdate'
INTEGER,SAVE :: ncid ! netCDF-ID
INTEGER :: dimid, varid
 
! OPEN FILE
CALL NFERR( status, &
nf90_open(TRIM(fname), NF90_NOWRITE, ncid) &
,21)
 
CALL NFERR( status, &
nf90_inq_varid(ncid, TRIM(varname), varid ) &
,22)
IF (status.ne.0) then
write(*,*) "variable not found in NetCDF file, skipping"
STOP
ENDIF
 
CALL NFERR( status, &
nf90_get_att(ncid, varid, TRIM(attname), str ) &
,23)
 
!CLOSE FILE
CALL NFERR( status, &
nf90_close(ncid) &
,24)
 
! RETURN
status = 0
END SUBROUTINE read_att
! ------------------------------------------------------------------------
 
! ------------------------------------------------------------------------
SUBROUTINE nc_dump(fname, x_data, y_data,z_data,t_data, &
x_units,y_units, z_units,t_units, aps, ak, bk, &
label,fold, sfold, mfold, dfold, tp, dp, pmin, pmax)
 
USE netcdf
 
IMPLICIT NONE
 
CHARACTER(LEN=*), INTENT(IN) :: fname
REAL, DIMENSION(:), INTENT(IN) :: x_data, y_data,z_data,t_data
CHARACTER(LEN=*), INTENT(IN) :: x_units,y_units, z_units,t_units
REAL, DIMENSION(:,:,:), INTENT(IN) :: aps
REAL, DIMENSION(:), INTENT(IN) :: ak,bk
REAL, DIMENSION(:,:,:,:), INTENT(IN) :: label
REAL, DIMENSION(:,:,:), INTENT(IN) :: fold
REAL, DIMENSION(:,:,:), INTENT(IN) :: sfold
REAL, DIMENSION(:,:,:), INTENT(IN) :: mfold
REAL, DIMENSION(:,:,:), INTENT(IN) :: dfold
REAL, DIMENSION(:,:,:), INTENT(IN) :: dp
REAL, DIMENSION(:,:,:), INTENT(IN) :: tp
REAL, DIMENSION(:,:,:), INTENT(IN) :: pmin
REAL, DIMENSION(:,:,:), INTENT(IN) :: pmax
 
! LOCAL
INTEGER :: status
CHARACTER(LEN=*), PARAMETER :: substr = 'nc_dump'
INTEGER :: ncid ! netCDF-ID
INTEGER :: nlon,nlat,nlev,ntime
INTEGER :: dimid_lat, dimid_lon, dimid_lev, dimid_ilev, dimid_time
INTEGER :: varid_lat, varid_lon, varid_lev, varid_ilev, varid_time
INTEGER :: varid_aps, varid_hyam, varid_hybm, varid_label
INTEGER :: varid_fold, varid_sfold, varid_mfold, varid_dfold, varid_tp, varid_dp
INTEGER :: varid_pmin, varid_pmax
 
!
CHARACTER(LEN=8) :: date
CHARACTER(LEN=10) :: time
CHARACTER(LEN=5) :: zone
 
 
nlon = SIZE(x_data)
nlat = SIZE(y_data)
nlev = SIZE(z_data)
ntime = SIZE(t_data)
 
! CREATE NEW FILE
CALL NFERR(status, &
nf90_create(TRIM(fname), NF90_CLOBBER, ncid) &
,51)
 
! ADD GLOBALE ATTRIBUTES
! - VERSION
CALL NFERR(status, &
nf90_put_att(ncid, NF90_GLOBAL, 'contact:', &
'Andrea Pozzer, MPIC, Mainz') &
,52)
! - DATE AND TIME
CALL DATE_AND_TIME(date, time, zone)
CALL NFERR(status, &
nf90_put_att(ncid, NF90_GLOBAL, 'date', date) &
,53)
CALL NFERR(status, &
nf90_put_att(ncid, NF90_GLOBAL, 'time', TRIM(time)//TRIM(zone)) &
,54)
! DEFINE DIMENSIONS
CALL NFERR(status, &
nf90_def_dim(ncid, 'lon', nlon, dimid_lon) &
,56)
CALL NFERR(status, &
nf90_def_dim(ncid, 'lat', nlat, dimid_lat) &
,57)
CALL NFERR(status, &
nf90_def_dim(ncid, 'lev', nlev, dimid_lev) &
,58)
CALL NFERR(status, &
nf90_def_dim(ncid, 'time', NF90_UNLIMITED, dimid_time) &
,59)
 
! DEFINE COORDINATE VARIABLES WITH ATTRIBUTES
CALL NFERR(status, &
nf90_def_var(ncid, 'lon', NF90_FLOAT, (/ dimid_lon /), varid_lon) &
,60)
CALL NFERR(status, &
nf90_put_att(ncid, varid_lon, 'long_name', 'longitude') &
,61)
CALL NFERR(status, &
nf90_put_att(ncid, varid_lon, 'units', x_units) &
,62)
 
CALL NFERR(status, &
nf90_def_var(ncid, 'lat', NF90_FLOAT, (/ dimid_lat /), varid_lat) &
,63)
CALL NFERR(status, &
nf90_put_att(ncid, varid_lat, 'long_name', 'latitude') &
,63)
CALL NFERR(status, &
nf90_put_att(ncid, varid_lat, 'units', y_units) &
,64)
 
CALL NFERR(status, &
nf90_def_var(ncid, 'lev', NF90_FLOAT, (/ dimid_lev /), varid_lev) &
,65)
CALL NFERR(status, &
nf90_put_att(ncid, varid_lev, 'long_name', 'level index') &
,66)
CALL NFERR(status, &
nf90_put_att(ncid, varid_lev, 'units', z_units) &
,67)
CALL NFERR(status, &
nf90_def_var(ncid, 'time', NF90_FLOAT, (/ dimid_time /), varid_time) &
,68)
CALL NFERR(status, &
nf90_put_att(ncid, varid_time, 'long_name', 'time') &
,69)
CALL NFERR(status, &
nf90_put_att(ncid, varid_time, 'units', t_units) &
,70)
 
! DEFINE VARIABLES
! - aps
CALL NFERR(status, &
nf90_def_var(ncid, 'APS', NF90_FLOAT &
, (/ dimid_lon, dimid_lat, dimid_time /), varid_aps) &
,78)
CALL NFERR(status, &
nf90_put_att(ncid, varid_aps, 'long_name' &
,'Surface Pressure') &
,79)
! - ak
CALL NFERR(status, &
nf90_def_var(ncid, 'hyam', NF90_FLOAT &
, (/ dimid_lev /), varid_hyam) &
,78)
CALL NFERR(status, &
nf90_put_att(ncid, varid_hyam, 'long_name' &
,'"hybrid A coefficient at layer midpoints') &
,79)
! - bk
CALL NFERR(status, &
nf90_def_var(ncid, 'hybm', NF90_FLOAT &
, (/ dimid_lev /), varid_hybm) &
,78)
CALL NFERR(status, &
nf90_put_att(ncid, varid_hybm, 'long_name' &
,'hybrid B coefficient at layer midpoints') &
,79)
! - label
CALL NFERR(status, &
nf90_def_var(ncid, 'label', NF90_FLOAT &
, (/ dimid_lon, dimid_lat, dimid_lev, dimid_time /), varid_label) &
,78)
CALL NFERR(status, &
nf90_put_att(ncid, varid_label, 'long_name' &
,'TR=1,ST=2,ST.CUTOFF=3,TR.TR.CUTOFF=4,PV.BLOB=5') &
,79)
! - fold
CALL NFERR(status, &
nf90_def_var(ncid, 'fold', NF90_FLOAT &
, (/ dimid_lon, dimid_lat, dimid_time /), varid_fold) &
,78)
CALL NFERR(status, &
nf90_put_att(ncid, varid_fold, 'long_name' &
, 'Tropopause folding: 1=YES, 2=NO') &
,79)
! - sfold
CALL NFERR(status, &
nf90_def_var(ncid, 'sfold', NF90_FLOAT &
, (/ dimid_lon, dimid_lat, dimid_time /), varid_sfold) &
,78)
CALL NFERR(status, &
nf90_put_att(ncid, varid_sfold, 'long_name' &
, 'Shallow tropopause folding (50< >200 hPa) : 1=YES, 2=NO') &
,79)
! - mfold
CALL NFERR(status, &
nf90_def_var(ncid, 'mfold', NF90_FLOAT &
, (/ dimid_lon, dimid_lat, dimid_time /), varid_mfold) &
,78)
CALL NFERR(status, &
nf90_put_att(ncid, varid_mfold, 'long_name' &
, 'Medium Tropopause folding (200< >350 hPa): 1=YES, 2=NO') &
,79)
! - dfold
CALL NFERR(status, &
nf90_def_var(ncid, 'dfold', NF90_FLOAT &
, (/ dimid_lon, dimid_lat, dimid_time /), varid_dfold) &
,78)
CALL NFERR(status, &
nf90_put_att(ncid, varid_dfold, 'long_name' &
, 'Deep Tropopause folding (>=350 hPa): 1=YES, 2=NO') &
,79)
! - tropopause
CALL NFERR(status, &
nf90_def_var(ncid, 'tp', NF90_FLOAT &
, (/ dimid_lon, dimid_lat, dimid_time /), varid_tp) &
,78)
CALL NFERR(status, &
nf90_put_att(ncid, varid_tp, 'long_name' &
, 'Tropopause') &
,79)
CALL NFERR(status, &
nf90_put_att(ncid, varid_tp, 'units' &
, 'hPa') &
,79)
! - folding extension
CALL NFERR(status, &
nf90_def_var(ncid, 'dp', NF90_FLOAT &
, (/ dimid_lon, dimid_lat, dimid_time /), varid_dp) &
,78)
CALL NFERR(status, &
nf90_put_att(ncid, varid_dp, 'long_name' &
, 'vertical extend of folding') &
,79)
CALL NFERR(status, &
nf90_put_att(ncid, varid_dp, 'units' &
, 'hPa') &
,79)
CALL NFERR(status, &
nf90_def_var(ncid, 'pmin', NF90_FLOAT &
, (/ dimid_lon, dimid_lat, dimid_time /), varid_pmin) &
,78)
CALL NFERR(status, &
nf90_put_att(ncid, varid_pmin, 'long_name' &
, 'vertical extend of folding') &
,79)
CALL NFERR(status, &
nf90_put_att(ncid, varid_pmin, 'units' &
, 'hPa') &
,79)
CALL NFERR(status, &
nf90_def_var(ncid, 'pmax', NF90_FLOAT &
, (/ dimid_lon, dimid_lat, dimid_time /), varid_pmax) &
,78)
CALL NFERR(status, &
nf90_put_att(ncid, varid_pmax, 'long_name' &
, 'vertical extend of folding') &
,79)
CALL NFERR(status, &
nf90_put_att(ncid, varid_pmax, 'units' &
, 'hPa') &
,79)
 
! SWITCH MODUS
CALL NFERR(status, &
nf90_enddef(ncid) &
,82)
 
! SAVE COORDINATE VARIBLES
CALL NFERR(status, &
nf90_put_var(ncid, varid_lon, x_data) &
,83)
CALL NFERR(status, &
nf90_put_var(ncid, varid_lat, y_data) &
,84)
CALL NFERR(status, &
nf90_put_var(ncid, varid_lev, z_data) &
,85)
CALL NFERR(status, &
nf90_put_var(ncid, varid_time, t_data) &
,86)
 
CALL NFERR(status, &
nf90_put_var(ncid, varid_aps, aps) &
,89)
CALL NFERR(status, &
nf90_put_var(ncid, varid_hyam, ak) &
,89)
CALL NFERR(status, &
nf90_put_var(ncid, varid_hybm, bk) &
,89)
CALL NFERR(status, &
nf90_put_var(ncid, varid_label, label) &
,89)
CALL NFERR(status, &
nf90_put_var(ncid, varid_fold, fold) &
,89)
CALL NFERR(status, &
nf90_put_var(ncid, varid_sfold, sfold) &
,89)
CALL NFERR(status, &
nf90_put_var(ncid, varid_mfold, mfold) &
,89)
CALL NFERR(status, &
nf90_put_var(ncid, varid_dfold, dfold) &
,89)
CALL NFERR(status, &
nf90_put_var(ncid, varid_tp, tp) &
,89)
CALL NFERR(status, &
nf90_put_var(ncid, varid_dp, dp) &
,89)
CALL NFERR(status, &
nf90_put_var(ncid, varid_pmin, pmin) &
,89)
CALL NFERR(status, &
nf90_put_var(ncid, varid_pmax, pmax) &
,89)
 
! CLOSE FILE
CALL NFERR(status, &
nf90_close(ncid) &
,90)
 
END SUBROUTINE nc_dump
 
 
! ------------------------------------------------------------------
SUBROUTINE NFERR(status,command,pos)
 
IMPLICIT NONE
 
! I/O
INTEGER, INTENT(OUT) :: status
INTEGER, INTENT(IN) :: command
INTEGER, INTENT(IN) :: pos
 
status=command
IF (status /= NF90_NOERR) THEN
WRITE(*,*) 'netCDF ERROR at position: ', pos
WRITE(*,*) 'netCDF ERROR status : ',status
WRITE(*,*) 'netCDF ERROR : ',nf90_strerror(status)
END IF
 
END SUBROUTINE NFERR
! ------------------------------------------------------------------
 
 
 
end module netcdf_tools