691if (
associated(this%voldati))
write(unit=lunit)this%voldati
693if (.not.
present(unit))
close(unit=lunit)
695end subroutine volgrid6d_write_on_file
704subroutine volgrid6d_read_from_file (this,unit,filename,description,tarray,filename_auto)
707integer,
intent(inout),
optional :: unit
708character(len=*),
INTENT(in),
optional :: filename
709character(len=*),
intent(out),
optional :: filename_auto
710character(len=*),
INTENT(out),
optional :: description
711integer,
intent(out),
optional :: tarray(8)
713integer :: ntime, ntimerange, nlevel, nvar
715character(len=254) :: ldescription,lfilename,arg
716integer :: ltarray(8),lunit
717logical :: opened,exist
725if (.not.
present(unit))
then
736lfilename=trim(arg)//
".vg6d"
737if (
index(arg,
'/',back=.true.) > 0) lfilename=lfilename(
index(arg,
'/',back=.true.)+1 : )
739if (
present(filename))
then
740 if (filename /=
"")
then
745if (
present(filename_auto))filename_auto=lfilename
748inquire(unit=lunit,opened=opened)
749if (.not. opened)
then
750 inquire(file=lfilename,exist=exist)
751 IF (.NOT. exist)
CALL raise_fatal_error(
'file '//trim(lfilename)//
' does not exist, cannot open')
752 open (unit=lunit,file=lfilename,form=
"UNFORMATTED")
755read(unit=lunit)ldescription
756read(unit=lunit)ltarray
758call l4f_log(l4f_info,
"Info: reading volgrid6d from file: "//trim(lfilename))
759call l4f_log(l4f_info,
"Info: description: "//trim(ldescription))
762if (
present(description))description=ldescription
763if (
present(tarray))tarray=ltarray
767read(unit=lunit) ntime, ntimerange, nlevel, nvar
770call volgrid6d_alloc (this, &
771 ntime=ntime, ntimerange=ntimerange, nlevel=nlevel, nvar=nvar)
773call volgrid6d_alloc_vol (this)
775if (
associated(this%time))
call read_unit(this%time, lunit)
776if (
associated(this%level))
read(unit=lunit)this%level
777if (
associated(this%timerange))
read(unit=lunit)this%timerange
778if (
associated(this%var))
read(unit=lunit)this%var
783if (
associated(this%voldati))
read(unit=lunit)this%voldati
785if (.not.
present(unit))
close(unit=lunit)
787end subroutine volgrid6d_read_from_file
809SUBROUTINE import_from_gridinfo(this, gridinfo, force, dup_mode, clone, &
813LOGICAL,
INTENT(in),
OPTIONAL :: force
814INTEGER,
INTENT(in),
OPTIONAL :: dup_mode
815LOGICAL ,
INTENT(in),
OPTIONAL :: clone
816LOGICAL,
INTENT(IN),
OPTIONAL :: isanavar
818CHARACTER(len=255) :: type
819INTEGER :: itime0, itimerange0, itime1, itimerange1, itime, itimerange, &
820 ilevel, ivar, ldup_mode
823REAL,
ALLOCATABLE :: tmpgrid(:,:)
825IF (
PRESENT(dup_mode))
THEN
831call get_val(this%griddim,proj_type=type)
834call l4f_category_log(this%category,l4f_debug,
"import_from_gridinfo: "//trim(type))
837if (.not.
c_e(type))
then
838 call copy(gridinfo%griddim, this%griddim)
842 CALL volgrid6d_alloc_vol(this, ini=.true.)
844else if (.not. (this%griddim == gridinfo%griddim ))
then
847 "volgrid and gridinfo grid type or size are different, gridinfo rejected")
854ilevel =
index(this%level, gridinfo%level)
855IF (ilevel == 0 .AND. optio_log(force))
THEN
856 ilevel =
index(this%level, vol7d_level_miss)
857 IF (ilevel /= 0) this%level(ilevel) = gridinfo%level
862 "volgrid6d: level not valid for volume, gridinfo rejected")
867IF (optio_log(isanavar))
THEN
869 itime1 =
SIZE(this%time)
871 itimerange1 =
SIZE(this%timerange)
873 correctedtime = gridinfo%time
874 IF (this%time_definition == 1) correctedtime = correctedtime + &
875 timedelta_new(sec=gridinfo%timerange%p1)
876 itime0 =
index(this%time, correctedtime)
877 IF (itime0 == 0 .AND. optio_log(force))
THEN
878 itime0 =
index(this%time, datetime_miss)
879 IF (itime0 /= 0) this%time(itime0) = correctedtime
881 IF (itime0 == 0)
THEN
883 "volgrid6d: time not valid for volume, gridinfo rejected")
889 itimerange0 =
index(this%timerange,gridinfo%timerange)
890 IF (itimerange0 == 0 .AND. optio_log(force))
THEN
891 itimerange0 =
index(this%timerange, vol7d_timerange_miss)
892 IF (itimerange0 /= 0) this%timerange(itimerange0) = gridinfo%timerange
894 IF (itimerange0 == 0)
THEN
896 "volgrid6d: timerange not valid for volume, gridinfo rejected")
900 itimerange1 = itimerange0
903ivar =
index(this%var, gridinfo%var)
904IF (ivar == 0 .AND. optio_log(force))
THEN
905 ivar =
index(this%var, volgrid6d_var_miss)
906 IF (ivar /= 0) this%var(ivar) = gridinfo%var
910 "volgrid6d: var not valid for volume, gridinfo rejected")
915DO itimerange = itimerange0, itimerange1
916 DO itime = itime0, itime1
917 IF (
ASSOCIATED(this%gaid))
THEN
919 IF (
c_e(this%gaid(ilevel,itime,itimerange,ivar)))
THEN
923 IF (optio_log(
clone))
CALL delete(this%gaid(ilevel,itime,itimerange,ivar))
926 IF (optio_log(
clone))
THEN
927 CALL copy(gridinfo%gaid, this%gaid(ilevel,itime,itimerange,ivar))
932 this%gaid(ilevel,itime,itimerange,ivar) = gridinfo%gaid
935 IF (
ASSOCIATED(this%voldati))
THEN
936 IF (.NOT.dup .OR. ldup_mode == 0)
THEN
937 this%voldati(:,:,ilevel,itime,itimerange,ivar) =
decode_gridinfo(gridinfo)
938 ELSE IF (ldup_mode == 1)
THEN
941 this%voldati(:,:,ilevel,itime,itimerange,ivar) = tmpgrid(:,:)
943 ELSE IF (ldup_mode == 2)
THEN
944 WHERE(.NOT.
c_e(this%voldati(:,:,ilevel,itime,itimerange,ivar)))
945 this%voldati(:,:,ilevel,itime,itimerange,ivar) =
decode_gridinfo(gridinfo)
952 "gaid not allocated, you probably need to call volgrid6d_alloc_vol first")
960END SUBROUTINE import_from_gridinfo
967SUBROUTINE export_to_gridinfo(this, gridinfo, itime, itimerange, ilevel, ivar, &
968 gaid_template, clone)
975TYPE(
grid_id),
INTENT(in),
OPTIONAL :: gaid_template
976LOGICAL,
INTENT(in),
OPTIONAL :: clone
979LOGICAL :: usetemplate
980REAL,
POINTER :: voldati(:,:)
987IF (.NOT.
c_e(this%gaid(ilevel,itime,itimerange,ivar)))
THEN
989 CALL l4f_category_log(this%category,l4f_debug,
"empty gaid found, skipping export")
995IF (
PRESENT(gaid_template))
THEN
996 CALL copy(gaid_template, gaid)
998 CALL l4f_category_log(this%category,l4f_debug,
"template cloned to a new gaid")
1000 usetemplate =
c_e(gaid)
1003IF (.NOT.usetemplate)
THEN
1004 IF (optio_log(
clone))
THEN
1005 CALL copy(this%gaid(ilevel,itime,itimerange,ivar), gaid)
1007 CALL l4f_category_log(this%category,l4f_debug,
"original gaid cloned to a new one")
1010 gaid = this%gaid(ilevel,itime,itimerange,ivar)
1014IF (this%time_definition == 1)
THEN
1015 correctedtime = this%time(itime) - &
1016 timedelta_new(sec=this%timerange(itimerange)%p1)
1018 correctedtime = this%time(itime)
1021CALL init(gridinfo,gaid, this%griddim, correctedtime, this%timerange(itimerange), &
1022 this%level(ilevel), this%var(ivar))
1025CALL export(gridinfo%griddim, gridinfo%gaid)
1027IF (
ASSOCIATED(this%voldati))
THEN
1028 CALL encode_gridinfo(gridinfo, this%voldati(:,:,ilevel,itime,itimerange,ivar))
1029ELSE IF (usetemplate)
THEN
1031 CALL volgrid_get_vol_2d(this, ilevel, itime, itimerange, ivar, voldati)
1036END SUBROUTINE export_to_gridinfo
1056SUBROUTINE import_from_gridinfovv(this, gridinfov, dup_mode, clone, decode, &
1057 time_definition, anavar, categoryappend)
1060INTEGER,
INTENT(in),
OPTIONAL :: dup_mode
1061LOGICAL ,
INTENT(in),
OPTIONAL :: clone
1062LOGICAL,
INTENT(in),
OPTIONAL :: decode
1063INTEGER,
INTENT(IN),
OPTIONAL :: time_definition
1064CHARACTER(len=*),
INTENT(IN),
OPTIONAL :: anavar(:)
1065CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
1067INTEGER :: i, j, stallo
1068INTEGER :: ngrid, ntime, ntimerange, nlevel, nvar
1070CHARACTER(len=512) :: a_name
1071TYPE(
datetime),
ALLOCATABLE :: correctedtime(:)
1072LOGICAL,
ALLOCATABLE :: isanavar(:)
1073TYPE(vol7d_var) :: lvar
1076if (
present(categoryappend))
then
1077 call l4f_launcher(a_name,a_name_append=trim(subcategory)//
"."//trim(categoryappend))
1079 call l4f_launcher(a_name,a_name_append=trim(subcategory))
1081category=l4f_category_get(a_name)
1087ngrid=count_distinct(gridinfov%array(1:gridinfov%arraysize)%griddim,back=.true.)
1089 ' different grid definition(s) found in input data')
1091ALLOCATE(this(ngrid),stat=stallo)
1094 CALL raise_fatal_error()
1097 IF (
PRESENT(categoryappend))
THEN
1098 CALL init(this(i), time_definition=time_definition, categoryappend=trim(categoryappend)//
"-vol"//
t2c(i))
1100 CALL init(this(i), time_definition=time_definition, categoryappend=
"vol"//
t2c(i))
1104this(:)%griddim=pack_distinct(gridinfov%array(1:gridinfov%arraysize)%griddim, &
1108ALLOCATE(isanavar(gridinfov%arraysize))
1109isanavar(:) = .false.
1110IF (
PRESENT(anavar))
THEN
1111 DO i = 1, gridinfov%arraysize
1112 DO j = 1,
SIZE(anavar)
1113 lvar =
convert(gridinfov%array(i)%var)
1114 IF (lvar%btable == anavar(j))
THEN
1115 isanavar(i) = .true.
1121 t2c(gridinfov%arraysize)//
' constant-data messages found in input data')
1125ALLOCATE(correctedtime(gridinfov%arraysize))
1126correctedtime(:) = gridinfov%array(1:gridinfov%arraysize)%time
1127IF (
PRESENT(time_definition))
THEN
1128 IF (time_definition == 1)
THEN
1129 DO i = 1, gridinfov%arraysize
1130 correctedtime(i) = correctedtime(i) + &
1131 timedelta_new(sec=gridinfov%array(i)%timerange%p1)
1137 IF (
PRESENT(anavar))
THEN
1138 j = count((this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim) &
1139 .AND. .NOT.isanavar(:))
1142 ' has only constant data, this is not allowed')
1144 CALL raise_fatal_error()
1147 ntime = count_distinct(correctedtime, &
1148 mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim) &
1149 .AND. .NOT.isanavar(:), back=.true.)
1150 ntimerange = count_distinct(gridinfov%array(1:gridinfov%arraysize)%timerange, &
1151 mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim) &
1152 .AND. .NOT.isanavar(:), back=.true.)
1153 nlevel = count_distinct(gridinfov%array(1:gridinfov%arraysize)%level, &
1154 mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim), &
1156 nvar = count_distinct(gridinfov%array(1:gridinfov%arraysize)%var, &
1157 mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim), &
1164 CALL volgrid6d_alloc(this(i),this(i)%griddim%dim,ntime=ntime, &
1165 ntimerange=ntimerange,nlevel=nlevel,nvar=nvar)
1167 this(i)%time = pack_distinct(correctedtime, ntime, &
1168 mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim) &
1169 .AND. .NOT.isanavar(:), back=.true.)
1170 CALL sort(this(i)%time)
1172 this(i)%timerange = pack_distinct(gridinfov%array( &
1173 1:gridinfov%arraysize)%timerange, ntimerange, &
1174 mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim) &
1175 .AND. .NOT.isanavar(:), back=.true.)
1176 CALL sort(this(i)%timerange)
1178 this(i)%level=pack_distinct(gridinfov%array(1:gridinfov%arraysize)%level, &
1179 nlevel,mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim), &
1181 CALL sort(this(i)%level)
1183 this(i)%var=pack_distinct(gridinfov%array(1:gridinfov%arraysize)%var, nvar, &
1184 mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim), &
1190 CALL volgrid6d_alloc_vol(this(i), decode=decode)
1194DEALLOCATE(correctedtime)
1196DO i = 1, gridinfov%arraysize
1201 "to volgrid6d index: "//
t2c(
index(this%griddim, gridinfov%array(i)%griddim)))
1204 CALL import(this(
index(this%griddim, gridinfov%array(i)%griddim)), &
1205 gridinfov%array(i), dup_mode=dup_mode,
clone=
clone, isanavar=isanavar(i))
1210CALL l4f_category_delete(category)
1212END SUBROUTINE import_from_gridinfovv
1220SUBROUTINE export_to_gridinfov(this, gridinfov, gaid_template, clone)
1223TYPE(
grid_id),
INTENT(in),
OPTIONAL :: gaid_template
1224LOGICAL,
INTENT(in),
OPTIONAL :: clone
1226INTEGER :: i ,itime, itimerange, ilevel, ivar
1227INTEGER :: ntime, ntimerange, nlevel, nvar
1237CALL dealloc(this%griddim%dim)
1240ntime=
size(this%time)
1241ntimerange=
size(this%timerange)
1242nlevel=
size(this%level)
1246 DO itimerange=1,ntimerange
1250 CALL init(gridinfol)
1251 CALL export(this, gridinfol, itime, itimerange, ilevel, ivar, &
1253 IF (
c_e(gridinfol%gaid))
THEN
1254 CALL insert(gridinfov, gridinfol)
1264END SUBROUTINE export_to_gridinfov
1272SUBROUTINE export_to_gridinfovv(this, gridinfov, gaid_template, clone)
1277TYPE(
grid_id),
INTENT(in),
OPTIONAL :: gaid_template
1278LOGICAL,
INTENT(in),
OPTIONAL :: clone
1285 "export_to_gridinfovv grid index: "//
t2c(i))
1290END SUBROUTINE export_to_gridinfovv
1302SUBROUTINE volgrid6d_import_from_file(this, filename, dup_mode, decode, &
1303 time_definition, anavar, categoryappend)
1305CHARACTER(len=*),
INTENT(in) :: filename
1306INTEGER,
INTENT(in),
OPTIONAL :: dup_mode
1307LOGICAL,
INTENT(in),
OPTIONAL :: decode
1308INTEGER,
INTENT(IN),
OPTIONAL :: time_definition
1309CHARACTER(len=*),
INTENT(IN),
OPTIONAL :: anavar(:)
1310character(len=*),
INTENT(in),
OPTIONAL :: categoryappend
1314CHARACTER(len=512) :: a_name
1318IF (
PRESENT(categoryappend))
THEN
1319 CALL l4f_launcher(a_name,a_name_append= &
1320 trim(subcategory)//
"."//trim(categoryappend))
1322 CALL l4f_launcher(a_name,a_name_append=trim(subcategory))
1324category=l4f_category_get(a_name)
1326CALL import(gridinfo, filename=filename, categoryappend=categoryappend)
1328IF (gridinfo%arraysize > 0)
THEN
1330 CALL import(this, gridinfo, dup_mode=dup_mode,
clone=.true., decode=decode, &
1331 time_definition=time_definition, anavar=anavar, &
1332 categoryappend=categoryappend)
1338 CALL l4f_category_log(category,l4f_info,
"file does not contain gridded data")
1342CALL l4f_category_delete(category)
1344END SUBROUTINE volgrid6d_import_from_file
1354SUBROUTINE volgrid6d_export_to_file(this, filename, gaid_template, categoryappend)
1356CHARACTER(len=*),
INTENT(in) :: filename
1357TYPE(
grid_id),
INTENT(in),
OPTIONAL :: gaid_template
1358character(len=*),
INTENT(in),
OPTIONAL :: categoryappend
1362CHARACTER(len=512) :: a_name
1364IF (
PRESENT(categoryappend))
THEN
1365 CALL l4f_launcher(a_name,a_name_append=trim(subcategory)//
"."//trim(categoryappend))
1367 CALL l4f_launcher(a_name,a_name_append=trim(subcategory))
1369category=l4f_category_get(a_name)
1375CALL l4f_category_log(category,l4f_info,
"writing volgrid6d to grib file: "//trim(filename))
1378 CALL export(this, gridinfo, gaid_template=gaid_template,
clone=.true.)
1379 IF (gridinfo%arraysize > 0)
THEN
1380 CALL export(gridinfo, filename)
1388CALL l4f_category_delete(category)
1390END SUBROUTINE volgrid6d_export_to_file
1396SUBROUTINE volgrid6dv_delete(this)
1401IF (
ASSOCIATED(this))
THEN
1402 DO i = 1,
SIZE(this)
1405 "delete volgrid6d vector index: "//trim(
to_char(i)))
1412END SUBROUTINE volgrid6dv_delete
1416SUBROUTINE volgrid6d_transform_compute(this, volgrid6d_in, volgrid6d_out, &
1417 lev_out, var_coord_vol, clone)
1419type(
volgrid6d),
INTENT(in) :: volgrid6d_in
1420type(
volgrid6d),
INTENT(inout) :: volgrid6d_out
1421TYPE(
vol7d_level),
INTENT(in),
OPTIONAL :: lev_out(:)
1422INTEGER,
INTENT(in),
OPTIONAL :: var_coord_vol
1423LOGICAL,
INTENT(in),
OPTIONAL :: clone
1425INTEGER :: ntime, ntimerange, inlevel, onlevel, nvar, &
1426 itime, itimerange, ilevel, ivar, levshift, levused, lvar_coord_vol, spos
1427REAL,
POINTER :: voldatiin(:,:,:), voldatiout(:,:,:), coord_3d_in(:,:,:)
1432call l4f_category_log(volgrid6d_in%category,l4f_debug,
"start volgrid6d_transform_compute")
1440lvar_coord_vol = optio_i(var_coord_vol)
1442if (
associated(volgrid6d_in%time))
then
1443 ntime=
size(volgrid6d_in%time)
1444 volgrid6d_out%time=volgrid6d_in%time
1447if (
associated(volgrid6d_in%timerange))
then
1448 ntimerange=
size(volgrid6d_in%timerange)
1449 volgrid6d_out%timerange=volgrid6d_in%timerange
1452IF (
ASSOCIATED(volgrid6d_in%level))
THEN
1453 inlevel=
SIZE(volgrid6d_in%level)
1455IF (
PRESENT(lev_out))
THEN
1456 onlevel=
SIZE(lev_out)
1457 volgrid6d_out%level=lev_out
1458ELSE IF (
ASSOCIATED(volgrid6d_in%level))
THEN
1459 onlevel=
SIZE(volgrid6d_in%level)
1460 volgrid6d_out%level=volgrid6d_in%level
1463if (
associated(volgrid6d_in%var))
then
1464 nvar=
size(volgrid6d_in%var)
1465 volgrid6d_out%var=volgrid6d_in%var
1468IF (.NOT.
ASSOCIATED(volgrid6d_in%voldati))
THEN
1469 ALLOCATE(voldatiin(volgrid6d_in%griddim%dim%nx, volgrid6d_in%griddim%dim%ny, &
1472IF (.NOT.
ASSOCIATED(volgrid6d_out%voldati))
THEN
1473 ALLOCATE(voldatiout(volgrid6d_out%griddim%dim%nx, volgrid6d_out%griddim%dim%ny, &
1477CALL get_val(this, levshift=levshift, levused=levused)
1479IF (
c_e(lvar_coord_vol))
THEN
1480 CALL get_val(this%trans, output_levtype=output_levtype)
1481 IF (output_levtype%level1 == 103 .OR. output_levtype%level1 == 108)
THEN
1482 spos = firsttrue(volgrid6d_in%level(:) == vol7d_level_new(1))
1485 'output level '//
t2c(output_levtype%level1)// &
1486 ' requested, but height/press of surface not provided in volume')
1488 IF (.NOT.
c_e(levshift) .AND. .NOT.
c_e(levused))
THEN
1490 'internal inconsistence, levshift and levused undefined when they should be')
1499 DO itimerange=1,ntimerange
1502 IF (
c_e(levshift) .AND.
c_e(levused))
THEN
1504 volgrid6d_in%gaid(levshift+1:levshift+levused,itime,itimerange,ivar) &
1507 DO ilevel = 1, min(inlevel,onlevel)
1509 IF (
c_e(volgrid6d_in%gaid(ilevel,itime,itimerange,ivar)) .AND. .NOT. &
1510 c_e(volgrid6d_out%gaid(ilevel,itime,itimerange,ivar)))
THEN
1512 IF (optio_log(
clone))
THEN
1513 CALL copy(volgrid6d_in%gaid(ilevel,itime,itimerange,ivar),&
1514 volgrid6d_out%gaid(ilevel,itime,itimerange,ivar))
1517 "cloning gaid, level "//
t2c(ilevel))
1520 volgrid6d_out%gaid(ilevel,itime,itimerange,ivar) = &
1521 volgrid6d_in%gaid(ilevel,itime,itimerange,ivar)
1526 DO ilevel = min(inlevel,onlevel) + 1, onlevel
1527 IF (
c_e(volgrid6d_in%gaid(inlevel,itime,itimerange,ivar)) .AND. .NOT. &
1528 c_e(volgrid6d_out%gaid(ilevel,itime,itimerange,ivar)))
then
1530 CALL copy(volgrid6d_in%gaid(inlevel,itime,itimerange,ivar),&
1531 volgrid6d_out%gaid(ilevel,itime,itimerange,ivar))
1534 "forced cloning gaid, level "//
t2c(inlevel)//
"->"//
t2c(ilevel))
1539 IF (
c_e(lvar_coord_vol))
THEN
1540 NULLIFY(coord_3d_in)
1541 CALL volgrid_get_vol_3d(volgrid6d_in, itime, itimerange, lvar_coord_vol, &
1545 coord_3d_in(:,:,levshift+1:levshift+levused) = rmiss
1547 DO ilevel = levshift+1, levshift+levused
1548 WHERE(
c_e(coord_3d_in(:,:,ilevel)) .AND.
c_e(coord_3d_in(:,:,spos)))
1549 coord_3d_in(:,:,ilevel) = coord_3d_in(:,:,ilevel) - &
1550 coord_3d_in(:,:,spos)
1552 coord_3d_in(:,:,ilevel) = rmiss
1558 CALL volgrid_get_vol_3d(volgrid6d_in, itime, itimerange, ivar, &
1560 IF (
ASSOCIATED(volgrid6d_out%voldati)) &
1561 CALL volgrid_get_vol_3d(volgrid6d_out, itime, itimerange, ivar, &
1563 IF (
c_e(lvar_coord_vol))
THEN
1564 CALL compute(this, voldatiin, voldatiout,
convert(volgrid6d_in%var(ivar)), &
1565 coord_3d_in(:,:,levshift+1:levshift+levused))
1567 CALL compute(this, voldatiin, voldatiout,
convert(volgrid6d_in%var(ivar)))
1569 CALL volgrid_set_vol_3d(volgrid6d_out, itime, itimerange, ivar, &
1575IF (
c_e(lvar_coord_vol))
THEN
1576 DEALLOCATE(coord_3d_in)
1578IF (.NOT.
ASSOCIATED(volgrid6d_in%voldati))
THEN
1579 DEALLOCATE(voldatiin)
1581IF (.NOT.
ASSOCIATED(volgrid6d_out%voldati))
THEN
1582 DEALLOCATE(voldatiout)
1586END SUBROUTINE volgrid6d_transform_compute
1595SUBROUTINE volgrid6d_transform(this, griddim, volgrid6d_in, volgrid6d_out, &
1596 lev_out, volgrid6d_coord_in, maskgrid, maskbounds, clone, decode, categoryappend)
1599TYPE(
volgrid6d),
INTENT(inout) :: volgrid6d_in
1600TYPE(
volgrid6d),
INTENT(out) :: volgrid6d_out
1601TYPE(
vol7d_level),
INTENT(in),
OPTIONAL,
TARGET :: lev_out(:)
1602TYPE(
volgrid6d),
INTENT(in),
OPTIONAL :: volgrid6d_coord_in
1603REAL,
INTENT(in),
OPTIONAL :: maskgrid(:,:)
1604REAL,
INTENT(in),
OPTIONAL :: maskbounds(:)
1605LOGICAL,
INTENT(in),
OPTIONAL :: clone
1606LOGICAL,
INTENT(in),
OPTIONAL :: decode
1607CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
1612TYPE(vol7d_var) :: vcoord_var
1613INTEGER :: i, k, ntime, ntimerange, nlevel, nvar, var_coord_in, var_coord_vol, &
1614 cf_out, nxc, nyc, nxi, nyi, i3, i4, i5, i6, &
1615 ulstart, ulend, spos
1616REAL,
ALLOCATABLE :: coord_3d_in(:,:,:)
1617TYPE(geo_proj) :: proj_in, proj_out
1618CHARACTER(len=80) :: trans_type
1620LOGICAL,
ALLOCATABLE :: mask_in(:)
1623call l4f_category_log(volgrid6d_in%category, l4f_debug,
"start volgrid6d_transform")
1631if (
associated(volgrid6d_in%time)) ntime=
size(volgrid6d_in%time)
1632if (
associated(volgrid6d_in%timerange)) ntimerange=
size(volgrid6d_in%timerange)
1633if (
associated(volgrid6d_in%level)) nlevel=
size(volgrid6d_in%level)
1634if (
associated(volgrid6d_in%var)) nvar=
size(volgrid6d_in%var)
1636IF (ntime == 0 .OR. ntimerange == 0 .OR. nlevel == 0 .OR. nvar == 0)
THEN
1638 "trying to transform an incomplete volgrid6d object, ntime="//
t2c(ntime)// &
1639 ' ntimerange='//
t2c(ntimerange)//
' nlevel='//
t2c(nlevel)//
' nvar='//
t2c(nvar))
1640 CALL init(volgrid6d_out)
1645CALL get_val(this, trans_type=trans_type)
1649IF (
PRESENT(griddim) .AND. (trans_type ==
'inter' .OR. trans_type ==
'boxinter' &
1650 .OR. trans_type ==
'stencilinter'))
THEN
1652 CALL get_val(griddim, component_flag=cf_out,
proj=proj_out)
1654 IF (proj_in /= proj_out)
CALL vg6d_wind_unrot(volgrid6d_in)
1655ELSE IF (
PRESENT(griddim))
THEN
1656 CALL get_val(griddim, component_flag=cf_out)
1661var_coord_vol = imiss
1662IF (trans_type ==
'vertint')
THEN
1663 IF (
PRESENT(lev_out))
THEN
1666 IF (
PRESENT(volgrid6d_coord_in))
THEN
1667 IF (
ASSOCIATED(volgrid6d_coord_in%voldati))
THEN
1670 IF (
SIZE(volgrid6d_coord_in%voldati,4) /= 1 .OR. &
1671 SIZE(volgrid6d_coord_in%voldati,5) /= 1)
THEN
1673 'volume providing constant input vertical coordinate must have &
1674 &only 1 time and 1 timerange')
1675 CALL init(volgrid6d_out)
1681 CALL get_val(this, output_levtype=output_levtype)
1683 IF (.NOT.
c_e(vcoord_var))
THEN
1685 'requested output level type '//
t2c(output_levtype%level1)// &
1686 ' does not correspond to any known physical variable for &
1687 &providing vertical coordinate')
1688 CALL init(volgrid6d_out)
1693 DO i = 1,
SIZE(volgrid6d_coord_in%var)
1694 IF (
convert(volgrid6d_coord_in%var(i)) == vcoord_var)
THEN
1700 IF (.NOT.
c_e(var_coord_in))
THEN
1702 'volume providing constant input vertical coordinate contains no &
1703 &variables matching output level type '//
t2c(output_levtype%level1))
1704 CALL init(volgrid6d_out)
1709 'Coordinate for vertint found in coord volume at position '// &
1715 CALL get_val(volgrid6d_coord_in%griddim, nx=nxc, ny=nyc)
1716 CALL get_val(volgrid6d_in%griddim, nx=nxi, ny=nyi)
1717 IF (nxc /= nxi .OR. nyc /= nyi)
THEN
1719 'volume providing constant input vertical coordinate must have &
1720 &the same grid as the input')
1722 'vertical coordinate: '//
t2c(nxc)//
'x'//
t2c(nyc)// &
1723 ', input volume: '//
t2c(nxi)//
'x'//
t2c(nyi))
1724 CALL init(volgrid6d_out)
1730 CALL get_val(this, input_levtype=input_levtype)
1732 (volgrid6d_coord_in%level(:)%level1 == input_levtype%level1) .AND. &
1733 (volgrid6d_coord_in%level(:)%level2 == input_levtype%level2)
1734 ulstart = firsttrue(mask_in)
1735 ulend = lasttrue(mask_in)
1736 IF (ulstart == 0 .OR. ulend == 0)
THEN
1738 'coordinate file does not contain levels of type '// &
1739 t2c(input_levtype%level1)//
'/'//
t2c(input_levtype%level2)// &
1740 ' specified for input data')
1741 CALL init(volgrid6d_out)
1746 coord_3d_in = volgrid6d_coord_in%voldati(:,:,ulstart:ulend,1,1,var_coord_in)
1748 IF (output_levtype%level1 == 103 .OR. &
1749 output_levtype%level1 == 108)
THEN
1750 spos = firsttrue(volgrid6d_coord_in%level(:) == vol7d_level_new(1))
1753 'output level '//
t2c(output_levtype%level1)// &
1754 ' requested, but height/press of surface not provided in coordinate file')
1755 CALL init(volgrid6d_out)
1759 DO k = 1,
SIZE(coord_3d_in,3)
1760 WHERE(
c_e(coord_3d_in(:,:,k)) .AND. &
1761 c_e(volgrid6d_coord_in%voldati(:,:,spos,1,1,var_coord_in)))
1762 coord_3d_in(:,:,k) = coord_3d_in(:,:,k) - &
1763 volgrid6d_coord_in%voldati(:,:,spos,1,1,var_coord_in)
1765 coord_3d_in(:,:,k) = rmiss
1773 IF (.NOT.
c_e(var_coord_in))
THEN
1775 CALL get_val(this, output_levtype=output_levtype)
1777 IF (
c_e(vcoord_var))
THEN
1778 DO i = 1,
SIZE(volgrid6d_in%var)
1779 IF (
convert(volgrid6d_in%var(i)) == vcoord_var)
THEN
1785 IF (
c_e(var_coord_vol))
THEN
1787 'Coordinate for vertint found in input volume at position '// &
1794 CALL init(volgrid6d_out, griddim=volgrid6d_in%griddim, &
1795 time_definition=volgrid6d_in%time_definition, categoryappend=categoryappend)
1796 IF (
c_e(var_coord_in))
THEN
1797 CALL init(grid_trans, this, lev_in=volgrid6d_in%level, lev_out=lev_out, &
1798 coord_3d_in=coord_3d_in, categoryappend=categoryappend)
1800 CALL init(grid_trans, this, lev_in=volgrid6d_in%level, lev_out=lev_out, &
1801 categoryappend=categoryappend)
1804 CALL get_val(grid_trans, output_level_auto=llev_out)
1805 IF (.NOT.
ASSOCIATED(llev_out)) llev_out => lev_out
1806 nlevel =
SIZE(llev_out)
1809 'volgrid6d_transform: vertint requested but lev_out not provided')
1810 CALL init(volgrid6d_out)
1816 CALL init(volgrid6d_out, griddim=griddim, &
1817 time_definition=volgrid6d_in%time_definition, categoryappend=categoryappend)
1818 CALL init(grid_trans, this, in=volgrid6d_in%griddim, out=volgrid6d_out%griddim, &
1819 maskgrid=maskgrid, maskbounds=maskbounds, categoryappend=categoryappend)
1823IF (
c_e(grid_trans))
THEN
1825 CALL volgrid6d_alloc(volgrid6d_out, ntime=ntime, nlevel=nlevel, &
1826 ntimerange=ntimerange, nvar=nvar)
1828 IF (
PRESENT(decode))
THEN
1831 ldecode =
ASSOCIATED(volgrid6d_in%voldati)
1834 decode_loop:
DO i6 = 1,nvar
1835 DO i5 = 1, ntimerange
1838 IF (
c_e(volgrid6d_in%gaid(i3,i4,i5,i6)))
THEN
1839 ldecode = ldecode .OR. grid_id_readonly(volgrid6d_in%gaid(i3,i4,i5,i6))
1847 IF (
PRESENT(decode))
THEN
1848 IF (ldecode.NEQV.decode)
THEN
1850 'volgrid6d_transform: decode status forced to .TRUE. because driver does not allow copy')
1854 CALL volgrid6d_alloc_vol(volgrid6d_out, decode=ldecode)
1859 IF (trans_type ==
'vertint')
THEN
1862 "volgrid6d_transform: vertint to "//
t2c(nlevel)//
" levels")
1864 CALL compute(grid_trans, volgrid6d_in, volgrid6d_out, lev_out=llev_out, &
1867 CALL compute(grid_trans, volgrid6d_in, volgrid6d_out,
clone=
clone)
1870 IF (cf_out == 0)
THEN
1871 CALL wind_unrot(volgrid6d_out)
1872 ELSE IF (cf_out == 1)
THEN
1873 CALL wind_rot(volgrid6d_out)
1879 'volgrid6d_transform: transformation not valid')
1885END SUBROUTINE volgrid6d_transform
1896SUBROUTINE volgrid6dv_transform(this, griddim, volgrid6d_in, volgrid6d_out, &
1897 lev_out, volgrid6d_coord_in, maskgrid, maskbounds, clone, decode, categoryappend)
1900TYPE(
volgrid6d),
INTENT(inout) :: volgrid6d_in(:)
1901TYPE(
volgrid6d),
POINTER :: volgrid6d_out(:)
1902TYPE(
vol7d_level),
INTENT(in),
OPTIONAL :: lev_out(:)
1903TYPE(
volgrid6d),
INTENT(in),
OPTIONAL :: volgrid6d_coord_in
1904REAL,
INTENT(in),
OPTIONAL :: maskgrid(:,:)
1905REAL,
INTENT(in),
OPTIONAL :: maskbounds(:)
1906LOGICAL,
INTENT(in),
OPTIONAL :: clone
1907LOGICAL,
INTENT(in),
OPTIONAL :: decode
1908CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
1913allocate(volgrid6d_out(
size(volgrid6d_in)),stat=stallo)
1915 call l4f_log(l4f_fatal,
"allocating memory")
1916 call raise_fatal_error()
1919do i=1,
size(volgrid6d_in)
1920 call transform(this, griddim, volgrid6d_in(i), volgrid6d_out(i), &
1921 lev_out=lev_out, volgrid6d_coord_in=volgrid6d_coord_in, &
1922 maskgrid=maskgrid, maskbounds=maskbounds, &
1923 clone=
clone, decode=decode, categoryappend=categoryappend)
1926END SUBROUTINE volgrid6dv_transform
1930SUBROUTINE volgrid6d_v7d_transform_compute(this, volgrid6d_in, vol7d_out, &
1931 networkname, noconvert)
1933type(
volgrid6d),
INTENT(in) :: volgrid6d_in
1934type(
vol7d),
INTENT(inout) :: vol7d_out
1935CHARACTER(len=*),
OPTIONAL,
INTENT(in) :: networkname
1936LOGICAL,
OPTIONAL,
INTENT(in) :: noconvert
1938INTEGER :: nntime, nana, ntime, ntimerange, nlevel, nvar, stallo
1939INTEGER :: itime, itimerange, ivar, inetwork
1940REAL,
ALLOCATABLE :: voldatir_out(:,:,:)
1942TYPE(
datetime),
ALLOCATABLE :: validitytime(:,:)
1943REAL,
POINTER :: voldatiin(:,:,:)
1946call l4f_category_log(volgrid6d_in%category,l4f_debug,
"start volgrid6d_v7d_transform_compute")
1955if (
present(networkname))
then
1956 call init(vol7d_out%network(1),name=networkname)
1958 call init(vol7d_out%network(1),name=
'generic')
1961if (
associated(volgrid6d_in%timerange))
then
1962 ntimerange=
size(volgrid6d_in%timerange)
1963 vol7d_out%timerange=volgrid6d_in%timerange
1966if (
associated(volgrid6d_in%time))
then
1967 ntime=
size(volgrid6d_in%time)
1969 if (vol7d_out%time_definition == volgrid6d_in%time_definition)
then
1972 vol7d_out%time=volgrid6d_in%time
1976 allocate (validitytime(ntime,ntimerange),stat=stallo)
1979 call raise_fatal_error()
1983 do itimerange=1,ntimerange
1984 if (vol7d_out%time_definition > volgrid6d_in%time_definition)
then
1985 validitytime(itime,itimerange) = &
1986 volgrid6d_in%time(itime) + timedelta_new(sec=volgrid6d_in%timerange(itimerange)%p1)
1988 validitytime(itime,itimerange) = &
1989 volgrid6d_in%time(itime) - timedelta_new(sec=volgrid6d_in%timerange(itimerange)%p1)
1994 nntime = count_distinct(reshape(validitytime,(/ntime*ntimerange/)), back=.true.)
1995 vol7d_out%time=pack_distinct(reshape(validitytime,(/ntime*ntimerange/)), nntime,back=.true.)
2000IF (
ASSOCIATED(volgrid6d_in%level))
THEN
2001 nlevel =
SIZE(volgrid6d_in%level)
2002 vol7d_out%level=volgrid6d_in%level
2005IF (
ASSOCIATED(volgrid6d_in%var))
THEN
2006 nvar =
SIZE(volgrid6d_in%var)
2007 IF (.NOT. optio_log(noconvert))
THEN
2008 CALL vargrib2varbufr(volgrid6d_in%var, vol7d_out%dativar%r, c_func)
2012nana =
SIZE(vol7d_out%ana)
2015IF (.NOT.
ASSOCIATED(volgrid6d_in%voldati))
THEN
2016 ALLOCATE(voldatiin(volgrid6d_in%griddim%dim%nx, volgrid6d_in%griddim%dim%ny, &
2020ALLOCATE(voldatir_out(nana,1,nlevel),stat=stallo)
2021IF (stallo /= 0)
THEN
2023 CALL raise_fatal_error()
2028 do itimerange=1,ntimerange
2040 CALL volgrid_get_vol_3d(volgrid6d_in, itime, itimerange, ivar, &
2043 CALL compute(this, voldatiin, voldatir_out, vol7d_out%dativar%r(ivar))
2045 if (vol7d_out%time_definition == volgrid6d_in%time_definition)
then
2046 vol7d_out%voldatir(:,itime,:,itimerange,ivar,inetwork) = &
2049 vol7d_out%voldatir(:,
index(vol7d_out%time,validitytime(itime,itimerange)),:,itimerange,ivar,inetwork)=&
2050 reshape(voldatir_out,(/nana,nlevel/))
2065deallocate(voldatir_out)
2066IF (.NOT.
ASSOCIATED(volgrid6d_in%voldati))
THEN
2067 DEALLOCATE(voldatiin)
2069if (
allocated(validitytime))
deallocate(validitytime)
2072IF (
ASSOCIATED(c_func))
THEN
2074 CALL compute(c_func(ivar), vol7d_out%voldatir(:,:,:,:,ivar,:))
2079end SUBROUTINE volgrid6d_v7d_transform_compute
2088SUBROUTINE volgrid6d_v7d_transform(this, volgrid6d_in, vol7d_out, v7d, &
2089 maskgrid, maskbounds, networkname, noconvert, find_index, categoryappend)
2091TYPE(
volgrid6d),
INTENT(inout) :: volgrid6d_in
2092TYPE(
vol7d),
INTENT(out) :: vol7d_out
2093TYPE(
vol7d),
INTENT(in),
OPTIONAL :: v7d
2094REAL,
INTENT(in),
OPTIONAL :: maskgrid(:,:)
2095REAL,
INTENT(in),
OPTIONAL :: maskbounds(:)
2096CHARACTER(len=*),
OPTIONAL,
INTENT(in) :: networkname
2097LOGICAL,
OPTIONAL,
INTENT(in) :: noconvert
2098PROCEDURE(basic_find_index),
POINTER,
OPTIONAL :: find_index
2099CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
2102INTEGER :: ntime, ntimerange, nlevel, nvar, nana, time_definition, nnetwork, stallo
2103INTEGER :: itime, itimerange, inetwork
2104TYPE(
datetime),
ALLOCATABLE :: validitytime(:,:)
2105INTEGER,
ALLOCATABLE :: point_index(:)
2106TYPE(
vol7d) :: v7d_locana
2109call l4f_category_log(volgrid6d_in%category,l4f_debug,
"start volgrid6d_v7d_transform")
2112call vg6d_wind_unrot(volgrid6d_in)
2120call get_val(this,time_definition=time_definition)
2121if (.not.
c_e(time_definition))
then
2125IF (
PRESENT(v7d))
THEN
2126 CALL vol7d_copy(v7d, v7d_locana)
2128 CALL init(v7d_locana)
2131if (
associated(volgrid6d_in%timerange)) ntimerange=
size(volgrid6d_in%timerange)
2133if (
associated(volgrid6d_in%time))
then
2135 ntime=
size(volgrid6d_in%time)
2137 if (time_definition /= volgrid6d_in%time_definition)
then
2140 allocate (validitytime(ntime,ntimerange),stat=stallo)
2143 call raise_fatal_error()
2147 do itimerange=1,ntimerange
2148 if (time_definition > volgrid6d_in%time_definition)
then
2149 validitytime(itime,itimerange) = &
2150 volgrid6d_in%time(itime) + timedelta_new(sec=volgrid6d_in%timerange(itimerange)%p1)
2152 validitytime(itime,itimerange) = &
2153 volgrid6d_in%time(itime) - timedelta_new(sec=volgrid6d_in%timerange(itimerange)%p1)
2158 ntime = count_distinct(reshape(validitytime,(/ntime*ntimerange/)), back=.true.)
2159 deallocate (validitytime)
2165if (
associated(volgrid6d_in%level)) nlevel=
size(volgrid6d_in%level)
2166if (
associated(volgrid6d_in%var)) nvar=
size(volgrid6d_in%var)
2168CALL init(grid_trans, this, volgrid6d_in%griddim, v7d_locana, &
2169 maskgrid=maskgrid, maskbounds=maskbounds, find_index=find_index, &
2170 categoryappend=categoryappend)
2171CALL init (vol7d_out,time_definition=time_definition)
2173IF (
c_e(grid_trans))
THEN
2175 nana=
SIZE(v7d_locana%ana)
2176 CALL vol7d_alloc(vol7d_out, nana=nana, ntime=ntime, nlevel=nlevel, &
2177 ntimerange=ntimerange, ndativarr=nvar, nnetwork=nnetwork)
2178 vol7d_out%ana = v7d_locana%ana
2180 CALL get_val(grid_trans, output_point_index=point_index)
2181 IF (
ALLOCATED(point_index))
THEN
2183 CALL vol7d_alloc(vol7d_out, nanavari=1)
2184 CALL init(vol7d_out%anavar%i(1),
'B01192')
2187 CALL vol7d_alloc_vol(vol7d_out)
2189 IF (
ALLOCATED(point_index))
THEN
2190 DO inetwork = 1, nnetwork
2191 vol7d_out%volanai(:,1,inetwork) = point_index(:)
2194 CALL compute(grid_trans, volgrid6d_in, vol7d_out, networkname, noconvert)
2196 CALL l4f_log(l4f_error,
'vg6d_v7d_transform: transformation not valid')
2204CALL vol7d_dballe_set_var_du(vol7d_out)
2209END SUBROUTINE volgrid6d_v7d_transform
2220SUBROUTINE volgrid6dv_v7d_transform(this, volgrid6d_in, vol7d_out, v7d, &
2221 maskgrid, maskbounds, networkname, noconvert, find_index, categoryappend)
2223TYPE(
volgrid6d),
INTENT(inout) :: volgrid6d_in(:)
2224TYPE(
vol7d),
INTENT(out) :: vol7d_out
2225TYPE(
vol7d),
INTENT(in),
OPTIONAL :: v7d
2226REAL,
INTENT(in),
OPTIONAL :: maskgrid(:,:)
2227REAL,
INTENT(in),
OPTIONAL :: maskbounds(:)
2228CHARACTER(len=*),
OPTIONAL,
INTENT(in) :: networkname
2229LOGICAL,
OPTIONAL,
INTENT(in) :: noconvert
2230PROCEDURE(basic_find_index),
POINTER,
OPTIONAL :: find_index
2231CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
2234TYPE(
vol7d) :: v7dtmp
2240DO i=1,
SIZE(volgrid6d_in)
2241 CALL transform(this, volgrid6d_in(i), v7dtmp, v7d=v7d, &
2242 maskgrid=maskgrid, maskbounds=maskbounds, &
2243 networkname=networkname, noconvert=noconvert, find_index=find_index, &
2244 categoryappend=categoryappend)
2245 CALL vol7d_append(vol7d_out, v7dtmp)
2248END SUBROUTINE volgrid6dv_v7d_transform
2252SUBROUTINE v7d_volgrid6d_transform_compute(this, vol7d_in, volgrid6d_out, networkname, gaid_template)
2254type(
vol7d),
INTENT(in) :: vol7d_in
2255type(
volgrid6d),
INTENT(inout) :: volgrid6d_out
2256CHARACTER(len=*),
OPTIONAL,
INTENT(in) :: networkname
2257TYPE(
grid_id),
OPTIONAL,
INTENT(in) :: gaid_template
2259integer :: nana, ntime, ntimerange, nlevel, nvar
2260INTEGER :: ilevel, itime, itimerange, ivar, inetwork
2262REAL,
POINTER :: voldatiout(:,:,:)
2263type(vol7d_network) :: network
2268 'start v7d_volgrid6d_transform_compute')
2276IF (
PRESENT(networkname))
THEN
2277 CALL init(network,name=networkname)
2278 inetwork =
index(vol7d_in%network,network)
2279 IF (inetwork <= 0)
THEN
2281 'network '//trim(networkname)//
' not found, first network will be transformed')
2289if (
associated(vol7d_in%time))
then
2290 ntime=
size(vol7d_in%time)
2291 volgrid6d_out%time=vol7d_in%time
2294if (
associated(vol7d_in%timerange))
then
2295 ntimerange=
size(vol7d_in%timerange)
2296 volgrid6d_out%timerange=vol7d_in%timerange
2299if (
associated(vol7d_in%level))
then
2300 nlevel=
size(vol7d_in%level)
2301 volgrid6d_out%level=vol7d_in%level
2304if (
associated(vol7d_in%dativar%r))
then
2305 nvar=
size(vol7d_in%dativar%r)
2306 CALL varbufr2vargrib(vol7d_in%dativar%r, volgrid6d_out%var, c_func, gaid_template)
2309nana=
SIZE(vol7d_in%voldatir, 1)
2311IF (.NOT.
ASSOCIATED(volgrid6d_out%voldati))
THEN
2312 ALLOCATE(voldatiout(volgrid6d_out%griddim%dim%nx, volgrid6d_out%griddim%dim%ny, &
2317 DO itimerange=1,ntimerange
2321 IF (
PRESENT(gaid_template))
THEN
2322 DO ilevel = 1, nlevel
2323 IF (any(
c_e(vol7d_in%voldatir(:,itime,ilevel,itimerange,ivar,inetwork))))
THEN
2324 CALL copy(gaid_template, volgrid6d_out%gaid(ilevel,itime,itimerange,ivar))
2326 volgrid6d_out%gaid(ilevel,itime,itimerange,ivar) = grid_id_new()
2332 IF (
ASSOCIATED(volgrid6d_out%voldati)) &
2333 CALL volgrid_get_vol_3d(volgrid6d_out, itime, itimerange, ivar, &
2336 CALL compute(this, &
2337 vol7d_in%voldatir(:,itime,:,itimerange,ivar,inetwork), voldatiout, &
2338 vol7d_in%dativar%r(ivar))
2340 IF (
ASSOCIATED(c_func))
THEN
2341 CALL compute(c_func(ivar), voldatiout(:,:,:))
2344 CALL volgrid_set_vol_3d(volgrid6d_out, itime, itimerange, ivar, &
2351IF (.NOT.
ASSOCIATED(volgrid6d_out%voldati))
THEN
2352 DEALLOCATE(voldatiout)
2354IF (
ASSOCIATED(c_func))
THEN
2358END SUBROUTINE v7d_volgrid6d_transform_compute
2367SUBROUTINE v7d_volgrid6d_transform(this, griddim, vol7d_in, volgrid6d_out, &
2368 networkname, gaid_template, categoryappend)
2372TYPE(
vol7d),
INTENT(inout) :: vol7d_in
2373TYPE(
volgrid6d),
INTENT(out) :: volgrid6d_out
2374CHARACTER(len=*),
OPTIONAL,
INTENT(in) :: networkname
2375TYPE(
grid_id),
OPTIONAL,
INTENT(in) :: gaid_template
2376CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
2379integer :: ntime, ntimerange, nlevel, nvar
2385CALL vol7d_alloc_vol(vol7d_in)
2386ntime=
SIZE(vol7d_in%time)
2387ntimerange=
SIZE(vol7d_in%timerange)
2388nlevel=
SIZE(vol7d_in%level)
2390if (
associated(vol7d_in%dativar%r)) nvar=
size(vol7d_in%dativar%r)
2393 CALL l4f_log(l4f_error, &
2394 "trying to transform a vol7d object incomplete or without real variables")
2395 CALL init(volgrid6d_out)
2400CALL init(grid_trans, this, vol7d_in, griddim, categoryappend=categoryappend)
2401CALL init(volgrid6d_out, griddim, time_definition=vol7d_in%time_definition, &
2402 categoryappend=categoryappend)
2404IF (
c_e(grid_trans))
THEN
2406 CALL volgrid6d_alloc(volgrid6d_out, griddim%dim, ntime=ntime, nlevel=nlevel, &
2407 ntimerange=ntimerange, nvar=nvar)
2409 CALL volgrid6d_alloc_vol(volgrid6d_out, decode=.true.)
2411 CALL compute(grid_trans, vol7d_in, volgrid6d_out, networkname, gaid_template)
2413 CALL vg6d_wind_rot(volgrid6d_out)
2416 CALL l4f_log(l4f_error,
'v7d_vg6d_transform: transformation not valid')
2422END SUBROUTINE v7d_volgrid6d_transform
2426SUBROUTINE v7d_v7d_transform_compute(this, vol7d_in, vol7d_out, lev_out, &
2429type(
vol7d),
INTENT(in) :: vol7d_in
2430type(
vol7d),
INTENT(inout) :: vol7d_out
2431TYPE(
vol7d_level),
INTENT(in),
OPTIONAL :: lev_out(:)
2432INTEGER,
INTENT(in),
OPTIONAL :: var_coord_vol
2434INTEGER :: itime, itimerange, ilevel, ivar, inetwork, &
2435 levshift, levused, lvar_coord_vol, spos
2436REAL,
ALLOCATABLE :: coord_3d_in(:,:,:)
2439lvar_coord_vol = optio_i(var_coord_vol)
2440vol7d_out%time(:) = vol7d_in%time(:)
2441vol7d_out%timerange(:) = vol7d_in%timerange(:)
2442IF (
PRESENT(lev_out))
THEN
2443 vol7d_out%level(:) = lev_out(:)
2445 vol7d_out%level(:) = vol7d_in%level(:)
2447vol7d_out%network(:) = vol7d_in%network(:)
2448IF (
ASSOCIATED(vol7d_in%dativar%r))
THEN
2449 vol7d_out%dativar%r(:) = vol7d_in%dativar%r(:)
2451 CALL get_val(this, levshift=levshift, levused=levused)
2453 IF (
c_e(lvar_coord_vol))
THEN
2454 CALL get_val(this%trans, output_levtype=output_levtype)
2455 IF (output_levtype%level1 == 103 .OR. output_levtype%level1 == 108)
THEN
2456 spos = firsttrue(vol7d_in%level(:) == vol7d_level_new(1))
2458 CALL l4f_log(l4f_error, &
2459 'output level '//
t2c(output_levtype%level1)// &
2460 ' requested, but height/press of surface not provided in volume')
2462 IF (.NOT.
c_e(levshift) .AND. .NOT.
c_e(levused))
THEN
2463 CALL l4f_log(l4f_error, &
2464 'internal inconsistence, levshift and levused undefined when they should be')
2466 ALLOCATE(coord_3d_in(
SIZE(vol7d_in%ana),1,
SIZE(vol7d_in%level)))
2471 DO inetwork = 1,
SIZE(vol7d_in%network)
2472 DO ivar = 1,
SIZE(vol7d_in%dativar%r)
2473 DO itimerange = 1,
SIZE(vol7d_in%timerange)
2474 DO itime = 1,
SIZE(vol7d_in%time)
2477 IF (
c_e(lvar_coord_vol))
THEN
2480 coord_3d_in(:,:,levshift+1:levshift+levused) = rmiss
2482 DO ilevel = levshift+1, levshift+levused
2483 WHERE(
c_e(vol7d_in%voldatir(:,itime:itime,ilevel,itimerange,lvar_coord_vol,inetwork)) .AND. &
2484 c_e(vol7d_in%voldatir(:,itime:itime,spos,itimerange,lvar_coord_vol,inetwork)))
2485 coord_3d_in(:,:,ilevel) = vol7d_in%voldatir(:,itime:itime,ilevel,itimerange,lvar_coord_vol,inetwork) - &
2486 vol7d_in%voldatir(:,itime:itime,spos,itimerange,lvar_coord_vol,inetwork)
2488 coord_3d_in(:,:,ilevel) = rmiss
2492 CALL compute(this, &
2493 vol7d_in%voldatir(:,itime,:,itimerange,ivar,inetwork), &
2494 vol7d_out%voldatir(:,itime:itime,:,itimerange,ivar,inetwork), &
2495 var=vol7d_in%dativar%r(ivar), &
2496 coord_3d_in=coord_3d_in)
2498 CALL compute(this, &
2499 vol7d_in%voldatir(:,itime,:,itimerange,ivar,inetwork), &
2500 vol7d_out%voldatir(:,itime:itime,:,itimerange,ivar,inetwork), &
2501 var=vol7d_in%dativar%r(ivar), &
2502 coord_3d_in=vol7d_in%voldatir(:,itime:itime,:,itimerange, &
2503 lvar_coord_vol,inetwork))
2506 CALL compute(this, &
2507 vol7d_in%voldatir(:,itime,:,itimerange,ivar,inetwork), &
2508 vol7d_out%voldatir(:,itime:itime,:,itimerange,ivar,inetwork), &
2509 var=vol7d_in%dativar%r(ivar))
2518END SUBROUTINE v7d_v7d_transform_compute
2528SUBROUTINE v7d_v7d_transform(this, vol7d_in, vol7d_out, v7d, maskbounds, &
2529 lev_out, vol7d_coord_in, categoryappend)
2531TYPE(
vol7d),
INTENT(inout) :: vol7d_in
2532TYPE(
vol7d),
INTENT(out) :: vol7d_out
2533TYPE(
vol7d),
INTENT(in),
OPTIONAL :: v7d
2534REAL,
INTENT(in),
OPTIONAL :: maskbounds(:)
2535TYPE(
vol7d_level),
INTENT(in),
OPTIONAL,
TARGET :: lev_out(:)
2536TYPE(
vol7d),
INTENT(in),
OPTIONAL :: vol7d_coord_in
2537CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
2539INTEGER :: nvar, inetwork
2543TYPE(vol7d_var) :: vcoord_var
2544REAL,
ALLOCATABLE :: coord_3d_in(:,:,:)
2545INTEGER :: var_coord_in, var_coord_vol, i, k, ulstart, ulend, spos
2546INTEGER,
ALLOCATABLE :: point_index(:)
2547TYPE(
vol7d) :: v7d_locana, vol7d_tmpana
2548CHARACTER(len=80) :: trans_type
2549LOGICAL,
ALLOCATABLE :: mask_in(:), point_mask(:)
2551CALL vol7d_alloc_vol(vol7d_in)
2553IF (
ASSOCIATED(vol7d_in%dativar%r)) nvar=
SIZE(vol7d_in%dativar%r)
2556IF (
PRESENT(v7d)) v7d_locana = v7d
2557CALL init(vol7d_out, time_definition=vol7d_in%time_definition)
2559CALL get_val(this, trans_type=trans_type)
2561var_coord_vol = imiss
2562IF (trans_type ==
'vertint')
THEN
2564 IF (
PRESENT(lev_out))
THEN
2568 IF (
PRESENT(vol7d_coord_in))
THEN
2569 IF (
ASSOCIATED(vol7d_coord_in%voldatir) .AND. &
2570 ASSOCIATED(vol7d_coord_in%dativar%r))
THEN
2573 IF (
SIZE(vol7d_coord_in%voldatir,2) /= 1 .OR. &
2574 SIZE(vol7d_coord_in%voldatir,4) /= 1 .OR. &
2575 SIZE(vol7d_coord_in%voldatir,6) /= 1)
THEN
2576 CALL l4f_log(l4f_error, &
2577 'volume providing constant input vertical coordinate must have &
2578 &only 1 time, 1 timerange and 1 network')
2584 CALL get_val(this, output_levtype=output_levtype)
2586 IF (.NOT.
c_e(vcoord_var))
THEN
2587 CALL l4f_log(l4f_error, &
2588 'requested output level type '//
t2c(output_levtype%level1)// &
2589 ' does not correspond to any known physical variable for &
2590 &providing vertical coordinate')
2595 var_coord_in =
index(vol7d_coord_in%dativar%r, vcoord_var)
2597 IF (var_coord_in <= 0)
THEN
2598 CALL l4f_log(l4f_error, &
2599 'volume providing constant input vertical coordinate contains no &
2600 &real variables matching output level type '//
t2c(output_levtype%level1))
2604 CALL l4f_log(l4f_info, &
2605 'Coordinate for vertint found in coord volume at position '// &
2609 CALL get_val(this, input_levtype=input_levtype)
2611 (vol7d_coord_in%level(:)%level1 == input_levtype%level1) .AND. &
2612 (vol7d_coord_in%level(:)%level2 == input_levtype%level2)
2613 ulstart = firsttrue(mask_in)
2614 ulend = lasttrue(mask_in)
2615 IF (ulstart == 0 .OR. ulend == 0)
THEN
2616 CALL l4f_log(l4f_error, &
2617 'coordinate file does not contain levels of type '// &
2618 t2c(input_levtype%level1)//
'/'//
t2c(input_levtype%level2)// &
2619 ' specified for input data')
2624 coord_3d_in = vol7d_coord_in%voldatir(:,1:1,ulstart:ulend,1,var_coord_in,1)
2626 IF (output_levtype%level1 == 103 &
2627 .OR. output_levtype%level1 == 108)
THEN
2628 spos = firsttrue(vol7d_coord_in%level(:) == vol7d_level_new(1))
2630 CALL l4f_log(l4f_error, &
2631 'output level '//
t2c(output_levtype%level1)// &
2632 ' requested, but height/press of surface not provided in coordinate file')
2636 DO k = 1,
SIZE(coord_3d_in,3)
2637 WHERE(
c_e(coord_3d_in(:,:,k)) .AND. &
2638 c_e(vol7d_coord_in%voldatir(:,1:1,spos,1,var_coord_in,1)))
2639 coord_3d_in(:,:,k) = coord_3d_in(:,:,k) - &
2640 vol7d_coord_in%voldatir(:,1:1,spos,1,var_coord_in,1)
2642 coord_3d_in(:,:,k) = rmiss
2650 IF (var_coord_in <= 0)
THEN
2652 CALL get_val(this, output_levtype=output_levtype)
2654 IF (
c_e(vcoord_var))
THEN
2655 DO i = 1,
SIZE(vol7d_in%dativar%r)
2656 IF (vol7d_in%dativar%r(i) == vcoord_var)
THEN
2662 IF (
c_e(var_coord_vol))
THEN
2663 CALL l4f_log(l4f_info, &
2664 'Coordinate for vertint found in input volume at position '// &
2671 IF (var_coord_in > 0)
THEN
2672 CALL init(grid_trans, this, lev_in=vol7d_in%level, lev_out=lev_out, &
2673 coord_3d_in=coord_3d_in, categoryappend=categoryappend)
2675 CALL init(grid_trans, this, lev_in=vol7d_in%level, lev_out=lev_out, &
2676 categoryappend=categoryappend)
2679 CALL get_val(grid_trans, output_level_auto=llev_out)
2680 IF (.NOT.
associated(llev_out)) llev_out => lev_out
2682 IF (
c_e(grid_trans)) then
2684 CALL vol7d_alloc(vol7d_out, nana=
SIZE(vol7d_in%ana), &
2685 ntime=
SIZE(vol7d_in%time), ntimerange=
SIZE(vol7d_in%timerange), &
2686 nlevel=
SIZE(llev_out), nnetwork=
SIZE(vol7d_in%network), ndativarr=nvar)
2687 vol7d_out%ana(:) = vol7d_in%ana(:)
2689 CALL vol7d_alloc_vol(vol7d_out)
2694 CALL compute(grid_trans, vol7d_in, vol7d_out, llev_out, &
2695 var_coord_vol=var_coord_vol)
2697 CALL l4f_log(l4f_error,
'v7d_v7d_transform: transformation not valid')
2701 CALL l4f_log(l4f_error, &
2702 'v7d_v7d_transform: vertint requested but lev_out not provided')
2708 CALL init(grid_trans, this, vol7d_in, v7d_locana, maskbounds=maskbounds, &
2709 categoryappend=categoryappend)
2712 IF (
c_e(grid_trans)) then
2714 CALL vol7d_alloc(vol7d_out, nana=
SIZE(v7d_locana%ana), &
2715 ntime=
SIZE(vol7d_in%time), ntimerange=
SIZE(vol7d_in%timerange), &
2716 nlevel=
SIZE(vol7d_in%level), nnetwork=
SIZE(vol7d_in%network), ndativarr=nvar)
2717 vol7d_out%ana = v7d_locana%ana
2719 CALL get_val(grid_trans, point_mask=point_mask, output_point_index=point_index)
2721 IF (
ALLOCATED(point_index))
THEN
2722 CALL vol7d_alloc(vol7d_out, nanavari=1)
2723 CALL init(vol7d_out%anavar%i(1),
'B01192')
2726 CALL vol7d_alloc_vol(vol7d_out)
2728 IF (
ALLOCATED(point_index))
THEN
2729 DO inetwork = 1,
SIZE(vol7d_in%network)
2730 vol7d_out%volanai(:,1,inetwork) = point_index(:)
2733 CALL compute(grid_trans, vol7d_in, vol7d_out)
2735 IF (
ALLOCATED(point_mask))
THEN
2736 IF (
SIZE(point_mask) /=
SIZE(vol7d_in%ana))
THEN
2737 CALL l4f_log(l4f_warn, &
2738 'v7d_v7d_transform: inconsistency in point size: '//
t2c(
SIZE(point_mask)) &
2739 //
':'//
t2c(
SIZE(vol7d_in%ana)))
2742 CALL l4f_log(l4f_debug,
'v7d_v7d_transform: merging ana from in to out')
2744 CALL vol7d_copy(vol7d_in, vol7d_tmpana, &
2745 lana=point_mask, lnetwork=(/.true./), &
2746 ltime=(/.false./), ltimerange=(/.false./), llevel=(/.false./))
2747 CALL vol7d_append(vol7d_out, vol7d_tmpana)
2752 CALL l4f_log(l4f_error,
'v7d_v7d_transform: transformation not valid')
2759IF (.NOT.
PRESENT(v7d))
CALL delete(v7d_locana)
2761END SUBROUTINE v7d_v7d_transform
2771subroutine vg6d_wind_unrot(this)
2774integer :: component_flag
2776call get_val(this%griddim,component_flag=component_flag)
2778if (component_flag == 1)
then
2780 "unrotating vector components")
2781 call vg6d_wind__un_rot(this,.false.)
2782 call set_val(this%griddim,component_flag=0)
2785 "no need to unrotate vector components")
2788end subroutine vg6d_wind_unrot
2796subroutine vg6d_wind_rot(this)
2799integer :: component_flag
2801call get_val(this%griddim,component_flag=component_flag)
2803if (component_flag == 0)
then
2805 "rotating vector components")
2806 call vg6d_wind__un_rot(this,.true.)
2807 call set_val(this%griddim,component_flag=1)
2810 "no need to rotate vector components")
2813end subroutine vg6d_wind_rot
2817SUBROUTINE vg6d_wind__un_rot(this,rot)
2821INTEGER :: i, j, k, l, a11, a12, a21, a22, stallo
2822double precision,
pointer :: rot_mat(:,:,:)
2823real,
allocatable :: tmp_arr(:,:)
2824REAL,
POINTER :: voldatiu(:,:), voldativ(:,:)
2825INTEGER,
POINTER :: iu(:), iv(:)
2827IF (.NOT.
ASSOCIATED(this%var))
THEN
2829 "trying to unrotate an incomplete volgrid6d object")
2830 CALL raise_fatal_error()
2834CALL volgrid6d_var_hor_comp_index(this%var, iu, iv)
2835IF (.NOT.
ASSOCIATED(iu))
THEN
2837 "unrotation impossible")
2838 CALL raise_fatal_error()
2843ALLOCATE(tmp_arr(this%griddim%dim%nx, this%griddim%dim%ny),stat=stallo)
2844IF (stallo /= 0)
THEN
2846 CALL raise_fatal_error()
2849IF (.NOT.
ASSOCIATED(this%voldati))
THEN
2850 ALLOCATE(voldatiu(this%griddim%dim%nx, this%griddim%dim%ny), &
2851 voldativ(this%griddim%dim%nx, this%griddim%dim%ny))
2854CALL griddim_unproj(this%griddim)
2855CALL wind_unrot(this%griddim, rot_mat)
2868 DO k = 1,
SIZE(this%timerange)
2869 DO j = 1,
SIZE(this%time)
2870 DO i = 1,
SIZE(this%level)
2872 CALL volgrid_get_vol_2d(this, i, j, k, iu(l), voldatiu)
2873 CALL volgrid_get_vol_2d(this, i, j, k, iv(l), voldativ)
2879 WHERE(voldatiu /= rmiss .AND. voldativ /= rmiss)
2880 tmp_arr(:,:) = real(voldatiu(:,:)*rot_mat(:,:,a11) + &
2881 voldativ(:,:)*rot_mat(:,:,a12))
2882 voldativ(:,:) = real(voldatiu(:,:)*rot_mat(:,:,a21) + &
2883 voldativ(:,:)*rot_mat(:,:,a22))
2884 voldatiu(:,:) = tmp_arr(:,:)
2890 CALL volgrid_set_vol_2d(this, i, j, k, iu(l), voldatiu)
2891 CALL volgrid_set_vol_2d(this, i, j, k, iv(l), voldativ)
2897IF (.NOT.
ASSOCIATED(this%voldati))
THEN
2898 DEALLOCATE(voldatiu, voldativ)
2900DEALLOCATE(rot_mat, tmp_arr, iu, iv)
2902END SUBROUTINE vg6d_wind__un_rot
2948subroutine vg6d_c2a (this)
2952integer :: ngrid,igrid,jgrid,ugrid,vgrid,tgrid
2953doubleprecision :: xmin, xmax, ymin, ymax
2954doubleprecision :: xmin_t, xmax_t, ymin_t, ymax_t
2955doubleprecision :: step_lon_t,step_lat_t
2956character(len=80) :: type_t,type
2963 call init(griddim_t)
2965 call get_val(this(igrid)%griddim,xmin=xmin_t, xmax=xmax_t, ymin=ymin_t, ymax=ymax_t,proj_type=type_t)
2966 step_lon_t=(xmax_t-xmin_t)/dble(this(igrid)%griddim%dim%nx-1)
2967 step_lat_t=(ymax_t-ymin_t)/dble(this(igrid)%griddim%dim%ny-1)
2980 if (this(igrid)%griddim == this(jgrid)%griddim ) cycle
2982 if (this(igrid)%griddim%dim%nx == this(jgrid)%griddim%dim%nx .and. &
2983 this(igrid)%griddim%dim%ny == this(jgrid)%griddim%dim%ny )
then
2985 call get_val(this(jgrid)%griddim,xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,proj_type=type)
2987 if (type_t /=
type )cycle
2993 call l4f_category_log(this(igrid)%category,l4f_debug,
"diff coordinate lon"//&
2994 to_char(
abs(xmin - (xmin_t+step_lon_t/2.d0)))//&
2995 to_char(
abs(xmax - (xmax_t+step_lon_t/2.d0))))
2996 call l4f_category_log(this(igrid)%category,l4f_debug,
"diff coordinate lat"//&
2997 to_char(
abs(ymin - (ymin_t+step_lat_t/2.d0)))//&
2998 to_char(
abs(ymax - (ymax_t+step_lat_t/2.d0))))
3001 if (
abs(xmin - (xmin_t+step_lon_t/2.d0)) < 1.d-3 .and.
abs(xmax - (xmax_t+step_lon_t/2.d0)) < 1.d-3 )
then
3002 if (
abs(ymin - ymin_t) < 1.d-3 .and.
abs(ymax - ymax_t) < 1.d-3 )
then
3018 if (
abs(ymin - (ymin_t+step_lat_t/2.d0)) < 1.d-3 .and.
abs(ymax - (ymax_t+step_lat_t/2.d0)) < 1.d-3 )
then
3019 if (
abs(xmin - xmin_t) < 1.d-3 .and.
abs(xmax - xmax_t) < 1.d-3 )
then
3034 call l4f_category_log(this(igrid)%category,l4f_debug,
"C grid: test U and V"//&
3038 call l4f_category_log(this(igrid)%category,l4f_debug,
"UV diff coordinate lon"//&
3039 to_char(
abs(xmin_t - xmin)-step_lon_t/2.d0)//&
3041 call l4f_category_log(this(igrid)%category,l4f_debug,
"UV diff coordinate lat"//&
3042 to_char(
abs(ymin_t - ymin) -step_lat_t/2.d0)//&
3046 if (
abs(ymin - (ymin_t+step_lat_t/2.d0)) < 2.d-3 .and.
abs(ymax - (ymax_t+step_lat_t/2.d0)) < 2.d-3 )
then
3047 if (
abs(xmin_t - (xmin+step_lon_t/2.d0)) < 2.d-3 .and.
abs(xmax_t - (xmax+step_lon_t/2.d0)) < 2.d-3 )
then
3050 call l4f_category_log(this(igrid)%category,l4f_debug,
"C grid: found U and V case up and right")
3056 call init(griddim_t,xmin=xmin, xmax=xmax, ymin=ymin_t, ymax=ymax_t)
3063 if (
c_e(ugrid))
then
3068 call vg6d_c2a_grid(this(ugrid),this(tgrid)%griddim,cgrid=1)
3070 call vg6d_c2a_grid(this(ugrid),griddim_t,cgrid=1)
3072 call vg6d_c2a_mat(this(ugrid),cgrid=1)
3075 if (
c_e(vgrid))
then
3080 call vg6d_c2a_grid(this(vgrid),this(tgrid)%griddim,cgrid=2)
3082 call vg6d_c2a_grid(this(vgrid),griddim_t,cgrid=2)
3084 call vg6d_c2a_mat(this(vgrid),cgrid=2)
3094end subroutine vg6d_c2a
3098subroutine vg6d_c2a_grid(this,griddim_t,cgrid)
3102integer,
intent(in) :: cgrid
3104doubleprecision :: xmin, xmax, ymin, ymax
3105doubleprecision :: step_lon,step_lat
3108if (
present(griddim_t))
then
3110 call get_val(griddim_t,xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
3111 call set_val(this%griddim,xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
3113 CALL griddim_setsteps(this%griddim)
3122 call l4f_category_log(this%category,l4f_debug,
"C grid: T points, nothing to do")
3129 call l4f_category_log(this%category,l4f_debug,
"C grid: U points, we need interpolation")
3132 call get_val(this%griddim, xmin=xmin, xmax=xmax)
3133 step_lon=(xmax-xmin)/dble(this%griddim%dim%nx-1)
3134 xmin=xmin-step_lon/2.d0
3135 xmax=xmax-step_lon/2.d0
3138 CALL griddim_setsteps(this%griddim)
3143 call l4f_category_log(this%category,l4f_debug,
"C grid: V points, we need interpolation")
3146 call get_val(this%griddim, ymin=ymin, ymax=ymax)
3147 step_lat=(ymax-ymin)/dble(this%griddim%dim%ny-1)
3148 ymin=ymin-step_lat/2.d0
3149 ymax=ymax-step_lat/2.d0
3150 call set_val(this%griddim, ymin=ymin, ymax=ymax)
3152 CALL griddim_setsteps(this%griddim)
3157 call raise_fatal_error ()
3164call griddim_unproj(this%griddim)
3167end subroutine vg6d_c2a_grid
3170subroutine vg6d_c2a_mat(this,cgrid)
3173integer,
intent(in) :: cgrid
3175INTEGER :: i, j, k, iv, stallo
3176REAL,
ALLOCATABLE :: tmp_arr(:,:)
3177REAL,
POINTER :: voldatiuv(:,:)
3180IF (cgrid == 0)
RETURN
3181IF (cgrid /= 1 .AND. cgrid /= 2)
THEN
3183 trim(
to_char(cgrid))//
" not known")
3189ALLOCATE(tmp_arr(this%griddim%dim%nx, this%griddim%dim%ny),stat=stallo)
3191 call l4f_log(l4f_fatal,
"allocating memory")
3192 call raise_fatal_error()
3196IF (.NOT.
ASSOCIATED(this%voldati))
THEN
3197 ALLOCATE(voldatiuv(this%griddim%dim%nx, this%griddim%dim%ny), stat=stallo)
3198 IF (stallo /= 0)
THEN
3199 CALL l4f_log(l4f_fatal,
"allocating memory")
3200 CALL raise_fatal_error()
3205 DO iv = 1,
SIZE(this%var)
3206 DO k = 1,
SIZE(this%timerange)
3207 DO j = 1,
SIZE(this%time)
3208 DO i = 1,
SIZE(this%level)
3209 tmp_arr(:,:) = rmiss
3210 CALL volgrid_get_vol_2d(this, i, j, k, iv, voldatiuv)
3213 WHERE(voldatiuv(1,:) /= rmiss .AND. voldatiuv(2,:) /= rmiss)
3214 tmp_arr(1,:) = voldatiuv(1,:) - (voldatiuv(2,:) - voldatiuv(1,:)) / 2.
3218 WHERE(voldatiuv(1:this%griddim%dim%nx-1,:) /= rmiss .AND. &
3219 voldatiuv(2:this%griddim%dim%nx,:) /= rmiss)
3220 tmp_arr(2:this%griddim%dim%nx,:) = &
3221 (voldatiuv(1:this%griddim%dim%nx-1,:) + &
3222 voldatiuv(2:this%griddim%dim%nx,:)) / 2.
3225 voldatiuv(:,:) = tmp_arr
3226 CALL volgrid_set_vol_2d(this, i, j, k, iv, voldatiuv)
3232ELSE IF (cgrid == 2)
THEN
3233 DO iv = 1,
SIZE(this%var)
3234 DO k = 1,
SIZE(this%timerange)
3235 DO j = 1,
SIZE(this%time)
3236 DO i = 1,
SIZE(this%level)
3237 tmp_arr(:,:) = rmiss
3238 CALL volgrid_get_vol_2d(this, i, j, k, iv, voldatiuv)
3241 WHERE(voldatiuv(:,1) /= rmiss .AND. voldatiuv(:,2) /= rmiss)
3242 tmp_arr(:,1) = voldatiuv(:,1) - (voldatiuv(:,2) - voldatiuv(:,1)) / 2.
3246 WHERE(voldatiuv(:,1:this%griddim%dim%ny-1) /= rmiss .AND. &
3247 voldatiuv(:,2:this%griddim%dim%ny) /= rmiss)
3248 tmp_arr(:,2:this%griddim%dim%ny) = &
3249 (voldatiuv(:,1:this%griddim%dim%ny-1) + &
3250 voldatiuv(:,2:this%griddim%dim%ny)) / 2.
3253 voldatiuv(:,:) = tmp_arr
3254 CALL volgrid_set_vol_2d(this, i, j, k, iv, voldatiuv)
3261IF (.NOT.
ASSOCIATED(this%voldati))
THEN
3262 DEALLOCATE(voldatiuv)
3266end subroutine vg6d_c2a_mat
3272subroutine display_volgrid6d (this)
3281print*,
"----------------------- volgrid6d display ---------------------"
3284IF (
ASSOCIATED(this%time))
then
3285 print*,
"---- time vector ----"
3286 print*,
"elements=",
size(this%time)
3287 do i=1,
size(this%time)
3292IF (
ASSOCIATED(this%timerange))
then
3293 print*,
"---- timerange vector ----"
3294 print*,
"elements=",
size(this%timerange)
3295 do i=1,
size(this%timerange)
3296 call display(this%timerange(i))
3300IF (
ASSOCIATED(this%level))
then
3301 print*,
"---- level vector ----"
3302 print*,
"elements=",
size(this%level)
3303 do i=1,
size(this%level)
3308IF (
ASSOCIATED(this%var))
then
3309 print*,
"---- var vector ----"
3310 print*,
"elements=",
size(this%var)
3311 do i=1,
size(this%var)
3316IF (
ASSOCIATED(this%gaid))
then
3317 print*,
"---- gaid vector (present mask only) ----"
3318 print*,
"elements=",shape(this%gaid)
3319 print*,
c_e(reshape(this%gaid,(/
SIZE(this%gaid)/)))
3322print*,
"--------------------------------------------------------------"
3325end subroutine display_volgrid6d
3331subroutine display_volgrid6dv (this)
3336print*,
"----------------------- volgrid6d vector ---------------------"
3338print*,
"elements=",
size(this)
3343 call l4f_category_log(this(i)%category,l4f_debug,
"ora mostro il vettore volgrid6d" )
3349print*,
"--------------------------------------------------------------"
3351end subroutine display_volgrid6dv
3356subroutine vg6dv_rounding(vg6din,vg6dout,level,timerange,nostatproc,merge)
3358type(
volgrid6d),
intent(out),
pointer :: vg6dout(:)
3361logical,
intent(in),
optional :: merge
3362logical,
intent(in),
optional :: nostatproc
3366allocate(vg6dout(
size(vg6din)))
3368do i = 1,
size(vg6din)
3369 call vg6d_rounding(vg6din(i),vg6dout(i),level,timerange,nostatproc,merge)
3372end subroutine vg6dv_rounding
3385subroutine vg6d_rounding(vg6din,vg6dout,level,timerange,nostatproc,merge)
3390logical,
intent(in),
optional :: merge
3392logical,
intent(in),
optional :: nostatproc
3394integer :: ilevel,itimerange
3395type(
vol7d_level) :: roundlevel(size(vg6din%level))
3398roundlevel=vg6din%level
3400if (
present(level))
then
3401 do ilevel = 1,
size(vg6din%level)
3402 if ((any(vg6din%level(ilevel) .almosteq. level)))
then
3403 roundlevel(ilevel)=level(1)
3408roundtimerange=vg6din%timerange
3410if (
present(timerange))
then
3411 do itimerange = 1,
size(vg6din%timerange)
3412 if ((any(vg6din%timerange(itimerange) .almosteq. timerange)))
then
3413 roundtimerange(itimerange)=timerange(1)
3420if (optio_log(nostatproc))
then
3421 roundtimerange(:)%timerange=254
3422 roundtimerange(:)%p2=0
3426call vg6d_reduce(vg6din,vg6dout,roundlevel,roundtimerange,merge)
3428end subroutine vg6d_rounding
3438subroutine vg6d_reduce(vg6din,vg6dout,roundlevel,roundtimerange,merge)
3443logical,
intent(in),
optional :: merge
3445integer :: nlevel,ntime,ntimerange,nvar,ilevel,itimerange,ivar,indl,indt,itime,nx,ny
3446real,
allocatable :: vol2d(:,:)
3448nx=vg6din%griddim%dim%nx
3449ny=vg6din%griddim%dim%ny
3450nlevel=count_distinct(roundlevel,back=.true.)
3451ntime=
size(vg6din%time)
3452ntimerange=count_distinct(roundtimerange,back=.true.)
3453nvar=
size(vg6din%var)
3455call init(vg6dout, vg6din%griddim, vg6din%time_definition, categoryappend=
"generated by vg6d_reduce")
3456call volgrid6d_alloc(vg6dout, vg6din%griddim%dim, ntime, nlevel, ntimerange, nvar)
3458if (
ASSOCIATED(vg6din%voldati) .or. optio_log(merge))
then
3459 call volgrid6d_alloc_vol(vg6dout,inivol=.true.,decode=.true.)
3460 allocate(vol2d(nx,ny))
3462 call volgrid6d_alloc_vol(vg6dout,inivol=.true.,decode=.false.)
3465vg6dout%time=vg6din%time
3466vg6dout%var=vg6din%var
3467vg6dout%timerange=pack_distinct(roundtimerange,ntimerange,back=.true.)
3468vg6dout%level=pack_distinct(roundlevel,nlevel,back=.true.)
3470CALL sort(vg6dout%timerange)
3471CALL sort(vg6dout%level)
3473do ilevel=1,
size(vg6din%level)
3474 indl=
index(vg6dout%level,roundlevel(ilevel))
3475 do itimerange=1,
size(vg6din%timerange)
3476 indt=
index(vg6dout%timerange,roundtimerange(itimerange))
3480 if (
ASSOCIATED(vg6din%voldati))
then
3481 vol2d=vg6din%voldati(:,:,ilevel,itime,itimerange,ivar)
3484 if (optio_log(merge))
then
3486 if ( .not.
ASSOCIATED(vg6din%voldati))
then
3487 CALL grid_id_decode_data(vg6din%gaid(ilevel,itime,itimerange,ivar), vol2d)
3491 where (.not.
c_e(vg6dout%voldati(:,:,indl,itime,indt,ivar)))
3493 vg6dout%voldati(:,:,indl,itime,indt,ivar)=vol2d
3496 else if (
ASSOCIATED(vg6din%voldati))
then
3497 if (.not. any(
c_e(vg6dout%voldati(:,:,indl,itime,indt,ivar))))
then
3498 vg6dout%voldati(:,:,indl,itime,indt,ivar)=vol2d
3502 if (
c_e(vg6din%gaid(ilevel,itime,itimerange,ivar)).and. .not.
c_e(vg6dout%gaid(indl,itime,indt,ivar)))
then
3503 call copy (vg6din%gaid(ilevel,itime,itimerange,ivar), vg6dout%gaid(indl,itime,indt,ivar))
3510if (
ASSOCIATED(vg6din%voldati) .or. optio_log(merge))
then
3514end subroutine vg6d_reduce