!
module addfld_module
!
! This software is part of the NCAR TIE-GCM. Use is governed by the
! Open Source Academic Research License Agreement contained in the file
! tiegcmlicense.txt.
!
! Save fields to secondary history (B. Foster: November, 2005)
! This module replaces the old addfsech routines.
!
! 7/13/06 btf: Using fsechist%data associated status to determine
! firstcall, rather than unames.
!
use params_module
implicit none
contains
!-----------------------------------------------------------------------
subroutine addfld(name,long_name,units,
| f,dname1,lb1,ub1,dname2,lb2,ub2,idx)
!
! Purpose:
! Save user-provided data f to secondary history field structure.
!
! Basically, the user provides a 2d array in any dimension combination
! (geo or mag) containing all or part of the global domain, along with
! the field name, long_name, and units. This data is saved in either
! 2d (lon,lat), or 3d (lon,lat,lev) data arrays (fsechist(i)%data) for
! the final history variable at the current model time.
!
! If all 3 spatial dimensions are desired on the history, idx is the
! index of the 3rd dimension (i.e., the dimension not specified by
! dname1 and dname2). In this case, addfld is called from inside a
! loop over the 3rd dimension. If only 2 spatial dimensions are desired
! (lon,lat) then idx can be zero. However, only the levs dimension can
! be excluded (i.e., idx can be zero only if saving a horizontal slice,
! see table below)
!
! This routine saves the data and field information in the derived type
! fsechist(n) (see fields.F). The final variable on the history will be
! either (lon,lat,lev,time), or (lon,lat,time), geo or mag (see subs
! def_fsech and wrfsech in nchist.F).
!
! On Input:
! dname1 Name (char string) of 1st dimension of f (see table below)
! lb1,ub1 Lower and upper integer bounds of 1st dimension of f
! (must be within full domain dimension sizes, see table below)
! dname2 Name (char string) of 2nd dimension of f (see table below)
! lb2,ub2 Lower and upper integer bounds of 2nd dimension of f
! (must be within full domain dimension sizes, see table below)
! f(lb1:ub1,lb2:ub2)
! Full or sub domain data to be saved for secondary history
! idx Index of 3rd dimension (or 0 if writing horizontal slice only)
! (idx can be 0 only if providing f(lon,lat) or f(lat,lon))
! name Name of field (should be in namelist SECFLDS) (<= 16 chars)
! long_name Optional long name of the field (<= 80 chars)
! units Optional units of the field (<= 16 chars)
!
! On Output:
! fsechist(n) is defined (derived type fields_sech, see fields.F). Data
! component has been allocated, and data from f (full or sub-domain)
! has been saved to fsechist(n)%data(lon,lat,lev)
!
! Caveat: If k > 0, midpoints is assumed.
!
! Note 1: nlevp1==nilevp1==nmlevp1==nimlevp1
! Note 2: For geo fields, only i=3:nlonp2 are written to histories (see sub wrfsech)
! Note 3: If the input lon dimension bounds are the full domain dimension
! (e.g., lb1,ub1==1,nlonp4), then it is assumed that only the master
! task is making the call for the current field. In this case (task0_only),
! the data are not gathered to the master task when writing the history
! (see sub mp_gather2root_fsech in mpi.F)
!
! Following is a table of valid combinations for dname1 and dname2, with the
! corresponding full-domain dimension sizes. Final fields on the secondary
! histories will be one of two types: f(lon,lat,lev,time) (when idx > 0),
! or f(lon,lat,time) (when idx==0).
!
! dname1 dname2 idx dim sizes Description
! (full domain) (of input f)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! lev lon j nlevp1,nlonp4 geo lat slice at j (midpoint levs)
! ilev lon j nilevp1,nlonp4 geo lat slice at j (interface levs)
! lon lev j nlonp4,nlevp1 geo lat slice at j (midpoint levs)
! lon ilev j nlonp4,nilevp1 geo lat slice at j (interface levs)
!
! lev lat i nlevp1,nlat geo lon slice at i (midpoint levs)
! ilev lat i nilevp1,nlat geo lon slice at i (interface levs)
! lat lev i nlat,nlevp1 geo lon slice at i (midpoint levs)
! lat ilev i nlat,nilevp1 geo lon slice at i (interface levs)
!
! lon lat k nlonp4,nlat geo horizontal slice at k
! lat lon k nlat,nlonp4 geo horizontal slice at k
! lon lat 0 nlonp4,nlat geo horizontal slice (2d)
! lat lon 0 nlat,nlonp4 geo horizontal slice (2d)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! mlev mlon j nmlevp1,nmlonp1 mag lat slice at j (midpoint levs)
! imlev mlon j nimlevp1,nmlonp1 mag lat slice at j (interface levs)
! mlon mlev j nmlonp1,nmlevp1 mag lat slice at j (midpoint levs)
! mlon imlev j nmlonp1,nimlevp1 mag lat slice at j (interface levs)
!
! mlev mlat i nmlevp1,nmlat mag lon slice at i (midpoint levs)
! imlev mlat i nimlevp1,nmlat mag lon slice at i (interface levs)
! mlat mlev i nmlat,nmlevp1 mag lon slice at i (midpoint levs)
! mlat imlev i nmlat,nimlevp1 mag lon slice at i (interface levs)
!
! mlon mlat k nmlonp1,nmlat mag horizontal slice at k
! mlat mlon k nmlat,nmlonp1 mag horizontal slice at k
! mlon mlat 0 nmlonp1,nmlat mag horizontal slice (2d)
! mlat mlon 0 nmlat,nmlonp1 mag horizontal slice (2d)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
! Example 1: save (lev,lon) array at j for 3d+time field on the history:
! real :: f(lev0:lev1,lon0:lon1) ! field on mpi task subdomains
! do j=lat0,lat1
! do i=lon0,lon1
! do k=lev0,lev1
! f(k,i) = [expression, probably involving k,i,j indices in other arrays]
! enddo
! enddo
! call addfld('MY_3D_DIAG','My long name','my units',
! f,'lev',lev0,lev1,'lon',lon0,lon1,j)
! enddo
!
! Example 2: save (lon,lat) array for 2d+time field on the history:
! real :: f(lon0:lon1,lat0:lat1) ! field on mpi task subdomains
! do j=lat0,lat1
! do i=lon0,lon1
! f(i,j) = [expression, probably involving i,j indices in other arrays]
! enddo
! enddo
! call addfld('MY_2D_DIAG','My long name','my units',
! f,'lon',lon0,lon1,'lat',lat0,lat1,0)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
! Use modules:
use fields_module,only: fsechist,f4d,nf4d,shortname_len
use hist_module,only: isechist,modeltime
use init_module,only: istep
use mpi_module,only: lat0,lat1,lon0,lon1
implicit none
!
! Args:
character(len=*),intent(in) :: dname1,dname2,name,long_name,units
integer,intent(in) :: lb1,ub1,lb2,ub2,idx
real,intent(in) :: f(lb1:ub1,lb2:ub2)
!
! Local:
integer :: ier,i,j,k,ii,ix,ixf4d,iseries_p,iseries_s
integer :: idim1,idim2,idim3,ihdim1,ihdim2,ihdim3
logical :: time2write,wrprim,wrsech,newseries_prim,
| newseries_sech,found
logical :: mag,geo,firstcall
real :: fmin,fmax
character(len=8) :: dname3
integer :: istep1=-1
!
! Valid dimension names on history (dname1 and dname2 must be one of these)
integer,parameter :: nvalid=8
character(len=8) :: valid_dnames(nvalid) =
| (/'lev ','ilev ','lon ','lat ', ! geo dim names
| 'mlev ','imlev ','mlon ','mlat '/) ! mag dim names
!
! Dimensions of global fsechist(ix)%data corresponding to dim names:
integer :: valid_dims(nvalid) =
| (/nlevp1 ,nilevp1 ,nlonp4 ,nlat , ! geo dim sizes
| nmlevp1 ,nimlevp1 ,nmlonp1 ,nmlat /) ! mag dim sizes
!
! External:
integer,external :: strloc
logical,external :: wrhist
!
! Return silently if not writing secondary histories:
if (isechist==0) return
!
! Save only if we are writing a secondary history this step:
time2write = wrhist(istep,modeltime,
| wrprim, newseries_prim, iseries_p,
| wrsech, newseries_sech, iseries_s)
if (.not.wrsech) return
if (istep1 < 0) istep1 = istep
!
! Check for name in fsechist(mxfsech) (namelist SECFLDS)
! If field name was not included in namelist SECFLDS, then the data
! pointer is not allocated, and the field is not defined on secondary
! history. See checks for data allocation (.not.associated(fsechist(ix)%data))
! in nchist.F and mpi.F.
!
ix = strloc(fsechist%short_name,mxfsech,name) ! mxfsech in params.F
if (ix <= 0) then ! field was not requested by SECFLDS in namelist read file
!
! Return silently to avoid cluttering stdout with warnings.
! (often addfld calls are left open in the code, but the field was not
! requested SECFLDS in the namelist read file)
!
return
endif
!
firstcall = .false.
if (.not.associated(fsechist(ix)%data)) firstcall = .true.
!
! Return with warning if field requested is already saved as 4d field
! (the 4d field will be written to the secondary history by nchist module).
! That is, the user added a call to addfld in which the field name (1st arg)
! is already in use as a 4d field (e.g., call addfld('TN'...))
!
ixf4d = strloc(f4d%short_name,nf4d,name)
if (ixf4d > 0) return
!
! Call is for a diagnostic secondary history field (2d or 3d, mag or geo).
! write(6,"('addfld: ix=',i3,' name=',a)") ix,name
!
! Dim types 1 and 2 must be different:
if ((index(dname1,'lev') > 0 .and. index(dname2,'lev') > 0).or.
| (index(dname1,'lon') > 0 .and. index(dname2,'lon') > 0).or.
| (index(dname1,'lat') > 0 .and. index(dname2,'lat') > 0)) then
write(6,"('>>> addfld: redundant dimension types: dname1,2=',
| 2a)") dname1,dname2
call shutdown('addfld')
endif
!
! Set idim1,idim2 (also validate dname1,2):
idim1 = 0 ; idim2 = 0
i = strloc(valid_dnames,nvalid,dname1)
if (i > 0) idim1 = valid_dims(i)
i = strloc(valid_dnames,nvalid,dname2)
if (i > 0) idim2 = valid_dims(i)
if (idim1==0) then
write(6,"(/,'>>> addfld: invalid dname1 = ',a,' for field ',a)")
| dname1,name
write(6,"('Valid dim names: ',8a8)") valid_dnames
call shutdown('addfld')
endif
if (idim2==0) then
write(6,"(/,'>>> addfld: invalid dname2 = ',a,' for field ',a)")
| dname2,name
write(6,"('Valid dim names: ',8a8)") valid_dnames
call shutdown('addfld')
endif
!
! Validate dimension bounds:
if (lb1 < 1 .or. ub1 > idim1) then
write(6,"(/,'>>> addfld bad bounds: dname1=',a,' idim1=',i4,
| ' lb1=',i4,' ub1=',i4)") dname1,idim1,lb1,ub1
call shutdown('addfld')
endif
if (lb2 < 1 .or. ub2 > idim2) then
write(6,"(/,'>>> addfld bad bounds: dname2=',a,' idim2=',i4,
| ' lb2=',i4,' ub2=',i4)") dname2,idim2,lb2,ub2
call shutdown('addfld')
endif
!
! If idx==0, it must refer to lev dimension:
if (idx==0.and.
| ((index(dname1,'lev') > 0.and.index(dname2,'lon') > 0).or.
| (index(dname1,'lev') > 0.and.index(dname2,'lat') > 0).or.
| (index(dname1,'lon') > 0.and.index(dname2,'lev') > 0).or.
| (index(dname1,'lat') > 0.and.index(dname2,'lev') > 0))) then
write(6,"(/,'>>> addfld: idx can be zero only when saving',
| ' f(lon,lat) or f(lat,lon): dname1,2=',2a)") dname1,dname2
call shutdown('addfld')
endif
!
! If longitude dimension is global, assume task0_only, i.e., assume
! field is presented by master task only -- in this case gather2root
! must not be called on this field (sub mp_gather2root_sech in mpi.F)
!
fsechist(ix)%task0_only = .false.
if ((index(dname1,'lon') > 0 .and. lb1==1.and.ub1==idim1) .or.
| (index(dname2,'lon') > 0 .and. lb2==1.and.ub2==idim2))
| fsechist(ix)%task0_only = .true.
!
! Mag or geo:
geo = .true. ; mag = .false.
i = strloc(valid_dnames(5:8),4,dname1)
if (i > 0) then
geo = .false.
mag = .true.
endif
!
! Dims must be all mag or all geo:
! (valid_dnames(1:4) are geo, valid_dnames(5:8) are mag)
!
if (geo) then
i = strloc(valid_dnames(5:8),4,dname1)
ii = strloc(valid_dnames(5:8),4,dname2)
if (i > 0.or.ii > 0) then
write(6,"(/,'>>> addfld: cannot mix geo and mag dimensions:',
| ' dname1,2=',2a)") dname1,dname2
endif
else
i = strloc(valid_dnames(1:4),4,dname1)
ii = strloc(valid_dnames(1:4),4,dname2)
if (i > 0.or.ii > 0) then
write(6,"(/,'>>> addfld: cannot mix geo and mag dimensions:',
| ' dname1,2=',2a)") dname1,dname2
endif
endif
!
! Set full-domain secondary output array dimensions from dname1,2:
! (geographic or magnetic, lev,lat,or lon)
fsechist(ix)%dimnames(:) = ' '
fsechist(ix)%dimsizes(:) = 0
if (index(dname1,'lon') > 0) fsechist(ix)%dimsizes(1) = idim1
if (index(dname1,'lat') > 0) fsechist(ix)%dimsizes(2) = idim1
if (index(dname1,'lev') > 0) fsechist(ix)%dimsizes(3) = idim1
if (index(dname2,'lon') > 0) fsechist(ix)%dimsizes(1) = idim2
if (index(dname2,'lat') > 0) fsechist(ix)%dimsizes(2) = idim2
if (index(dname2,'lev') > 0) fsechist(ix)%dimsizes(3) = idim2
!
if (index(dname1,'lon') > 0) fsechist(ix)%dimnames(1) = dname1
if (index(dname1,'lat') > 0) fsechist(ix)%dimnames(2) = dname1
if (index(dname1,'lev') > 0) fsechist(ix)%dimnames(3) = dname1
if (index(dname2,'lon') > 0) fsechist(ix)%dimnames(1) = dname2
if (index(dname2,'lat') > 0) fsechist(ix)%dimnames(2) = dname2
if (index(dname2,'lev') > 0) fsechist(ix)%dimnames(3) = dname2
!
! Set third dimension of full-domain output array according to idx:
!
! 11/05 btf: If the 3rd dimension is lev, we don't know
! if idx > 0 refers to midpoints or interfaces. This will
! happen only when addfld is called from a k-loop and idx > 0.
!
idim3 = 0 ! 2d array
fsechist(ix)%ndims = 2
! 3d field:
if (idx > 0) then ! 3d array
fsechist(ix)%ndims = 3
if (index(dname1,'lev')==0.and.index(dname2,'lev')==0) then
idim3 = nlevp1 ; if (mag) idim3 = nmlevp1 ! but could be nilevp1 or nimlevp1
elseif (index(dname1,'lat')==0.and.index(dname2,'lat')==0) then
idim3 = nlat ; if (mag) idim3 = nmlat
elseif (index(dname1,'lon')==0.and.index(dname2,'lon')==0) then
idim3 = nlonp4 ; if (mag) idim3 = nmlonp1
endif
if (idim3==0) then
write(6,"('>>> addfld: cannot determine idim3: dname1,2=',a,',
| ',a)") trim(dname1),trim(dname2)
call shutdown('addfld')
endif
dname3 = ' '
do i=1,nvalid
if (idim3==valid_dims(i)) dname3 = valid_dnames(i)
enddo
do i=1,3
if (len_trim(fsechist(ix)%dimnames(i))==0) then
fsechist(ix)%dimsizes(i) = idim3
fsechist(ix)%dimnames(i) = dname3
endif
enddo
if (idx > idim3) then
write(6,"(/,'>>> addfld: field ',a,': bad idx=',i4,
| ' dname3=',a,' idim3=',i4,' (idx should not be > idim3)')")
| trim(name),idx,dname3,idim3
call shutdown('addfld')
endif
endif
!
do i=1,fsechist(ix)%ndims
if (len_trim(fsechist(ix)%dimnames(i))==0.or.
| fsechist(ix)%dimsizes(i)==0) then
write(6,"(/,'>>> addfld: could not determine dimension ',i2,
| ': field ',a,': dimname=',a,' dimsize=',i3)") i,
| trim(name),fsechist(ix)%dimnames(i),fsechist(ix)%dimsizes(i)
call shutdown('addfld')
endif
enddo
!
! write(6,"('addfld: field ',a,' fsechist(ix)%dimnames=',a,',',a,
! | ',',a,' dimsizes=',3i4)") name,trim(fsechist(ix)%dimnames(1)),
! | trim(fsechist(ix)%dimnames(2)),trim(fsechist(ix)%dimnames(3)),
! | fsechist(ix)%dimsizes
!
! Allocate fsechist(ix)%data, either 2d (lon,lat) or 3d (lon,lat,lev):
! (field is mag or geo)
!
! 3d (i,j,k):
if (.not.associated(fsechist(ix)%data)) then
ihdim1=fsechist(ix)%dimsizes(1)
ihdim2=fsechist(ix)%dimsizes(2)
ihdim3=fsechist(ix)%dimsizes(3)
if (idx > 0) then
allocate(fsechist(ix)%data(ihdim1,ihdim2,ihdim3),stat=ier)
if (ier /= 0) then
write(6,"(/,'>>> Error allocating 3d (ihdim1=',i3,
| ',ihdim2=',i3,',ihdim3=',i3,') for field ',a)")
| ihdim1,ihdim2,ihdim3,trim(fsechist(ix)%short_name)
call shutdown('addfld')
else
write(6,"(/,'Allocated 3d sech field ',a,'(',a,'=',i3,
| ',',a,'=',i3,',',a,'=',i3,')')")
| trim(fsechist(ix)%short_name),
| trim(fsechist(ix)%dimnames(1)),ihdim1,
| trim(fsechist(ix)%dimnames(2)),ihdim2,
| trim(fsechist(ix)%dimnames(3)),ihdim3
endif
!
! 2d (i,j):
else ! idx==0 (no lev dimension)
allocate(fsechist(ix)%data(ihdim1,ihdim2,1),stat=ier)
if (ier /= 0) then
write(6,"(/,'>>> Error allocating 2d (ihdim1=',i3,
| ' ihdim2=',i3,' for field ',a)") ihdim1,ihdim2,
| trim(fsechist(ix)%short_name)
call shutdown('addfld')
else
write(6,"(/,'Allocated 2d sech field ',a,'(',a,'=',i3,
| ',',a,'=',i3,')')") trim(fsechist(ix)%short_name),
| trim(fsechist(ix)%dimnames(1)),ihdim1,
| trim(fsechist(ix)%dimnames(2)),ihdim2
endif
endif
fsechist(ix)%data = spval ! init whole array to missing value
endif ! this field data not associated
!
! Assign data to fsech(ix)%data(i,j,k) from f, according to dim names:
!
! f(lev,lon), idx=j
if (index(dname1,'lev') > 0.and.index(dname2,'lon') > 0) then
do k=lb1,ub1
do i=lb2,ub2
fsechist(ix)%data(i,idx,k) = f(k,i)
enddo
! if (trim(fsechist(ix)%short_name)=='RP')
! | write(6,"('addfld: idx=',i3,' k=',i3,' name=',a,
! | ' data(:,idx,k)=',/,(6e12.4))") idx,k,
! | fsechist(ix)%short_name,fsechist(ix)%data(:,idx,k)
enddo
!
! f(lev,lat), idx=i
elseif (index(dname1,'lev') > 0.and.index(dname2,'lat') > 0) then
do k=lb1,ub1
do j=lb2,ub2
fsechist(ix)%data(idx,j,k) = f(k,j)
enddo
enddo
!
! f(lon,lev), idx=j
elseif (index(dname1,'lon') > 0.and.index(dname2,'lev') > 0) then
do k=lb2,ub2
do i=lb1,ub1
fsechist(ix)%data(i,idx,k) = f(i,k)
enddo
enddo
!
! f(lon,lat), idx = k or 0
elseif (index(dname1,'lon') > 0.and.index(dname2,'lat') > 0) then
if (idx==0) then
do j=lb2,ub2
do i=lb1,ub1
fsechist(ix)%data(i,j,1) = f(i,j)
enddo
enddo
else ! 3d
do j=lb2,ub2
do i=lb1,ub1
fsechist(ix)%data(i,j,idx) = f(i,j)
enddo
enddo
endif
!
! f(lat,lev), idx = i
elseif (index(dname1,'lat') > 0.and.index(dname2,'lev') > 0) then
do i=lb2,ub2
do j=lb1,ub1
fsechist(ix)%data(idx,j,k) = f(j,k)
enddo
enddo
!
! f(lat,lon), idx = k or 0
elseif (index(dname1,'lat') > 0.and.index(dname2,'lon') > 0) then
if (idx==0) then
do j=lb1,ub1
do i=lb2,ub2
fsechist(ix)%data(i,j,1) = f(j,i)
enddo
enddo
else
do j=lb1,ub1
do i=lb2,ub2
fsechist(ix)%data(i,j,idx) = f(j,i)
enddo
enddo
endif
endif
!
! Set other structure components:
fsechist(ix)%long_name = long_name
fsechist(ix)%units = units
fsechist(ix)%geo = geo
fsechist(ix)%mag = mag
!
! Report to stdout:
if (firstcall) then
write(6,"(/,'Initialized diagnostic secondary history field ',
| a,' (ix=',i3,'):')") trim(fsechist(ix)%short_name),ix
write(6,"(' short_name = ',a)") trim(fsechist(ix)%short_name)
write(6,"(' long_name = ',a)") trim(fsechist(ix)%long_name)
write(6,"(' units = ',a)") trim(fsechist(ix)%units)
write(6,"(' geo = ',l1)") fsechist(ix)%geo
write(6,"(' mag = ',l1)") fsechist(ix)%mag
write(6,"(' dimnames = ',3a)") fsechist(ix)%dimnames
write(6,"(' dimsizes = ',3i4)") fsechist(ix)%dimsizes
write(6,"(' ndims = ', i2)") fsechist(ix)%ndims
write(6,"(' task0_only = ', l5)") fsechist(ix)%task0_only
endif ! first call for this field
!
! Print non-spval min,max, i.e., portion contributed by current task:
! if (fsechist(ix)%ndims==3) then
! call fminmaxspv(fsechist(ix)%data(:,:,:),
! | fsechist(ix)%dimsizes(1)*fsechist(ix)%dimsizes(2)*
! | fsechist(ix)%dimsizes(3),fmin,fmax,spval)
! write(6,"(' 3-d min,max (non-spval portion so far)= ',
! | 2e12.4)") fmin,fmax
! else
! call fminmaxspv(fsechist(ix)%data(:,:,1),
! | fsechist(ix)%dimsizes(1)*fsechist(ix)%dimsizes(2),
! | fmin,fmax,spval)
! write(6,"(' 2-d min,max (non-spval portion so far)= ',
! | 2e12.4)") fmin,fmax
! endif
end subroutine addfld
!-----------------------------------------------------------------------
end module addfld_module