mo_grib.f90
MODULE mo_grib
! Descritpion:
!
5: ! This module contains the output routines based on GRIB edition 1.
!
! For a detailed description of the contained subroutines look on their
! header.
!
10: ! Most important new features:
! - Use of a code table in section 1 which allows the use of 128 tables
! each containg 128 different variables.
! - Uses the sub center entry (255 for ECHAM). Center is still 98 for
! ECMWF
15: ! - The century parameter is correctly used. Now, in 360 day mode,
! 25599 years can be simulated before an overflow in the date occures.
! - Uses correct dates in 365/366 day mode. Allows for clean forecast and
! nudging runs.
! - Either the intrinsic codegb5 C function can be used for encoding, or
20: ! the ECMWF EMOS library. For this define for compile -DEMOS for CFLAGS
! and FFLAGS.
! - The output is now real standard !!!
!
! Authors:
25: !
! L. Kornblueh, MPI, December 1998
! L. Kornblueh, MPI, April 1998, added NWP forecast mode
!
30: IMPLICIT NONE
INTEGER :: kleng ! size of kgrib
INTEGER :: klengo ! byte size of kgrib after codegb5
INTEGER, ALLOCATABLE :: kgrib(:) ! grib data set
35:
INTEGER :: ksec0(2), ksec1(43)
INTEGER :: ksec2_gp(22), ksec2_sp(22)
INTEGER :: ksec3(2), ksec4(42)
40: REAL, ALLOCATABLE :: psec2(:), psec3(:), psec4(:)
INTEGER :: kword ! output size in words
INTEGER :: klenp ! size of psec4
45:
INTEGER, PARAMETER :: nint = 0
INTEGER, PARAMETER :: nbit = BIT_SIZE(nint)
INTEGER, PARAMETER :: iobyte = nbit/8
50: INTEGER, PARAMETER :: local_table = 128 ! local code table
INTEGER, PARAMETER :: nudging_table = 129 ! nudging code table
INTEGER, PARAMETER :: g4x_table = 130 ! G4X code table
INTEGER, PARAMETER :: center_id = 98 ! identification of centre
INTEGER :: model_id ! model identification
55: INTEGER, PARAMETER :: grid_type = 255 ! grid definition
INTEGER, PARAMETER :: nflag = 128 ! flag(see code table 1)
INTEGER :: code_parameter ! data field code
! (see code table 2 )
INTEGER :: level_type ! indicator of type of level
60: ! (see code table 3)
INTEGER :: level_p1 ! data level1 (see code table 3)
INTEGER :: level_p2 ! data level2 (see code table 3)
! reference time of data
65:
INTEGER :: year ! year of century
INTEGER :: month ! month
INTEGER :: day ! day
INTEGER :: hour ! hour
70: INTEGER :: minute ! minute
INTEGER :: time_unit = 0 ! unit of time range
! (see code table 4)
INTEGER :: time_p1 = 0 ! time range 1
75: INTEGER :: time_p2 = 0 ! time range 2
INTEGER :: range_flag = 10 ! time range flag
! (see code table 5)
INTEGER :: century ! century
INTEGER, PARAMETER :: subcenter_id = 255 ! subcenter
80: INTEGER, PARAMETER :: decimal_scale = 0 ! decimal scale
INTEGER, PARAMETER :: local_extension = 0 ! local extension flag
INTEGER :: forecast_hours = 0 ! number of forecast hours in
! NWP mode
85:
! Variables for block 2
INTEGER :: npvct ! no. of vertical ccordinate parameters
90: ! type gaussian latitudes ( data_type = 4)
INTEGER, PARAMETER :: ngptype = 4 ! data representation type (code table 6)
INTEGER :: npalat ! number of points along a latitude
INTEGER :: npamer ! number of points along a meridional
95: INTEGER :: nglaor ! latitude of origin in degree*1000
INTEGER, PARAMETER :: ngloor = 0 ! longitude of origin in degree*1000
INTEGER, PARAMETER :: neinfl = 128 ! resolution flag(code table 7)
INTEGER :: nglaex ! latitude of extreme point degree*1000
INTEGER :: ngloex ! longitude of extreme point degree*1000
100: INTEGER :: ngdi ! latitude increment in degree*1000
INTEGER :: nglpe ! numder of latitude lines between a pole
! and the equator
INTEGER, PARAMETER :: nscmfl = 0 ! scanning mode(flags see table 8)
105: ! type ( spherical harmonics data_type = 50)
INTEGER, PARAMETER :: nsptype = 50 ! data representation type (code table 6)
INTEGER :: nfasmj ! j - pentagonal resolution parameter
INTEGER :: nfasmk ! k - pentagonal resolution parameter
110: INTEGER :: nfasmm ! m - pentagonal resolution parameter
INTEGER, PARAMETER :: nrtsh = 1 ! representation type (see table 9)
INTEGER, PARAMETER :: nrmsh = 1 ! representation mode (see table 10)
CONTAINS
115:
SUBROUTINE set_output_time
USE mo_exception
USE mo_start_dataset, ONLY : nstep
120: USE mo_control, ONLY : ncbase, ntbase, dtime, lnwp
USE mo_year, ONLY : ic2ymd
IMPLICIT NONE
125: INTEGER, SAVE :: century0, year0, month0, day0, hour0, minute0
! Local scalars:
REAL :: zsecs
130: INTEGER :: itbase, date
LOGICAL, SAVE :: l_nwp_started = .FALSE.
! Intrinsic functions
135: INTRINSIC MOD
IF (.NOT. l_nwp_started .OR. .NOT. lnwp ) THEN
zsecs = (nstep+1)*dtime + ntbase
date = ic2ymd(ncbase+INT(zsecs/86400))
140: itbase = MOD(zsecs,86400.)
year = date/10000
month = MOD(date/100,100)
day = MOD(date,100)
hour = itbase/3600
145: minute = MOD(itbase/60,60)
century = year/100+1
year = MOD(year,100)
150: IF (century > 255 .OR. century < 0) THEN
CALL finish ('set_output_time', &
& 'This century is not supported in GRIB 1')
END IF
END IF
155:
IF (lnwp .AND. .NOT. l_nwp_started) THEN
century0 = century
year0 = year
month0 = month
160: day0 = day
hour0 = hour
minute0 = minute-ANINT(dtime/60) ! adjust for start of NWP forecast
ENDIF
165: IF (lnwp) THEN ! NWP mode ..., not analysis
l_nwp_started = .TRUE.
range_flag = 0
time_unit = 1
century = century0
170: year = year0
month = month0
day = day0
hour = hour0
minute = minute0
175: ENDIF
IF (l_nwp_started) THEN
forecast_hours = ANINT((nstep+1)*dtime/3600)
END IF
180:
END SUBROUTINE set_output_time
SUBROUTINE close_grib_file
185: USE mo_post, ONLY: nunitdf
IMPLICIT NONE
INTEGER :: iret
190:
! External subroutines
EXTERNAL pbclose
CALL pbclose(nunitdf,iret)
195:
END SUBROUTINE close_grib_file
SUBROUTINE open_grib_file
200: USE mo_doctor, ONLY: nout, nerr
USE mo_exception, ONLY: finish
USE mo_post, ONLY: nunitdf
USE mo_filename, ONLY: path_limit, standard_grib_file
USE mo_control, ONLY: lnwp
205:
IMPLICIT NONE
INTEGER :: iret
210: CHARACTER (3) :: fcst
CHARACTER (path_limit+5) :: filename
! External subroutines
EXTERNAL pbopen
215:
IF (lnwp) THEN
CALL set_output_time
IF (forecast_hours > 744) THEN
WRITE (nerr,*) 'NWP mode makes no sense for this time range.'
220: WRITE (nerr,*) 'Please change to the climate mode.'
CALL finish('open_grib_file','Run terminated.')
END IF
WRITE (fcst,'(i3.3)') forecast_hours
filename = TRIM(standard_grib_file)//'+'//fcst
225: WRITE (nout,*) 'NWP forecast (', forecast_hours, &
'hr) GRIB output: ', &
TRIM(filename(INDEX(filename,'/',.TRUE.)+1:))
ELSE
filename = TRIM(standard_grib_file)
230: END IF
CALL pbopen(nunitdf,filename,'w',iret)
IF (iret /= 0) THEN
WRITE (nerr,*) 'Could not open file: ', filename
235: CALL finish('open_grib_file','Run terminated.')
END IF
END SUBROUTINE open_grib_file
240: SUBROUTINE init_grib
USE mo_constants, ONLY : api
USE mo_control, ONLY : nvclev, nlon, ngl, nhgl, nn, nm, nk, vct, lanalysis
USE mo_gaussgrid, ONLY : gmu
245: USE mo_doctor, ONLY : nvers
USE mo_post, ONLY : nlalo
IMPLICIT NONE
250: ! allocate grib work array
nlalo = nlon*ngl
kleng = nlalo + 1000 ! only for max 32 (64)bit packing
! dependend on generic INTEGER size,
255: ! nlalo - packed field size
! 1000 - upper limit for header size
IF (.NOT. ALLOCATED(kgrib)) ALLOCATE(kgrib(kleng))
! set all GRIB block entries to 0
260:
ksec0(:) = 0
ksec1(:) = 0
ksec2_gp(:) = 0
ksec2_sp(:) = 0
265: ksec3(:) = 0
ksec4(:) = 0
! set model version and reset century
270: model_id = nvers
century = 0
level_p1 = 0
! set fixed values in GRIB block 1, common to all GRIB sets
275:
! to set standard output to analysis mode lnwp must be .false. and
! lanalysis .true.. Later with working adiabatic NMI this must be
! changed to range_flag = 1. range_flag = 0 is default for output
! of the OI or 3DVAR analysis, this are given external. In case a
280: ! 4DVAR is running inside ECHAM this value must be adjusted to
! range_flag = 0, as it is now.
! For usual climate mode runs range_flag is set to 10, which means
! The time given in the GRIB header in section 1 is the valid time
! and time_p1 and time_p2 do not mean anything.
285:
IF (lanalysis) range_flag = 0
ksec1(2) = center_id
ksec1(3) = model_id
290: ksec1(4) = grid_type ! No WMO predefined
ksec1(5) = nflag
ksec1(8) = level_p1
ksec1(15) = time_unit
ksec1(16) = time_p1
295: ksec1(17) = time_p2
ksec1(18) = range_flag
ksec1(21) = century ! preliminary value
ksec1(22) = subcenter_id
ksec1(23) = decimal_scale
300: ksec1(24) = local_extension ! Until now no extension, might be used
! for paleoclimate runs centuries
npalat = nlon
npamer = ngl
305: nglaor = ASIN(gmu(1))*180000.0/api
nglaex = -nglaor
ngdi = 360000/npalat
ngloex = 360000-ngdi
nglpe = nhgl
310:
ksec2_gp(1) = ngptype
ksec2_gp(2) = npalat
ksec2_gp(3) = npamer
ksec2_gp(4) = nglaor
315: ksec2_gp(5) = ngloor
ksec2_gp(6) = neinfl
ksec2_gp(7) = nglaex
ksec2_gp(8) = ngloex
ksec2_gp(9) = ngdi
320: ksec2_gp(10) = nglpe
ksec2_gp(11) = nscmfl
ksec2_gp(12) = 2*nvclev
nfasmj = nn
325: nfasmk = nk
nfasmm = nm
ksec2_sp(1) = nsptype
ksec2_sp(2) = nfasmj
330: ksec2_sp(3) = nfasmk
ksec2_sp(4) = nfasmm
ksec2_sp(5) = nrtsh
ksec2_sp(6) = nrmsh
ksec2_sp(12) = 2*nvclev
335:
ksec3(1:2) = 0
ksec4(1:42) = 0
340: ksec4(2) = 24 ! bits, predefined standard in ECHAM4
#ifdef EMOS
ALLOCATE (psec2(10+2*nvclev))
ALLOCATE (psec3(0))
345: psec2(1:10) = 0.0
psec2(11:10+2*nvclev) = vct(1:2*nvclev)
! ECMWF GRIB library does not except centuries < 20,
! so checking has to be disabled (ly365 seems to be
350: ! correct with checking)
CALL grsvck (0)
#endif
355:
END SUBROUTINE init_grib
SUBROUTINE outsp
360: ! Description:
!
! Control postprocessing of spectral fields
!
! Method:
365: !
! Pack and write
!
! *outsp* is called from subroutine *stepon*
!
370: ! Results:
! Packed spectral fields are written out
!
! Externals:
! *codegb5* - pack to grib code
375: !
! Authors:
!
! U. Schlese, DKRZ, July 1994, original source
! U. Schulzweida, MPI, May 1998, f90 rewrite
380: ! L. Kornblueh, MPI, November 1998, f90 rewrite, extended to GRIB 1
!
! for more details see file AUTHORS
!
385: USE mo_exception
USE mo_memory_sp
USE mo_parameters
USE mo_start_dataset
USE mo_control
390: USE mo_post
USE mo_year
USE mo_doctor
USE mo_mpi, ONLY: p_pe, p_io
USE mo_decomposition, ONLY: dcg => global_decomposition
395: USE mo_transpose, ONLY: gather_sp
IMPLICIT NONE
! Local scalars:
400: INTEGER :: ierr, iret, jlev
! Local arrays:
REAL, POINTER :: ptr3d(:,:,:)
REAL, POINTER :: z3d(:,:,:)
405: REAL :: psec4(2,nsp)
TYPE (memory_info) :: info
! External subroutines
410: EXTERNAL codegb5, pbwrite
! Executable statements
NULLIFY (z3d)
415:
! Finish definition GRIB block 1
ksec1(1) = local_table
420: CALL set_output_time
IF (lnwp) THEN
ksec1(15) = time_unit
ksec1(16) = forecast_hours
425: ksec1(18) = range_flag
END IF
ksec1(10) = year
ksec1(11) = month
ksec1(12) = day
430: ksec1(13) = hour
ksec1(14) = minute
ksec1(21) = century
!#ifdef DEBUG
435: IF (lnwp) THEN
WRITE (nout,'(a,i2.2,a1,i2.2,2x,i2.2,a1,i2.2,a1,i4.4,a,i3.3,a)') &
& ' Output time: ', hour, ':', minute, &
& day, '.', month, '.', year+100*(century-1), &
& ', forecast at +', forecast_hours, ' hr'
440: ELSE
WRITE (nout,'(a,i2.2,a1,i2.2,2x,i2.2,a1,i2.2,a1,i4.4)') &
& ' Output time: ', hour, ':', minute, &
& day, '.', month, '.', year+100*(century-1)
END IF
445: !#endif
IF (lppp .OR. lppt) THEN
CALL get_entry (sp, 'STP', ptr3d)
450: CALL get_info (sp, 'STP', info)
ALLOCATE (z3d(info%gdim_1, info%gdim_2, info%gdim_3))
CALL gather_sp (z3d, ptr3d, dcg)
END IF
455:
! Pressure
IF (lppp) THEN
460: code_parameter = 152
level_type = 1
level_p2 = 0
ksec1(6) = code_parameter
465: ksec1(7) = level_type
ksec1(9) = level_p2
IF (p_pe == p_io) THEN
psec4(:,:) = z3d(nlevp1,:,:)
470: #ifdef EMOS
ksec4(1) = n2sp
CALL gribex (ksec0, ksec1, ksec2_sp, psec2, ksec3, psec3, &
ksec4, psec4, n2sp, kgrib, kleng, kword, 'C', ierr)
#else
475: CALL codegb5 (psec4,n2sp,16,nbit,ksec1,ksec2_sp,vct,2*nvclev, &
kgrib,klengo,kword,0,ierr)
#endif
CALL pbwrite(nunitdf,kgrib,kword*iobyte,iret)
480: IF (iret /= kword*iobyte) CALL errmsg
END IF
END IF
! Temperature
485:
IF (lppt) THEN
code_parameter = 130
level_type = 109
490:
ksec1(6) = code_parameter
ksec1(7) = level_type
DO jlev = 1, nlev
495:
level_p2 = jlev
ksec1(9) = level_p2
IF (p_pe == p_io) THEN
500: psec4(:,:) = z3d(jlev,:,:)
#ifdef EMOS
ksec4(1) = n2sp
CALL gribex (ksec0, ksec1, ksec2_sp, psec2, ksec3, psec3, &
ksec4, psec4, n2sp, kgrib, kleng, kword, 'C', ierr)
505:
#else
CALL codegb5(psec4,n2sp,16,nbit,ksec1,ksec2_sp,vct,2*nvclev, &
kgrib,klengo,kword,0,ierr)
#endif
510: CALL pbwrite(nunitdf,kgrib,kword*iobyte,iret)
IF (iret /= kword*iobyte) CALL errmsg
END IF
END DO
515: END IF
IF (ASSOCIATED(z3d)) DEALLOCATE(z3d)
! Divergence
520:
IF (lppd) THEN
CALL get_entry (sp, 'SD', ptr3d)
CALL get_info (sp, 'SD', info)
525: ALLOCATE (z3d(info%gdim_1, info%gdim_2, info%gdim_3))
CALL gather_sp (z3d, ptr3d, dcg)
code_parameter = 155
level_type = 109
530:
ksec1(6) = code_parameter
ksec1(7) = level_type
DO jlev = 1, nlev
535:
level_p2 = jlev
ksec1(9) = level_p2
IF (p_pe == p_io) THEN
540: psec4(:,:) = z3d(jlev,:,:)
#ifdef EMOS
ksec4(1) = n2sp
CALL gribex (ksec0, ksec1, ksec2_sp, psec2, ksec3, psec3, &
ksec4, psec4, n2sp, kgrib, kleng, kword, 'C', ierr)
545: #else
CALL codegb5(psec4,n2sp,16,nbit,ksec1,ksec2_sp,vct,2*nvclev, &
kgrib,klengo,kword,0,ierr)
#endif
CALL pbwrite(nunitdf,kgrib,kword*iobyte,iret)
550:
IF (iret /= kword*iobyte) CALL errmsg
END IF
END DO
END IF
555:
IF (ASSOCIATED(z3d)) DEALLOCATE(z3d)
! Vorticity
560: IF (lppvo) THEN
CALL get_entry (sp, 'SVO', ptr3d)
CALL get_info (sp, 'SVO', info)
ALLOCATE (z3d(info%gdim_1, info%gdim_2, info%gdim_3))
565: CALL gather_sp (z3d, ptr3d, dcg)
code_parameter = 138
level_type = 109
570: ksec1(6) = code_parameter
ksec1(7) = level_type
DO jlev = 1, nlev
575: level_p2 = jlev
ksec1(9) = level_p2
IF (p_pe == p_io) THEN
psec4(:,:) = z3d(jlev,:,:)
580: #ifdef EMOS
ksec4(1) = n2sp
CALL gribex (ksec0, ksec1, ksec2_sp, psec2, ksec3, psec3, &
ksec4, psec4, n2sp, kgrib, kleng, kword, 'C', ierr)
585: #else
CALL codegb5(psec4,n2sp,16,nbit,ksec1,ksec2_sp,vct,2*nvclev, &
kgrib,klengo,kword,0,ierr)
#endif
CALL pbwrite(nunitdf,kgrib,kword*iobyte,iret)
590:
IF (iret /= kword*iobyte) CALL errmsg
END IF
END DO
END IF
595:
IF (ASSOCIATED(z3d)) DEALLOCATE(z3d)
RETURN
600: CONTAINS
SUBROUTINE errmsg
CALL finish('outsp', &
'I/O error on spectral harmonics output - disk full?')
605: END SUBROUTINE errmsg
END SUBROUTINE outsp
SUBROUTINE outgpi
610:
! Description:
!
! Control postprocessing of gridpoint fields.
!
615: ! Method:
!
! 2-dimensional fields are allocated
! all values are collected in these field sorted north to south
! at completion packing is done by codegb5.
620: !
! *outgpi* is called from subroutine *scan1sl*
!
! Results:
! The selected output fields are written out.
625: !
! Externals:
! *codegb5* - pack grid point field to grib code
!
! Authors:
630: !
! U. Schlese, DKRZ, July 1994, original source
! L. Kornblueh, MPI, May 1998, f90 rewrite
! U. Schulzweida, MPI, May 1998, f90 rewrite
!
635: ! for more details see file AUTHORS
!
USE mo_exception
USE mo_parameters
640: USE mo_start_dataset
USE mo_control
USE mo_post
USE mo_year
USE mo_doctor
645: USE mo_memory_g3b
USE mo_mpi, ONLY: p_pe, p_io
USE mo_decomposition, ONLY: dcg => global_decomposition
USE mo_transpose, ONLY: gather_gp
650: IMPLICIT NONE
! Local scalars:
REAL :: quot
INTEGER :: ierr, ilev, ipbits, iret, jfield, jlev
655: CHARACTER (8) :: yname
! Local arrays:
! pointers for grid point space
660: REAL, POINTER :: ptr3d(:,:,:)
REAL, POINTER :: z3d(:,:,:)
REAL :: psec4(nlon,ngl)
TYPE (memory_info) :: info
665:
! External subroutines
EXTERNAL codegb5, pbwrite
! Intrinsic functions
670: INTRINSIC MOD
! Executable statements
675: ! Definition grib blocks
quot = 1.0/(nptime*dtime)
! Define grib block 1
680:
ksec1(1) = local_table
CALL set_output_time
685: IF (lnwp) THEN
ksec1(15) = time_unit
ksec1(16) = forecast_hours
ksec1(18) = range_flag
END IF
690: ksec1(10) = year
ksec1(11) = month
ksec1(12) = day
ksec1(13) = hour
ksec1(14) = minute
695: ksec1(21) = century
! Loop over all fields
DO jfield = ncdmin, ncdmax
700:
yname = yn(jfield)
IF (yname(1:1)==' ') CYCLE
IF (yname(1:3)=='G4X') CYCLE
705: ipbits = npbits(jfield)
code_parameter = jfield
ksec1(6) = code_parameter
710:
CALL get_entry (g3b, yname, ptr3d)
CALL get_info (g3b, yname, info)
ALLOCATE (z3d(info%gdim_1, info%gdim_2, info%gdim_3))
call gather_gp (z3d, ptr3d, dcg)
715:
IF (p_pe == p_io) THEN
DO jlev = 1, nlevp1
ilev = npplev(jlev,jfield)
720:
IF (ilev == 0) EXIT
! Find first puzzle piece
725: IF (ilev == -100) THEN
IF (SIZE(z3d,3) == 1) THEN
psec4(1:nlon,:) = z3d(1:nlon,:,1)
ELSE IF (SIZE(z3d,2) == 1) THEN
psec4(1:nlon,:) = z3d(1:nlon,1,:)
730: END IF
ELSE
psec4(1:nlon,:) = z3d(1:nlon,ilev,:)
END IF
IF (laccu(jfield)) psec4(:,:) = psec4(:,:)*quot
735:
! Define variable parts of grib block 1
! Due to the limited layer information possibilities in
! GRIB 1 the soil depth is given as the height of the
740: ! layer center (level_type 111).
IF (ilev == -100) THEN
IF (code_parameter == 167 .OR. & ! 2 m temperature
code_parameter == 168 .OR. & ! 2 m dew point
745: code_parameter == 201 .OR. & ! 2 m maximum temperature
code_parameter == 202 .OR. & ! 2 m minimum temperature
code_parameter == 216) THEN ! 2 m maximum wind speed
level_type = 105
level_p2 = 2
750: ELSE IF (code_parameter == 165 .OR. & ! 10 m u
code_parameter == 166 .OR. & ! 10 m v
code_parameter == 171) THEN ! 10 m wind speed
level_type = 105
level_p2 = 10
755: ELSE IF (code_parameter == 207) THEN ! TD3
level_type = 111
level_p2 = 3
ELSE IF (code_parameter == 208) THEN ! TD4
level_type = 111
760: level_p2 = 19
ELSE IF (code_parameter == 209) THEN ! TD5
level_type = 111
level_p2 = 78
ELSE IF (code_parameter == 170) THEN ! TD
765: level_type = 111
level_p1 = 268/256
level_p2 = MOD(268, 256)
ELSE IF (code_parameter == 183) THEN ! TDCL
level_type = 111
770: level_p1 = 698/256
level_p2 = MOD(698, 256)
ELSE
level_type = 1
level_p2 = 0
775: ENDIF
ELSE
level_type = 109
level_p2 = ilev
END IF
780:
ksec1(7) = level_type
ksec1(8) = level_p1
ksec1(9) = level_p2
785: ! Pack and write
#ifdef EMOS
ksec4(1) = nlalo
ksec4(2) = ipbits
CALL gribex (ksec0, ksec1, ksec2_gp, psec2, ksec3, psec3, &
790: ksec4, psec4, nlalo, kgrib, kleng, kword, 'C', ierr)
#else
CALL codegb5(psec4,nlalo,ipbits,nbit,ksec1,ksec2_gp,vct, &
2*nvclev, kgrib,klengo,kword,0,ierr)
795: #endif
CALL pbwrite(nunitdf,kgrib,kword*iobyte,iret)
IF (iret /= kword*iobyte) THEN
CALL finish('outgpi','I/O error - disks full?')
800: END IF
level_p1 = 0
ksec1(8) = level_p1
805: END DO
END IF
DEALLOCATE (z3d)
810: END DO
END SUBROUTINE outgpi
SUBROUTINE outgpli
815:
! Description:
!
! Controls postprocessing of semi lagrangian gridpoint fields.
!
820: ! Method:
!
! 2-dimensional fields are allocated
! all values are collected in these field sorted north to south
! at completion packing is done by codegb
825: !
! *outgpli* is called from subroutine *scan1sl*
!
! Results:
! The selected output fields are written out
830: !
! Externals:
! *codegb5* - pack grid point field to grib code
!
! Authors:
835: !
! U. Schlese, DKRZ, July 1994, original source
! L. Kornblueh, MPI, May 1998, f90 rewrite
! U. Schulzweida, MPI, May 1998, f90 rewrite
!
840: ! for more details see file AUTHORS
!
USE mo_exception
USE mo_parameters
845: USE mo_start_dataset
USE mo_control
USE mo_post
USE mo_tracer
USE mo_year
850: USE mo_memory_gl
USE mo_mpi, ONLY: p_pe, p_io
USE mo_decomposition, ONLY: dcg => global_decomposition
USE mo_transpose, ONLY: gather_gp
855: IMPLICIT NONE
! Local scalars:
INTEGER :: ierr, iret, ivar, ivbase, jfield, jlev, jt
CHARACTER (8) :: yname
860:
! Local arrays:
! pointers for grid point space
REAL, POINTER :: ptr3d(:,:,:)
865: REAL, POINTER :: z3d(:,:,:)
REAL :: psec4(nlon,ngl)
TYPE (memory_info) :: info
870: INTEGER :: icodem(jptrac+jps), idispv(jptrac+jps), illen(jptrac+jps)
! External subroutines
EXTERNAL codegb5, pbwrite
875: ! Executable statements
! Create tables
ivar = 0
880:
! q-field
IF (lppq) THEN
ivar = ivar + 1
885: icodem(ivar) = 133
idispv(ivar) = 0
illen(ivar) = nlp2
END IF
890: ! x-field
IF (lppx) THEN
ivar = ivar + 1
icodem(ivar) = 153
895: idispv(ivar) = nlp2*nlev
illen(ivar) = nlp2
END IF
! Additional tracers
900:
DO jt = 1, ntrac
IF (lppxt(jt)) THEN
ivar = ivar + 1
icodem(ivar) = ntcode(jt)
905: idispv(ivar) = 2*nlp2*nlev + (jt-1)*nlon*nlev
illen(ivar) = nlon
END IF
END DO
910: ! Define grib blocks
! Define grib block 1
ksec1(1) = local_table
915:
level_type = 109
ksec1(7) = level_type
CALL set_output_time
920:
IF (lnwp) THEN
ksec1(15) = time_unit
ksec1(16) = forecast_hours
ksec1(18) = range_flag
925: END IF
ksec1(10) = year
ksec1(11) = month
ksec1(12) = day
ksec1(13) = hour
930: ksec1(14) = minute
ksec1(21) = century
! Loop over all fields
935: DO jfield = 1, ivar
ivbase = idispv(jfield) + 1
code_parameter = icodem(jfield)
ksec1(6) = code_parameter
940:
IF (jfield == 1) yname = 'Q'
IF (jfield == 2) yname = 'X'
IF (jfield >= 3) yname = 'XT'
945: CALL get_entry (gl, yname, ptr3d)
CALL get_info (gl, yname, info)
ALLOCATE (z3d(info%gdim_1, info%gdim_2, info%gdim_3))
call gather_gp (z3d, ptr3d, dcg)
950: IF (p_pe == p_io) THEN
DO jlev = 1, nlev
level_p2 = jlev
955: ksec1(9) = level_p2
psec4(1:nlon,:) = z3d(1:nlon,jlev,:)
! Pack and write
960: #ifdef EMOS
ksec4(1) = nlalo
CALL gribex (ksec0, ksec1, ksec2_gp, psec2, ksec3, psec3, &
ksec4, psec4, nlalo, kgrib, kleng, kword, 'C', ierr)
#else
965: CALL codegb5(psec4,nlalo,16,nbit,ksec1,ksec2_gp,vct,2*nvclev, &
kgrib,klengo,kword,0,ierr)
#endif
CALL pbwrite(nunitdf,kgrib,kword*iobyte,iret)
970: IF (iret /= kword*iobyte) THEN
CALL finish('outgpli','I/O error - disks full?')
END IF
END DO
975: END IF
DEALLOCATE (z3d)
END DO
END SUBROUTINE outgpli
980:
SUBROUTINE outgpx
! Description:
!
985: ! Control postprocessing of g4x-gridpoint fields.
!
! Method:
!
! 2-dimensional fields are allocated
990: ! all values are collected in these field sorted north to south
! at completion packing is done by codegb
!
! *outgpx* is called from subroutine *scan1sl*
!
995: ! Results:
! The selected output fields are written out
!
! Externals:
! *codegb* - pack grid point field to grib code
1000: !
! Authors:
!
! U. Schlese, DKRZ, July 1994, original source
! L. Kornblueh, MPI, May 1998, f90 rewrite
1005: ! U. Schulzweida, MPI, May 1998, f90 rewrite
!
! for more details see file AUTHORS
!
1010: USE mo_exception
USE mo_parameters
USE mo_start_dataset
USE mo_control
USE mo_post
1015: USE mo_year
!!$ IMPLICIT NONE
!!$
!!$ ! Local scalars:
1020: !!$ REAL :: quot
!!$ INTEGER :: iend, ierr, ig4base, ilev, &
!!$& ipbits, iret, irow, istart, ix, jfield, jlev
!!$ CHARACTER (8) :: yname
!!$
1025: !!$ ! Local arrays:
!!$ ! pointers for grid point space
!!$ REAL, POINTER :: z1dim(:)
!!$ REAL :: psec4(nlon,ngl)
!!$
1030: !!$ ! External subroutines
!!$ EXTERNAL codegb5, pbwrite
!!$
!!$ ! Intrinsic functions
!!$ INTRINSIC INT, MOD
1035: !!$
!!$
!!$ ! Executable statements
!!$
!!$ ! Define grib block 1
1040: !!$
!!$ ksec1(1) = g4x_table
!!$
!!$ CALL set_output_time
!!$
1045: !!$ IF (lnwp) THEN
!!$ ksec1(15) = time_unit
!!$ ksec1(16) = forecast_hours
!!$ ksec1(18) = range_flag
!!$ END IF
1050: !!$ ksec1(10) = year
!!$ ksec1(11) = month
!!$ ksec1(12) = day
!!$ ksec1(13) = hour
!!$ ksec1(14) = minute
1055: !!$ ksec1(21) = century
!!$
!!$ ! Loop over all g4x-fields
!!$
!!$ DO jfield = 1, 128
1060: !!$ yname = yn(jfield)
!!$ IF (yname(1:1)==' ') CYCLE
!!$ IF (yname(1:3)/='G4X') CYCLE
!!$ ipbits = npbits(jfield)
!!$
1065: !!$ DO jlev = 1, nlevp1
!!$ ilev = npplev(jlev,jfield)
!!$ IF (ilev==0) EXIT
!!$
!!$ ! Find first puzzle piece
1070: !!$
!!$ ig4base = m_offset_g4b
!!$ istart = ig4base
!!$ IF (ilev==0) THEN
!!$ istart = istart + (ilev-1)*nlon
1075: !!$ END IF
!!$ iend = istart + nlon - 1
!!$ z1dim => mptglo(ng1a) %vp(istart:iend)
!!$
!!$ ! Collect data
1080: !!$
!!$ DO irow = 1, ngl
!!$ IF (MOD(irow,2)==1) THEN
!!$ ix = (irow+1)/2
!!$ ELSE
1085: !!$ ix = ngl + 1 - irow/2
!!$ END IF
!!$
!!$ IF (laccu(jfield)) THEN
!!$ quot = 1.0/(n4ptime*dtime)
1090: !!$ psec4(1:nlon,ix) = z1dim(1:nlon)*quot
!!$ ELSE
!!$ psec4(1:nlon,ix) = z1dim(1:nlon)
!!$ END IF
!!$ istart = istart + nlon
1095: !!$ iend = iend + nlon - 1
!!$ z1dim => mptglo(ng1a) %vp(istart:iend)
!!$ END DO
!!$
!!$ ! Pack and write
1100: !!$
!!$ code_parameter = jfield
!!$ IF (ilev==-100) THEN
!!$ level_type = 1
!!$ level_p2 = 0
1105: !!$ ELSE
!!$ level_type = 109
!!$ level_p2 = ilev
!!$ END IF
!!$
1110: !!$ ksec1(7) = level_type
!!$ ksec1(9) = level_p2
!!$
!!$#ifdef EMOS
!!$ ksec4(1) = nlalo
1115: !!$ CALL gribex (ksec0, ksec1, ksec2_gp, psec2, ksec3, psec3, &
!!$& ksec4, psec4, nlalo, kgrib, kleng, kword, 'C', ierr)
!!$
!!$#else
!!$ CALL codegb5(psec4,nlalo,ipbits,nbit,ksec1,ksec2_gp,vct,2*nvclev,&
1120: !!$& kgrib,klengo,kword,0,ierr)
!!$#endif
!!$ CALL pbwrite(nunitg4x,kgrib,kword*iobyte,iret)
!!$
!!$ IF (iret /= kword*iobyte) THEN
1125: !!$ CALL finish('outgpx','I/O error - disks full?')
!!$ END IF
!!$
!!$ END DO
!!$ END DO
1130:
END SUBROUTINE outgpx
END MODULE mo_grib
Info Section
uses: mo_constants, mo_control, mo_decomposition, mo_doctor, mo_exception
mo_filename, mo_gaussgrid, mo_memory_g3b, mo_memory_gl, mo_memory_sp
mo_mpi, mo_parameters, mo_post, mo_start_dataset, mo_tracer
mo_transpose, mo_year
calls: codegb5, errmsg, finish, gather_gp, gather_sp
get_entry, get_info, gribex, grsvck, pbclose
pbopen, pbwrite, set_output_time
back to top
ECHAM 4 vf90 (C) 1998 Max-Planck-Institut für Meteorologie, Hamburg
Wed Nov 24 01:25:21 CST 1999
HTML derived from FORTRAN source by f2html.pl v0.3 (C) 1997,98 Beroud Jean-Marc.